fipy.terms.advectionTerm¶

Classes

AdvectionTerm([coeff])

The AdvectionTerm object constructs the b vector contribution for the advection term given by

class fipy.terms.advectionTerm.AdvectionTerm(coeff=None)¶

Bases: FirstOrderAdvectionTerm

The AdvectionTerm object constructs the b vector contribution for the advection term given by

\[u \abs{\nabla \phi}\]

from the advection equation given by:

\[\frac{\partial \phi}{\partial t} + u \abs{\nabla \phi} = 0\]

The construction of the gradient magnitude term requires upwinding as in the standard FirstOrderAdvectionTerm. The higher order terms are incorporated as follows. The formula used here is given by:

\[u_P \abs{\nabla \phi}_P = \max \left( u_P , 0 \right) \left[ \sum_A \min \left( D_{AP}, 0 \right)^2 \right]^{1/2} + \min \left( u_P , 0 \right) \left[ \sum_A \max \left( D_{AP}, 0 \right)^2 \right]^{1/2}\]

where,

\[D_{AP} = \frac{ \phi_A - \phi_P } { d_{AP}} - \frac{ d_{AP} } {2} m \left(L_A, L_P \right)\]

and

\[\begin{split}m\left(x, y\right) &= x \qquad \text{if $\abs{x} \le \abs{y} \forall xy \ge 0$} \\ m\left(x, y\right) &= y \qquad \text{if $\abs{x} > \abs{y} \forall xy \ge 0$} \\ m\left(x, y\right) &= 0 \qquad \text{if $xy < 0$}\end{split}\]

also,

\[\begin{split}L_A &= \frac{\phi_{AA} + \phi_P - 2 \phi_A}{d_{AP}^2} \\ L_P &= \frac{\phi_{A} + \phi_{PP} - 2 \phi_P}{d_{AP}^2}\end{split}\]

Here are some simple test cases for this problem:

>>> from fipy.meshes import Grid1D
>>> from fipy.solvers import *
>>> SparseMatrix = LinearPCGSolver()._matrixClass
>>> mesh = Grid1D(dx = 1., nx = 3)

Trivial test:

>>> from fipy.variables.cellVariable import CellVariable
>>> coeff = CellVariable(mesh = mesh, value = numerix.zeros(3, 'd'))
>>> v, L, b = AdvectionTerm(0.)._buildMatrix(coeff, SparseMatrix)
>>> print(numerix.allclose(b, numerix.zeros(3, 'd'), atol = 1e-10)) 
True

Less trivial test:

>>> coeff = CellVariable(mesh = mesh, value = numerix.arange(3))
>>> v, L, b = AdvectionTerm(1.)._buildMatrix(coeff, SparseMatrix)
>>> print(numerix.allclose(b, numerix.array((0., -1., -1.)), atol = 1e-10)) 
True

Even less trivial

>>> coeff = CellVariable(mesh = mesh, value = numerix.arange(3))
>>> v, L, b = AdvectionTerm(-1.)._buildMatrix(coeff, SparseMatrix)
>>> print(numerix.allclose(b, numerix.array((1., 1., 0.)), atol = 1e-10)) 
True

Another trivial test case (more trivial than a trivial test case standing on a harpsichord singing “trivial test cases are here again”)

>>> vel = numerix.array((-1, 2, -3))
>>> coeff = CellVariable(mesh = mesh, value = numerix.array((4, 6, 1)))
>>> v, L, b = AdvectionTerm(vel)._buildMatrix(coeff, SparseMatrix)
>>> print(numerix.allclose(b, -vel * numerix.array((2, numerix.sqrt(5**2 + 2**2), 5)), atol = 1e-10)) 
True

Somewhat less trivial test case:

>>> from fipy.meshes import Grid2D
>>> mesh = Grid2D(dx = 1., dy = 1., nx = 2, ny = 2)
>>> vel = numerix.array((3, -5, -6, -3))
>>> coeff = CellVariable(mesh = mesh, value = numerix.array((3, 1, 6, 7)))
>>> v, L, b = AdvectionTerm(vel)._buildMatrix(coeff, SparseMatrix)
>>> answer = -vel * numerix.array((2, numerix.sqrt(2**2 + 6**2), 1, 0))
>>> print(numerix.allclose(b, answer, atol = 1e-10)) 
True

For the above test cases the AdvectionTerm gives the same result as the AdvectionTerm. The following test imposes a quadratic field. The higher order term can resolve this field correctly.

\[\phi = x^2\]

The returned vector b should have the value:

\[-\abs{\nabla \phi} = -\left|\frac{\partial \phi}{\partial x}\right| = - 2 \abs{x}\]

