examples.convection.exponential1D package

Submodules

examples.convection.exponential1D.cylindricalMesh1D module

This example solves the steady-state cylindrical convection-diffusion equation given by

\nabla \cdot \left(D \nabla \phi + \vec{u} \phi \right) = 0

with coefficients D = 1 and \vec{u} = (10,), or

>>> diffCoeff = 1.
>>> convCoeff = (10.,)

We define a 1D cylindrical mesh representing an annulus

>>> from fipy import CellVariable, CylindricalGrid1D, DiffusionTerm, ExponentialConvectionTerm, Viewer
>>> from fipy.tools import numerix

>>> r0 = 1.
>>> r1 = 2.
>>> nr = 100
>>> mesh = CylindricalGrid1D(dr=(r1 - r0) / nr, nr=nr) + ((r0,),)

The solution variable is initialized to valueLeft:

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

>>> var = CellVariable(mesh=mesh, name = "variable")

and impose the boundary conditions

\phi = \begin{cases} 0& \text{at $r = r_0$,} \\ 1& \text{at $r = r_1$,} \end{cases}

with

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

The equation is created with the DiffusionTerm and ExponentialConvectionTerm.

>>> eq = (DiffusionTerm(coeff=diffCoeff)
...       + ExponentialConvectionTerm(coeff=convCoeff))

More details of the benefits and drawbacks of each type of convection term can be found in Numerical Schemes. Essentially, the ExponentialConvectionTerm and PowerLawConvectionTerm will both handle most types of convection-diffusion cases, with the PowerLawConvectionTerm being more efficient.

We solve the equation

>>> eq.solve(var=var)

and test the solution against the analytical result

or

>>> axis = 0
>>> try:
...     from scipy.special import expi 
...     r = mesh.cellCenters[axis]
...     AA = numerix.exp(convCoeff[axis] / diffCoeff * (r1 - r))
...     BB = expi(convCoeff[axis] * r0 / diffCoeff) - expi(convCoeff[axis] * r / diffCoeff) 
...     CC = expi(convCoeff[axis] * r0 / diffCoeff) - expi(convCoeff[axis] * r1 / diffCoeff) 
...     analyticalArray = AA * BB / CC 
... except ImportError:
...     print("The SciPy library is unavailable. It is required for testing purposes.")

>>> print(var.allclose(analyticalArray, atol=1e-3)) 
1

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

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

examples.convection.exponential1D.cylindricalMesh1DNonUniform module

This example solves the steady-state cylindrical convection-diffusion equation given by

\nabla \cdot \left(D \nabla \phi + \vec{u} \phi \right) = 0

with coefficients D = 1 and \vec{u} = (10,), or

>>> diffCoeff = 1.
>>> convCoeff = ((10.,),)

We define a 1D cylindrical mesh representing an annulus. The mesh has a non-constant cell spacing.

>>> from fipy import CellVariable, CylindricalGrid1D, DiffusionTerm, ExponentialConvectionTerm, Viewer
>>> from fipy.tools import numerix

>>> r0 = 1.
>>> r1 = 2.
>>> nr = 100
>>> Rratio = (r1 / r0)**(1 / float(nr))
>>> dr = r0 * (Rratio - 1) * Rratio**numerix.arange(nr)
>>> mesh = CylindricalGrid1D(dr=dr) + ((r0,),)

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

The solution variable is initialized to valueLeft:

>>> var = CellVariable(mesh=mesh, name = "variable")

and impose the boundary conditions

with

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

The equation is created with the DiffusionTerm and ExponentialConvectionTerm.

>>> eq = (DiffusionTerm(coeff=diffCoeff)
...       + ExponentialConvectionTerm(coeff=convCoeff))

More details of the benefits and drawbacks of each type of convection term can be found in Numerical Schemes. Essentially, the ExponentialConvectionTerm and PowerLawConvectionTerm will both handle most types of convection-diffusion cases, with the PowerLawConvectionTerm being more efficient.

We solve the equation

>>> eq.solve(var=var)

and test the solution against the analytical result

or

>>> axis = 0

