Skip to content

Commit 85f0714

Browse files
paveltominpaveltomin
authored
fix: clean up set sparsity pattern logic + matrix setup for poromechanics with contact and wells (#3834)
1 parent 28da7fb commit 85f0714

22 files changed

+234
-390
lines changed

‎src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1197,7 +1197,7 @@ void PhysicsSolverBase::setupSystem( DomainPartition & domain,
11971197
if( setSparsity )
11981198
{
11991199
SparsityPattern< globalIndex > pattern;
1200-
dofManager.setSparsityPattern( pattern );
1200+
setSparsityPattern( domain, dofManager, localMatrix, pattern );
12011201
localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) );
12021202
}
12031203
localMatrix.setName( this->getName() + "/matrix" );
@@ -1209,6 +1209,14 @@ void PhysicsSolverBase::setupSystem( DomainPartition & domain,
12091209
solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS );
12101210
}
12111211

1212+
void PhysicsSolverBase::setSparsityPattern( DomainPartition & GEOS_UNUSED_PARAM( domain ),
1213+
DofManager & dofManager,
1214+
CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),
1215+
SparsityPattern< globalIndex > & pattern )
1216+
{
1217+
dofManager.setSparsityPattern( pattern );
1218+
}
1219+
12121220
void PhysicsSolverBase::setSystemSetupTimestamp( Timestamp timestamp )
12131221
{
12141222
m_systemSetupTimestamp = timestamp;

‎src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -402,6 +402,19 @@ class PhysicsSolverBase : public ExecutableGroup
402402
ParallelVector & solution,
403403
bool const setSparsity = true );
404404

405+
/**
406+
* @brief Set the sparsity pattern of the linear system matrix
407+
* @param domain the domain containing the mesh and fields
408+
* @param dofManager degree-of-freedom manager associated with the linear system
409+
* @param localMatrix the system matrix
410+
* @param pattern the sparsity pattern to be filled
411+
*/
412+
virtual void
413+
setSparsityPattern( DomainPartition & domain,
414+
DofManager & dofManager,
415+
CRSMatrix< real64, globalIndex > & localMatrix,
416+
SparsityPattern< globalIndex > & pattern );
417+
405418
/**
406419
* @brief Create a preconditioner for this solver's linear system.
407420
* @param domain the domain containing the mesh and fields

‎src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp‎

Lines changed: 11 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -106,73 +106,44 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL
106106
setInputFlag( dataRepository::InputFlags::FALSE );
107107
}
108108

109-
/**
110-
* @brief default destructor
111-
*/
112-
virtual ~CoupledReservoirAndWellsBase () override {}
113-
114109
/**
115110
* @defgroup Solver Interface Functions
116111
*
117112
* These functions provide the primary interface that is required for derived classes
118113
*/
119114
/**@{*/
120115

121-
virtual void
122-
setupSystem( DomainPartition & domain,
123-
DofManager & dofManager,
124-
CRSMatrix< real64, globalIndex > & localMatrix,
125-
ParallelVector & rhs,
126-
ParallelVector & solution,
127-
bool const setSparsity = true ) override
116+
virtual void setSparsityPattern( DomainPartition & domain,
117+
DofManager & dofManager,
118+
CRSMatrix< real64, globalIndex > & localMatrix,
119+
SparsityPattern< globalIndex > & pattern ) override
128120
{
129-
GEOS_MARK_FUNCTION;
130-
131-
// call reservoir solver setup (needed in case of SinglePhasePoromechanicsConformingFractures)
132-
reservoirSolver()->setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );
133-
134-
dofManager.setDomain( domain );
135-
136-
Base::setupDofs( domain, dofManager );
137-
dofManager.reorderByRank();
138-
139-
// Set the sparsity pattern without reservoir-well coupling
121+
// Set the reservoir sparsity pattern without reservoir-well coupling
140122
SparsityPattern< globalIndex > patternDiag;
141-
dofManager.setSparsityPattern( patternDiag );
123+
reservoirSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternDiag );
142124

143125
// Get the original row lengths (diagonal blocks only)
144-
array1d< localIndex > rowLengths( patternDiag.numRows() );
126+
array1d< localIndex > rowLengths( patternDiag.numRows());
145127
for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow )
146128
{
147129
rowLengths[localRow] = patternDiag.numNonZeros( localRow );
148130
}
149131

150132
// Add the number of nonzeros induced by coupling on perforations
151-
addCouplingNumNonzeros( domain, dofManager, rowLengths.toView() );
133+
addCouplingNumNonzeros( domain, dofManager, rowLengths.toView());
152134

153135
// Create a new pattern with enough capacity for coupled matrix
154-
SparsityPattern< globalIndex > pattern;
155-
pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() );
136+
pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data());
156137

157138
// Copy the original nonzeros
158139
for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow )
159140
{
160141
globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous();
161-
pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) );
142+
pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ));
162143
}
163144

