examples.elphf.diffusion package¶

Submodules¶

examples.elphf.diffusion.mesh1D module¶

A simple 1D example to test the setup of the multicomponent diffusion equations. The diffusion equation for each species in single-phase multicomponent system can be expressed as

$\frac{\partial C_j}{\partial t} = D_{jj}\nabla^2 C_j + D_{j}\nabla\cdot \frac{C_j}{1 - \sum_{\substack{k=2\\ k \neq j}}^{n-1} C_k} \sum_{\substack{i=2\\ i \neq j}}^{n-1} \nabla C_i$

where $C_j$ is the concentration of the $j^\text{th}$ species, $t$ is time, $D_{jj}$ is the self-diffusion coefficient of the $j^\text{th}$ species, and $\sum_{\substack{i=2\\ i \neq j}}^{n-1}$ represents the summation over all substitutional species in the system, excluding the solvent and the component of interest.

We solve the problem on a 1D mesh

>>> nx = 400
>>> dx = 0.01
>>> L = nx * dx

>>> from fipy import CellVariable, FaceVariable, Grid1D, TransientTerm, DiffusionTerm, PowerLawConvectionTerm, DefaultAsymmetricSolver, Viewer
>>> mesh = Grid1D(dx = dx, nx = nx)

One component in this ternary system will be designated the “solvent”

>>> class ComponentVariable(CellVariable):
...     def __init__(self, mesh, value = 0., name = '',
...                  standardPotential = 0., barrier = 0.,
...                  diffusivity = None, valence = 0, equation = None):
...         CellVariable.__init__(self, mesh = mesh, value = value,
...                               name = name)
...         self.standardPotential = standardPotential
...         self.barrier = barrier
...         self.diffusivity = diffusivity
...         self.valence = valence
...         self.equation = equation
...
...     def copy(self):
...         return self.__class__(mesh = self.mesh,
...                               value = self.value,
...                               name = self.name,
...                               standardPotential =
...                                   self.standardPotential,
...                               barrier = self.barrier,
...                               diffusivity = self.diffusivity,
...                               valence = self.valence,
...                               equation = self.equation)

>>> solvent = ComponentVariable(mesh = mesh, name = 'Cn', value = 1.)

We can create an arbitrary number of components, simply by providing a tuple or list of components

>>> substitutionals = [
...     ComponentVariable(mesh = mesh, name = 'C1', diffusivity = 1.,
...                       standardPotential = 1., barrier = 1.),
...     ComponentVariable(mesh = mesh, name = 'C2', diffusivity = 1.,
...                       standardPotential = 1., barrier = 1.),
...     ]

>>> interstitials = []

>>> for component in substitutionals:
...     solvent -= component

We separate the solution domain into two different concentration regimes

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

We create one diffusion equation for each substitutional component

>>> for Cj in substitutionals:
...     CkSum = ComponentVariable(mesh = mesh, value = 0.)
...     CkFaceSum = FaceVariable(mesh = mesh, value = 0.)
...     for Ck in [Ck for Ck in substitutionals if Ck is not Cj]:
...         CkSum += Ck
...         CkFaceSum += Ck.harmonicFaceValue
...
...     convectionCoeff = CkSum.faceGrad \
...                       * (Cj.diffusivity / (1. - CkFaceSum))
...
...     Cj.equation = (TransientTerm()
...                    == DiffusionTerm(coeff=Cj.diffusivity)
...                    + PowerLawConvectionTerm(coeff=convectionCoeff))
...     Cj.solver = DefaultAsymmetricSolver(precon=None, iterations=3200)

If we are running interactively, we create a viewer to see the results

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

Now, we iterate the problem to equilibrium, plotting as we go

>>> from builtins import range
>>> for i in range(40):
...     for Cj in substitutionals:
...         Cj.equation.solve(var=Cj,
...                           dt=10000.,
...                           solver=Cj.solver)
...     if __name__ == '__main__':
...         viewer.plot()

Since there is nothing to maintain the concentration separation in this problem, we verify that the concentrations have become uniform

Note

Between petsc=3.13.2=h82b89f7_0 and petsc=3.13.4=h82b89f7_0, PETSc ceased achieving 1e-7 tolerance when solving on 2 processors on Linux. Solving on macOS is OK. Solving on 1, 3, or 4 processors is OK.

>>> print(substitutionals[0].allclose(0.45, rtol = 2e-7, atol = 2e-7))
True
>>> print(substitutionals[1].allclose(0.45, rtol = 2e-7, atol = 2e-7))
True

examples.elphf.diffusion.mesh1Ddimensional module¶

In this example, we present the same three-component diffusion problem introduced in examples/elphf/diffusion/mesh1D.py but we demonstrate FiPy’s facility to use dimensional quantities.

>>> import warnings
>>> warnings.warn("\n\n\tSupport for physical dimensions is incomplete.\n\tIt is not possible to solve dimensional equations.\n")

>>> from fipy import CellVariable, FaceVariable, PhysicalField, Grid1D, TransientTerm, DiffusionTerm, PowerLawConvectionTerm, LinearLUSolver, Viewer
>>> from fipy.tools import numerix

We solve the problem on a 40 mm long 1D mesh

>>> nx = 40
>>> dx = PhysicalField(1., "mm")
>>> L = nx * dx
>>> mesh = Grid1D(dx = dx, nx = nx)

Again, one component in this ternary system will be designated the “solvent”

>>> class ComponentVariable(CellVariable):
...     def __init__(self, mesh, value = 0., name = '',
...                  standardPotential = 0., barrier = 0.,
...                  diffusivity = None, valence = 0, equation = None):
...         CellVariable.__init__(self, mesh = mesh, value = value,
...                               name = name)
...         self.standardPotential = Variable(standardPotential)
...         self.barrier = Variable(barrier)
...         self.diffusivity = Variable(diffusivity)
...         self.valence = valence
...         self.equation = equation
...
...     def copy(self):
...         return self.__class__(mesh = self.mesh,
...                               value = self.value,
...                               name = self.name,
...                               standardPotential =
...                                   self.standardPotential,
...                               barrier = self.barrier,
...                               diffusivity = self.diffusivity,
...                               valence = self.valence,
...                               equation = self.equation)

>>> solvent = ComponentVariable(mesh = mesh, name = 'Cn', value = "1 mol/m**3")

We can create an arbitrary number of components, simply by providing a Tuple or list of components

>>> substitutionals = [
...     ComponentVariable(mesh = mesh, name = 'C1', diffusivity = "1e-9 m**2/s",
...                       standardPotential = 1., barrier = 1., value = "0.3 mol/m**3"),
...     ComponentVariable(mesh = mesh, name = 'C2', diffusivity = "1e-9 m**2/s",
...                       standardPotential = 1., barrier = 1., value = "0.6 mol/m**3"),
...     ]

>>> interstitials = []

>>> for component in substitutionals:
...     solvent -= component

We separate the solution domain into two different concentration regimes

>>> x = mesh.cellCenters[0]
>>> substitutionals[0].setValue("0.3 mol/m**3")
>>> substitutionals[0].setValue("0.6 mol/m**3", where=x > L / 2)
>>> substitutionals[1].setValue("0.6 mol/m**3")
>>> substitutionals[1].setValue("0.3 mol/m**3", where=x > L / 2)

We create one diffusion equation for each substitutional component

>>> for Cj in substitutionals:
...     CkSum = ComponentVariable(mesh = mesh, value = 0.)
...     CkFaceSum = FaceVariable(mesh = mesh, value = 0.)
...     for Ck in [Ck for Ck in substitutionals if Ck is not Cj]:
...         CkSum += Ck
...         CkFaceSum += Ck.harmonicFaceValue
...
...     convectionCoeff = CkSum.faceGrad \
...                       * (Cj.diffusivity / (1. - CkFaceSum))
...
...     Cj.equation = (TransientTerm()
...                    == DiffusionTerm(coeff=Cj.diffusivity)
...                    + PowerLawConvectionTerm(coeff = convectionCoeff))

If we are running interactively, we create a viewer to see the results

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

Now, we iterate the problem to equilibrium, plotting as we go

>>> solver = LinearLUSolver()

>>> from builtins import range
>>> for i in range(40):
...     for Cj in substitutionals:
...         Cj.updateOld()
...     for Cj in substitutionals:
...         Cj.equation.solve(var = Cj,
...                           dt = "1000 s",
...                           solver = solver)
...     if __name__ == '__main__':
...         viewer.plot()

Since there is nothing to maintain the concentration separation in this problem, we verify that the concentrations have become uniform

>>> print(substitutionals[0].scaled.allclose("0.45 mol/m**3",
...     atol = "1e-7 mol/m**3", rtol = 1e-7))
1
>>> print(substitutionals[1].scaled.allclose("0.45 mol/m**3",
...     atol = "1e-7 mol/m**3", rtol = 1e-7))
1

Note

The absolute tolerance atol must be in units compatible with the value to be checked, but the relative tolerance rtol is dimensionless.

examples.elphf.diffusion.mesh2D module¶

The same three-component diffusion problem as introduced in examples/elphf/diffusion/mesh1D.py but in 2D:

>>> from fipy import CellVariable, FaceVariable, Grid2D, TransientTerm, DiffusionTerm, PowerLawConvectionTerm, Viewer

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

One component in this ternary system will be designated the “solvent”

>>> class ComponentVariable(CellVariable):
...     def __init__(self, mesh, value = 0., name = '', standardPotential = 0.,
...                  barrier = 0., diffusivity = None, valence = 0, equation = None):
...         CellVariable.__init__(self, mesh = mesh, value = value, name = name)
...         self.standardPotential = standardPotential
...         self.barrier = barrier
...         self.diffusivity = diffusivity
...         self.valence = valence
...         self.equation = equation
...
...     def copy(self):
...         return self.__class__(mesh = self.mesh, value = self.value,
...                               name = self.name,
...                               standardPotential = self.standardPotential,
...                               barrier = self.barrier,
...                               diffusivity = self.diffusivity,
...                               valence = self.valence,
...                               equation = self.equation)

>>> solvent = ComponentVariable(mesh = mesh, name = 'Cn', value = 1.)

We can create an arbitrary number of components, simply by providing a Tuple or list of components

>>> substitutionals = [
...     ComponentVariable(mesh = mesh, name = 'C1', diffusivity = 1.,
...                       standardPotential = 1., barrier = 1.),
...     ComponentVariable(mesh = mesh, name = 'C2', diffusivity = 1.,
...                       standardPotential = 1., barrier = 1.),
...     ]

>>> interstitials = []

>>> for component in substitutionals:
...     solvent -= component

Although we are not interested in them for this problem, we create one field to represent the “phase” (1 everywhere)

>>> phase = CellVariable(mesh = mesh, name = 'xi', value = 1.)

and one field to represent the electrostatic potential (0 everywhere)

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

Although it is constant in this problem, in later problems we will need the following functions of the phase field

>>> def pPrime(xi):
...     return 30. * (xi * (1 - xi))**2

>>> def gPrime(xi):
...     return 2 * xi * (1 - xi) * (1 - 2 * xi)

We separate the solution domain into two different concentration regimes

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

We create one diffusion equation for each substitutional component

>>> for Cj in substitutionals:
...     CkSum = ComponentVariable(mesh = mesh, value = 0.)
...     CkFaceSum = FaceVariable(mesh = mesh, value = 0.)
...     for Ck in [Ck for Ck in substitutionals if Ck is not Cj]:
...         CkSum += Ck
...         CkFaceSum += Ck.harmonicFaceValue
...
...     counterDiffusion = CkSum.faceGrad
...     phaseTransformation = \
...         (pPrime(phase.harmonicFaceValue) * Cj.standardPotential \
...         + gPrime(phase.harmonicFaceValue) * Cj.barrier) \
...             * phase.faceGrad
...     electromigration = Cj.valence * potential.faceGrad
...     convectionCoeff = counterDiffusion \
...         + solvent.harmonicFaceValue \
...             * (phaseTransformation + electromigration)
...     convectionCoeff *= (Cj.diffusivity / (1. - CkFaceSum))
...
...     Cj.equation = (TransientTerm()
...                    == DiffusionTerm(coeff=Cj.diffusivity)
...                    + PowerLawConvectionTerm(coeff=convectionCoeff))

If we are running interactively, we create a viewer to see the results

>>> if __name__ == '__main__':
...     viewers = [Viewer(vars=field, datamin=0, datamax=1)
...                for field in [solvent] + substitutionals]
...     for viewer in viewers:
...         viewer.plot()
...     steps = 40
...     tol = 1e-7
... else:
...     steps = 20
...     tol = 1e-4

Now, we iterate the problem to equilibrium, plotting as we go

>>> from builtins import range
>>> for i in range(steps):
...     for Cj in substitutionals:
...         Cj.equation.solve(var = Cj,
...                           dt = 10000)
...     if __name__ == '__main__':
...         for viewer in viewers:
...             viewer.plot()

Since there is nothing to maintain the concentration separation in this problem, we verify that the concentrations have become uniform

>>> substitutionals[0].allclose(0.45, rtol = tol, atol = tol).value
1
>>> substitutionals[1].allclose(0.45, rtol = tol, atol = tol).value
1

We now rerun the problem with an initial condition that only has a concentration step in one corner.

>>> x, y = mesh.cellCenters
>>> substitutionals[0].setValue(0.3)
>>> substitutionals[0].setValue(0.6, where=(x > L / 2.) & (y > L / 2.))
>>> substitutionals[1].setValue(0.6)
>>> substitutionals[1].setValue(0.3, where=(x > L / 2.) & (y > L / 2.))

We iterate the problem to equilibrium again

>>> from builtins import range
>>> for i in range(steps):
...     for Cj in substitutionals:
...         Cj.equation.solve(var = Cj,
...                           dt = 10000)
...     if __name__ == '__main__':
...         for viewer in viewers:
...             viewer.plot()

and verify that the correct uniform concentrations are achieved

>>> substitutionals[0].allclose(0.375, rtol = tol, atol = tol).value
1
>>> substitutionals[1].allclose(0.525, rtol = tol, atol = tol).value
1

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