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

Class LinearPDE

   object --+    
            |    
LinearProblem --+
                |
               LinearPDE
Known Subclasses:

This class is used to define a general linear, steady, second order PDE for an unknown function u on a given domain defined through a Domain object.

For a single PDE having a solution with a single component the linear PDE is defined in the following form:

-(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 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 A, B, C, D, X and Y have to be specified through Data objects in Function and the coefficients 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. 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)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y

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 PDE. The coefficients d and y are each a scalar in FunctionOnBoundary and the coefficients d_reduced and y_reduced are each a scalar in ReducedFunctionOnBoundary.

Constraints for the solution prescribe the value of the solution at certain locations in the domain. They have the form

u=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 PDE or the boundary condition.

The PDE 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 of PDEs and a solution with several components the PDE has the form

-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, 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])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]

The coefficient d is of rank two and y is of rank one both in FunctionOnBoundary. The coefficients d_reduced is of rank two and y_reduced is of rank one both in ReducedFunctionOnBoundary.

Constraints take the form

u[i]=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 system of PDEs is symmetrical if

LinearPDE 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 PDEs 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 = LinearPDE(dom)
   p.setValue(A=kronecker(dom), D=1, Y=0.5)
   u = p.getSolution()
Instance Methods [hide private]
 
__init__(self, domain, numEquations=None, numSolutions=None, debug=False)
Initializes a new linear PDE.
str
__str__(self)
Returns the string representation of the PDE.
bool
checkSymmetry(self, verbose=True)
Tests the PDE for symmetry.
 
createOperator(self)
Returns an instance of a new operator.
Data
getFlux(self, u=None)
Returns the flux J for a given u.
 
getRequiredOperatorType(self)
Returns the system type which needs to be used by the current set up.
Data
getResidual(self, u=None)
Returns the residual of u or the current solution if u is not present.
Data
getSolution(self)
Returns the solution of the PDE.
tuple of Operator, and Data
getSystem(self)
Returns the operator and right hand side of the PDE.
 
insertConstraint(self, rhs_only=False)
Applies the constraints defined by q and r to the PDE.
 
setValue(self, **coefficients)
Sets new values to coefficients.

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, setDebug, setDebugOff, setDebugOn, 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, debug=False)
(Constructor)

 

Initializes a new linear PDE.

Parameters:
  • domain (Domain) - domain of the PDE
  • numEquations - number of equations. If None the number of equations is extracted from the PDE coefficients.
  • numSolutions - number of solution components. If None the number of solution components is extracted from the PDE coefficients.
  • debug - if True debug information is printed
Overrides: object.__init__

__str__(self)
(Informal representation operator)

 

Returns the string representation of the PDE.

Returns: str
a simple representation of the PDE
Overrides: object.__str__

checkSymmetry(self, verbose=True)

 

Tests the PDE 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 operator.

Overrides: LinearProblem.createOperator

getFlux(self, u=None)

 

Returns the flux J for a given u.

J[i,j]=(A[i,j,k,l]+A_reduced[A[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]

or

J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]

Parameters:
  • u (Data or None) - argument in the flux. If u is not present or equals None the current solution is used.
Returns: Data
flux

getRequiredOperatorType(self)

 

Returns the system type which needs to be used by the current set up.

Overrides: LinearProblem.getRequiredOperatorType

getResidual(self, u=None)

 

Returns the residual of u or the current solution if u is not present.

Parameters:
Returns: Data
residual of u

getSolution(self)

 

Returns the solution of the PDE.

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

insertConstraint(self, rhs_only=False)

 

Applies the constraints defined by q and r to the PDE.

Parameters:
  • rhs_only (bool) - if True only the right hand side is altered by the constraint

setValue(self, **coefficients)

 

Sets new values to coefficients.

Parameters:
Raises:
Overrides: LinearProblem.setValue