Skip to content

Commit 61a8c9f

Browse files
frankfeifanfrankfeifan
authored
MoMaS Medium Test Case with a Kinetic Reaction (#12)
* cleaned up momas and added medium test case * removed ad-hoc parts for momasMedium
1 parent 2f94ce2 commit 61a8c9f

File tree

5 files changed

+176
-12
lines changed

5 files changed

+176
-12
lines changed

‎src/reactions/exampleSystems/MoMasBenchmark.hpp‎

Lines changed: 85 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -15,13 +15,14 @@
1515

1616
namespace hpcReact
1717
{
18-
namespace MomMasBenchmark
18+
namespace MoMasBenchmark
1919
{
2020
// *****UNCRUSTIFY-OFF******
2121

22-
using simpleSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >;
22+
using easyCaseType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >;
23+
using mediumCaseType = reactionsSystems::MixedReactionsParameters< double, int, int, 14, 10, 9 >;
2324

24-
constexpr simpleSystemType simpleSystemParams =
25+
constexpr easyCaseType easyCaseParams =
2526
{
2627
// Stoichiometric matrix [7 rows × 12 columns]
2728
// Columns 0–6 are secondary species (must be -1 on the diagonal)
@@ -30,7 +31,7 @@ namespace MomMasBenchmark
3031
// C1 C2 C3 C4 C5 CS1 CS2 | X1 X2 X3 X4 S
3132
{ -1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0 }, // C1 = -X2
3233
{ 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
34+
{ 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 1, 0 }, // C3 = -X2 + X4
3435
{ 0, 0, 0, -1, 0, 0, 0, 0, -4, 1, 3, 0 }, // C4 = -4X2 + X3 + 3X4
3536
{ 0, 0, 0, 0, -1, 0, 0, 0, 4, 3, 1, 0 }, // C5 = 4X2 + 3X3 + X4
3637
{ 0, 0, 0, 0, 0, -1, 0, 0, 3, 1, 0, 1 }, // CS1 = 3X2 + X3 + S
@@ -39,9 +40,9 @@ namespace MomMasBenchmark
3940

4041
// Equilibrium constants K
4142
{
42-
1.0e12, // C1 + X2 = ??
43+
1.0e12, // C1 + X2 = inf
4344
1.0, // C2 = X2 + X3
44-
1.0, // C3 = X2 + X4
45+
1.0, // C3 = -X2 + X4
4546
1.0e1, // C4 + 4X2 = X3 + 3X4
4647
1.0e-35, // C5 = 4X2 + 3X3 + X4
4748
1.0e-6, // CS1 = 3X2 + X3 + S
@@ -52,7 +53,7 @@ namespace MomMasBenchmark
5253
{
5354
0.0, // C1 = -X2
5455
0.0, // C2 = X2 + X3
55-
0.0, // C3 = X2 + X4
56+
0.0, // C3 = -X2 + X4
5657
0.0, // C4 = -4X2 + X3 + 3X4
5758
0.0, // C5 = 4X2 + 3X3 + X4
5859
0.0, // CS1 = 3X2 + X3 + S
@@ -63,7 +64,7 @@ namespace MomMasBenchmark
6364
{
6465
0.0, // C1 = -X2
6566
0.0, // C2 = X2 + X3
66-
0.0, // C3 = X2 + X4
67+
0.0, // C3 = -X2 + X4
6768
0.0, // C4 = -4X2 + X3 + 3X4
6869
0.0, // C5 = 4X2 + 3X3 + X4
6970
0.0, // CS1 = 3X2 + X3 + S
@@ -81,6 +82,81 @@ namespace MomMasBenchmark
8182
}
8283
};
8384

85+
constexpr mediumCaseType mediumCaseParams =
86+
{
87+
// Stoichiometric matrix [10 rows × 14 columns]
88+
// Columns 0–8 are secondary species (must be -1 on the diagonal)
89+
// Columns 9–13 are primary species: {X1, X2, X3, X4, S}
90+
{
91+
// C1 C2 C3 C4 C5 C6 C7 CS1 CS2 | X1 X2 X3 X4 S
92+
{ -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0 }, // C1 = -X2
93+
{ 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0 }, // C2 = X2 + X3
94+
{ 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0 }, // C3 = -X2 + X4
95+
{ 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -4, 1, 3, 0 }, // C4 = -4X2 + X3 + 3X4
96+
{ 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 4, 3, 1, 0 }, // C5 = 4X2 + 3X3 + X4
97+
{ 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 10, 3, 0, 0 }, // C6 = 10X2 + 3X3
98+
{ 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -8, 0, 2, 0 }, // C7 = -8X2 + 2X4
99+
{ 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 3, 1, 0, 1 }, // CS1 = 3X2 + X3 + S
100+
{ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -3, 0, 1, 2 }, // CS2 = -3X2 + X4 + 2S
101+
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 1, 0 }, // Cc = -3X2 + X4 (kinetic)
102+
},
103+
104+
// Equilibrium constants K
105+
{
106+
1.0e12, // C1 + X2 = inf
107+
1.0, // C2 = X2 + X3
108+
1.0, // C3 = -X2 + X4
109+
1.0e1, // C4 + 4X2 = X3 + 3X4
110+
1.0e-35, // C5 = 4X2 + 3X3 + X4
111+
1.0e-32, // C6 = 10X2 + 3X3
112+
1.0e4, // C7 + 8X2 = 2X4
113+
1.0e-6, // CS1 = 3X2 + X3 + S
114+
1.0e1, // CS2 + 3X2 = X4 + 2S
115+
5 // Cc + 3X2 = X4 (kinetic)
116+
},
117+
118+
// Forward rate constants
119+
{
120+
0.0, // C1 = -X2
121+
0.0, // C2 = X2 + X3
122+
0.0, // C3 = -X2 + X4
123+
0.0, // C4 = -4X2 + X3 + 3X4
124+
0.0, // C5 = 4X2 + 3X3 + X4
125+
0.0, // C6 = 10X2 + 3X3
126+
0.0, // C7 = -8X2 + 2X4
127+
0.0, // CS1 = 3X2 + X3 + S
128+
0.0, // CS2 = -3X2 + X4 + 2S
129+
10.0 // Cc = -3X2 + X4 (kinetic)
130+
},
131+
132+
// Reverse rate constants
133+
{
134+
0.0, // C1 = -X2
135+
0.0, // C2 = X2 + X3
136+
0.0, // C3 = -X2 + X4
137+
0.0, // C4 = -4X2 + X3 + 3X4
138+
0.0, // C5 = 4X2 + 3X3 + X4
139+
0.0, // C6 = 10X2 + 3X3
140+
0.0, // C7 = -8X2 + 2X4
141+
0.0, // CS1 = 3X2 + X3 + S
142+
0.0, // CS2 = -3X2 + X4 + 2S
143+
2.0 // Cc = -3X2 + X4 (kinetic)
144+
},
145+
146+
// Flag of mobile secondary species
147+
{ 1, // C1 = -X2
148+
1, // C2 = X2 + X3
149+
1, // C3 = -X2 + X4
150+
1, // C4 = -4X2 + X3 + 3X4
151+
1, // C5 = 4X2 + 3X3 + X4
152+
1, // C6 = 10X2 + 3X3
153+
1, // C7 = -8X2 + 2X4
154+
0, // CS1 = 3X2 + X3 + S
155+
0, // CS2 = -3X2 + X4 + 2S
156+
1 // Cc = -3X2 + X4 (kinetic)
157+
}
158+
};
159+
84160
// *****UNCRUSTIFY-ON******
85-
} // namespace MomMasBenchmark
161+
} // namespace MoMasBenchmark
86162
} // namespace hpcReact

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ set( testSourceFiles
33
testEquilibriumReactions.cpp
44
testKineticReactions.cpp
55
testMomasEasyCase.cpp
6+
testMomasMediumCase.cpp
67
)
78

89
set( dependencyList hpcReact gtest )

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

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@
1414
#include "../MoMasBenchmark.hpp"
1515

1616
using namespace hpcReact;
17-
using namespace hpcReact::MomMasBenchmark;
17+
using namespace hpcReact::MoMasBenchmark;
1818
using namespace hpcReact::unitTest_utilities;
1919

2020
//******************************************************************************
@@ -26,7 +26,7 @@ TEST( testEquilibriumReactions, testMoMasAllEquilibrium )
2626
int,
2727
int >;
2828

29-
constexpr int numPrimarySpecies = hpcReact::MomMasBenchmark::simpleSystemParams.numPrimarySpecies();
29+
constexpr int numPrimarySpecies = hpcReact::MoMasBenchmark::easyCaseParams.numPrimarySpecies();
3030

3131
double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] =
3232
{
@@ -57,7 +57,7 @@ TEST( testEquilibriumReactions, testMoMasAllEquilibrium )
5757

5858
double logPrimarySpeciesConcentration[numPrimarySpecies];
5959
EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
60-
hpcReact::MomMasBenchmark::simpleSystemParams.equilibriumReactionsParameters(),
60+
hpcReact::MoMasBenchmark::easyCaseParams.equilibriumReactionsParameters(),
6161
targetAggregatePrimarySpeciesConcentration,
6262
logInitialPrimarySpeciesConcentration,
6363
logPrimarySpeciesConcentration );

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

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
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::MoMasBenchmark;
18+
using namespace hpcReact::unitTest_utilities;
19+
20+
//******************************************************************************
21+
22+
TEST( testEquilibriumReactions, testMoMasMediumEquilibrium )
23+
{
24+
using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double,
25+
int,
26+
int >;
27+
28+
constexpr int numPrimarySpecies = hpcReact::MoMasBenchmark::mediumCaseParams.numPrimarySpecies();
29+
30+
double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] =
31+
{
32+
1.0e-20, // X1
33+
-3.0, // X2
34+
1.0e-20, // X3
35+
1.0, // X4
36+
1.0 // S
37+
};
38+
39+
double const initialPrimarySpeciesConcentration[numPrimarySpecies] =
40+
{
41+
1.0e-20, // X1
42+
0.02, // X2
43+
1.0e-20, // X3
44+
1.0, // X4
45+
1.0 // S
46+
};
47+
48+
double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
49+
{
50+
log( initialPrimarySpeciesConcentration[0] ),
51+
log( initialPrimarySpeciesConcentration[1] ),
52+
log( initialPrimarySpeciesConcentration[2] ),
53+
log( initialPrimarySpeciesConcentration[3] ),
54+
log( initialPrimarySpeciesConcentration[4] )
55+
};
56+
57+
double logPrimarySpeciesConcentration[numPrimarySpecies];
58+
EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
59+
hpcReact::MoMasBenchmark::mediumCaseParams.equilibriumReactionsParameters(),
60+
targetAggregatePrimarySpeciesConcentration,
61+
logInitialPrimarySpeciesConcentration,
62+
logPrimarySpeciesConcentration );
63+
64+
double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
65+
{
66+
9.9999999999999919e-21, // X1
67+
0.14796989521717838, // X2
68+
5.7165444793692536e-24, // X3
69+
0.025616412699749774, // X4
70+
0.53958559521499294 // S
71+
};
72+
73+
for( int r=0; r<numPrimarySpecies; ++r )
74+
{
75+
EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
76+
}
77+
78+
79+
}
80+
81+
int main( int argc, char * * argv )
82+
{
83+
::testing::InitGoogleTest( &argc, argv );
84+
int const result = RUN_ALL_TESTS();
85+
return result;
86+
}

‎src/reactions/reactionsSystems/KineticReactions_impl.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -271,6 +271,7 @@ KineticReactions< REAL_TYPE,
271271
RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], s_ri ) : 0.0;
272272
quotient *= productTerm_i;
273273
}
274+
274275
if constexpr( CALCULATE_DERIVATIVES )
275276
{
276277
for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )

0 commit comments

Comments
 (0)