>>> try:
...     U = convCoeff[0][0]
...     from scipy.special import expi 
...     r = mesh.cellCenters[axis]
...     AA = numerix.exp(U / diffCoeff * (r1 - r))
...     BB = expi(U * r0 / diffCoeff) - expi(U * r / diffCoeff) 
...     CC = expi(U * r0 / diffCoeff) - expi(U * r1 / diffCoeff) 
...     analyticalArray = AA * BB / CC 
... except ImportError:
...     print("The SciPy library is unavailable. It is required for testing purposes.")

>>> print(var.allclose(analyticalArray, atol=1e-3)) 
1

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

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

examples.convection.exponential1D.mesh1D module

Solve the steady-state convection-diffusion equation in one dimension.

This example solves the steady-state convection-diffusion equation given by

\nabla \cdot \left(D \nabla \phi + \vec{u} \phi \right) = 0

with coefficients D = 1 and \vec{u} = 10\hat{\i}, or

>>> diffCoeff = 1.
>>> convCoeff = (10.,)

We define a 1D mesh

>>> from fipy import CellVariable, Grid1D, DiffusionTerm, ExponentialConvectionTerm, Viewer
>>> from fipy.tools import numerix

>>> L = 10.
>>> nx = 10
>>> mesh = Grid1D(dx=L / nx, nx=nx)

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

The solution variable is initialized to valueLeft:

>>> var = CellVariable(mesh=mesh, name="variable")

and impose the boundary conditions

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

with

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

The equation is created with the DiffusionTerm and ExponentialConvectionTerm. The scheme used by the convection term needs to calculate a Péclet number and thus the diffusion term instance must be passed to the convection term.

>>> eq = (DiffusionTerm(coeff=diffCoeff)
...       + ExponentialConvectionTerm(coeff=convCoeff))

More details of the benefits and drawbacks of each type of convection term can be found in Numerical Schemes. Essentially, the ExponentialConvectionTerm and PowerLawConvectionTerm will both handle most types of convection-diffusion cases, with the PowerLawConvectionTerm being more efficient.

We solve the equation

>>> eq.solve(var=var)

and test the solution against the analytical result

\phi = \frac{1 - \exp(-u_x x / D)}{1 - \exp(-u_x L / D)}

or

>>> axis = 0
>>> x = mesh.cellCenters[axis]
>>> CC = 1. - numerix.exp(-convCoeff[axis] * x / diffCoeff)
>>> DD = 1. - numerix.exp(-convCoeff[axis] * L / diffCoeff)
>>> analyticalArray = CC / DD
>>> print(var.allclose(analyticalArray))
1

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

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

examples.convection.exponential1D.tri2D module

This example solves the steady-state convection-diffusion equation as described in examples.diffusion.convection.exponential1D.input but uses a Tri2D mesh.

Here the axes are reversed (nx = 1, ny = 1000) and

\vec{u} = (0, 10)

>>> from fipy import CellVariable, Tri2D, DiffusionTerm, ExponentialConvectionTerm, DefaultAsymmetricSolver, Viewer
>>> from fipy.tools import numerix

>>> L = 10.
>>> nx = 1
>>> ny = 1000
>>> mesh = Tri2D(dx = L / ny, dy = L / ny, nx = nx, ny = ny)

>>> valueBottom = 0.
>>> valueTop = 1.

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

>>> var.constrain(valueBottom, mesh.facesBottom)
>>> var.constrain(valueTop, mesh.facesTop)

>>> diffCoeff = 1.
>>> convCoeff = numerix.array(((0.,), (10.,)))

>>> eq = (DiffusionTerm(coeff=diffCoeff)
...       + ExponentialConvectionTerm(coeff=convCoeff))

>>> eq.solve(var = var,
...          solver=DefaultAsymmetricSolver(iterations=10000))

The analytical solution test for this problem is given by:

>>> axis = 1
>>> y = mesh.cellCenters[axis]
>>> CC = 1. - numerix.exp(-convCoeff[axis] * y / diffCoeff)
>>> DD = 1. - numerix.exp(-convCoeff[axis] * L / diffCoeff)
>>> analyticalArray = CC / DD

>>> print(var.allclose(analyticalArray, rtol = 1e-6, atol = 1e-6))
1

>>> if __name__ == '__main__':
...     viewer = Viewer(vars = var)
...     viewer.plot()
Last updated on Jun 27, 2023. Created using Sphinx 6.2.1.