Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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 .vscode/settings.json
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
"dm",
"drucker",
"elastoplastic",
"Gmsh",
"Knepley",
"neumann",
"nodeset",
Expand Down
15 changes: 9 additions & 6 deletions CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,16 +13,19 @@ The version numbers are in the form `MAJOR.MINOR.PATCH`, where major releases in
* Add displacement scale and improve the names of the other scales.
* The time scale for poroelasticity is derived from the other scales and nominal material properties.
* Add discussion of nondimensionalization for each of the governing equations.
* Update default solver tolerances consistent with new nondimensionalization scales.
* Update default solver tolerances consistent with the new nondimensionalization scales.
* Add option to use adaptive time stepping and update examples `reverse-2d` Steps 7 and 8 and `magma-2d` Step 1 to demonstrate its use.
* Update solver settings for elasticity for better scalability with problem size.
* ASCII mesh format
* Changed `group` to `vertex-group`; remove `group` `type`.
* Renamed `VertexGroup` to `BoundaryGroup` in `meshio.gmsh_utils` and changed default behavior to not be recursive (generate "face" groups, not "vertex" groups).
* Update multi-grid preconditioner solver settings for elasticity for better scalability with problem size.
* Mesh formats
* Add support for use of Cubit side sets and Gmsh physical groups with edges (2D) and surfaces (3D). Use of Cubit vertex sets and Gmsh physical groups with vertices are deprecated and support for them will be removed in version 6.0.0.
* Use of Cubit side sets and Gmsh physical groups with curves (2D) and surfaces (3D) are using
* ASCII mesh format
* Changed `group` to `vertex-group`; remove `group` `type`.
* Renamed `VertexGroup` to `BoundaryGroup` in `meshio.gmsh_utils` and changed default behavior to not be recursive (generate "face" groups, not "vertex" groups).
* Use the VTU (XML) format for VTK files instead of the legacy ASCII format.
* Material `description` property is no longer used; a deprecation warning is printed to stdout if it is specified. This feature will be removed in v6.0.
* Internal labels for meshes and fields are now derived from the Pyre component name and identifier.
* Update PETSc to 3.23.1
* Update PETSc to 3.24.0
* **Added**
* Improved documentation for `pylith_eqinfo` and illustrate use in `examples/strikdslip-2d`, `examples/crustal-strikeslip-2d`, and `examples/crustal-strikeslip-3d`.
* Added `pylith_convertmesh` for converting Cubit and Gmsh files to the PETSc HDF5 mesh format.
Expand Down
2 changes: 1 addition & 1 deletion libsrc/pylith/fekernels/IsotropicLinearGenMaxwell.cc
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,6 @@

#include "pylith/fekernels/IsotropicLinearGenMaxwell.hh" // Implementation of object methods.

const PylithInt pylith::fekernels::IsotropicLinearGenMaxwell::numParallel = 3;
// const PylithInt pylith::fekernels::IsotropicLinearGenMaxwell::numParallel = 3;

