Navigation

Bookmark and Share FiPy: A Finite Volume PDE Solver Using Python
Version 3.0.1-dev139-ge5d2233

Table Of Contents

Previous topic

steppers Package

Next topic

tests Package

This Page

Contact

FiPy developers
Jonathan Guyer
Daniel Wheeler
James Warren

Join our mailing list

100 Bureau Drive, M/S 6555
Gaithersburg, MD 20899

301-975-5329 Telephone
301-975-4553 Facsimile

terms Package

terms Package

exception fipy.terms.ExplicitVariableError(s='Terms with explicit Variables cannot mix with Terms with implicit Variables.')

Bases: exceptions.Exception

exception fipy.terms.TermMultiplyError(s='Must multiply terms by int or float.')

Bases: exceptions.Exception

exception fipy.terms.AbstractBaseClassError(s="can't instantiate abstract base class")

Bases: exceptions.NotImplementedError

exception fipy.terms.VectorCoeffError(s='The coefficient must be a vector value.')

Bases: exceptions.TypeError

exception fipy.terms.SolutionVariableNumberError(s='Different number of solution variables and equations.')

Bases: exceptions.Exception

exception fipy.terms.SolutionVariableRequiredError(s='The solution variable needs to be specified.')

Bases: exceptions.Exception

exception fipy.terms.IncorrectSolutionVariable(s='The solution variable is incorrect.')

Bases: exceptions.Exception

fipy.terms.ConvectionTerm

alias of PowerLawConvectionTerm

class fipy.terms.FirstOrderAdvectionTerm(coeff=None)

Bases: fipy.terms.nonDiffusionTerm._NonDiffusionTerm

The FirstOrderAdvectionTerm 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. The formula used here is given by:

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

Here are some simple test cases for this problem:

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

Trivial test:

>>> var = CellVariable(value = numerix.zeros(3, 'd'), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(0.)._buildMatrix(var, SparseMatrix)
>>> print numerix.allclose(b, numerix.zeros(3, 'd'), atol = 1e-10) 
True

Less trivial test:

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

Even less trivial

>>> var = CellVariable(value = numerix.arange(3), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(-1.)._buildMatrix(var, 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))
>>> var = CellVariable(value = numerix.array((4,6,1)), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(vel)._buildMatrix(var, 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))
>>> var = CellVariable(value = numerix.array((3 , 1, 6, 7)), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(vel)._buildMatrix(var, SparseMatrix)
>>> answer = -vel * numerix.array((2, numerix.sqrt(2**2 + 6**2), 1, 0))
>>> print numerix.allclose(b, answer, atol = 1e-10) 
True
class fipy.terms.AdvectionTerm(coeff=None)

Bases: fipy.terms.firstOrderAdvectionTerm.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 incorperated 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

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$}

also,

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}

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.

class fipy.terms.TransientTerm(coeff=1.0, var=None)

Bases: fipy.terms.cellTerm.CellTerm

The TransientTerm represents

\int_V \frac{\partial (\rho \phi)}{\partial t} dV \simeq \frac{(\rho_{P} \phi_{P} - \rho_{P}^\text{old} \phi_P^\text{old}) V_P}{\Delta t}

where \rho is the coeff value.

The following test case verifies that variable coefficients and old coefficient values work correctly. We will solve the following equation

\frac{ \partial \phi^2 } { \partial t } = k.

The analytic solution is given by

\phi = \sqrt{ \phi_0^2 + k t },

where \phi_0 is the initial value.

>>> phi0 = 1.
>>> k = 1.
>>> dt = 1.
>>> relaxationFactor = 1.5
>>> steps = 2
>>> sweeps = 8

>>> from fipy.meshes import Grid1D
>>> mesh = Grid1D(nx = 1)
>>> from fipy.variables.cellVariable import CellVariable
>>> var = CellVariable(mesh = mesh, value = phi0, hasOld = 1)
>>> from fipy.terms.transientTerm import TransientTerm
>>> from fipy.terms.implicitSourceTerm import ImplicitSourceTerm

