Generic Chain Reaction for Kinetic Reaction Verification (#13)

frankfeifan · web-flow · commit 2b2aeaf270ef · 2025-10-29T17:09:43.000-07:00
* cleaned up momas and added medium test case

* added 3 species full kinetic serial chain reaction

* added an option variable for user-defined kinetic rate calculation method

* modified computeReactionRatesTest and added a new test for chain reaction

* added a kinetic reaction test to call quotient reaction rate compute
diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp
@@ -63,7 +63,9 @@ simpleKineticTestRateParams =
   // reverse rate constants
   { 1.0, 0.5 },
   // flag of mobile secondary species
-  { 1, 1 }
+  { 1, 1 },
+  // Use the forward and reverse to calculate the kinetic reaction rates
+  0
 };
 
 using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >;
@@ -84,7 +86,9 @@ simpleTestRateParams =
   // reverse rate constants
   { 1.0, 0.5 },
   // flag of mobile secondary species
-  { 1, 1 }
+  { 1, 1 },
+  // Use the forward and reverse to calculate the kinetic reaction rates
+  0
 };
 
 // *****UNCRUSTIFY-ON******
diff --git a/src/reactions/exampleSystems/ChainGeneric.hpp b/src/reactions/exampleSystems/ChainGeneric.hpp
@@ -0,0 +1,68 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: (BSD-3-Clause)
+ *
+ * Copyright (c) 2025- Lawrence Livermore National Security LLC
+ * All rights reserved
+ *
+ * See top level LICENSE files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+#pragma once
+
+#include "../reactionsSystems/Parameters.hpp"
+
+namespace hpcReact
+{
+namespace ChainGeneric
+{
+// *****UNCRUSTIFY-OFF******
+
+  using serialAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 3, 3, 0 >;
+
+  constexpr serialAllKineticType serialAllKineticParams =
+  {
+    // Stoichiometric matrix [3 rows × 3 columns]
+    // Columns 0–3 are primary species: {C1, C2, C3 }
+    {
+      // C1   C2   C3 
+      {  -1,   1,   0 },  // C1 = C2
+      {   0,  -1,   1 },  // C2 = C3
+      {   0,   0,  -1 },  // C3 = 
+    },
+
+    // Equilibrium constants K
+    {
+      0, // C1 = C2
+      0, // C2 = C3
+      0  // C3
+    },
+
+    // Forward rate constants 
+    { 
+      0.05, // C1 = C2
+      0.03, // C2 = C3
+      0.02, // C3
+    },
+
+    // Reverse rate constants 
+    { 
+      0.0, // C1 = C2
+      0.0, // C2 = C3 
+      0.0  // C3
+    },
+
+    // Flag of mobile secondary species
+    { 
+      1, // C1 = C2
+      1, // C2 = C3 
+      1  // C3
+    },
+
+    0 // Use the forward and reverse to calculate the kinetic reaction rates
+  };
+
+// *****UNCRUSTIFY-ON******
+} // namespace ChainGeneric
+} // namespace hpcReact
diff --git a/src/reactions/exampleSystems/unitTests/CMakeLists.txt b/src/reactions/exampleSystems/unitTests/CMakeLists.txt
@@ -4,6 +4,7 @@ set( testSourceFiles
      testKineticReactions.cpp
      testMomasEasyCase.cpp
      testMomasMediumCase.cpp
+     testGenericChainReactions.cpp
    )
 
 set( dependencyList hpcReact gtest )