Build the test case in the following way,

>>> mesh = Grid1D(dx = 1., nx = 5)
>>> vel = 1.
>>> coeff = CellVariable(mesh = mesh, value = mesh.cellCenters[0]**2)
>>> v, L, b = __AdvectionTerm(vel)._buildMatrix(coeff, SparseMatrix)

The first order term is not accurate. The first and last element are ignored because they don’t have any neighbors for higher order evaluation

>>> print(numerix.allclose(CellVariable(mesh=mesh,
... value=b).globalValue[1:-1], -2 * mesh.cellCenters.globalValue[0][1:-1]))
False

The higher order term is spot on.

>>> v, L, b = AdvectionTerm(vel)._buildMatrix(coeff, SparseMatrix)
>>> print(numerix.allclose(CellVariable(mesh=mesh,
... value=b).globalValue[1:-1], -2 * mesh.cellCenters.globalValue[0][1:-1]))
True

The AdvectionTerm will also resolve a circular field with more accuracy,

\[\phi = \left( x^2 + y^2 \right)^{1/2}\]

Build the test case in the following way,

>>> mesh = Grid2D(dx = 1., dy = 1., nx = 10, ny = 10)
>>> vel = 1.
>>> x, y = mesh.cellCenters
>>> r = numerix.sqrt(x**2 + y**2)
>>> coeff = CellVariable(mesh = mesh, value = r)
>>> v, L, b = __AdvectionTerm(1.)._buildMatrix(coeff, SparseMatrix)
>>> error = CellVariable(mesh=mesh, value=b + 1)
>>> ans = CellVariable(mesh=mesh, value=b + 1)
>>> ans[(x > 2) & (x < 8) & (y > 2) & (y < 8)] = 0.123105625618
>>> print((error <= ans).all())
True

The maximum error is large (about 12 %) for the first order advection.

>>> v, L, b = AdvectionTerm(1.)._buildMatrix(coeff, SparseMatrix)
>>> error = CellVariable(mesh=mesh, value=b + 1)
>>> ans = CellVariable(mesh=mesh, value=b + 1)
>>> ans[(x > 2) & (x < 8) & (y > 2) & (y < 8)] = 0.0201715476598
>>> print((error <= ans).all())
True

The maximum error is 2 % when using a higher order contribution.

Create a Term.

Parameters:: coeff (float or CellVariable or FaceVariable) – Coefficient for the term. FaceVariable objects are only acceptable for diffusion or convection terms.

property RHSvector¶: Return the RHS vector calculated in solve() or sweep(). The cacheRHSvector() method should be called before solve() or sweep() to cache the vector.

__eq__(other)¶: Return self==value.

__hash__()¶: Return hash(self).

__mul__(other)¶

Multiply a term

>>> 2. * __NonDiffusionTerm(coeff=0.5)
__NonDiffusionTerm(coeff=1.0)

Test for ticket:291.

>>> from fipy import PowerLawConvectionTerm
>>> PowerLawConvectionTerm(coeff=[[1], [0]]) * 1.0
PowerLawConvectionTerm(coeff=array([[ 1.],
       [ 0.]]))

__neg__()¶

Negate a Term.

>>> -__NonDiffusionTerm(coeff=1.)
__NonDiffusionTerm(coeff=-1.0)

__repr__()¶

The representation of a Term object is given by,

>>> print(__UnaryTerm(123.456))
__UnaryTerm(coeff=123.456)

__rmul__(other)¶

Multiply a term

>>> 2. * __NonDiffusionTerm(coeff=0.5)
__NonDiffusionTerm(coeff=1.0)

Test for ticket:291.

>>> from fipy import PowerLawConvectionTerm
>>> PowerLawConvectionTerm(coeff=[[1], [0]]) * 1.0
PowerLawConvectionTerm(coeff=array([[ 1.],
       [ 0.]]))

cacheMatrix()¶: Informs solve() and sweep() to cache their matrix so that matrix can return the matrix.

cacheRHSvector()¶: Informs solve() and sweep() to cache their right hand side vector so that getRHSvector() can return it.

justErrorVector(var=None, solver=None, boundaryConditions=(), dt=1.0, underRelaxation=None, residualFn=None)¶

Builds the Term’s linear system once.

This method also recalculates and returns the error as well as applying under-relaxation.

justErrorVector returns the overlapping local value in parallel (not the non-overlapping value).

>>> from fipy.solvers import DummySolver
>>> from fipy import *
>>> m = Grid1D(nx=10)
>>> v = CellVariable(mesh=m)
>>> len(DiffusionTerm().justErrorVector(v, solver=DummySolver())) == m.numberOfCells
True

