Class TransportPDE
object --+
|
LinearProblem --+
|
TransportPDE
- Known Subclasses:
-
This class is used to define a transport problem given by a general
linear, time dependent, second order PDE for an unknown, non-negative,
time-dependent function u on a given domain defined
through a Domain object.
For a single equation with a solution with a single component the
transport problem is defined in the following form:
(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)
where u_t denotes the time derivative of u and grad(F) denotes the spatial
derivative of F. Einstein's summation convention,
ie. summation over indexes appearing twice in a term of a sum performed,
is used. The coefficients M, A,
B, C, D, X and Y have to be specified through
Data
objects in Function and the coefficients M_reduced, A_reduced, B_reduced, C_reduced, D_reduced, X_reduced and Y_reduced have to be specified through Data
objects in ReducedFunction. It is also allowed to use objects that
can be converted into such Data
objects. M and M_reduced are
scalar, A and A_reduced are rank
two, B, C, X,
B_reduced, C_reduced and X_reduced are rank one and D, D_reduced, Y and Y_reduced are scalar.
The following natural boundary conditions are considered:
n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t
where n is the outer normal field. Notice that the
coefficients A, A_reduced, B, B_reduced, X
and X_reduced are defined in the transport problem.
The coefficients m, d and y are each a scalar in FunctionOnBoundary and the coefficients m_reduced, d_reduced and y_reduced are each a scalar in ReducedFunctionOnBoundary.
Constraints for the solution prescribing the value of the solution at
certain locations in the domain have the form
u_t=r where q>0
r and q are each scalar where
q is the characteristic function defining where the
constraint is applied. The constraints override any other condition set
by the transport problem or the boundary condition.
The transport problem is symmetrical if
A[i,j]=A[j,i] and B[j]=C[j]
and A_reduced[i,j]=A_reduced[j,i] and B_reduced[j]=C_reduced[j]
For a system and a solution with several components the transport
problem has the form
(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i]
A and A_reduced are of rank
four, B, B_reduced, C and C_reduced are each of rank
three, M, M_reduced, D, D_reduced, X_reduced and X are each of rank two
and Y and Y_reduced are of rank
one. The natural boundary conditions take the form:
n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]=
(m[i,k]+m_reduced[i,k])*u[k]_t
The coefficient d and m are of
rank two and y is of rank one with all in FunctionOnBoundary. The coefficients d_reduced and m_reduced are of rank
two and y_reduced is of rank one all in ReducedFunctionOnBoundary.
Constraints take the form
u[i]_t=r[i] where q[i]>0
r and q are each rank one.
Notice that at some locations not necessarily all components must have a
constraint.
The transport problem is symmetrical if
-
M[i,k]=M[i,k]
-
M_reduced[i,k]=M_reduced[i,k]
-
A[i,j,k,l]=A[k,l,i,j]
-
A_reduced[i,j,k,l]=A_reduced[k,l,i,j]
-
B[i,j,k]=C[k,i,j]
-
B_reduced[i,j,k]=C_reduced[k,i,j]
-
D[i,k]=D[i,k]
-
D_reduced[i,k]=D_reduced[i,k]
-
m[i,k]=m[k,i]
-
m_reduced[i,k]=m_reduced[k,i]
-
d[i,k]=d[k,i]
-
d_reduced[i,k]=d_reduced[k,i]
TransportPDE also supports solution discontinuities over
a contact region in the domain. To specify the conditions across the
discontinuity we are using the generalised flux J
which, in the case of a system of PDEs and several components of the
solution, is defined as
J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]
For the case of single solution component and single PDE J is defined as
J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]
In the context of discontinuities n denotes the
normal on the discontinuity pointing from side 0 towards side 1
calculated from getNormal of FunctionOnContactZero. For a system of transport
problems the contact condition takes the form
n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])-
(d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]
where J0 and J1 are the fluxes
on side 0 and side 1 of the discontinuity, respectively. jump(u), which is the difference of the solution at side
1 and at side 0, denotes the jump of u across
discontinuity along the normal calculated by jump. The
coefficient d_contact is of rank two and y_contact is of rank one both in FunctionOnContactZero or FunctionOnContactOne. The coefficient d_contact_reduced is of rank two and y_contact_reduced is of rank one both in ReducedFunctionOnContactZero or ReducedFunctionOnContactOne. In case of a single PDE and
a single component solution the contact condition takes the form
n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)
In this case the coefficient d_contact and y_contact are each scalar both in FunctionOnContactZero or FunctionOnContactOne and the coefficient d_contact_reduced and y_contact_reduced are each scalar both in ReducedFunctionOnContactZero or ReducedFunctionOnContactOne.
Typical usage:
p = TransportPDE(dom)
p.setValue(M=1., C=[-1.,0.])
p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
t = 0
dt = 0.1
while (t < 1.):
u = p.solve(dt)
|
__init__(self,
domain,
numEquations=None,
numSolutions=None,
useBackwardEuler=False,
debug=False)
Initializes a transport problem. |
|
|
str
|
__str__(self)
Returns the string representation of the transport problem. |
|
|
bool
|
checkSymmetry(self,
verbose=True)
Tests the transport problem for symmetry. |
|
|
|
|
float
|
|
float
|
|
float
|
|
Data
|
getSolution(self,
dt)
Returns the solution of the problem. |
|
|
tuple of Operator, and Data
|
getSystem(self)
Returns the operator and right hand side of the PDE. |
|
|
float
|
getUnlimitedTimeStepSize(self)
Returns the value returned by the getSafeTimeStepSize
method to indicate no limit on the safe time step size. |
|
|
|
|
|
setDebug(self,
flag)
Switches debug output on if flag is True, otherwise it
is switched off. |
|
|
|
|
|
|
|
|
|
setValue(self,
**coefficients)
Sets new values to coefficients. |
|
|
bool
|
|
Inherited from LinearProblem :
alteredCoefficient ,
checkReciprocalSymmetry ,
checkSymmetricTensor ,
createCoefficient ,
createRightHandSide ,
createSolution ,
getCoefficient ,
getCurrentOperator ,
getCurrentRightHandSide ,
getCurrentSolution ,
getDim ,
getDomain ,
getDomainStatus ,
getFunctionSpaceForCoefficient ,
getFunctionSpaceForEquation ,
getFunctionSpaceForSolution ,
getNumEquations ,
getNumSolutions ,
getOperator ,
getOperatorType ,
getRightHandSide ,
getShapeOfCoefficient ,
getSolverOptions ,
getSystemStatus ,
hasCoefficient ,
initializeSystem ,
introduceCoefficients ,
invalidateOperator ,
invalidateRightHandSide ,
invalidateSolution ,
invalidateSystem ,
isOperatorValid ,
isRightHandSideValid ,
isSolutionValid ,
isSymmetric ,
isSystemValid ,
isUsingLumping ,
reduceEquationOrder ,
reduceSolutionOrder ,
resetAllCoefficients ,
resetOperator ,
resetRightHandSide ,
resetSolution ,
setReducedOrderForEquationOff ,
setReducedOrderForEquationOn ,
setReducedOrderForEquationTo ,
setReducedOrderForSolutionOff ,
setReducedOrderForSolutionOn ,
setReducedOrderForSolutionTo ,
setReducedOrderOff ,
setReducedOrderOn ,
setReducedOrderTo ,
setSolution ,
setSolverOptions ,
setSymmetry ,
setSymmetryOff ,
setSymmetryOn ,
setSystemStatus ,
trace ,
validOperator ,
validRightHandSide ,
validSolution
Inherited from object :
__delattr__ ,
__getattribute__ ,
__hash__ ,
__new__ ,
__reduce__ ,
__reduce_ex__ ,
__repr__ ,
__setattr__
|
Inherited from object :
__class__
|
__init__(self,
domain,
numEquations=None,
numSolutions=None,
useBackwardEuler=False,
debug=False)
(Constructor)
|
|
Initializes a transport problem.
- Parameters:
domain (Domain) - domain of the PDE
numEquations - number of equations. If None the number of equations
is extracted from the coefficients.
numSolutions - number of solution components. If None the number of
solution components is extracted from the coefficients.
debug - if True debug information is printed
useBackwardEuler (bool ) - if set the backward Euler scheme is used. Otherwise the
Crank-Nicholson scheme is applied. Note that backward Euler
scheme will return a safe time step size which is practically
infinity as the scheme is unconditional unstable. The
Crank-Nicholson scheme provides a higher accuracy but requires to
limit the time step size to be stable.
- Overrides:
object.__init__
|
__str__(self)
(Informal representation operator)
|
|
Returns the string representation of the transport problem.
- Returns:
str
- a simple representation of the transport problem
- Overrides:
object.__str__
|
checkSymmetry(self,
verbose=True)
|
|
Tests the transport problem for symmetry.
- Parameters:
verbose (bool ) - if set to True or not present a report on coefficients which
break the symmetry is printed.
- Returns:
bool
- True if the PDE is symmetric
- Overrides:
LinearProblem.checkSymmetry
Note:
This is a very expensive operation. It should be used for degugging
only! The symmetry flag is not altered.
|
getConstraintWeightingFactor(self)
|
|
returns the weighting factor used to insert the constraints into the
problem
- Returns:
float
- value for the weighting factor
|
getRequiredOperatorType(self)
|
|
Returns the system type which needs to be used by the current set
up.
- Returns:
float
- a code to indicate the type of transport problem scheme used
- Overrides:
LinearProblem.getRequiredOperatorType
|
getSafeTimeStepSize(self)
|
|
Returns a safe time step size to do the next time step.
- Returns:
float
- safe time step size
Note:
If not getSafeTimeStepSize() <
getUnlimitedTimeStepSize() any time step size can be
used.
|
getUnlimitedTimeStepSize(self)
|
|
Returns the value returned by the getSafeTimeStepSize
method to indicate no limit on the safe time step size.
- Returns:
float
- the value used to indicate that no limit is set to the time step
size
Note:
Typically the biggest positive float is returned
|
setConstraintWeightingFactor(self,
value=67108864.0)
|
|
Sets the weighting factor used to insert the constraints into the
problem
- Parameters:
value (large positive float ) - value for the weighting factor
|
Switches debug output on if flag is True, otherwise it is
switched off.
- Parameters:
flag (bool ) - desired debug status
- Overrides:
LinearProblem.setDebug
|
setInitialSolution(self,
u)
|
|
Sets the initial solution.
- Parameters:
Note:
u must be non-negative
|
setValue(self,
**coefficients)
|
|
Sets new values to coefficients.
- Parameters:
- Raises:
- Overrides:
LinearProblem.setValue
|
Returns true if backward Euler is used. Otherwise false is
returned.
- Returns: bool
|