From b652c037ebf60d7bee086048a3fc0f6d5d2261d8 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Mon, 9 Feb 2026 14:40:34 -0600 Subject: [PATCH 1/2] fix seq coupling for multiphase poromechanics --- .../multiphysics/MultiphasePoromechanics.hpp | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp index de5de39883a..e74a6c066dd 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp @@ -107,6 +107,31 @@ class MultiphasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHANI GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); } + virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override + { + GEOS_MARK_FUNCTION; + + Base::mapSolutionBetweenSolvers( domain, solverType ); + + /// After the solid mechanics solver + if( solverType == static_cast< integer >( Base::SolverType::SolidMechanics ) + && !this->m_performStressInitialization ) // do not update during poromechanics initialization + { + this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // update mass after porosity change due to mechanics solve + this->flowSolver()->updateMass( subRegion ); + } ); + } ); + } + } + /** * @brief Helper function to recompute the bulk density * @param[in] subRegion the element subRegion From f08d67d19536291e38e0ed0e31b8d2fe9c142c60 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 10 Feb 2026 12:04:34 +0100 Subject: [PATCH 2/2] discard Well / test update compAmount --- .../fluidFlow/CompositionalMultiphaseBase.cpp | 9 +- .../multiphysics/MultiphasePoromechanics.cpp | 99 ++++++++++--------- .../multiphysics/MultiphasePoromechanics.hpp | 4 +- ...iphasePoromechanicsConformingFractures.cpp | 6 +- ...asePoromechanicsConformingFracturesALM.cpp | 6 +- .../PoromechanicsInitialization.cpp | 16 +-- 6 files changed, 78 insertions(+), 62 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index ec8268e0079..aee38af565f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -780,6 +780,7 @@ void CompositionalMultiphaseBase::updateCompAmount( ElementSubRegionBase & subRe arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< flow::compAmount >(); arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView1d< real64 const > const deltaVolume = subRegion.getField< fields::flow::deltaVolume >(); string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() ); CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName ); arrayView2d< real64 const > const porosity = porousMaterial.getPorosity(); @@ -799,7 +800,7 @@ void CompositionalMultiphaseBase::updateCompAmount( ElementSubRegionBase & subRe for( integer ic = 0; ic < numComp; ++ic ) { // m = phi*V*z_c*rho_T - compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compFrac[ei][ic] * totalDens[ei][0]; + compAmount[ei][ic] = porosity[ei][0] * (volume[ei]+deltaVolume[ei]) * compFrac[ei][ic] * totalDens[ei][0]; } } ); } @@ -811,7 +812,11 @@ void CompositionalMultiphaseBase::updateCompAmount( ElementSubRegionBase & subRe { for( integer ic = 0; ic < numComp; ++ic ) { - compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compDens[ei][ic]; + compAmount[ei][ic] = porosity[ei][0] * (volume[ei]+deltaVolume[ei]) * compDens[ei][ic]; + GEOS_LOG_RANK(GEOS_FMT("Attempted Correction @{}/{} : {} -> {} ", ei, ic, + porosity[ei][0]*volume[ei]* compDens[ei][ic], + porosity[ei][0]*deltaVolume[ei]* compDens[ei][ic] + )); } } ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp index 545ebfcfc43..6f57204d32e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp @@ -49,22 +49,31 @@ MultiphasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::MultiphasePoromechanic : Base( name, parent ) { } +// template<> +// void MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsLagrangianFEM >::assembleSystem( real64 const time, +// real64 const dt, +// DomainPartition & domain, +// DofManager const & dofManager, +// CRSMatrixView< real64, globalIndex const > const & localMatrix, +// arrayView1d< real64 > const & localRhs ) +// { +// GEOS_MARK_FUNCTION; + +// Base::assembleSystem( time, dt, domain, dofManager, localMatrix, localRhs ); + +// // assemble well contributions +// flowSolver()->wellSolver()->assembleSystem( time, dt, domain, dofManager, localMatrix, localRhs ); +// flowSolver()->assembleCouplingTerms( time, dt, domain, dofManager, localMatrix, localRhs ); +// } template<> -void MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsLagrangianFEM >::assembleSystem( real64 const time, - real64 const dt, - DomainPartition & domain, - DofManager const & dofManager, - CRSMatrixView< real64, globalIndex const > const & localMatrix, - arrayView1d< real64 > const & localRhs ) -{ - GEOS_MARK_FUNCTION; - - Base::assembleSystem( time, dt, domain, dofManager, localMatrix, localRhs ); +void MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsLagrangeContact >::assembleElementBasedTerms( real64 const time_n, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{} - // assemble well contributions - flowSolver()->wellSolver()->assembleSystem( time, dt, domain, dofManager, localMatrix, localRhs ); - flowSolver()->assembleCouplingTerms( time, dt, domain, dofManager, localMatrix, localRhs ); -} template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > void MultiphasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementBasedTerms( real64 const time_n, @@ -204,32 +213,32 @@ void MultiphasePoromechanics<>::setMGRStrategy() EnumStrings< LinearSolverParameters::MGR::StrategyType >::toString( linearSolverParameters.mgr.strategy ))); } -template<> -void MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsLagrangianFEM >::setMGRStrategy() -{ - // same as CompositionalMultiphaseReservoirAndWells< MultiphasePoromechanics<> >::setMGRStrategy - - LinearSolverParameters & linearSolverParameters = m_linearSolverParameters.get(); - - if( linearSolverParameters.preconditionerType != LinearSolverParameters::PreconditionerType::mgr ) - return; - - linearSolverParameters.mgr.separateComponents = true; - linearSolverParameters.dofsPerNode = 3; - - if( dynamic_cast< CompositionalMultiphaseHybridFVM * >( this->flowSolver() ) ) - { - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}/{}", - this->getName(), this->getCatalogName(), this->flowSolver()->getCatalogName() ) ); - } - else - { - linearSolverParameters.mgr.strategy = LinearSolverParameters::MGR::StrategyType::multiphasePoromechanicsReservoirFVM; - } - GEOS_LOG_LEVEL_RANK_0( logInfo::LinearSolver, - GEOS_FMT( "{}: MGR strategy set to {}", getName(), - EnumStrings< LinearSolverParameters::MGR::StrategyType >::toString( linearSolverParameters.mgr.strategy ))); -} +// template<> +// void MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsLagrangianFEM >::setMGRStrategy() +// { +// // same as CompositionalMultiphaseReservoirAndWells< MultiphasePoromechanics<> >::setMGRStrategy + +// LinearSolverParameters & linearSolverParameters = m_linearSolverParameters.get(); + +// if( linearSolverParameters.preconditionerType != LinearSolverParameters::PreconditionerType::mgr ) +// return; + +// linearSolverParameters.mgr.separateComponents = true; +// linearSolverParameters.dofsPerNode = 3; + +// if( dynamic_cast< CompositionalMultiphaseHybridFVM * >( this->flowSolver() ) ) +// { +// GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}/{}", +// this->getName(), this->getCatalogName(), this->flowSolver()->getCatalogName() ) ); +// } +// else +// { +// linearSolverParameters.mgr.strategy = LinearSolverParameters::MGR::StrategyType::multiphasePoromechanicsReservoirFVM; +// } +// GEOS_LOG_LEVEL_RANK_0( logInfo::LinearSolver, +// GEOS_FMT( "{}: MGR strategy set to {}", getName(), +// EnumStrings< LinearSolverParameters::MGR::StrategyType >::toString( linearSolverParameters.mgr.strategy ))); +// } template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > void MultiphasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::updateBulkDensity( ElementSubRegionBase & subRegion ) @@ -254,14 +263,14 @@ void MultiphasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::updateBulkDensity template class MultiphasePoromechanics<>; template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsLagrangeContact >; template class MultiphasePoromechanics< CompositionalMultiphaseBase, SolidMechanicsAugmentedLagrangianContact >; -template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> >; -template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsLagrangeContact >; -template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsAugmentedLagrangianContact >; +// template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> >; +// template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsLagrangeContact >; +// template class MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<>, SolidMechanicsAugmentedLagrangianContact >; namespace { -typedef MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > MultiphaseReservoirPoromechanics; -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphaseReservoirPoromechanics, string const &, Group * const ) +// typedef MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > MultiphaseReservoirPoromechanics; +// REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphaseReservoirPoromechanics, string const &, Group * const ) typedef MultiphasePoromechanics<> MultiphasePoromechanics; REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphasePoromechanics, string const &, Group * const ) } diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp index e74a6c066dd..9977fcc3973 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp @@ -126,7 +126,9 @@ class MultiphasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHANI auto & subRegion ) { // update mass after porosity change due to mechanics solve - this->flowSolver()->updateMass( subRegion ); + GEOS_LOG_RANK("------ START ---------------------------------------------------"); + this->flowSolver()->updateCompAmount( subRegion ); + GEOS_LOG_RANK("------ END ---------------------------------------------------"); } ); } ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp index 0ef1cb73f36..b3e1aadb9c6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp @@ -267,12 +267,12 @@ assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, } template class MultiphasePoromechanicsConformingFractures<>; -template class MultiphasePoromechanicsConformingFractures< CompositionalMultiphaseReservoirAndWells<> >; +// template class MultiphasePoromechanicsConformingFractures< CompositionalMultiphaseReservoirAndWells<> >; namespace { -typedef MultiphasePoromechanicsConformingFractures< CompositionalMultiphaseReservoirAndWells<> > MultiphaseReservoirPoromechanicsConformingFractures; -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphaseReservoirPoromechanicsConformingFractures, string const &, Group * const ) +// typedef MultiphasePoromechanicsConformingFractures< CompositionalMultiphaseReservoirAndWells<> > MultiphaseReservoirPoromechanicsConformingFractures; +// REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphaseReservoirPoromechanicsConformingFractures, string const &, Group * const ) typedef MultiphasePoromechanicsConformingFractures<> MultiphasePoromechanicsConformingFractures; REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphasePoromechanicsConformingFractures, string const &, Group * const ) } diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp index ef17e087d1d..54166073786 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp @@ -175,12 +175,12 @@ assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & mesh, template class MultiphasePoromechanicsConformingFracturesALM<>; -template class MultiphasePoromechanicsConformingFracturesALM< CompositionalMultiphaseReservoirAndWells<> >; +// template class MultiphasePoromechanicsConformingFracturesALM< CompositionalMultiphaseReservoirAndWells<> >; namespace { -typedef MultiphasePoromechanicsConformingFracturesALM< CompositionalMultiphaseReservoirAndWells<> > MultiphaseReservoirPoromechanicsConformingFracturesALM; -REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphaseReservoirPoromechanicsConformingFracturesALM, string const &, Group * const ) +// typedef MultiphasePoromechanicsConformingFracturesALM< CompositionalMultiphaseReservoirAndWells<> > MultiphaseReservoirPoromechanicsConformingFracturesALM; +// REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphaseReservoirPoromechanicsConformingFracturesALM, string const &, Group * const ) typedef MultiphasePoromechanicsConformingFracturesALM<> MultiphasePoromechanicsConformingFracturesALM; REGISTER_CATALOG_ENTRY( PhysicsSolverBase, MultiphasePoromechanicsConformingFracturesALM, string const &, Group * const ) } diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index 288e282dae1..7ebecb32add 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -141,11 +141,11 @@ namespace typedef PoromechanicsInitialization< MultiphasePoromechanics<> > MultiphasePoromechanicsInitialization; typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFractures<> > MultiphasePoromechanicsConformingFracturesInitialization; typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFracturesALM<> > MultiphasePoromechanicsConformingFracturesALMInitialization; -typedef PoromechanicsInitialization< MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > > MultiphaseReservoirPoromechanicsInitialization; -typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFractures< CompositionalMultiphaseReservoirAndWells<> > > - CompositionalMultiphaseReservoirPoromechanicsConformingFracturesInitialization; -typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFracturesALM< CompositionalMultiphaseReservoirAndWells<> > > - CompositionalMultiphaseReservoirPoromechanicsConformingFracturesALMInitialization; +// typedef PoromechanicsInitialization< MultiphasePoromechanics< CompositionalMultiphaseReservoirAndWells<> > > MultiphaseReservoirPoromechanicsInitialization; +// typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFractures< CompositionalMultiphaseReservoirAndWells<> > > + // CompositionalMultiphaseReservoirPoromechanicsConformingFracturesInitialization; +// typedef PoromechanicsInitialization< MultiphasePoromechanicsConformingFracturesALM< CompositionalMultiphaseReservoirAndWells<> > > + // CompositionalMultiphaseReservoirPoromechanicsConformingFracturesALMInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanics<> > SinglePhasePoromechanicsInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures<> > SinglePhasePoromechanicsConformingFracturesInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsConformingFracturesInitialization; @@ -157,9 +157,9 @@ typedef PoromechanicsInitialization< HydrofractureSolver< SinglePhasePoromechani REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsConformingFracturesALMInitialization, string const &, Group * const ) -REGISTER_CATALOG_ENTRY( TaskBase, MultiphaseReservoirPoromechanicsInitialization, string const &, Group * const ) -REGISTER_CATALOG_ENTRY( TaskBase, CompositionalMultiphaseReservoirPoromechanicsConformingFracturesInitialization, string const &, Group * const ) -REGISTER_CATALOG_ENTRY( TaskBase, CompositionalMultiphaseReservoirPoromechanicsConformingFracturesALMInitialization, string const &, Group * const ) +// REGISTER_CATALOG_ENTRY( TaskBase, MultiphaseReservoirPoromechanicsInitialization, string const &, Group * const ) +// REGISTER_CATALOG_ENTRY( TaskBase, CompositionalMultiphaseReservoirPoromechanicsConformingFracturesInitialization, string const &, Group * const ) +// REGISTER_CATALOG_ENTRY( TaskBase, CompositionalMultiphaseReservoirPoromechanicsConformingFracturesALMInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhaseReservoirPoromechanicsConformingFracturesInitialization, string const &, Group * const )