feat: add mixed equilibrium and kinetic reactions

hostconfigs/apple/macOS_base.cmake

@@ -1,5 +1,5 @@
 message( "this hostconfig assumes you are using homebrew")
-message( "brew install bison cmake gcc git-lfs open-mpi openblas python3")
+message( "brew install bison cmake gcc git-lfs open-mpi openblas python3 llvm cppcheck lcov")
 
 message( "CMAKE_SYSTEM_PROCESSOR = ${CMAKE_SYSTEM_PROCESSOR}" )
 message("CONFIG_NAME = ${CONFIG_NAME}")

src/CMakeLists.txt

@@ -7,13 +7,14 @@ set( hpcReact_headers
      reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp
      reactions/bulkGeneric/KineticReactions.hpp
      reactions/bulkGeneric/KineticReactions_impl.hpp
+     reactions/bulkGeneric/MixedEquilibriumKineticReactions.hpp
+     reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp
      reactions/bulkGeneric/Parameters.hpp
      reactions/bulkGeneric/ParametersPredefined.hpp
      reactions/bulkGeneric/SpeciesUtilities.hpp
      )
 
-set( hpcReact_sources
-     )
+set( hpcReact_sources)
 
 find_package(LAPACK REQUIRED)
 find_package(BLAS REQUIRED)
