feat: add mixed system.

Author: CusiniM

src/reactions/bulkGeneric/MixedEquilibriumKineticReactions.hpp

@@ -0,0 +1,91 @@
+#pragma once
+
+#include "common/macros.hpp"
+
+/** @file MixedEquilibriumKineticReactions.hpp
+ *  @brief Header file for the MixedEquilibriumKineticReactions class.
+ *  @author HPC-REACT Team
+ *  @date 2025
+ */
+
+namespace hpcReact
+{
+namespace bulkGeneric
+{
+
+/**
+ * @brief Class for computing reaction rates and species rates for a given set of reactions.
+ * @tparam REAL_TYPE The type of the real numbers used in the class.
+ * @tparam INT_TYPE The type of the integer numbers used in the class.
+ * @tparam INDEX_TYPE The type of the index used in the class.
+ * @tparam LOGE_CONCENTRATION Whether to use logarithm of concentration for the calculations.
+ * @details
+ *   This class provides the ablity to compute kinetic reactions.
+ */
+template< typename REAL_TYPE,
+          typename INT_TYPE,
+          typename INDEX_TYPE,
+          bool LOGE_CONCENTRATION >
+class MixedEquilibriumKineticReactions
+{
+public:
+
+  /// Type alias for the real type used in the class.
+  using RealType = REAL_TYPE;
+
+  /// Type alias for the integer type used in the class.
+  using IntType = INT_TYPE;
+
+  /// Type alias for the index type used in the class.
+  using IndexType = INDEX_TYPE;
+
+  using kineticReactions = KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, LOGE_CONCENTRATION >;
+
+  template< typename PARAMS_DATA,
+            typename ARRAY_1D_TO_CONST,
+            typename ARRAY_1D,
+            typename ARRAY_2D >
+  static HPCREACT_HOST_DEVICE inline void
+  updateMixedSystem( RealType const & temperature,
+                     PARAMS_DATA const & params,
+                     ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
+                     ARRAY_1D & logSecondarySpeciesConcentrations,
+                     ARRAY_1D & aggregatePrimarySpeciesConcentrations,
+                     ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations,
+                     ARRAY_1D & reactionRates,
+                     ARRAY_2D & reactionRatesDerivatives )
+  {
+    constexpr int numSpecies = PARAMS_DATA::numSpecies;
+    constexpr int numSecondarySpecies = PARAMS_DATA::numReactions;
+    constexpr int numPrimarySpecies = numSpecies - numSecondarySpecies;
+
+    // 1. Compute new aggregate species from primary species
+    calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE,
+                                                    INT_TYPE,
+                                                    INDEX_TYPE >( params, 
+                                                                  logPrimarySpeciesConcentrations,
+                                                                  logSecondarySpeciesConcentrations,
+                                                                  aggregatePrimarySpeciesConcentrations,
+                                                                  dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations )
+
+    // 2. Compute the reaction rates for all kinetic reactions
+    kineticReactions::computeReactionRates( temperature,
+                                            params,
+                                            speciesConcentration,
+                                            reactionRates,
+                                            reactionRatesDerivatives );
+
+    // 3. Compute aggregate species rates
+    
+    
+  }
+
+private:
+
+};
+
+} // namespace bulkGeneric
+} // namespace hpcReact
+
+#include "MixedEquilibriumKineticReactions_impl.hpp"
+#include "common/macrosCleanup.hpp"

src/reactions/bulkGeneric/SpeciesUtilities.hpp

@@ -89,6 +89,51 @@ void calculateLogSecondarySpeciesConcentrationWrtLogC( PARAMS_DATA const & param
   } );
 }
 
+template< typename REAL_TYPE,
+          typename INT_TYPE,
+          typename INDEX_TYPE,
+          typename PARAMS_DATA,
+          typename ARRAY_1D_TO_CONST,
+          typename ARRAY_1D,
+          typename ARRAY_2D >
+HPCREACT_HOST_DEVICE
+inline
+void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
+                                                     ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
+                                                     ARRAY_1D & logSecondarySpeciesConcentrations,
+                                                     ARRAY_1D & aggregatePrimarySpeciesConcentrations,
+                                                     ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations )
+{
+  constexpr int numSpecies = PARAMS_DATA::numSpecies;
+  constexpr int numSecondarySpecies = PARAMS_DATA::numReactions;
+  constexpr int numPrimarySpecies = numSpecies - numSecondarySpecies;
+
+  calculateLogSecondarySpeciesConcentration< REAL_TYPE,
+                                             INT_TYPE,
+                                             INDEX_TYPE >( params,
+                                                           logPrimarySpeciesConcentrations,
+                                                           logSecondarySpeciesConcentrations );
+
+  for( int i = 0; i < numPrimarySpecies; ++i )
+  {
+    REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] );
+    aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
+    dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
+    for( int j = 0; j < numSecondarySpecies; ++j )
+    {
+      REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] );
+      aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
+      for( int k=0; k<numPrimarySpecies; ++k )
+      {
+        REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = params.stoichiometricMatrix( j, k+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
+        dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j,
+                                                                                                                                   i+numSecondarySpecies ) *
+                                                                                                      dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration;
+      }
+    }
+  }
+}
+
 template< typename REAL_TYPE,
           typename INT_TYPE,
           typename INDEX_TYPE,