Class TransportPDE

   object --+    
            |    
LinearProblem --+
                |
               TransportPDE

Known Subclasses:

• heat.TemperatureCartesian

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)