Package esys :: Package escript :: Module linearPDEs :: Class TransportPDE
[hide private]
[frames] | no frames]

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

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)
Instance Methods [hide private]
 
__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.
 
createOperator(self)
Returns an instance of a new transport operator.
float
getConstraintWeightingFactor(self)
returns the weighting factor used to insert the constraints into the problem
float
getRequiredOperatorType(self)
Returns the system type which needs to be used by the current set up.
float
getSafeTimeStepSize(self)
Returns a safe time step size to do the next time step.
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.
 
setConstraintWeightingFactor(self, value=67108864.0)
Sets the weighting factor used to insert the constraints into the problem
 
setDebug(self, flag)
Switches debug output on if flag is True, otherwise it is switched off.
 
setDebugOff(self)
Switches debug output off.
 
setDebugOn(self)
Switches debug output on.
 
setInitialSolution(self, u)
Sets the initial solution.
 
setValue(self, **coefficients)
Sets new values to coefficients.
bool
useBackwardEuler(self)
Returns true if backward Euler is used.

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__

Properties [hide private]

Inherited from object: __class__

Method Details [hide private]

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

createOperator(self)

 

Returns an instance of a new transport operator.

Overrides: LinearProblem.createOperator

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.

getSolution(self, dt)

 

Returns the solution of the problem.

Returns: Data
the solution
Overrides: LinearProblem.getSolution

getSystem(self)

 

Returns the operator and right hand side of the PDE.

Returns: tuple of Operator, and Data
the discrete version of the PDE
Overrides: LinearProblem.getSystem

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

setDebug(self, flag)

 

Switches debug output on if flag is True, otherwise it is switched off.

Parameters:
  • flag (bool) - desired debug status
Overrides: LinearProblem.setDebug

setDebugOff(self)

 

Switches debug output off.

Overrides: LinearProblem.setDebugOff

setDebugOn(self)

 

Switches debug output on.

Overrides: LinearProblem.setDebugOn

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

useBackwardEuler(self)

 

Returns true if backward Euler is used. Otherwise false is returned.

Returns: bool