164145
// Add the nonzeros from coupling
165-
addCouplingSparsityPattern( domain, dofManager, pattern.toView() );
166-
167-
// Finally, steal the pattern into a CRS matrix
168-
localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) );
169-
localMatrix.setName( this->getName() + "/localMatrix" );
170-
171-
rhs.setName( this->getName() + "/rhs" );
172-
rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS );
173-
174-
solution.setName( this->getName() + "/solution" );
175-
solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS );
146+
addCouplingSparsityPattern( domain, dofManager, pattern.toView());
176147
}
177148

178149
/**@}*/

‎src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp‎

Lines changed: 14 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -264,11 +264,11 @@ real64 HydrofractureSolver< POROMECHANICS_SOLVER >::fullyCoupledSolverStep( real
264264
// Only build the sparsity pattern if the mesh has changed
265265
if( meshModificationTimestamp > getSystemSetupTimestamp() )
266266
{
267-
setupSystem( domain,
268-
m_dofManager,
269-
m_localMatrix,
270-
m_rhs,
271-
m_solution );
267+
this->setupSystem( domain,
268+
m_dofManager,
269+
m_localMatrix,
270+
m_rhs,
271+
m_solution );
272272
setSystemSetupTimestamp( meshModificationTimestamp );
273273
}
274274

@@ -486,61 +486,38 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setupCoupling( DomainPartition
486486
}
487487

488488
template< typename POROMECHANICS_SOLVER >
489-
void HydrofractureSolver< POROMECHANICS_SOLVER >::setupSystem( DomainPartition & domain,
490-
DofManager & dofManager,
491-
CRSMatrix< real64, globalIndex > & localMatrix,
492-
ParallelVector & rhs,
493-
ParallelVector & solution,
494-
bool const setSparsity )
489+
void HydrofractureSolver< POROMECHANICS_SOLVER >::setSparsityPattern( DomainPartition & domain,
490+
DofManager & dofManager,
491+
CRSMatrix< real64, globalIndex > & localMatrix,
492+
SparsityPattern< globalIndex > & pattern )
495493
{
496-
GEOS_MARK_FUNCTION;
497-
498-
GEOS_UNUSED_VAR( setSparsity );
499-
500-
dofManager.setDomain( domain );
501-
502-
setupDofs( domain, dofManager );
503-
dofManager.reorderByRank();
504-
505-
localIndex const numLocalRows = dofManager.numLocalDofs();
506-
507494
SparsityPattern< globalIndex > patternOriginal;
508495
dofManager.setSparsityPattern( patternOriginal );
509496

510497
// Get the original row lengths (diagonal blocks only)
511-
array1d< localIndex > rowLengths( patternOriginal.numRows() );
498+
array1d< localIndex > rowLengths( patternOriginal.numRows());
512499
for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
513500
{
514501
rowLengths[localRow] = patternOriginal.numNonZeros( localRow );
515502
}
516503

517504
// Add the number of nonzeros induced by coupling
518-
addFluxApertureCouplingNNZ( domain, dofManager, rowLengths.toView() );
505+
addFluxApertureCouplingNNZ( domain, dofManager, rowLengths.toView());
519506

520507
// Create a new pattern with enough capacity for coupled matrix
521-
SparsityPattern< globalIndex > pattern;
522508
pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(),
523509
patternOriginal.numColumns(),
524-
rowLengths.data() );
510+
rowLengths.data());
525511

526512
// Copy the original nonzeros
527513
for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
528514
{
529515
globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous();
530-
pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) );
516+
pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ));
531517
}
532518

533519
// Add the nonzeros from coupling
534-
addFluxApertureCouplingSparsityPattern( domain, dofManager, pattern.toView() );
535-
536-
localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) );
537-
localMatrix.setName( this->getName() + "/matrix" );
538-
539-
rhs.setName( this->getName() + "/rhs" );
540-
rhs.create( numLocalRows, MPI_COMM_GEOS );
541-
542-
solution.setName( this->getName() + "/solution" );
543-
solution.create( numLocalRows, MPI_COMM_GEOS );
520+
addFluxApertureCouplingSparsityPattern( domain, dofManager, pattern.toView());
544521

545522
setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );
546523
}

‎src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp‎

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -102,12 +102,10 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER
102102
virtual void setupCoupling( DomainPartition const & domain,
103103
DofManager & dofManager ) const override final;
104104

105-
virtual void setupSystem( DomainPartition & domain,
106-
DofManager & dofManager,
107-
CRSMatrix< real64, globalIndex > & localMatrix,
108-
ParallelVector & rhs,
109-
ParallelVector & solution,
110-
bool const setSparsity = true ) override;
105+
virtual void setSparsityPattern( DomainPartition & domain,
106+
DofManager & dofManager,
107+
CRSMatrix< real64, globalIndex > & localMatrix,
108+
SparsityPattern< globalIndex > & pattern ) override;
111109