Relaxation, given by relaxationFactor, is required for a converged solution.

>>> eq = TransientTerm(var) == ImplicitSourceTerm(-relaxationFactor) \
...                            + var * relaxationFactor + k

A number of sweeps at each time step are required to let the relaxation take effect.

>>> for step in range(steps):
...     var.updateOld()
...     for sweep in range(sweeps):
...         eq.solve(var, dt = dt)

Compare the final result with the analytical solution.

>>> from fipy.tools import numerix
>>> print var.allclose(numerix.sqrt(k * dt * steps + phi0**2))
1
class fipy.terms.DiffusionTerm(coeff=(1.0, ), var=None)

Bases: fipy.terms.diffusionTermNoCorrection.DiffusionTermNoCorrection

This term represents a higher order diffusion term. The order of the term is determined by the number of coeffs, such that:

DiffusionTerm(D1)

represents a typical 2nd-order diffusion term of the form

\nabla\cdot\left(D_1 \nabla \phi\right)

and:

DiffusionTerm((D1,D2))

represents a 4th-order Cahn-Hilliard term of the form

\nabla \cdot \left\{ D_1 \nabla \left[ \nabla\cdot\left( D_2 \nabla \phi\right) \right] \right\}

and so on.

class fipy.terms.DiffusionTermCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

class fipy.terms.DiffusionTermNoCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

class fipy.terms.DiffusionTermCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

class fipy.terms.DiffusionTermNoCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

class fipy.terms.ExplicitDiffusionTerm(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

The discretization for the ExplicitDiffusionTerm is given by

\int_V \nabla \cdot (\Gamma\nabla\phi) dV \simeq \sum_f \Gamma_f \frac{\phi_A^\text{old}-\phi_P^\text{old}}{d_{AP}} A_f

where \phi_A^\text{old} and \phi_P^\text{old} are the old values of the variable. The term is added to the RHS vector and makes no contribution to the solution matrix.

fipy.terms.ImplicitDiffusionTerm

alias of DiffusionTerm

class fipy.terms.ImplicitSourceTerm(coeff=0.0, var=None)

Bases: fipy.terms.sourceTerm.SourceTerm

The ImplicitSourceTerm represents

\int_V \phi S \,dV \simeq \phi_P S_P V_P

where S is the coeff value.

class fipy.terms.ResidualTerm(equation, underRelaxation=1.0)

Bases: fipy.terms.explicitSourceTerm._ExplicitSourceTerm

The ResidualTerm is a special form of explicit SourceTerm that adds the residual of one equation to another equation. Useful for Newton’s method.

class fipy.terms.CentralDifferenceConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.abstractConvectionTerm._AbstractConvectionTerm

This Term represents

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the central differencing scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.
class fipy.terms.ExplicitUpwindConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.abstractUpwindConvectionTerm._AbstractUpwindConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P^\text{old} +(1-\alpha_f)\phi_A^\text{old} and \alpha_f is calculated using the upwind scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.
class fipy.terms.ExponentialConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.asymmetricConvectionTerm._AsymmetricConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the exponential scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.
class fipy.terms.HybridConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.asymmetricConvectionTerm._AsymmetricConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the hybrid scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.
class fipy.terms.PowerLawConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.asymmetricConvectionTerm._AsymmetricConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the power law scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.
class fipy.terms.UpwindConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.abstractUpwindConvectionTerm._AbstractUpwindConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the upwind convection scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.
class fipy.terms.VanLeerConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.explicitUpwindConvectionTerm.ExplicitUpwindConvectionTerm

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.
class fipy.terms.FirstOrderAdvectionTerm(coeff=None)

Bases: fipy.terms.nonDiffusionTerm._NonDiffusionTerm

The FirstOrderAdvectionTerm 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. The formula used here is given by:

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

Here are some simple test cases for this problem:

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

Trivial test:

>>> var = CellVariable(value = numerix.zeros(3, 'd'), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(0.)._buildMatrix(var, SparseMatrix)
>>> print numerix.allclose(b, numerix.zeros(3, 'd'), atol = 1e-10) 
True

Less trivial test:

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

Even less trivial

>>> var = CellVariable(value = numerix.arange(3), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(-1.)._buildMatrix(var, 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))
>>> var = CellVariable(value = numerix.array((4,6,1)), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(vel)._buildMatrix(var, 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))
>>> var = CellVariable(value = numerix.array((3 , 1, 6, 7)), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(vel)._buildMatrix(var, SparseMatrix)
>>> answer = -vel * numerix.array((2, numerix.sqrt(2**2 + 6**2), 1, 0))
>>> print numerix.allclose(b, answer, atol = 1e-10) 
True
class fipy.terms.AdvectionTerm(coeff=None)

Bases: fipy.terms.firstOrderAdvectionTerm.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 incorperated 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

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$}

also,

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}

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.

abstractBinaryTerm Module

abstractConvectionTerm Module

abstractDiffusionTerm Module

abstractUpwindConvectionTerm Module

advectionTerm Module

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

Bases: fipy.terms.firstOrderAdvectionTerm.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 incorperated 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

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$}

also,

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}

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.

asymmetricConvectionTerm Module

binaryTerm Module

cellTerm Module

class fipy.terms.cellTerm.CellTerm(coeff=1.0, var=None)

Bases: fipy.terms.nonDiffusionTerm._NonDiffusionTerm

Attention

This class is abstract. Always create one of its subclasses.

centralDiffConvectionTerm Module

class fipy.terms.centralDiffConvectionTerm.CentralDifferenceConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.abstractConvectionTerm._AbstractConvectionTerm

This Term represents

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the central differencing scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.

coupledBinaryTerm Module

diffusionTerm Module

class fipy.terms.diffusionTerm.DiffusionTerm(coeff=(1.0, ), var=None)

Bases: fipy.terms.diffusionTermNoCorrection.DiffusionTermNoCorrection

This term represents a higher order diffusion term. The order of the term is determined by the number of coeffs, such that:

DiffusionTerm(D1)

represents a typical 2nd-order diffusion term of the form

\nabla\cdot\left(D_1 \nabla \phi\right)

and:

DiffusionTerm((D1,D2))

represents a 4th-order Cahn-Hilliard term of the form

\nabla \cdot \left\{ D_1 \nabla \left[ \nabla\cdot\left( D_2 \nabla \phi\right) \right] \right\}

and so on.

class fipy.terms.diffusionTerm.DiffusionTermCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

class fipy.terms.diffusionTerm.DiffusionTermNoCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

diffusionTermCorrection Module

class fipy.terms.diffusionTermCorrection.DiffusionTermCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

diffusionTermNoCorrection Module

