proteus.Transport module

This module contains methods for solving

\[m_t^i + \nabla \cdot (f^i - \sum_k a^{ik} \nabla \phi^k) + H^i(\nabla u) + r^i = 0 \quad i=0,\ldots,nc-1\]

The solution is the vector of components \(u=u_0,\ldots,u_{nc-1}\), and the nonlinear coefficients are \(m^i,f^i,a^{ik},\phi^k, H^i\) and \(r^i\).

class proteus.Transport.StorageSet(initializer=[], shape=(0,), storageType='d')

Bases: set

allocate(storageDict)

class proteus.Transport.OneLevelTransport(uDict, phiDict, testSpaceDict, matType, dofBoundaryConditionsDict, dofBoundaryConditionsSetterDict, coefficients, elementQuadrature, elementBoundaryQuadrature, fluxBoundaryConditionsDict=None, advectiveFluxBoundaryConditionsSetterDict=None, diffusiveFluxBoundaryConditionsSetterDictDict=None, stressFluxBoundaryConditionsSetterDict=None, stabilization=None, shockCapturing=None, conservativeFluxDict=None, numericalFluxType=None, TimeIntegrationClass=None, massLumping=False, reactionLumping=False, options=None, name='defaultName', reuse_trial_and_test_quadrature=False, sd=True, movingDomain=False, bdyNullSpace=False, hasCutCells=False)

Bases: proteus.NonlinearSolvers.NonlinearEquation

A class for finite element discretizations of multicomponent advective-diffusive-reactive transport on a single spatial mesh.

Objects of this type take the initial-boundary value problem for

\[m_t^i + \nabla \cdot (\mathbf{f}^i - \sum_k \mathbf{a}^{ik} \nabla \phi^k) + H^i(\nabla u) + r^i = 0\]

and turn it into the discrete (nonlinear or linear) algebraic problem

\[F(U) = 0\]

where F and U have dimension self.dim. The Jacobian of F or an approximation for it may also be provided.

The NonlinearEquation interface is

• self.dim

• getResidual(u,r)

• getJacobian(jacobian).

The rest of the functions in this class are either private functions or return various other pieces of information.

Variables

ebq_global[('velocityAverage',0)] (array) – This attribute stores the average velocity along an edge given a discontinous velocity field.

Allocate storage and initialize some variables.

Parameters

• uDict (dict) – Dictionary of proteus.FemTools.FiniteElementFunction objects.

• phiDict (dict) – Dictionary of proteus.FemTools.FiniteElementFunction objects.

• testSpaceDict (dict) – Dictionary of FiniteElementSpace objects

• dofBoundaryConditionsDict (dict) – Dictionary of DOFBoundaryConditions objects for the Dirichlet conditions.

• coefficients (proteus.TransportCoefficients.TC_base) – Problem’s Transport Coefficients class.

• elementQuadratureDict (dict) – Dictionary of dictionaries of quadrature rules for each element integral in each component equation.

• elementBoundaryQuadratureDict (dict) – Dictionary of dictionaries of quadrature rules for each element boundary integral in each component equation

• stabilization (bool) –

• shockCapturing (bool) –

• numericalFlux (bool) –

• bdyNullSpace (bool) – Indicates whether the boundary conditions create a global null space.

Notes

The constructor sets the input arguments, calculates dimensions, and allocates storage. The meanings of variable suffixes are

• global – per physical domain

• element – per element

• elementBoundary – per element boundary

The prefix n means ‘number of’.

Storage is divided into quantities required at different sets of points or geometric entities. Each type of storage has a dictionary for all the quantities of that type. The names and dimensions of the storage dictionaries are

• e – at element

• q – at element quadrature, unique to elements

• ebq – at element boundary quadrature, unique to elements

• ebq_global – at element boundary quadrature, unique to element boundary

• ebqe – at element boundary quadrature, unique to global, exterior element boundary

• phi_ip – at the generalized interpolation points required to build a nonlinear phi

setupFieldStrides(interleaved=True)

Set up the stride/offset layout for this object. The default layout is interleaved.

setInitialConditions(getInitialConditionsDict, T=0.0)

setInitialConditionsTPF(getInitialConditionsDict, idxDict, T=0.0)

archiveAnalyticalSolutions(archive, analyticalSolutionsDict, T=0.0, tCount=0)

setInitialConditionsFromDofs(getInitialConditionDofs, T=0.0)

initializeTimeIntegration(t0, tOut, u0, r0, DTSET=None)

initializeTimeHistory()

updateTimeHistory(T, resetFromDOF=False)

put a version on each level rather than just multilevel?

calculateAuxiliaryQuantitiesAfterStep()

getResidual(u, r)

Calculate the element residuals and add in to the global residual

getJacobian(jacobian, skipMassTerms=False)

getJacobian_dense(jacobian)

getJacobian_CSR(jacobian)

Add in the element jacobians to the global jacobian

calculateElementResidual()

Calculate all the element residuals

estimate_mt()

calculateElementJacobian(skipMassTerms=False)

calculateElementMassJacobian()

calculate just the mass matrix terms for element jacobian (i.e., those that multiply the accumulation term) does not include dt

calculateElementBoundaryJacobian()

calculateExteriorElementBoundaryJacobian()

calculateCoefficients()

calculateSolutionAtQuadrature()

calculateStrongResidualAndAdjoint(cq)

break off computation of strong residual and adjoint so that we can do this at different quadrature locations?

calculateElementCoefficients()

calculate the nonlinear coefficients at the quadrature points and nodes

calculateElementBoundaryCoefficients()

Calculate the nonlinear coefficients at the element boundary quadrature points

calculateExteriorElementBoundaryCoefficients()

Calculate the nonlinear coefficients at global exterior element boundary quadrature points

calculateQuadrature()

updateAfterMeshMotion()

calculateElementQuadrature()

Calculate the physical location and weights of the quadrature rules and the shape information at the quadrature points.

This function should be called only when the mesh changes.

calculateElementBoundaryQuadrature()

Calculate the physical location and weights of the quadrature rules and the shape information at the quadrature points on element boundaries.

This function should be called only when the mesh changes.

calculateExteriorElementBoundaryQuadrature()

Calculate the physical location and weights of the quadrature rules and the shape information at the quadrature points on global element boundaries.

This function should be called only when the mesh changes.

setFreeDOF(free_u)

Set the free(non-Dirichlet) DOF from the global DOF

updateLocal2Global()

Build a mapping between local DOF numbers and global free DOF numbers, that is, EXCLUDING Dirichlet DOF. We have to use a list of local dof because of the exclusion of Dirichlet nodes (otherwise we could just loop over range(self.nDOF_element).

setInflowFlux()

setUnknowns(free_u)

initializeJacobian()

viewSolution(plotOffSet=None, titleModifier='', dgridnx=50, dgridny=50, dgridp=16.0, pause=False)

viewSolutionVTK(plotOffSet=None, titleModifier='', dgridnx=50, dgridny=50, dgridp=16.0, pause=False)

saveSolution()

write_ebq_geo_Ensight(filename)

write_ebq_velocity_Ensight(filename, nOutput, append=False, firstVariable=True, case_filename=None)

archiveFiniteElementSolutions(archive, t, tCount, initialPhase=False, writeVectors=True, meshChanged=False, femSpaceWritten={}, writeVelocityPostProcessor=True)

write finite element solutions to archive at time t, tries to group finite element functions by their space

archiveElementQuadratureValues(archive, t, tCount, scalarKeys=None, vectorKeys=None, tensorKeys=None, initialPhase=False, meshChanged=False)

write element quadrature dictionary values to archive at time t, groups all entries under same mesh at a given time level

archiveElementBoundaryQuadratureValues(archive, t, tCount, scalarKeys=None, vectorKeys=None, tensorKeys=None, initialPhase=False, meshChanged=False)

write elementBoundary quadrature dictionary values to archive at time t, groups all entries under same mesh at a given time level

archiveExteriorElementBoundaryQuadratureValues(archive, t, tCount, scalarKeys=None, vectorKeys=None, tensorKeys=None, initialPhase=False, meshChanged=False)

write exteriorElementBoundary quadrature dictionary values to archive at time t, groups all entries under same mesh at a given time level

archiveFiniteElementResiduals(archive, t, tCount, res_dict, res_name_base='spatial_residual')

write assembled spatial residual stored in r at time t

ASSUMES archiveFiniteElementSolutions has already been called for t and tCount!!!

initializeMassJacobian()

Setup the storage for the mass jacobian and return as a `SparseMat` or `Mat` based on self.matType

initializeSpatialJacobian()

Setup the storage for the spatial jacobian and return as a `SparseMat` or `Mat` based on self.matType

calculateElementLoadCoefficients_inhomogeneous()

Calculate the non-solution dependent portion of the Reaction term, ‘r’ Temporary fix for model reduction for linear problems Just Zero’s the solution and calls the usual update

calculateElementLoad_inhomogeneous()

Calculate all the portion of the weak element residual corresponding to terms that the ‘inhomogeneous’ or constant portion of the traditional load vector. This includes purely spatio-temporal portions of the reaction term, ‘r’, and boundary condition terms This is a temporary fix for linear model reduction.

calculateExteriorElementBoundaryCoefficients_inhomogeneous()

Calculate the coefficients at global exterior element boundary quadrature points for terms that are independent of the solution, i.e., space and time dependent portions of boundary integrals Just sets the solution to zero and evaluates the boundary integrals This is a temporary fix for linear model reduction

getLoadVector(f)

Return the non-solution dependent portion of the Reaction term, ‘r’ and flux boundary conditions

getMassJacobian(jacobian)

assemble the portion of the jacobian coming from the time derivative terms

getSpatialJacobian(jacobian)

getSpatialResidual(u, r)

Calculate the element spatial residuals and add in to the global residual This is being used right now to test nonlinear pod and deim

getMassResidual(u, r)

Calculate the portion of the element residuals associated with the temporal term and assemble into a global residual This is being used right now to test nonlinear pod and deim and is inefficient

attachLaplaceOperator(nu=1.0)

Attach a Discrete Laplace Operator to the Transport Class.

Parameters

nu (float) – Viscosity parameter for the Laplace operator.

Notes

This function creates a Laplace matrix and stores the result in the proteus.OneLevelTransport class to assist in the construction of physics based preconditione

class proteus.Transport.MultilevelTransport(problem, numerics, mlMesh, OneLevelTransportType=<class 'proteus.Transport.OneLevelTransport'>)

Bases: object

Nonlinear ADR on a multilevel mesh

initialize(nd, mlMesh, TrialSpaceTypeDict, TestSpaceTypeDict, matType, dirichletConditionsSetterDict, coefficients, quadratureDict, elementBoundaryQuadratureDict, fluxBoundaryConditionsDict='noFlow', advectiveFluxBoundaryConditionsSetterDict={}, diffusiveFluxBoundaryConditionsSetterDictDict={}, stressFluxBoundaryConditionsSetterDict={}, stabilization=None, shockCapturing=None, conservativeFlux=None, numericalFluxType=None, TimeIntegrationClass=<class 'proteus.TimeIntegration.BackwardEuler'>, massLumping=False, reactionLumping=False, weakDirichletConditions=None, options=None, useSparseDiffusion=True, movingDomain=False, PhiSpaceTypeDict=None)

setInitialConditions(getInitialConditionsDict, T=0.0)

setInitialConditionsFromDofs(T)

calculateAuxiliaryQuantitiesAfterStep()

attachModels(modelList)

viewJacobian(file_prefix='./dump_', b=None)

Save parallel Jacobian list and (optionally) right-hand side b to disk if they possess save methods, otherwise, do nothing.