diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index c5de93442f2..f303a54d8d7 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr2125-10859-1db8623 + baseline: integratedTests/baseline_integratedTests-pr3396-10875-13b594f allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 26c0bdfc084..6a66008eb18 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3396 (2024-03-21) +===================== +Use solid mechanics solver directly to perform poromechanics initialization. + PR #2125 (2024-03-20) ===================== Phase-field nucleation model. diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 5e71b0d87ed..d52bb3716d4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -168,6 +168,8 @@ class FlowSolverBase : public PhysicsSolverBase GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." ); } + void initializeState( DomainPartition & domain ); + virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } virtual void initializeThermalState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); } @@ -236,8 +238,6 @@ class FlowSolverBase : public PhysicsSolverBase virtual void initializePostInitialConditionsPreSubGroups() override; - void initializeState( DomainPartition & domain ); - virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); } void initializePorosityAndPermeability( MeshLevel & mesh, string_array const & regionNames ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index 9f6ebe22d02..4cdb07fad86 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -231,6 +231,8 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL } } + void initializeState( DomainPartition & domain ) const { return reservoirSolver()->initializeState( domain ); } + void assembleFluxTerms( real64 const dt, DomainPartition const & domain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp index 30259c11b43..922c27319d6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp @@ -41,7 +41,7 @@ class MultiphasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHANI using Base::m_stabilizationType; using Base::m_stabilizationRegionNames; using Base::m_stabilizationMultiplier; - using Base::getLogLevel; + using Base::updateBulkDensity; /** * @brief main constructor for MultiphasePoromechanics Objects diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index 955b745bcfd..6791bee3f49 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -21,7 +21,6 @@ #include "events/tasks/TasksManager.hpp" #include "physicsSolvers/PhysicsSolverManager.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsStatistics.hpp" #include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp" #include "physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp" @@ -87,8 +86,9 @@ postInputInitialization() TasksManager & tasksManager = problemManager.getGroup< TasksManager >( "Tasks" ); GEOS_THROW_IF( !tasksManager.hasGroup( m_solidMechanicsStatisticsName ), - GEOS_FMT( "{}: statistics task named {} not found", + GEOS_FMT( "{}: {} task named {} not found", getWrapperDataContext( viewKeyStruct::solidMechanicsStatisticsNameString() ), + SolidMechanicsStatistics::catalogName(), m_solidMechanicsStatisticsName ), InputError ); @@ -116,7 +116,9 @@ execute( real64 const time_n, m_solidMechanicsStateResetTask.execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); - m_poromechanicsSolver->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); + m_poromechanicsSolver->flowSolver()->initializeState( domain ); + m_poromechanicsSolver->updateBulkDensity( domain ); + m_poromechanicsSolver->solidMechanicsSolver()->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain ); GEOS_LOG_LEVEL_RANK_0( logInfo::SolverInitialization, GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` has completed stress initialization",//1 getName(), time_n + dt, m_poromechanicsSolverName ) ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp index 60e8535a2b7..9a6d06ed454 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp @@ -91,8 +91,6 @@ class PoromechanicsInitialization : public TaskBase void postInputInitialization() override; -// void registerDataOnMesh( Group & meshBodies ) override; - /// Name of the poromechanics solver string m_poromechanicsSolverName; diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 099f3fb41c9..8c19b341664 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -81,18 +81,14 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER PoromechanicsSolver( const string & name, dataRepository::Group * const parent ) : Base( name, parent ), - m_isThermal( 0 ) + m_isThermal( 0 ), + m_performStressInitialization( false ) { this->registerWrapper( viewKeyStruct::isThermalString(), &m_isThermal ). setApplyDefaultValue( 0 ). setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "Flag indicating whether the problem is thermal or not. Set isThermal=\"1\" to enable the thermal coupling" ); - this->registerWrapper( viewKeyStruct::performStressInitializationString(), &m_performStressInitialization ). - setApplyDefaultValue( false ). - setInputFlag( dataRepository::InputFlags::FALSE ). - setDescription( "Flag to indicate that the solver is going to perform stress initialization" ); - this->registerWrapper( viewKeyStruct::stabilizationTypeString(), &m_stabilizationType ). setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "StabilizationType. Options are:\n" + @@ -220,12 +216,9 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER setRestartFlags( dataRepository::RestartFlags::NO_WRITE ). setSizedFromParent( 0 ); - if( this->getNonlinearSolverParameters().m_couplingType == NonlinearSolverParameters::CouplingType::Sequential ) - { - // register the bulk density for use in the solid mechanics solver - // ideally we would resize it here as well, but the solid model name is not available yet (see below) - subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() ); - } + // register the bulk density for use in the solid mechanics solver + // ideally we would resize it here as well, but the solid model name is not available yet (see below) + subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() ); if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local ) { @@ -302,9 +295,10 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER * @brief Utility function to set the stress initialization flag * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise */ - void setStressInitialization( integer const performStressInitialization ) + void setStressInitialization( bool const performStressInitialization ) { m_performStressInitialization = performStressInitialization; + solidMechanicsSolver()->setStressInitialization( performStressInitialization ); } struct viewKeyStruct : Base::viewKeyStruct @@ -315,9 +309,6 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER /// Flag to indicate that the simulation is thermal constexpr static char const * isThermalString() { return "isThermal"; } - /// Flag to indicate that the solver is going to perform stress initialization - constexpr static char const * performStressInitializationString() { return "performStressInitialization"; } - /// Type of pressure stabilization constexpr static char const * stabilizationTypeString() {return "stabilizationType"; } @@ -394,6 +385,23 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER } + void updateBulkDensity( DomainPartition & domain ) + { + this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // update the bulk density + // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the + // pressure/temperature terms + updateBulkDensity( subRegion ); + } ); + } ); + } + protected: template< typename CONSTITUTIVE_BASE, @@ -561,20 +569,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER /// After the flow solver if( solverType == static_cast< integer >( SolverType::Flow ) ) { - this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - - mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, - auto & subRegion ) - { - // update the bulk density - // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the - // pressure/temperature terms - updateBulkDensity( subRegion ); - } ); - } ); + updateBulkDensity( domain ); } /// After the solid mechanics solver @@ -668,7 +663,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER integer m_isThermal; /// Flag to indicate that the solver is going to perform stress initialization - integer m_performStressInitialization; + bool m_performStressInitialization; /// Type of stabilization used stabilization::StabilizationType m_stabilizationType; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index 553740546f9..66d12a7da9c 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -41,6 +41,7 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN using Base::m_stabilizationType; using Base::m_stabilizationRegionNames; using Base::m_stabilizationMultiplier; + using Base::updateBulkDensity; /** * @brief main constructor for SinglePhasePoromechanics objects diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index dd1623ab5dd..d5b672efeaa 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -60,7 +60,8 @@ SolidMechanicsLagrangianFEM::SolidMechanicsLagrangianFEM( const string & name, m_maxForce( 0.0 ), m_maxNumResolves( 10 ), m_strainTheory( 0 ), - m_isFixedStressPoromechanicsUpdate( false ) + m_isFixedStressPoromechanicsUpdate( false ), + m_performStressInitialization( false ) { registerWrapper( viewKeyStruct::newmarkGammaString(), &m_newmarkGamma ). @@ -1006,6 +1007,8 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain, ParallelVector & solution, bool const setSparsity ) { + GEOS_LOG( "SolidMechanicsLagrangianFEM::setupSystem" ); + GEOS_MARK_FUNCTION; PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); @@ -1021,25 +1024,6 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain, arrayView1d< globalIndex const > const dofNumber = nodeManager.getReference< globalIndex_array >( dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); - if( m_contactRelationName != viewKeyStruct::noContactRelationNameString() ) - { - ElementRegionManager const & elemManager = mesh.getElemManager(); - string_array allFaceElementRegions; - elemManager.forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion const & elemRegion ) - { - allFaceElementRegions.emplace_back( elemRegion.getName() ); - } ); - - finiteElement:: - fillSparsity< FaceElementSubRegion, - solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, - allFaceElementRegions, - this->getDiscretizationName(), - dofNumber, - dofManager.rankOffset(), - sparsityPattern ); - - } finiteElement:: fillSparsity< CellElementSubRegion, solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, @@ -1072,7 +1056,7 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM MeshLevel & mesh, string_array const & regionNames ) { - if( m_isFixedStressPoromechanicsUpdate ) + if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization ) { set< string > poromechanicsRegions; set< string > mechanicsRegions; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index 5d783ac67a9..56b78476aa1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -278,6 +278,15 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase return m_rigidBodyModes; } + /* + * @brief Utility function to set the stress initialization flag + * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise + */ + void setStressInitialization( bool const performStressInitialization ) + { + m_performStressInitialization = performStressInitialization; + } + protected: virtual void postInputInitialization() override; @@ -293,8 +302,11 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase real64 m_maxForce = 0.0; integer m_maxNumResolves; integer m_strainTheory; -// MPI_iCommData m_iComm; + + /// Flag to indicate that the solver is running in fixed stress (sequential) mode bool m_isFixedStressPoromechanicsUpdate; + /// Flag to indicate that the solver is going to perform stress initialization + bool m_performStressInitialization; /// Rigid body modes array1d< ParallelVector > m_rigidBodyModes; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index e4045f4ef5c..aac23546d9c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -640,7 +640,7 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time, assembleContact( domain, dofManager, localMatrix, localRhs ); // for sequential: add (fixed) pressure force contribution into residual (no derivatives) - if( m_isFixedStressPoromechanicsUpdate ) + if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization ) { forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh,