modified computeReactionRatesTest and added a new test for chain reaction

Author: frankfeifan

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

src/reactions/exampleSystems/unitTests/CMakeLists.txt

@@ -4,6 +4,7 @@ set( testSourceFiles
      testKineticReactions.cpp
      testMomasEasyCase.cpp
      testMomasMediumCase.cpp
+     testGenericChainReactions.cpp
    )
 
 set( dependencyList hpcReact gtest )

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;
+}

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

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

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,10 +83,12 @@ 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 );
 }

src/reactions/massActions/MassActions.hpp

@@ -76,7 +76,7 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params,
   {
     return;
   }
-  
+
   massActions_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE,
                                                                INT_TYPE,
                                                                INDEX_TYPE >( params,

src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp

@@ -115,9 +115,9 @@ EquilibriumReactions< REAL_TYPE,
 {
   if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 )
   {
-    return; 
+    return;
   }
-  
+
   HPCREACT_UNUSED_VAR( temperature );
   constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
 

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,27 @@ 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 );
+    }
+    else
+    {
+      throw std::runtime_error( "Error: No such option implemented so far for reaction rate calculation." );
+    }
   }
 
 

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 )