class fipy.terms.diffusionTermNoCorrection.DiffusionTermNoCorrection(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

explicitDiffusionTerm Module

class fipy.terms.explicitDiffusionTerm.ExplicitDiffusionTerm(coeff=(1.0, ), var=None)

Bases: fipy.terms.abstractDiffusionTerm._AbstractDiffusionTerm

The discretization for the ExplicitDiffusionTerm is given by

\int_V \nabla \cdot (\Gamma\nabla\phi) dV \simeq \sum_f \Gamma_f \frac{\phi_A^\text{old}-\phi_P^\text{old}}{d_{AP}} A_f

where \phi_A^\text{old} and \phi_P^\text{old} are the old values of the variable. The term is added to the RHS vector and makes no contribution to the solution matrix.

explicitSourceTerm Module

explicitUpwindConvectionTerm Module

class fipy.terms.explicitUpwindConvectionTerm.ExplicitUpwindConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.abstractUpwindConvectionTerm._AbstractUpwindConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P^\text{old} +(1-\alpha_f)\phi_A^\text{old} and \alpha_f is calculated using the upwind scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.

exponentialConvectionTerm Module

class fipy.terms.exponentialConvectionTerm.ExponentialConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.asymmetricConvectionTerm._AsymmetricConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the exponential scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.

faceTerm Module

class fipy.terms.faceTerm.FaceTerm(coeff=1.0, var=None)

Bases: fipy.terms.nonDiffusionTerm._NonDiffusionTerm

Attention

This class is abstract. Always create one of its subclasses.

firstOrderAdvectionTerm Module

class fipy.terms.firstOrderAdvectionTerm.FirstOrderAdvectionTerm(coeff=None)

Bases: fipy.terms.nonDiffusionTerm._NonDiffusionTerm

The FirstOrderAdvectionTerm 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. The formula used here is given by:

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

Here are some simple test cases for this problem:

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

Trivial test:

>>> var = CellVariable(value = numerix.zeros(3, 'd'), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(0.)._buildMatrix(var, SparseMatrix)
>>> print numerix.allclose(b, numerix.zeros(3, 'd'), atol = 1e-10) 
True

Less trivial test:

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

Even less trivial

>>> var = CellVariable(value = numerix.arange(3), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(-1.)._buildMatrix(var, 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))
>>> var = CellVariable(value = numerix.array((4,6,1)), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(vel)._buildMatrix(var, 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))
>>> var = CellVariable(value = numerix.array((3 , 1, 6, 7)), mesh = mesh)
>>> v, L, b = FirstOrderAdvectionTerm(vel)._buildMatrix(var, SparseMatrix)
>>> answer = -vel * numerix.array((2, numerix.sqrt(2**2 + 6**2), 1, 0))
>>> print numerix.allclose(b, answer, atol = 1e-10) 
True

hybridConvectionTerm Module

class fipy.terms.hybridConvectionTerm.HybridConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.asymmetricConvectionTerm._AsymmetricConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the hybrid scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.

implicitDiffusionTerm Module

fipy.terms.implicitDiffusionTerm.ImplicitDiffusionTerm

alias of DiffusionTerm

implicitSourceTerm Module

class fipy.terms.implicitSourceTerm.ImplicitSourceTerm(coeff=0.0, var=None)

Bases: fipy.terms.sourceTerm.SourceTerm

The ImplicitSourceTerm represents

\int_V \phi S \,dV \simeq \phi_P S_P V_P

where S is the coeff value.

nonDiffusionTerm Module

powerLawConvectionTerm Module

class fipy.terms.powerLawConvectionTerm.PowerLawConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.asymmetricConvectionTerm._AsymmetricConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the power law scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.

residualTerm Module

class fipy.terms.residualTerm.ResidualTerm(equation, underRelaxation=1.0)

Bases: fipy.terms.explicitSourceTerm._ExplicitSourceTerm

The ResidualTerm is a special form of explicit SourceTerm that adds the residual of one equation to another equation. Useful for Newton’s method.

sourceTerm Module

class fipy.terms.sourceTerm.SourceTerm(coeff=0.0, var=None)

Bases: fipy.terms.cellTerm.CellTerm

Attention

This class is abstract. Always create one of its subclasses.

term Module

class fipy.terms.term.Term(coeff=1.0, var=None)

Bases: object

Attention

This class is abstract. Always create one of its subclasses.

Create a Term.

Parameters :
  • coeff: The coefficient for the term. A CellVariable or number. FaceVariable objects are also acceptable for diffusion or convection terms.
RHSvector

Return the RHS vector caculated in solve() or sweep(). The cacheRHSvector() method should be called before solve() or sweep() to cache the vector.

cacheMatrix()

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

cacheRHSvector()

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

copy()
getDefaultSolver(var=None, solver=None, *args, **kwargs)
getMatrix(*args, **kwds)

Deprecated since version 3.0: use the matrix property instead

getRHSvector(*args, **kwds)

Deprecated since version 3.0: use the rHSvector property instead

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.

Parameters :
  • var: The variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
  • solver: The iterative solver to be used to solve the linear system of equations. Defaults to LinearPCGSolver for Pysparse and LinearLUSolver for Trilinos.
  • boundaryConditions: A tuple of boundaryConditions.
  • dt: The time step size.
  • underRelaxation: Usually a value between 0 and 1 or None in the case of no under-relaxation
  • residualFn: A function that takes var, matrix, and RHSvector arguments used to customize the residual calculation.

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
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.

Parameters :
  • var: The variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
  • solver: The iterative solver to be used to solve the linear system of equations. Defaults to LinearPCGSolver for Pysparse and LinearLUSolver for Trilinos.
  • boundaryConditions: A tuple of boundaryConditions.
  • dt: The time step size.
  • underRelaxation: Usually a value between 0 and 1 or None in the case of no under-relaxation
  • residualFn: A function that takes var, matrix, and RHSvector arguments used to customize the residual calculation.

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(DiffusionTerm().justResidualVector(v)) == m.numberOfCells
True
matrix

Return the matrix caculated 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: The variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
  • solver: The iterative solver to be used to solve the linear system of equations. Defaults to LinearPCGSolver for Pysparse and LinearLUSolver for Trilinos.
  • boundaryConditions: A tuple of boundaryConditions.
  • dt: The time step size.
  • underRelaxation: Usually a value between 0 and 1 or None in the case of no under-relaxation
  • residualFn: A function that takes var, matrix, and RHSvector arguments used to customize the residual calculation.
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: The variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.
  • solver: The iterative solver to be used to solve the linear system of equations. Defaults to LinearPCGSolver for Pysparse and LinearLUSolver for Trilinos.
  • boundaryConditions: A tuple of boundaryConditions.
  • dt: The time step 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: The variable to be solved for. Provides the initial condition, the old value and holds the solution on completion.

  • solver: The iterative solver to be used to solve the linear system of equations. Defaults to LinearPCGSolver for Pysparse and LinearLUSolver for Trilinos.

  • boundaryConditions: A tuple of boundaryConditions.

  • dt: The time step size.

  • underRelaxation: Usually a value between 0 and 1 or None in the case of no under-relaxation

  • residualFn: A function that takes var, matrix, and RHSvector arguments, used to customize the residual calculation.

  • cacheResidual: If True, calculate and store the residual vector

    \vec{r}=\mathsf{L}\vec{x} - \vec{b} in the residualVector member of Term

  • cacheError: 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

test Module

transientTerm Module

class fipy.terms.transientTerm.TransientTerm(coeff=1.0, var=None)

Bases: fipy.terms.cellTerm.CellTerm

The TransientTerm represents

\int_V \frac{\partial (\rho \phi)}{\partial t} dV \simeq \frac{(\rho_{P} \phi_{P} - \rho_{P}^\text{old} \phi_P^\text{old}) V_P}{\Delta t}

where \rho is the coeff value.

The following test case verifies that variable coefficients and old coefficient values work correctly. We will solve the following equation

\frac{ \partial \phi^2 } { \partial t } = k.

The analytic solution is given by

\phi = \sqrt{ \phi_0^2 + k t },

where \phi_0 is the initial value.

>>> phi0 = 1.
>>> k = 1.
>>> dt = 1.
>>> relaxationFactor = 1.5
>>> steps = 2
>>> sweeps = 8

>>> from fipy.meshes import Grid1D
>>> mesh = Grid1D(nx = 1)
>>> from fipy.variables.cellVariable import CellVariable
>>> var = CellVariable(mesh = mesh, value = phi0, hasOld = 1)
>>> from fipy.terms.transientTerm import TransientTerm
>>> from fipy.terms.implicitSourceTerm import ImplicitSourceTerm

Relaxation, given by relaxationFactor, is required for a converged solution.

>>> eq = TransientTerm(var) == ImplicitSourceTerm(-relaxationFactor) \
...                            + var * relaxationFactor + k

A number of sweeps at each time step are required to let the relaxation take effect.

>>> for step in range(steps):
...     var.updateOld()
...     for sweep in range(sweeps):
...         eq.solve(var, dt = dt)

Compare the final result with the analytical solution.

>>> from fipy.tools import numerix
>>> print var.allclose(numerix.sqrt(k * dt * steps + phi0**2))
1

unaryTerm Module

upwindConvectionTerm Module

class fipy.terms.upwindConvectionTerm.UpwindConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.abstractUpwindConvectionTerm._AbstractUpwindConvectionTerm

The discretization for this Term is given by

\int_V \nabla \cdot (\vec{u} \phi)\,dV \simeq \sum_{f} (\vec{n} \cdot \vec{u})_f \phi_f A_f

where \phi_f=\alpha_f \phi_P +(1-\alpha_f)\phi_A and \alpha_f is calculated using the upwind convection scheme. For further details see Numerical Schemes.

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.

vanLeerConvectionTerm Module

class fipy.terms.vanLeerConvectionTerm.VanLeerConvectionTerm(coeff=1.0, var=None)

Bases: fipy.terms.explicitUpwindConvectionTerm.ExplicitUpwindConvectionTerm

Create a _AbstractConvectionTerm object.

>>> from fipy import *
>>> m = Grid1D(nx = 2)
>>> cv = CellVariable(mesh = m)
>>> fv = FaceVariable(mesh = m)
>>> vcv = CellVariable(mesh=m, rank=1)
>>> vfv = FaceVariable(mesh=m, rank=1)
>>> __ConvectionTerm(coeff = cv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = fv) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> __ConvectionTerm(coeff = vcv)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = vfv)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.]]), mesh=UniformGrid1D(dx=1.0, nx=2)))
>>> __ConvectionTerm(coeff = (1,))
__ConvectionTerm(coeff=(1,))
>>> ExplicitUpwindConvectionTerm(coeff = (0,)).solve(var=cv, solver=DummySolver()) 
Traceback (most recent call last):
    ...
