From 5408b7d3abbcb29ab1bc5d0ff125a2ed3b2d6728 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 19 Apr 2024 13:29:51 -0500 Subject: [PATCH] Freeze porosity when doing poromechanics initialization (#3056) --- integratedTests | 2 +- src/coreComponents/constitutive/solid/PorousSolid.hpp | 10 ++++++++++ .../multiphysics/MultiphasePoromechanics.cpp | 2 ++ .../multiphysics/PoromechanicsSolver.hpp | 3 ++- .../multiphysics/SinglePhasePoromechanics.cpp | 2 ++ .../SinglePhasePoromechanicsEmbeddedFractures.hpp | 1 + .../poromechanicsKernels/MultiphasePoromechanics.hpp | 4 ++++ .../MultiphasePoromechanics_impl.hpp | 5 ++++- .../poromechanicsKernels/SinglePhasePoromechanics.hpp | 3 +++ .../SinglePhasePoromechanics_impl.hpp | 5 ++++- .../ThermalMultiphasePoromechanics.hpp | 3 +++ .../ThermalMultiphasePoromechanics_impl.hpp | 3 +++ .../ThermalSinglePhasePoromechanics.hpp | 3 +++ .../ThermalSinglePhasePoromechanics_impl.hpp | 3 +++ 14 files changed, 45 insertions(+), 4 deletions(-) diff --git a/integratedTests b/integratedTests index 671feaaf28d..268498a406a 160000 --- a/integratedTests +++ b/integratedTests @@ -1 +1 @@ -Subproject commit 671feaaf28d51f101144145afa00a6c535669ff1 +Subproject commit 268498a406af2a44fdf7a6952d6703ca7dab30e8 diff --git a/src/coreComponents/constitutive/solid/PorousSolid.hpp b/src/coreComponents/constitutive/solid/PorousSolid.hpp index 2ed3b1fb4a3..078ca36ce7d 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.hpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.hpp @@ -82,6 +82,7 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, real64 ( & dTotalStress_dPressure )[6], real64 ( & dTotalStress_dTemperature )[6], DiscretizationOps & stiffness, + integer const performStressInitialization, real64 & porosity, real64 & porosity_n, real64 & dPorosity_dVolStrain, @@ -116,6 +117,15 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, dPorosity_dPressure, dPorosity_dTemperature ); + // skip porosity update when doing poromechanics initialization + if( performStressInitialization ) + { + porosity = porosityInit; + dPorosity_dVolStrain = 0.0; + dPorosity_dPressure = 0.0; + dPorosity_dTemperature = 0.0; + } + // Save the derivative of solid density wrt pressure for the computation of the body force dSolidDensity_dPressure = m_porosityUpdate.dGrainDensity_dPressure(); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp index ac12d92bbf7..4874ee13e30 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp @@ -224,6 +224,7 @@ void MultiphasePoromechanics< FLOW_SOLVER >::assembleElementBasedTerms( real64 c this->flowSolver()->numFluidComponents(), this->flowSolver()->numFluidPhases(), this->flowSolver()->useTotalMassEquation(), + this->m_performStressInitialization, FlowSolverBase::viewKeyStruct::fluidNamesString() ); } else @@ -242,6 +243,7 @@ void MultiphasePoromechanics< FLOW_SOLVER >::assembleElementBasedTerms( real64 c this->flowSolver()->numFluidPhases(), this->flowSolver()->useSimpleAccumulation(), this->flowSolver()->useTotalMassEquation(), + this->m_performStressInitialization, FlowSolverBase::viewKeyStruct::fluidNamesString() ); } } ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index e0b93cb71bd..6d85889df85 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -389,7 +389,8 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER } /// After the solid mechanics solver - if( solverType == static_cast< integer >( SolverType::SolidMechanics ) ) + if( solverType == static_cast< integer >( SolverType::SolidMechanics ) + && !m_performStressInitialization ) // do not update during poromechanics initialization { // compute the average of the mean total stress increment over quadrature points averageMeanTotalStressIncrement( domain ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp index d21eef612ae..445661f7daf 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp @@ -223,6 +223,7 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementB localRhs, dt, flowDofKey, + this->m_performStressInitialization, FlowSolverBase::viewKeyStruct::fluidNamesString() ); } else @@ -237,6 +238,7 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementB localRhs, dt, flowDofKey, + this->m_performStressInitialization, FlowSolverBase::viewKeyStruct::fluidNamesString() ); } } ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp index e34b3289a7a..5cd29456ab4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -185,6 +185,7 @@ real64 SinglePhasePoromechanicsEmbeddedFractures::assemblyLaunch( MeshLevel & me dt, gravityVectorData, flowDofKey, + m_performStressInitialization, FlowSolverBase::viewKeyStruct::fluidNamesString() ); real64 const maxForce = diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp index 0367fb42b67..c6dabc3a35a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp @@ -102,6 +102,7 @@ class MultiphasePoromechanics : localIndex const numPhases, integer const useSimpleAccumulation, integer const useTotalMassEquation, + integer const performStressInitialization, string const fluidModelKey ); //***************************************************************************** @@ -354,6 +355,8 @@ class MultiphasePoromechanics : /// Use total mass equation flag integer const m_useTotalMassEquation; + + integer const m_performStressInitialization; }; using MultiphasePoromechanicsKernelFactory = @@ -369,6 +372,7 @@ using MultiphasePoromechanicsKernelFactory = localIndex const, integer const, integer const, + integer const, string const >; /** diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp index 51645867ddc..b7104c38fb8 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp @@ -54,6 +54,7 @@ MultiphasePoromechanics( NodeManager const & nodeManager, localIndex const numPhases, integer const useSimpleAccumulation, integer const useTotalMassEquation, + integer const performStressInitialization, string const fluidModelKey ): Base( nodeManager, edgeManager, @@ -79,7 +80,8 @@ MultiphasePoromechanics( NodeManager const & nodeManager, m_numComponents( numComponents ), m_numPhases( numPhases ), m_useSimpleAccumulation( useSimpleAccumulation ), - m_useTotalMassEquation( useTotalMassEquation ) + m_useTotalMassEquation( useTotalMassEquation ), + m_performStressInitialization( performStressInitialization ) { GEOS_ERROR_IF_GT_MSG( m_numComponents, maxNumComponents, "MultiphasePoromechanics solver allows at most " << @@ -162,6 +164,7 @@ smallStrainUpdate( localIndex const k, stack.dTotalStress_dPressure, stack.dTotalStress_dTemperature, stack.stiffness, + m_performStressInitialization, porosity, porosity_n, dPorosity_dVolStrain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp index 08b9c07e5d7..f8902f25f4e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp @@ -94,6 +94,7 @@ class SinglePhasePoromechanics : real64 const inputDt, real64 const (&gravityVector)[3], string const inputFlowDofKey, + integer const performStressInitialization, string const fluidModelKey ); //***************************************************************************** @@ -257,6 +258,7 @@ class SinglePhasePoromechanics : /// Derivative of fluid density wrt pressure arrayView2d< real64 const > const m_dFluidDensity_dPressure; + integer const m_performStressInitialization; }; using SinglePhasePoromechanicsKernelFactory = @@ -268,6 +270,7 @@ using SinglePhasePoromechanicsKernelFactory = real64 const, real64 const (&)[3], string const, + integer const, string const >; /** diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp index 023aeb08bff..e76c112220a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp @@ -50,6 +50,7 @@ SinglePhasePoromechanics( NodeManager const & nodeManager, real64 const inputDt, real64 const (&gravityVector)[3], string const inputFlowDofKey, + integer const performStressInitialization, string const fluidModelKey ): Base( nodeManager, edgeManager, @@ -68,7 +69,8 @@ SinglePhasePoromechanics( NodeManager const & nodeManager, fluidModelKey ), m_fluidDensity( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ), m_fluidDensity_n( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ), - m_dFluidDensity_dPressure( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() ) + m_dFluidDensity_dPressure( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() ), + m_performStressInitialization( performStressInitialization ) {} template< typename SUBREGION_TYPE, @@ -100,6 +102,7 @@ smallStrainUpdate( localIndex const k, stack.dTotalStress_dPressure, stack.dTotalStress_dTemperature, stack.stiffness, + m_performStressInitialization, porosity, porosity_n, dPorosity_dVolStrain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp index 5c60e6820da..62833506841 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp @@ -76,6 +76,7 @@ class ThermalMultiphasePoromechanics : using Base::m_numComponents; using Base::m_numPhases; using Base::m_useTotalMassEquation; + using Base::m_performStressInitialization; using Base::m_solidDensity; using Base::m_fluidPhaseMassDensity; using Base::m_dFluidPhaseMassDensity; @@ -110,6 +111,7 @@ class ThermalMultiphasePoromechanics : localIndex const numComponents, localIndex const numPhases, integer const useTotalMassEquation, + integer const performStressInitialization, string const fluidModelKey ); /** @@ -341,6 +343,7 @@ using ThermalMultiphasePoromechanicsKernelFactory = localIndex const, localIndex const, integer const, + integer const, string const >; } // namespace thermalporomechanicsKernels diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp index 0c19141489f..fe18bd5f75c 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp @@ -52,6 +52,7 @@ ThermalMultiphasePoromechanics( NodeManager const & nodeManager, localIndex const numComponents, localIndex const numPhases, integer const useTotalMassEquation, + integer const performStressInitialization, string const fluidModelKey ): Base( nodeManager, edgeManager, @@ -71,6 +72,7 @@ ThermalMultiphasePoromechanics( NodeManager const & nodeManager, numPhases, false, // do not use simple accumulation form useTotalMassEquation, + performStressInitialization, fluidModelKey ), m_rockInternalEnergy_n( inputConstitutiveType.getInternalEnergy_n() ), m_rockInternalEnergy( inputConstitutiveType.getInternalEnergy() ), @@ -137,6 +139,7 @@ smallStrainUpdate( localIndex const k, stack.dTotalStress_dPressure, stack.dTotalStress_dTemperature, stack.stiffness, + m_performStressInitialization, porosity, porosity_n, dPorosity_dVolStrain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp index 5fe11040d59..178da2ea3bc 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp @@ -70,6 +70,7 @@ class ThermalSinglePhasePoromechanics : using Base::m_solidDensity; using Base::m_flowDofNumber; using Base::m_dt; + using Base::m_performStressInitialization; /** * @brief Constructor @@ -90,6 +91,7 @@ class ThermalSinglePhasePoromechanics : real64 const inputDt, real64 const (&gravityVector)[3], string const inputFlowDofKey, + integer const performStressInitialization, string const fluidModelKey ); //***************************************************************************** @@ -305,6 +307,7 @@ using ThermalSinglePhasePoromechanicsKernelFactory = real64 const, real64 const (&)[3], string const, + integer const, string const >; } // namespace thermalPoromechanicsKernels diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp index a3f3d4db3a9..76584fb52dd 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp @@ -46,6 +46,7 @@ ThermalSinglePhasePoromechanics( NodeManager const & nodeManager, real64 const inputDt, real64 const (&gravityVector)[3], string const inputFlowDofKey, + integer const performStressInitialization, string const fluidModelKey ): Base( nodeManager, edgeManager, @@ -61,6 +62,7 @@ ThermalSinglePhasePoromechanics( NodeManager const & nodeManager, inputDt, gravityVector, inputFlowDofKey, + performStressInitialization, fluidModelKey ), m_dFluidDensity_dTemperature( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dTemperature() ), @@ -121,6 +123,7 @@ smallStrainUpdate( localIndex const k, stack.dTotalStress_dPressure, stack.dTotalStress_dTemperature, stack.stiffness, + m_performStressInitialization, porosity, porosity_n, dPorosity_dVolStrain,