Skip to content

Fix poroelasticity with gravity #851

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 7 commits into from
Aug 15, 2025
Merged

Fix poroelasticity with gravity #851

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
1cbeb1d
Updated PETSc to 3.23.5.
baagaard-usgs Aug 6, 2025
14cdb39
FIX: Use scale and magnitude of values to determine tolerance for ful…
baagaard-usgs Aug 12, 2025
38ab7de
FIX: Fix typo in setting context for poroelasticity with gravity.
baagaard-usgs Aug 6, 2025
df37fa4
FIX: Correct solution for full-scale test of poroelasticity with pres…
baagaard-usgs Aug 14, 2025
7e21b68
Add journal debug output of auxiliary field for elasticity and poroel…
baagaard-usgs Aug 6, 2025
74b7b58
Add full-scale test of poroelasticity with gravity and/or body forces.
baagaard-usgs Aug 6, 2025
ec91c85
Update CHANGES.
baagaard-usgs Aug 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ The version numbers are in the form `MAJOR.MINOR.PATCH`, where major releases in
* Account for processes without cells in initial condition patch when verifying configuration.
* Avoid divide by zero for `KinSrcRamp` when final slip is zero.
* Remove `solid_bulk_modulus` as a spatial database field for poroelasticity; compute it from the other fields.
* Fixed typos in setting up gravity with poroelasticity.
* Fix pythia import in `pylith_eqinfo`.

## Version 4.2.0 (2025-01-15)
Expand Down
3 changes: 2 additions & 1 deletion configure.ac
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])])

dnl PETSC
AC_LANG(C)
CIT_PATH_PETSC([3.23.4])
CIT_PATH_PETSC([3.23.5])
CIT_HEADER_PETSC
CIT_CHECK_LIB_PETSC

