Skip to content

feat: generalize use of sparse matrices in kernels #3989

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Open
rrsettgast wants to merge 7 commits into
base: develop
Choose a base branch
from
Open

feat: generalize use of sparse matrices in kernels #3989

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
1c660e4
generalize use of sparse matricies for mechanics kernels
rrsettgast Mar 3, 2026
0ff890f
fix some issues on macOS
rrsettgast Mar 4, 2026
636459e
more generalization
rrsettgast Mar 15, 2026
f799518
more generalization
rrsettgast Mar 15, 2026
4eb0808
Merge branch 'develop' into feature/generalizeSparseMatrixUsage
rrsettgast Mar 27, 2026
0bae70e
Merge branch 'develop' into feature/generalizeSparseMatrixUsage
rrsettgast Apr 8, 2026
022ac84
Merge branch 'develop' into feature/generalizeSparseMatrixUsage
rrsettgast Apr 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/coreComponents/codingUtilities/UnitTestUtilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,8 +190,8 @@ void compareMatrices( MATRIX const & matrix1,
ASSERT_EQ( matrix1.numLocalRows(), matrix2.numLocalRows() );
ASSERT_EQ( matrix1.numLocalCols(), matrix2.numLocalCols() );

CRSMatrix< real64, globalIndex > const mat1 = matrix1.extract();
CRSMatrix< real64, globalIndex > const mat2 = matrix2.extract();
DefaultGlobalMatrix const mat1 = matrix1.extract();
DefaultGlobalMatrix const mat2 = matrix2.extract();

compareLocalMatrices( mat1.toViewConst(), mat2.toViewConst(), relTol, absTol, matrix1.ilower() );
}
Expand Down
9 changes: 9 additions & 0 deletions src/coreComponents/common/DataTypes.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,10 +304,19 @@ using SparsityPatternView = LvArray::SparsityPatternView< COL_INDEX, INDEX_TYPE
template< typename T, typename COL_INDEX=globalIndex, typename INDEX_TYPE=localIndex >
using CRSMatrix = LvArray::CRSMatrix< T, COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer >;

/// Default matrix type used for global system assembly and extraction.
using DefaultGlobalMatrix = CRSMatrix< real64, globalIndex >;

/// Alias for CRS Matrix View.
template< typename T, typename COL_INDEX=globalIndex, typename INDEX_TYPE=localIndex >
using CRSMatrixView = LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer >;

/// Default matrix view used for global system assembly kernels and solvers.
using DefaultGlobalMatrixView = CRSMatrixView< real64, globalIndex const >;

/// Default const matrix view used for read-only access to assembled system.
using DefaultGlobalConstMatrixView = CRSMatrixView< real64 const, globalIndex const >;

///@}

//END_SPHINX_INCLUDE_00
Expand Down
6 changes: 4 additions & 2 deletions src/coreComponents/common/FieldSpecificationOps.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -550,11 +550,12 @@ struct FieldSpecificationEqual : public FieldSpecificationOp< OpEqual >
* negate the rhs vector upon assembly. Thus, it sets the value to negative of the desired
* update for the field. For a linear problem, this may lead to unexpected results.
*/
template< typename MATRIX_VIEW = DefaultGlobalMatrixView >
GEOS_HOST_DEVICE
static inline void
SpecifyFieldValue( globalIndex const dof,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
real64 & rhs,
real64 const bcValue,
real64 const fieldValue )
Expand Down Expand Up @@ -640,11 +641,12 @@ struct FieldSpecificationAdd : public FieldSpecificationOp< OpAdd >
* @param[in] fieldValue unused.
*
*/
template< typename MATRIX_VIEW = DefaultGlobalMatrixView >
GEOS_HOST_DEVICE
static inline void
SpecifyFieldValue( globalIndex const dof,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
real64 & rhs,
real64 const bcValue,
real64 const fieldValue )
Expand Down
4 changes: 2 additions & 2 deletions src/coreComponents/dataRepository/python/PyGroup.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,8 +296,8 @@ static PyObject * PyGroup_register( PyGroup * const self, PyObject * const args
return nullptr;
}

std::function< void( CRSMatrix< real64, globalIndex >, array1d< real64 > ) > wrapedCallback =
LvArray::python::PythonFunction< CRSMatrix< real64, globalIndex >, array1d< real64 > > { callback };
std::function< void( DefaultGlobalMatrix, array1d< real64 > ) > wrapedCallback =
LvArray::python::PythonFunction< DefaultGlobalMatrix, array1d< real64 > > { callback };

if( self->group->registerCallback( &wrapedCallback, typeid( wrapedCallback ) ) )
{
Expand Down
3 changes: 2 additions & 1 deletion src/coreComponents/dataRepository/unitTests/testErrorHandling.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include "common/initializeEnvironment.hpp"

#include <gtest/gtest.h>
#include <csignal>
#include <filesystem>
#include <signal.h>

Expand Down Expand Up @@ -445,7 +446,7 @@ TEST( ErrorHandling, VerifySignalHandlerLogs )

ErrorLogger::global().enableFileOutput( true );

raise( SIGINT );
std::raise( SIGINT );

EXPECT_TRUE( signalHappened );

Expand Down
54 changes: 30 additions & 24 deletions src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,13 +206,14 @@ class FieldSpecificationBase : public dataRepository::Group
*
* @note This function is rarely used directly. More often it is called by other ApplyBoundaryCondition functions.
*/
template< typename FIELD_OP, typename POLICY, typename T, int NDIM, int USD >
template< typename FIELD_OP, typename POLICY, typename T, int NDIM, int USD,
typename MATRIX_VIEW = DefaultGlobalMatrixView >
void applyBoundaryConditionToSystemKernel( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs,
ArrayView< T const, NDIM, USD > const & fieldView ) const;

Expand All @@ -235,14 +236,15 @@ class FieldSpecificationBase : public dataRepository::Group
* This function applies the boundary condition to a linear system of equations. This function is
* typically called from within the lambda to a call to BoundaryConditionManager::ApplyBoundaryCondition().
*/
template< typename FIELD_OP, typename POLICY >
template< typename FIELD_OP, typename POLICY,
typename MATRIX_VIEW = DefaultGlobalMatrixView >
void applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
dataRepository::Group const & dataGroup,
string const & fieldName,
string const & dofMapName,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs ) const;

/**
Expand All @@ -267,14 +269,15 @@ class FieldSpecificationBase : public dataRepository::Group
* typically called from within the lambda to a call to
* BoundaryConditionManager::ApplyBoundaryCondition().
*/
template< typename FIELD_OP, typename POLICY, typename LAMBDA >
template< typename FIELD_OP, typename POLICY, typename LAMBDA,
typename MATRIX_VIEW = DefaultGlobalMatrixView >
void
applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs,
LAMBDA && lambda ) const;

Expand All @@ -301,15 +304,16 @@ class FieldSpecificationBase : public dataRepository::Group
* typically called from within the lambda to a call to
* BoundaryConditionManager::ApplyBoundaryCondition().
*/
template< typename FIELD_OP, typename POLICY, typename LAMBDA >
template< typename FIELD_OP, typename POLICY, typename LAMBDA,
typename MATRIX_VIEW = DefaultGlobalMatrixView >
void
applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
real64 const dt,
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs,
LAMBDA && lambda ) const;

Expand Down Expand Up @@ -342,15 +346,16 @@ class FieldSpecificationBase : public dataRepository::Group
* Therefore, the compositional solvers do not call applyBoundaryConditionToSystem, but instead call computeRhsContribution directly, and
* apply these rhs contributions "manually" according to the equation layout used in the solver
*/
template< typename FIELD_OP, typename POLICY, typename LAMBDA >
template< typename FIELD_OP, typename POLICY, typename LAMBDA,
typename MATRIX_VIEW = DefaultGlobalMatrixView >
void
computeRhsContribution( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
real64 const dt,
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< globalIndex > const & dof,
arrayView1d< real64 > const & rhsContribution,
LAMBDA && lambda ) const;
Expand All @@ -366,10 +371,11 @@ class FieldSpecificationBase : public dataRepository::Group
*
* This function zeroes the rows of the matrix that correspond to boundary conditions.
*/
template< typename POLICY >
template< typename POLICY,
typename MATRIX_VIEW = DefaultGlobalMatrixView >
void zeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet,
arrayView1d< globalIndex const > const & dofMap,
CRSMatrixView< real64, globalIndex const > const & matrix ) const;
MATRIX_VIEW const & matrix ) const;

/**
* @brief View keys
Expand Down Expand Up @@ -693,13 +699,13 @@ void FieldSpecificationBase::applyFieldValue( SortedArrayView< localIndex const
}, wrapper );
}

template< typename FIELD_OP, typename POLICY, typename T, int NDIM, int USD >
template< typename FIELD_OP, typename POLICY, typename T, int NDIM, int USD, typename MATRIX_VIEW >
void FieldSpecificationBase::applyBoundaryConditionToSystemKernel( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs,
ArrayView< T const, NDIM, USD > const & fieldView ) const
{
Expand All @@ -713,14 +719,14 @@ void FieldSpecificationBase::applyBoundaryConditionToSystemKernel( SortedArrayVi
} );
}

template< typename FIELD_OP, typename POLICY >
template< typename FIELD_OP, typename POLICY, typename MATRIX_VIEW >
void FieldSpecificationBase::applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
dataRepository::Group const & dataGroup,
string const & fieldName,
string const & dofMapName,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs ) const
{
dataRepository::WrapperBase const & wrapper = dataGroup.getWrapperBase( fieldName );
Expand All @@ -743,15 +749,15 @@ void FieldSpecificationBase::applyBoundaryConditionToSystem( SortedArrayView< lo
}, wrapper );
}

template< typename FIELD_OP, typename POLICY, typename LAMBDA >
template< typename FIELD_OP, typename POLICY, typename LAMBDA, typename MATRIX_VIEW >
void
FieldSpecificationBase::
applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
real64 const time,
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs,
LAMBDA && lambda ) const
{
Expand All @@ -766,7 +772,7 @@ FieldSpecificationBase::
std::forward< LAMBDA >( lambda ) );
}

template< typename FIELD_OP, typename POLICY, typename LAMBDA >
template< typename FIELD_OP, typename POLICY, typename LAMBDA, typename MATRIX_VIEW >
void
FieldSpecificationBase::
applyBoundaryConditionToSystem( SortedArrayView< localIndex const > const & targetSet,
Expand All @@ -775,7 +781,7 @@ FieldSpecificationBase::
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< real64 > const & rhs,
LAMBDA && lambda ) const
{
Expand All @@ -799,7 +805,7 @@ FieldSpecificationBase::
FIELD_OP::template prescribeRhsValues< POLICY >( rhs, dof, dofRankOffset, rhsContribution );
}

template< typename FIELD_OP, typename POLICY, typename LAMBDA >
template< typename FIELD_OP, typename POLICY, typename LAMBDA, typename MATRIX_VIEW >
void
FieldSpecificationBase::
computeRhsContribution( SortedArrayView< localIndex const > const & targetSet,
Expand All @@ -808,7 +814,7 @@ FieldSpecificationBase::
dataRepository::Group const & dataGroup,
arrayView1d< globalIndex const > const & dofMap,
globalIndex const dofRankOffset,
CRSMatrixView< real64, globalIndex const > const & matrix,
MATRIX_VIEW const & matrix,
arrayView1d< globalIndex > const & dof,
arrayView1d< real64 > const & rhsContribution,
LAMBDA && lambda ) const
Expand Down Expand Up @@ -868,10 +874,10 @@ FieldSpecificationBase::
}


template< typename POLICY >
template< typename POLICY, typename MATRIX_VIEW >
void FieldSpecificationBase::zeroSystemRowsForBoundaryCondition( SortedArrayView< localIndex const > const & targetSet,
arrayView1d< globalIndex const > const & dofMap,
CRSMatrixView< real64, globalIndex const > const & matrix ) const
MATRIX_VIEW const & matrix ) const

{
integer const component = ( getComponent() >=0 ) ? getComponent() : 0;
Expand Down
7 changes: 4 additions & 3 deletions src/coreComponents/finiteElement/kernelInterface/ImplicitKernelBase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,8 @@ template< typename SUBREGION_TYPE,
typename CONSTITUTIVE_TYPE,
typename FE_TYPE,
int NUM_DOF_PER_TEST_SP,
int NUM_DOF_PER_TRIAL_SP >
int NUM_DOF_PER_TRIAL_SP,
typename MATRIX_VIEW = DefaultGlobalMatrixView >
class ImplicitKernelBase : public KernelBase< SUBREGION_TYPE,
CONSTITUTIVE_TYPE,
FE_TYPE,
Expand Down Expand Up @@ -92,7 +93,7 @@ class ImplicitKernelBase : public KernelBase< SUBREGION_TYPE,
CONSTITUTIVE_TYPE & inputConstitutiveType,
arrayView1d< globalIndex const > const & inputDofNumber,
globalIndex const rankOffset,
CRSMatrixView< real64, globalIndex const > const & inputMatrix,
MATRIX_VIEW const & inputMatrix,
arrayView1d< real64 > const & inputRhs,
real64 const inputDt ):
Base( elementSubRegion,
Expand Down Expand Up @@ -217,7 +218,7 @@ class ImplicitKernelBase : public KernelBase< SUBREGION_TYPE,
globalIndex const m_dofRankOffset;

/// The global Jacobian matrix.
CRSMatrixView< real64, globalIndex const > const m_matrix;
MATRIX_VIEW const m_matrix;

/// The global residaul vector.
arrayView1d< real64 > const m_rhs;
Expand Down
Loading
Loading