examples.diffusion package¶

Subpackages¶

Submodules¶

examples.diffusion.anisotropy module¶

Solve the diffusion equation with an anisotropic diffusion coefficient.

We wish to solve the problem

$\frac{\partial \phi}{\partial t} = \partial_j \Gamma_{ij} \partial_i \phi$

on a circular domain centered at $(0, 0)$ . We can choose an anisotropy ratio of 5 such that

$\Gamma' = \begin{pmatrix} 0.2 & 0 \\ 0 & 1 \end{pmatrix}$

A new matrix is formed by rotating $\Gamma'$ such that

$R = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix}$

and

$\Gamma = R \Gamma' R^T$

In the case of a point source at $(0, 0)$ a reference solution is given by,

$\phi \left( X, Y, t \right) = Q \frac{ \exp \left( -\frac{1}{4 t} \left( \frac{ X^2 }{ \Gamma'_{00}} + \frac{ Y^2 }{ \Gamma'_{11}} \right) \right) }{ 4 \pi t \sqrt{\Gamma'_{00} \Gamma'_{11}} }$

where $\left(X, Y \right)^T = R \left(x, y \right)^T$ and $Q$ is the initial mass.

>>> from fipy import CellVariable, Gmsh2D, Viewer, TransientTerm, DiffusionTermCorrection
>>> from fipy.tools import serialComm, numerix

Import a mesh previously created using Gmsh.

>>> import os
>>> mesh = Gmsh2D(os.path.splitext(__file__)[0] + '.msh', communicator=serialComm) 

Set the center-most cell to have a value.

>>> var = CellVariable(mesh=mesh, hasOld=1) 
>>> x, y = mesh.cellCenters 
>>> var[numerix.argmin(x**2 + y**2)] = 1. 

Choose an orientation for the anisotropy.

>>> theta = numerix.pi / 4.
>>> rotationMatrix = numerix.array(((numerix.cos(theta), numerix.sin(theta)), \
...                                 (-numerix.sin(theta), numerix.cos(theta))))
>>> gamma_prime = numerix.array(((0.2, 0.), (0., 1.)))
>>> DOT = numerix.NUMERIX.dot
>>> gamma = DOT(DOT(rotationMatrix, gamma_prime), numerix.transpose(rotationMatrix))

Make the equation, viewer and solve.

>>> eqn = TransientTerm() == DiffusionTermCorrection((gamma,))

>>> if __name__ == '__main__':
...     viewer = Viewer(var, datamin=0.0, datamax=0.001)

>>> mass = float(var.cellVolumeAverage * numerix.sum(mesh.cellVolumes)) 
>>> time = 0
>>> dt=0.00025

>>> from builtins import range
>>> for i in range(20):
...     var.updateOld() 
...     res = 1.
...
...     while res > 1e-2:
...         res = eqn.sweep(var, dt=dt) 
...
...     if __name__ == '__main__':
...         viewer.plot()
...     time += dt

Compare with the analytical solution (within 5% accuracy).

>>> X, Y = numerix.dot(mesh.cellCenters, CellVariable(mesh=mesh, rank=2, value=rotationMatrix)) 
>>> solution = mass * numerix.exp(-(X**2 / gamma_prime[0][0] + Y**2 / gamma_prime[1][1]) / (4 * time)) / (4 * numerix.pi * time * numerix.sqrt(gamma_prime[0][0] * gamma_prime[1][1])) 
>>> print(max(abs((var - solution) / max(solution))) < 0.08) 
True

examples.diffusion.circle module¶

Solve the diffusion equation in a circular domain meshed with triangles.

This example demonstrates how to solve a simple diffusion problem on a non-standard mesh with varying boundary conditions. The Gmsh package is used to create the mesh. Firstly, define some parameters for the creation of the mesh,

>>> cellSize = 0.05
>>> radius = 1.

The cellSize is the preferred edge length of each mesh element and the radius is the radius of the circular mesh domain. In the following code section a file is created with the geometry that describes the mesh. For details of how to write such geometry files for Gmsh, see the gmsh manual.

The mesh created by Gmsh is then imported into FiPy using the Gmsh2D object.

>>> from fipy import CellVariable, Gmsh2D, TransientTerm, DiffusionTerm, Viewer
>>> from fipy.tools import numerix

>>> mesh = Gmsh2D('''
...               cellSize = %(cellSize)g;
...               radius = %(radius)g;
...               Point(1) = {0, 0, 0, cellSize};
...               Point(2) = {-radius, 0, 0, cellSize};
...               Point(3) = {0, radius, 0, cellSize};
...               Point(4) = {radius, 0, 0, cellSize};
...               Point(5) = {0, -radius, 0, cellSize};
...               Circle(6) = {2, 1, 3};
...               Circle(7) = {3, 1, 4};
...               Circle(8) = {4, 1, 5};
...               Circle(9) = {5, 1, 2};
...               Line Loop(10) = {6, 7, 8, 9};
...               Plane Surface(11) = {10};
...               ''' % locals()) 

Using this mesh, we can construct a solution variable

>>> phi = CellVariable(name = "solution variable",
...                    mesh = mesh,
...                    value = 0.) 

We can now create a Viewer to see the mesh

>>> viewer = None
>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer = Viewer(vars=phi, datamin=-1, datamax=1.)
...     viewer.plotMesh()

We set up a transient diffusion equation

>>> D = 1.
>>> eq = TransientTerm() == DiffusionTerm(coeff=D)

The following line extracts the $x$ coordinate values on the exterior faces. These are used as the boundary condition fixed values.

>>> X, Y = mesh.faceCenters 

>>> phi.constrain(X, mesh.exteriorFaces) 

We first step through the transient problem

>>> timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
>>> steps = 10
>>> from builtins import range
>>> for step in range(steps):
...     eq.solve(var=phi,
...              dt=timeStepDuration) 
...     if viewer is not None:
...         viewer.plot() 

evolution of diffusion problem on a circular mesh

If we wanted to plot or analyze the results of this calculation with another application, we could export tab-separated-values with

TSVViewer(vars=(phi, phi.grad)).plot(filename="myTSV.tsv")

The values are listed at the cell centers. Particularly for irregular meshes, no specific ordering should be relied upon. Vector quantities are listed in multiple columns, one for each mesh dimension.

This problem again has an analytical solution that depends on the error function, but it’s a bit more complicated due to the varying boundary conditions and the different horizontal diffusion length at different vertical positions

>>> x, y = mesh.cellCenters 
>>> t = timeStepDuration * steps

>>> phiAnalytical = CellVariable(name="analytical value",
...                              mesh=mesh) 

>>> x0 = radius * numerix.cos(numerix.arcsin(y)) 
>>> try:
...     from scipy.special import erf 
...     ## This function can sometimes throw nans on OS X
...     ## see http://projects.scipy.org/scipy/scipy/ticket/325
...     phiAnalytical.setValue(x0 * (erf((x0+x) / (2 * numerix.sqrt(D * t)))
...                                  - erf((x0-x) / (2 * numerix.sqrt(D * t))))) 
... except ImportError:
...     print("The SciPy library is not available to test the solution to \
... the transient diffusion equation")

>>> print(phi.allclose(phiAnalytical, atol = 7e-2)) 
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Transient diffusion. Press <return> to proceed...")

As in the earlier examples, we can also directly solve the steady-state diffusion problem.

>>> DiffusionTerm(coeff=D).solve(var=phi) 

The values at the elements should be equal to their x coordinate

>>> print(phi.allclose(x, atol = 0.03)) 
1

Display the results if run as a script.

>>> from fipy import input
>>> if viewer is not None:
...     viewer.plot()
...     input("Steady-state diffusion. Press <return> to proceed...")

steady-state solution to diffusion on a circular mesh

examples.diffusion.circleQuad module¶

Solve the diffusion equation in a circular domain meshed with quadrangles.

This example demonstrates how to solve a simple diffusion problem on a non-standard mesh with varying boundary conditions. The Gmsh package is used to create the mesh. Firstly, define some parameters for the creation of the mesh,

>>> cellSize = 0.05
>>> radius = 1.

The cellSize is the preferred edge length of each mesh element and the radius is the radius of the circular mesh domain. In the following code section a file is created with the geometry that describes the mesh. For details of how to write such geometry files for Gmsh, see the gmsh manual.

The mesh created by Gmsh is then imported into FiPy using the Gmsh2D object.

>>> from fipy import CellVariable, Gmsh2D, TransientTerm, DiffusionTerm, Viewer
>>> from fipy.tools import numerix

>>> mesh = Gmsh2D('''
...               cellSize = %(cellSize)g;
...               radius = %(radius)g;
...               Point(1) = {0, 0, 0, cellSize};
...               Point(2) = {-radius, 0, 0, cellSize};
...               Point(3) = {0, radius, 0, cellSize};
...               Point(4) = {radius, 0, 0, cellSize};
...               Point(5) = {0, -radius, 0, cellSize};
...               Circle(6) = {2, 1, 3};
...               Circle(7) = {3, 1, 4};
...               Circle(8) = {4, 1, 5};
...               Circle(9) = {5, 1, 2};
...               Line Loop(10) = {6, 7, 8, 9};
...               Plane Surface(11) = {10};
...               Recombine Surface{11};
...               ''' % locals()) 

Using this mesh, we can construct a solution variable

>>> phi = CellVariable(name = "solution variable",
...                    mesh = mesh,
...                    value = 0.) 

We can now create a Viewer to see the mesh

>>> viewer = None
>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer = Viewer(vars=phi, datamin=-1, datamax=1.)
...     viewer.plotMesh()
...     input("Irregular circular mesh. Press <return> to proceed...")

documentation/generated/examples/circleMesh.*

We set up a transient diffusion equation

>>> D = 1.
>>> eq = TransientTerm() == DiffusionTerm(coeff=D)

The following line extracts the $x$ coordinate values on the exterior faces. These are used as the boundary condition fixed values.

>>> X, Y = mesh.faceCenters 
>>> phi.constrain(X, mesh.exteriorFaces) 

We first step through the transient problem

>>> timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
>>> steps = 10
>>> from builtins import range
>>> for step in range(steps):
...     eq.solve(var=phi,
...              dt=timeStepDuration) 
...     if viewer is not None:
...         viewer.plot()

documentation/generated/examples/circleTransient.*

If we wanted to plot or analyze the results of this calculation with another application, we could export tab-separated-values with

TSVViewer(vars=(phi, phi.grad)).plot(filename="myTSV.tsv")

The values are listed at the Cell centers. Particularly for irregular meshes, no specific ordering should be relied upon. Vector quantities are listed in multiple columns, one for each mesh dimension.

This problem again has an analytical solution that depends on the error function, but it’s a bit more complicated due to the varying boundary conditions and the different horizontal diffusion length at different vertical positions

>>> x, y = mesh.cellCenters 
>>> t = timeStepDuration * steps

>>> phiAnalytical = CellVariable(name="analytical value",
...                              mesh=mesh) 

>>> x0 = radius * numerix.cos(numerix.arcsin(y)) 
>>> try:
...     from scipy.special import erf 
...     ## This function can sometimes throw nans on OS X
...     ## see http://projects.scipy.org/scipy/scipy/ticket/325
...     phiAnalytical.setValue(x0 * (erf((x0+x) / (2 * numerix.sqrt(D * t)))
...                                  - erf((x0-x) / (2 * numerix.sqrt(D * t))))) 
... except ImportError:
...     print("The SciPy library is not available to test the solution to \
... the transient diffusion equation")

>>> print(phi.allclose(phiAnalytical, atol = 7e-2)) 
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Transient diffusion. Press <return> to proceed...")

As in the earlier examples, we can also directly solve the steady-state diffusion problem.

>>> DiffusionTerm(coeff=D).solve(var=phi) 

The values at the elements should be equal to their x coordinate

>>> print(phi.allclose(x, atol = 0.035)) 
1

Display the results if run as a script.

>>> from fipy import input
>>> if viewer is not None:
...     viewer.plot()
...     input("Steady-state diffusion. Press <return> to proceed...")

documentation/generated/examples/circleSteadyState.*

examples.diffusion.coupled module¶

Solve the biharmonic equation as a coupled pair of diffusion equations.

FiPy has only first order time derivatives so equations such as the biharmonic wave equation written as

$\frac{\partial^4 v}{\partial x^4} + \frac{\partial^2 v}{\partial t^2} &= 0$

cannot be represented as a single equation. We need to decompose the biharmonic equation into two equations that are first order in time in the following way,

$\frac{\partial^2 v_0}{\partial x^2} + \frac{\partial v_1}{\partial t} &= 0 \\ \frac{\partial^2 v_1}{\partial x^2} - \frac{\partial v_0}{\partial t} &= 0$

