GEOS-DEV · rrsettgast · Sep 3, 2025 · Jul 14, 2025 · Jul 25, 2025 · Jul 28, 2025
diff --git a/cmake/Macros.cmake b/cmake/Macros.cmake
@@ -42,6 +42,7 @@ macro(hpcReact_add_code_checks)
                         --suppress=unusedFunction 
                         --suppress=constStatement 
                         --suppress=unusedStructMember
+                        -DHPCREACT_HOST_DEVICE= 
                         -I../hpcReact/src )
 
     if( UNCRUSTIFY_FOUND )

diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp
@@ -61,7 +61,9 @@ simpleKineticTestRateParams =
   // forward rate constants
   { 1.0, 0.5 },
   // reverse rate constants
-  { 1.0, 0.5 }
+  { 1.0, 0.5 },
+  // flag of mobile secondary species
+  { 1, 1 }
 };
 
 using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >;
@@ -80,7 +82,9 @@ simpleTestRateParams =
   // forward rate constants
   { 1.0, 0.5 },
   // reverse rate constants
-  { 1.0, 0.5 }
+  { 1.0, 0.5 },
+  // flag of mobile secondary species
+  { 1, 1 }
 };
 
 // *****UNCRUSTIFY-ON******

diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp
@@ -28,9 +28,9 @@ namespace MomMasBenchmark
   // Columns 7–11 are primary species: {X1, X2, X3, X4, S}
   {
     // C1   C2   C3   C4   C5   CS1  CS2  |  X1  X2  X3   X4  S
-    { -1,   0,   0,   0,   0,   0,   0,      0,  -1,  0,  0,  0 },  // C1 = X2
-    {  0,  -1,   0,   0,   0,   0,   0,      0,   1,  1,  0,  0 },  // C2 = X2 * X3
-    {  0,   0,  -1,   0,   0,   0,   0,      0,  -1,  0,  1,  0 },  // C3 = X2 * X4
+    { -1,   0,   0,   0,   0,   0,   0,      0,  -1,  0,  0,  0 },  // C1 = -X2
+    {  0,  -1,   0,   0,   0,   0,   0,      0,   1,  1,  0,  0 },  // C2 = X2 + X3
+    {  0,   0,  -1,   0,   0,   0,   0,      0,  -1,  0,  1,  0 },  // C3 = X2 + X4
     {  0,   0,   0,  -1,   0,   0,   0,      0,  -4,  1,  3,  0 },  // C4 = -4X2 + X3 + 3X4
     {  0,   0,   0,   0,  -1,   0,   0,      0,   4,  3,  1,  0 },  // C5 = 4X2 + 3X3 + X4
     {  0,   0,   0,   0,   0,  -1,   0,      0,   3,  1,  0,  1 },  // CS1 = 3X2 + X3 + S
@@ -39,34 +39,46 @@ namespace MomMasBenchmark
 
   // Equilibrium constants K
   {
-    1.0e-12,   // C1
-    1.0,       // C2
-    1.0,       // C3
-    0.1,       // C4
-    1.0e35,    // C5
-    1.0e6,     // CS1
-    1.0e-1     // CS2
+    1.0e12,   //   C1 + X2 = ??
+    1.0,      //        C2 = X2 + X3
+    1.0,      //        C3 = X2 + X4
+    1.0e1,    //  C4 + 4X2 = X3 + 3X4
+    1.0e-35,  //        C5 = 4X2 + 3X3 + X4
+    1.0e-6,   //       CS1 = 3X2 + X3 + S
+    1.0e1     // CS2 + 3X2 = + X4 + 2S
   },
 
   // Forward rate constants 
-  { 0.0,
-    0.0, 
-    0.0, 
-    0.0, 
-    0.0, 
-    0.0, 
-    0.0 
+  { 
+    0.0,   // C1 = -X2
+    0.0,   // C2 = X2 + X3
+    0.0,   // C3 = X2 + X4
+    0.0,   // C4 = -4X2 + X3 + 3X4
+    0.0,   // C5 = 4X2 + 3X3 + X4
+    0.0,   // CS1 = 3X2 + X3 + S
+    0.0    // CS2 = -3X2 + X4 + 2S
   },
 
   // Reverse rate constants 
-  { 0.0,
-    0.0, 
-    0.0, 
-    0.0, 
-    0.0, 
-    0.0, 
-    0.0 
+  { 
+    0.0,   // C1 = -X2
+    0.0,   // C2 = X2 + X3
+    0.0,   // C3 = X2 + X4
+    0.0,   // C4 = -4X2 + X3 + 3X4
+    0.0,   // C5 = 4X2 + 3X3 + X4
+    0.0,   // CS1 = 3X2 + X3 + S
+    0.0    // CS2 = -3X2 + X4 + 2S
   },
+
+  // Flag of mobile secondary species
+  { 1, // C1 = -X2
+    1, // C2 = X2 + X3
+    1, // C3 = -X2 + X4
+    1, // C4 = -4X2 + X3 + 3X4
+    1, // C5 = 4X2 + 3X3 + X4
+    0, // CS1 = 3X2 + X3 + S
+    0  // CS2 = -3X2 + X4 + 2S
+  }
 };
 
 // *****UNCRUSTIFY-ON******

diff --git a/src/reactions/exampleSystems/unitTests/CMakeLists.txt b/src/reactions/exampleSystems/unitTests/CMakeLists.txt
@@ -2,6 +2,7 @@
 set( testSourceFiles
      testEquilibriumReactions.cpp
      testKineticReactions.cpp
+     testMomasEasyCase.cpp
    )
 
 set( dependencyList hpcReact gtest )

diff --git a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp
@@ -0,0 +1,87 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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/equilibriumReactionsTestUtilities.hpp"
+#include "../MoMasBenchmark.hpp"
+
+using namespace hpcReact;
+using namespace hpcReact::MomMasBenchmark;
+using namespace hpcReact::unitTest_utilities;
+
+//******************************************************************************
+
+TEST( testEquilibriumReactions, testMoMasAllEquilibrium )
+{
+
+  using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double,
+                                                                           int,
+                                                                           int >;
+
+  constexpr int numPrimarySpecies = hpcReact::MomMasBenchmark::simpleSystemParams.numPrimarySpecies();
+
+  double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] =
+  {
+    1.0e-20, // X1
+    -2.0, // X2
+    1.0e-20, // X3
+    2.0, // X4
+    1.0 // S
+  };
+
+  double const initialPrimarySpeciesConcentration[numPrimarySpecies] =
+  {
+    1.0e-20, // X1
+    0.02, // X2
+    1.0e-20, // X3
+    1.0, // X4
+    1.00 // S
+  };
+
+  double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
+  {
+    log( initialPrimarySpeciesConcentration[0] ),
+    log( initialPrimarySpeciesConcentration[1] ),
+    log( initialPrimarySpeciesConcentration[2] ),
+    log( initialPrimarySpeciesConcentration[3] ),
+    log( initialPrimarySpeciesConcentration[4] )
+  };
+
+  double logPrimarySpeciesConcentration[numPrimarySpecies];
+  EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
+                                                          hpcReact::MomMasBenchmark::simpleSystemParams.equilibriumReactionsParameters(),
+                                                          targetAggregatePrimarySpeciesConcentration,
+                                                          logInitialPrimarySpeciesConcentration,
+                                                          logPrimarySpeciesConcentration );
+
+  double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
+  {
+    9.9999999999999919e-21, // X1
+    0.25971841330881928, // X2
+    1.4603613417111526e-24, // X3
+    0.3495378685828045, // X4
+    0.39074371811222675 // S
+  };
+
+  for( int r=0; r<numPrimarySpecies; ++r )
+  {
+    EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
+  }
+
+
+}
+
+int main( int argc, char * * argv )
+{
+  ::testing::InitGoogleTest( &argc, argv );
+  int const result = RUN_ALL_TESTS();
+  return result;
+}
diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp
@@ -101,17 +101,31 @@ constexpr CArrayWrapper<double, 11> reverseRates =
     2.51E+06,   //      CaCl2 = Ca+2 + 2Cl- 
     2.69E+07,   //      MgSO4 = Mg+2 + SO4-2
     6.62E+07,   //     NaSO4- = Na+ + SO4-2
