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",
122 this->
template setConstitutiveName< constitutive::CoupledSolidBase >( subRegion,
126 this->
template setConstitutiveName< constitutive::PorosityBase >( subRegion,
127 constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString(),
"porosity" );
131 this->
template setConstitutiveName< constitutive::HydraulicApertureBase >( subRegion,
142 ": Local stabilization has been temporarily disabled",
145 DomainPartition & domain = this->
template getGroupByPath< DomainPartition >(
"/Problem/domain" );
147 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&] (
string const &,
151 ElementRegionManager & elementRegionManager = mesh.getElemManager();
152 elementRegionManager.forElementSubRegions< CellElementSubRegion >( regionNames,
153 [&]( localIndex const,
154 ElementSubRegionBase & subRegion )
156 if( subRegion.hasField< fields::poromechanics::bulkDensity >() )
159 constitutive::CoupledSolidBase const & solid =
160 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
161 subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
162 subRegion.getField< fields::poromechanics::bulkDensity >().resizeDimension< 1 >( solid.getDensity().size( 1 ) );
170 Base::registerDataOnMesh( meshBodies );
175 solidMechanicsSolver()->enableFixedStressPoromechanicsUpdate();
177 flowSolver()->enableFixedStressPoromechanicsUpdate();
180 if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
182 flowSolver()->enableJumpStabilization();
185 this->
template forDiscretizationOnMeshTargets( meshBodies, [&] (
string const &,
195 subRegion.registerWrapper<
string >( viewKeyStruct::porousMaterialNamesString() ).
198 setSizedFromParent( 0 );
201 subRegion.registerWrapper<
string >( constitutive::CoupledSolidBase::viewKeyStruct::porosityModelNameString() ).
204 setSizedFromParent( 0 );
208 subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() );
210 if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
212 subRegion.registerField< fields::flow::macroElementIndex >( this->getName());
213 subRegion.registerField< fields::flow::elementStabConstant >( this->getName());
223 flowSolver()->setKeepVariablesConstantDuringInitStep( m_performStressInitialization );
225 if( this->m_stabilizationType == stabilization::StabilizationType::Global || this->m_stabilizationType == stabilization::StabilizationType::Local )
227 this->updateStabilizationParameters( domain );
230 Base::implicitStepSetup( time_n, dt, domain );
239 solidMechanicsSolver()->setupDofs( domain, dofManager );
240 flowSolver()->setupDofs( domain, dofManager );
241 this->setupCoupling( domain, dofManager );
244 virtual bool checkSequentialConvergence(
int const & iter,
250 auto & subcycling = this->getNonlinearSolverParameters().m_subcyclingOption;
251 auto const subcycling_orig = subcycling;
252 if( m_performStressInitialization )
255 bool isConverged = Base::checkSequentialConvergence( iter, time_n, dt, domain );
258 subcycling = subcycling_orig;
269 return std::get< toUnderlying( SolverType::SolidMechanics ) >( m_solvers );
278 return std::get< toUnderlying( SolverType::Flow ) >( m_solvers );
285 void setStressInitialization(
bool const performStressInitialization )
287 m_performStressInitialization = performStressInitialization;
288 solidMechanicsSolver()->setStressInitialization( performStressInitialization );
317 for(
string const & regionName : m_stabilizationRegionNames )
319 regionFilter.insert( regionName );
323 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&] (
string const &,
328 string_array filteredTargetRegionNames;
329 filteredTargetRegionNames.reserve( targetRegionNames.size() );
331 for( string const & targetRegionName : targetRegionNames )
333 if( regionFilter.count( targetRegionName ) )
335 filteredTargetRegionNames.emplace_back( targetRegionName );
340 mesh.getElemManager().forElementSubRegions( filteredTargetRegionNames, [&](
localIndex const,
341 ElementSubRegionBase & subRegion )
343 arrayView1d< integer > const macroElementIndex = subRegion.getField< fields::flow::macroElementIndex >();
344 arrayView1d< real64 > const elementStabConstant = subRegion.getField< fields::flow::elementStabConstant >();
346 constitutive::CoupledSolidBase const & porousSolid =
347 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion,
348 subRegion.getReference< string >( viewKeyStruct::porousMaterialNamesString() ) );
350 arrayView1d< real64 const > const bulkModulus = porousSolid.getBulkModulus();
351 arrayView1d< real64 const > const shearModulus = porousSolid.getShearModulus();
352 arrayView1d< real64 const > const biotCoefficient = porousSolid.getBiotCoefficient();
354 real64 const stabilizationMultiplier = m_stabilizationMultiplier;
356 forAll< parallelDevicePolicy<> >( subRegion.size(), [bulkModulus,
359 stabilizationMultiplier,
361 elementStabConstant] GEOS_HOST_DEVICE ( localIndex const ei )
364 real64 const bM = bulkModulus[ei];
365 real64 const sM = shearModulus[ei];
366 real64 const bC = biotCoefficient[ei];
368 macroElementIndex[ei] = 1;
369 elementStabConstant[ei] = stabilizationMultiplier * 9.0 * (bC * bC) / (32.0 * (10.0 * sM / 3.0 + bM));
376 void updateBulkDensity( DomainPartition & domain )
378 this->
template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&](
string const &,
382 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
388 updateBulkDensity( subRegion );
395 template<
typename TYPE_LIST,
396 typename KERNEL_WRAPPER,
397 typename ... PARAMS >
398 real64 assemblyLaunch( MeshLevel & mesh,
399 DofManager
const & dofManager,
401 string const & materialNamesString,
402 CRSMatrixView< real64, globalIndex const >
const & localMatrix,
403 arrayView1d< real64 >
const & localRhs,
405 PARAMS && ... params )
409 NodeManager
const & nodeManager = mesh.getNodeManager();
411 string const dofKey = dofManager.getKey( fields::solidMechanics::totalDisplacement::key() );
412 arrayView1d< globalIndex const >
const & dofNumber = nodeManager.getReference<
globalIndex_array >( dofKey );
414 real64 const gravityVectorData[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( this->gravityVector() );
416 KERNEL_WRAPPER kernelWrapper( dofNumber,
417 dofManager.rankOffset(),
422 std::forward< PARAMS >( params )... );
424 return finiteElement::
425 regionBasedKernelApplication< parallelDevicePolicy< >,
428 this->solidMechanicsSolver()->getDiscretizationName(),
435 void recordAverageMeanTotalStressIncrement( DomainPartition & domain,
436 array1d< real64 > & averageMeanTotalStressIncrement )
438 averageMeanTotalStressIncrement.resize( 0 );
439 this->
template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&](
string const &,
443 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
447 string const & solidName =
448 subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
449 constitutive::CoupledSolidBase & solid =
450 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
452 arrayView1d< const real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
453 for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
455 averageMeanTotalStressIncrement.emplace_back( averageMeanTotalStressIncrement_k[k] );
461 void applyAcceleratedAverageMeanTotalStressIncrement( DomainPartition & domain,
462 array1d< real64 > & averageMeanTotalStressIncrement )
465 this->
template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&](
string const &,
469 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
473 string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
474 constitutive::CoupledSolidBase & solid =
475 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
476 auto & porosityModel = dynamic_cast< constitutive::BiotPorosity const & >( solid.getBasePorosityModel());
477 arrayView1d< real64 > const & averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
478 for( localIndex k = 0; k < localIndex( averageMeanTotalStressIncrement_k.size()); k++ )
480 porosityModel.updateAverageMeanTotalStressIncrement( k, averageMeanTotalStressIncrement[i] );
487 real64 computeAitkenRelaxationFactor( array1d< real64 >
const & s0,
488 array1d< real64 >
const & s1,
489 array1d< real64 >
const & s1_tilde,
490 array1d< real64 >
const & s2_tilde,
493 array1d< real64 > r1 = axpy( s1_tilde, s0, -1.0 );
494 array1d< real64 > r2 = axpy( s2_tilde, s1, -1.0 );
497 array1d< real64 > diff = axpy( r2, r1, -1.0 );
499 real64 const denom = dot( diff, diff );
500 real64 const numer = dot( r1, diff );
503 if( !isZero( denom ))
505 omega1 = -1.0 * omega0 * numer / denom;
510 array1d< real64 > computeUpdate( array1d< real64 >
const & s1,
511 array1d< real64 >
const & s2_tilde,
514 return axpy( scale( s1, 1.0 - omega1 ),
515 scale( s2_tilde, omega1 ),
519 void startSequentialIteration(
integer const & iter,
520 DomainPartition & domain )
override
522 if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
526 recordAverageMeanTotalStressIncrement( domain, m_s1 );
532 m_s1_tilde = m_s2_tilde;
538 void finishSequentialIteration(
integer const & iter,
539 DomainPartition & domain )
override
541 if( this->getNonlinearSolverParameters().m_nonlinearAccelerationType == NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
550 m_omega1 = computeAitkenRelaxationFactor( m_s0, m_s1, m_s1_tilde, m_s2_tilde, m_omega0 );
551 m_s2 = computeUpdate( m_s1, m_s2_tilde, m_omega1 );
552 applyAcceleratedAverageMeanTotalStressIncrement( domain, m_s2 );
562 if( solverType ==
static_cast< integer >( SolverType::Flow ) )
564 updateBulkDensity( domain );
568 if( solverType ==
static_cast< integer >( SolverType::SolidMechanics )
569 && !m_performStressInitialization )
572 averageMeanTotalStressIncrement( domain );
574 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&](
string const &,
579 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
584 flowSolver()->updatePorosityAndPermeability( subRegion );
586 updateBulkDensity( subRegion );
592 if( solverType ==
static_cast< integer >( SolverType::SolidMechanics ) &&
593 this->getNonlinearSolverParameters().m_nonlinearAccelerationType== NonlinearSolverParameters::NonlinearAccelerationType::Aitken )
595 recordAverageMeanTotalStressIncrement( domain, m_s2_tilde );
605 this->
template forDiscretizationOnMeshTargets<>( domain.
getMeshBodies(), [&](
string const &,
609 mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
613 string const & solidName = subRegion.template getReference< string >( viewKeyStruct::porousMaterialNamesString() );
614 constitutive::CoupledSolidBase & solid =
615 this->template getConstitutiveModel< constitutive::CoupledSolidBase >( subRegion, solidName );
617 arrayView2d< real64 const > const meanTotalStressIncrement_k = solid.getMeanTotalStressIncrement_k();
618 arrayView1d< real64 > const averageMeanTotalStressIncrement_k = solid.getAverageMeanTotalStressIncrement_k();
620 finiteElement::FiniteElementBase & subRegionFE =
621 subRegion.template getReference< finiteElement::FiniteElementBase >( solidMechanicsSolver()->getDiscretizationName() );
624 finiteElement::FiniteElementDispatchHandler< BASE_FE_TYPES >::
625 dispatch3D( subRegionFE, [&] ( auto const finiteElement )
627 using FE_TYPE = decltype( finiteElement );
630 AverageOverQuadraturePoints1DKernelFactory::
631 createAndLaunch< CellElementSubRegion,
633 parallelDevicePolicy<> >( mesh.getNodeManager(),
634 mesh.getEdgeManager(),
635 mesh.getFaceManager(),
638 meanTotalStressIncrement_k,
639 averageMeanTotalStressIncrement_k );
645 virtual void updateBulkDensity( ElementSubRegionBase & subRegion ) = 0;
647 virtual void validateNonlinearAcceleration()
override
651 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.
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.
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 initializePostInitialConditionsPreSubGroups() override
Called by InitializePostInitialConditions() prior to initializing sub-Groups.
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.
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
virtual void initializePostInitialConditionsPreSubGroups()
Called by InitializePostInitialConditions() prior to initializing sub-Groups.
@ 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.
ENUM_STRINGS(LinearSolverParameters::SolverType, "direct", "cg", "gmres", "fgmres", "bicgstab", "preconditioner")
Declare strings associated with enumeration values.
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).
std::int32_t integer
Signed integer type.
@ None
No Schur complement - just block-GS/block-Jacobi preconditioner.
Array< T, 1 > array1d
Alias for 1D array.
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.