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
rrsettgast 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 3 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
0.01 // 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;
}
8 changes: 7 additions & 1 deletion src/reactions/reactionsSystems/KineticReactions_impl.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -236,8 +236,9 @@ KineticReactions< REAL_TYPE,
}

// get/calculate the forward and reverse rate constants for this reaction
RealType const rateConstant = params.rateConstantForward( r ); //* exp( -params.m_activationEnergy[r] / ( constants::R *
// RealType const rateConstant = params.rateConstantForward( r ); //* exp( -params.m_activationEnergy[r] / ( constants::R *
// temperature ) );
RealType rateConstant;
RealType const equilibriumConstant = params.equilibriumConstant( r );

RealType quotient = 1.0;
Expand All @@ -253,6 +254,8 @@ KineticReactions< REAL_TYPE,
}
quotient = exp( logQuotient );

rateConstant = quotient < equilibriumConstant ? params.rateConstantForward( r ) : params.rateConstantReverse( r );

if constexpr( CALCULATE_DERIVATIVES )
{
for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
Expand All @@ -271,6 +274,9 @@ KineticReactions< REAL_TYPE,
RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], s_ri ) : 0.0;
quotient *= productTerm_i;
}

rateConstant = quotient < equilibriumConstant ? params.rateConstantForward( r ) : params.rateConstantReverse( r );

if constexpr( CALCULATE_DERIVATIVES )
{
for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i )
Expand Down
3 changes: 2 additions & 1 deletion src/reactions/reactionsSystems/Parameters.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,8 @@ struct KineticReactionsParameters
RealType stoichiometricMatrix( IndexType const r, int const i ) const { return m_stoichiometricMatrix[r][i]; }
RealType rateConstantForward( IndexType const r ) const { return m_rateConstantForward[r]; }
RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; }
RealType equilibriumConstant( IndexType const r ) const { return m_rateConstantForward[r] / m_rateConstantReverse[r]; }
// RealType equilibriumConstant( IndexType const r ) const { return m_rateConstantForward[r] / m_rateConstantReverse[r]; }
RealType equilibriumConstant( IndexType const r ) const { return m_equilibiriumConstant[r]; } // Just for MoMaS medium now


CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix;
Expand Down