// End of file
2 changes: 1 addition & 1 deletion libsrc/pylith/fekernels/IsotropicLinearGenMaxwell.hh
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ class pylith::fekernels::IsotropicLinearGenMaxwell {
// PUBLIC MEMBERS /////////////////////////////////////////////////////////////////////////////
public:

static const PylithInt numParallel;
static const PylithInt numParallel = 3;

struct Context {
PylithReal shearModulus;
Expand Down
40 changes: 20 additions & 20 deletions libsrc/pylith/materials/Elasticity.cc
Original file line number Diff line number Diff line change
Expand Up @@ -472,30 +472,30 @@ pylith::materials::Elasticity::_setKernelsResidual(pylith::feassemble::Integrato
default:
PYLITH_COMPONENT_LOGICERROR("Unknown case (bitUse=" << bitUse << ") for residual kernels.");
} // switch
const PetscPointFn* r1 = _rheology->getKernelf1v(coordsys);
PetscPointFn* r1 = _rheology->getKernelf1v(coordsys);

std::vector<ResidualKernels> kernels;
switch (_formulation) {
case QUASISTATIC: {
const PetscPointFn* f0u = r0;
const PetscPointFn* f1u = r1;
PetscPointFn* f0u = r0;
PetscPointFn* f1u = r1;

kernels.resize(1);
kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);
break;
} // QUASISTATIC
case DYNAMIC: {
// Displacement
const PetscPointFn* f0u = pylith::fekernels::DispVel::f0u;
const PetscPointFn* f1u = NULL;
const PetscPointFn* g0u = pylith::fekernels::DispVel::g0u;
const PetscPointFn* g1u = NULL;
PetscPointFn* f0u = pylith::fekernels::DispVel::f0u;
PetscPointFn* f1u = NULL;
PetscPointFn* g0u = pylith::fekernels::DispVel::g0u;
PetscPointFn* g1u = NULL;

// Velocity
const PetscPointFn* f0v = pylith::fekernels::Elasticity::f0v;
const PetscPointFn* f1v = NULL;
const PetscPointFn* g0v = r0;
const PetscPointFn* g1v = r1;
PetscPointFn* f0v = pylith::fekernels::Elasticity::f0v;
PetscPointFn* f1v = NULL;
PetscPointFn* g0v = r0;
PetscPointFn* g1v = r1;

kernels.resize(4);
kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);
Expand All @@ -506,16 +506,16 @@ pylith::materials::Elasticity::_setKernelsResidual(pylith::feassemble::Integrato
} // DYNAMIC
case DYNAMIC_IMEX: {
// Displacement
const PetscPointFn* f0u = pylith::fekernels::DispVel::f0u;
const PetscPointFn* f1u = NULL;
const PetscPointFn* g0u = pylith::fekernels::DispVel::g0u;
const PetscPointFn* g1u = NULL;
PetscPointFn* f0u = pylith::fekernels::DispVel::f0u;
PetscPointFn* f1u = NULL;
PetscPointFn* g0u = pylith::fekernels::DispVel::g0u;
PetscPointFn* g1u = NULL;

// Velocity
const PetscPointFn* f0v = pylith::fekernels::DispVel::f0v;
const PetscPointFn* f1v = NULL;
const PetscPointFn* g0v = r0;
const PetscPointFn* g1v = r1;
PetscPointFn* f0v = pylith::fekernels::DispVel::f0v;
PetscPointFn* f1v = NULL;
PetscPointFn* g0v = r0;
PetscPointFn* g1v = r1;

kernels.resize(4);
kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);
Expand Down Expand Up @@ -638,7 +638,7 @@ pylith::materials::Elasticity::_setKernelsDerivedField(pylith::feassemble::Integ
kernels[0] = ProjectKernels("cauchy_stress", _rheology->getKernelCauchyStressVector(coordsys));

const int spaceDim = coordsys->getSpaceDim();
const PetscPointFn* strainKernel =
PetscPointFn* strainKernel =
(3 == spaceDim) ? pylith::fekernels::Elasticity3D::infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::ElasticityPlaneStrain::infinitesimalStrain_asVector :
NULL;
Expand Down
10 changes: 5 additions & 5 deletions libsrc/pylith/materials/IncompressibleElasticity.cc
Original file line number Diff line number Diff line change
Expand Up @@ -328,12 +328,12 @@ pylith::materials::IncompressibleElasticity::_setKernelsResidual(pylith::feassem
} // switch

// Displacement
const PetscPointFn* f0u = r0;
const PetscPointFn* f1u = _rheology->getKernelf1u(coordsys);
PetscPointFn* f0u = r0;
PetscPointFn* f1u = _rheology->getKernelf1u(coordsys);

// Pressure
const PetscPointFn* f0p = _rheology->getKernelf0p(coordsys);
const PetscPointFn* f1p = NULL;
PetscPointFn* f0p = _rheology->getKernelf0p(coordsys);
PetscPointFn* f1p = NULL;

std::vector<ResidualKernels> kernels(2);
kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);
Expand Down Expand Up @@ -429,7 +429,7 @@ pylith::materials::IncompressibleElasticity::_setKernelsDerivedField(pylith::fea
kernels[0] = ProjectKernels("cauchy_stress", _rheology->getKernelCauchyStressVector(coordsys));

const int spaceDim = coordsys->getSpaceDim();
const PetscPointFn* strainKernel =
PetscPointFn* strainKernel =
(3 == spaceDim) ? pylith::fekernels::Elasticity3D::infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::ElasticityPlaneStrain::infinitesimalStrain_asVector :
NULL;
Expand Down
4 changes: 2 additions & 2 deletions libsrc/pylith/materials/IsotropicLinearGenMaxwell.cc
Original file line number Diff line number Diff line change
Expand Up @@ -224,11 +224,11 @@ pylith::materials::IsotropicLinearGenMaxwell::addKernelsUpdateStateVars(std::vec
PYLITH_COMPONENT_DEBUG("addKernelsUpdateStateVars(kernels="<<kernels<<", coordsys="<<coordsys<<")");

const int spaceDim = coordsys->getSpaceDim();
const PetscPointFn* funcViscousStrain =
PetscPointFn* funcViscousStrain =
(3 == spaceDim) ? pylith::fekernels::IsotropicLinearGenMaxwell3D::viscousStrain_infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearGenMaxwellPlaneStrain::viscousStrain_infinitesimalStrain_asVector :
NULL;
const PetscPointFn* funcTotalStrain =
PetscPointFn* funcTotalStrain =
(3 == spaceDim) ? pylith::fekernels::Elasticity3D::infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::ElasticityPlaneStrain::infinitesimalStrain_asVector :
NULL;
Expand Down
4 changes: 2 additions & 2 deletions libsrc/pylith/materials/IsotropicLinearMaxwell.cc
Original file line number Diff line number Diff line change
Expand Up @@ -223,11 +223,11 @@ pylith::materials::IsotropicLinearMaxwell::addKernelsUpdateStateVars(std::vector
PYLITH_COMPONENT_DEBUG("addKernelsUpdateStateVars(kernels="<<kernels<<", coordsys="<<coordsys<<")");

const int spaceDim = coordsys->getSpaceDim();
const PetscPointFn* funcViscousStrain =
PetscPointFn* funcViscousStrain =
(3 == spaceDim) ? pylith::fekernels::IsotropicLinearMaxwell3D::viscousStrain_infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearMaxwellPlaneStrain::viscousStrain_infinitesimalStrain_asVector :
NULL;
const PetscPointFn* funcTotalStrain =
PetscPointFn* funcTotalStrain =
(3 == spaceDim) ? pylith::fekernels::Elasticity3D::infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::ElasticityPlaneStrain::infinitesimalStrain_asVector :
NULL;
Expand Down
4 changes: 2 additions & 2 deletions libsrc/pylith/materials/IsotropicLinearPoroelasticity.cc
Original file line number Diff line number Diff line change
Expand Up @@ -634,7 +634,7 @@ pylith::materials::IsotropicLinearPoroelasticity::addKernelsUpdateStateVarsImpli
if (_useStateVars) {
const int spaceDim = coordsys->getSpaceDim();

const PetscPointFn* funcPorosity =
PetscPointFn* funcPorosity =
(3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::updatePorosityImplicit :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::updatePorosityImplicit :
NULL;
Expand All @@ -660,7 +660,7 @@ pylith::materials::IsotropicLinearPoroelasticity::addKernelsUpdateStateVarsExpli
if (_useStateVars) {
const int spaceDim = coordsys->getSpaceDim();

const PetscPointFn* funcPorosity =
PetscPointFn* funcPorosity =
(3 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticity3D::updatePorosityExplicit :
(2 == spaceDim) ? pylith::fekernels::IsotropicLinearPoroelasticityPlaneStrain::updatePorosityExplicit :
NULL;
Expand Down
4 changes: 2 additions & 2 deletions libsrc/pylith/materials/IsotropicPowerLaw.cc
Original file line number Diff line number Diff line change
Expand Up @@ -227,13 +227,13 @@ pylith::materials::IsotropicPowerLaw::addKernelsUpdateStateVars(std::vector<Proj
PYLITH_COMPONENT_DEBUG("addKernelsUpdateStateVars(kernels="<<kernels<<", coordsys="<<coordsys<<")");

const int spaceDim = coordsys->getSpaceDim();
const PetscPointFn* funcViscousStrain =
PetscPointFn* funcViscousStrain =
(!_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicPowerLaw3D::viscousStrain_infinitesimalStrain_asVector :
(!_useReferenceState && 2 == spaceDim) ? pylith::fekernels::IsotropicPowerLawPlaneStrain::viscousStrain_infinitesimalStrain_asVector :
(_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicPowerLaw3D::viscousStrain_infinitesimalStrain_refState_asVector :
(_useReferenceState && 2 == spaceDim) ? pylith::fekernels::IsotropicPowerLawPlaneStrain::viscousStrain_infinitesimalStrain_refState_asVector :
NULL;
const PetscPointFn* funcDeviatoricStress =
PetscPointFn* funcDeviatoricStress =
(!_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicPowerLaw3D::deviatoricStress_infinitesimalStrain_asVector :
(!_useReferenceState && 2 == spaceDim) ? pylith::fekernels::IsotropicPowerLawPlaneStrain::deviatoricStress_infinitesimalStrain_asVector :
(_useReferenceState && 3 == spaceDim) ? pylith::fekernels::IsotropicPowerLaw3D::deviatoricStress_infinitesimalStrain_refState_asVector :
Expand Down
54 changes: 27 additions & 27 deletions libsrc/pylith/materials/Poroelasticity.cc
Original file line number Diff line number Diff line change
Expand Up @@ -421,17 +421,17 @@ pylith::materials::Poroelasticity::_setKernelsResidual(pylith::feassemble::Integ
case QUASISTATIC: {
if (!_useStateVars) {
// Displacement
const PetscPointFn* f0u = r0;
const PetscPointFn* f1u = _rheology->getKernelf1u_implicit(coordsys);
PetscPointFn* f0u = r0;
PetscPointFn* f1u = _rheology->getKernelf1u_implicit(coordsys);

// Pressure
PetscPointFn* f0p = _rheology->getKernelf0p_implicit(coordsys, _useBodyForce, _gravityField, _useSourceDensity);
PetscPointFn* f1p = _rheology->getKernelf1p_implicit(coordsys, _useBodyForce, _gravityField); // darcy
// velocity

// Volumetric Strain
const PetscPointFn* f0e = pylith::fekernels::Poroelasticity::f0e;
const PetscPointFn* f1e = NULL;
PetscPointFn* f0e = pylith::fekernels::Poroelasticity::f0e;
PetscPointFn* f1e = NULL;

kernels.resize(3);
kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);
Expand All @@ -440,27 +440,27 @@ pylith::materials::Poroelasticity::_setKernelsResidual(pylith::feassemble::Integ
} else {
// Displacement
PetscPointFn* f0u = r0;
const PetscPointFn* f1u = _rheology->getKernelf1u_implicit(coordsys);
PetscPointFn* f1u = _rheology->getKernelf1u_implicit(coordsys);

// Pressure
PetscPointFn* f0p = _rheology->getKernelf0p_implicit(coordsys, _useBodyForce, _gravityField, _useSourceDensity);
PetscPointFn* f1p = _rheology->getKernelf1p_implicit(coordsys, _useBodyForce, _gravityField);

// Volumetric Strain
const PetscPointFn* f0e = pylith::fekernels::Poroelasticity::f0e;
const PetscPointFn* f1e = NULL;
PetscPointFn* f0e = pylith::fekernels::Poroelasticity::f0e;
PetscPointFn* f1e = NULL;

// Velocity
const PetscPointFn* f0v = pylith::fekernels::Poroelasticity::f0v_implicit;
const PetscPointFn* f1v = NULL;
PetscPointFn* f0v = pylith::fekernels::Poroelasticity::f0v_implicit;
PetscPointFn* f1v = NULL;

// Time derivative of pressure
const PetscPointFn* f0pdot = pylith::fekernels::Poroelasticity::f0pdot;
const PetscPointFn* f1pdot = NULL;
PetscPointFn* f0pdot = pylith::fekernels::Poroelasticity::f0pdot;
PetscPointFn* f1pdot = NULL;

// Time derivative of Volumetric Strain
const PetscPointFn* f0edot = pylith::fekernels::Poroelasticity::f0edot;
const PetscPointFn* f1edot = NULL;
PetscPointFn* f0edot = pylith::fekernels::Poroelasticity::f0edot;
PetscPointFn* f1edot = NULL;

kernels.resize(6);
kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);
Expand All @@ -475,22 +475,22 @@ pylith::materials::Poroelasticity::_setKernelsResidual(pylith::feassemble::Integ
case DYNAMIC_IMEX:
case DYNAMIC: {
// Displacement
const PetscPointFn* f0u = pylith::fekernels::DispVel::f0u;
const PetscPointFn* f1u = NULL;
const PetscPointFn* g0u = pylith::fekernels::Poroelasticity::g0u;
const PetscPointFn* g1u = NULL;
PetscPointFn* f0u = pylith::fekernels::DispVel::f0u;
PetscPointFn* f1u = NULL;
PetscPointFn* g0u = pylith::fekernels::Poroelasticity::g0u;
PetscPointFn* g1u = NULL;

// Pressure
const PetscPointFn* f0p = _rheology->getKernelf0p_explicit(coordsys);
const PetscPointFn* f1p = NULL;
const PetscPointFn* g0p = _rheology->getKernelg0p(coordsys, _useBodyForce, _gravityField, _useSourceDensity);
const PetscPointFn* g1p = _rheology->getKernelg1p_explicit(coordsys, _gravityField); // Darcy velocity
PetscPointFn* f0p = _rheology->getKernelf0p_explicit(coordsys);
PetscPointFn* f1p = NULL;
PetscPointFn* g0p = _rheology->getKernelg0p(coordsys, _useBodyForce, _gravityField, _useSourceDensity);
PetscPointFn* g1p = _rheology->getKernelg1p_explicit(coordsys, _gravityField); // Darcy velocity

// Velocity
const PetscPointFn* f0v = pylith::fekernels::Poroelasticity::f0v_explicit;
const PetscPointFn* f1v = NULL;
const PetscPointFn* g0v = r0;
const PetscPointFn* g1v = _rheology->getKernelg1v_explicit(coordsys);
PetscPointFn* f0v = pylith::fekernels::Poroelasticity::f0v_explicit;
PetscPointFn* f1v = NULL;
PetscPointFn* g0v = r0;
PetscPointFn* g1v = _rheology->getKernelg1v_explicit(coordsys);

kernels.resize(6);
kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);
Expand Down Expand Up @@ -775,13 +775,13 @@ pylith::materials::Poroelasticity::_setKernelsDerivedField(pylith::feassemble::I
kernels[0] = ProjectKernels("cauchy_stress", _rheology->getKernelCauchyStressVector(coordsys));

const int spaceDim = coordsys->getSpaceDim();
const PetscPointFn* strainKernel =
PetscPointFn* strainKernel =
(3 == spaceDim) ? pylith::fekernels::Elasticity3D::infinitesimalStrain_asVector :
(2 == spaceDim) ? pylith::fekernels::ElasticityPlaneStrain::infinitesimalStrain_asVector :
NULL;
kernels[1] = ProjectKernels("cauchy_strain", strainKernel);

const PetscPointFn* bulkDensity = pylith::fekernels::Poroelasticity::bulkDensity_asScalar;
PetscPointFn* bulkDensity = pylith::fekernels::Poroelasticity::bulkDensity_asScalar;

kernels[2] = ProjectKernels("bulk_density", bulkDensity);

Expand Down