@@ -35,11 +36,13 @@ blt_add_library( NAME             hpcReact
 
 target_include_directories( hpcReact
                             INTERFACE
-                            $<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/include>
-                            $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src>
+                            $<BUILD_INTERFACE:${CMAKE_BINARY_DIR}../include>
+                            $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
                             $<INSTALL_INTERFACE:include> )
 
 blt_print_target_properties( TARGET hpcReact )
+message(STATUS "HPCReact/src CMAKE_CURRENT_SOURCE_DIR: ${CMAKE_CURRENT_SOURCE_DIR}")
+
 # install( FILES ${hpcReact_headers} 
 #          DESTINATION include )
 

src/common/CArrayWrapper.hpp

@@ -1,5 +1,6 @@
 #pragma once
 #include "macros.hpp"
+#include <initializer_list>
 
 
 /**
@@ -26,37 +27,72 @@ struct CArrayWrapper;
 template< typename T, int DIM0 >
 struct CArrayWrapper< T, DIM0 >
 {
+  // default constructor
+  constexpr CArrayWrapper() = default;
+
+  /**
+   * @brief Construct a CArrayWrapper from an initializer list.
+   *
+   * Allows brace-initialization with a list of values:
+   * @code
+   * CArrayWrapper< double, 3 > arr = {1.0, 2.0, 3.0};
+   * @endcode
+   *
+   * @param init An initializer list with exactly DIM0 elements.
+   *
+   * @note No runtime bounds checking is performed on the initializer size.
+   */
+  constexpr CArrayWrapper( std::initializer_list< T > init )
+  {
+    // static_assert(init.size() == DIM0, "Size mismatch"); // needs c++20
+    int i = 0;
+    for( auto const & val : init )
+    {
+      data[i++] = val;
+    }
+  }
+  /**
+   * @brief Copy constructor.
+   * @param src The source CArrayWrapper to copy from.
+   */ 
+  constexpr CArrayWrapper( CArrayWrapper const & src )
+  {
+    for( std::size_t i = 0; i < DIM0; i++ )
+    {
+      data[i] = src.data[i];
+    }
+  }
 
   /**
    * @brief Read/write access to an element by index.
    * @param dim The index (must be in range [0, DIM0)).
    * @return Reference to the element at the specified index.
    */
-  HPCREACT_HOST_DEVICE inline T & operator()( int const dim ) { return data[dim]; }
+  HPCREACT_HOST_DEVICE constexpr inline T & operator()( int const dim ) { return data[dim]; }
 
   /**
    * @brief Read-only access to an element by index (const overload).
    * @param dim The index (must be in range [0, DIM0)).
    * @return Const reference to the element at the specified index.
    */
-  HPCREACT_HOST_DEVICE inline T const & operator()( int const dim ) const { return data[dim]; }
+  HPCREACT_HOST_DEVICE constexpr inline T const & operator()( int const dim ) const { return data[dim]; }
 
   /**
    * @brief Subscript operator for read/write access.
    * @param dim The index (must be in range [0, DIM0)).
    * @return Reference to the element at the specified index.
    */
-  HPCREACT_HOST_DEVICE inline T & operator[]( int const dim ) { return data[dim]; }
+  HPCREACT_HOST_DEVICE constexpr inline T & operator[]( int const dim ) { return data[dim]; }
 
   /**
    * @brief Subscript operator for read-only access (const overload).
    * @param dim The index (must be in range [0, DIM0)).
    * @return Const reference to the element at the specified index.
    */
-  HPCREACT_HOST_DEVICE inline T const & operator[]( int const dim ) const { return data[dim]; }
+  HPCREACT_HOST_DEVICE constexpr inline T const & operator[]( int const dim ) const { return data[dim]; }
 
   /// The underlying 1D C-style array.
-  T data[DIM0];
+  T data[DIM0]{};
 };
 
 /**
@@ -72,13 +108,60 @@ struct CArrayWrapper< T, DIM0 >
 template< typename T, int DIM0, int DIM1 >
 struct CArrayWrapper< T, DIM0, DIM1 >
 {
+  // default constructor
+  constexpr CArrayWrapper() = default;
+
+  /**
+   * @brief Copy constructor.
+   * @param src The source CArrayWrapper to copy from.
+   */
+  constexpr CArrayWrapper( CArrayWrapper const & src )
+  {
+    for( std::size_t i = 0; i < DIM0; i++ )
+    {
+      for ( std::size_t j = 0; j < DIM1; j++)
+        data[i][j] = src.data[i][j];
+    }
+  }
+
+  /**
+   * @brief Construct a 2D CArrayWrapper from nested initializer lists.
+   *
+   * Allows brace-initialization with a matrix-like structure:
+   * @code
+   * CArrayWrapper<double, 2, 3> mat = {
+   *   {1.0, 2.0, 3.0},
+   *   {4.0, 5.0, 6.0}
+   * };
+   * @endcode
+   *
+   * @param init A nested initializer list with exactly D0 rows and D1 elements per row.
+   *
+   * @note No runtime bounds checking is performed on the initializer dimensions.
+   */
+  constexpr CArrayWrapper( std::initializer_list< std::initializer_list< T > > init )
+  {
+    // static_assert(init.size() == DIM0, "Size mismatch"); // needs c++20
+    int i = 0;
+    for( auto const & row : init )
+    {
+      // static_assert(row.size() == DIM1, "Size mismatch"); // needs c++20
+      int j = 0;
+      for( auto const & val : row )
+      {
+        data[i][j++] = val;
+      }
+      ++i;
+    }
+  }
+
   /**
    * @brief Read/write access to an element by 2D indices.
    * @param dim0 Index in the first dimension (range [0, DIM0)).
    * @param dim1 Index in the second dimension (range [0, DIM1)).
    * @return Reference to the element at the specified 2D location.
    */
-  HPCREACT_HOST_DEVICE inline T & operator()( int const dim0, int const dim1 )
+  HPCREACT_HOST_DEVICE constexpr inline T & operator()( int const dim0, int const dim1 )
   {
     return data[dim0][dim1];
   }
@@ -89,7 +172,7 @@ struct CArrayWrapper< T, DIM0, DIM1 >
    * @param dim1 Index in the second dimension (range [0, DIM1)).
    * @return Const reference to the element at the specified 2D location.
    */
-  HPCREACT_HOST_DEVICE inline T const & operator()( int const dim0, int const dim1 ) const
+  HPCREACT_HOST_DEVICE constexpr inline T const & operator()( int const dim0, int const dim1 ) const
   {
     return data[dim0][dim1];
   }
@@ -101,7 +184,7 @@ struct CArrayWrapper< T, DIM0, DIM1 >
    *
    * This allows usage like `obj[dim0][dim1]`.
    */
-  HPCREACT_HOST_DEVICE inline T       ( & operator[]( int const dim0 ))[DIM1]
+  HPCREACT_HOST_DEVICE constexpr inline T       ( & operator[]( int const dim0 ))[DIM1]
   {
     return data[dim0];
   }
@@ -111,13 +194,13 @@ struct CArrayWrapper< T, DIM0, DIM1 >
    * @param dim0 The row index (range [0, DIM0)).
    * @return Const reference to an array of type T[DIM1].
    */
-  HPCREACT_HOST_DEVICE inline T const (&operator[]( int const dim0 ) const)[DIM1]
+  HPCREACT_HOST_DEVICE constexpr inline T const (&operator[]( int const dim0 ) const)[DIM1]
   {
     return data[dim0];
   }
 
   /// The underlying 2D C-style array of size DIM0 x DIM1.
-  T data[DIM0][DIM1];
+  T data[DIM0][DIM1]{};
 };
 
 /**
@@ -134,6 +217,54 @@ struct CArrayWrapper< T, DIM0, DIM1 >
 template< typename T, int DIM0, int DIM1, int DIM2 >
 struct CArrayWrapper< T, DIM0, DIM1, DIM2 >
 {
+  // default constructor
+  constexpr CArrayWrapper() = default;
+
+  /**
+   * @brief Construct a 3D CArrayWrapper from nested initializer lists.
+   *
+   * Enables tensor-like initialization using triple-nested braces:
+   * @code
+   * CArrayWrapper<double, 2, 2, 2> cube = {
+   *   {
+   *     {1.0, 2.0},
+   *     {3.0, 4.0}
+   *   },
+   *   {
+   *     {5.0, 6.0},
+   *     {7.0, 8.0}
+   *   }
+   * };
+   * @endcode
+   *
+   * @param init A three-level nested initializer list with D0 planes, D1 rows per plane,
+   *             and D2 elements per row.
+   *
+   * @note This constructor does not perform size validation. Incorrect initializer sizes
+   *       may lead to undefined behavior.
+   */
+  constexpr CArrayWrapper( std::initializer_list< std::initializer_list< std::initializer_list< T > > > init )
+  {
+    // static_assert(init.size() == DIM0, "Size mismatch"); // needs c++20
+    int i = 0;
+    for( auto const & plane : init )
+    {
+      // static_assert(plane.size() == DIM1, "Size mismatch"); // needs c++20
+      int j = 0;
+      for( auto const & row : plane )
+      {
+        // static_assert(row.size() == DIM2, "Size mismatch"); // needs c++20
+        int k = 0;
+        for( auto const & val : row )
+        {
+          data[i][j][k++] = val;
+        }
+        ++j;
+      }
+      ++i;
+    }
+  }
+
   /**
    * @brief Read/write access to an element by 3D indices.
    * @param dim0 Index in the first dimension (range [0, DIM0)).
@@ -144,7 +275,7 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 >
    * @note Currently, this function incorrectly indexes data[dim0][dim1], missing dim2.
    *       It should be `data[dim0][dim1][dim2]`. Please correct if intended.
    */
-  HPCREACT_HOST_DEVICE inline T & operator()( int const dim0, int const dim1, int const dim2 )
+  HPCREACT_HOST_DEVICE constexpr inline T & operator()( int const dim0, int const dim1, int const dim2 )
   {
     // NOTE: This looks like a bug in your original code. Should be data[dim0][dim1][dim2].
     return data[dim0][dim1][dim2];
@@ -157,7 +288,7 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 >
    * @param dim2 Index in the third dimension (range [0, DIM2)).
    * @return Const reference to the element at the specified 3D location.
    */
-  HPCREACT_HOST_DEVICE inline T const & operator()( int const dim0, int const dim1, int const dim2 ) const
+  HPCREACT_HOST_DEVICE constexpr inline T const & operator()( int const dim0, int const dim1, int const dim2 ) const
   {
     // NOTE: Same potential bug as above. Should be data[dim0][dim1][dim2].
     return data[dim0][dim1][dim2];
@@ -170,7 +301,7 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 >
    *
    * This allows usage like `obj[dim0][dim1][dim2]`.
    */
-  HPCREACT_HOST_DEVICE inline T       ( & operator[]( int const dim0 ))[DIM1][DIM2]
+  HPCREACT_HOST_DEVICE constexpr inline T ( & operator[]( int const dim0 ))[DIM1][DIM2]
   {
     return data[dim0];
   }
@@ -180,11 +311,11 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 >
    * @param dim0 The slice index (range [0, DIM0)).
    * @return Const reference to an array of type T[DIM1][DIM2].
    */
-  HPCREACT_HOST_DEVICE inline T const (&operator[]( int const dim0 ) const)[DIM1][DIM2]
+  HPCREACT_HOST_DEVICE constexpr inline T const (&operator[]( int const dim0 ) const)[DIM1][DIM2]
   {
     return data[dim0];
   }
 
   /// The underlying 3D C-style array of size DIM0 x DIM1 x DIM2.
-  T data[DIM0][DIM1][DIM2];
+  T data[DIM0][DIM1][DIM2]{};
 };

src/reactions/bulkGeneric/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp

@@ -30,9 +30,7 @@ EquilibriumReactions< REAL_TYPE,
                                                                                               ARRAY_2D & jacobian )
 {
   HPCREACT_UNUSED_VAR( temperature );
-  constexpr int numSpecies = PARAMS_DATA::numSpecies;
-  constexpr int numSecondarySpecies = PARAMS_DATA::numReactions;
-  constexpr int numPrimarySpecies = numSpecies - numSecondarySpecies;
+  constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
 
   RealType aggregatePrimaryConcentrations[numPrimarySpecies] = {0.0};
   ARRAY_2D dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations = {{{0.0}}};
@@ -68,9 +66,7 @@ EquilibriumReactions< REAL_TYPE,
                                                                   ARRAY_1D & logPrimarySpeciesConcentration )
 {
   HPCREACT_UNUSED_VAR( temperature );
-  constexpr int numSpecies = PARAMS_DATA::numSpecies;
-  constexpr int numReactions = PARAMS_DATA::numReactions;
-  constexpr int numPrimarySpecies = numSpecies - numReactions;
+  constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
 
   double residual[numPrimarySpecies] = { 0.0 };
 //  double aggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 };

src/reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp

@@ -32,8 +32,8 @@ EquilibriumReactions< REAL_TYPE,
 {
 
   HPCREACT_UNUSED_VAR( temperature );
-  constexpr int numSpecies = PARAMS_DATA::numSpecies;
-  constexpr int numReactions = PARAMS_DATA::numReactions;
+  constexpr int numSpecies = PARAMS_DATA::numSpecies();
+  constexpr int numReactions = PARAMS_DATA::numReactions();
 
   // initialize the species concentration
   RealType speciesConcentration[numSpecies];
@@ -115,8 +115,8 @@ EquilibriumReactions< REAL_TYPE,
                                                                 ARRAY_1D & speciesConcentration )
 {
   HPCREACT_UNUSED_VAR( temperature );
-  constexpr int numSpecies = PARAMS_DATA::numSpecies;
-  constexpr int numReactions = PARAMS_DATA::numReactions;
+  constexpr int numSpecies = PARAMS_DATA::numSpecies();
+  constexpr int numReactions = PARAMS_DATA::numReactions();
   double residual[numReactions] = { 0.0 };
   double xi[numReactions] = { 0.0 };
   double dxi[numReactions] = { 0.0 };

src/reactions/bulkGeneric/KineticReactions.hpp

@@ -79,7 +79,7 @@ class KineticReactions
                         ARRAY_1D_TO_CONST const & speciesConcentration,
                         ARRAY_1D & reactionRates )
   {
-    REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions][PARAMS_DATA::numSpecies] = { {0.0} };
+    REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions()][PARAMS_DATA::numSpecies()] = { {0.0} };
     computeReactionRates_impl< PARAMS_DATA, false >( temperature,
                                                      params,
                                                      speciesConcentration,