Skip to content

Commit 28eae7a

Browse files
baagaard-usgsbaagaard-usgs
authored
Merge pull request #851 from baagaard-usgs/fix-poroelasticity-gravity
Fix poroelasticity with gravity
2 parents 843c054 + ec91c85 commit 28eae7a

37 files changed

+3026
-167
lines changed

‎CHANGES.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ The version numbers are in the form `MAJOR.MINOR.PATCH`, where major releases in
2929
* Account for processes without cells in initial condition patch when verifying configuration.
3030
* Avoid divide by zero for `KinSrcRamp` when final slip is zero.
3131
* Remove `solid_bulk_modulus` as a spatial database field for poroelasticity; compute it from the other fields.
32+
* Fixed typos in setting up gravity with poroelasticity.
3233
* Fix pythia import in `pylith_eqinfo`.
3334

3435
## Version 4.2.0 (2025-01-15)

‎configure.ac‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -114,7 +114,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])])
114114

115115
dnl PETSC
116116
AC_LANG(C)
117-
CIT_PATH_PETSC([3.23.4])
117+
CIT_PATH_PETSC([3.23.5])
118118
CIT_HEADER_PETSC
119119
CIT_CHECK_LIB_PETSC
120120

@@ -266,6 +266,7 @@ AC_CONFIG_FILES([Makefile
266266
tests/fullscale/viscoelasticity/nofaults-2d/Makefile
267267
tests/fullscale/viscoelasticity/nofaults-3d/Makefile
268268
tests/fullscale/poroelasticity/Makefile
269+
tests/fullscale/poroelasticity/nofaults-2d/Makefile
269270
tests/fullscale/poroelasticity/faults-2d/Makefile
270271
tests/fullscale/poroelasticity/terzaghi/Makefile
271272
tests/fullscale/poroelasticity/mandel/Makefile

‎libsrc/pylith/fekernels/IsotropicLinearPoroelasticity.hh‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1181,7 +1181,7 @@ public:
11811181
pylith::fekernels::Poroelasticity::Context poroelasticContext;
11821182
pylith::fekernels::Poroelasticity::setContextQuasistatic(
11831183
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
1184-
pylith::fekernels::Poroelasticity::setContextBodyForce(
1184+
pylith::fekernels::Poroelasticity::setContextGravity(
11851185
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
11861186

11871187
// Rheology Context
@@ -1233,7 +1233,7 @@ public:
12331233
pylith::fekernels::Poroelasticity::Context poroelasticContext;
12341234
pylith::fekernels::Poroelasticity::setContextQuasistatic(
12351235
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
1236-
pylith::fekernels::Poroelasticity::setContextBodyForce(
1236+
pylith::fekernels::Poroelasticity::setContextGravity(
12371237
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
12381238

12391239
// Rheology Context
@@ -3287,7 +3287,7 @@ public:
32873287
pylith::fekernels::Poroelasticity::Context poroelasticContext;
32883288
pylith::fekernels::Poroelasticity::setContextQuasistatic(
32893289
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
3290-
pylith::fekernels::Poroelasticity::setContextGravity(
3290+
pylith::fekernels::Poroelasticity::setContextGravityBodyForce(
32913291
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
32923292

32933293
// Rheology Context
@@ -3337,7 +3337,7 @@ public:
33373337
pylith::fekernels::Poroelasticity::Context poroelasticContext;
33383338
pylith::fekernels::Poroelasticity::setContextQuasistatic(
33393339
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
3340-
pylith::fekernels::Poroelasticity::setContextGravity(
3340+
pylith::fekernels::Poroelasticity::setContextGravityBodyForce(
33413341
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
33423342

33433343
// Rheology Context
@@ -3895,7 +3895,7 @@ public:
38953895
pylith::fekernels::Poroelasticity::Context poroelasticContext;
38963896
pylith::fekernels::Poroelasticity::setContextDynamic(
38973897
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
3898-
pylith::fekernels::Poroelasticity::setContextGravitySourceDensity(
3898+
pylith::fekernels::Poroelasticity::setContextGravityBodyForceSourceDensity(
38993899
&poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x);
39003900

39013901
// Rheology Context

‎libsrc/pylith/fekernels/Poroelasticity.hh‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -435,8 +435,8 @@ public:
435435
assert(context);
436436
assert(numS >= 3);
437437

438-
const PylithInt i_gravityField = 4;
439-
const PylithInt i_bodyForce = 5;
438+
const PylithInt i_bodyForce = 4;
439+
const PylithInt i_gravityField = 5;
440440

441441
assert(aOff[i_gravityField] >= 0);
442442
assert(aOff[i_bodyForce] >= 0);
@@ -513,8 +513,8 @@ public:
513513
assert(context);
514514
assert(numS >= 3);
515515

516-
const PylithInt i_gravityField = 4;
517-
const PylithInt i_bodyForce = 5;
516+
const PylithInt i_bodyForce = 4;
517+
const PylithInt i_gravityField = 5;
518518
const PylithInt i_sourceDensity = 6;
519519

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

890-
// Poroelastic Auxiliariies
891890
const PylithScalar bulkDensity = context.bulkDensity;
892891
const PylithScalar* gravityField = context.gravityField;
893892

‎libsrc/pylith/materials/Elasticity.cc‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -249,6 +249,11 @@ pylith::materials::Elasticity::createAuxiliaryField(const pylith::topology::Fiel
249249
assert(auxiliaryFactory);
250250
auxiliaryFactory->setValuesFromDB();
251251

252+
pythia::journal::debug_t debug("elasticity.view_auxiliary_field");
253+
if (debug.state()) {
254+
auxiliaryField->view("Elasticity auxiliary field");
255+
} // if
256+
252257
_Elasticity::Events::logger.eventEnd(_Elasticity::Events::createAuxiliaryField);
253258
PYLITH_METHOD_RETURN(auxiliaryField);
254259
} // createAuxiliaryField

‎libsrc/pylith/materials/IsotropicLinearPoroelasticity.cc‎

Lines changed: 59 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -147,47 +147,39 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelg0p(const spatialdata
147147

148148
switch (bitUse) {
149149
case 0x0:
150-
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
151-
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
152-
NULL;
153-
break;
154150
case 0x1:
155-
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
156-
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
157-
NULL;
158-
break;
159151
case 0x2:
152+
case 0x3:
160153
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
161154
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
162155
NULL;
163156
break;
164157
case 0x4:
158+
// aOff for sourceDensity is 3
165159
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source :
166160
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source :
167-
NULL; // aOff for sourceDensity is 3
168-
break;
169-
case 0x3:
170-
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p :
171-
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p :
172161
NULL;
173162
break;
174163
case 0x5:
164+
// aOff for sourceDensity is 4
175165
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source_body :
176166
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source_body :
177-
NULL; // aOff for sourceDensity is 4
167+
NULL;
178168
break;
179169
case 0x6:
170+
// aOff for sourceDensity is 4
180171
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source_grav :
181172
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source_grav :
182-
NULL; // aOff for sourceDensity is 4
173+
NULL;
183174
break;
184175
case 0x7:
176+
// aOff for sourceDensity is 5
185177
g0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::g0p_source_grav_body :
186178
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::g0p_source_grav_body :
187-
NULL; // aOff for sourceDensity is 5
179+
NULL;
188180
break;
189181
default:
190-
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ") for Poroelasticity RHS residual kernels.");
182+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
191183
} // switch
192184

193185
PYLITH_METHOD_RETURN(g0p);
@@ -231,8 +223,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelg1p_explicit(const sp
231223
NULL;
232224
break;
233225
default:
234-
PYLITH_COMPONENT_ERROR("Unknown combination of flags for _useTensorPermeability="<<_useTensorPermeability<<", _gravityField="<<_gravityField<<").");
235-
throw std::logic_error("Unknown combination of flags.");
226+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
236227
} // switch
237228

238229
PYLITH_METHOD_RETURN(g1p);
@@ -264,8 +255,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelg1v_explicit(const sp
264255
NULL;
265256
break;
266257
default:
267-
PYLITH_COMPONENT_ERROR("Unknown combination of flags for getKernelg1v_explicit.");
268-
throw std::logic_error("Unknown combination of flags.");
258+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
269259
} // switch
270260

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

312302
switch (bitUse) {
313303
case 0x0:
314-
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
315-
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
316-
NULL;
317-
break;
318304
case 0x1:
319-
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
320-
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
321-
NULL;
322-
break;
323305
case 0x2:
306+
case 0x3:
324307
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
325308
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
326309
NULL;
327310
break;
328311
case 0x4:
312+
// aOff for sourceDensity is 3
329313
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source :
330314
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source :
331-
NULL; // aOff for sourceDensity is 3
332-
break;
333-
case 0x3:
334-
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit :
335-
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit :
336315
NULL;
337316
break;
338317
case 0x5:
318+
// aOff for sourceDensity is 4
339319
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source_body :
340320
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source_body :
341-
NULL; // aOff for sourceDensity is 4
321+
NULL;
342322
break;
343323
case 0x6:
324+
// aOff for sourceDensity is 4
344325
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source_grav :
345326
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source_grav :
346-
NULL; // aOff for sourceDensity is 4
327+
NULL;
347328
break;
348329
case 0x7:
330+
// aOff for sourceDensity is 5
349331
f0p = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::f0p_implicit_source_grav_body :
350332
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::f0p_implicit_source_grav_body :
351-
NULL; // aOff for sourceDensity is 5
333+
NULL;
352334
break;
353335
default:
354-
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ") for Poroelasticity LHS residual kernels.");
336+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
355337
} // switch
356338

357339
PYLITH_METHOD_RETURN(f0p);
@@ -383,8 +365,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelf1u_implicit(const sp
383365
NULL;
384366
break;
385367
default:
386-
PYLITH_COMPONENT_ERROR("Unknown combination of flags for getKernelf1u_implicit.");
387-
throw std::logic_error("Unknown combination of flags.");
368+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
388369
} // switch
389370

390371
PYLITH_METHOD_RETURN(f1u);
@@ -451,8 +432,7 @@ pylith::materials::IsotropicLinearPoroelasticity::getKernelf1p_implicit(const sp
451432
break;
452433

453434
default:
454-
PYLITH_COMPONENT_ERROR("Unknown combination of flags for _useTensorPermeability="<<_useTensorPermeability<<", _useBodyForce="<<_useBodyForce<<", _gravityField="<<_gravityField<<").");
455-
throw std::logic_error("Unknown combination of flags.");
435+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
456436
} // switch
457437

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

537517
const int spaceDim = coordsys->getSpaceDim();
538-
PetscPointJacFn* Jf3pp =
539-
(!_useTensorPermeability && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp :
540-
(!_useTensorPermeability && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp :
541-
(_useTensorPermeability && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp_tensor_permeability :
542-
(_useTensorPermeability && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp_tensor_permeability :
543-
NULL;
518+
const int bitTensorPermeability = _useTensorPermeability ? 0x1 : 0x0;
519+
const int bitUse = bitTensorPermeability;
520+
521+
PetscPointJacFn* Jf3pp = NULL;
522+
switch (bitUse) {
523+
case 0x0:
524+
Jf3pp = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp :
525+
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp :
526+
NULL;
527+
break;
528+
case 0x1:
529+
Jf3pp = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::Jf3pp_tensor_permeability :
530+
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::Jf3pp_tensor_permeability :
531+
NULL;
532+
break;
533+
default:
534+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
535+
} // switch
544536

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

575567
const int spaceDim = coordsys->getSpaceDim();
576-
PetscPointFn* kernel =
577-
(!_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_asVector :
578-
(!_useReferenceState && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_asVector :
579-
(_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_refState_asVector :
580-
(_useReferenceState && 2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_refState_asVector :
581-
NULL;
568+
const int bitReferenceState = _useReferenceState ? 0x1 : 0x0;
569+
const int bitUse = bitReferenceState;
570+
571+
PetscPointFn* kernel = NULL;
572+
573+
switch (bitUse) {
574+
case 0x0:
575+
kernel = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_asVector :
576+
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_asVector :
577+
NULL;
578+
break;
579+
case 0x1:
580+
kernel = (3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::cauchyStress_infinitesimalStrain_refState_asVector :
581+
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::cauchyStress_infinitesimalStrain_refState_asVector :
582+
NULL;
583+
break;
584+
default:
585+
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ").");
586+
} // switch
582587

583588
PYLITH_METHOD_RETURN(kernel);
584589
} // getKernelCauchyStressVector

‎libsrc/pylith/materials/Poroelasticity.cc‎

Lines changed: 6 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -254,6 +254,11 @@ pylith::materials::Poroelasticity::createAuxiliaryField(const pylith::topology::
254254
assert(auxiliaryFactory);
255255
auxiliaryFactory->setValuesFromDB();
256256

257+
pythia::journal::debug_t debug("poroelasticity.view_auxiliary_field");
258+
if (debug.state()) {
259+
auxiliaryField->view("Poroelasticity auxiliary field");
260+
} // if
261+
257262
PYLITH_METHOD_RETURN(auxiliaryField);
258263
} // createAuxiliaryField
259264

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

396401
const int bitBodyForce = _useBodyForce ? 0x1 : 0x0;
397402
const int bitGravity = _gravityField ? 0x2 : 0x0;
398-
const int bitSourceDensity = _useSourceDensity ? 0x4 : 0x0;
399-
const int bitUse = bitBodyForce | bitGravity | bitSourceDensity;
403+
const int bitUse = bitBodyForce | bitGravity;
400404

401405
PetscPointFn* r0 = NULL;
402406
switch (bitUse) {
@@ -408,20 +412,9 @@ pylith::materials::Poroelasticity::_setKernelsResidual(pylith::feassemble::Integ
408412
case 0x2:
409413
r0 = pylith::fekernels::Poroelasticity::g0v_grav;
410414
break;
411-
case 0x4:
412-
break;
413415
case 0x3:
414416
r0 = pylith::fekernels::Poroelasticity::g0v_grav_bodyforce;
415417
break;
416-
case 0x5:
417-
r0 = pylith::fekernels::Poroelasticity::g0v_bodyforce;
418-
break;
419-
case 0x6:
420-
r0 = pylith::fekernels::Poroelasticity::g0v_grav;
421-
break;
422-
case 0x7:
423-
r0 = pylith::fekernels::Poroelasticity::g0v_grav_bodyforce;
424-
break;
425418
default:
426419
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ") for residual kernels.");
427420
} // switch

0 commit comments

Comments
 (0)