Skip to content

MoMaS Medium Test Case with a Kinetic Reaction #12

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Oct 23, 2025
Merged

MoMaS Medium Test Case with a Kinetic Reaction #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
3c38401
cleaned up momas and added medium test case
frankfeifan Sep 10, 2025
2aeda9c
corrections
frankfeifan Sep 28, 2025
7192d47
corrected comments
frankfeifan Oct 1, 2025
4abdabb
Merge remote-tracking branch 'origin/main' into feat/frankfei/momasMe…
frankfeifan Oct 22, 2025
96bcdde
removed ad-hoc parts for momasMedium
frankfeifan Oct 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
94 changes: 85 additions & 9 deletions src/reactions/exampleSystems/MoMasBenchmark.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,13 +15,14 @@

namespace hpcReact
{
namespace MomMasBenchmark
namespace MoMasBenchmark
{
// *****UNCRUSTIFY-OFF******

using simpleSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >;
using easyCaseType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >;
using mediumCaseType = reactionsSystems::MixedReactionsParameters< double, int, int, 14, 10, 9 >;

constexpr simpleSystemType simpleSystemParams =
constexpr easyCaseType easyCaseParams =
{
// Stoichiometric matrix [7 rows × 12 columns]
// Columns 0–6 are secondary species (must be -1 on the diagonal)
Expand All @@ -30,7 +31,7 @@ namespace MomMasBenchmark
// 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
{ 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
Expand All @@ -39,9 +40,9 @@ namespace MomMasBenchmark

// Equilibrium constants K
{
1.0e12, // C1 + X2 = ??
1.0e12, // C1 + X2 = inf
1.0, // C2 = X2 + X3
1.0, // C3 = X2 + X4
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
Expand All @@ -52,7 +53,7 @@ namespace MomMasBenchmark
{
0.0, // C1 = -X2
0.0, // C2 = X2 + X3
0.0, // C3 = X2 + X4
0.0, // C3 = -X2 + X4
0.0, // C4 = -4X2 + X3 + 3X4
0.0, // C5 = 4X2 + 3X3 + X4
0.0, // CS1 = 3X2 + X3 + S
Expand All @@ -63,7 +64,7 @@ namespace MomMasBenchmark
{
0.0, // C1 = -X2
0.0, // C2 = X2 + X3
0.0, // C3 = X2 + X4
0.0, // C3 = -X2 + X4
0.0, // C4 = -4X2 + X3 + 3X4
0.0, // C5 = 4X2 + 3X3 + X4
0.0, // CS1 = 3X2 + X3 + S
Expand All @@ -81,6 +82,81 @@ namespace MomMasBenchmark
}
};

constexpr mediumCaseType mediumCaseParams =
{
// Stoichiometric matrix [10 rows × 14 columns]
// Columns 0–8 are secondary species (must be -1 on the diagonal)
// Columns 9–13 are primary species: {X1, X2, X3, X4, S}
{
// C1 C2 C3 C4 C5 C6 C7 CS1 CS2 | X1 X2 X3 X4 S
{ -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0 }, // C1 = -X2
{ 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0 }, // C2 = X2 + X3
{ 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0 }, // C3 = -X2 + X4
{ 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -4, 1, 3, 0 }, // C4 = -4X2 + X3 + 3X4
{ 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 4, 3, 1, 0 }, // C5 = 4X2 + 3X3 + X4
{ 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 10, 3, 0, 0 }, // C6 = 10X2 + 3X3
{ 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -8, 0, 2, 0 }, // C7 = -8X2 + 2X4
{ 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 3, 1, 0, 1 }, // CS1 = 3X2 + X3 + S
{ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -3, 0, 1, 2 }, // CS2 = -3X2 + X4 + 2S
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 1, 0 }, // Cc = -3X2 + X4 (kinetic)
},

// Equilibrium constants K
{
1.0e12, // C1 + X2 = inf
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-32, // C6 = 10X2 + 3X3
1.0e4, // C7 + 8X2 = 2X4
1.0e-6, // CS1 = 3X2 + X3 + S
1.0e1, // CS2 + 3X2 = X4 + 2S
5 // Cc + 3X2 = X4 (kinetic)
},

// Forward rate constants
{
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, // C6 = 10X2 + 3X3
0.0, // C7 = -8X2 + 2X4
0.0, // CS1 = 3X2 + X3 + S
0.0, // CS2 = -3X2 + X4 + 2S
10.0 // Cc = -3X2 + X4 (kinetic)
},

// Reverse rate constants
{
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, // C6 = 10X2 + 3X3
0.0, // C7 = -8X2 + 2X4
0.0, // CS1 = 3X2 + X3 + S
0.0, // CS2 = -3X2 + X4 + 2S
2.0 // Cc = -3X2 + X4 (kinetic)
},

// 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
1, // C6 = 10X2 + 3X3
1, // C7 = -8X2 + 2X4
0, // CS1 = 3X2 + X3 + S
0, // CS2 = -3X2 + X4 + 2S
1 // Cc = -3X2 + X4 (kinetic)
}
};

// *****UNCRUSTIFY-ON******
} // namespace MomMasBenchmark
} // namespace MoMasBenchmark
} // namespace hpcReact
1 change: 1 addition & 0 deletions src/reactions/exampleSystems/unitTests/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ set( testSourceFiles
testEquilibriumReactions.cpp
testKineticReactions.cpp
testMomasEasyCase.cpp
testMomasMediumCase.cpp
)