Parameters:

var (CellVariable) – Variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
solver (Solver) – Iterative solver to be used to solve the linear system of equations. The default sovler depends on the solver package selected.
boundaryConditions (tuple of BoundaryCondition) –
dt (float) – Timestep size.
underRelaxation (float) – Usually a value between 0 and 1 or None in the case of no under-relaxation
residualFn (function) – Takes var, matrix, and RHSvector arguments, used to customize the residual calculation.

Returns:

error – The residual vector $\vec{e}=\mathsf{L}\vec{x}_\text{old} - \vec{b}$

Return type:

CellVariable

justResidualVector(var=None, solver=None, boundaryConditions=(), dt=None, underRelaxation=None, residualFn=None)¶

Builds the Term’s linear system once.

This method also recalculates and returns the residual as well as applying under-relaxation.

justResidualVector returns the overlapping local value in parallel (not the non-overlapping value).

>>> from fipy import *
>>> m = Grid1D(nx=10)
>>> v = CellVariable(mesh=m)
>>> len(numerix.asarray(DiffusionTerm().justResidualVector(v))) == m.numberOfCells
True

Parameters:

var (CellVariable) – Variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
solver (Solver) – Iterative solver to be used to solve the linear system of equations. The default sovler depends on the solver package selected.
boundaryConditions (tuple of BoundaryCondition) –
dt (float) – Timestep size.
underRelaxation (float) – Usually a value between 0 and 1 or None in the case of no under-relaxation
residualFn (function) – Takes var, matrix, and RHSvector arguments, used to customize the residual calculation.

Returns:

residual – The residual vector $\vec{r}=\mathsf{L}\vec{x} - \vec{b}$

Return type:

CellVariable

property matrix¶: Return the matrix calculated in solve() or sweep(). The cacheMatrix() method should be called before solve() or sweep() to cache the matrix.

residualVectorAndNorm(var=None, solver=None, boundaryConditions=(), dt=None, underRelaxation=None, residualFn=None)¶

Builds the Term’s linear system once.

This method also recalculates and returns the residual as well as applying under-relaxation.

Parameters:

var (CellVariable) – Variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
solver (Solver) – Iterative solver to be used to solve the linear system of equations. The default sovler depends on the solver package selected.
boundaryConditions (tuple of BoundaryCondition) –
dt (float) – Timestep size.
underRelaxation (float) – Usually a value between 0 and 1 or None in the case of no under-relaxation
residualFn (function) – Takes var, matrix, and RHSvector arguments, used to customize the residual calculation.

Returns:

residual (~fipy.variables.cellVariable.CellVariable) – The residual vector $\vec{r}=\mathsf{L}\vec{x} - \vec{b}$
norm (float) – The L2 norm of residual, $\|\vec{r}\|_2$

solve(var=None, solver=None, boundaryConditions=(), dt=None)¶

Builds and solves the Term’s linear system once. This method does not return the residual. It should be used when the residual is not required.

Parameters:

var (CellVariable) – Variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
solver (Solver) – Iterative solver to be used to solve the linear system of equations. The default sovler depends on the solver package selected.
boundaryConditions (tuple of BoundaryCondition) –
dt (float) – Timestep size.

sweep(var=None, solver=None, boundaryConditions=(), dt=None, underRelaxation=None, residualFn=None, cacheResidual=False, cacheError=False)¶

Builds and solves the Term’s linear system once. This method also recalculates and returns the residual as well as applying under-relaxation.

Parameters:

var (CellVariable) – Variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
solver (Solver) – Iterative solver to be used to solve the linear system of equations. The default sovler depends on the solver package selected.
boundaryConditions (tuple of BoundaryCondition) –
dt (float) – Timestep size.
underRelaxation (float) – Usually a value between 0 and 1 or None in the case of no under-relaxation
residualFn (function) – Takes var, matrix, and RHSvector arguments, used to customize the residual calculation.
cacheResidual (bool) – If True, calculate and store the residual vector $\vec{r}=\mathsf{L}\vec{x} - \vec{b}$ in the residualVector member of Term
cacheError (bool) – If True, use the residual vector $\vec{r}$ to solve $\mathsf{L}\vec{e}=\vec{r}$ for the error vector $\vec{e}$ and store it in the errorVector member of Term

Returns:

residual – The residual vector $\vec{r}=\mathsf{L}\vec{x} - \vec{b}$

Return type:

CellVariable

Last updated on Jun 26, 2024. Created using Sphinx 7.1.2.