From 3c384017ef1f5e4bb2d2b5f10be53feb168736a8 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Wed, 10 Sep 2025 13:19:24 -0700 Subject: [PATCH 01/10] cleaned up momas and added medium test case --- .../exampleSystems/MoMasBenchmark.hpp | 86 +++++++++++++++++-- .../exampleSystems/unitTests/CMakeLists.txt | 1 + .../unitTests/testMomasEasyCase.cpp | 6 +- .../unitTests/testMomasMediumCase.cpp | 86 +++++++++++++++++++ .../KineticReactions_impl.hpp | 8 +- 5 files changed, 178 insertions(+), 9 deletions(-) create mode 100644 src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 4eae0d3..9937658 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) @@ -39,7 +40,7 @@ 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.0e1, // C4 + 4X2 = X3 + 3X4 @@ -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, 2, 0 }, // Cc = -3X2 + 2X4 (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 = 2X4 (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 + 2X4 (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 + 2X4 (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 + 2X4 (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; } + + rateConstant = quotient < equilibriumConstant ? params.rateConstantForward( r ) : params.rateConstantReverse( r ); + if constexpr( CALCULATE_DERIVATIVES ) { for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) From 2aeda9cd1c59e548a9ef8b3c804cfcee2ab75b32 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Sun, 28 Sep 2025 14:28:06 -0700 Subject: [PATCH 02/10] corrections --- src/reactions/exampleSystems/MoMasBenchmark.hpp | 10 +++++----- src/reactions/reactionsSystems/Parameters.hpp | 3 ++- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 9937658..1af6608 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -98,7 +98,7 @@ constexpr mediumCaseType mediumCaseParams = { 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, 2, 0 }, // Cc = -3X2 + 2X4 (kinetic) + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 1, 0 }, // Cc = -3X2 + X4 (kinetic) }, // Equilibrium constants K @@ -112,7 +112,7 @@ constexpr mediumCaseType mediumCaseParams = 1.0e4, // C7 + 8X2 = 2X4 1.0e-6, // CS1 = 3X2 + X3 + S 1.0e1, // CS2 + 3X2 = X4 + 2S - 5, // Cc + 3X2 = 2X4 (kinetic) + 5 // Cc + 3X2 = X4 (kinetic) }, // Forward rate constants @@ -126,7 +126,7 @@ constexpr mediumCaseType mediumCaseParams = 0.0, // C7 = -8X2 + 2X4 0.0, // CS1 = 3X2 + X3 + S 0.0, // CS2 = -3X2 + X4 + 2S - 10.0 // Cc = -3X2 + 2X4 (kinetic) + 10.0 // Cc = -3X2 + X4 (kinetic) }, // Reverse rate constants @@ -140,7 +140,7 @@ constexpr mediumCaseType mediumCaseParams = 0.0, // C7 = -8X2 + 2X4 0.0, // CS1 = 3X2 + X3 + S 0.0, // CS2 = -3X2 + X4 + 2S - 0.01 // Cc = -3X2 + 2X4 (kinetic) + 0.01 // Cc = -3X2 + X4 (kinetic) }, // Flag of mobile secondary species @@ -153,7 +153,7 @@ constexpr mediumCaseType mediumCaseParams = 1, // C7 = -8X2 + 2X4 0, // CS1 = 3X2 + X3 + S 0, // CS2 = -3X2 + X4 + 2S - 1 // Cc = -3X2 + 2X4 (kinetic) + 1 // Cc = -3X2 + X4 (kinetic) } }; diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index d2ca7ba..03e7f1a 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -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; From 7192d4717be45d72beabeb7c972e7dcee190f352 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Wed, 1 Oct 2025 14:16:27 -0700 Subject: [PATCH 03/10] corrected comments --- src/reactions/exampleSystems/MoMasBenchmark.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 1af6608..16072c0 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -31,7 +31,7 @@ namespace MoMasBenchmark // 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 @@ -42,7 +42,7 @@ namespace MoMasBenchmark { 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 @@ -53,7 +53,7 @@ namespace MoMasBenchmark { 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 @@ -64,7 +64,7 @@ namespace MoMasBenchmark { 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 @@ -91,7 +91,7 @@ constexpr mediumCaseType mediumCaseParams = // 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, -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 @@ -105,7 +105,7 @@ constexpr mediumCaseType mediumCaseParams = { 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-32, // C6 = 10X2 + 3X3 @@ -119,7 +119,7 @@ constexpr mediumCaseType mediumCaseParams = { 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, // C6 = 10X2 + 3X3 @@ -133,7 +133,7 @@ constexpr mediumCaseType mediumCaseParams = { 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, // C6 = 10X2 + 3X3 From c95d4caea34ebeecfdf2d3e899aa436cd8c53134 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Wed, 22 Oct 2025 12:22:39 -0700 Subject: [PATCH 04/10] added 3 species full kinetic serial chain reaction --- src/reactions/exampleSystems/ChainGeneric.hpp | 66 +++++++++++++++++++ src/reactions/massActions/MassActions.hpp | 28 +++++--- ...ionsAggregatePrimaryConcentration_impl.hpp | 5 ++ .../MixedEquilibriumKineticReactions_impl.hpp | 40 +++++++---- 4 files changed, 120 insertions(+), 19 deletions(-) create mode 100644 src/reactions/exampleSystems/ChainGeneric.hpp diff --git a/src/reactions/exampleSystems/ChainGeneric.hpp b/src/reactions/exampleSystems/ChainGeneric.hpp new file mode 100644 index 0000000..8e37bf1 --- /dev/null +++ b/src/reactions/exampleSystems/ChainGeneric.hpp @@ -0,0 +1,66 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: (BSD-3-Clause) + * + * Copyright (c) 2025- Lawrence Livermore National Security LLC + * All rights reserved + * + * See top level LICENSE files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#pragma once + +#include "../reactionsSystems/Parameters.hpp" + +namespace hpcReact +{ +namespace ChainGeneric +{ +// *****UNCRUSTIFY-OFF****** + + using serialAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 3, 3, 0 >; + + constexpr serialAllKineticType serialAllKineticParams = +{ + // Stoichiometric matrix [3 rows × 3 columns] + // Columns 0–3 are primary species: {C1, C2, C3 } + { + // C1 C2 C3 + { -1, 1, 0 }, // C1 = C2 + { 0, -1, 1 }, // C2 = C3 + { 0, 0, -1 }, // C3 = + }, + + // Equilibrium constants K + { + 0, // C1 = C2 + 0, // C2 = C3 + 0 // C3 + }, + + // Forward rate constants + { + 0.05, // C1 = C2 + 0.03, // C2 = C3 + 0.02, // C3 + }, + + // Reverse rate constants + { + 0.0, // C1 = C2 + 0.0, // C2 = C3 + 0.0 // C3 + }, + + // Flag of mobile secondary species + { + 1, // C1 = C2 + 1, // C2 = C3 + 1 // C3 + } +}; + +// *****UNCRUSTIFY-ON****** +} // namespace ChainGeneric +} // namespace hpcReact diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index f8d6b96..e73ea6b 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -72,6 +72,11 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, ARRAY_1D & logSecondarySpeciesConcentrations ) { + if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) + { + return; + } + massActions_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, @@ -173,15 +178,22 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, { constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); - REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0}; + if constexpr( numSecondarySpecies > 0 ) + { + REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0}; - calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations, - aggregatePrimarySpeciesConcentrations, - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } + else + { + GEOS_UNUSED_VAR( logPrimarySpeciesConcentrations, aggregatePrimarySpeciesConcentrations, dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } } diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index 07e99f3..5cf2ca5 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -113,6 +113,11 @@ EquilibriumReactions< REAL_TYPE, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, ARRAY_1D & logPrimarySpeciesConcentration ) { + if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) + { + return; + } + HPCREACT_UNUSED_VAR( temperature ); constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp index f6dfd18..e898670 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp @@ -58,16 +58,32 @@ MixedEquilibriumKineticReactions< REAL_TYPE, ARRAY_1D_PRIMARY & aggregateSpeciesRates, ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ) { - // 1. Compute new aggregate species from primary species - massActions::calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params.equilibriumReactionsParameters(), - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations, - aggregatePrimarySpeciesConcentrations, - mobileAggregatePrimarySpeciesConcentrations, - dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, - dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + if constexpr( PARAMS_DATA::numEquilibriumReactions() > 0 ) + { + // 1. Compute new aggregate species from primary species + massActions::calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params.equilibriumReactionsParameters(), + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + mobileAggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, + dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + } + else + { + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); + aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; + mobileAggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; + dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; + dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; + } + } if constexpr( PARAMS_DATA::numKineticReactions() > 0 ) { @@ -119,6 +135,8 @@ MixedEquilibriumKineticReactions< REAL_TYPE, ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations ) { + HPCREACT_UNUSED_VAR( surfaceArea ); + constexpr IntType numSpecies = PARAMS_DATA::numSpecies(); constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); @@ -147,7 +165,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, kineticReactions::computeReactionRates( temperature, params.kineticReactionsParameters(), logSpeciesConcentration, - surfaceArea, + // surfaceArea, reactionRates, reactionRatesDerivatives ); From 96bcddea1f8134aee52b2270fd0ec68262970f86 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Wed, 22 Oct 2025 12:29:44 -0700 Subject: [PATCH 05/10] removed ad-hoc parts for momasMedium --- src/reactions/exampleSystems/MoMasBenchmark.hpp | 2 +- src/reactions/reactionsSystems/KineticReactions_impl.hpp | 7 +------ src/reactions/reactionsSystems/Parameters.hpp | 3 +-- 3 files changed, 3 insertions(+), 9 deletions(-) diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 16072c0..c0c8f06 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -140,7 +140,7 @@ constexpr mediumCaseType mediumCaseParams = 0.0, // C7 = -8X2 + 2X4 0.0, // CS1 = 3X2 + X3 + S 0.0, // CS2 = -3X2 + X4 + 2S - 0.01 // Cc = -3X2 + X4 (kinetic) + 2.0 // Cc = -3X2 + X4 (kinetic) }, // Flag of mobile secondary species diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp index 2c97275..00debf7 100644 --- a/src/reactions/reactionsSystems/KineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/KineticReactions_impl.hpp @@ -236,9 +236,8 @@ 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; @@ -254,8 +253,6 @@ 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 ) @@ -274,8 +271,6 @@ 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 ) { diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index 03e7f1a..d2ca7ba 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -100,8 +100,7 @@ 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_equilibiriumConstant[r]; } // Just for MoMaS medium now + RealType equilibriumConstant( IndexType const r ) const { return m_rateConstantForward[r] / m_rateConstantReverse[r]; } CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; From aa697d02e34528c6c68c81223297eb69d39f9917 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Thu, 23 Oct 2025 17:17:30 -0700 Subject: [PATCH 06/10] added an option variable for user-defined kinetic rate calculation method --- src/reactions/exampleSystems/ChainGeneric.hpp | 68 ++++++++++--------- .../MixedEquilibriumKineticReactions_impl.hpp | 31 ++++++--- src/reactions/reactionsSystems/Parameters.hpp | 9 ++- 3 files changed, 64 insertions(+), 44 deletions(-) diff --git a/src/reactions/exampleSystems/ChainGeneric.hpp b/src/reactions/exampleSystems/ChainGeneric.hpp index 8e37bf1..8b0acd8 100644 --- a/src/reactions/exampleSystems/ChainGeneric.hpp +++ b/src/reactions/exampleSystems/ChainGeneric.hpp @@ -22,44 +22,46 @@ namespace ChainGeneric using serialAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 3, 3, 0 >; constexpr serialAllKineticType serialAllKineticParams = -{ - // Stoichiometric matrix [3 rows × 3 columns] - // Columns 0–3 are primary species: {C1, C2, C3 } { - // C1 C2 C3 - { -1, 1, 0 }, // C1 = C2 - { 0, -1, 1 }, // C2 = C3 - { 0, 0, -1 }, // C3 = - }, + // Stoichiometric matrix [3 rows × 3 columns] + // Columns 0–3 are primary species: {C1, C2, C3 } + { + // C1 C2 C3 + { -1, 1, 0 }, // C1 = C2 + { 0, -1, 1 }, // C2 = C3 + { 0, 0, -1 }, // C3 = + }, - // Equilibrium constants K - { - 0, // C1 = C2 - 0, // C2 = C3 - 0 // C3 - }, + // Equilibrium constants K + { + 0, // C1 = C2 + 0, // C2 = C3 + 0 // C3 + }, + + // Forward rate constants + { + 0.05, // C1 = C2 + 0.03, // C2 = C3 + 0.02, // C3 + }, - // Forward rate constants - { - 0.05, // C1 = C2 - 0.03, // C2 = C3 - 0.02, // C3 - }, + // Reverse rate constants + { + 0.0, // C1 = C2 + 0.0, // C2 = C3 + 0.0 // C3 + }, - // Reverse rate constants - { - 0.0, // C1 = C2 - 0.0, // C2 = C3 - 0.0 // C3 - }, + // Flag of mobile secondary species + { + 1, // C1 = C2 + 1, // C2 = C3 + 1 // C3 + }, - // Flag of mobile secondary species - { - 1, // C1 = C2 - 1, // C2 = C3 - 1 // C3 - } -}; + 0 // Use the forward and reverse to calculate the kinetic reaction rates + }; // *****UNCRUSTIFY-ON****** } // namespace ChainGeneric diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp index e898670..fb13116 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp @@ -134,9 +134,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, ARRAY_1D & reactionRates, ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations ) -{ - HPCREACT_UNUSED_VAR( surfaceArea ); - +{ constexpr IntType numSpecies = PARAMS_DATA::numSpecies(); constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); @@ -162,12 +160,27 @@ MixedEquilibriumKineticReactions< REAL_TYPE, CArrayWrapper< RealType, numKineticReactions, numSpecies > reactionRatesDerivatives; - kineticReactions::computeReactionRates( temperature, - params.kineticReactionsParameters(), - logSpeciesConcentration, - // surfaceArea, - reactionRates, - reactionRatesDerivatives ); + if( params.reactionRatesUpdateOption() == 0 ) + { + kineticReactions::computeReactionRates( temperature, + params.kineticReactionsParameters(), + logSpeciesConcentration, + reactionRates, + reactionRatesDerivatives ); + } + else if( params.reactionRatesUpdateOption() == 1 ) + { + kineticReactions::computeReactionRates( temperature, + params.kineticReactionsParameters(), + logSpeciesConcentration, + surfaceArea, + reactionRates, + reactionRatesDerivatives ); + } + else + { + throw std::runtime_error( "Error: No such option implemented so far for reaction rate calculation."); + } // Compute the reaction rates derivatives w.r.t. log primary species concentrations for( IntType i = 0; i < numKineticReactions; ++i ) diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index d2ca7ba..0956d9b 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -129,12 +129,14 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > const & equilibriumConstant, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse, - CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag ): + CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag, + IntType const reactionRatesUpdateOption = 1 ): m_stoichiometricMatrix( stoichiometricMatrix ), m_equilibriumConstant( equilibriumConstant ), m_rateConstantForward( rateConstantForward ), m_rateConstantReverse( rateConstantReverse ), - m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ) + m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ), + m_reactionRatesUpdateOption( reactionRatesUpdateOption ) {} static constexpr IndexType numReactions() { return NUM_REACTIONS; } @@ -232,12 +234,15 @@ struct MixedReactionsParameters RealType equilibriumConstant( IndexType const r ) const { return m_equilibriumConstant[r]; } RealType rateConstantForward( IndexType const r ) const { return m_rateConstantForward[r]; } RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; } + IntType reactionRatesUpdateOption() const { return m_reactionRatesUpdateOption; } CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; CArrayWrapper< RealType, NUM_REACTIONS > m_equilibriumConstant; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag; + + IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. }; From 74925cf9492876ccbde1794e9ab8d52130509e27 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 24 Oct 2025 18:13:49 -0700 Subject: [PATCH 07/10] modified computeReactionRatesTest and added a new test for chain reaction --- src/reactions/exampleSystems/BulkGeneric.hpp | 8 ++- .../exampleSystems/unitTests/CMakeLists.txt | 1 + .../unitTests/testGenericChainReactions.cpp | 64 +++++++++++++++++++ .../unitTests/testKineticReactions.cpp | 3 + src/reactions/geochemistry/Carbonate.hpp | 2 +- .../testGeochemicalKineticReactions.cpp | 14 ++++ src/reactions/massActions/MassActions.hpp | 2 +- ...ionsAggregatePrimaryConcentration_impl.hpp | 4 +- .../reactionsSystems/KineticReactions.hpp | 29 +++++++-- .../KineticReactions_impl.hpp | 2 +- .../MixedEquilibriumKineticReactions_impl.hpp | 29 ++------- src/reactions/reactionsSystems/Parameters.hpp | 14 ++-- .../kineticReactionsTestUtilities.hpp | 2 + 13 files changed, 134 insertions(+), 40 deletions(-) create mode 100644 src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp index 036a6cd..77547ce 100644 --- a/src/reactions/exampleSystems/BulkGeneric.hpp +++ b/src/reactions/exampleSystems/BulkGeneric.hpp @@ -63,7 +63,9 @@ simpleKineticTestRateParams = // reverse rate constants { 1.0, 0.5 }, // flag of mobile secondary species - { 1, 1 } + { 1, 1 }, + // Use the forward and reverse to calculate the kinetic reaction rates + 0 }; using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >; @@ -84,7 +86,9 @@ simpleTestRateParams = // reverse rate constants { 1.0, 0.5 }, // flag of mobile secondary species - { 1, 1 } + { 1, 1 }, + // Use the forward and reverse to calculate the kinetic reaction rates + 0 }; // *****UNCRUSTIFY-ON****** diff --git a/src/reactions/exampleSystems/unitTests/CMakeLists.txt b/src/reactions/exampleSystems/unitTests/CMakeLists.txt index 710de67..09ea7a0 100644 --- a/src/reactions/exampleSystems/unitTests/CMakeLists.txt +++ b/src/reactions/exampleSystems/unitTests/CMakeLists.txt @@ -4,6 +4,7 @@ set( testSourceFiles testKineticReactions.cpp testMomasEasyCase.cpp testMomasMediumCase.cpp + testGenericChainReactions.cpp ) set( dependencyList hpcReact gtest ) diff --git a/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp new file mode 100644 index 0000000..6e0f9a5 --- /dev/null +++ b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp @@ -0,0 +1,64 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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/kineticReactionsTestUtilities.hpp" +#include "../ChainGeneric.hpp" + + +using namespace hpcReact; +using namespace hpcReact::unitTest_utilities; + +//****************************************************************************** +TEST( testChainGenericKineticReactions, computeReactionRatesTest_chainReactionParams ) +{ + using namespace hpcReact::ChainGeneric; + + double const initialSpeciesConcentration[] = + { + 1.0, // C1 + 1e-8, // C2 + 1e-8 // C3 + }; + double const surfaceArea[] = + { + 0.0, // C1 -> C2 + 0.0, // C2 -> C3 + 0.0 // C3 -> + }; + + double const expectedReactionRates[] = + { + 0.05, // C1 -> C2 + 3e-10, // C2 -> C3 + 2e-10 // C3 -> + }; + double const expectedReactionRatesDerivatives[][3] = + { { 0.05, 0.0, 0.0 }, + { 0.0, 0.03, 0.0 }, + { 0.0, 0.0, 0.02 } }; + computeReactionRatesTest< double, false >( serialAllKineticParams.kineticReactionsParameters(), + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); + computeReactionRatesTest< double, true >( serialAllKineticParams.kineticReactionsParameters(), + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} diff --git a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp index 7baeb0e..1286169 100644 --- a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp +++ b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp @@ -24,16 +24,19 @@ using namespace hpcReact::unitTest_utilities; TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams ) { double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; + double const surfaceArea[] = { 0.0, 0.0 }; double const expectedReactionRates[] = { 1.0, 0.25 }; double const expectedReactionRatesDerivatives[][5] = { { 2.0, -0.5, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.5, 0.25, 0.0 } }; computeReactionRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); } diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp index f36fe0a..52ae3ab 100644 --- a/src/reactions/geochemistry/Carbonate.hpp +++ b/src/reactions/geochemistry/Carbonate.hpp @@ -113,7 +113,7 @@ using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParame using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 10, 10 >; using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 16, 10, 9 >; -constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag ); +constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag, 0 ); constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag ); constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag ); diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp index b84f1bd..dd36fb7 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp @@ -43,6 +43,18 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic ) 1.09 // Na+1 }; + double const surfaceArea[10] = { 0.0, // OH- + H+ = H2O + 0.0, // CO2 + H2O = H+ + HCO3- + 0.0, // CO3-2 + H+ = HCO3- + 0.0, // CaHCO3+ = Ca+2 + HCO3- + 0.0, // CaSO4 = Ca+2 + SO4-2 + 0.0, // CaCl+ = Ca+2 + Cl- + 0.0, // CaCl2 = Ca+2 + 2Cl- + 0.0, // MgSO4 = Mg+2 + SO4-2 + 0.0, // NaSO4- = Na+ + SO4-2 + 0.0, // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + }; + double const expectedReactionRates[10] = { -0.001424736, // OH- + H+ = H2O -12610.7392, // CO2 + H2O = H+ + HCO3- -0.175591624, // CO3-2 + H+ = HCO3- @@ -71,10 +83,12 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic ) computeReactionRatesTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); computeReactionRatesTest< double, true >( carbonateSystemAllKinetic.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); } diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index e73ea6b..a72fe03 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -76,7 +76,7 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, { return; } - + massActions_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index 5cf2ca5..2cd0a13 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -115,9 +115,9 @@ EquilibriumReactions< REAL_TYPE, { if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) { - return; + return; } - + HPCREACT_UNUSED_VAR( temperature ); constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp index d66f1e9..0481e44 100644 --- a/src/reactions/reactionsSystems/KineticReactions.hpp +++ b/src/reactions/reactionsSystems/KineticReactions.hpp @@ -13,6 +13,8 @@ #include "common/macros.hpp" +#include + /** @file KineticReactions.hpp * @brief Header file for the KineticReactions class. * @author HPC-REACT Team @@ -130,12 +132,27 @@ class KineticReactions ARRAY_1D & reactionRates, ARRAY_2D & reactionRatesDerivatives ) { - computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature, - params, - speciesConcentration, - surfaceArea, - reactionRates, - reactionRatesDerivatives ); + if( params.reactionRatesUpdateOption() == 0 ) + { + computeReactionRates_impl< PARAMS_DATA, true >( temperature, + params, + speciesConcentration, + reactionRates, + reactionRatesDerivatives ); + } + else if( params.reactionRatesUpdateOption() == 1 ) + { + computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature, + params, + speciesConcentration, + surfaceArea, + reactionRates, + reactionRatesDerivatives ); + } + else + { + throw std::runtime_error( "Error: No such option implemented so far for reaction rate calculation." ); + } } diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp index 00debf7..59e8aa1 100644 --- a/src/reactions/reactionsSystems/KineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/KineticReactions_impl.hpp @@ -271,7 +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 ) diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp index fb13116..f24e4d6 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp @@ -134,7 +134,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, ARRAY_1D & reactionRates, ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations ) -{ +{ constexpr IntType numSpecies = PARAMS_DATA::numSpecies(); constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); @@ -160,27 +160,12 @@ MixedEquilibriumKineticReactions< REAL_TYPE, CArrayWrapper< RealType, numKineticReactions, numSpecies > reactionRatesDerivatives; - if( params.reactionRatesUpdateOption() == 0 ) - { - kineticReactions::computeReactionRates( temperature, - params.kineticReactionsParameters(), - logSpeciesConcentration, - reactionRates, - reactionRatesDerivatives ); - } - else if( params.reactionRatesUpdateOption() == 1 ) - { - kineticReactions::computeReactionRates( temperature, - params.kineticReactionsParameters(), - logSpeciesConcentration, - surfaceArea, - reactionRates, - reactionRatesDerivatives ); - } - else - { - throw std::runtime_error( "Error: No such option implemented so far for reaction rate calculation."); - } + kineticReactions::computeReactionRates( temperature, + params.kineticReactionsParameters(), + logSpeciesConcentration, + surfaceArea, + reactionRates, + reactionRatesDerivatives ); // Compute the reaction rates derivatives w.r.t. log primary species concentrations for( IntType i = 0; i < numKineticReactions; ++i ) diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index 0956d9b..0e8e95f 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -89,11 +89,13 @@ struct KineticReactionsParameters constexpr KineticReactionsParameters( CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > const & stoichiometricMatrix, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse, - CArrayWrapper< RealType, NUM_REACTIONS > const & equilibriumConstant ): + CArrayWrapper< RealType, NUM_REACTIONS > const & equilibriumConstant, + IntType const reactionRatesUpdateOption ): m_stoichiometricMatrix( stoichiometricMatrix ), m_rateConstantForward( rateConstantForward ), m_rateConstantReverse( rateConstantReverse ), - m_equilibiriumConstant( equilibriumConstant ) // Initialize to empty array + m_equilibiriumConstant( equilibriumConstant ), // Initialize to empty array + m_reactionRatesUpdateOption( reactionRatesUpdateOption ) {} @@ -102,11 +104,14 @@ struct KineticReactionsParameters RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; } RealType equilibriumConstant( IndexType const r ) const { return m_rateConstantForward[r] / m_rateConstantReverse[r]; } + IntType reactionRatesUpdateOption() const { return m_reactionRatesUpdateOption; } CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; CArrayWrapper< RealType, NUM_REACTIONS > m_equilibiriumConstant; + + IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. }; @@ -192,7 +197,7 @@ struct MixedReactionsParameters equilibriumConstant( i ) = m_equilibriumConstant( numEquilibriumReactions() + i ); } - return { kineticMatrix, rateConstantForward, rateConstantReverse, equilibriumConstant }; + return { kineticMatrix, rateConstantForward, rateConstantReverse, equilibriumConstant, m_reactionRatesUpdateOption }; } void verifyParameterConsistency() @@ -234,7 +239,6 @@ struct MixedReactionsParameters RealType equilibriumConstant( IndexType const r ) const { return m_equilibriumConstant[r]; } RealType rateConstantForward( IndexType const r ) const { return m_rateConstantForward[r]; } RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; } - IntType reactionRatesUpdateOption() const { return m_reactionRatesUpdateOption; } CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; CArrayWrapper< RealType, NUM_REACTIONS > m_equilibriumConstant; @@ -242,7 +246,7 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag; - IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. + IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. }; diff --git a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp index 1c650a3..789bd50 100644 --- a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp @@ -36,6 +36,7 @@ template< typename REAL_TYPE, typename PARAMS_DATA > void computeReactionRatesTest( PARAMS_DATA const & params, REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], + REAL_TYPE const (&surfaceArea)[PARAMS_DATA::numReactions()], REAL_TYPE const (&expectedReactionRates)[PARAMS_DATA::numReactions()], REAL_TYPE const (&expectedReactionRatesDerivatives)[PARAMS_DATA::numReactions()][PARAMS_DATA::numSpecies()] ) { @@ -79,6 +80,7 @@ void computeReactionRatesTest( PARAMS_DATA const & params, KineticReactionsType::computeReactionRates( temperature, params, speciesConcentration, + surfaceArea, reactionRates, reactionRatesDerivatives ); From 56c7efc7bd874234b1250ad2dbd64af22378b849 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Mon, 27 Oct 2025 09:55:45 -0700 Subject: [PATCH 08/10] codecov --- src/reactions/reactionsSystems/KineticReactions.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp index 0481e44..a5f5575 100644 --- a/src/reactions/reactionsSystems/KineticReactions.hpp +++ b/src/reactions/reactionsSystems/KineticReactions.hpp @@ -149,10 +149,6 @@ class KineticReactions reactionRates, reactionRatesDerivatives ); } - else - { - throw std::runtime_error( "Error: No such option implemented so far for reaction rate calculation." ); - } } From f811c0a45abd11f798f0315debbc632e3a1a5ee3 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Tue, 28 Oct 2025 11:12:16 -0700 Subject: [PATCH 09/10] added a kinetic reaction test to call quotient reaction rate compute --- .../testGeochemicalKineticReactions.cpp | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp index dd36fb7..d8d2415 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp @@ -93,7 +93,48 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic ) expectedReactionRatesDerivatives ); } +TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem ) +{ + double const initialSpeciesConcentration[16] = + { + 1.0e-16, // OH- + 1.0e-16, // CO2 + 1.0e-16, // CO3-2 + 1.0e-16, // CaHCO3+ + 1.0e-16, // CaSO4 + 1.0e-16, // CaCl+ + 1.0e-16, // CaCl2 + 1.0e-16, // MgSO4 + 1.0e-16, // NaSO4- + 3.76e-1, // H+ + 3.76e-1, // HCO3- + 3.87e-2, // Ca+2 + 3.21e-2, // SO4-2 + 1.89, // Cl- + 1.65e-2, // Mg+2 + 1.09 // Na+1 + }; + + double const surfaceArea[1] = { 1e6 }; // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + + double const expectedReactionRates[1] = { 1.5488389999999999 }; // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + + double const expectedReactionRatesDerivatives[1][16] = + { + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.0877659574468075e-03, -3.0877659574468075e-03, -2.9999999999999997e-02, 0, 0, 0, 0 } + }; + computeReactionRatesTest< double, false >( carbonateSystem.kineticReactionsParameters(), + initialSpeciesConcentration, + surfaceArea, + expectedReactionRates, + expectedReactionRatesDerivatives ); + computeReactionRatesTest< double, true >( carbonateSystem.kineticReactionsParameters(), + initialSpeciesConcentration, + surfaceArea, + expectedReactionRates, + expectedReactionRatesDerivatives ); +} //****************************************************************************** From adf47412bc3b0185372970a554b28190939c73ed Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Wed, 29 Oct 2025 16:57:24 -0700 Subject: [PATCH 10/10] fixed cuda --- src/reactions/reactionsSystems/Parameters.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index a883ce8..54cddaa 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -105,7 +105,7 @@ struct KineticReactionsParameters HPCREACT_HOST_DEVICE RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; } HPCREACT_HOST_DEVICE RealType equilibriumConstant( IndexType const r ) const { return m_rateConstantForward[r] / m_rateConstantReverse[r]; } - IntType reactionRatesUpdateOption() const { return m_reactionRatesUpdateOption; } + HPCREACT_HOST_DEVICE IntType reactionRatesUpdateOption() const { return m_reactionRatesUpdateOption; } CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward;