Skip to content

fix: clean up set sparsity pattern logic + matrix setup for poromechanics with contact and wells #3834

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

Merged
merged 17 commits into from
Sep 30, 2025
Merged

fix: clean up set sparsity pattern logic + matrix setup for poromechanics with contact and wells #3834

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
8209ece
refactor: clean up set sparsity pattern logic
Sep 25, 2025
2cb7e06
code style and doxygen
Sep 25, 2025
61dab0f
missing override
Sep 25, 2025
6840707
fix
Sep 25, 2025
9c242ec
Merge branch 'develop' into pt/set-pattern
paveltomin Sep 25, 2025
87ee761
Merge branch 'develop' into pt/set-pattern
paveltomin Sep 25, 2025
19e1d17
bug fix
Sep 26, 2025
52370d3
Merge branch 'pt/set-pattern' of https://github.com/GEOS-DEV/GEOS int…
Sep 26, 2025
000e415
Fix formatting in Group.hpp
paveltomin Sep 26, 2025
5392697
Disable SCCACHE in CI workflow
paveltomin Sep 26, 2025
8619b94
Remove USE_SCCACHE from CI tests configuration
paveltomin Sep 26, 2025
810defa
Remove unnecessary blank line in Group.hpp
paveltomin Sep 26, 2025
ec946b0
Merge branch 'develop' into pt/set-pattern
paveltomin Sep 26, 2025
137c6e9
fix: matrix setup for poromechanics with contact and wells (#3839)
paveltomin Sep 29, 2025
2001883
Merge branch 'develop' into pt/set-pattern
paveltomin Sep 29, 2025
1ae2748
Document localMatrix parameter in setSparsityPattern
paveltomin Sep 30, 2025
ee1e2ac
Merge branch 'develop' into pt/set-pattern
paveltomin Sep 30, 2025
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
Prev Previous commit
Next Next commit
fix: matrix setup for poromechanics with contact and wells (#3839)
  • Loading branch information
@paveltomin
paveltomin authored Sep 29, 2025
commit 137c6e97aa48cf1898da8ee2bdf345a6b4d29f06
3 changes: 2 additions & 1 deletion src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1197,7 +1197,7 @@ void PhysicsSolverBase::setupSystem( DomainPartition & domain,
if( setSparsity )
{
SparsityPattern< globalIndex > pattern;
setSparsityPattern( domain, dofManager, pattern );
setSparsityPattern( domain, dofManager, localMatrix, pattern );
localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) );
}
localMatrix.setName( this->getName() + "/matrix" );
Expand All @@ -1211,6 +1211,7 @@ void PhysicsSolverBase::setupSystem( DomainPartition & domain,

void PhysicsSolverBase::setSparsityPattern( DomainPartition & GEOS_UNUSED_PARAM( domain ),
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),
SparsityPattern< globalIndex > & pattern )
{
dofManager.setSparsityPattern( pattern );
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -411,6 +411,7 @@ class PhysicsSolverBase : public ExecutableGroup
virtual void
setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern );

/**
Expand Down
27 changes: 3 additions & 24 deletions src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,42 +106,21 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL
setInputFlag( dataRepository::InputFlags::FALSE );
}

/**
* @brief default destructor
*/
virtual ~CoupledReservoirAndWellsBase () override {}

/**
* @defgroup Solver Interface Functions
*
* These functions provide the primary interface that is required for derived classes
*/
/**@{*/

virtual void
setupSystem( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
ParallelVector & rhs,
ParallelVector & solution,
bool const setSparsity = true ) override
{
GEOS_MARK_FUNCTION;

// call reservoir solver setup (needed in case of SinglePhasePoromechanicsConformingFractures)
// TODO this logic does not really work - the pattern can be overwritten below
reservoirSolver()->setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );

PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );
}

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override
{
// Set the sparsity pattern without reservoir-well coupling
// Set the reservoir sparsity pattern without reservoir-well coupling
SparsityPattern< globalIndex > patternDiag;
dofManager.setSparsityPattern( patternDiag );
reservoirSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternDiag );

// Get the original row lengths (diagonal blocks only)
array1d< localIndex > rowLengths( patternDiag.numRows());
Expand Down
28 changes: 8 additions & 20 deletions src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -264,11 +264,11 @@ real64 HydrofractureSolver< POROMECHANICS_SOLVER >::fullyCoupledSolverStep( real
// Only build the sparsity pattern if the mesh has changed
if( meshModificationTimestamp > getSystemSetupTimestamp() )
{
setupSystem( domain,
m_dofManager,
m_localMatrix,
m_rhs,
m_solution );
this->setupSystem( domain,
m_dofManager,
m_localMatrix,
m_rhs,
m_solution );
setSystemSetupTimestamp( meshModificationTimestamp );
}

Expand Down Expand Up @@ -485,24 +485,10 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setupCoupling( DomainPartition
}
}

template< typename POROMECHANICS_SOLVER >
void HydrofractureSolver< POROMECHANICS_SOLVER >::setupSystem( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
ParallelVector & rhs,
ParallelVector & solution,
bool const setSparsity )
{
GEOS_MARK_FUNCTION;

PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );

setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );
}

template< typename POROMECHANICS_SOLVER >
void HydrofractureSolver< POROMECHANICS_SOLVER >::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern )
{
SparsityPattern< globalIndex > patternOriginal;
Expand Down Expand Up @@ -532,6 +518,8 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setSparsityPattern( DomainPart

// Add the nonzeros from coupling
addFluxApertureCouplingSparsityPattern( domain, dofManager, pattern.toView());

setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );
}

template< typename POROMECHANICS_SOLVER >
Expand Down
8 changes: 1 addition & 7 deletions src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,15 +102,9 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER
virtual void setupCoupling( DomainPartition const & domain,
DofManager & dofManager ) const override final;

