FiPyfipy.terms.advectionTerm¶
Classes
|
The AdvectionTerm object constructs the b vector contribution for the advection term given by |
- class fipy.terms.advectionTerm.AdvectionTerm(coeff=None)¶
Bases:
FirstOrderAdvectionTermThe 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
bshould 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 (
floatorCellVariableorFaceVariable) – 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 (
tupleofBoundaryCondition) –dt (
float) – Timestep size.underRelaxation (
float) – Usually a value between 0 and 1 or None in the case of no under-relaxationresidualFn (
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:
- 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 (
tupleofBoundaryCondition) –dt (
float) – Timestep size.underRelaxation (
float) – Usually a value between 0 and 1 or None in the case of no under-relaxationresidualFn (
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:
- 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 (
tupleofBoundaryCondition) –dt (
float) – Timestep size.underRelaxation (
float) – Usually a value between 0 and 1 or None in the case of no under-relaxationresidualFn (
function) – Takes var, matrix, and RHSvector arguments, used to customize the residual calculation.
- Returns:
residual (
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 (
tupleofBoundaryCondition) –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 (
tupleofBoundaryCondition) –dt (
float) – Timestep size.underRelaxation (
float) – Usually a value between 0 and 1 or None in the case of no under-relaxationresidualFn (
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 TermcacheError (
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: