diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 4eae0d3..c0c8f06 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/reactions/exampleSystems/unitTests/CMakeLists.txt b/src/reactions/exampleSystems/unitTests/CMakeLists.txt index a5174b1..710de67 100644 --- a/src/reactions/exampleSystems/unitTests/CMakeLists.txt +++ b/src/reactions/exampleSystems/unitTests/CMakeLists.txt @@ -3,6 +3,7 @@ set( testSourceFiles testEquilibriumReactions.cpp testKineticReactions.cpp testMomasEasyCase.cpp + testMomasMediumCase.cpp ) set( dependencyList hpcReact gtest ) diff --git a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp index 2f77886..3fe3eb6 100644 --- a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp +++ b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp @@ -14,7 +14,7 @@ #include "../MoMasBenchmark.hpp" using namespace hpcReact; -using namespace hpcReact::MomMasBenchmark; +using namespace hpcReact::MoMasBenchmark; using namespace hpcReact::unitTest_utilities; //****************************************************************************** @@ -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] = { @@ -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 ); diff --git a/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp new file mode 100644 index 0000000..dffb5da --- /dev/null +++ b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp @@ -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 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 )