21 #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_
22 #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSSOLVER_HPP_
28 #include "constitutive/solid/CoupledSolidBase.hpp"
29 #include "constitutive/contact/HydraulicApertureBase.hpp"
32 #include "codingUtilities/Utilities.hpp"
37 namespace stabilization
39 enum class StabilizationType :
integer
53 template<
typename FLOW_SOLVER,
typename MECHANICS_SOLVER = Sol
idMechanicsLagrangianFEM >
81 :
Base( name, parent ),
86 setApplyDefaultValue( 0 ).
88 setDescription(
"Flag indicating whether the problem is thermal or not. Set isThermal=\"1\" to enable the thermal coupling" );
92 setDescription(
"StabilizationType. Options are:\n" +
93 toString( stabilization::StabilizationType::None ) +
"- Add no stabilization to mass equation \n" +
94 toString( stabilization::StabilizationType::Global ) +
"- Add jump stabilization to all faces \n" +
95 toString( stabilization::StabilizationType::Local ) +
"- Add jump stabilization on interior of macro elements" );
98 setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ).
100 setDescription(
"Regions where stabilization is applied." );
103 setApplyDefaultValue( 1.0 ).
105 setDescription(
"Constant multiplier of stabilization strength" );
113 GEOS_FMT(
"{} {}: The attribute `{}` of the flow solver must be thermal since the poromechanics solver is thermal",
118 GEOS_FMT(
"{} {}: The attribute `{}` of solid mechanics solver `{}` must be `{}`",
120 SolidMechanicsLagrangianFEM::viewKeyStruct::timeIntegrationOptionString(),
130 this->
template setConstitutiveName< constitutive::CoupledSolidBase >( subRegion,
134 this->
template setConstitutiveName< constitutive::PorosityBase >( subRegion,
135 constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString(),
"porosity" );
139 this->
template setConstitutiveName< constitutive::HydraulicApertureBase >( subRegion,
150 ": Local stabilization has been temporarily disabled",
153 DomainPartition & domain = this->
template getGroupByPath< DomainPartition >(
"/Problem/domain" );
155 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&] (
string const &,
159 ElementRegionManager & elementRegionManager = mesh.getElemManager();
160 elementRegionManager.forElementSubRegions< CellElementSubRegion >( regionNames,
161 [&]( localIndex const,
162 ElementSubRegionBase & subRegion )
164 if( subRegion.hasField< fields::poromechanics::bulkDensity >() )
167 constitutive::CoupledSolidBase const & solid =
168 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
169 subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
170 subRegion.getField< fields::poromechanics::bulkDensity >().resizeDimension< 1 >( solid.getDensity().size( 1 ) );
178 Base::registerDataOnMesh( meshBodies );
183 solidMechanicsSolver()->enableFixedStressPoromechanicsUpdate();
185 flowSolver()->enableFixedStressPoromechanicsUpdate();
188 if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
190 flowSolver()->enableJumpStabilization();
193 this->forDiscretizationOnMeshTargets( meshBodies, [&] (
string const &,
203 subRegion.registerWrapper<
string >( viewKeyStruct::porousMaterialNamesString() ).
206 setSizedFromParent( 0 );
209 subRegion.registerWrapper<
string >( constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString() ).
212 setSizedFromParent( 0 );
216 subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() );
218 if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
220 subRegion.registerField< fields::flow::macroElementIndex >( this->getName());
221 subRegion.registerField< fields::flow::elementStabConstant >( this->getName());
231 flowSolver()->setKeepVariablesConstantDuringInitStep( m_performStressInitialization );
233 if( this->m_stabilizationType == stabilization::StabilizationType::Global || this->m_stabilizationType == stabilization::StabilizationType::Local )
235 this->updateStabilizationParameters( domain );
238 Base::implicitStepSetup( time_n, dt, domain );
247 solidMechanicsSolver()->setupDofs( domain, dofManager );
248 flowSolver()->setupDofs( domain, dofManager );
249 this->setupCoupling( domain, dofManager );
252 virtual bool checkSequentialConvergence(
int const & iter,
258 auto & subcycling = this->getNonlinearSolverParameters().m_subcyclingOption;
259 auto const subcycling_orig = subcycling;
260 if( m_performStressInitialization )
263 bool isConverged = Base::checkSequentialConvergence( iter, time_n, dt, domain );
266 subcycling = subcycling_orig;
277 return std::get< toUnderlying( SolverType::SolidMechanics ) >( m_solvers );
286 return std::get< toUnderlying( SolverType::Flow ) >( m_solvers );
293 void setStressInitialization(
bool const performStressInitialization )
295 m_performStressInitialization = performStressInitialization;
296 solidMechanicsSolver()->setStressInitialization( performStressInitialization );
325 for(
string const & regionName : m_stabilizationRegionNames )
327 regionFilter.insert( regionName );
331 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&] (
string const &,
336 string_array filteredTargetRegionNames;
337 filteredTargetRegionNames.reserve( targetRegionNames.size() );
339 for( string const & targetRegionName : targetRegionNames )
341 if( regionFilter.count( targetRegionName ) )
343 filteredTargetRegionNames.emplace_back( targetRegionName );
348 mesh.getElemManager().forElementSubRegions( filteredTargetRegionNames, [&](
localIndex const,
349 ElementSubRegionBase & subRegion )
351 arrayView1d< integer > const macroElementIndex = subRegion.getField< fields::flow::macroElementIndex >();
352 arrayView1d< real64 > const elementStabConstant = subRegion.getField< fields::flow::elementStabConstant >();
354 constitutive::CoupledSolidBase const & porousSolid =
355 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
356 subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
358 arrayView1d< real64 const > const bulkModulus = porousSolid.getBulkModulus();
359 arrayView1d< real64 const > const shearModulus = porousSolid.getShearModulus();
360 arrayView1d< real64 const > const biotCoefficient = porousSolid.getBiotCoefficient();
362 real64 const stabilizationMultiplier = m_stabilizationMultiplier;
364 forAll< parallelDevicePolicy<> >( subRegion.size(), [bulkModulus,
367 stabilizationMultiplier,
369 elementStabConstant] GEOS_HOST_DEVICE ( localIndex const ei )
372 real64 const bM = bulkModulus[ei];
373 real64 const sM = shearModulus[ei];
374 real64 const bC = biotCoefficient[ei];
376 macroElementIndex[ei] = 1;
377 elementStabConstant[ei] = stabilizationMultiplier * 9.0 * (bC * bC) / (32.0 * (10.0 * sM / 3.0 + bM));
384 void updateBulkDensity( DomainPartition & domain )
386 this->
template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&](
string const &,
390 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
396 updateBulkDensity( subRegion );
403 template<
typename TYPE_LIST,
404 typename KERNEL_WRAPPER,
405 typename ... PARAMS >
406 real64 assemblyLaunch( MeshLevel & mesh,
407 DofManager
const & dofManager,
409 string const & materialNamesString,
410 CRSMatrixView< real64, globalIndex const >
const & localMatrix,
411 arrayView1d< real64 >
const & localRhs,
413 PARAMS && ... params )
417 NodeManager
const & nodeManager = mesh.getNodeManager();
419 string const dofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
420 arrayView1d< globalIndex const >
const & dofNumber = nodeManager.getReference<
globalIndex_array >( dofKey );
422 real64 const gravityVectorData[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( this->gravityVector() );
424 KERNEL_WRAPPER kernelWrapper( dofNumber,
425 dofManager.rankOffset(),
430 std::forward< PARAMS >( params )... );
432 return finiteElement::
433 regionBasedKernelApplication< parallelDevicePolicy< >,
436 this->solidMechanicsSolver()->getDiscretizationName(),
443 void recordAverageMeanTotalStressIncrement( DomainPartition & domain,
444 array1d< real64 > & averageMeanTotalStressIncrement )
446 averageMeanTotalStressIncrement.resize( 0 );
447 this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&](
string const &,
451 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
455 string const & solidName =
456 subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
457 constitutive::CoupledSolidBase & solid =
458 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
460 arrayView1d< const real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
461 for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
463 averageMeanTotalStressIncrement.emplace_back( averageMeanTotalStressIncrement_k[k] );
469 void applyAcceleratedAverageMeanTotalStressIncrement( DomainPartition & domain,
470 array1d< real64 > & averageMeanTotalStressIncrement )
473 this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&](
string const &,
477 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
481 string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
482 constitutive::CoupledSolidBase & solid =
483 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
484 auto & porosityModel = dynamic_cast< constitutive::BiotPorosity const & >( solid.getBasePorosityModel());
485 arrayView1d< real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
486 for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
488 porosityModel.updateAverageMeanTotalStressIncrement( k, averageMeanTotalStressIncrement[i] );
495 real64 computeAitkenRelaxationFactor( array1d< real64 >
const & s0,
496 array1d< real64 >
const & s1,
497 array1d< real64 >
const & s1_tilde,
498 array1d< real64 >
const & s2_tilde,
501 array1d< real64 > r1 = axpy( s1_tilde, s0, -1.0 );
502 array1d< real64 > r2 = axpy( s2_tilde, s1, -1.0 );
505 array1d< real64 > diff = axpy( r2, r1, -1.0 );
507 real64 const denom = dot( diff, diff );
508 real64 const numer = dot( r1, diff );
511 if( !isZero( denom ))
513 omega1 = -1.0 * omega0 * numer / denom;
518 array1d< real64 > computeUpdate( array1d< real64 >
const & s1,
519 array1d< real64 >
const & s2_tilde,
522 return axpy( scale( s1, 1.0 - omega1 ),
523 scale( s2_tilde, omega1 ),
527 void startSequentialIteration(
integer const & iter,
528 DomainPartition & domain )
override
530 if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
534 recordAverageMeanTotalStressIncrement( domain, m_s1 );
540 m_s1_tilde = m_s2_tilde;
546 void finishSequentialIteration(
integer const & iter,
547 DomainPartition & domain )
override
549 if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
558 m_omega1 = computeAitkenRelaxationFactor( m_s0, m_s1, m_s1_tilde, m_s2_tilde, m_omega0 );
559 m_s2 = computeUpdate( m_s1, m_s2_tilde, m_omega1 );
560 applyAcceleratedAverageMeanTotalStressIncrement( domain, m_s2 );
570 if( solverType ==
static_cast< integer >( SolverType::Flow ) )
572 updateBulkDensity( domain );
576 if( solverType ==
static_cast< integer >( SolverType::SolidMechanics )
577 && !m_performStressInitialization )
580 averageMeanTotalStressIncrement( domain );
582 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&](
string const &,
587 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
592 flowSolver()->updatePorosityAndPermeability( subRegion );
594 updateBulkDensity( subRegion );
600 if( solverType ==
static_cast< integer >( SolverType::SolidMechanics ) &&
601 this->getNonlinearSolverParameters().m_nonlinearAccelerationType== NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
603 recordAverageMeanTotalStressIncrement( domain, m_s2_tilde );
613 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&](
string const &,
617 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
621 string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
622 constitutive::CoupledSolidBase & solid =
623 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
625 arrayView2d< real64 const > const meanTotalStressIncrement_k = solid.getMeanTotalStressIncrement_k();
626 arrayView1d< real64 > const averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
628 finiteElement::FiniteElementBase & subRegionFE =
629 subRegion.template getReference< finiteElement::FiniteElementBase >( solidMechanicsSolver()->getDiscretizationName() );
632 finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::
633 dispatch3D( subRegionFE, [&] ( auto const finiteElement )
635 using FE_TYPE = decltype( finiteElement );
638 AverageOverQuadraturePoints1DKernelFactory::
639 createAndLaunch< CellElementSubRegion,
641 parallelDevicePolicy<> >( mesh.getNodeManager(),
642 mesh.getEdgeManager(),
643 mesh.getFaceManager(),
646 meanTotalStressIncrement_k,
647 averageMeanTotalStressIncrement_k );
653 virtual void updateBulkDensity( ElementSubRegionBase & subRegion ) = 0;
655 virtual void validateNonlinearAcceleration()
override
659 GEOS_ERROR(
"Nonlinear acceleration is not implemented for MPI runs" );
#define GEOS_ERROR(msg)
Raise a hard error and terminate the program.
#define GEOS_THROW_IF(EXP, msg, TYPE)
Conditionally throw an exception.
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
virtual void postInputInitialization() override
std::tuple< SOLVERS *... > m_solvers
Pointers of the single-physics solvers.
The DoFManager is responsible for allocating global dofs, constructing sparsity patterns,...
Partition of the decomposed physical domain. It also manages the connexion information to its neighbo...
Group const & getMeshBodies() const
Get the mesh bodies, const version.
The ElementRegionManager class provides an interface to ObjectManagerBase in order to manage ElementR...
void forElementSubRegions(LAMBDA &&lambda)
This function is used to launch kernel function over the element subregions of all the subregion type...
Class facilitating the representation of a multi-level discretization of a MeshBody.
ElementRegionManager const & getElemManager() const
Get the element region manager.
@ Sequential
Sequential coupling.
virtual string getCatalogName() const =0
virtual void setConstitutiveNamesCallSuper(ElementSubRegionBase &subRegion) const
This function sets constitutive name fields on an ElementSubRegionBase, and calls the base function i...
CRSMatrix< real64, globalIndex > m_localMatrix
Local system matrix and rhs.
DofManager m_dofManager
Data structure to handle degrees of freedom.
ParallelVector m_solution
System solution vector.
ParallelVector m_rhs
System right-hand side vector.
array1d< real64 > m_s0
Member variables needed for Nonlinear Acceleration ( Aitken ). Naming convention follows ( Jiang & Tc...
virtual void implicitStepSetup(real64 const &time_n, real64 const &dt, DomainPartition &domain) override
function to perform setup for implicit timestep
real64 m_stabilizationMultiplier
Stabilization Multiplier.
virtual void initializePreSubGroups() override
Called by Initialize() prior to initializing sub-Groups.
stabilization::StabilizationType m_stabilizationType
Type of stabilization used.
virtual void postInputInitialization() override
FLOW_SOLVER * flowSolver() const
accessor for the pointer to the flow solver
PoromechanicsSolver(const string &name, dataRepository::Group *const parent)
main constructor for CoupledSolver Objects
virtual void setConstitutiveNamesCallSuper(ElementSubRegionBase &subRegion) const override final
This function sets constitutive name fields on an ElementSubRegionBase, and calls the base function i...
MECHANICS_SOLVER * solidMechanicsSolver() const
accessor for the pointer to the solid mechanics solver
virtual void mapSolutionBetweenSolvers(DomainPartition &domain, integer const solverType) override
Maps the solution obtained from one solver to the fields used by the other solver(s)
integer m_isThermal
Flag to determine whether or not this is a thermal simulation.
bool m_performStressInitialization
Flag to indicate that the solver is going to perform stress initialization.
string_array m_stabilizationRegionNames
Names of regions where stabilization applied.
virtual void setupDofs(DomainPartition const &domain, DofManager &dofManager) const override
Populate degree-of-freedom manager with fields relevant to this solver.
virtual void registerDataOnMesh(dataRepository::Group &meshBodies) override
Register wrappers that contain data on the mesh objects.
void averageMeanTotalStressIncrement(DomainPartition &domain)
Helper function to average the mean total stress increment over quadrature points.
static string coupledSolverAttributePrefix()
String used to form the solverName used to register solvers in CoupledSolver.
@ QuasiStatic
QuasiStatic.
virtual void initializePreSubGroups()
Called by Initialize() prior to initializing sub-Groups.
Wrapper< TBASE > & registerWrapper(string const &name, wrapperMap::KeyIndex::index_type *const rkey=nullptr)
Create and register a Wrapper around a new object.
string const & getName() const
Get group name.
DataContext const & getWrapperDataContext(KEY key) const
@ NOPLOT
Do not ever write to plot file.
@ OPTIONAL
Optional in input.
@ NO_WRITE
Do not write into restart.
stdVector< string > string_array
A 1-dimensional array of geos::string types.
int MPI_COMM_GEOS
Global MPI communicator used by GEOSX.
array1d< globalIndex > globalIndex_array
A 1-dimensional array of geos::globalIndex types.
std::set< T > set
A set of local indices.
double real64
64-bit floating point type.
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
int integer
Signed integer type.
Array< T, 1 > array1d
Alias for 1D array.
ENUM_STRINGS(LinearSolverParameters::SolverType, "direct", "cg", "gmres", "fgmres", "bicgstab", "richardson", "preconditioner")
Declare strings associated with enumeration values.
Provides enum <-> string conversion facilities.
constexpr static char const * porousMaterialNamesString()
Names of the porous materials.
constexpr static const char * stabilizationRegionNamesString()
Name of regions on which stabilization applied.
static constexpr char const * hydraulicApertureRelationNameString()
Name of the hydraulicApertureRelationName.
constexpr static char const * isThermalString()
Flag to indicate that the simulation is thermal.
constexpr static char const * stabilizationTypeString()
Type of pressure stabilization.
constexpr static const char * stabilizationMultiplierString()
Multiplier on stabilization strength.
Structure to hold scoped key names.