examples.diffusion.explicit package¶

Submodules¶

examples.diffusion.explicit.mesh1D module¶

This input file again solves a 1D diffusion problem as in examples.diffusion.steadyState.mesh1D, the difference being that this transient example is solved explicitly.

We create a 1D mesh:

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

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

and we initialize a CellVariable to initialValue:

>>> valueLeft = 0.
>>> initialValue = 1.
>>> var = CellVariable(
...     name = "concentration",
...     mesh = mesh,
...     value = initialValue)

The transient diffusion equation

$\frac{\partial \phi}{\partial t} = \nabla \cdot (D \nabla \phi)$

is represented by a TransientTerm and an ExplicitDiffusionTerm.

We take the diffusion coefficient $D = 1$

>>> diffusionCoeff = 1.

We build the equation:

>>> eq = TransientTerm() == ExplicitDiffusionTerm(coeff = diffusionCoeff)

and the boundary conditions:

>>> var.constrain(valueLeft, mesh.facesLeft)

In this case, many steps have to be taken to reach equilibrium. A loop is required to execute the necessary time steps:

>>> timeStepDuration = 0.1
>>> steps = 100
>>> from builtins import range
>>> for step in range(steps):
...     eq.solve(var=var, dt=timeStepDuration)

The analytical solution for this transient diffusion problem is given by $\phi = \erf(x/2\sqrt{D t})$ .

The result is tested against the expected profile:

>>> Lx = nx * dx
>>> x = mesh.cellCenters[0]
>>> t = timeStepDuration * steps
>>> epsi = x / numerix.sqrt(t * diffusionCoeff)
>>> from scipy.special import erf 
>>> analyticalArray = erf(epsi/2) 
>>> print(var.allclose(analyticalArray, atol = 2e-3)) 
1

If the problem is run interactively, we can view the result:

>>> if __name__ == '__main__':
...     viewer = Viewer(vars = (var,))
...     viewer.plot()

examples.diffusion.explicit.mixedelement module¶

This input file again solves an explicit 1D diffusion problem as in ./examples/diffusion/mesh1D.py but on a mesh with both square and triangular elements. The term used is the ExplicitDiffusionTerm. In this case many steps have to be taken to reach equilibrium. The timeStepDuration parameter specifies the size of each time step and steps is the number of time steps.

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

>>> dx = 1.
>>> dy = 1.
>>> nx = 10

>>> valueLeft = 0.
>>> valueRight = 1.
>>> D = 1.
>>> timeStepDuration = 0.01 * dx**2 / (2 * D)
>>> if __name__ == '__main__':
...     steps = 10000
... else:
...     steps = 10

>>> gridMesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=2)
>>> triMesh = Tri2D(dx=dx, dy=dy, nx=nx, ny=1) + ((dx*nx,), (0,))
>>> otherGridMesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=1) + ((dx*nx,), (1,))
>>> bigMesh = gridMesh + triMesh + otherGridMesh

>>> L = dx * nx * 2

>>> var = CellVariable(name="concentration",
...                    mesh=bigMesh,
...                    value=valueLeft)

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

>>> var.constrain(valueLeft, where=bigMesh.facesLeft)
>>> var.constrain(valueRight, where=bigMesh.facesRight)

In a semi-infinite domain, the analytical solution for this transient diffusion problem is given by $\phi = 1 - \erf((L - x)/2\sqrt{D t})$ , which is a reasonable approximation at early times. At late times, the solution is just a straight line. If the SciPy library is available, the result is tested against the expected profile:

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

>>> if __name__ == '__main__':
...     varAnalytical = valueLeft + (valueRight - valueLeft) * x / L
...     atol = 0.2
... else:
...     try:
...         from scipy.special import erf 
...         varAnalytical = 1 - erf((L - x) / (2 * numerix.sqrt(D * t))) 
...         atol = 0.03
...     except ImportError:
...         print("The SciPy library is not available to test the solution to ...           the transient diffusion equation")

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=var)

>>> from builtins import range
>>> for step in range(steps):
...     eqn.solve(var, dt=timeStepDuration)
...     if (step % 100) == 0 and (__name__ == '__main__'):
...         viewer.plot()

We check the answer against the analytical result

>>> print(var.allclose(varAnalytical, atol=atol)) 
1

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input('finished')

examples.diffusion.explicit.test module¶

examples.diffusion.explicit.tri2D module¶

This input file again solves a 1D diffusion problem as in ./examples/diffusion/steadyState/mesh1D.py. The difference in this example is that the solution method is explicit. The equation used is the ExplicitDiffusionEquation. In this case many steps have to be taken to reach equilibrium. The timeStepDuration parameter specifies the size of each time step and steps is the number of time steps.

>>> dx = 1.
>>> dy = 1.
>>> nx = 10
>>> ny = 1
>>> valueLeft = 0.
>>> valueRight = 1.
>>> timeStepDuration = 0.02
>>> steps = 10

A loop is required to execute the necessary time steps:

>>> from builtins import range
>>> for step in range(steps):
...     eq.solve(var, solver=solver, dt=timeStepDuration)

The result is again tested in the same way:

>>> Lx = nx * dx
>>> x = mesh.cellCenters[0]
>>> print(var.allclose(answer, rtol = 1e-8))
1

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