set( dependencyList hpcReact gtest )
Expand Down
6 changes: 3 additions & 3 deletions src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
#include "../MoMasBenchmark.hpp"

using namespace hpcReact;
using namespace hpcReact::MomMasBenchmark;
using namespace hpcReact::MoMasBenchmark;
using namespace hpcReact::unitTest_utilities;

//******************************************************************************
Expand All @@ -26,7 +26,7 @@ TEST( testEquilibriumReactions, testMoMasAllEquilibrium )
int,
int >;

constexpr int numPrimarySpecies = hpcReact::MomMasBenchmark::simpleSystemParams.numPrimarySpecies();
constexpr int numPrimarySpecies = hpcReact::MoMasBenchmark::easyCaseParams.numPrimarySpecies();

double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] =
{
Expand Down Expand Up @@ -57,7 +57,7 @@ TEST( testEquilibriumReactions, testMoMasAllEquilibrium )

double logPrimarySpeciesConcentration[numPrimarySpecies];
EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0,
hpcReact::MomMasBenchmark::simpleSystemParams.equilibriumReactionsParameters(),
hpcReact::MoMasBenchmark::easyCaseParams.equilibriumReactionsParameters(),
targetAggregatePrimarySpeciesConcentration,
logInitialPrimarySpeciesConcentration,
logPrimarySpeciesConcentration );
Expand Down
86 changes: 86 additions & 0 deletions src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,86 @@
/*
* ------------------------------------------------------------------------------------------------------------
* 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::MoMasBenchmark;
using namespace hpcReact::unitTest_utilities;

//******************************************************************************

TEST( testEquilibriumReactions, testMoMasMediumEquilibrium )
{
using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double,
int,
int >;

constexpr int numPrimarySpecies = hpcReact::MoMasBenchmark::mediumCaseParams.numPrimarySpecies();

double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] =
{
1.0e-20, // X1
-3.0, // X2
1.0e-20, // X3
1.0, // X4
1.0 // S
};

double const initialPrimarySpeciesConcentration[numPrimarySpecies] =
{
1.0e-20, // X1
0.02, // X2
1.0e-20, // X3
1.0, // X4
1.0 // 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::MoMasBenchmark::mediumCaseParams.equilibriumReactionsParameters(),
targetAggregatePrimarySpeciesConcentration,
logInitialPrimarySpeciesConcentration,
logPrimarySpeciesConcentration );

double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
{
9.9999999999999919e-21, // X1
0.14796989521717838, // X2
5.7165444793692536e-24, // X3
0.025616412699749774, // X4
0.53958559521499294 // 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;
}
1 change: 1 addition & 0 deletions src/reactions/reactionsSystems/KineticReactions_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,6 +271,7 @@ KineticReactions< REAL_TYPE,
RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], s_ri ) : 0.0;
quotient *= productTerm_i;
}

if constexpr( CALCULATE_DERIVATIVES )
{
for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
Expand Down