Adding New Governing Equations and/or Bulk Rheologies#
There are four basic tasks for adding new physics in the form of a governing equation:
Select the fields for the solution. This will control the form of the partial differential equation and the terms in the residuals and Jacobians.
Derive the pointwise functions for the residuals and Jacobians. Determine flags that will be used to indicate which terms to include.
Determine which parameters in the pointwise functions could vary in space as well as any state variables. We bundle all state variables and spatially varying parameters into a field called the auxiliary field. Each material has a separate auxiliary field.
Parameters that are spatially uniform are treated separately from the parameters in the auxiliary field.
Material is responsible for the terms in the governing equations associated with the domain (i.e., volume integrals in a 3D domain and surface integrals in a 2D domain).
A separate object implements the bulk rheology for a specific governing equation.
Fig. 144 shows the objects used to implement multiple rheologies for the elasticity equation: an isotropic, linear elastic rheology for incompressible elasticity, and an isotropic, linear elastic rheology for poroelasticity.
Elasticity object describes the physics for the elasticity equation, including the pointwise functions and flags for turning on optional terms (such as inertia) in the governing equation, and
RheologyElasticity defines the interface for bulk elastic rheologies.
Define solution subfields.
All subfields in the solution field are
SolutionSubfieldobjects (see Fig. 145). PyLith already includes several solution subfields:
SubfieldDisplacementDisplacement vector field.
SubfieldVelocityVelocity vector field.
SubfieldLagrangeFaultLagrange multiplier field for fault constraints.
SubfieldPressureFluid pressure or mean stress scalar field.
SubfieldTemperatureTemperature scalar field.
PyLith includes solution field containers with predefined subfields:
SolnDispSolution composed of a displacement field.
SolnDispVelSolution composed of displacement and velocity fields.
SolnDispPresSolution composed of displacement and mean stress (pressure) fields.
SolnDispLagrangeSolution composed of displacement and Lagrange multiplier fields.
SolnDispPresLagrangeSolution composed of displacement, mean stress (pressure), and Lagrange multiplier subfields.
SolnDispVelLagrangeSolution composed of displacement, velocity, and Lagrange multiplier subfields.
Define auxiliary subfields.
The auxiliary subfields for a governing equation are defined as facilities in a Pyre Component. For example, the ones for
AuxFieldsElasticity. The order of the subfields is defined not by the order they are listed in the Pyre component, but by the order they are added to the auxiliary field in the C++ object. The auxiliary subfields bulk rheologies are defined in the same way.
A single auxiliary field will be created for each material; it contains the auxiliary subfields from both the governing equation and the bulk rheology.
Flags to turn on/off terms in governing equation.
For the elasticity equation, we sometimes do not include body forces or inertial terms in our simulations. Rather than implement these cases as separate materials, we simply include flags in the material to turn these terms on/off. The flags are implemented as Pyre properties in our material component.
Define auxiliary subfields.
We build the auxiliary field using classes derived from
pylith::feassemble::AuxiliaryFactory. The method corresponding to each subfield specifies the name of the subfield, its components, and scale for nondimensionalizing. We generally create a single auxiliary factory object for each governing equation, but not each bulk constitutive model, because constitutive models for the same governing equation often have many of the same subfields. For example, most of our bulk constitutive models for the elasticity contain density, bulk modulus, and shear modulus auxiliary subfields.
Within the concrete implementation of the material and bulk rheology objects, we add the subfields to the auxiliary field. The order in which they are added determines the order they will be in the auxiliary field. You will need to use know this order when you implement the pointwise functions. See Fig. 146 for more information on the layout of the auxiliary field.
Implement the pointwise functions.
The pointwise functions for the residuals, Jacobians, and projections follow nearly identical interfaces. In most cases you should implement the terms in the governing equation in 3D with 2D plane strain using the more general 3D version. While this is slightly less efficient, it results in less code to maintain.
Set the pointwise functions.
We set the pointwise functions for the RHS and LHS residuals and Jacobians, taking into consideration which optional terms of the governing equation have been selected by the user.