Skip to content

fix: SEQ Coupling for Multiphase Poromechanics #3969

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
jhuang2601 wants to merge 2 commits into
base: develop
Choose a base branch
from
Draft

fix: SEQ Coupling for Multiphase Poromechanics #3969

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b652c03
fix seq coupling for multiphase poromechanics
jhuang2601 Feb 9, 2026
f08d67d
discard Well / test update compAmount
jafranc Feb 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 7 additions & 2 deletions src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand All @@ -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];
}
} );
}
Expand All @@ -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]
));
}
} );
}
Expand Down
99 changes: 54 additions & 45 deletions src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 )
Expand All @@ -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 )
}
Expand Down
27 changes: 27 additions & 0 deletions src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,33 @@ 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
GEOS_LOG_RANK("------ START ---------------------------------------------------");
this->flowSolver()->updateCompAmount( subRegion );
GEOS_LOG_RANK("------ END ---------------------------------------------------");
} );
} );
}
}

/**
* @brief Helper function to recompute the bulk density
* @param[in] subRegion the element subRegion
Expand Down
6 changes: 3 additions & 3 deletions ...coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
}
Expand Down
6 changes: 3 additions & 3 deletions ...eComponents/physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
}
Expand Down
16 changes: 8 additions & 8 deletions src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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 )
Expand Down
Loading