examples.convection package¶

Subpackages¶

Submodules¶

examples.convection.peclet module¶

This example tests diffusion-convection for increasing Péclet numbers. This test case has been introduced because LinearCGSSolver was not working with Péclet numbers over 1. LinearLUSolver is now the default for ConvectionTerm. For nx = 1000 the LinearGMRESSolver does not work.

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

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

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

>>> var = CellVariable(name = "solution variable", mesh=mesh, value=valueLeft)

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

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

>>> convCoeff = 1.0
>>> peclet = 1e-3
>>> allcloseList = []
>>> from fipy import input
>>> from builtins import str
>>> while peclet < 1e4:
...     var[:] = valueLeft
...     diffCoeff = convCoeff * dx / peclet
...     eq = (TransientTerm(1e-4)
...           == DiffusionTerm(coeff=diffCoeff)
...           + PowerLawConvectionTerm(coeff=(convCoeff,)))
...     eq.solve(var=var, dt=1.)
...     x = mesh.cellCenters[0]
...     arg0 = -convCoeff * x / diffCoeff
...     arg0 = numerix.where(arg0 < -200, -200, arg0)
...     arg1 = -convCoeff * L / diffCoeff
...     arg1 = (arg1 >= -200) * (arg1 + 200) - 200
...     CC = 1. - numerix.exp(arg0)
...     DD = 1. - numerix.exp(arg1)
...     analyticalArray = CC / DD
...     allcloseList.append(var.allclose(CC / DD, rtol = 1e-2, atol = 1e-2).value)
...     if __name__ == '__main__':
...         viewer.plot()
...         input("Peclet number: " + str(peclet) + ", press key")
...     peclet *= 10

>>> print(allcloseList)
[True, True, True, True, True, True, True]

examples.convection.robin module¶

Solve an advection-diffusion equation with a Robin boundary condition.

This example demonstrates how to apply a Robin boundary condition to an advection-diffusion equation. The equation we wish to solve is given by,

$0 &=\frac{\partial^2 C}{\partial x^2} - P \frac{\partial C}{\partial x} -D C \qquad 0 < x <1\\ x=0: P &= -\frac{\partial C}{\partial x} + P C \\ x=1: \frac{\partial C}{\partial x} & = 0$

The analytical solution for this equation is given by,

$C \left( x \right) = \frac{ 2 P \exp{\left(\frac{P x}{2}\right)} \left[ \left(P + A \right) \exp{\left(\frac{A}{2} \left(1 - x\right)\right)} - \left(P - A \right) \exp{\left(-\frac{A}{2} \left(1 - x\right)\right)} \right]} { \left(P + A \right)^2 \exp{\left(\frac{A}{2}\right)} - \left(P - A \right)^2 \exp{\left(-\frac{A}{2}\right)}}$

where

$A = \sqrt{P^2 + 4D}$

>>> from fipy import CellVariable, FaceVariable, Grid1D, DiffusionTerm, PowerLawConvectionTerm, ImplicitSourceTerm, Viewer
>>> from fipy.tools import numerix
>>> nx = 100
>>> dx = 1.0 / nx

>>> mesh = Grid1D(nx=nx, dx=dx)
>>> C = CellVariable(mesh=mesh)

>>> D = 2.0
>>> P = 3.0

>>> C.faceGrad.constrain([0], mesh.facesRight)

We note that the Robin condition exactly defines the flux on the left, so we introduce a corresponding divergence source to the equation.

Note

Zeroing out the coefficients of the equation at this boundary is probably not necessary due to the default no-flux boundary conditions of cell-centered finite volume, but it’s a safe precaution.

>>> convectionCoeff = FaceVariable(mesh=mesh, value=[P])
>>> convectionCoeff[..., mesh.facesLeft.value] = 0.
>>> diffusionCoeff = FaceVariable(mesh=mesh, value=1.)
>>> diffusionCoeff[..., mesh.facesLeft.value] = 0.

>>> eq = (PowerLawConvectionTerm(coeff=convectionCoeff)
...       == DiffusionTerm(coeff=diffusionCoeff) - ImplicitSourceTerm(coeff=D)
...       - (P * mesh.facesLeft).divergence)

>>> A = numerix.sqrt(P**2 + 4 * D)

>>> x = mesh.cellCenters[0]
>>> CAnalytical = CellVariable(mesh=mesh)
>>> CAnalytical.setValue(2 * P * numerix.exp(P * x / 2)
...                      * ((P + A) * numerix.exp(A / 2 * (1 - x))
...                         - (P - A) * numerix.exp(-A / 2 *(1 - x)))
...                      / ((P + A)**2*numerix.exp(A / 2)
...                         - (P - A)**2 * numerix.exp(-A / 2)))

>>> if __name__ == '__main__':
...     C.name = '$C$'
...     CAnalytical.name = '$C_{analytical}$'
...     viewer = Viewer(vars=(C, CAnalytical))

>>> if __name__ == '__main__':
...     restol = 1e-5
...     anstol = 1e-3
... else:
...     restol = 0.5
...     anstol = 0.15

>>> res = 1e+10

>>> while res > restol:
...     res = eq.sweep(var=C)
...     if __name__ == '__main__':
...         viewer.plot()

>>> print(C.allclose(CAnalytical, rtol=anstol, atol=anstol))
True

examples.convection.source module¶

Solve a convection problem with a source.

This example solves the equation

$\frac{\partial \phi}{\partial x} + \alpha \phi = 0$

with $\phi \left( 0 \right) = 1$ at $x = 0$ . The boundary condition at $x = L$ is an outflow boundary condition requiring the use of an artificial constraint to be set on the right hand side faces. Exterior faces without constraints are considered to have zero outflow. An ImplicitSourceTerm object will be used to represent this term. The derivative of $\phi$ can be represented by a ConvectionTerm with a constant unitary velocity field from left to right. The following is an example code that includes a test against the analytical result.

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

>>> L = 10.
>>> nx = 5000
>>> dx =  L / nx
>>> mesh = Grid1D(dx=dx, nx=nx)
>>> phi0 = 1.0
>>> alpha = 1.0
>>> phi = CellVariable(name=r"$\phi$", mesh=mesh, value=phi0)
>>> solution = CellVariable(name=r"solution", mesh=mesh, value=phi0 * numerix.exp(-alpha * mesh.cellCenters[0]))

>>> from fipy import input
>>> if __name__ == "__main__":
...     viewer = Viewer(vars=(phi, solution))
...     viewer.plot()
...     input("press key to continue")

>>> phi.constrain(phi0, mesh.facesLeft)
>>> ## fake outflow condition
>>> phi.faceGrad.constrain([0], mesh.facesRight)

>>> eq = PowerLawConvectionTerm((1,)) + ImplicitSourceTerm(alpha)
>>> eq.solve(phi)
>>> print(numerix.allclose(phi, phi0 * numerix.exp(-alpha * mesh.cellCenters[0]), atol=1e-3))
True

>>> from fipy import input
>>> if __name__ == "__main__":
...     viewer = Viewer(vars=(phi, solution))
...     viewer.plot()
...     input("finished")

examples.convection.test module¶

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