-    8.55E-03    // CaCO3 + H+ = Ca+2 + HCO3-
+    8.55E-03   // CaCO3 + H+ = Ca+2 + HCO3-
     //  1           //   Ca(OH)2(s) + 2H+ = Ca2+ + 2H2O (kinetic)
   };
+
+constexpr CArrayWrapper<int, 11> mobileSpeciesFlag = 
+  { 1,   //   OH- + H+ = H2O         
+    1,   //  CO2 + H2O = H+ + HCO3-  
+    1,   // CO3-2 + H+ = HCO3-       
+    1,   //      H2CO3 = H+ + HCO3-  
+    1,   //    CaHCO3+ = Ca+2 + HCO3-
+    1,   //      CaSO4 = Ca+2 + SO4-2
+    1,   //      CaCl+ = Ca+2 + Cl-  
+    1,   //      CaCl2 = Ca+2 + 2Cl- 
+    1,   //      MgSO4 = Mg+2 + SO4-2
+    1,   //     NaSO4- = Na+ + SO4-2
+    1   // CaCO3 + H+ = Ca+2 + HCO3-
+  };
   }
   using carbonateSystemAllKineticType     = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 0 >;
   using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 11 >;
   using carbonateSystemType               = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 11, 10 >;
 
-  constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
-  constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
-  constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
+  constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
+  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 );
 
 // *****UNCRUSTIFY-ON******
 } // namespace geochemistry

