examples.diffusion.mesh1D¶
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)¶
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 at ,
>>> phi = CellVariable(name="solution variable",
... mesh=mesh,
... value=0.)
We’ll let
>>> D = 1.
for now.
The boundary conditions
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, .
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 .
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 . 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...")
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...")
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
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 .
>>> 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...")
Often, boundary conditions may be functions of another variable in the system or of time.
For example, to have
we will need to declare time 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...")
Many interesting problems do not have simple, uniform diffusivities. We consider a steady-state diffusion problem
with a spatially varying diffusion coefficient
and with boundary conditions at and at , where is the length of the solution domain. Exact numerical answers to this problem are found when the mesh has cell centers that lie at and , or when the number of cells in the mesh satisfies , where is an integer. The mesh we’ve been using thus far is satisfactory, with and .
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
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...")
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 representing the temperature,
With constant and uniform density , heat capacity and thermal conductivity , this is often written like Eq. (1), but replacing with . 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 and
then we have
However, using a DiffusionTerm
with the same coefficient as that in the
section above is incorrect, as the steady state governing equation reduces to
, 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,
>>> 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...")
Often, the diffusivity is not only non-uniform, but also depends on the value of the variable, such that
(2)¶
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,
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
>>> 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...")
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)
We assign no explicit boundary conditions, leaving the default no-flux boundary conditions, and solve
>>> 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...")
and see that 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
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 ,
the initial condition no longer appears and 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
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.