Skip to content

Commit b172c96

Browse files
jkgollajkgolla
committed
Implement unit test for MoMaS benchmark easy case
1 parent 952ba2e commit b172c96

File tree

3 files changed

+189
-0
lines changed

3 files changed

+189
-0
lines changed

‎src/reactions/geochemistry/Momas.hpp‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
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+
#pragma once
13+
14+
#include "../reactionsSystems/Parameters.hpp"
15+
16+
namespace hpcReact
17+
{
18+
19+
namespace geochemistry
20+
{
21+
// turn off uncrustify to allow for better readability of the parameters
22+
// *****UNCRUSTIFY-OFF******
23+
24+
namespace momas
25+
{
26+
27+
constexpr CArrayWrapper<double, 7, 12> stoichMatrix =
28+
{
29+
// C1 C2 C3 C4 C5 CS1 CS2 | X1 X2 X3 X4 S
30+
{ -1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0 }, // C1 = -X2
31+
{ 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0 }, // C2 = X2 + X3
32+
{ 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 1, 0 }, // C3 = X2 + X4
33+
{ 0, 0, 0, -1, 0, 0, 0, 0, -4, 1, 3, 0 }, // C4 = -4X2 + X3 + 3X4
34+
{ 0, 0, 0, 0, -1, 0, 0, 0, 4, 3, 1, 0 }, // C5 = 4X2 + 3X3 + X4
35+
{ 0, 0, 0, 0, 0, -1, 0, 0, 3, 1, 0, 1 }, // CS1 = 3X2 + X3 + S
36+
{ 0, 0, 0, 0, 0, 0, -1, 0, -3, 0, 1, 2 } // CS2 = -3X2 + X4 + 2S
37+
};
38+
39+
constexpr CArrayWrapper<double, 7> equilibriumConstants =
40+
{
41+
1.0e-12, // C1 = -X2
42+
1.0, // C2 = X2 + X3
43+
1.0, // C3 = X2 + X4
44+
0.1, // C4 = -4X2 + X3 + 3X4
45+
1.0e35, // C5 = 4X2 + 3X3 + X4
46+
1.0e6, // CS1 = 3X2 + X3 + S
47+
1.0e-1 // CS2 = -3X2 + X4 + 2S
48+
};
49+
50+
constexpr CArrayWrapper<double, 7> forwardRates =
51+
{
52+
0.0, // C1 = -X2
53+
0.0, // C2 = X2 + X3
54+
0.0, // C3 = X2 + X4
55+
0.0, // C4 = -4X2 + X3 + 3X4
56+
0.0, // C5 = 4X2 + 3X3 + X4
57+
0.0, // CS1 = 3X2 + X3 + S
58+
0.0 // CS2 = -3X2 + X4 + 2S
59+
};
60+
61+
constexpr CArrayWrapper<double, 7> reverseRates =
62+
{
63+
0.0, // C1 = -X2
64+
0.0, // C2 = X2 + X3
65+
0.0, // C3 = X2 + X4
66+
0.0, // C4 = -4X2 + X3 + 3X4
67+
0.0, // C5 = 4X2 + 3X3 + X4
68+
0.0, // CS1 = 3X2 + X3 + S
69+
0.0 // CS2 = -3X2 + X4 + 2S
70+
};
71+
72+
constexpr CArrayWrapper<int, 7> mobileSpeciesFlag =
73+
{
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+
};
82+
}
83+
using momasSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 0 >;
84+
using momasSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >;
85+
using momasSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >;
86+
87+
constexpr momasSystemAllKineticType momasSystemAllKinetic( momas::stoichMatrix, momas::equilibriumConstants, momas::forwardRates, momas::reverseRates, momas::mobileSpeciesFlag );
88+
constexpr momasSystemAllEquilibriumType momasSystemAllEquilibrium( momas::stoichMatrix, momas::equilibriumConstants, momas::forwardRates, momas::reverseRates, momas::mobileSpeciesFlag );
89+
constexpr momasSystemType carbonateSystem( momas::stoichMatrix, momas::equilibriumConstants, momas::forwardRates, momas::reverseRates, momas::mobileSpeciesFlag );
90+
91+
// *****UNCRUSTIFY-ON******
92+
} // namespace geochemistry
93+
} // namespace hpcReact

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

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

89
set( dependencyList hpcReact gtest )

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

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
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 "../Momas.hpp"
15+
16+
using namespace hpcReact;
17+
using namespace hpcReact::geochemistry;
18+
using namespace hpcReact::unitTest_utilities;
19+
20+
21+
// //******************************************************************************
22+
// TEST( testEquilibriumReactions, testEnforceEquilibrium )
23+
// {
24+
// double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
25+
// double const expectedSpeciesConcentrations[5] = { 3.92138294e-01, 3.03930853e-01, 5.05945481e-01, 7.02014628e-01, 5.95970745e-01 };
26+
27+
28+
// std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
29+
// testEnforceEquilibrium< double, 2 >( simpleTestRateParams.equilibriumReactionsParameters(),
30+
// initialSpeciesConcentration,
31+
// expectedSpeciesConcentrations );
32+
33+
// }
34+
35+
36+
//******************************************************************************
37+
38+
TEST( testEquilibriumReactions, testMoMasAllEquilibrium )
39+
{
40+
41+
using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double,
42+
int,
43+
int >;
44+
45+
constexpr int numPrimarySpecies = hpcReact::geochemistry::momasSystemAllEquilibrium.numPrimarySpecies();
46+
47+
double const initialPrimarySpeciesConcentration[numPrimarySpecies] =
48+
{
49+
1.00e-20, // X1
50+
2.00, // X2
51+
1.00e-20, // X3
52+
2.00, // X4
53+
1.00 // S
54+
};
55+
56+
57+
58+
double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] =
59+
{
60+
log( initialPrimarySpeciesConcentration[0] ),
61+
log( initialPrimarySpeciesConcentration[1] ),
62+
log( initialPrimarySpeciesConcentration[2] ),
63+
log( initialPrimarySpeciesConcentration[3] ),
64+
log( initialPrimarySpeciesConcentration[4] )
65+
};
66+
67+
double logPrimarySpeciesConcentration[numPrimarySpecies];
68+
EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
69+
hpcReact::geochemistry::momasSystemAllEquilibrium.equilibriumReactionsParameters(),
70+
logInitialPrimarySpeciesConcentration,
71+
logPrimarySpeciesConcentration );
72+
73+
double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
74+
{
75+
9.7051090442170804E-21, // X1
76+
5.0023298955833342E-12, // X2
77+
1.9327372426296357E-33, // X3
78+
7.3929274619958745E-12, // X4
79+
9.8708294125907346E-13 // S
80+
};
81+
82+
for( int r=0; r<numPrimarySpecies; ++r )
83+
{
84+
EXPECT_NEAR( exp( logPrimarySpeciesConcentration[r] ), expectedPrimarySpeciesConcentrations[r], 1.0e-8 * expectedPrimarySpeciesConcentrations[r] );
85+
}
86+
87+
88+
}
89+
90+
int main( int argc, char * * argv )
91+
{
92+
::testing::InitGoogleTest( &argc, argv );
93+
int const result = RUN_ALL_TESTS();
94+
return result;
95+
}

0 commit comments

Comments
 (0)