feat: generalize use of sparse matrices in kernels

src/coreComponents/codingUtilities/UnitTestUtilities.hpp

@@ -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() );
 }

src/coreComponents/common/DataTypes.hpp

@@ -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

src/coreComponents/common/FieldSpecificationOps.hpp

@@ -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 )
@@ -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 )

src/coreComponents/dataRepository/python/PyGroup.cpp

@@ -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 ) ) )
   {

src/coreComponents/dataRepository/unitTests/testErrorHandling.cpp

@@ -25,6 +25,7 @@
 #include "common/initializeEnvironment.hpp"
 
 #include <gtest/gtest.h>
+#include <csignal>
 #include <filesystem>
 #include <signal.h>
 
@@ -445,7 +446,7 @@ TEST( ErrorHandling, VerifySignalHandlerLogs )
 
   ErrorLogger::global().enableFileOutput( true );
 
-  raise( SIGINT );
+  std::raise( SIGINT );
 
   EXPECT_TRUE( signalHappened );
 

src/coreComponents/fieldSpecification/FieldSpecificationBase.hpp

@@ -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;
 
@@ -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;
 
   /**
@@ -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;
 
@@ -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;
 
@@ -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;
@@ -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
@@ -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
 {
@@ -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 );
@@ -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
 {
@@ -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,
@@ -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
 {
@@ -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,
@@ -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
@@ -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;

src/coreComponents/finiteElement/kernelInterface/ImplicitKernelBase.hpp

@@ -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,
@@ -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,
@@ -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;