In this first example we want to solve the Laplace Equation (2) a special case of the Poisson Equation (1) for the absence of any charges.
Prior to actually solving the PDE we have to define a mesh (or grid), on which the equation shall be solved, and a couple of boundary conditions. Lines 6-9 define some support variables and a 2D mesh. Lines 10-12 define further support variables, which are used later to define boundary conditions and save the result in NumPy arrays.
Lines 14-21 define the physical quantitiy (potential in the code, \phi in the above equaitons), the equation with this quantity, that needs to be solved, and the boundary conditions. In this example the boundary conditions represent three capacitor plates thus that the mesh shows the space between the plates. One plate is located at the top and extends along the entire x-range. The two other plates are at the left and right hand side, respectively, cover only the y<y0 range.
As shown in this example, boundary conditions (BCs) are objects of type FixedValue or FixedFlux (not shown). If more than one BC need to be defined, this can be done in a Python list or tuple. But notice: FiPy takes the term boundary condition very seriously. Boundary conditions can only be defined on so called exterior faces. Thus, if boundary conditions for a given problem cannot be applied onto the surfaces of a "common" mesh (line, square or cube), one has to construct a mesh suitable for the problem using a tool called gmsh .
Lines 26-27 store the solution in a file. The specified file extensions determines the format of the output. Omitting arguments to viewer.plot() starts an interactive display module (e.g. Matplotlib, Mayavi2, etc.).
In case of "common" meshes it is rather easy to convert the CellVariable into a NumPy array as is shown in lines 31-39. In case of a uncommon mesh one would have to call the __call__() method of the CellVariable and evaluate each point for a given coordinate (will be shown in another example).
#!/usr/bin/env python from fipy import * import sys nx = 300 dx = 0.05 L = nx * dx mesh = Grid2D(dx=dx,dy=dx,nx=nx,ny=nx) x,y = mesh.getCellCenters() x0,y0 = L/2,L/2 X,Y = mesh.getFaceCenters() potential = CellVariable(mesh=mesh, name='potential', value=0.) potential.equation = (DiffusionTerm(coeff = 1.) == 0.) bcs = ( FixedValue(value=5,faces=mesh.getFacesLeft() & (Y<y0) ), FixedValue(value=0,faces=mesh.getFacesRight() & (Y<y0) ), FixedValue(value=2,faces=mesh.getFacesTop() ), ) potential.equation.solve(var=potential, boundaryConditions=bcs) # The default visualization method viewer = viewers.Viewer(vars=(potential,)) viewer.plot("/path/to/output.vtk") # The follow evaluation of the solution is only # possible with "common" meshes. result = pylab.array(potential) result = result.reshape((nx,nx)) xx,yy = pylab.array(x), pylab.array(y) xx,yy = xx.reshape((nx,nx)), yy.reshape((nx,nx)) pylab.title("Capacitor: +5/+0 V, Top Plate: +2 V") pylab.contourf(xx,yy,result,levels=pylab.linspace(result.min(),result.max(),64)) pylab.colorbar() pylab.show()