Reimplement fault constitutive models using new formulation

[baagaard-usgs]

FaultCohesiveDyn

• fault_rheology
• traction_perturbation
• open_free_surface

Friction ISA FaultRheology

• FaultRheology implements T(x,t) without assuming friction coefficient
• Friction implements T_f = C + \mu_f T_n
• Friction models implement computing \mu(f) as a function of slip, slip rate, and state variables

Friction Models

• StaticFriction
• SlipWeakening

TractionPerturbation [DONE]

Implements T_initial with same T(x,t) as boundary conditions

Issues

• How do we set the direction of friction in 3D?
  • If sliding, then we can use the current direction of slip rate (Lagrange multiplier is 0).
  • If not sliding, then friction should be opposite of Lagrange multiplier (need to account for initial traction).

[baagaard-usgs]

plith::faults

FaultCohesiveDyn

setKernelsResidual()

• f0u_neg = faultRheology->getF0uNegKernel()
• f0u_pos = faultRheology->getF0uPosKernel()
• f0l kernel that enforces traction * slip_opening = 0; traction only = 0 if slip opening > 0

FrictionStatic

• Implements getF0NegKernel() [returns pylith::fekernels::FrictionStatic::f0u_neg]
• Implements getF0PosKernel() [returns pylith::fekernels::FrictionStatic::f0u_pos]

pylith::fekernels

FaultFriction

Context

• normal (argument n) and reference directions (from constants)
• tractions in global coordinates (transform them to local in local function)
• dispN (disp on negative side of the fault) (spaceDim)
• dispP (disp on positive side fo teh fault) (spaceDim)
• velN (disp on negative side of the fault) (spaceDim)
• velP (disp on positive side fo teh fault) (spaceDim)
• cohesion
• openFree surface flag
• dim
• t

• General kernels f0u (separate kernels for 2D and 3D)
  • Calls frictionCoefficientFn() passing contexts as arguments
  • Use contexts to compute lambda^* = traction_friction - lambda (all components in global Cartesian coordinate system)
  • traction_friction = C + friction_coefficient * traction_normal (makes sense in fault coordinate system; shear only)
  • Transforming between global and fault coordinate systems - use normal passed in and reference directions in constants -> BoundaryDirections to get tangential directions (similar to transforming slip to XYZ)
• FrictionContext: cohesion, normalTraction(scalar), lagrangeMultiplier(dim)

FrictionStatic

• Context - staticCoefficient
• Has f0u_neg and f0u_pos "entry" function that calls FaultFriction general kernels
  • Create context and then call FaultFriction::f0u(), passing its own FrictionStatic context and any other general contexts it created up front from auxiliary, solution fields, constants
  • Analogous to IsotropicLinearElasticity::f1v_infinitessimalStrain; instead of calling Elasticity::f1v we will call FaultFriction::f0u
• Implements a function that computes coefficient of friction using the FrictionStatic context passed down by FaultFriction::f0u(); analogous to cauchyStress() function in IsotropicLinearElasticity
  • Arguments: scalar t, slip, slipRate; FrictionStatic context

[baagaard-usgs]

We probably want to consider a zero tolerance (with a potentially separate tolerance for the fault normal direction) to detect "zero" slip and slip rate.

As in v2.x, these would be set in faults::FaultCohesiveDyn at the Python level and passed down to C++ with reasonable default values. These would be passed to the kernels in constants and set via _updateKernelConstants().

[baagaard-usgs]

@knepley The residual kernels for the contributions to the displacement field for the fault faces (which is the Lagrange multiplier) only get the solution for one side of the fault. For example, pylith::fekernels::FaultCohesiveKin::f0u_pos() gets passed in sOff with sOff[0] == 2. sOff is being set by PetscDSGetComponentOffsetsCohesive() which takes the side of the fault as an argument.

For fault friction we want to compute slip and slip rate, which are parameters to the friction model, to compute the traction term which gets combined with the Lagrange multiplier. I think we want to pass in the solution for both sides of the fault for all kernels involving the cohesive cell, even though the kernel "output" is only for one side in the case of the functions for the fault faces.

[knepley]

@baagaard-usgs Oh, that is right. We wanted to mimic a regular kernel. If we now want it to take all values, we can just pass 0 instead of s to the OffsetCohesive(), or even have a switch to indicate what kind of kernel input you want.

[baagaard-usgs]

We decided to change the PETSc code so that we always get all solution fields for both sides of the fault when integrating the fault faces or cohesive cells.

Matt has created a PETSc branch knepley/fix-cohesive-offsets that updates the solution field offsets so that we will get solution fields for both sides of the fault in the kernels involving cohesive cells (fault faces and cohesive cell). I am working on a PyLith branch update-fault-kernels-solution-offsets in my fork that updates the fault kernels accordingly. We hope to get these branches integrated into PETSc knepley/pylith and PyLIth main by the end of this week (Aug 4).