Skip to content

Commit f20e43a

Browse files
frankfeifanjkgollarrsettgast
authored
feat: MoMaS Benchmark (#11)
* Implement unit test for MoMaS benchmark easy case * cleanup of massActions...and change definition of K in the momas example * add new interface function to specify target aggregate concentrations for equilibrium. Fix momas test * move the momas example to exampleSystems * changed back the carbonate parameters and leave it to another PR --------- Co-authored-by: jkgolla <[email protected]> Co-authored-by: Randolph Settgast <[email protected]>
1 parent 057bb05 commit f20e43a

16 files changed

+387
-69
lines changed

‎cmake/Macros.cmake‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,7 @@ macro(hpcReact_add_code_checks)
4242
--suppress=unusedFunction
4343
--suppress=constStatement
4444
--suppress=unusedStructMember
45+
-DHPCREACT_HOST_DEVICE=
4546
-I../hpcReact/src )
4647

4748
if( UNCRUSTIFY_FOUND )

‎src/reactions/exampleSystems/BulkGeneric.hpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,9 @@ simpleKineticTestRateParams =
6161
// forward rate constants
6262
{ 1.0, 0.5 },
6363
// reverse rate constants
64-
{ 1.0, 0.5 }
64+
{ 1.0, 0.5 },
65+
// flag of mobile secondary species
66+
{ 1, 1 }
6567
};
6668

6769
using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >;
@@ -80,7 +82,9 @@ simpleTestRateParams =
8082
// forward rate constants
8183
{ 1.0, 0.5 },
8284
// reverse rate constants
83-
{ 1.0, 0.5 }
85+
{ 1.0, 0.5 },
86+
// flag of mobile secondary species
87+
{ 1, 1 }
8488
};
8589

8690
// *****UNCRUSTIFY-ON******

‎src/reactions/exampleSystems/MoMasBenchmark.hpp‎

Lines changed: 36 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -28,9 +28,9 @@ namespace MomMasBenchmark
2828
// Columns 7–11 are primary species: {X1, X2, X3, X4, S}
2929
{
3030
// C1 C2 C3 C4 C5 CS1 CS2 | X1 X2 X3 X4 S
31-
{ -1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0 }, // C1 = X2
32-
{ 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0 }, // C2 = X2 * X3
33-
{ 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 1, 0 }, // C3 = X2 * X4
31+
{ -1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0 }, // C1 = -X2
32+
{ 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0 }, // C2 = X2 + X3
33+
{ 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 1, 0 }, // C3 = X2 + X4
3434
{ 0, 0, 0, -1, 0, 0, 0, 0, -4, 1, 3, 0 }, // C4 = -4X2 + X3 + 3X4
3535
{ 0, 0, 0, 0, -1, 0, 0, 0, 4, 3, 1, 0 }, // C5 = 4X2 + 3X3 + X4
3636
{ 0, 0, 0, 0, 0, -1, 0, 0, 3, 1, 0, 1 }, // CS1 = 3X2 + X3 + S
@@ -39,34 +39,46 @@ namespace MomMasBenchmark
3939

4040
// Equilibrium constants K
4141
{
42-
1.0e-12, // C1
43-
1.0, // C2
44-
1.0, // C3
45-
0.1, // C4
46-
1.0e35, // C5
47-
1.0e6, // CS1
48-
1.0e-1 // CS2
42+
1.0e12, // C1 + X2 = ??
43+
1.0, // C2 = X2 + X3
44+
1.0, // C3 = X2 + X4
45+
1.0e1, // C4 + 4X2 = X3 + 3X4
46+
1.0e-35, // C5 = 4X2 + 3X3 + X4
47+
1.0e-6, // CS1 = 3X2 + X3 + S
48+
1.0e1 // CS2 + 3X2 = + X4 + 2S
4949
},
5050

5151
// Forward rate constants
52-
{ 0.0,
53-
0.0,
54-
0.0,
55-
0.0,
56-
0.0,
57-
0.0,
58-
0.0
52+
{
53+
0.0, // C1 = -X2
54+
0.0, // C2 = X2 + X3
55+
0.0, // C3 = X2 + X4
56+
0.0, // C4 = -4X2 + X3 + 3X4
57+
0.0, // C5 = 4X2 + 3X3 + X4
58+
0.0, // CS1 = 3X2 + X3 + S
59+
0.0 // CS2 = -3X2 + X4 + 2S
5960
},
6061

6162
// Reverse rate constants
62-
{ 0.0,
63-
0.0,
64-
0.0,
65-
0.0,
66-
0.0,
67-
0.0,
68-
0.0
63+
{
64+
0.0, // C1 = -X2
65+
0.0, // C2 = X2 + X3
66+
0.0, // C3 = X2 + X4
67+
0.0, // C4 = -4X2 + X3 + 3X4
68+
0.0, // C5 = 4X2 + 3X3 + X4
69+
0.0, // CS1 = 3X2 + X3 + S
70+
0.0 // CS2 = -3X2 + X4 + 2S
6971
},
72+
73+
// Flag of mobile secondary species
74+
{ 1, // C1 = -X2
75+
1, // C2 = X2 + X3
76+
1, // C3 = -X2 + X4
77+
1, // C4 = -4X2 + X3 + 3X4
78+
1, // C5 = 4X2 + 3X3 + X4
79+
0, // CS1 = 3X2 + X3 + S
80+
0 // CS2 = -3X2 + X4 + 2S
81+
}
7082
};
7183

7284
// *****UNCRUSTIFY-ON******

‎src/reactions/exampleSystems/unitTests/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
set( testSourceFiles
33
testEquilibriumReactions.cpp
44
testKineticReactions.cpp
5+
testMomasEasyCase.cpp
56
)
67

78
set( dependencyList hpcReact gtest )

‎src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
/*
2+
* ------------------------------------------------------------------------------------------------------------
3+
* SPDX-License-Identifier: (BSD-3-Clause)
4+
*
5+
* Copyright (c) 2025- Lawrence Livermore National Security LLC
6+
* All rights reserved
7+
*
8+
* See top level LICENSE files for details.
9+
* ------------------------------------------------------------------------------------------------------------
10+
*/
11+
12+
13+
#include "reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp"
14+
#include "../MoMasBenchmark.hpp"
15+
16+
using namespace hpcReact;
17+
using namespace hpcReact::MomMasBenchmark;
18+
using namespace hpcReact::unitTest_utilities;
19+
20+
//******************************************************************************
21+
22+
TEST( testEquilibriumReactions, testMoMasAllEquilibrium )
23+
{
24+
25+
using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double,
26+
int,
27+
int >;
28+
29+
constexpr int numPrimarySpecies = hpcReact::MomMasBenchmark::simpleSystemParams.numPrimarySpecies();
30+
31+
double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] =
32+
{
33+
1.0e-20, // X1
34+
-2.0, // X2
35+
1.0e-20, // X3
36+
2.0, // X4
37+
1.0 // S
38+
};
39+
40+
double const initialPrimarySpeciesConcentration[numPrimarySpecies] =
41+
{
42+
1.0e-20, // X1
43+
0.02, // X2
44+
1.0e-20, // X3
45+
1.0, // X4
46+
1.00 // S
47+
};
48+
49+
double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
50+
{
51+
log( initialPrimarySpeciesConcentration[0] ),
52+
log( initialPrimarySpeciesConcentration[1] ),
53+
log( initialPrimarySpeciesConcentration[2] ),
54+
log( initialPrimarySpeciesConcentration[3] ),
55+
log( initialPrimarySpeciesConcentration[4] )
56+
};
57+
58+
double logPrimarySpeciesConcentration[numPrimarySpecies];
59+
EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
60+
hpcReact::MomMasBenchmark::simpleSystemParams.equilibriumReactionsParameters(),
61+
targetAggregatePrimarySpeciesConcentration,
62+
logInitialPrimarySpeciesConcentration,
63+
logPrimarySpeciesConcentration );
64+
65+
double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
66+
{
67+
9.9999999999999919e-21, // X1
68+
0.25971841330881928, // X2
69+
1.4603613417111526e-24, // X3
70+
0.3495378685828045, // X4
71+
0.39074371811222675 // S
72+
};
73+
74+
for( int r=0; r<numPrimarySpecies; ++r )
75+
{
76+
EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
77+
}
78+
79+
80+
}
81+
82+
int main( int argc, char * * argv )
83+
{
84+
::testing::InitGoogleTest( &argc, argv );
85+
int const result = RUN_ALL_TESTS();
86+
return result;
87+
}

‎src/reactions/geochemistry/Carbonate.hpp‎

Lines changed: 18 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -101,17 +101,31 @@ constexpr CArrayWrapper<double, 11> reverseRates =
101101
2.51E+06, // CaCl2 = Ca+2 + 2Cl-
102102
2.69E+07, // MgSO4 = Mg+2 + SO4-2
103103
6.62E+07, // NaSO4- = Na+ + SO4-2
104-
8.55E-03 // CaCO3 + H+ = Ca+2 + HCO3-
104+
8.55E-03 // CaCO3 + H+ = Ca+2 + HCO3-
105105
// 1 // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic)
106106
};
107+
108+
constexpr CArrayWrapper<int, 11> mobileSpeciesFlag =
109+
{ 1, // OH- + H+ = H2O
110+
1, // CO2 + H2O = H+ + HCO3-
111+
1, // CO3-2 + H+ = HCO3-
112+
1, // H2CO3 = H+ + HCO3-
113+
1, // CaHCO3+ = Ca+2 + HCO3-
114+
1, // CaSO4 = Ca+2 + SO4-2
115+
1, // CaCl+ = Ca+2 + Cl-
116+
1, // CaCl2 = Ca+2 + 2Cl-
117+
1, // MgSO4 = Mg+2 + SO4-2
118+
1, // NaSO4- = Na+ + SO4-2
119+
1 // CaCO3 + H+ = Ca+2 + HCO3-
120+
};
107121
}
108122
using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 0 >;
109123
using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 11 >;
110124
using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 11, 10 >;
111125