virtual void setupSystem( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
ParallelVector & rhs,
ParallelVector & solution,
bool const setSparsity = true ) override;

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override;

virtual void implicitStepSetup( real64 const & time_n,
Expand Down
40 changes: 17 additions & 23 deletions src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,32 +61,14 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
DofManager::Connector::Elem );
}

virtual void setupSystem( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
ParallelVector & rhs,
ParallelVector & solution,
bool const setSparsity = true ) override
{
GEOS_MARK_FUNCTION;

PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );

setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );

// if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct )
// {
// m_precond = createPreconditioner( domain );
// }
}

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override
{
/// 2. Add coupling terms not added by the DofManager.
// start with the flow solver sparsity pattern (it could be reservoir + wells)
SparsityPattern< globalIndex > patternOriginal;
dofManager.setSparsityPattern( patternOriginal );
this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal );

// Get the original row lengths (diagonal blocks only)
array1d< localIndex > rowLengths( patternOriginal.numRows());
Expand All @@ -112,6 +94,8 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,

// Add the nonzeros from coupling
addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView());

setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );
}

virtual void assembleSystem( real64 const time_n,
Expand Down Expand Up @@ -183,6 +167,8 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
{
GEOS_MARK_FUNCTION;

integer const numComp = numFluidComponents();

this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, // meshBodyName,
MeshLevel const & mesh,
string_array const & ) // regionNames
Expand Down Expand Up @@ -228,7 +214,10 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
if( k1 != k0 )
{
localIndex const numNodesPerElement = elemsToNodes[sei[iconn][k1]].size();
rowLengths[rowNumber] += 3*numNodesPerElement;
for( integer ic = 0; ic < numComp; ic++ )
{
rowLengths[rowNumber + ic] += 3*numNodesPerElement;
}
}
}
}
Expand All @@ -251,6 +240,8 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
{
GEOS_MARK_FUNCTION;

integer const numComp = numFluidComponents();

this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel const & mesh,
string_array const & )
Expand Down Expand Up @@ -325,7 +316,10 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
for( localIndex i = 0; i < 3; ++i )
{
globalIndex const colIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i );
pattern.insertNonZero( rowIndex, colIndex );
for( integer ic = 0; ic < numComp; ic++ )
{
pattern.insertNonZero( rowIndex + ic, colIndex );
}
}
}
}
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,7 @@ void SinglePhasePoromechanicsEmbeddedFractures::setupCoupling( DomainPartition c

void SinglePhasePoromechanicsEmbeddedFractures::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),
SparsityPattern< globalIndex > & pattern )
{
// Set the sparsity pattern without the Kwu and Kuw blocks.
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ class SinglePhasePoromechanicsEmbeddedFractures : public SinglePhasePoromechanic

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override;

virtual void
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,7 @@ LaplaceFEM::LaplaceFEM( const string & name,
*/
void LaplaceFEM::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),
SparsityPattern< globalIndex > & pattern )
{
pattern.resize( dofManager.numLocalDofs(),
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ class LaplaceFEM : public LaplaceBaseH1

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override;

virtual void
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -989,6 +989,7 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain,

void SolidMechanicsLagrangianFEM::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),
SparsityPattern< globalIndex > & pattern )
{
pattern.resize( dofManager.numLocalDofs(),
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,7 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase
virtual void
setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override;

virtual std::unique_ptr< PreconditionerBase< LAInterface > >
Expand Down
1 change: 1 addition & 0 deletions ...onents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -217,6 +217,7 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do

void SolidMechanicsAugmentedLagrangianContact::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),
SparsityPattern< globalIndex > & pattern )
{
// Set the sparsity pattern without the Abu and Aub blocks.
Expand Down
1 change: 1 addition & 0 deletions ...onents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override final;

virtual void implicitStepSetup( real64 const & time_n,
Expand Down
3 changes: 2 additions & 1 deletion src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -230,11 +230,12 @@ void SolidMechanicsEmbeddedFractures::setupDofs( DomainPartition const & domain,

void SolidMechanicsEmbeddedFractures::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern )
{
if( m_useStaticCondensation )
{
SolidMechanicsLagrangianFEM::setSparsityPattern( domain, dofManager, pattern );
SolidMechanicsLagrangianFEM::setSparsityPattern( domain, dofManager, localMatrix, pattern );
return;
}

Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override;

virtual void implicitStepComplete( real64 const & time_n,
Expand Down
3 changes: 2 additions & 1 deletion src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -215,10 +215,11 @@ void SolidMechanicsLagrangeContact::initializePreSubGroups()

void SolidMechanicsLagrangeContact::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern )
{
// avoid calling SolidMechanicsLagrangianFEM::setSparsityPattern
PhysicsSolverBase::setSparsityPattern( domain, dofManager, pattern );
PhysicsSolverBase::setSparsityPattern( domain, dofManager, localMatrix, pattern );
}

void SolidMechanicsLagrangeContact::implicitStepSetup( real64 const & time_n,
Expand Down
1 change: 1 addition & 0 deletions src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override;

virtual std::unique_ptr< PreconditionerBase< LAInterface > >
Expand Down
1 change: 1 addition & 0 deletions ...ponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -217,6 +217,7 @@ void SolidMechanicsLagrangeContactBubbleStab::setupSystem( DomainPartition & dom

void SolidMechanicsLagrangeContactBubbleStab::setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),
SparsityPattern< globalIndex > & pattern )
{
// Set the sparsity pattern without the Abu and Aub blocks.
Expand Down
1 change: 1 addition & 0 deletions ...ponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ class SolidMechanicsLagrangeContactBubbleStab : public ContactSolverBase

virtual void setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) override final;

virtual void
Expand Down
Loading