diff --git a/src/reactions/geochemistry/Forge.hpp b/src/reactions/geochemistry/Forge.hpp
@@ -119,12 +119,35 @@ constexpr CArrayWrapper< double, 19 > reverseRateConstant =
    1.0    // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq)
 };
 
+constexpr CArrayWrapper< int, 19 > mobileSpeciesFlag =
+{
+   1,   // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻
+   1,   // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻
+   1,   // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻
+   1,   // CaCl⁺ ⇌ Ca²⁺ + Cl⁻
+   1,   // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source)
+   1,   // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻
+   1,   // MgCO₃(aq) + H⁺ ⇌ Mg²⁺ + HCO₃⁻
+   1,   // MgCl⁺ ⇌ Mg²⁺ + Cl⁻
+   1,   // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻
+   1,   // HSO₄⁻ ⇌ H⁺ + SO₄²⁻
+   1,   // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻
+   1,   // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq)
+   1,   // NaHSilO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq)
+   1,   // NaCl ⇌ Na⁺ + Cl⁻
+   1,   // KCl ⇌ K⁺ + Cl⁻
+   1,   // KSO₄⁻ ⇌ K⁺ + SO₄²⁻
+   1,   // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻
+   1,   // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq)
+   1    // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq)
+};
+
 }
 
 using forgeSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 26, 19, 16 >;
 
 
-constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant );
+constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant, forge::mobileSpeciesFlag );
 
 
 // *****UNCRUSTIFY-ON******

diff --git a/src/reactions/geochemistry/Ultramafics.hpp b/src/reactions/geochemistry/Ultramafics.hpp
@@ -126,15 +126,40 @@ constexpr CArrayWrapper<double, 21> reverseRates =
     2.83E-44,   //  Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O
     2.10E-25    //  Mg(OH)2 + 2H+ = Mg++ + 2H2O
   };
-};
+
+constexpr CArrayWrapper<int, 21> mobileSpeciesFlag = 
+  { 
+    1,   //  OH- + H+ = H2O         
+    1,   //  CO2(aq) + H2O = HCO3- + H+  
+    1,   //  CO3-- + H+ = HCO3-       
+    1,   //  Mg2OH+++ + H+ = 2Mg++ + H2O
+    1,   //  Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O
+    1,   //  MgOH+ + H+ = Mg++ + H2O
+    1,   //  Mg2CO3++ + H+ = 2Mg++ + HCO3-
+    1,   //  MgCO3 + H+ = Mg++ + HCO3-
+    1,   //  MgHCO3+ = Mg++ + HCO3-
+    1,   //  Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O
+    1,   //  MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O
+    1,   //  MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O
+    1,   //  H2SiO4-- + 2H+ = SiO2(aq) + 2H2O
+    1,   //  H3SiO4- + H+ = SiO2(aq) + 2H2O
+    1,   //  H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O
+    1,   //  H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O
+    1,   //  Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O
+    1,   //  MgCO3 + H+ = Mg++ + HCO3-
+    1,   //  SiO2 = SiO2(aq)
+    1,   //  Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O
+    1    //  Mg(OH)2 + 2H+ = Mg++ + 2H2O
+  };
+}
 
   using ultramaficSystemAllKineticType     = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 0 >;
   using ultramaficSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 21 >;
   using ultramaficSystemType               = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 16 >;
 
-  constexpr ultramaficSystemAllKineticType     ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
-  constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
-  constexpr ultramaficSystemType               ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
+  constexpr ultramaficSystemAllKineticType     ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
+  constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
+  constexpr ultramaficSystemType               ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
 
 // *****UNCRUSTIFY-ON******
 } // namespace geochemistry

diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp
@@ -133,10 +133,10 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 )
   };
 
   double logPrimarySpeciesConcentration[numPrimarySpecies];
-  EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
-                                                          hpcReact::geochemistry::carbonateSystemAllEquilibrium.equilibriumReactionsParameters(),
-                                                          logInitialPrimarySpeciesConcentration,
-                                                          logPrimarySpeciesConcentration );
+  EquilibriumReactionsType::enforceEquilibrium_LogAggregate( 0,
+                                                             hpcReact::geochemistry::carbonateSystemAllEquilibrium.equilibriumReactionsParameters(),
+                                                             logInitialPrimarySpeciesConcentration,
+                                                             logPrimarySpeciesConcentration );
 
   double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
   {