Expand Down Expand Up @@ -266,6 +266,7 @@ AC_CONFIG_FILES([Makefile
tests/fullscale/viscoelasticity/nofaults-2d/Makefile
tests/fullscale/viscoelasticity/nofaults-3d/Makefile
tests/fullscale/poroelasticity/Makefile
tests/fullscale/poroelasticity/nofaults-2d/Makefile
tests/fullscale/poroelasticity/faults-2d/Makefile
tests/fullscale/poroelasticity/terzaghi/Makefile
tests/fullscale/poroelasticity/mandel/Makefile
Expand Down
10 changes: 5 additions & 5 deletions libsrc/pylith/fekernels/IsotropicLinearPoroelasticity.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1181,7 +1181,7 @@ public:
pylith::fekernels::Poroelasticity::Context poroelasticContext;
pylith::fekernels::Poroelasticity::setContextQuasistatic(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
pylith::fekernels::Poroelasticity::setContextBodyForce(
pylith::fekernels::Poroelasticity::setContextGravity(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);

// Rheology Context
Expand Down Expand Up @@ -1233,7 +1233,7 @@ public:
pylith::fekernels::Poroelasticity::Context poroelasticContext;
pylith::fekernels::Poroelasticity::setContextQuasistatic(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
pylith::fekernels::Poroelasticity::setContextBodyForce(
pylith::fekernels::Poroelasticity::setContextGravity(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);

// Rheology Context
Expand Down Expand Up @@ -3287,7 +3287,7 @@ public:
pylith::fekernels::Poroelasticity::Context poroelasticContext;
pylith::fekernels::Poroelasticity::setContextQuasistatic(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
pylith::fekernels::Poroelasticity::setContextGravity(
pylith::fekernels::Poroelasticity::setContextGravityBodyForce(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);

// Rheology Context
Expand Down Expand Up @@ -3337,7 +3337,7 @@ public:
pylith::fekernels::Poroelasticity::Context poroelasticContext;
pylith::fekernels::Poroelasticity::setContextQuasistatic(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
pylith::fekernels::Poroelasticity::setContextGravity(
pylith::fekernels::Poroelasticity::setContextGravityBodyForce(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);

// Rheology Context
Expand Down Expand Up @@ -3895,7 +3895,7 @@ public:
pylith::fekernels::Poroelasticity::Context poroelasticContext;
pylith::fekernels::Poroelasticity::setContextDynamic(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
pylith::fekernels::Poroelasticity::setContextGravitySourceDensity(
pylith::fekernels::Poroelasticity::setContextGravityBodyForceSourceDensity(
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);

// Rheology Context
Expand Down
9 changes: 4 additions & 5 deletions libsrc/pylith/fekernels/Poroelasticity.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -435,8 +435,8 @@ public:
assert(context);
assert(numS >= 3);

const PylithInt i_gravityField = 4;
const PylithInt i_bodyForce = 5;
const PylithInt i_bodyForce = 4;
const PylithInt i_gravityField = 5;

assert(aOff[i_gravityField] >= 0);
assert(aOff[i_bodyForce] >= 0);
Expand Down Expand Up @@ -513,8 +513,8 @@ public:
assert(context);
assert(numS >= 3);

const PylithInt i_gravityField = 4;
const PylithInt i_bodyForce = 5;
const PylithInt i_bodyForce = 4;
const PylithInt i_gravityField = 5;
const PylithInt i_sourceDensity = 6;

assert(aOff[i_gravityField] >= 0);
Expand Down Expand Up @@ -887,7 +887,6 @@ public:
pylith::fekernels::Poroelasticity::setContextGravity(
&context, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);

// Poroelastic Auxiliariies
const PylithScalar bulkDensity = context.bulkDensity;
const PylithScalar* gravityField = context.gravityField;

Expand Down
5 changes: 5 additions & 0 deletions libsrc/pylith/materials/Elasticity.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,6 +249,11 @@ pylith::materials::Elasticity::createAuxiliaryField(const pylith::topology::Fiel
assert(auxiliaryFactory);
auxiliaryFactory->setValuesFromDB();

pythia::journal::debug_t debug("elasticity.view_auxiliary_field");
if (debug.state()) {
auxiliaryField->view("Elasticity auxiliary field");
} // if

_Elasticity::Events::logger.eventEnd(_Elasticity::Events::createAuxiliaryField);
PYLITH_METHOD_RETURN(auxiliaryField);
} // createAuxiliaryField
Expand Down
113 changes: 59 additions & 54 deletions libsrc/pylith/materials/IsotropicLinearPoroelasticity.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,47 +147,39 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelg0p(const spatialdata

switch (bitUse) {
case 0x0:
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
NULL;
break;
case 0x1:
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
NULL;
break;
case 0x2:
case 0x3:
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
NULL;
break;
case 0x4:
// aOff for sourceDensity is 3
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source :
NULL; // aOff for sourceDensity is 3
break;
case 0x3:
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
NULL;
break;
case 0x5:
// aOff for sourceDensity is 4
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source_body :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source_body :
NULL; // aOff for sourceDensity is 4
NULL;
break;
case 0x6:
// aOff for sourceDensity is 4
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source_grav :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source_grav :
NULL; // aOff for sourceDensity is 4
NULL;
break;
case 0x7:
// aOff for sourceDensity is 5
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source_grav_body :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source_grav_body :
NULL; // aOff for sourceDensity is 5
NULL;
break;
default:
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ") for Poroelasticity RHS residual kernels.");
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(g0p);
Expand Down Expand Up @@ -231,8 +223,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelg1p_explicit(const sp
NULL;
break;
default:
PYLITH_COMPONENT_ERROR("Unknown combination of flags for _useTensorPermeability="<<_useTensorPermeability<<", _gravityField="<<_gravityField<<").");
throw std::logic_error("Unknown combination of flags.");
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(g1p);
Expand Down Expand Up @@ -264,8 +255,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelg1v_explicit(const sp
NULL;
break;
default:
PYLITH_COMPONENT_ERROR("Unknown combination of flags for getKernelg1v_explicit.");
throw std::logic_error("Unknown combination of flags.");
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(g1v);
Expand Down Expand Up @@ -311,47 +301,39 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelf0p_implicit(const sp

switch (bitUse) {
case 0x0:
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
NULL;
break;
case 0x1:
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
NULL;
break;
case 0x2:
case 0x3:
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
NULL;
break;
case 0x4:
// aOff for sourceDensity is 3
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source :
NULL; // aOff for sourceDensity is 3
break;
case 0x3:
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
NULL;
break;
case 0x5:
// aOff for sourceDensity is 4
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source_body :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source_body :
NULL; // aOff for sourceDensity is 4
NULL;
break;
case 0x6:
// aOff for sourceDensity is 4
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source_grav :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source_grav :
NULL; // aOff for sourceDensity is 4
NULL;
break;
case 0x7:
// aOff for sourceDensity is 5
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source_grav_body :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source_grav_body :
NULL; // aOff for sourceDensity is 5
NULL;
break;
default:
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ") for Poroelasticity LHS residual kernels.");
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(f0p);
Expand Down Expand Up @@ -383,8 +365,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelf1u_implicit(const sp
NULL;
break;
default:
PYLITH_COMPONENT_ERROR("Unknown combination of flags for getKernelf1u_implicit.");
throw std::logic_error("Unknown combination of flags.");
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(f1u);
Expand Down Expand Up @@ -451,8 +432,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelf1p_implicit(const sp
break;

default:
PYLITH_COMPONENT_ERROR("Unknown combination of flags for _useTensorPermeability="<<_useTensorPermeability<<", _useBodyForce="<<_useBodyForce<<", _gravityField="<<_gravityField<<").");
throw std::logic_error("Unknown combination of flags.");
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(f1p);
Expand Down Expand Up @@ -535,12 +515,24 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelJf3pp(const spatialda
PYLITH_COMPONENT_DEBUG("getKernelJf3pp(coordsys="<<typeid(coordsys).name()<<")");

const int spaceDim = coordsys->getSpaceDim();
PetscPointJacFn* Jf3pp =
(!_useTensorPermeability && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp :
(!_useTensorPermeability && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp :
(_useTensorPermeability && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp_tensor_permeability :
(_useTensorPermeability && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp_tensor_permeability :
NULL;
const int bitTensorPermeability = _useTensorPermeability ? 0x1 : 0x0;
const int bitUse = bitTensorPermeability;

PetscPointJacFn* Jf3pp = NULL;
switch (bitUse) {
case 0x0:
Jf3pp = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp :
NULL;
break;
case 0x1:
Jf3pp = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp_tensor_permeability :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp_tensor_permeability :
NULL;
break;
default:
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(Jf3pp);
} // getKerneJf3pp
Expand Down Expand Up @@ -573,12 +565,25 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelCauchyStressVector(co
PYLITH_COMPONENT_DEBUG("getKernelCauchyStressVector(coordsys="<<typeid(coordsys).name()<<")");

const int spaceDim = coordsys->getSpaceDim();
PetscPointFn* kernel =
(!_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_asVector :
(!_useReferenceState && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_asVector :
(_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_refState_asVector :
(_useReferenceState && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_refState_asVector :
NULL;
const int bitReferenceState = _useReferenceState ? 0x1 : 0x0;
const int bitUse = bitReferenceState;

PetscPointFn* kernel = NULL;

switch (bitUse) {
case 0x0:
kernel = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_asVector :
NULL;
break;
case 0x1:
kernel = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_refState_asVector :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_refState_asVector :
NULL;
break;
default:
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
} // switch

PYLITH_METHOD_RETURN(kernel);
} // getKernelCauchyStressVector
Expand Down
19 changes: 6 additions & 13 deletions libsrc/pylith/materials/Poroelasticity.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,11 @@ pylith::materials::Poroelasticity::createAuxiliaryField(const pylith::topology::
assert(auxiliaryFactory);
auxiliaryFactory->setValuesFromDB();

pythia::journal::debug_t debug("poroelasticity.view_auxiliary_field");
if (debug.state()) {
auxiliaryField->view("Poroelasticity auxiliary field");
} // if

PYLITH_METHOD_RETURN(auxiliaryField);
} // createAuxiliaryField

Expand Down Expand Up @@ -395,8 +400,7 @@ pylith::materials::Poroelasticity::_setKernelsResidual(pylith::feassemble::Integ

const int bitBodyForce = _useBodyForce ? 0x1 : 0x0;
const int bitGravity = _gravityField ? 0x2 : 0x0;
const int bitSourceDensity = _useSourceDensity ? 0x4 : 0x0;
const int bitUse = bitBodyForce | bitGravity | bitSourceDensity;
const int bitUse = bitBodyForce | bitGravity;

PetscPointFn* r0 = NULL;
switch (bitUse) {
Expand All @@ -408,20 +412,9 @@ pylith::materials::Poroelasticity::_setKernelsResidual(pylith::feassemble::Integ
case 0x2:
r0 = pylith::fekernels::Poroelasticity::g0v_grav;
break;
case 0x4:
break;
case 0x3:
r0 = pylith::fekernels::Poroelasticity::g0v_grav_bodyforce;
break;
case 0x5:
r0 = pylith::fekernels::Poroelasticity::g0v_bodyforce;
break;
case 0x6:
r0 = pylith::fekernels::Poroelasticity::g0v_grav;
break;
case 0x7:
r0 = pylith::fekernels::Poroelasticity::g0v_grav_bodyforce;
break;
default:
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ") for residual kernels.");
} // switch
Expand Down
Loading