diff --git a/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp
@@ -0,0 +1,64 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: (BSD-3-Clause)
+ *
+ * Copyright (c) 2025- Lawrence Livermore National Security LLC
+ * All rights reserved
+ *
+ * See top level LICENSE files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+#include "reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp"
+#include "../ChainGeneric.hpp"
+
+
+using namespace hpcReact;
+using namespace hpcReact::unitTest_utilities;
+
+//******************************************************************************
+TEST( testChainGenericKineticReactions, computeReactionRatesTest_chainReactionParams )
+{
+  using namespace hpcReact::ChainGeneric;
+
+  double const initialSpeciesConcentration[] =
+  {
+    1.0, // C1
+    1e-8, // C2
+    1e-8 // C3
+  };
+  double const surfaceArea[] =
+  {
+    0.0, // C1 -> C2
+    0.0, // C2 -> C3
+    0.0 // C3 ->
+  };
+
+  double const expectedReactionRates[] =
+  {
+    0.05, // C1 -> C2
+    3e-10, // C2 -> C3
+    2e-10 // C3 ->
+  };
+  double const expectedReactionRatesDerivatives[][3] =
+  { { 0.05, 0.0, 0.0 },
+    { 0.0, 0.03, 0.0 },
+    { 0.0, 0.0, 0.02 } };
+  computeReactionRatesTest< double, false >( serialAllKineticParams.kineticReactionsParameters(),
+                                             initialSpeciesConcentration,
+                                             surfaceArea, // No use. Just to pass something here
+                                             expectedReactionRates,
+                                             expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double, true >( serialAllKineticParams.kineticReactionsParameters(),
+                                            initialSpeciesConcentration,
+                                            surfaceArea, // No use. Just to pass something here
+                                            expectedReactionRates,
+                                            expectedReactionRatesDerivatives );
+}
+
+int main( int argc, char * * argv )
+{
+  ::testing::InitGoogleTest( &argc, argv );
+  int const result = RUN_ALL_TESTS();
+  return result;
+}
diff --git a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp
@@ -24,16 +24,19 @@ using namespace hpcReact::unitTest_utilities;
 TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams )
 {
   double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
+  double const surfaceArea[] = { 0.0, 0.0 };
   double const expectedReactionRates[] = { 1.0, 0.25 };
   double const expectedReactionRatesDerivatives[][5] =
   { { 2.0, -0.5, 0.0, 0.0, 0.0 },
     { 0.0, 0.0, 0.5, 0.25, 0.0 } };
   computeReactionRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
                                              initialSpeciesConcentration,
+                                             surfaceArea, // No use. Just to pass something here
                                              expectedReactionRates,
                                              expectedReactionRatesDerivatives );
   computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(),
                                             initialSpeciesConcentration,
+                                            surfaceArea, // No use. Just to pass something here
                                             expectedReactionRates,
                                             expectedReactionRatesDerivatives );
 }
diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp
@@ -113,7 +113,7 @@ using carbonateSystemAllKineticType     = reactionsSystems::MixedReactionsParame
 using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 10, 10 >;
 using carbonateSystemType               = reactionsSystems::MixedReactionsParameters< double, int, int, 16, 10, 9 >;
 
-constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
+constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag, 0 );
 constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
 constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
 
diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp
@@ -43,6 +43,18 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic )
     1.09 // Na+1
   };
 
+  double const surfaceArea[10] = { 0.0, // OH- + H+ = H2O
+                                   0.0, // CO2 + H2O = H+ + HCO3-
+                                   0.0, // CO3-2 + H+ = HCO3-
+                                   0.0, // CaHCO3+ = Ca+2 + HCO3-
+                                   0.0, // CaSO4 = Ca+2 + SO4-2
+                                   0.0, // CaCl+ = Ca+2 + Cl-
+                                   0.0, // CaCl2 = Ca+2 + 2Cl-
+                                   0.0, // MgSO4 = Mg+2 + SO4-2
+                                   0.0, // NaSO4- = Na+ + SO4-2
+                                   0.0, // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
+  };
+
   double const expectedReactionRates[10] = { -0.001424736, //             OH- + H+ = H2O
                                              -12610.7392, //              CO2 + H2O = H+ + HCO3-
                                              -0.175591624, //             CO3-2 + H+ = HCO3-
@@ -71,15 +83,58 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic )
 
   computeReactionRatesTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(),
                                              initialSpeciesConcentration,
+                                             surfaceArea, // No use. Just to pass something here
                                              expectedReactionRates,
                                              expectedReactionRatesDerivatives );
   computeReactionRatesTest< double, true >( carbonateSystemAllKinetic.kineticReactionsParameters(),
                                             initialSpeciesConcentration,
+                                            surfaceArea, // No use. Just to pass something here
                                             expectedReactionRates,
                                             expectedReactionRatesDerivatives );
 }
 
+TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem )
+{
+  double const initialSpeciesConcentration[16] =
+  {
+    1.0e-16, // OH-
+    1.0e-16, // CO2
+    1.0e-16, // CO3-2
+    1.0e-16, // CaHCO3+
+    1.0e-16, // CaSO4
+    1.0e-16, // CaCl+
+    1.0e-16, // CaCl2
+    1.0e-16, // MgSO4
+    1.0e-16, // NaSO4-
+    3.76e-1, // H+
+    3.76e-1, // HCO3-
+    3.87e-2, // Ca+2
+    3.21e-2, // SO4-2
+    1.89, // Cl-
+    1.65e-2, // Mg+2
+    1.09 // Na+1
+  };
+
+  double const surfaceArea[1] = { 1e6 }; // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
+
+  double const expectedReactionRates[1] = { 1.5488389999999999 }; // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
+
+  double const expectedReactionRatesDerivatives[1][16] =
+  {
+    { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.0877659574468075e-03, -3.0877659574468075e-03, -2.9999999999999997e-02, 0, 0, 0, 0 }
+  };
 
+  computeReactionRatesTest< double, false >( carbonateSystem.kineticReactionsParameters(),
+                                             initialSpeciesConcentration,
+                                             surfaceArea,
+                                             expectedReactionRates,
+                                             expectedReactionRatesDerivatives );
+  computeReactionRatesTest< double, true >( carbonateSystem.kineticReactionsParameters(),
+                                            initialSpeciesConcentration,
+                                            surfaceArea,
+                                            expectedReactionRates,
+                                            expectedReactionRatesDerivatives );
+}
 
 //******************************************************************************
 
diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp
@@ -72,6 +72,11 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params,
                                                 ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
                                                 ARRAY_1D & logSecondarySpeciesConcentrations )
 {
+  if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 )
+  {
+    return;
+  }
+
   massActions_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE,
                                                                INT_TYPE,
                                                                INDEX_TYPE >( params,
@@ -173,15 +178,22 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
 {
   static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
 
-  REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0};
+  if constexpr( numSecondarySpecies > 0 )
+  {
+    REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0};
 
-  calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE,
-                                                  INT_TYPE,
-                                                  INDEX_TYPE >( params,
-                                                                logPrimarySpeciesConcentrations,
-                                                                logSecondarySpeciesConcentrations,
-                                                                aggregatePrimarySpeciesConcentrations,
-                                                                dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations );
+    calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE,
+                                                    INT_TYPE,
+                                                    INDEX_TYPE >( params,
+                                                                  logPrimarySpeciesConcentrations,
+                                                                  logSecondarySpeciesConcentrations,
+                                                                  aggregatePrimarySpeciesConcentrations,
+                                                                  dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations );
+  }
+  else
+  {
+    GEOS_UNUSED_VAR( logPrimarySpeciesConcentrations, aggregatePrimarySpeciesConcentrations, dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations );
+  }
 
 }
 
diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp
@@ -113,6 +113,11 @@ EquilibriumReactions< REAL_TYPE,
                                                                   ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0,
                                                                   ARRAY_1D & logPrimarySpeciesConcentration )
 {
+  if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 )
+  {
+    return;
+  }
+
   HPCREACT_UNUSED_VAR( temperature );
   static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
 
diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp
@@ -13,6 +13,8 @@
 
 #include "common/macros.hpp"
 
+#include <stdexcept>
+
 /** @file KineticReactions.hpp
  *  @brief Header file for the KineticReactions class.
  *  @author HPC-REACT Team
@@ -130,12 +132,23 @@ class KineticReactions
                         ARRAY_1D & reactionRates,
                         ARRAY_2D & reactionRatesDerivatives )
   {
-    computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature,
-                                                            params,
-                                                            speciesConcentration,
-                                                            surfaceArea,
-                                                            reactionRates,
-                                                            reactionRatesDerivatives );
+    if( params.reactionRatesUpdateOption() == 0 )
+    {
+      computeReactionRates_impl< PARAMS_DATA, true >( temperature,
+                                                      params,
+                                                      speciesConcentration,
+                                                      reactionRates,
+                                                      reactionRatesDerivatives );
+    }
+    else if( params.reactionRatesUpdateOption() == 1 )
+    {
+      computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature,
+                                                              params,
+                                                              speciesConcentration,
+                                                              surfaceArea,
+                                                              reactionRates,
+                                                              reactionRatesDerivatives );
+    }
   }
 
 
diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp
@@ -271,7 +271,7 @@ KineticReactions< REAL_TYPE,
         RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], s_ri ) : 0.0;
         quotient *= productTerm_i;
       }
-      
+
       if constexpr( CALCULATE_DERIVATIVES )
       {
         for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp
diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp
diff --git a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp

Original file line number	Diff line number	Diff line change
`@@ -4,6 +4,7 @@ set( testSourceFiles`
`4`	`4`	`testKineticReactions.cpp`
`5`	`5`	`testMomasEasyCase.cpp`
`6`	`6`	`testMomasMediumCase.cpp`
	`7`	`+ testGenericChainReactions.cpp`
`7`	`8`	`)`
`8`	`9`
`9`	`10`	`set( dependencyList hpcReact gtest )`