add system from Tiras

Author: rrsettgast

src/reactions/bulkGeneric/KineticReactions_impl.hpp

@@ -3,6 +3,10 @@
 #include "common/macros.hpp"
 
 #include <cmath>
+#include <string>
+#include <iostream>
+
+#define OUTPUT_RATE_EQUATIONS 0
 
 namespace hpcReact
 {
@@ -27,45 +31,88 @@ KineticReactions< REAL_TYPE,
                                            RealConstDataArrayView1d const & primarySpeciesConcentration,
                                            RealDataArrayView1d & reactionRates )
 {
+#if OUTPUT_RATE_EQUATIONS==1
+  std::string const c[] = { "A", "B", "C", "D", "E" };
+#endif
+
   for( IntType i = 0; i < PARAMS_DATA::numPrimarySpecies; ++i )
   {
+#if OUTPUT_RATE_EQUATIONS==1
+    std::string rateString;
+#endif
     reactionRates[i] = 0.0;
     for( IntType r=0; r<PARAMS_DATA::numKineticReactions; ++r )
     {
       RealType const forwardRateConstant = params.m_rateConstant[r] * exp( -params.m_activationEnergy[r] / ( constants::R * temperature ) );
-      RealType const reverseRateConstant = params.m_equilibriumConstant[r] / forwardRateConstant;
+      RealType const reverseRateConstant = forwardRateConstant / params.m_equilibriumConstant[r] ;
+
+//      printf( "forwardRateConstant = %8.4e, reverseRateConstant = %8.4e\n", forwardRateConstant, reverseRateConstant );
       RealType const s_ir = params.m_stoichiometricMatrix[r][i];
 
       RealType productConcPlus = 1.0;
       RealType productConcMinus = 1.0;
       bool productConcPlusFlag = false;
       bool productConcMinusFlag = false;
+#if OUTPUT_RATE_EQUATIONS==1
+      std::string productConcPlusString;
+      std::string productConcMinusString;
+#endif
       for( IntType j = 0; j < PARAMS_DATA::numPrimarySpecies; ++j )
       {
         RealType const s_jr = params.m_stoichiometricMatrix[r][j];
         if( s_jr < 0.0 )
         {
           productConcPlus *= pow( primarySpeciesConcentration[j], -s_jr );
+#if OUTPUT_RATE_EQUATIONS==1
+          productConcPlusString += " c(" +c[j]+")" ;//+ "^" + std::to_string(-s_jr) + " ";
+#endif
           productConcPlusFlag = true;
         }
         else if( s_jr > 0.0 )
         {
           productConcMinus *= pow( primarySpeciesConcentration[j], s_jr );
+#if OUTPUT_RATE_EQUATIONS==1
+          productConcMinusString += " c("+c[j]+")" ;//+ "^" + std::to_string(s_jr) + " ";
+#endif
           productConcMinusFlag = true;
         }
       }
 
       RealType reactionRate_n = 0.0;
+#if OUTPUT_RATE_EQUATIONS==1
+      std::string reactionRateString = "{ ";
+#endif
       if( productConcPlusFlag )
       {
         reactionRate_n += forwardRateConstant * productConcPlus;
+#if OUTPUT_RATE_EQUATIONS==1
+        reactionRateString += " k_" + std::to_string(r+1) + "f" + productConcPlusString;
+#endif
       }
       if( productConcMinusFlag )
       {
         reactionRate_n -= reverseRateConstant * productConcMinus;
+#if OUTPUT_RATE_EQUATIONS==1
+        reactionRateString += " - k_" + std::to_string(r+1) + "r" + productConcMinusString;
+#endif
       }
+
+#if OUTPUT_RATE_EQUATIONS==1
+      reactionRateString += " }";
+      if( s_ir > 0.1 || s_ir < -0.1 )
+      {
+        std::string s_irString;
+        s_irString.resize(10);
+        snprintf(s_irString.data(), 5, "%+.1f", s_ir);
+        rateString += s_irString + " " + reactionRateString;
+      }
+#endif
       reactionRates[i] += s_ir * reactionRate_n;
+
     }
+#if OUTPUT_RATE_EQUATIONS==1
+    std::cout<<rateString<<std::endl;
+#endif
   }
 }
 

src/reactions/bulkGeneric/Parameters.hpp

@@ -14,6 +14,7 @@ namespace bulkGeneric
 
 template< typename REAL_TYPE,
           typename INT_TYPE,
+          typename INDEX_TYPE,
           int NUM_PRIMARY_SPECIES,
           int NUM_KINETIC_REACTIONS >
 struct KineticParameters
@@ -26,11 +27,29 @@ struct KineticParameters
 
   RealType m_activationEnergy[numKineticReactions];
   RealType m_equilibriumConstant[numKineticReactions];