Historically, FiPy required systems of coupled equations to be solved successively, “sweeping” the equations to convergence. As a practical example, we use the following system

$\frac{\partial v_0}{\partial t} &= 0.01 \nabla^2 v_0 - \nabla^2 v_1 \\ \frac{\partial v_1}{\partial t} &= \nabla^2 v_0 + 0.01 \nabla^2 v_1$

subject to the boundary conditions

$\begin{align*} v_0|_{x=0} &= 0 & v_0|_{x=1} &= 1 \\ v_1|_{x=0} &= 1 & v_1|_{x=1} &= 0 \end{align*}$

This system closely resembles the pure biharmonic equation, but has an additional diffusion contribution to improve numerical stability. The example system is solved with the following block of code using explicit coupling for the cross-coupled terms.

>>> from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer

>>> m = Grid1D(nx=100, Lx=1.)

>>> v0 = CellVariable(mesh=m, hasOld=True, value=0.5)
>>> v1 = CellVariable(mesh=m, hasOld=True, value=0.5)

>>> v0.constrain(0, m.facesLeft)
>>> v0.constrain(1, m.facesRight)

>>> v1.constrain(1, m.facesLeft)
>>> v1.constrain(0, m.facesRight)

>>> eq0 = TransientTerm() == DiffusionTerm(coeff=0.01) - v1.faceGrad.divergence
>>> eq1 = TransientTerm() == v0.faceGrad.divergence + DiffusionTerm(coeff=0.01)

>>> vi = Viewer((v0, v1))

>>> from builtins import range
>>> for t in range(100):
...     v0.updateOld()
...     v1.updateOld()
...     res0 = res1 = 1e100
...     while max(res0, res1) > 0.1:
...         res0 = eq0.sweep(var=v0, dt=1e-5)
...         res1 = eq1.sweep(var=v1, dt=1e-5)
...     if t % 10 == 0:
...         vi.plot()

The uncoupled method still works, but it can be advantageous to solve the two equations simultaneously. In this case, by coupling the equations, we can eliminate the explicit sources and dramatically increase the time steps:

>>> v0.value = 0.5
>>> v1.value = 0.5

>>> eqn0 = TransientTerm(var=v0) == DiffusionTerm(0.01, var=v0) - DiffusionTerm(1, var=v1)
>>> eqn1 = TransientTerm(var=v1) == DiffusionTerm(1, var=v0) + DiffusionTerm(0.01, var=v1)

>>> eqn = eqn0 & eqn1

>>> from builtins import range
>>> for t in range(1):
...     v0.updateOld()
...     v1.updateOld()
...     eqn.solve(dt=1.e-3)
...     vi.plot()

It is also possible to pose the same equations in vector form:

>>> v = CellVariable(mesh=m, hasOld=True, value=[[0.5], [0.5]], elementshape=(2,))

>>> v.constrain([[0], [1]], m.facesLeft)
>>> v.constrain([[1], [0]], m.facesRight)

>>> eqn = TransientTerm([[1, 0],
...                      [0, 1]]) == DiffusionTerm([[[0.01, -1],
...                                                  [1, 0.01]]])

>>> vi = Viewer((v[0], v[1]))

>>> from builtins import range
>>> for t in range(1):
...     v.updateOld()
...     eqn.solve(var=v, dt=1.e-3)
...     vi.plot()

Whether you pose your problem in coupled or vector form should be dictated by the underlying physics. If $v_0$ and $v_1$ represent the concentrations of two conserved species, then it is natural to write two separate governing equations and to couple them. If they represent two components of a vector field, then the vector formulation is obviously more natural. FiPy will solve the same matrix system either way.

examples.diffusion.electrostatics module¶

Solve the Poisson equation in one dimension.

The Poisson equation is a particular example of the steady-state diffusion equation. We examine a few cases in one dimension.

>>> from fipy import CellVariable, Grid1D, Viewer, DiffusionTerm

>>> nx = 200
>>> dx = 0.01
>>> L = nx * dx
>>> mesh = Grid1D(dx = dx, nx = nx)

Given the electrostatic potential $\phi$ ,

>>> potential = CellVariable(mesh=mesh, name='potential', value=0.)

the permittivity $\epsilon$ ,

>>> permittivity = 1

the concentration $C_j$ of the $j^\text{th}$ component with valence $z_j$ (we consider only a single component $C_\text{e}^{-}$ with valence with $z_{\text{e}^{-}} = -1$ )

>>> electrons = CellVariable(mesh=mesh, name='e-')
>>> electrons.valence = -1

and the charge density $\rho$ ,

>>> charge = electrons * electrons.valence
>>> charge.name = "charge"

The dimensionless Poisson equation is

$\nabla\cdot\left(\epsilon\nabla\phi\right) = -\rho = -\sum_{j=1}^n z_j C_j$

>>> potential.equation = (DiffusionTerm(coeff = permittivity)
...                       + charge == 0)

Because this equation admits an infinite number of potential profiles, we must constrain the solution by fixing the potential at one point:

>>> potential.constrain(0., mesh.facesLeft)

First, we obtain a uniform charge distribution by setting a uniform concentration of electrons $C_{\text{e}^{-}} = 1$ .

>>> electrons.setValue(1.)

and we solve for the electrostatic potential

>>> potential.equation.solve(var=potential)

This problem has the analytical solution

$\psi(x) = \frac{x^2}{2} - 2x$

>>> x = mesh.cellCenters[0]
>>> analytical = CellVariable(mesh=mesh, name="analytical solution",
...                           value=(x**2)/2 - 2*x)

which has been satisfactorily obtained

>>> print(potential.allclose(analytical, rtol = 2e-5, atol = 2e-5))
1

If we are running the example interactively, we view the result

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer = Viewer(vars=(charge, potential, analytical))
...     viewer.plot()
...     input("Press any key to continue...")

electrostatic potential generated by uniform charge

Next, we segregate all of the electrons to right side of the domain

$C_{\text{e}^{-}} = \begin{cases} 0& \text{for $x \le L/2$,} \\ 1& \text{for $x > L/2$.} \end{cases}$

>>> x = mesh.cellCenters[0]
>>> electrons.setValue(0.)
>>> electrons.setValue(1., where=x > L / 2.)

and again solve for the electrostatic potential

>>> potential.equation.solve(var=potential)

which now has the analytical solution

$\psi(x) = \begin{cases} -x& \text{for $x \le L/2$,} \\ \frac{(x-1)^2}{2} - x& \text{for $x > L/2$.} \end{cases}$

>>> analytical.setValue(-x)
>>> analytical.setValue(((x-1)**2)/2 - x, where=x > L/2)

>>> print(potential.allclose(analytical, rtol = 2e-5, atol = 2e-5))
1

and again view the result

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Press any key to continue...")

electrostatic potential generated by negative charge on right

Finally, we segregate all of the electrons to the left side of the domain

$C_{\text{e}^{-}} = \begin{cases} 1& \text{for $x \le L/2$,} \\ 0& \text{for $x > L/2$.} \end{cases}$

>>> electrons.setValue(1.)
>>> electrons.setValue(0., where=x > L / 2.)

and again solve for the electrostatic potential

>>> potential.equation.solve(var=potential)

which has the analytical solution

$\psi(x) = \begin{cases} \frac{x^2}{2} - x& \text{for $x \le L/2$,} \\ -\frac{1}{2}& \text{for $x > L/2$.} \end{cases}$

We again verify that the correct equilibrium is attained

>>> analytical.setValue((x**2)/2 - x)
>>> analytical.setValue(-0.5, where=x > L / 2)

>>> print(potential.allclose(analytical, rtol = 2e-5, atol = 2e-5))
1

and once again view the result

>>> if __name__ == '__main__':
...     viewer.plot()

electrostatic potential generated by negative charge on left

examples.diffusion.mesh1D module¶

Solve a one-dimensional diffusion equation under different conditions.

To run this example from the base FiPy directory, type:

$ python examples/diffusion/mesh1D.py

at the command line. Different stages of the example should be displayed, along with prompting messages in the terminal.

This example takes the user through assembling a simple problem with FiPy. It describes different approaches to a 1D diffusion problem with constant diffusivity and fixed value boundary conditions such that,

(1)¶ $\frac{\partial \phi}{\partial t} = D \nabla^2 \phi.$

The first step is to define a one dimensional domain with 50 solution points. The Grid1D object represents a linear structured grid. The parameter dx refers to the grid spacing (set to unity here).

>>> from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer
>>> from fipy.tools import numerix

>>> nx = 50
>>> dx = 1.
>>> mesh = Grid1D(nx=nx, dx=dx)

FiPy solves all equations at the centers of the cells of the mesh. We thus need a CellVariable object to hold the values of the solution, with the initial condition $\phi = 0$ at $t = 0$ ,

>>> phi = CellVariable(name="solution variable",
...                    mesh=mesh,
...                    value=0.)

We’ll let

>>> D = 1.

for now.

The boundary conditions

$\phi = \begin{cases} 0& \text{at $x = 1$,} \\ 1& \text{at $x = 0$.} \end{cases}$

are formed with a value

>>> valueLeft = 1
>>> valueRight = 0

and a set of faces over which they apply.

Note

Only faces around the exterior of the mesh can be used for boundary conditions.

For example, here the exterior faces on the left of the domain are extracted by mesh.facesLeft. The boundary conditions are applied using phi. constrain() with these faces and a value (valueLeft).

>>> phi.constrain(valueRight, mesh.facesRight)
>>> phi.constrain(valueLeft, mesh.facesLeft)

Note

If no boundary conditions are specified on exterior faces, the default boundary condition is no-flux, $\vec{n} \cdot (D \nabla \phi) \rvert_\text{someFaces} = 0$ .

If you have ever tried to numerically solve Eq. (1), you most likely attempted “explicit finite differencing” with code something like:

for step in range(steps):
    for j in range(cells):
        phi_new[j] = phi_old[j] \
          + (D * dt / dx**2) * (phi_old[j+1] - 2 * phi_old[j] + phi_old[j-1])
    time += dt

plus additional code for the boundary conditions. In FiPy, you would write

>>> eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D)

The largest stable timestep that can be taken for this explicit 1D diffusion problem is $\Delta t \le \Delta x^2 / (2 D)$ .

We limit our steps to 90% of that value for good measure

>>> timeStepDuration = 0.9 * dx**2 / (2 * D)
>>> steps = 100

If we’re running interactively, we’ll want to view the result, but not if this example is being run automatically as a test. We accomplish this by having Python check if this script is the “__main__” script, which will only be true if we explicitly launched it and not if it has been imported by another script such as the automatic tester. The factory function Viewer() returns a suitable viewer depending on available viewers and the dimension of the mesh.

>>> phiAnalytical = CellVariable(name="analytical value",
...                              mesh=mesh)

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=(phi, phiAnalytical),
...                     datamin=0., datamax=1.)
...     viewer.plot()

In a semi-infinite domain, the analytical solution for this transient diffusion problem is given by $\phi = 1 - \erf(x/2\sqrt{D t})$ . If the SciPy library is available, the result is tested against the expected profile:

>>> x = mesh.cellCenters[0]
>>> t = timeStepDuration * steps

>>> try:
...     from scipy.special import erf 
...     phiAnalytical.setValue(1 - erf(x / (2 * numerix.sqrt(D * t)))) 
... except ImportError:
...     print("The SciPy library is not available to test the solution to \
... the transient diffusion equation")

We then solve the equation by repeatedly looping in time:

>>> from builtins import range
>>> for step in range(steps):
...     eqX.solve(var=phi,
...               dt=timeStepDuration)
...     if __name__ == '__main__':
...         viewer.plot()

>>> print(phi.allclose(phiAnalytical, atol = 7e-4)) 
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Explicit transient diffusion. Press <return> to proceed...")

solution to diffusion problem evolved by explicit time steps

Although explicit finite differences are easy to program, we have just seen that this 1D transient diffusion problem is limited to taking rather small time steps. If, instead, we represent Eq. (1) as:

phi_new[j] = phi_old[j] \
  + (D * dt / dx**2) * (phi_new[j+1] - 2 * phi_new[j] + phi_new[j-1])

it is possible to take much larger time steps. Because phi_new appears on both the left and right sides of the equation, this form is called “implicit”. In general, the “implicit” representation is much more difficult to program than the “explicit” form that we just used, but in FiPy, all that is needed is to write

>>> eqI = TransientTerm() == DiffusionTerm(coeff=D)

reset the problem

>>> phi.setValue(valueRight)

and rerun with much larger time steps

>>> timeStepDuration *= 10
>>> steps //= 10
>>> from builtins import range
>>> for step in range(steps):
...     eqI.solve(var=phi,
...               dt=timeStepDuration)
...     if __name__ == '__main__':
...         viewer.plot()

>>> print(phi.allclose(phiAnalytical, atol = 2e-2)) 
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Implicit transient diffusion. Press <return> to proceed...")

solution to diffusion problem evolved by implicit time steps

Note that although much larger stable timesteps can be taken with this implicit version (there is, in fact, no limit to how large an implicit timestep you can take for this particular problem), the solution is less accurate. One way to achieve a compromise between stability and accuracy is with the Crank-Nicholson scheme, represented by:

phi_new[j] = phi_old[j] + (D * dt / (2 * dx**2)) * \
                ((phi_new[j+1] - 2 * phi_new[j] + phi_new[j-1])
                 + (phi_old[j+1] - 2 * phi_old[j] + phi_old[j-1]))

which is essentially an average of the explicit and implicit schemes from above. This can be rendered in FiPy as easily as

>>> eqCN = eqX + eqI

We again reset the problem

>>> phi.setValue(valueRight)

and apply the Crank-Nicholson scheme until the end, when we apply one step of the fully implicit scheme to drive down the error (see, e.g., section 19.2 of [22]).

>>> from builtins import range
>>> for step in range(steps - 1):
...     eqCN.solve(var=phi,
...                dt=timeStepDuration)
...     if __name__ == '__main__':
...         viewer.plot()
>>> eqI.solve(var=phi,
...           dt=timeStepDuration)
>>> if __name__ == '__main__':
...     viewer.plot()

>>> print(phi.allclose(phiAnalytical, atol = 3e-3)) 
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Crank-Nicholson transient diffusion. Press <return> to proceed...")

As mentioned above, there is no stable limit to how large a time step can be taken for the implicit diffusion problem. In fact, if the time evolution of the problem is not interesting, it is possible to eliminate the time step altogether by omitting the TransientTerm. The steady-state diffusion equation

$D \nabla^2 \phi = 0$

is represented in FiPy by

>>> DiffusionTerm(coeff=D).solve(var=phi)

>>> if __name__ == '__main__':
...     viewer.plot()

The analytical solution to the steady-state problem is no longer an error function, but simply a straight line, which we can confirm to a tolerance of $10^{-10}$ .

>>> L = nx * dx
>>> print(phi.allclose(valueLeft + (valueRight - valueLeft) * x / L,
...                    rtol = 1e-10, atol = 1e-10))
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Implicit steady-state diffusion. Press <return> to proceed...")

steady-state solution to diffusion problem

Often, boundary conditions may be functions of another variable in the system or of time.

For example, to have

$\phi = \begin{cases} (1 + \sin t) / 2 &\text{on $ x = 0 $} \\ 0 &\text{on $ x = L $} \\ \end{cases}$

we will need to declare time $t$ as a Variable

>>> time = Variable()

and then declare our boundary condition as a function of this Variable

>>> del phi.faceConstraints
>>> valueLeft = 0.5 * (1 + numerix.sin(time))
>>> phi.constrain(valueLeft, mesh.facesLeft)
>>> phi.constrain(0., mesh.facesRight)

>>> eqI = TransientTerm() == DiffusionTerm(coeff=D)

When we update time at each timestep, the left-hand boundary condition will automatically update,

>>> dt = .1
>>> while time() < 15:
...     time.setValue(time() + dt)
...     eqI.solve(var=phi, dt=dt)
...     if __name__ == '__main__':
...         viewer.plot()

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Time-dependent boundary condition. Press <return> to proceed...")

solution to diffusion problem with a time-dependent Dirichlet boundary condition

Many interesting problems do not have simple, uniform diffusivities. We consider a steady-state diffusion problem

$\nabla \cdot ( D \nabla \phi) = 0,$

with a spatially varying diffusion coefficient

$D = \begin{cases} 1& \text{for $ 0 < x < L / 4 $,} \\ 0.1& \text{for $ L / 4 \le x < 3 L / 4 $,} \\ 1& \text{for $ 3 L / 4 \le x < L $,} \end{cases}$

and with boundary conditions $\phi = 0$ at $x = 0$ and $D \frac{\partial \phi}{\partial x} = 1$ at $x = L$ , where $L$ is the length of the solution domain. Exact numerical answers to this problem are found when the mesh has cell centers that lie at $L / 4$ and $3 L / 4$ , or when the number of cells in the mesh $N_i$ satisfies $N_i = 4 i + 2$ , where $i$ is an integer. The mesh we’ve been using thus far is satisfactory, with $N_i = 50$ and $i = 12$ .

Because FiPy considers diffusion to be a flux from one cell to the next, through the intervening face, we must define the non-uniform diffusion coefficient on the mesh faces

>>> D = FaceVariable(mesh=mesh, value=1.0)
>>> X = mesh.faceCenters[0]
>>> D.setValue(0.1, where=(L / 4. <= X) & (X < 3. * L / 4.))

The boundary conditions are a fixed value of

>>> valueLeft = 0.

to the left and a fixed flux of

>>> fluxRight = 1.

to the right:

>>> phi = CellVariable(mesh=mesh)
>>> phi.faceGrad.constrain([fluxRight], mesh.facesRight)
>>> phi.constrain(valueLeft, mesh.facesLeft)

We re-initialize the solution variable

>>> phi.setValue(0)

and obtain the steady-state solution with one implicit solution step

>>> DiffusionTerm(coeff = D).solve(var=phi)

The analytical solution is simply

$\phi = \begin{cases} x & \text{for $ 0 < x < L/4 $,} \\ 10 x - 9L/4 & \text{for $ L/4 \le x < 3 L / 4 $,} \\ x + 18 L / 4 & \text{for $ 3 L / 4 \le x < L $,} \end{cases}$

or

>>> x = mesh.cellCenters[0]
>>> phiAnalytical.setValue(x)
>>> phiAnalytical.setValue(10 * x - 9. * L / 4.,
...                        where=(L / 4. <= x) & (x < 3. * L / 4.))
>>> phiAnalytical.setValue(x + 18. * L / 4.,
...                        where=3. * L / 4. <= x)
>>> print(phi.allclose(phiAnalytical, atol = 1e-8, rtol = 1e-8))
1

And finally, we can plot the result

>>> from fipy import input
>>> if __name__ == '__main__':
...     Viewer(vars=(phi, phiAnalytical)).plot()
...     input("Non-uniform steady-state diffusion. Press <return> to proceed...")

steady-state solution to diffusion problem with a non-uniform diffusivity

Note that for problems involving heat transfer and other similar conservation equations, it is important to ensure that we begin with the correct form of the equation. For example, for heat transfer with $\phi$ representing the temperature,

$\frac{\partial}{\partial t} \left(\rho \hat{C}_p \phi\right) = \nabla \cdot [ k \nabla \phi ].$

With constant and uniform density $\rho$ , heat capacity $\hat{C}_p$ and thermal conductivity $k$ , this is often written like Eq. (1), but replacing $D$ with $\alpha = \frac{k}{\rho \hat{C}_p}$ . However, when these parameters vary either in position or time, it is important to be careful with the form of the equation used. For example, if $k = 1$ and

$\rho \hat{C}_p = \begin{cases} 1& \text{for $ 0 < x < L / 4 $,} \\ 10& \text{for $ L / 4 \le x < 3 L / 4 $,} \\ 1& \text{for $ 3 L / 4 \le x < L $,} \end{cases},$

then we have

$\alpha = \begin{cases} 1& \text{for $ 0 < x < L / 4 $,} \\ 0.1& \text{for $ L / 4 \le x < 3 L / 4 $,} \\ 1& \text{for $ 3 L / 4 \le x < L $,} \end{cases}.$

However, using a DiffusionTerm with the same coefficient as that in the section above is incorrect, as the steady state governing equation reduces to $0 = \nabla^2\phi$ , which results in a linear profile in 1D, unlike that for the case above with spatially varying diffusivity. Similar care must be taken if there is time dependence in the parameters in transient problems.

We can illustrate the differences with an example. We define field variables for the correct and incorrect solution

>>> phiT = CellVariable(name="correct", mesh=mesh)
>>> phiF = CellVariable(name="incorrect", mesh=mesh)
>>> phiT.faceGrad.constrain([fluxRight], mesh.facesRight)
>>> phiF.faceGrad.constrain([fluxRight], mesh.facesRight)
>>> phiT.constrain(valueLeft, mesh.facesLeft)
>>> phiF.constrain(valueLeft, mesh.facesLeft)
>>> phiT.setValue(0)
>>> phiF.setValue(0)

The relevant parameters are

>>> k = 1.
>>> alpha_false = FaceVariable(mesh=mesh, value=1.0)
>>> X = mesh.faceCenters[0]
>>> alpha_false.setValue(0.1, where=(L / 4. <= X) & (X < 3. * L / 4.))
>>> eqF = 0 == DiffusionTerm(coeff=alpha_false)
>>> eqT = 0 == DiffusionTerm(coeff=k)
>>> eqF.solve(var=phiF)
>>> eqT.solve(var=phiT)

Comparing to the correct analytical solution, $\phi = x$

>>> x = mesh.cellCenters[0]
>>> phiAnalytical.setValue(x)
>>> print(phiT.allclose(phiAnalytical, atol = 1e-8, rtol = 1e-8)) 
1

and finally, plot

>>> from fipy import input
>>> if __name__ == '__main__':
...     Viewer(vars=(phiT, phiF)).plot()
...     input("Non-uniform thermal conductivity. Press <return> to proceed...")

representation of difference between non-uniform alpha and D

Often, the diffusivity is not only non-uniform, but also depends on the value of the variable, such that

(2)¶ $\frac{\partial \phi}{\partial t} = \nabla \cdot [ D(\phi) \nabla \phi].$

With such a non-linearity, it is generally necessary to “sweep” the solution to convergence. This means that each time step should be calculated over and over, using the result of the previous sweep to update the coefficients of the equation, without advancing in time. In FiPy, this is accomplished by creating a solution variable that explicitly retains its “old” value by specifying hasOld when you create it. The variable does not move forward in time until it is explicitly told to updateOld(). In order to compare the effects of different numbers of sweeps, let us create a list of variables: phi[0] will be the variable that is actually being solved and phi[1] through phi[4] will display the result of taking the corresponding number of sweeps (phi[1] being equivalent to not sweeping at all).

>>> valueLeft = 1.
>>> valueRight = 0.
>>> phi = [
...     CellVariable(name="solution variable",
...                  mesh=mesh,
...                  value=valueRight,
...                  hasOld=1),
...     CellVariable(name="1 sweep",
...                  mesh=mesh),
...     CellVariable(name="2 sweeps",
...                  mesh=mesh),
...     CellVariable(name="3 sweeps",
...                  mesh=mesh),
...     CellVariable(name="4 sweeps",
...                  mesh=mesh)
... ]

If, for example,

$D = D_0 (1 - \phi)$

we would simply write Eq. (2) as

>>> D0 = 1.
>>> eq = TransientTerm() == DiffusionTerm(coeff=D0 * (1 - phi[0]))

Note

Because of the non-linearity, the Crank-Nicholson scheme does not work for this problem.

We apply the same boundary conditions that we used for the uniform diffusivity cases

>>> phi[0].constrain(valueRight, mesh.facesRight)
>>> phi[0].constrain(valueLeft, mesh.facesLeft)

Although this problem does not have an exact transient solution, it can be solved in steady-state, with

$\phi(x) = 1 - \sqrt{\frac{x}{L}}$

>>> x = mesh.cellCenters[0]
>>> phiAnalytical.setValue(1. - numerix.sqrt(x/L))

We create a viewer to compare the different numbers of sweeps with the analytical solution from before.

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=phi + [phiAnalytical],
...                     datamin=0., datamax=1.)
...     viewer.plot()

As described above, an inner “sweep” loop is generally required for the solution of non-linear or multiple equation sets. Often a conditional is required to exit this “sweep” loop given some convergence criteria. Instead of using the solve() method equation, when sweeping, it is often useful to call sweep() instead. The sweep() method behaves the same way as solve(), but returns the residual that can then be used as part of the exit condition.

We now repeatedly run the problem with increasing numbers of sweeps.

>>> from fipy import input
>>> from builtins import range
>>> for sweeps in range(1, 5):
...     phi[0].setValue(valueRight)
...     for step in range(steps):
...         # only move forward in time once per time step
...         phi[0].updateOld()
...
...         # but "sweep" many times per time step
...         for sweep in range(sweeps):
...             res = eq.sweep(var=phi[0],
...                            dt=timeStepDuration)
...         if __name__ == '__main__':
...             viewer.plot()
...
...     # copy the final result into the appropriate display variable
...     phi[sweeps].setValue(phi[0])
...     if __name__ == '__main__':
...         viewer.plot()
...         input("Implicit variable diffusivity. %d sweep(s). \
... Residual = %f. Press <return> to proceed..." % (sweeps, (abs(res))))

As can be seen, sweeping does not dramatically change the result, but the “residual” of the equation (a measure of how accurately it has been solved) drops about an order of magnitude with each additional sweep.

Attention

Choosing an optimal balance between the number of time steps, the number of sweeps, the number of solver iterations, and the solver tolerance is more art than science and will require some experimentation on your part for each new problem.

Finally, we can increase the number of steps to approach equilibrium, or we can just solve for it directly

>>> eq = DiffusionTerm(coeff=D0 * (1 - phi[0]))

>>> phi[0].setValue(valueRight)
>>> res = 1e+10
>>> while res > 1e-6:
...     res = eq.sweep(var=phi[0],
...                    dt=timeStepDuration)

>>> print(phi[0].allclose(phiAnalytical, atol = 1e-1))
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Implicit variable diffusivity - steady-state. \
... Press <return> to proceed...")

solution to a diffusion problem a non-linear diffusivity

Fully implicit solutions are not without their pitfalls, particularly in steady state. Consider a localized block of material diffusing in a closed box.

>>> phi = CellVariable(mesh=mesh, name=r"$\phi$")

>>> phi.value = 0.
>>> phi.setValue(1., where=(x > L/2. - L/10.) & (x < L/2. + L/10.))
>>> if __name__ == '__main__':
...     viewer = Viewer(vars=phi, datamin=-0.1, datamax=1.1)

initial condition for no-flux boundary conditions

We assign no explicit boundary conditions, leaving the default no-flux boundary conditions, and solve

$\partial\phi/\partial t = \nabla\cdot(D\nabla\phi)$

>>> D = 1.
>>> eq = TransientTerm() == DiffusionTerm(D)

>>> dt = 10. * dx**2 / (2 * D)
>>> steps = 200

>>> from builtins import range
>>> for step in range(steps):
...     eq.solve(var=phi, dt=dt)
...     if __name__ == '__main__':
...         viewer.plot()
>>> from fipy import input
>>> if __name__ == '__main__':
...     input("No-flux - transient. \
... Press <return> to proceed...")

long-time solution for no-flux boundary conditions

and see that $\phi$ dissipates to the expected average value of 0.2 with reasonable accuracy.

>>> print(numerix.allclose(phi, 0.2, atol=1e-5))
True

If we reset the initial condition

>>> phi.value = 0.
>>> phi.setValue(1., where=(x > L/2. - L/10.) & (x < L/2. + L/10.))
>>> if __name__ == '__main__':
...     viewer.plot()

and solve the steady-state problem

>>> DiffusionTerm(coeff=D).solve(var=phi) 
>>> if __name__ == '__main__':
...     viewer.plot()
>>> from fipy import input
>>> if __name__ == '__main__':
...     input("No-flux - stead-state failure. \
... Press <return> to proceed...")

>>> print(numerix.allclose(phi, 0.0)) 
True

(failed) steady-state solution for no-flux boundary conditions

we find that the value is uniformly zero! What happened to our no-flux boundary conditions?

The problem is that in the implicit discretization of $\nabla\cdot(D\nabla\phi) = 0$ ,

$\begin{vmatrix} \frac{D}{{\Delta x}^2} & -\frac{D}{{\Delta x}^2} & & & & & \\\\[1em] \ddots & \ddots & \ddots & & & & \\[1em] & -\frac{D}{{\Delta x}^2} & \frac{2 D}{{\Delta x}^2} & -\frac{D}{{\Delta x}^2} & & & \\\\[1em] & & -\frac{D}{{\Delta x}^2} & \frac{2 D}{{\Delta x}^2} & -\frac{D}{{\Delta x}^2} & & \\\\[1em] & & & -\frac{D}{{\Delta x}^2} & \frac{2 D}{{\Delta x}^2} & -\frac{D}{{\Delta x}^2} & \\\\[1em] & & & & \ddots & \ddots & \ddots \\\\[1em] & & & & & -\frac{D}{{\Delta x}^2} & \frac{D}{{\Delta x}^2} \\\\[1em] \end{vmatrix} \begin{vmatrix} \phi^\text{new}_{0} \\\\[1em] \vdots \\[1em] \phi^\text{new}_{j-1} \\\\[1em] \phi^\text{new}_{j} \\\\[1em] \phi^\text{new}_{j+1} \\\\[1em] \vdots \\\\[1em] \phi^\text{new}_{N-1} \end{vmatrix} = \begin{vmatrix} 0 \\\\[1em] \vdots \\\\[1em] 0 \\\\[1em] 0 \\\\[1em] 0 \\\\[1em] \vdots \\\\[1em] 0 \end{vmatrix}$

the initial condition $\phi^\text{old}$ no longer appears and $\phi = 0$ is a perfectly legitimate solution to this matrix equation.

The solution is to run the transient problem and to take one enormous time step

>>> phi.value = 0.
>>> phi.setValue(1., where=(x > L/2. - L/10.) & (x < L/2. + L/10.))
>>> if __name__ == '__main__':
...     viewer.plot()

>>> (TransientTerm() == DiffusionTerm(D)).solve(var=phi, dt=1e6*dt)
>>> if __name__ == '__main__':
...     viewer.plot()
>>> from fipy import input
>>> if __name__ == '__main__':
...     input("No-flux - steady-state. \
... Press <return> to proceed...")

>>> print(numerix.allclose(phi, 0.2, atol=1e-5))
True

steady-state solution for no-flux boundary conditions

If this example had been written primarily as a script, instead of as documentation, we would delete every line that does not begin with either “>>>” or “...”, and then delete those prefixes from the remaining lines, leaving:

## This script was derived from
## 'examples/diffusion/mesh1D.py'

nx = 50
dx = 1.
mesh = Grid1D(nx = nx, dx = dx)
phi = CellVariable(name="solution variable",
                   mesh=mesh,
                   value=0)