TransientTermError: The equation requires a TransientTerm with explicit convection.
>>> (TransientTerm(0.) - ExplicitUpwindConvectionTerm(coeff = (0,))).solve(var=cv, solver=DummySolver(), dt=1.)

>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = 1)).solve(var=cv, solver=DummySolver(), dt=1.) 
Traceback (most recent call last):
    ...
VectorCoeffError: The coefficient must be a vector value.
>>> m2 = Grid2D(nx=2, ny=1)
>>> cv2 = CellVariable(mesh=m2)
>>> vcv2 = CellVariable(mesh=m2, rank=1)
>>> vfv2 = FaceVariable(mesh=m2, rank=1)
>>> __ConvectionTerm(coeff=vcv2)
__ConvectionTerm(coeff=_ArithmeticCellToFaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> __ConvectionTerm(coeff=vfv2)
__ConvectionTerm(coeff=FaceVariable(value=array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]]), mesh=UniformGrid2D(dx=1.0, nx=2, dy=1.0, ny=1)))
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = ((0,),(0,)))).solve(var=cv2, solver=DummySolver(), dt=1.)
>>> (TransientTerm() - ExplicitUpwindConvectionTerm(coeff = (0,0))).solve(var=cv2, solver=DummySolver(), dt=1.)
Parameters :
  • coeff : The Term‘s coefficient value.

The National Institute of Standards and Technology (NIST) is an agency of the U.S. Department of Commerce

Privacy Policy / Security Notice / Accessibility Statement / Disclaimer / Freedom of Information Act (FOIA) /
Environmental Policy Statement / No Fear Act Policy / NIST Information Quality Standards /
Scientific Integrity Summary

Last updated: Apr 02, 2013 | Contact: Ctcms Webmaster