-  CArrayWrapper< RealType, numKineticReactions, numPrimarySpecies>  m_stoichiometricMatrix;
+  RealType m_stoichiometricMatrix[numKineticReactions][numPrimarySpecies];
 
   RealType m_rateConstant[numKineticReactions];
 };
 
+template< typename REAL_TYPE,
+          typename INT_TYPE,
+          typename INDEX_TYPE,
+          int NUM_PRIMARY_SPECIES,
+          int NUM_SECONDARY_SPECIES,
+          int NUM_KINETIC_REACTIONS,
+          int NUM_EQUILIBRIUM_REACTIONS >
+struct BulkParameters
+{
+  using RealType = REAL_TYPE;
+  using IntType = INT_TYPE;
+
+  static constexpr IntType numPrimarySpecies = NUM_PRIMARY_SPECIES;
+  static constexpr IntType numSecondarySpecies = NUM_SECONDARY_SPECIES;
+  static constexpr IntType numKineticReactions = NUM_KINETIC_REACTIONS;
+  static constexpr IntType numEquilibriumReactions = NUM_EQUILIBRIUM_REACTIONS;
+
+};
 
 
 

src/reactions/bulkGeneric/ParametersPredefined.hpp

@@ -8,7 +8,7 @@ namespace bulkGeneric
 {
 
 constexpr
-KineticParameters< double, int, 5, 5 > bicarbonateBuffer =
+KineticParameters< double, int, int, 6, 5 > bicarbonateBuffer =
 {
   // taken from 
   // Simplified Models of the Bicarbonate Buffer for Scaled Simulations of CO2 Electrolyzers
@@ -22,15 +22,59 @@ KineticParameters< double, int, 5, 5 > bicarbonateBuffer =
   // m_equilibriumConstant
     { 4.27E+07, 4.27E-07, 2.09E-04, 2.09E+10, 1.00E+14 },
   // m_stoichiometricMatrix
-  { { {  1,  1, -1,  0,  0 },
-      {  1,  0, -1,  0, -1 },
-      {  0, -1, -1,  1,  0 },
-      {  0,  0, -1,  1,  1 },
-      {  0,  1,  0,  0,  1 } }
+  { { -1, -1,  1,  0,  0,  0 },
+    { -1,  0,  1,  0,  1, -1 },
+    {  0,  1,  1, -1,  0,  0 },
+    {  0,  0,  1, -1, -1,  0 },
+    {  0, -1,  0,  0, -1,  1 }
   },
   // m_rateConstant
   { 8.42E+03, 3.71E-02, 1.25E+06, 1.24E+12, 2.30E+10 }
 };
 
+
+constexpr
+KineticParameters< double, int, int, 6, 5 > bicarbonateBufferTest =
+{
+  // taken from 
+  // Simplified Models of the Bicarbonate Buffer for Scaled Simulations of CO2 Electrolyzers
+  // Thomas Moore, Tiras Y. Lin, Thomas Roy, Sarah E. Baker, Eric B. Duoss, Christopher Hahn, and Victor A. Beck
+  // Industrial & Engineering Chemistry Research 2023 62 (40), 16291-16301
+  // DOI: 10.1021/acs.iecr.3c02504
+
+  // m_activationEnergy
+    { 0.0, 0.0, 0.0, 0.0, 0.0 },
+  // m_equilibriumConstant
+    { 4.27E+07, 4.27E-07, 2.09E-04, 2.09E+10, 1.00E+14 },
+  // m_stoichiometricMatrix
+  // C02, OH-, HCO3-, CO3^2-, H+, H20
+  { { -1, -1,  1,  0,  0,  0 },
+    { -1,  0,  1,  0,  1, -1 },
+    {  0,  1,  1, -1,  0,  0 },
+    {  0,  0,  1, -1, -1,  0 },
+    {  0, -1,  0,  0, -1,  1 }
+  },
+  // m_rateConstant
+  { 8.42E+03, 3.71E-02, 1.25E+03, 1.24E+3, 2.30E+3 }
+};
+
+constexpr
+KineticParameters< double, int, int, 5, 2 > simpleTest =
+{
+  // m_activationEnergy
+    { 0.0, 0.0 },
+  // m_equilibriumConstant
+    { 1.0, 1.0 },
+  // m_stoichiometricMatrix
+  // C02, OH-, HCO3-, CO3^2-, H+, H20
+  { { -2,  1,  1,  0,  0 },
+    {  0,  0, -1, -1,  2 }
+  },
+  // m_rateConstant
+  { 1.0, 0.5 }
+};
+
+
+
 } // namespace bulkGeneric
 } // namespace hpcReact

src/reactions/bulkGeneric/unitTests/testKineticReactions.cpp