eq = DiffusionTerm(coeff=D0 * (1 - phi[0]))
phi[0].setValue(valueRight)
res = 1e+10
while res > 1e-6:
    res = eq.sweep(var=phi[0],
                   dt=timeStepDuration)

print phi[0].allclose(phiAnalytical, atol = 1e-1)
# Expect:
# 1
#
if __name__ == '__main__':
    viewer.plot()
    input("Implicit variable diffusivity - steady-state. \
Press <return> to proceed...")

Your own scripts will tend to look like this, although you can always write them as doctest scripts if you choose. You can obtain a plain script like this from some …/example.py by typing:

$ python setup.py copy_script --From .../example.py --To myExample.py

at the command line.

Most of the FiPy examples will be a mixture of plain scripts and doctest documentation/tests.

examples.diffusion.mesh20x20 module¶

Solve a two-dimensional diffusion problem in a square domain.

This example solves a diffusion problem and demonstrates the use of applying boundary condition patches.

>>> from fipy import CellVariable, Grid2D, Viewer, TransientTerm, DiffusionTerm
>>> from fipy.tools import numerix

>>> nx = 20
>>> ny = nx
>>> dx = 1.
>>> dy = dx
>>> L = dx * nx
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

We create a CellVariable and initialize it to zero:

>>> phi = CellVariable(name = "solution variable",
...                    mesh = mesh,
...                    value = 0.)

and then create a diffusion equation. This is solved by default with an iterative conjugate gradient solver.

>>> D = 1.
>>> eq = TransientTerm() == DiffusionTerm(coeff=D)

We apply Dirichlet boundary conditions

>>> valueTopLeft = 0
>>> valueBottomRight = 1

to the top-left and bottom-right corners. Neumann boundary conditions are automatically applied to the top-right and bottom-left corners.

>>> X, Y = mesh.faceCenters
>>> facesTopLeft = ((mesh.facesLeft & (Y > L / 2))
...                 | (mesh.facesTop & (X < L / 2)))
>>> facesBottomRight = ((mesh.facesRight & (Y < L / 2))
...                     | (mesh.facesBottom & (X > L / 2)))

>>> phi.constrain(valueTopLeft, facesTopLeft)
>>> phi.constrain(valueBottomRight, facesBottomRight)

We create a viewer to see the results

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=phi, datamin=0., datamax=1.)
...     viewer.plot()

and solve the equation by repeatedly looping in time:

>>> timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
>>> steps = 10
>>> from builtins import range
>>> for step in range(steps):
...     eq.solve(var=phi,
...              dt=timeStepDuration)
...     if __name__ == '__main__':
...         viewer.plot()

solution to diffusion problem on a 2D domain with some Dirichlet boundaries

We can test the value of the bottom-right corner cell.

>>> print(numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2))
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Implicit transient diffusion. Press <return> to proceed...")

We can also solve the steady-state problem directly

>>> DiffusionTerm().solve(var=phi)
>>> if __name__ == '__main__':
...     viewer.plot()

stead-state solution to diffusion problem on a 2D domain with some Dirichlet boundaries

and test the value of the bottom-right corner cell.

>>> print(numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2))
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Implicit steady-state diffusion. Press <return> to proceed...")

examples.diffusion.mesh20x20Coupled module¶

Solve a coupled set of diffusion equations in two dimensions.

This example solves a diffusion problem and demonstrates the use of applying boundary condition patches.

>>> from fipy import CellVariable, Grid2D, Viewer, TransientTerm, DiffusionTerm
>>> from fipy.tools import numerix

>>> nx = 20
>>> ny = nx
>>> dx = 1.
>>> dy = dx
>>> L = dx * nx
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

We create a CellVariable and initialize it to zero:

>>> phi = CellVariable(name = "solution variable",
...                    mesh = mesh,
...                    value = 0.)

and then create a diffusion equation. This is solved by default with an iterative conjugate gradient solver.

>>> D = 1.
>>> eq = TransientTerm(var=phi) == DiffusionTerm(coeff=D, var=phi)

We apply Dirichlet boundary conditions

>>> valueTopLeft = 0
>>> valueBottomRight = 1

to the top-left and bottom-right corners. Neumann boundary conditions are automatically applied to the top-right and bottom-left corners.

>>> x, y = mesh.faceCenters
>>> facesTopLeft = ((mesh.facesLeft & (y > L / 2))
...                 | (mesh.facesTop & (x < L / 2)))
>>> facesBottomRight = ((mesh.facesRight & (y < L / 2))
...                     | (mesh.facesBottom & (x > L / 2)))

>>> phi.constrain(valueTopLeft, facesTopLeft)
>>> phi.constrain(valueBottomRight, facesBottomRight)

We create a viewer to see the results

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=phi, datamin=0., datamax=1.)
...     viewer.plot()

and solve the equation by repeatedly looping in time:

>>> timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
>>> steps = 10
>>> from builtins import range
>>> for step in range(steps):
...     eq.solve(dt=timeStepDuration)
...     if __name__ == '__main__':
...         viewer.plot()

documentation/generated/examples/mesh20x20transient.*

We can test the value of the bottom-right corner cell.

>>> print(numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2))
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Implicit transient diffusion. Press <return> to proceed...")

We can also solve the steady-state problem directly

>>> DiffusionTerm(var=phi).solve()
>>> if __name__ == '__main__':
...     viewer.plot()

documentation/generated/examples/mesh20x20steadyState.*

and test the value of the bottom-right corner cell.

>>> print(numerix.allclose(phi(((L,), (0,))), valueBottomRight, atol = 1e-2))
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Implicit steady-state diffusion. Press <return> to proceed...")

examples.diffusion.test module¶

Run all the test cases in examples/diffusion/

examples.diffusion.variable module¶

This example is a 1D steady state diffusion test case as in ./examples/diffusion/variable/mesh2x1/input.py with then number of cells set to nx = 10.

A simple analytical answer can be used to test the result:

>>> DiffusionTerm(coeff = diffCoeff).solve(var)
>>> if __name__ == "__main__":
...     viewer = Viewer(vars = var)
...     viewer.plot()
>>> x = mesh.cellCenters[0]
>>> values = numerix.where(x < 3. * L / 4., 10 * x - 9. * L / 4., x + 18. * L / 4.)
>>> values = numerix.where(x < L / 4., x, values)
>>> print(var.allclose(values, atol = 1e-8, rtol = 1e-8))
1

Last updated on Jun 27, 2023. Created using Sphinx 6.2.1.