112-
constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
113-
constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
114-
constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
126+
constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
127+
constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
128+
constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
115129

116130
// *****UNCRUSTIFY-ON******
117131
} // namespace geochemistry

‎src/reactions/geochemistry/Forge.hpp‎

Lines changed: 24 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -119,12 +119,35 @@ constexpr CArrayWrapper< double, 19 > reverseRateConstant =
119119
1.0 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq)
120120
};
121121

122+
constexpr CArrayWrapper< int, 19 > mobileSpeciesFlag =
123+
{
124+
1, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻
125+
1, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻
126+
1, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻
127+
1, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻
128+
1, // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source)
129+
1, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻
130+
1, // MgCO₃(aq) + H⁺ ⇌ Mg²⁺ + HCO₃⁻
131+
1, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻
132+
1, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻
133+
1, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻
134+
1, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻
135+
1, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq)
136+
1, // NaHSilO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq)
137+
1, // NaCl ⇌ Na⁺ + Cl⁻
138+
1, // KCl ⇌ K⁺ + Cl⁻
139+
1, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻
140+
1, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻
141+
1, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq)
142+
1 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq)
143+
};
144+
122145
}
123146

124147
using forgeSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 26, 19, 16 >;
125148

126149

127-
constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant );
150+
constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant, forge::mobileSpeciesFlag );
128151

129152

130153
// *****UNCRUSTIFY-ON******

‎src/reactions/geochemistry/Ultramafics.hpp‎

Lines changed: 29 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -126,15 +126,40 @@ constexpr CArrayWrapper<double, 21> reverseRates =
126126
2.83E-44, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O
127127
2.10E-25 // Mg(OH)2 + 2H+ = Mg++ + 2H2O
128128
};
129-
};
129+
130+
constexpr CArrayWrapper<int, 21> mobileSpeciesFlag =
131+
{
132+
1, // OH- + H+ = H2O
133+
1, // CO2(aq) + H2O = HCO3- + H+
134+
1, // CO3-- + H+ = HCO3-
135+
1, // Mg2OH+++ + H+ = 2Mg++ + H2O
136+
1, // Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O
137+
1, // MgOH+ + H+ = Mg++ + H2O
138+
1, // Mg2CO3++ + H+ = 2Mg++ + HCO3-
139+
1, // MgCO3 + H+ = Mg++ + HCO3-
140+
1, // MgHCO3+ = Mg++ + HCO3-
141+
1, // Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O
142+
1, // MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O
143+
1, // MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O
144+
1, // H2SiO4-- + 2H+ = SiO2(aq) + 2H2O
145+
1, // H3SiO4- + H+ = SiO2(aq) + 2H2O
146+
1, // H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O
147+
1, // H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O
148+
1, // Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O
149+
1, // MgCO3 + H+ = Mg++ + HCO3-
150+
1, // SiO2 = SiO2(aq)
151+
1, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O
152+
1 // Mg(OH)2 + 2H+ = Mg++ + 2H2O
153+
};
154+
}
130155

131156
using ultramaficSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 0 >;
132157
using ultramaficSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 21 >;
133158
using ultramaficSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 16 >;
134159

135-
constexpr ultramaficSystemAllKineticType ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
136-
constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
137-
constexpr ultramaficSystemType ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
160+
constexpr ultramaficSystemAllKineticType ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
161+
constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
162+
constexpr ultramaficSystemType ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
138163

139164
// *****UNCRUSTIFY-ON******
140165
} // namespace geochemistry

‎src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -133,10 +133,10 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 )
133133
};
134134

135135
double logPrimarySpeciesConcentration[numPrimarySpecies];
136-
EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
137-
hpcReact::geochemistry::carbonateSystemAllEquilibrium.equilibriumReactionsParameters(),
138-
logInitialPrimarySpeciesConcentration,
139-
logPrimarySpeciesConcentration );
136+
EquilibriumReactionsType::enforceEquilibrium_LogAggregate( 0,
137+
hpcReact::geochemistry::carbonateSystemAllEquilibrium.equilibriumReactionsParameters(),
138+
logInitialPrimarySpeciesConcentration,
139+
logPrimarySpeciesConcentration );
140140

141141
double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
142142
{

0 commit comments

Comments
 (0)