112110
virtual void implicitStepSetup( real64 const & time_n,
113111
real64 const & dt,

‎src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp‎

Lines changed: 23 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -61,69 +61,41 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
6161
DofManager::Connector::Elem );
6262
}
6363

64-
virtual void setupSystem( DomainPartition & domain,
65-
DofManager & dofManager,
66-
CRSMatrix< real64, globalIndex > & localMatrix,
67-
ParallelVector & rhs,
68-
ParallelVector & solution,
69-
bool const setSparsity = true ) override
64+
virtual void setSparsityPattern( DomainPartition & domain,
65+
DofManager & dofManager,
66+
CRSMatrix< real64, globalIndex > & localMatrix,
67+
SparsityPattern< globalIndex > & pattern ) override
7068
{
71-
GEOS_MARK_FUNCTION;
72-
73-
GEOS_UNUSED_VAR( setSparsity );
74-
75-
/// 1. Add all coupling terms handled directly by the DofManager
76-
dofManager.setDomain( domain );
77-
this->setupDofs( domain, dofManager );
78-
dofManager.reorderByRank();
79-
80-
/// 2. Add coupling terms not added by the DofManager.
81-
localIndex const numLocalRows = dofManager.numLocalDofs();
82-
69+
// start with the flow solver sparsity pattern (it could be reservoir + wells)
8370
SparsityPattern< globalIndex > patternOriginal;
84-
dofManager.setSparsityPattern( patternOriginal );
71+
this->flowSolver()->setSparsityPattern( domain, dofManager, localMatrix, patternOriginal );
8572

8673
// Get the original row lengths (diagonal blocks only)
87-
array1d< localIndex > rowLengths( patternOriginal.numRows() );
74+
array1d< localIndex > rowLengths( patternOriginal.numRows());
8875
for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
8976
{
9077
rowLengths[localRow] = patternOriginal.numNonZeros( localRow );
9178
}
9279

9380
// Add the number of nonzeros induced by coupling
94-
addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() );
81+
addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView());
9582

9683
// Create a new pattern with enough capacity for coupled matrix
97-
SparsityPattern< globalIndex > pattern;
9884
pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(),
9985
patternOriginal.numColumns(),
100-
rowLengths.data() );
86+
rowLengths.data());
10187

10288
// Copy the original nonzeros
10389
for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow )
10490
{
10591
globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous();
106-
pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) );
92+
pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ));
10793
}
10894

10995
// Add the nonzeros from coupling
110-
addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() );
111-
112-
localMatrix.setName( this->getName() + "/matrix" );
113-
localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) );
114-
115-
rhs.setName( this->getName() + "/rhs" );
116-
rhs.create( numLocalRows, MPI_COMM_GEOS );
117-
118-
solution.setName( this->getName() + "/solution" );
119-
solution.create( numLocalRows, MPI_COMM_GEOS );
96+
addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView());
12097

12198
setUpDflux_dApertureMatrix( domain, dofManager, localMatrix );
122-
123-
// if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct )
124-
// {
125-
// m_precond = createPreconditioner( domain );
126-
// }
12799
}
128100

129101
virtual void assembleSystem( real64 const time_n,
@@ -195,6 +167,8 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
195167
{
196168
GEOS_MARK_FUNCTION;
197169

170+
integer const numComp = numFluidComponents();
171+
198172
this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, // meshBodyName,
199173
MeshLevel const & mesh,
200174
string_array const & ) // regionNames
@@ -240,7 +214,10 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
240214
if( k1 != k0 )
241215
{
242216
localIndex const numNodesPerElement = elemsToNodes[sei[iconn][k1]].size();
243-
rowLengths[rowNumber] += 3*numNodesPerElement;
217+
for( integer ic = 0; ic < numComp; ic++ )
218+
{
219+
rowLengths[rowNumber + ic] += 3*numNodesPerElement;
220+
}
244221
}
245222
}
246223
}
@@ -263,6 +240,8 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
263240
{
264241
GEOS_MARK_FUNCTION;
265242

243+
integer const numComp = numFluidComponents();
244+
266245
this->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
267246
MeshLevel const & mesh,
268247
string_array const & )
@@ -337,7 +316,10 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER,
337316
for( localIndex i = 0; i < 3; ++i )
338317
{
339318
globalIndex const colIndex = dispDofNumber[faceToNodeMap( faceIndex, a )] + LvArray::integerConversion< globalIndex >( i );
340-
pattern.insertNonZero( rowIndex, colIndex );
319+
for( integer ic = 0; ic < numComp; ic++ )
320+
{
321+
pattern.insertNonZero( rowIndex + ic, colIndex );
322+
}
341323
}
342324
}
343325
}

0 commit comments

Comments
 (0)