examples.diffusion.explicit.mesh1DΒΆ

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()
Last updated on Jun 26, 2024. Created using Sphinx 7.1.2.