From 8209ecebc9859f1dd2783977361348a47667bdbc Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 25 Sep 2025 12:26:37 -0500 Subject: [PATCH 01/11] refactor: clean up set sparsity pattern logic --- .../physicsSolvers/PhysicsSolverBase.cpp | 9 ++- .../physicsSolvers/PhysicsSolverBase.hpp | 5 ++ .../CoupledReservoirAndWellsBase.hpp | 34 +++----- .../multiphysics/HydrofractureSolver.cpp | 39 ++++----- .../multiphysics/HydrofractureSolver.hpp | 4 + .../PoromechanicsConformingFractures.hpp | 48 +++++------ ...glePhasePoromechanicsEmbeddedFractures.cpp | 42 +++------- ...glePhasePoromechanicsEmbeddedFractures.hpp | 9 +-- .../physicsSolvers/simplePDE/LaplaceFEM.cpp | 44 ++++------ .../physicsSolvers/simplePDE/LaplaceFEM.hpp | 14 +--- .../SolidMechanicsLagrangianFEM.cpp | 29 ++++--- .../SolidMechanicsLagrangianFEM.hpp | 5 ++ ...lidMechanicsAugmentedLagrangianContact.cpp | 50 +++++------- ...lidMechanicsAugmentedLagrangianContact.hpp | 4 + .../SolidMechanicsEmbeddedFractures.cpp | 80 +++++++------------ .../SolidMechanicsEmbeddedFractures.hpp | 9 +-- .../contact/SolidMechanicsLagrangeContact.cpp | 1 - ...olidMechanicsLagrangeContactBubbleStab.cpp | 41 ++++------ ...olidMechanicsLagrangeContactBubbleStab.hpp | 4 + .../contact/SolidMechanicsPenaltyContact.cpp | 31 +++---- .../contact/SolidMechanicsPenaltyContact.hpp | 12 +-- 21 files changed, 202 insertions(+), 312 deletions(-) diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 7c3da3cc8cb..718f935cb6a 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -1197,7 +1197,7 @@ void PhysicsSolverBase::setupSystem( DomainPartition & domain, if( setSparsity ) { SparsityPattern< globalIndex > pattern; - dofManager.setSparsityPattern( pattern ); + setSparsityPattern( domain, dofManager, pattern ); localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); } localMatrix.setName( this->getName() + "/matrix" ); @@ -1209,6 +1209,13 @@ void PhysicsSolverBase::setupSystem( DomainPartition & domain, solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); } +void PhysicsSolverBase::setSparsityPattern( DomainPartition & GEOS_UNUSED_PARAM( domain ), + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) +{ + dofManager.setSparsityPattern( pattern ); +} + void PhysicsSolverBase::setSystemSetupTimestamp( Timestamp timestamp ) { m_systemSetupTimestamp = timestamp; diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index 0ee65309d95..9b2fc51ebb6 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -402,6 +402,11 @@ class PhysicsSolverBase : public ExecutableGroup ParallelVector & solution, bool const setSparsity = true ); + virtual void + setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ); + /** * @brief Create a preconditioner for this solver's linear system. * @param domain the domain containing the mesh and fields diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index 1370432d37c..89a6ae679fc 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -129,50 +129,42 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL 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 ); - dofManager.setDomain( domain ); - - Base::setupDofs( domain, dofManager ); - dofManager.reorderByRank(); + PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + } + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override + { // Set the sparsity pattern without reservoir-well coupling SparsityPattern< globalIndex > patternDiag; dofManager.setSparsityPattern( patternDiag ); // Get the original row lengths (diagonal blocks only) - array1d< localIndex > rowLengths( patternDiag.numRows() ); + array1d< localIndex > rowLengths( patternDiag.numRows()); for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { rowLengths[localRow] = patternDiag.numNonZeros( localRow ); } // Add the number of nonzeros induced by coupling on perforations - addCouplingNumNonzeros( domain, dofManager, rowLengths.toView() ); + addCouplingNumNonzeros( domain, dofManager, rowLengths.toView()); // Create a new pattern with enough capacity for coupled matrix - SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() ); + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data()); // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) ); + globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } // Add the nonzeros from coupling - addCouplingSparsityPattern( domain, dofManager, pattern.toView() ); - - // Finally, steal the pattern into a CRS matrix - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); - localMatrix.setName( this->getName() + "/localMatrix" ); - - rhs.setName( this->getName() + "/rhs" ); - rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - - solution.setName( this->getName() + "/solution" ); - solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + addCouplingSparsityPattern( domain, dofManager, pattern.toView()); } /**@}*/ diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index f29995ba60a..2eb5d02ba55 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -495,54 +495,43 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setupSystem( DomainPartition & { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( setSparsity ); + PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); - dofManager.setDomain( domain ); - - setupDofs( domain, dofManager ); - dofManager.reorderByRank(); - - localIndex const numLocalRows = dofManager.numLocalDofs(); + setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); +} +template< typename POROMECHANICS_SOLVER > +void HydrofractureSolver< POROMECHANICS_SOLVER >::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) +{ SparsityPattern< globalIndex > patternOriginal; dofManager.setSparsityPattern( patternOriginal ); // Get the original row lengths (diagonal blocks only) - array1d< localIndex > rowLengths( patternOriginal.numRows() ); + array1d< localIndex > rowLengths( patternOriginal.numRows()); for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) { rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); } // Add the number of nonzeros induced by coupling - addFluxApertureCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addFluxApertureCouplingNNZ( domain, dofManager, rowLengths.toView()); // Create a new pattern with enough capacity for coupled matrix - SparsityPattern< globalIndex > pattern; pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), patternOriginal.numColumns(), - rowLengths.data() ); + rowLengths.data()); // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) { - globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + globalIndex const *cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow )); } // Add the nonzeros from coupling - addFluxApertureCouplingSparsityPattern( domain, dofManager, pattern.toView() ); - - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); - localMatrix.setName( this->getName() + "/matrix" ); - - rhs.setName( this->getName() + "/rhs" ); - rhs.create( numLocalRows, MPI_COMM_GEOS ); - - solution.setName( this->getName() + "/solution" ); - solution.create( numLocalRows, MPI_COMM_GEOS ); - - setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); + addFluxApertureCouplingSparsityPattern( domain, dofManager, pattern.toView()); } template< typename POROMECHANICS_SOLVER > diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp index a8f431e8d7a..9d9e5da29c9 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp @@ -109,6 +109,10 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER ParallelVector & solution, bool const setSparsity = true ) override; + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override; + virtual void implicitStepSetup( real64 const & time_n, real64 const & dt, DomainPartition & domain ) override final; diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index 7db990f9c29..7c2d407bc79 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -70,60 +70,48 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( setSparsity ); + PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); - /// 1. Add all coupling terms handled directly by the DofManager - dofManager.setDomain( domain ); - this->setupDofs( domain, dofManager ); - dofManager.reorderByRank(); + setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); - /// 2. Add coupling terms not added by the DofManager. - localIndex const numLocalRows = dofManager.numLocalDofs(); + // if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + // { + // m_precond = createPreconditioner( domain ); + // } + } + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override + { + /// 2. Add coupling terms not added by the DofManager. SparsityPattern< globalIndex > patternOriginal; dofManager.setSparsityPattern( patternOriginal ); // Get the original row lengths (diagonal blocks only) - array1d< localIndex > rowLengths( patternOriginal.numRows() ); + array1d< localIndex > rowLengths( patternOriginal.numRows()); for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) { rowLengths[localRow] = patternOriginal.numNonZeros( localRow ); } // Add the number of nonzeros induced by coupling - addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView() ); + addTransmissibilityCouplingNNZ( domain, dofManager, rowLengths.toView()); // Create a new pattern with enough capacity for coupled matrix - SparsityPattern< globalIndex > pattern; pattern.resizeFromRowCapacities< parallelHostPolicy >( patternOriginal.numRows(), patternOriginal.numColumns(), - rowLengths.data() ); + rowLengths.data()); // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) { - globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow ) ); + globalIndex const *cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow )); } // Add the nonzeros from coupling - addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView() ); - - localMatrix.setName( this->getName() + "/matrix" ); - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); - - rhs.setName( this->getName() + "/rhs" ); - rhs.create( numLocalRows, MPI_COMM_GEOS ); - - solution.setName( this->getName() + "/solution" ); - solution.create( numLocalRows, MPI_COMM_GEOS ); - - setUpDflux_dApertureMatrix( domain, dofManager, localMatrix ); - - // if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) - // { - // m_precond = createPreconditioner( domain ); - // } + addTransmissibilityCouplingPattern( domain, dofManager, pattern.toView()); } virtual void assembleSystem( real64 const time_n, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index 146fa3ad0a5..35d9f5a196e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -121,60 +121,36 @@ void SinglePhasePoromechanicsEmbeddedFractures::setupCoupling( DomainPartition c meshTargets ); } -void SinglePhasePoromechanicsEmbeddedFractures::setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity ) +void SinglePhasePoromechanicsEmbeddedFractures::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) { - // Add missing couplings ( matrix pressure with displacement jump and jump - displacement ) - - GEOS_MARK_FUNCTION; - - GEOS_UNUSED_VAR( setSparsity ); - - dofManager.setDomain( domain ); - setupDofs( domain, dofManager ); - dofManager.reorderByRank(); - // Set the sparsity pattern without the Kwu and Kuw blocks. SparsityPattern< globalIndex > patternDiag; dofManager.setSparsityPattern( patternDiag ); // Get the original row lengths (diagonal blocks only) - array1d< localIndex > rowLengths( patternDiag.numRows() ); + array1d< localIndex > rowLengths( patternDiag.numRows()); for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { rowLengths[localRow] = patternDiag.numNonZeros( localRow ); } // Add the number of nonzeros induced by coupling jump-pm - addCouplingNumNonzeros( domain, dofManager, rowLengths.toView() ); + addCouplingNumNonzeros( domain, dofManager, rowLengths.toView()); // Create a new pattern with enough capacity for coupled matrix - SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() ); + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data()); // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) ); + globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } // Add the nonzeros from coupling - addCouplingSparsityPattern( domain, dofManager, pattern.toView() ); - - // Finally, steal the pattern into a CRS matrix - localMatrix.setName( this->getName() + "/localMatrix" ); - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); - - rhs.setName( this->getName() + "/rhs" ); - rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - - solution.setName( this->getName() + "/solution" ); - solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + addCouplingSparsityPattern( domain, dofManager, pattern.toView()); } void SinglePhasePoromechanicsEmbeddedFractures::addCouplingNumNonzeros( DomainPartition & domain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp index 4c2732cd83d..ecfbbaaf7ba 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -52,12 +52,9 @@ class SinglePhasePoromechanicsEmbeddedFractures : public SinglePhasePoromechanic virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) 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, + SparsityPattern< globalIndex > & pattern ); virtual void setupCoupling( DomainPartition const & domain, diff --git a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp index 08e01c68ec2..28d1d687c96 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp +++ b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp @@ -83,36 +83,25 @@ LaplaceFEM::LaplaceFEM( const string & name, {} //END_SPHINX_INCLUDE_CONSTRUCTOR -LaplaceFEM::~LaplaceFEM() -{ - // TODO Auto-generated destructor stub -} - -/* SETUP SYSTEM - Setting up the system using the base class method +/* SETUP MATRIX PATTERN + This method sets up the sparsity pattern of the matrix that will be used to solve the + Laplace equation. */ -void LaplaceFEM::setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity ) +void LaplaceFEM::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) { - GEOS_MARK_FUNCTION; - PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + pattern.resize( dofManager.numLocalDofs(), + dofManager.numGlobalDofs(), + 8 * 8 * 3 ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) { NodeManager const & nodeManager = mesh.getNodeManager(); - string const dofKey = dofManager.getKey( m_fieldName ); arrayView1d< globalIndex const > const & - dofIndex = nodeManager.getReference< globalIndex_array >( dofKey ); - - SparsityPattern< globalIndex > sparsityPattern( dofManager.numLocalDofs(), - dofManager.numGlobalDofs(), - 8*8*3 ); + dofIndex = nodeManager.getReference< globalIndex_array >( dofManager.getKey( m_fieldName ) ); finiteElement::fillSparsity< CellElementSubRegion, LaplaceFEMKernel >( mesh, @@ -120,13 +109,12 @@ void LaplaceFEM::setupSystem( DomainPartition & domain, this->getDiscretizationName(), dofIndex, dofManager.rankOffset(), - sparsityPattern ); + pattern ); - sparsityPattern.compress(); - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( sparsityPattern ) ); } ); -} + pattern.compress(); +} /* ASSEMBLE SYSTEM diff --git a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp index 802d6cb1bca..5ca5d971a98 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp @@ -38,9 +38,6 @@ class LaplaceFEM : public LaplaceBaseH1 LaplaceFEM( const string & name, Group * const parent ); - /// Destructor - virtual ~LaplaceFEM() override; - /// "CatalogName()" return the string used as XML tag in the input file. It ties the XML tag with /// this C++ classes. This is important. static string catalogName() { return "LaplaceFEM"; } @@ -58,13 +55,10 @@ class LaplaceFEM : public LaplaceBaseH1 // /**@{*/ //START_SPHINX_INCLUDE_SOLVERINTERFACE - virtual void - setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity = false ) override; + + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override; virtual void assembleSystem( real64 const time, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 92fff235619..c2ee4de1214 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -978,11 +978,22 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain, bool const setSparsity ) { GEOS_MARK_FUNCTION; + PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); - SparsityPattern< globalIndex > sparsityPattern( dofManager.numLocalDofs(), - dofManager.numGlobalDofs(), - 8*8*3*1.2 ); + if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) + { + m_precond = createPreconditioner( domain ); + } +} + +void SolidMechanicsLagrangianFEM::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) +{ + pattern.resize( dofManager.numLocalDofs(), + dofManager.numGlobalDofs(), + 8*8*3*1.2 ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, @@ -996,21 +1007,15 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain, fillSparsity< CellElementSubRegion, solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, regionNames, - this->getDiscretizationName(), + getDiscretizationName(), dofNumber, dofManager.rankOffset(), - sparsityPattern ); + pattern ); } ); - sparsityPattern.compress(); - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( sparsityPattern ) ); - - if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) - { - m_precond = createPreconditioner( domain ); - } + pattern.compress(); } void SolidMechanicsLagrangianFEM::computeRigidBodyModes( DomainPartition & domain ) const diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index 89579b79d15..5872ca4b491 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -113,6 +113,11 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase ParallelVector & solution, bool const setSparsity = false ) override; + virtual void + setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override; + virtual std::unique_ptr< PreconditionerBase< LAInterface > > createPreconditioner( DomainPartition & domain ) const override; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 3076a3d0776..d760be84278 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -114,11 +114,11 @@ void SolidMechanicsAugmentedLagrangianContact::registerDataOnMesh( dataRepositor FaceManager & faceManager = meshLevel.getFaceManager(); // Register the total bubble displacement - faceManager.registerField< contact::totalBubbleDisplacement >( this->getName() ). + faceManager.registerField< contact::totalBubbleDisplacement >( getName() ). reference().resizeDimension< 1 >( 3 ); // Register the incremental bubble displacement - faceManager.registerField< contact::incrementalBubbleDisplacement >( this->getName() ). + faceManager.registerField< contact::incrementalBubbleDisplacement >( getName() ). reference().resizeDimension< 1 >( 3 ); } ); @@ -202,7 +202,6 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( setSparsity ); // Create the lists of interface elements that have same type. createFaceTypeList( domain ); @@ -213,48 +212,39 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do // Create the list of cell elements that they are enriched with bubble functions. createBubbleCellList( domain ); - dofManager.setDomain( domain ); - setupDofs( domain, dofManager ); - dofManager.reorderByRank(); + PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); +} +void SolidMechanicsAugmentedLagrangianContact::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) +{ // Set the sparsity pattern without the Abu and Aub blocks. SparsityPattern< globalIndex > patternDiag; dofManager.setSparsityPattern( patternDiag ); // Get the original row lengths (diagonal blocks only) - array1d< localIndex > rowLengths( patternDiag.numRows() ); + array1d< localIndex > rowLengths( patternDiag.numRows()); for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { rowLengths[localRow] = patternDiag.numNonZeros( localRow ); } // Add the number of nonzeros induced by coupling - this->addCouplingNumNonzeros( domain, dofManager, rowLengths.toView() ); + addCouplingNumNonzeros( domain, dofManager, rowLengths.toView()); // Create a new pattern with enough capacity for coupled matrix - SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() ); + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data()); // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) ); + globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } // Add the nonzeros from coupling - this->addCouplingSparsityPattern( domain, dofManager, pattern.toView() ); - - // Finally, steal the pattern into a CRS matrix - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); - localMatrix.setName( this->getName() + "/localMatrix" ); - - rhs.setName( this->getName() + "/rhs" ); - rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - - solution.setName( this->getName() + "/solution" ); - solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - + addCouplingSparsityPattern( domain, dofManager, pattern.toView()); } void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const & time_n, @@ -400,7 +390,7 @@ void SolidMechanicsAugmentedLagrangianContact::assembleSystem( real64 const time arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - string const & fractureRegionName = this->getUniqueFractureRegionName(); + string const & fractureRegionName = getUniqueFractureRegionName(); forFiniteElementOnStickFractureSubRegions( meshName, [&] ( string const &, finiteElement::FiniteElementBase const & subRegionFE, @@ -748,7 +738,7 @@ void SolidMechanicsAugmentedLagrangianContact::applySystemSolution( DofManager c arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - string const & fractureRegionName = this->getUniqueFractureRegionName(); + string const & fractureRegionName = getUniqueFractureRegionName(); CRSMatrix< real64, globalIndex > const voidMatrix; array1d< real64 > const voidRhs; @@ -1146,8 +1136,8 @@ void SolidMechanicsAugmentedLagrangianContact::updateStickSlipList( DomainPartit slipList_v[kfe] = vals_v[nStick+kfe]; } ); - this->m_faceTypesToFaceElementsStick[meshName][finiteElementName] = stickList; - this->m_faceTypesToFaceElementsSlip[meshName][finiteElementName] = slipList; + m_faceTypesToFaceElementsStick[meshName][finiteElementName] = stickList; + m_faceTypesToFaceElementsSlip[meshName][finiteElementName] = slipList; GEOS_LOG_LEVEL_RANK_0( logInfo::ConfigurationStatistics, GEOS_FMT( "# stick elements: {}, # slip elements: {}", nStick, nSlip )) } ); @@ -1230,8 +1220,8 @@ void SolidMechanicsAugmentedLagrangianContact::createFaceTypeList( DomainPartiti quadList_v[kfe] = vals_v[nTri+kfe]; } ); - this->m_faceTypesToFaceElements[meshName]["Quadrilateral"] = quadList; - this->m_faceTypesToFaceElements[meshName]["Triangle"] = triList; + m_faceTypesToFaceElements[meshName]["Quadrilateral"] = quadList; + m_faceTypesToFaceElements[meshName]["Triangle"] = triList; } ); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 65e934c1146..c9be177fbc3 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -59,6 +59,10 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase ParallelVector & solution, bool const setSparsity = true ) override final; + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override final; + virtual void implicitStepSetup( real64 const & time_n, real64 const & dt, DomainPartition & domain ) override final; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp index 87ddc22cd5f..ee1478aae2b 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp @@ -114,7 +114,7 @@ void SolidMechanicsEmbeddedFractures::registerDataOnMesh( dataRepository::Group void SolidMechanicsEmbeddedFractures::initializePostInitialConditionsPreSubGroups() { ContactSolverBase::initializePostInitialConditionsPreSubGroups(); - updateState( this->getGroupByPath< DomainPartition >( "/Problem/domain" ) ); + updateState( getGroupByPath< DomainPartition >( "/Problem/domain" ) ); } void SolidMechanicsEmbeddedFractures::resetStateToBeginningOfStep( DomainPartition & domain ) @@ -233,66 +233,42 @@ void SolidMechanicsEmbeddedFractures::setupDofs( DomainPartition const & domain, } } -void SolidMechanicsEmbeddedFractures::setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity ) +void SolidMechanicsEmbeddedFractures::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) { - GEOS_MARK_FUNCTION; - - if( !m_useStaticCondensation ) + if( m_useStaticCondensation ) { + SolidMechanicsLagrangianFEM::setSparsityPattern( domain, dofManager, pattern ); + return; + } - GEOS_UNUSED_VAR( setSparsity ); - - dofManager.setDomain( domain ); - setupDofs( domain, dofManager ); - dofManager.reorderByRank(); - - // Set the sparsity pattern without the Kwu and Kuw blocks. - SparsityPattern< globalIndex > patternDiag; - dofManager.setSparsityPattern( patternDiag ); - - // Get the original row lengths (diagonal blocks only) - array1d< localIndex > rowLengths( patternDiag.numRows() ); - for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) - { - rowLengths[localRow] = patternDiag.numNonZeros( localRow ); - } - - // Add the number of nonzeros induced by coupling - addCouplingNumNonzeros( domain, dofManager, rowLengths.toView() ); - - // Create a new pattern with enough capacity for coupled matrix - SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() ); - - // Copy the original nonzeros - for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) - { - globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) ); - } + // Set the sparsity pattern without the Kwu and Kuw blocks. + SparsityPattern< globalIndex > patternDiag; + dofManager.setSparsityPattern( patternDiag ); - // Add the nonzeros from coupling - addCouplingSparsityPattern( domain, dofManager, pattern.toView() ); + // Get the original row lengths (diagonal blocks only) + array1d< localIndex > rowLengths( patternDiag.numRows() ); + for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) + { + rowLengths[localRow] = patternDiag.numNonZeros( localRow ); + } - // Finally, steal the pattern into a CRS matrix - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); - localMatrix.setName( this->getName() + "/localMatrix" ); + // Add the number of nonzeros induced by coupling + addCouplingNumNonzeros( domain, dofManager, rowLengths.toView() ); - rhs.setName( this->getName() + "/rhs" ); - rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + // Create a new pattern with enough capacity for coupled matrix + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() ); - solution.setName( this->getName() + "/solution" ); - solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - } - else + // Copy the original nonzeros + for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - SolidMechanicsLagrangianFEM::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) ); } + + // Add the nonzeros from coupling + addCouplingSparsityPattern( domain, dofManager, pattern.toView() ); } void SolidMechanicsEmbeddedFractures::assembleSystem( real64 const time, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp index 5eee0928256..b4834f474e7 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp @@ -52,12 +52,9 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase virtual void setupDofs( DomainPartition const & domain, DofManager & dofManager ) const override; - 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, + SparsityPattern< globalIndex > & pattern ) override; virtual void implicitStepComplete( real64 const & time_n, real64 const & dt, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index fe70b6d2160..ba38fddf34d 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -232,7 +232,6 @@ void SolidMechanicsLagrangeContact::setupSystem( DomainPartition & domain, { createPreconditioner( domain ); } - } void SolidMechanicsLagrangeContact::implicitStepSetup( real64 const & time_n, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp index a00980e3979..8dcbc236484 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp @@ -197,10 +197,8 @@ void SolidMechanicsLagrangeContactBubbleStab::setupSystem( DomainPartition & dom CRSMatrix< real64, globalIndex > & localMatrix, ParallelVector & rhs, ParallelVector & solution, - bool const GEOS_UNUSED_PARAM( setSparsity ) ) + bool const setSparsity ) { - - // setup monolithic coupled system // Create the lists of interface elements that have same type. @@ -212,50 +210,41 @@ void SolidMechanicsLagrangeContactBubbleStab::setupSystem( DomainPartition & dom // Create the list of cell elements that they are enriched with bubble functions. createBubbleCellList( domain ); - dofManager.setDomain( domain ); - setupDofs( domain, dofManager ); - dofManager.reorderByRank(); + PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + + computeRotationMatrices( domain ); +} +void SolidMechanicsLagrangeContactBubbleStab::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) +{ // Set the sparsity pattern without the Abu and Aub blocks. SparsityPattern< globalIndex > patternDiag; dofManager.setSparsityPattern( patternDiag ); // Get the original row lengths (diagonal blocks only) - array1d< localIndex > rowLengths( patternDiag.numRows() ); + array1d< localIndex > rowLengths( patternDiag.numRows()); for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { rowLengths[localRow] = patternDiag.numNonZeros( localRow ); } // Add the number of nonzeros induced by coupling - this->addCouplingNumNonzeros( domain, dofManager, rowLengths.toView() ); + this->addCouplingNumNonzeros( domain, dofManager, rowLengths.toView()); // Create a new pattern with enough capacity for coupled matrix - SparsityPattern< globalIndex > pattern; - pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() ); + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data()); // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); - pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) ); + globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } // Add the nonzeros from coupling - this->addCouplingSparsityPattern( domain, dofManager, pattern.toView() ); - - // Finally, steal the pattern into a CRS matrix - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); - localMatrix.setName( this->getName() + "/localMatrix" ); - - rhs.setName( this->getName() + "/rhs" ); - rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - - solution.setName( this->getName() + "/solution" ); - solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); - - computeRotationMatrices( domain ); - + this->addCouplingSparsityPattern( domain, dofManager, pattern.toView()); } void SolidMechanicsLagrangeContactBubbleStab::computeRotationMatrices( DomainPartition & domain ) const diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp index bafdec52fac..c7ef79a835c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp @@ -69,6 +69,10 @@ class SolidMechanicsLagrangeContactBubbleStab : public ContactSolverBase ParallelVector & solution, bool const setSparsity = true ) override final; + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override final; + virtual void implicitStepSetup( real64 const & time_n, real64 const & dt, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.cpp index 852dbe48aca..b7cfd31b0c5 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.cpp @@ -44,24 +44,13 @@ SolidMechanicsPenaltyContact::SolidMechanicsPenaltyContact( const string & name, setDescription( "Value of the penetration penalty stiffness. Units of Pressure/length" ); } -SolidMechanicsPenaltyContact::~SolidMechanicsPenaltyContact() +void SolidMechanicsPenaltyContact::setSparsityPattern( DofManager & dofManager, + DomainPartition & domain, + SparsityPattern< globalIndex > & pattern ) { - // TODO Auto-generated destructor stub -} - -void SolidMechanicsPenaltyContact::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, false ); - - SparsityPattern< globalIndex > sparsityPattern( dofManager.numLocalDofs(), - dofManager.numGlobalDofs(), - 8*8*3*1.2 ); + pattern.resize( dofManager.numLocalDofs(), + dofManager.numGlobalDofs(), + 8*8*3*1.2 ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, @@ -71,7 +60,6 @@ void SolidMechanicsPenaltyContact::setupSystem( DomainPartition & domain, arrayView1d< globalIndex const > const dofNumber = nodeManager.getReference< globalIndex_array >( dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); - ElementRegionManager const & elemManager = mesh.getElemManager(); string_array allFaceElementRegions; elemManager.forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion const & elemRegion ) @@ -86,7 +74,7 @@ void SolidMechanicsPenaltyContact::setupSystem( DomainPartition & domain, this->getDiscretizationName(), dofNumber, dofManager.rankOffset(), - sparsityPattern ); + pattern ); finiteElement::fillSparsity< CellElementSubRegion, solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, @@ -94,13 +82,12 @@ void SolidMechanicsPenaltyContact::setupSystem( DomainPartition & domain, this->getDiscretizationName(), dofNumber, dofManager.rankOffset(), - sparsityPattern ); + pattern ); } ); - sparsityPattern.compress(); - localMatrix.assimilate< parallelDevicePolicy<> >( std::move( sparsityPattern ) ); + pattern.compress(); } void SolidMechanicsPenaltyContact::assembleSystem( real64 const time, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.hpp index d6c0e9ffc3e..0dff34a4b7a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsPenaltyContact.hpp @@ -37,8 +37,6 @@ class SolidMechanicsPenaltyContact : public ContactSolverBase SolidMechanicsPenaltyContact( const string & name, Group * const parent ); - ~SolidMechanicsPenaltyContact() override; - /** * @brief name of the node manager in the object catalog * @return string that contains the catalog name to generate a new NodeManager object through the object catalog. @@ -55,13 +53,9 @@ class SolidMechanicsPenaltyContact : public ContactSolverBase /// String used to form the solverName used to register single-physics solvers in CoupledSolver static string coupledSolverAttributePrefix() { return "PenaltyContact"; } - virtual void - setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity = true ) override final; + virtual void setSparsityPattern( DofManager & dofManager, + DomainPartition & domain, + SparsityPattern< globalIndex > & pattern ) override final; virtual void assembleSystem( real64 const time, From 2cb7e065a683c714197c3f816abf5b3a00dbf268 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 25 Sep 2025 12:32:56 -0500 Subject: [PATCH 02/11] code style and doxygen --- src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp | 6 ++++++ .../multiphysics/CoupledReservoirAndWellsBase.hpp | 2 +- .../physicsSolvers/multiphysics/HydrofractureSolver.cpp | 2 +- .../multiphysics/PoromechanicsConformingFractures.hpp | 2 +- .../SinglePhasePoromechanicsEmbeddedFractures.cpp | 2 +- .../contact/SolidMechanicsAugmentedLagrangianContact.cpp | 2 +- .../contact/SolidMechanicsEmbeddedFractures.hpp | 4 ++-- .../contact/SolidMechanicsLagrangeContactBubbleStab.cpp | 2 +- 8 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index 9b2fc51ebb6..e3116f65540 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -402,6 +402,12 @@ class PhysicsSolverBase : public ExecutableGroup ParallelVector & solution, bool const setSparsity = true ); + /** + * @brief Set the sparsity pattern of the linear system matrix + * @param domain the domain containing the mesh and fields + * @param dofManager degree-of-freedom manager associated with the linear system + * @param pattern the sparsity pattern to be filled + */ virtual void setSparsityPattern( DomainPartition & domain, DofManager & dofManager, diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index 89a6ae679fc..7af00b2d6d5 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -159,7 +159,7 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index 2eb5d02ba55..06d5f7a816e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -526,7 +526,7 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::setSparsityPattern( DomainPart // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) { - globalIndex const *cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow )); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index 7c2d407bc79..4a2f96f6bf6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -106,7 +106,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternOriginal.numRows(); ++localRow ) { - globalIndex const *cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); + globalIndex const * cols = patternOriginal.getColumns( localRow ).dataIfContiguous(); pattern.insertNonZeros( localRow, cols, cols + patternOriginal.numNonZeros( localRow )); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index 35d9f5a196e..c41067a82d2 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -145,7 +145,7 @@ void SinglePhasePoromechanicsEmbeddedFractures::setSparsityPattern( DomainPartit // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index d760be84278..2aa85c8b4f5 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -239,7 +239,7 @@ void SolidMechanicsAugmentedLagrangianContact::setSparsityPattern( DomainPartiti // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp index b4834f474e7..b2fec4d7b8f 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp @@ -53,8 +53,8 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase DofManager & dofManager ) const override; virtual void setSparsityPattern( DomainPartition & domain, - DofManager & dofManager, - SparsityPattern< globalIndex > & pattern ) override; + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override; virtual void implicitStepComplete( real64 const & time_n, real64 const & dt, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp index 8dcbc236484..bd980d76519 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp @@ -239,7 +239,7 @@ void SolidMechanicsLagrangeContactBubbleStab::setSparsityPattern( DomainPartitio // Copy the original nonzeros for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) { - globalIndex const *cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow )); } From 61dab0f16103a762e5e2fec120dd5a581fcddb9f Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 25 Sep 2025 12:54:40 -0500 Subject: [PATCH 03/11] missing override --- .../multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp index ecfbbaaf7ba..23ffc6b6c74 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -54,7 +54,7 @@ class SinglePhasePoromechanicsEmbeddedFractures : public SinglePhasePoromechanic virtual void setSparsityPattern( DomainPartition & domain, DofManager & dofManager, - SparsityPattern< globalIndex > & pattern ); + SparsityPattern< globalIndex > & pattern ) override; virtual void setupCoupling( DomainPartition const & domain, From 6840707b542c22e4ae4c22424730e6f2839e9359 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 25 Sep 2025 14:29:52 -0500 Subject: [PATCH 04/11] fix --- .../solidMechanics/SolidMechanicsLagrangianFEM.hpp | 2 +- .../contact/SolidMechanicsEmbeddedFractures.cpp | 5 ----- .../contact/SolidMechanicsEmbeddedFractures.hpp | 2 -- 3 files changed, 1 insertion(+), 8 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index 5872ca4b491..b0b2723ef63 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -111,7 +111,7 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase CRSMatrix< real64, globalIndex > & localMatrix, ParallelVector & rhs, ParallelVector & solution, - bool const setSparsity = false ) override; + bool setSparsity = true ) override; virtual void setSparsityPattern( DomainPartition & domain, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp index ee1478aae2b..4bdca0418d5 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp @@ -59,11 +59,6 @@ SolidMechanicsEmbeddedFractures::SolidMechanicsEmbeddedFractures( const string & setDescription( "Value of the penetration penalty stiffness. Units of Pressure/length" ); } -SolidMechanicsEmbeddedFractures::~SolidMechanicsEmbeddedFractures() -{ - // TODO Auto-generated destructor stub -} - void SolidMechanicsEmbeddedFractures::postInputInitialization() { ContactSolverBase::postInputInitialization(); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp index b2fec4d7b8f..fbaff010cf0 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp @@ -32,8 +32,6 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase SolidMechanicsEmbeddedFractures( const string & name, Group * const parent ); - ~SolidMechanicsEmbeddedFractures() override; - /** * @brief name of the node manager in the object catalog * @return string that contains the catalog name to generate a new NodeManager object through the object catalog. From 19e1d17bbc8b0ee296e3069cd7d40cde324d143d Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Thu, 25 Sep 2025 19:03:53 -0500 Subject: [PATCH 05/11] bug fix --- .../contact/SolidMechanicsLagrangeContact.cpp | 23 ++++--------------- .../contact/SolidMechanicsLagrangeContact.hpp | 10 +++----- 2 files changed, 8 insertions(+), 25 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index ba38fddf34d..f13cefe4f87 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -213,25 +213,12 @@ void SolidMechanicsLagrangeContact::initializePreSubGroups() } -void SolidMechanicsLagrangeContact::setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const GEOS_UNUSED_PARAM( setSparsity ) ) +void SolidMechanicsLagrangeContact::setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) { - if( m_precond ) - { - m_precond->clear(); - } - - // setup monolithic coupled system - PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, true ); // "true" is to force setSparsity - - if( !m_precond && m_linearSolverParameters.get().solverType != LinearSolverParameters::SolverType::direct ) - { - createPreconditioner( domain ); - } + // avoid calling SolidMechanicsLagrangianFEM::setSparsityPattern + PhysicsSolverBase::setSparsityPattern( domain, dofManager, pattern ); } void SolidMechanicsLagrangeContact::implicitStepSetup( real64 const & time_n, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp index 5fd6832deb1..b0089d86fcc 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp @@ -57,13 +57,9 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase setupDofs( DomainPartition const & domain, DofManager & dofManager ) const override; - virtual void - setupSystem( DomainPartition & domain, - DofManager & dofManager, - CRSMatrix< real64, globalIndex > & localMatrix, - ParallelVector & rhs, - ParallelVector & solution, - bool const setSparsity = true ) override final; + virtual void setSparsityPattern( DomainPartition & domain, + DofManager & dofManager, + SparsityPattern< globalIndex > & pattern ) override; virtual std::unique_ptr< PreconditionerBase< LAInterface > > createPreconditioner( DomainPartition & domain ) const override; From 000e41596334aac8e826af8e7860e65dc301c8e6 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 26 Sep 2025 09:50:11 -0500 Subject: [PATCH 06/11] Fix formatting in Group.hpp --- src/coreComponents/dataRepository/Group.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/coreComponents/dataRepository/Group.hpp b/src/coreComponents/dataRepository/Group.hpp index 0aef2b787b8..38659902b44 100644 --- a/src/coreComponents/dataRepository/Group.hpp +++ b/src/coreComponents/dataRepository/Group.hpp @@ -1752,6 +1752,7 @@ Group::addLogLevel() } } /* end namespace dataRepository */ + } /* end namespace geos */ #endif /* GEOS_DATAREPOSITORY_GROUP_HPP_ */ From 53926972b9f9d02b047c315a6a054c3a590c5da0 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 26 Sep 2025 10:29:27 -0500 Subject: [PATCH 07/11] Disable SCCACHE in CI workflow --- .github/workflows/ci_tests.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml index e1e7506540f..e8ab48b2596 100644 --- a/.github/workflows/ci_tests.yml +++ b/.github/workflows/ci_tests.yml @@ -236,6 +236,7 @@ jobs: GCP_BUCKET: ${{ matrix.GCP_BUCKET }} HOST_CONFIG: ${{ matrix.HOST_CONFIG }} RUNS_ON: ubuntu-22.04 + USE_SCCACHE: false secrets: inherit # If the 'ci: run integrated tests' PR label is found, the integrated tests will be run immediately after the cpu jobs. From 8619b944a5afc98da1007008f00b7d4862d8218a Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 26 Sep 2025 13:02:47 -0500 Subject: [PATCH 08/11] Remove USE_SCCACHE from CI tests configuration Removed USE_SCCACHE environment variable from CI workflow. --- .github/workflows/ci_tests.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml index e8ab48b2596..e1e7506540f 100644 --- a/.github/workflows/ci_tests.yml +++ b/.github/workflows/ci_tests.yml @@ -236,7 +236,6 @@ jobs: GCP_BUCKET: ${{ matrix.GCP_BUCKET }} HOST_CONFIG: ${{ matrix.HOST_CONFIG }} RUNS_ON: ubuntu-22.04 - USE_SCCACHE: false secrets: inherit # If the 'ci: run integrated tests' PR label is found, the integrated tests will be run immediately after the cpu jobs. From 810defa832c2d433cc132dbabbdb6b1e5bd5dd1e Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 26 Sep 2025 13:03:03 -0500 Subject: [PATCH 09/11] Remove unnecessary blank line in Group.hpp --- src/coreComponents/dataRepository/Group.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/coreComponents/dataRepository/Group.hpp b/src/coreComponents/dataRepository/Group.hpp index 38659902b44..0aef2b787b8 100644 --- a/src/coreComponents/dataRepository/Group.hpp +++ b/src/coreComponents/dataRepository/Group.hpp @@ -1752,7 +1752,6 @@ Group::addLogLevel() } } /* end namespace dataRepository */ - } /* end namespace geos */ #endif /* GEOS_DATAREPOSITORY_GROUP_HPP_ */ From 137c6e97aa48cf1898da8ee2bdf345a6b4d29f06 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 29 Sep 2025 18:39:16 -0500 Subject: [PATCH 10/11] fix: matrix setup for poromechanics with contact and wells (#3839) --- .../physicsSolvers/PhysicsSolverBase.cpp | 3 +- .../physicsSolvers/PhysicsSolverBase.hpp | 1 + .../CoupledReservoirAndWellsBase.hpp | 27 ++----------- .../multiphysics/HydrofractureSolver.cpp | 28 ++++--------- .../multiphysics/HydrofractureSolver.hpp | 8 +--- .../PoromechanicsConformingFractures.hpp | 40 ++++++++----------- ...glePhasePoromechanicsEmbeddedFractures.cpp | 1 + ...glePhasePoromechanicsEmbeddedFractures.hpp | 1 + .../physicsSolvers/simplePDE/LaplaceFEM.cpp | 1 + .../physicsSolvers/simplePDE/LaplaceFEM.hpp | 1 + .../SolidMechanicsLagrangianFEM.cpp | 1 + .../SolidMechanicsLagrangianFEM.hpp | 1 + ...lidMechanicsAugmentedLagrangianContact.cpp | 1 + ...lidMechanicsAugmentedLagrangianContact.hpp | 1 + .../SolidMechanicsEmbeddedFractures.cpp | 3 +- .../SolidMechanicsEmbeddedFractures.hpp | 1 + .../contact/SolidMechanicsLagrangeContact.cpp | 3 +- .../contact/SolidMechanicsLagrangeContact.hpp | 1 + ...olidMechanicsLagrangeContactBubbleStab.cpp | 1 + ...olidMechanicsLagrangeContactBubbleStab.hpp | 1 + 20 files changed, 48 insertions(+), 77 deletions(-) diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index 718f935cb6a..3f0d505bf3c 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -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" ); @@ -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 ); diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index e3116f65540..d98d261f7a4 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -411,6 +411,7 @@ class PhysicsSolverBase : public ExecutableGroup virtual void setSparsityPattern( DomainPartition & domain, DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, SparsityPattern< globalIndex > & pattern ); /** diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp index 7af00b2d6d5..31924cfa0f1 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp @@ -106,11 +106,6 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL setInputFlag( dataRepository::InputFlags::FALSE ); } - /** - * @brief default destructor - */ - virtual ~CoupledReservoirAndWellsBase () override {} - /** * @defgroup Solver Interface Functions * @@ -118,30 +113,14 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL */ /**@{*/ - 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()); diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index 06d5f7a816e..e36249bbc76 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -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 ); } @@ -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; @@ -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 > diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp index 9d9e5da29c9..02e2e70d830 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp @@ -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, diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index 4a2f96f6bf6..74d18973e67 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -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()); @@ -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, @@ -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 @@ -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; + } } } } @@ -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 & ) @@ -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 ); + } } } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index c41067a82d2..a327d48e569 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -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. diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp index 23ffc6b6c74..4b7d6fa3437 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp @@ -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 diff --git a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp index 28d1d687c96..009b4d62162 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp +++ b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.cpp @@ -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(), diff --git a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp index 5ca5d971a98..0b3727c5dc6 100644 --- a/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp +++ b/src/coreComponents/physicsSolvers/simplePDE/LaplaceFEM.hpp @@ -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 diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index c2ee4de1214..68b967bfb20 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -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(), diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index b0b2723ef63..12a98378146 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -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 > > diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 2aa85c8b4f5..ebe7df14c9a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -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. diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index c9be177fbc3..c0d87a34bf9 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -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, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp index 4bdca0418d5..0fe6810a514 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp @@ -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; } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp index fbaff010cf0..2834db10804 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp @@ -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, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index f13cefe4f87..5b67c829b29 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -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, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp index b0089d86fcc..674a08876c1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp @@ -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 > > diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp index bd980d76519..893334100a4 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp @@ -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. diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp index c7ef79a835c..9a30e208464 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp @@ -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 From 1ae2748b3d11113344690829340dbea556cf724e Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 30 Sep 2025 10:11:18 -0500 Subject: [PATCH 11/11] Document localMatrix parameter in setSparsityPattern Added parameter documentation for localMatrix in setSparsityPattern. --- src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index d98d261f7a4..9126296ffda 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -406,6 +406,7 @@ class PhysicsSolverBase : public ExecutableGroup * @brief Set the sparsity pattern of the linear system matrix * @param domain the domain containing the mesh and fields * @param dofManager degree-of-freedom manager associated with the linear system + * @param localMatrix the system matrix * @param pattern the sparsity pattern to be filled */ virtual void