@@ -8,7 +8,39 @@
 using namespace hpcReact;
 using namespace hpcReact::bulkGeneric;
 
-TEST( bulkGeneric, test_computeReactionRates )
+// TEST( bulkGeneric, test_computeReactionRates )
+// {
+//   using KineticReactionsType = KineticReactions< double, 
+//                                                  double * const,
+//                                                  double const * const, 
+//                                                  int, 
+//                                                  int >;
+
+//   double const temperature = 298.15;
+//   double primarySpeciesConcentration[6] = { 0.01, 0.01, 0.01, 0.01, 0.01, 1.0 };
+//   double reactionRates[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+  
+//   KineticReactionsType::computeReactionRates( temperature, 
+//                                               bicarbonateBuffer, 
+//                                               primarySpeciesConcentration, 
+//                                               reactionRates );
+
+//   printf( "reactionRates = {%8.4e, %8.4e, %8.4e, %8.4e, %8.4e, %8.4e }\n", 
+//           reactionRates[0],
+//           reactionRates[1],
+//           reactionRates[2],
+//           reactionRates[3],
+//           reactionRates[4],
+//           reactionRates[5] );
+//   // EXPECT_NEAR( reactionRates[0], 0.0, 1.0e-8 );
+//   // EXPECT_NEAR( reactionRates[1], 0.0, 1.0e-8 );
+//   // EXPECT_NEAR( reactionRates[2], 0.0, 1.0e-8 );
+//   // EXPECT_NEAR( reactionRates[3], 0.0, 1.0e-8 );
+//   // EXPECT_NEAR( reactionRates[4], 0.0, 1.0e-8 );
+  
+// }
+
+TEST( bulkGeneric, test_computeReactionRatesIntegral )
 {
   using KineticReactionsType = KineticReactions< double, 
                                                  double * const,
@@ -17,20 +49,44 @@ TEST( bulkGeneric, test_computeReactionRates )
                                                  int >;
 
   double const temperature = 298.15;
-  double const primarySpeciesConcentration[] = { 1.0, 1.0, 1.0, 1.0, 1.0 };
+  double primarySpeciesConcentration[] = { 1.0, 0.0, 0.5, 1.0, 0.0 };
   double reactionRates[] = { 0.0, 0.0, 0.0, 0.0, 0.0 };
-  
-  KineticReactionsType::computeReactionRates( temperature, 
-                                              bicarbonateBuffer, 
-                                              primarySpeciesConcentration, 
-                                              reactionRates );
-
-  EXPECT_NEAR( reactionRates[0], 0.0, 1.0e-8 );
-  EXPECT_NEAR( reactionRates[1], 0.0, 1.0e-8 );
-  EXPECT_NEAR( reactionRates[2], 0.0, 1.0e-8 );
-  EXPECT_NEAR( reactionRates[3], 0.0, 1.0e-8 );
-  EXPECT_NEAR( reactionRates[4], 0.0, 1.0e-8 );
-  
+
+  constexpr double dt = 1.0e-1;
+  constexpr int numPrimarySpecies = decltype(simpleTest)::numPrimarySpecies;
+//  constexpr int numKineticReactions = decltype(simpleTest)::numKineticReactions;
+
+  //double time = 0.0;
+  for( int i = 0; i < 1000; ++i )
+  {
+    KineticReactionsType::computeReactionRates( temperature, 
+                                                simpleTest, 
+                                                primarySpeciesConcentration, 
+                                                reactionRates );
+
+    for( int j = 0; j < numPrimarySpecies; ++j )
+    {
+      primarySpeciesConcentration[j] += dt * reactionRates[j];
+    }
+    //time += dt;
+    // if( i % 10 == 0 )
+    // {
+    //   printf( " {% 12.8e} { % 12.8e, % 12.8e, % 12.8e, % 12.8e, % 12.8e }\n", 
+    //           time,
+    //           primarySpeciesConcentration[0],
+    //           primarySpeciesConcentration[1],
+    //           primarySpeciesConcentration[2],
+    //           primarySpeciesConcentration[3],
+    //           primarySpeciesConcentration[4] );
+    // }
+  }
+
+  EXPECT_NEAR( primarySpeciesConcentration[0], 3.92138294e-01, 1.0e-8 );
+  EXPECT_NEAR( primarySpeciesConcentration[1], 3.03930853e-01, 1.0e-8 );
+  EXPECT_NEAR( primarySpeciesConcentration[2], 5.05945481e-01, 1.0e-8 );
+  EXPECT_NEAR( primarySpeciesConcentration[3], 7.02014628e-01, 1.0e-8 );
+  EXPECT_NEAR( primarySpeciesConcentration[4], 5.95970745e-01, 1.0e-8 );
+
 }
 
 int main( int argc, char * * argv )