From 6c89cc570fa725328ae802ccc585771feab6f015 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Wed, 26 Nov 2025 17:17:29 -0800 Subject: [PATCH 1/5] add initial constitutive relations api --- src/CMakeLists.txt | 4 +- src/constitutive/activity/Bdot.hpp | 63 +++++++++++++++++++ .../ionicStrength/SpeciatedIonicStrength.hpp | 59 +++++++++++++++++ .../StoichiometricIonicStrength.hpp | 42 +++++++++++++ src/constitutive/unitTests/CMakeLists.txt | 22 +++++++ src/constitutive/unitTests/testBdot.cpp | 49 +++++++++++++++ .../unitTests/testIonicStrength.cpp | 48 ++++++++++++++ src/reactions/reactionsSystems/Parameters.hpp | 14 +++++ 8 files changed, 300 insertions(+), 1 deletion(-) create mode 100644 src/constitutive/activity/Bdot.hpp create mode 100644 src/constitutive/ionicStrength/SpeciatedIonicStrength.hpp create mode 100644 src/constitutive/ionicStrength/StoichiometricIonicStrength.hpp create mode 100644 src/constitutive/unitTests/CMakeLists.txt create mode 100644 src/constitutive/unitTests/testBdot.cpp create mode 100644 src/constitutive/unitTests/testIonicStrength.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 096ba7c..0b7dbf8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,7 @@ set( hpcReact_headers common/macros.hpp common/CArrayWrapper.hpp + constitutive/activity/Bdot.hpp reactions/exampleSystems/BulkGeneric.hpp reactions/geochemistry/Carbonate.hpp reactions/geochemistry/Forge.hpp @@ -67,10 +68,11 @@ message(STATUS "HPCReact/src CMAKE_CURRENT_SOURCE_DIR: ${CMAKE_CURRENT_SOURCE_DI # hpcReact_add_code_checks( PREFIX hpcReact # EXCLUDES "blt/*" ) +add_subdirectory( common/unitTests ) +add_subdirectory( constitutive/unitTests) add_subdirectory( reactions/exampleSystems/unitTests ) add_subdirectory( reactions/geochemistry/unitTests ) add_subdirectory( reactions/massActions/unitTests ) -add_subdirectory( common/unitTests ) add_subdirectory( docs ) hpcReact_add_code_checks( PREFIX hpcReact diff --git a/src/constitutive/activity/Bdot.hpp b/src/constitutive/activity/Bdot.hpp new file mode 100644 index 0000000..b7dc360 --- /dev/null +++ b/src/constitutive/activity/Bdot.hpp @@ -0,0 +1,63 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "common/macros.hpp" + +namespace hpcReact +{ + +template< typename REAL_TYPE, + typename INDEX_TYPE, + typename IONIC_STRENGTH_TYPE > +class Bdot +{ +public: + using RealType = REAL_TYPE; + using IndexType = INDEX_TYPE; + +struct Params : public IONIC_STRENGTH_TYPE::PARAMS +{ + +}; + + +template< typename ARRAY_1D_TO_CONST, + typename ARRAY_1D > +static inline HPCREACT_HOST_DEVICE +void +calculateActivities( IONIC_STRENGTH_TYPE::PARAMS const & params, + ARRAY_1D_TO_CONST const & speciesConcentrations, + ARRAY_1D & activities ) +{ + + RealType const ionicStrength = IONIC_STRENGTH_TYPE::calculate( params, speciesConcentrations ); + RealType const sqrtI = sqrt(ionicStrength); + RealType const A_gamma = 2; + RealType const B_gamma = 1.6; + auto const & speciesCharge = params.m_speciesCharge; + auto const & a = params.m_ionSizeParameter; + auto const & b = params.m_bdotParameter; + + + + constexpr IndexType numSpecies = params.numSpecies(); + for( IndexType i=0; i +class SpeciatedIonicStrength +{ +public: + using RealType = REAL_TYPE; + using IndexType = INDEX_TYPE; + +struct Params +{ + + HPCREACT_HOST_DEVICE static constexpr IndexType numSpecies() { return NUM_SPECIES; } + HPCREACT_HOST_DEVICE constexpr CArrayWrapper< RealType, NUM_SPECIES > & speciesCharge() { return m_speciesCharge; } + + CArrayWrapper< RealType, NUM_SPECIES > m_speciesCharge; +}; + + +template< typename ARRAY_1D_TO_CONST > +static inline HPCREACT_HOST_DEVICE +REAL_TYPE +calculate( Params const & params, + ARRAY_1D_TO_CONST const & speciesConcentration ) +{ + REAL_TYPE I = 0.0; + auto const & numSpecies = params.numSpecies(); + auto const & speciesCharge = params.m_speciesCharge; + for( int i=0; i +class SpeciatedIonicStrength +{ +public: + +template< typename ARRAY_1D_TO_CONST > +static inline HPCREACT_HOST_DEVICE +REAL_TYPE +calculate( ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D_TO_CONST const & speciesCharge, + int const numSpecies ) +{ + REAL_TYPE I = 0.0; + for( int i=0; i + +using namespace hpcReact; + + + + +constexpr SpeciatedIonicStrength::Params<3> testParams +{ + // Species charge + { 1.0, -1.0, 2.0 } +}; + +TEST( testBdot, testIonicStrength ) +{ + double speciesConcentration[ testParams.numSpecies() ] = { 0.1, 0.2, 0.3 }; + + double I = SpeciatedIonicStrength< double, int >::calculate( testParams, + speciesConcentration ); + + double expectedI = 0.5 * ( 0.1 * 1.0 * 1.0 + 0.2 * (-1.0) * (-1.0) + 0.3 * 2.0 * 2.0 ); + EXPECT_DOUBLE_EQ( I, expectedI ); + +} + + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} diff --git a/src/constitutive/unitTests/testIonicStrength.cpp b/src/constitutive/unitTests/testIonicStrength.cpp new file mode 100644 index 0000000..c29ba04 --- /dev/null +++ b/src/constitutive/unitTests/testIonicStrength.cpp @@ -0,0 +1,48 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" +#include "common/CArrayWrapper.hpp" +#include "common/pmpl.hpp" + +#include + +using namespace hpcReact; + + +using IonicStrength = SpeciatedIonicStrength; + +constexpr IonicStrength::Params testParams +{ + // Species charge + { 1.0, -1.0, 2.0 } +}; + +TEST( testBdot, testIonicStrength ) +{ + double speciesConcentration[ testParams.numSpecies() ] = { 0.1, 0.2, 0.3 }; + + double I = IonicStrength::calculate( testParams, + speciesConcentration ); + + double expectedI = 0.5 * ( 0.1 * 1.0 * 1.0 + 0.2 * (-1.0) * (-1.0) + 0.3 * 2.0 * 2.0 ); + EXPECT_DOUBLE_EQ( I, expectedI ); + +} + + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index 54cddaa..ca1557c 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -116,6 +116,10 @@ struct KineticReactionsParameters }; + + + + template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, @@ -136,12 +140,18 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse, CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag, + CArrayWrapper< RealType, NUM_SPECIES > const & speciesCharge, + CArrayWrapper< RealType, NUM_SPECIES > const & ionSizeParameter, + CArrayWrapper< RealType, NUM_SPECIES > const & bdotParameter, IntType const reactionRatesUpdateOption = 1 ): m_stoichiometricMatrix( stoichiometricMatrix ), m_equilibriumConstant( equilibriumConstant ), m_rateConstantForward( rateConstantForward ), m_rateConstantReverse( rateConstantReverse ), m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ), + m_speciesCharge( speciesCharge ), + m_ionSizeParameter( ionSizeParameter ), + m_bdotParameter( bdotParameter ), m_reactionRatesUpdateOption( reactionRatesUpdateOption ) {} @@ -250,6 +260,10 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag; + CArrayWrapper< RealType, NUM_SPECIES > m_speciesCharge; + CArrayWrapper< RealType, NUM_SPECIES > m_ionSizeParameter; + CArrayWrapper< RealType, NUM_SPECIES > m_bdotParameter; + IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. }; From 1413a1bbf001a50d4db1b55a7527fd5d983dd277 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Thu, 8 Jan 2026 10:32:30 -0800 Subject: [PATCH 2/5] compiles...still need to deal with params and actual calls for c->activity --- src/constitutive/activity/Bdot.hpp | 40 ++++- src/constitutive/activity/DebyeHuckel.hpp | 155 ++++++++++++++++++ src/constitutive/activity/Identity.hpp | 45 +++++ .../StoichiometricIonicStrength.hpp | 10 +- src/constitutive/unitTests/testBdot.cpp | 4 +- src/reactions/exampleSystems/BulkGeneric.hpp | 9 +- .../unitTests/testKineticReactions.cpp | 10 +- .../unitTests/testMomasEasyCase.cpp | 9 +- .../unitTests/testMomasMediumCase.cpp | 9 +- .../testGeochemicalEquilibriumReactions.cpp | 9 +- .../testGeochemicalMixedReactions.cpp | 5 +- .../reactionsSystems/EquilibriumReactions.hpp | 3 +- ...ionsAggregatePrimaryConcentration_impl.hpp | 18 +- ...uilibriumReactionsReactionExtents_impl.hpp | 12 +- .../reactionsSystems/KineticReactions.hpp | 1 + .../KineticReactions_impl.hpp | 59 ++++++- .../MixedEquilibriumKineticReactions.hpp | 3 +- .../MixedEquilibriumKineticReactions_impl.hpp | 6 + src/reactions/reactionsSystems/Parameters.hpp | 10 -- .../equilibriumReactionsTestUtilities.hpp | 8 +- .../kineticReactionsTestUtilities.hpp | 5 + .../mixedReactionsTestUtilities.hpp | 5 +- 22 files changed, 363 insertions(+), 72 deletions(-) create mode 100644 src/constitutive/activity/DebyeHuckel.hpp create mode 100644 src/constitutive/activity/Identity.hpp diff --git a/src/constitutive/activity/Bdot.hpp b/src/constitutive/activity/Bdot.hpp index b7dc360..57c37da 100644 --- a/src/constitutive/activity/Bdot.hpp +++ b/src/constitutive/activity/Bdot.hpp @@ -10,7 +10,7 @@ */ #pragma once -#include "common/macros.hpp" +#include "DebyeHuckel.hpp" namespace hpcReact { @@ -24,39 +24,61 @@ class Bdot using RealType = REAL_TYPE; using IndexType = INDEX_TYPE; -struct Params : public IONIC_STRENGTH_TYPE::PARAMS -{ + +struct Params : public IONIC_STRENGTH_TYPE::Params +{ + RealType m_ionSizeParameter[IONIC_STRENGTH_TYPE::Params::numSpecies()]; + RealType m_bdotParameter[IONIC_STRENGTH_TYPE::Params::numSpecies()]; }; + template< typename ARRAY_1D_TO_CONST, typename ARRAY_1D > static inline HPCREACT_HOST_DEVICE void -calculateActivities( IONIC_STRENGTH_TYPE::PARAMS const & params, +calculateActivities( Params const & params, ARRAY_1D_TO_CONST const & speciesConcentrations, ARRAY_1D & activities ) { RealType const ionicStrength = IONIC_STRENGTH_TYPE::calculate( params, speciesConcentrations ); RealType const sqrtI = sqrt(ionicStrength); - RealType const A_gamma = 2; - RealType const B_gamma = 1.6; + RealType const rho_w = 997.0479; // kg/m3 + RealType const eps_r = 78.54; // dimensionless + RealType const T_K = 298.15; + RealType const A_gamma = DebyeHuckel::A_gamma( T_K, rho_w, eps_r ); + RealType const B_gamma = DebyeHuckel::B_gamma( T_K, rho_w, eps_r ); auto const & speciesCharge = params.m_speciesCharge; auto const & a = params.m_ionSizeParameter; auto const & b = params.m_bdotParameter; - constexpr IndexType numSpecies = params.numSpecies(); + const IndexType numSpecies = params.numSpecies(); for( IndexType i=0; i::log10_gamma( sqrtI, + speciesCharge[i], + a[i], + A_gamma, + B_gamma ); + activities[i] = DebyeHuckel_term + b[i] * ionicStrength; } } +private: + // ------------------------------------------------------------- + // Physical constants (CODATA 2018-ish), SI units + // ------------------------------------------------------------- + static constexpr double pi = 3.14159265358979323846; + static constexpr double e0 = 8.8541878128e-12; // vacuum permittivity [F/m] + static constexpr double eChg = 1.602176634e-19; // elementary charge [C] + static constexpr double kB = 1.380649e-23; // Boltzmann constant [J/K] + static constexpr double NA = 6.02214076e23; // Avogadro [1/mol] + static constexpr double ln10 = 2.3025850929940459; + }; diff --git a/src/constitutive/activity/DebyeHuckel.hpp b/src/constitutive/activity/DebyeHuckel.hpp new file mode 100644 index 0000000..b174ac7 --- /dev/null +++ b/src/constitutive/activity/DebyeHuckel.hpp @@ -0,0 +1,155 @@ +#pragma once + +#include "common/macros.hpp" + +#include + +/** + * @file DebyeHuckel.hpp + * @brief Debye–Hückel A^γ and B parameters for aqueous electrolytes. + * + * This header provides helper functions to compute the Debye–Hückel + * parameters A^γ and B in their "native" (natural-log) form for + * molal (mol/kg) activity-coefficient models. + * + * The functions are expressed in terms of fundamental physical constants + * and water properties (density and relative permittivity). They can be + * used directly in Debye–Hückel or extended Debye–Hückel/B-dot models. + */ + +template< typename REAL_TYPE > +class DebyeHuckel +{ +public: + using RealType = REAL_TYPE; + + /// π (pi). + static constexpr RealType pi = 3.141592653589793e+00; + + /// Vacuum permittivity ε₀ [F/m]. + static constexpr RealType e0 = 8.854187812800001e-12; + + /// Elementary charge e [C]. + static constexpr RealType eChg = 1.602176634000000e-19; + + /// Boltzmann constant k_B [J/K]. + static constexpr RealType kB = 1.380649000000000e-23; + + /// Avogadro constant N_A [1/mol]. + static constexpr RealType NA = 6.022140760000000e+23; + + /// Natural logarithm of 10, ln(10). + static constexpr RealType ln10 = 2.302585092994046e+00; + + /// Inverse of ln(10). + static constexpr RealType invln10 = 4.342944819032518e-01; + + + // ------------------------------------------------------------- + // Debye–Hückel A^γ (natural log, molal scale) + // ------------------------------------------------------------- + + /** + * @brief Debye–Hückel A^γ parameter in natural-log form. + * + * Computes the coefficient A^γ(T,ρ,ε_r) used in the Debye–Hückel + * expression for the natural logarithm of the activity coefficient: + * + * \f[ + * \ln \gamma_i = + * - A^\gamma_{\ln}(T,P) \, z_i^2 + * \frac{\sqrt{I}}{1 + B(T,P)\, a_i \sqrt{I}} + * \f] + * + * where: + * - \f$ I \f$ is ionic strength in mol/kg (molal), + * - \f$ z_i \f$ is the ionic charge, + * - \f$ a_i \f$ is the ion-size parameter (length), + * - A^γ is independent of the log base (this function is for ln). + * + * The implementation follows the "native" Debye–Hückel form, + * using fundamental physical constants without any 1/ln(10) factors. + * + * @param T_K Temperature in kelvin [K]. + * @param rho_w Density of water in g/L (≈ kg/m³ numerically). + * @param eps_r Relative permittivity (dielectric constant) of water. + * @return A^γ in units consistent with molal ionic strength, for use + * in ln(γ) expressions. + */ + static inline HPCREACT_HOST_DEVICE + RealType A_gamma( RealType const T_K, + RealType const rho_w, + RealType const eps_r ) + { + RealType const num = ::pow( eChg, 3.0 ) * ::sqrt( 2.0 * pi * NA * rho_w ); + RealType const den = ::pow( 4.0 * pi * e0 * eps_r * kB * T_K, 1.5 ); + return num / den; + } + + + // ------------------------------------------------------------- + // Debye–Hückel B (natural log, molal scale) + // ------------------------------------------------------------- + + /** + * @brief Debye–Hückel B parameter in natural-log form. + * + * Computes the Debye–Hückel length-scale parameter B(T,ρ,ε_r) used + * in the extended Debye–Hückel law: + * + * \f[ + * \ln \gamma_i = + * - A^\gamma_{\ln}(T,P) \, z_i^2 + * \frac{\sqrt{I}}{1 + B(T,P)\, a_i \sqrt{I}} \; , + * \f] + * + * where: + * - \f$ I \f$ is ionic strength in mol/kg, + * - \f$ a_i \f$ is an ion-size parameter (length). + * + * The combination \f$ B a_i \sqrt{I} \f$ is dimensionless; the + * absolute units of B therefore depend on the length units chosen + * for \f$ a_i \f$. + * + * @param T_K Temperature in kelvin [K]. + * @param rho_w Density of water in g/L (≈ kg/m³ numerically). + * @param eps_r Relative permittivity (dielectric constant) of water. + * @return B parameter for use in ln(γ) expressions. + */ + static inline HPCREACT_HOST_DEVICE + RealType B_gamma( RealType const T_K, + RealType const rho_w, + RealType const eps_r ) + { + RealType const num = 2.0 * NA * rho_w * eChg * eChg; + RealType const den = e0 * eps_r * kB * T_K; + return ::sqrt( num / den ); + } + + + static inline HPCREACT_HOST_DEVICE + RealType log10_gamma( RealType const I, + RealType const zi, + RealType const ai, + RealType const T_K ) + { + RealType const sqrtI = sqrt(I); + RealType const A = A_gamma(T_K); + RealType const B = B_gamma(T_K); + RealType const denom = 1 + B * ai * sqrtI; + return -A * zi * zi * sqrtI / denom; + } + + + static inline HPCREACT_HOST_DEVICE + RealType log10_gamma( RealType const sqrtI, + RealType const zi, + RealType const ai, + RealType const A, + RealType const B ) + { + RealType const denom = 1 + B * ai * sqrtI; + return -A * zi * zi * sqrtI / denom; + } + +}; diff --git a/src/constitutive/activity/Identity.hpp b/src/constitutive/activity/Identity.hpp new file mode 100644 index 0000000..ea61d33 --- /dev/null +++ b/src/constitutive/activity/Identity.hpp @@ -0,0 +1,45 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "common/macros.hpp" + +namespace hpcReact +{ + +template< typename REAL_TYPE, + typename INDEX_TYPE, + typename IONIC_STRENGTH_TYPE > +class Identity +{ +public: + using RealType = REAL_TYPE; + using IndexType = INDEX_TYPE; + + template< typename ARRAY_1D_TO_CONST, + typename ARRAY_1D > + static inline HPCREACT_HOST_DEVICE + void + calculateActivities( IONIC_STRENGTH_TYPE::PARAMS const & params, + ARRAY_1D_TO_CONST const & speciesConcentrations, + ARRAY_1D & activities ) + { + constexpr IndexType numSpecies = params.numSpecies(); + for( IndexType i=0; i -class SpeciatedIonicStrength +class StoichiometricIonicStrength { public: @@ -27,13 +27,7 @@ calculate( ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D_TO_CONST const & speciesCharge, int const numSpecies ) { - REAL_TYPE I = 0.0; - for( int i=0; i::Params<3> testParams +constexpr SpeciatedIonicStrength::Params testParams { // Species charge { 1.0, -1.0, 2.0 } @@ -32,7 +32,7 @@ TEST( testBdot, testIonicStrength ) { double speciesConcentration[ testParams.numSpecies() ] = { 0.1, 0.2, 0.3 }; - double I = SpeciatedIonicStrength< double, int >::calculate( testParams, + double I = SpeciatedIonicStrength< double, int, 3 >::calculate( testParams, speciesConcentration ); double expectedI = 0.5 * ( 0.1 * 1.0 * 1.0 + 0.2 * (-1.0) * (-1.0) + 0.3 * 2.0 * 2.0 ); diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp index 77547ce..aea179b 100644 --- a/src/reactions/exampleSystems/BulkGeneric.hpp +++ b/src/reactions/exampleSystems/BulkGeneric.hpp @@ -45,7 +45,7 @@ namespace bulkGeneric // um1Constants }; -using simpleKineticTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 0 >; +using simpleKineticTestType = reactionsSystems::KineticReactionsParameters< double, int, int, 5, 2 >; constexpr simpleKineticTestType @@ -56,16 +56,15 @@ simpleKineticTestRateParams = { -2, 1, 1, 0, 0 }, { 0, 0, -1, -1, 2 } }, - // equilibrium constants - { 1.0, 1.0 }, // forward rate constants { 1.0, 0.5 }, // reverse rate constants { 1.0, 0.5 }, - // flag of mobile secondary species - { 1, 1 }, + // equilibrium constants + { 1.0, 1.0 }, // Use the forward and reverse to calculate the kinetic reaction rates 0 + }; using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >; diff --git a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp index 1286169..7c5bbdc 100644 --- a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp +++ b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp @@ -29,12 +29,12 @@ TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams 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(), + computeReactionRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams, initialSpeciesConcentration, surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); - computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(), + computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams, initialSpeciesConcentration, surfaceArea, // No use. Just to pass something here expectedReactionRates, @@ -52,12 +52,12 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams { 0.0, 0.0, -0.5, -0.25, 0.0 }, { 0.0, 0.0, 1.0, 0.5, 0.0 } }; - computeSpeciesRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(), + computeSpeciesRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams, initialSpeciesConcentration, expectedSpeciesRates, expectedSpeciesRatesDerivatives ); - computeSpeciesRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(), + computeSpeciesRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams, initialSpeciesConcentration, expectedSpeciesRates, expectedSpeciesRatesDerivatives ); @@ -70,7 +70,7 @@ TEST( testKineticReactions, testTimeStep ) double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01, 7.02014627734060e-01, 5.95970744531880e-01 }; - timeStepTest< double, false >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(), + timeStepTest< double, false >( bulkGeneric::simpleKineticTestRateParams, 2.0, 10, initialSpeciesConcentration, diff --git a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp index 33b4669..6745a96 100644 --- a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp +++ b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp @@ -23,11 +23,14 @@ using namespace hpcReact::unitTest_utilities; void testMoMasAllEquilibriumHelper() { - using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, - int, - int >; static constexpr int numPrimarySpecies = hpcReact::MoMasBenchmark::easyCaseParams.numPrimarySpecies(); + static constexpr int numSpecies = hpcReact::MoMasBenchmark::easyCaseParams.numSpecies(); + + using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, + int, + int, + Bdot< double, int, SpeciatedIonicStrength< double, int, numSpecies > >>; double logPrimarySpeciesConcentration[numPrimarySpecies]; diff --git a/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp index 534cd03..ec86b19 100644 --- a/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp +++ b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp @@ -22,11 +22,14 @@ using namespace hpcReact::unitTest_utilities; void testMoMasMediumEquilibriumHelper() { - using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, - int, - int >; static constexpr int numPrimarySpecies = hpcReact::MoMasBenchmark::mediumCaseParams.numPrimarySpecies(); + static constexpr int numSpecies = hpcReact::MoMasBenchmark::mediumCaseParams.numSpecies(); + + using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, + int, + int, + Bdot< double, int, SpeciatedIonicStrength< double, int, numSpecies > >>; diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp index ee38825..9f9fd3e 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp @@ -101,11 +101,14 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium ) TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 ) { - using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, - int, - int >; static constexpr int numPrimarySpecies = hpcReact::geochemistry::carbonateSystemAllEquilibrium.numPrimarySpecies(); + static constexpr int numSpecies = hpcReact::geochemistry::carbonateSystemAllEquilibrium.numSpecies(); + + using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, + int, + int, + Bdot< double, int, SpeciatedIonicStrength< double, int, numSpecies > >>; double const initialPrimarySpeciesConcentration[numPrimarySpecies] = { diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp index ce3b244..59f90bf 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp @@ -11,6 +11,8 @@ #include "reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp" #include "../GeochemicalSystems.hpp" +#include "constitutive/activity/Bdot.hpp" +#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" using namespace hpcReact; @@ -50,7 +52,8 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem ) 1.070434904554991 // Na+1 }; - timeStepTest< double, true >( carbonateSystem, + using ActivityModelType = Bdot< double, int, SpeciatedIonicStrength< double, int, carbonateSystemType::numSpecies() > >; + timeStepTest< double, true, ActivityModelType >( carbonateSystem, 1.0, 10, initialAggregateSpeciesConcentration, diff --git a/src/reactions/reactionsSystems/EquilibriumReactions.hpp b/src/reactions/reactionsSystems/EquilibriumReactions.hpp index 447a219..8e0b79f 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactions.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactions.hpp @@ -35,7 +35,8 @@ namespace reactionsSystems */ template< typename REAL_TYPE, typename INT_TYPE, - typename INDEX_TYPE > + typename INDEX_TYPE, + typename ACTIVITY_MODEL > class EquilibriumReactions { public: diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index 5dad38c..b6a516d 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -23,7 +23,8 @@ namespace reactionsSystems template< typename REAL_TYPE, typename INT_TYPE, - typename INDEX_TYPE > + typename INDEX_TYPE, + typename ACTIVITY_MODEL > template< typename PARAMS_DATA, typename ARRAY_1D, typename ARRAY_1D_TO_CONST, @@ -34,7 +35,8 @@ inline void EquilibriumReactions< REAL_TYPE, INT_TYPE, - INDEX_TYPE >::computeResidualAndJacobianAggregatePrimaryConcentrations( RealType const & temperature, + INDEX_TYPE, + ACTIVITY_MODEL >::computeResidualAndJacobianAggregatePrimaryConcentrations( RealType const & temperature, PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & targetAggregatePrimaryConcentrations, ARRAY_1D_TO_CONST2 const & logPrimarySpeciesConcentration, @@ -64,7 +66,8 @@ EquilibriumReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, - typename INDEX_TYPE > + typename INDEX_TYPE, + typename ACTIVITY_MODEL > template< typename PARAMS_DATA, typename ARRAY_1D, typename ARRAY_1D_TO_CONST > @@ -72,7 +75,8 @@ HPCREACT_HOST_DEVICE inline void EquilibriumReactions< REAL_TYPE, INT_TYPE, - INDEX_TYPE >::enforceEquilibrium_LogAggregate( REAL_TYPE const & temperature, + INDEX_TYPE, + ACTIVITY_MODEL >::enforceEquilibrium_LogAggregate( REAL_TYPE const & temperature, PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, ARRAY_1D & logPrimarySpeciesConcentration ) @@ -99,7 +103,8 @@ EquilibriumReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, - typename INDEX_TYPE > + typename INDEX_TYPE, + typename ACTIVITY_MODEL > template< typename PARAMS_DATA, typename ARRAY_1D, typename ARRAY_1D_TO_CONST > @@ -107,7 +112,8 @@ HPCREACT_HOST_DEVICE inline void EquilibriumReactions< REAL_TYPE, INT_TYPE, - INDEX_TYPE >::enforceEquilibrium_Aggregate( REAL_TYPE const & temperature, + INDEX_TYPE, + ACTIVITY_MODEL >::enforceEquilibrium_Aggregate( REAL_TYPE const & temperature, PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & targetAggregatePrimarySpeciesConcentration, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp index d776644..d3742d1 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp @@ -24,7 +24,8 @@ constexpr bool debugPrinting = false; template< typename REAL_TYPE, typename INT_TYPE, - typename INDEX_TYPE > + typename INDEX_TYPE, + typename ACTIVITY_MODEL > template< typename PARAMS_DATA, typename ARRAY_1D, typename ARRAY_1D_TO_CONST, @@ -34,7 +35,8 @@ HPCREACT_HOST_DEVICE inline void EquilibriumReactions< REAL_TYPE, INT_TYPE, - INDEX_TYPE >::computeResidualAndJacobianReactionExtents( REAL_TYPE const & temperature, + INDEX_TYPE, + ACTIVITY_MODEL >::computeResidualAndJacobianReactionExtents( REAL_TYPE const & temperature, PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & speciesConcentration0, ARRAY_1D_TO_CONST2 const & xi, @@ -112,7 +114,8 @@ EquilibriumReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, - typename INDEX_TYPE > + typename INDEX_TYPE, + typename ACTIVITY_MODEL > template< typename PARAMS_DATA, typename ARRAY_1D, typename ARRAY_1D_TO_CONST > @@ -120,7 +123,8 @@ HPCREACT_HOST_DEVICE inline void EquilibriumReactions< REAL_TYPE, INT_TYPE, - INDEX_TYPE >::enforceEquilibrium_Extents( REAL_TYPE const & temperature, + INDEX_TYPE, + ACTIVITY_MODEL >::enforceEquilibrium_Extents( REAL_TYPE const & temperature, PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & speciesConcentration0, ARRAY_1D & speciesConcentration ) diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp index a5f5575..4c6a3d7 100644 --- a/src/reactions/reactionsSystems/KineticReactions.hpp +++ b/src/reactions/reactionsSystems/KineticReactions.hpp @@ -38,6 +38,7 @@ namespace reactionsSystems template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > class KineticReactions { diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp index b38a69e..8307036 100644 --- a/src/reactions/reactionsSystems/KineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/KineticReactions_impl.hpp @@ -34,6 +34,7 @@ namespace reactionsSystems template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > template< typename PARAMS_DATA, bool CALCULATE_DERIVATIVES, @@ -44,6 +45,7 @@ HPCREACT_HOST_DEVICE inline void KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, + ACTIVITY_MODEL, LOGE_CONCENTRATION >::computeReactionRates_impl( RealType const &, //temperature, PARAMS_DATA const & params, @@ -57,6 +59,23 @@ KineticReactions< REAL_TYPE, HPCREACT_UNUSED_VAR( reactionRatesDerivatives ); } + RealType activities[ PARAMS_DATA::numSpecies() ]; + + if constexpr( LOGE_CONCENTRATION ) + { + RealType expSpeciesConcentration[ PARAMS_DATA::numSpecies() ]; + for( INT_TYPE i=0; i 0.0 ) { - productConcReverse += s_ri * speciesConcentration[i]; + productConcReverse += s_ri * activities[i]; } } @@ -130,7 +149,7 @@ KineticReactions< REAL_TYPE, { RealType const s_ri = params.stoichiometricMatrix( r, i ); - RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], fabs( s_ri ) ) : 0.0; + RealType const productTerm_i = activities[i] > 1e-100 ? pow( activities[i], fabs( s_ri ) ) : 0.0; if( s_ri < 0.0 ) { @@ -150,7 +169,7 @@ KineticReactions< REAL_TYPE, { if( i==j ) { - dProductConcForward_dC[j] *= -s_ri * pow( speciesConcentration[i], -s_ri-1 ); + dProductConcForward_dC[j] *= -s_ri * pow( activities[i], -s_ri-1 ); dProductConcReverse_dC[j] = 0.0; } else @@ -165,7 +184,7 @@ KineticReactions< REAL_TYPE, { if( i==j ) { - dProductConcReverse_dC[j] *= s_ri * pow( speciesConcentration[i], s_ri-1 ); + dProductConcReverse_dC[j] *= s_ri * pow( activities[i], s_ri-1 ); dProductConcForward_dC[j] = 0.0; } else @@ -197,6 +216,7 @@ KineticReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > template< typename PARAMS_DATA, bool CALCULATE_DERIVATIVES, @@ -208,6 +228,7 @@ HPCREACT_HOST_DEVICE inline void KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, + ACTIVITY_MODEL, LOGE_CONCENTRATION >::computeReactionRatesQuotient_impl( RealType const &, //temperature, PARAMS_DATA const & params, @@ -221,6 +242,22 @@ KineticReactions< REAL_TYPE, HPCREACT_UNUSED_VAR( reactionRatesDerivatives ); } + RealType activities[ PARAMS_DATA::numSpecies() ]; + + if constexpr( LOGE_CONCENTRATION ) + { + RealType expSpeciesConcentration[ PARAMS_DATA::numSpecies() ]; + for( INT_TYPE i=0; i 1e-100 ? pow( speciesConcentration[i], s_ri ) : 0.0; + RealType const productTerm_i = activities[i] > 1e-100 ? pow( activities[i], s_ri ) : 0.0; quotient *= productTerm_i; } @@ -279,7 +316,7 @@ KineticReactions< REAL_TYPE, RealType const s_ri = params.stoichiometricMatrix( r, i ); if( s_ri > 0.0 || s_ri < 0.0 ) { - reactionRatesDerivatives( r, i ) = -rateConstant * surfaceArea[r] * s_ri * quotient / ( equilibriumConstant * speciesConcentration[i] ); + reactionRatesDerivatives( r, i ) = -rateConstant * surfaceArea[r] * s_ri * quotient / ( equilibriumConstant * activities[i] ); } else { @@ -296,6 +333,7 @@ KineticReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > template< typename PARAMS_DATA, bool CALCULATE_DERIVATIVES, @@ -306,6 +344,7 @@ HPCREACT_HOST_DEVICE inline void KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, + ACTIVITY_MODEL, LOGE_CONCENTRATION >::computeSpeciesRates_impl( RealType const & temperature, PARAMS_DATA const & params, @@ -321,6 +360,8 @@ KineticReactions< REAL_TYPE, HPCREACT_UNUSED_VAR( speciesRatesDerivatives ); } + + computeReactionRates< PARAMS_DATA >( temperature, params, speciesConcentration, reactionRates, reactionRatesDerivatives ); for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) @@ -351,6 +392,7 @@ KineticReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > template< typename PARAMS_DATA, typename ARRAY_1D, @@ -360,6 +402,7 @@ HPCREACT_HOST_DEVICE inline void KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, + ACTIVITY_MODEL, LOGE_CONCENTRATION >::timeStep( RealType const dt, RealType const & temperature, PARAMS_DATA const & params, diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp index 9f0f9de..0950cb6 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp @@ -37,6 +37,7 @@ namespace reactionsSystems template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > class MixedEquilibriumKineticReactions { @@ -52,7 +53,7 @@ class MixedEquilibriumKineticReactions using IndexType = INDEX_TYPE; /// Type alias for the Kinetic reactions type used in the class. - using kineticReactions = KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, LOGE_CONCENTRATION >; + using kineticReactions = KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, ACTIVITY_MODEL, LOGE_CONCENTRATION >; /** * @brief Update a mixed chemical system by computing secondary species concentrations, diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp index f24e4d6..92aabd9 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp @@ -30,6 +30,7 @@ namespace reactionsSystems template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > template< typename PARAMS_DATA, typename ARRAY_1D_TO_CONST, @@ -43,6 +44,7 @@ HPCREACT_HOST_DEVICE inline void MixedEquilibriumKineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, + ACTIVITY_MODEL, LOGE_CONCENTRATION >::updateMixedSystem_impl( RealType const & temperature, PARAMS_DATA const & params, @@ -114,6 +116,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > template< typename PARAMS_DATA, typename ARRAY_1D_TO_CONST, @@ -125,6 +128,7 @@ HPCREACT_HOST_DEVICE inline void MixedEquilibriumKineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, + ACTIVITY_MODEL, LOGE_CONCENTRATION >::computeReactionRates_impl( RealType const & temperature, PARAMS_DATA const & params, @@ -188,6 +192,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, + typename ACTIVITY_MODEL, bool LOGE_CONCENTRATION > template< typename PARAMS_DATA, typename ARRAY_1D_TO_CONST, @@ -200,6 +205,7 @@ HPCREACT_HOST_DEVICE inline void MixedEquilibriumKineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, + ACTIVITY_MODEL, LOGE_CONCENTRATION >::computeAggregateSpeciesRates_impl( PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & speciesConcentration, diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index ca1557c..ae87b3d 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -140,18 +140,12 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse, CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag, - CArrayWrapper< RealType, NUM_SPECIES > const & speciesCharge, - CArrayWrapper< RealType, NUM_SPECIES > const & ionSizeParameter, - CArrayWrapper< RealType, NUM_SPECIES > const & bdotParameter, IntType const reactionRatesUpdateOption = 1 ): m_stoichiometricMatrix( stoichiometricMatrix ), m_equilibriumConstant( equilibriumConstant ), m_rateConstantForward( rateConstantForward ), m_rateConstantReverse( rateConstantReverse ), m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ), - m_speciesCharge( speciesCharge ), - m_ionSizeParameter( ionSizeParameter ), - m_bdotParameter( bdotParameter ), m_reactionRatesUpdateOption( reactionRatesUpdateOption ) {} @@ -260,10 +254,6 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag; - CArrayWrapper< RealType, NUM_SPECIES > m_speciesCharge; - CArrayWrapper< RealType, NUM_SPECIES > m_ionSizeParameter; - CArrayWrapper< RealType, NUM_SPECIES > m_bdotParameter; - IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. }; diff --git a/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp index 58b0bdb..d259c0e 100644 --- a/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp @@ -13,6 +13,8 @@ #include "reactions/reactionsSystems/EquilibriumReactions.hpp" #include "common/macros.hpp" #include "common/pmpl.hpp" +#include "constitutive/activity/Bdot.hpp" +#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" #include @@ -58,7 +60,8 @@ void computeResidualAndJacobianTest( PARAMS_DATA const & params, { using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< REAL_TYPE, int, - int >; + int, + Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > > >; static constexpr int numSpecies = PARAMS_DATA::numSpecies(); static constexpr int numReactions = PARAMS_DATA::numReactions(); @@ -128,7 +131,8 @@ void testEnforceEquilibrium( PARAMS_DATA const & params, { using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< REAL_TYPE, int, - int >; + int, + Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > > >; static constexpr int numSpecies = PARAMS_DATA::numSpecies(); diff --git a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp index 3fe18ed..46addd3 100644 --- a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp @@ -15,6 +15,8 @@ #include "common/pmpl.hpp" #include "common/printers.hpp" #include "common/CArrayWrapper.hpp" +#include "constitutive/activity/Bdot.hpp" +#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" #include "reactions/reactionsSystems/KineticReactions.hpp" #include @@ -66,6 +68,7 @@ void computeReactionRatesTest( PARAMS_DATA const & params, using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE, int, int, + Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >, LOGE_CONCENTRATION >; static constexpr int numSpecies = PARAMS_DATA::numSpecies(); @@ -164,6 +167,7 @@ void computeSpeciesRatesTest( PARAMS_DATA const & params, using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE, int, int, + Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >, LOGE_CONCENTRATION >; static constexpr int numSpecies = PARAMS_DATA::numSpecies(); @@ -242,6 +246,7 @@ void timeStepTest( PARAMS_DATA const & params, using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE, int, int, + Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >, LOGE_CONCENTRATION >; static constexpr int numSpecies = PARAMS_DATA::numSpecies(); diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp index 7426c03..0bbd7d4 100644 --- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp @@ -29,6 +29,7 @@ namespace unitTest_utilities //****************************************************************************** template< typename REAL_TYPE, bool LOGE_CONCENTRATION, + typename ACTIVITY_MODEL, typename PARAMS_DATA > void timeStepTest( PARAMS_DATA const & params, REAL_TYPE const dt, @@ -53,10 +54,12 @@ void timeStepTest( PARAMS_DATA const & params, using MixedReactionsType = reactionsSystems::MixedEquilibriumKineticReactions< REAL_TYPE, int, int, + ACTIVITY_MODEL, LOGE_CONCENTRATION >; using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< REAL_TYPE, int, - int >; + int, + ACTIVITY_MODEL >; // constexpr int numSpecies = PARAMS_DATA::numSpecies(); static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); From 806547dc68f36e61309c124d5f53e7aaf84638da Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Thu, 8 Jan 2026 17:34:55 -0800 Subject: [PATCH 3/5] testKineticReactions is function with activity...need updated answers to compare with --- src/constitutive/activity/Identity.hpp | 12 +- src/reactions/exampleSystems/BulkGeneric.hpp | 28 ++++- .../unitTests/testKineticReactions.cpp | 114 ++++++++++++------ .../reactionsSystems/KineticReactions.hpp | 14 +++ .../KineticReactions_impl.hpp | 25 ++-- .../kineticReactionsTestUtilities.hpp | 57 +++++---- 6 files changed, 176 insertions(+), 74 deletions(-) diff --git a/src/constitutive/activity/Identity.hpp b/src/constitutive/activity/Identity.hpp index ea61d33..18a47bc 100644 --- a/src/constitutive/activity/Identity.hpp +++ b/src/constitutive/activity/Identity.hpp @@ -24,15 +24,21 @@ class Identity using RealType = REAL_TYPE; using IndexType = INDEX_TYPE; + struct Params + { + HPCREACT_HOST_DEVICE static constexpr IndexType numSpecies() { return IONIC_STRENGTH_TYPE::Params::numSpecies(); } + }; + template< typename ARRAY_1D_TO_CONST, - typename ARRAY_1D > + typename ARRAY_1D, + typename PARAMS > static inline HPCREACT_HOST_DEVICE void - calculateActivities( IONIC_STRENGTH_TYPE::PARAMS const & params, + calculateActivities( PARAMS const & , ARRAY_1D_TO_CONST const & speciesConcentrations, ARRAY_1D & activities ) { - constexpr IndexType numSpecies = params.numSpecies(); + constexpr IndexType numSpecies = PARAMS::numSpecies(); for( IndexType i=0; i; +using simpleKineticParamsType = reactionsSystems::KineticReactionsParameters< double, int, int, 5, 2 >; constexpr -simpleKineticTestType +simpleKineticParamsType simpleKineticTestRateParams = { // stoichiometric matrix @@ -64,13 +67,30 @@ simpleKineticTestRateParams = { 1.0, 1.0 }, // Use the forward and reverse to calculate the kinetic reaction rates 0 +}; + +using simpleIonicStrengthType = SpeciatedIonicStrength< double, int, 5 >; +using simpleActivityParamsType = Bdot< double, int, simpleIonicStrengthType >::Params; +constexpr +simpleActivityParamsType +simpleActivityTestParams = +{ + // species charge + {{ 2.0, -1.0, 1.0, 0.0, -1.0 }}, + // ion size parameter + { 4.0, 3.5, 3.5, 3.5, 3.5 }, + // bdot parameter + { 0.1, 0.1, 0.1, 0.0, 0.1 } }; -using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >; +Identity< double, int, simpleIonicStrengthType >::Params simpleIdentityActivityTestParams = {}; + + +using simpleMixedParamsType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >; constexpr -simpleTestType +simpleMixedParamsType simpleTestRateParams = { // stoichiometric matrix diff --git a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp index 7c5bbdc..ad720f8 100644 --- a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp +++ b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp @@ -13,6 +13,10 @@ #include "reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp" #include "../BulkGeneric.hpp" +#include "constitutive/activity/Bdot.hpp" +#include "constitutive/activity/Identity.hpp" +#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" + #include @@ -23,27 +27,69 @@ using namespace hpcReact::unitTest_utilities; //****************************************************************************** TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams ) { + using IonicStrengthType = bulkGeneric::simpleIonicStrengthType; + 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, - initialSpeciesConcentration, - surfaceArea, // No use. Just to pass something here - expectedReactionRates, - expectedReactionRatesDerivatives ); - computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams, - initialSpeciesConcentration, - surfaceArea, // No use. Just to pass something here - expectedReactionRates, - expectedReactionRatesDerivatives ); + + { + using ActivityType = Identity< double, int, IonicStrengthType >; + 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, + ActivityType >( bulkGeneric::simpleKineticTestRateParams, + bulkGeneric::simpleIdentityActivityTestParams, + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); + computeReactionRatesTest< double, + true, + ActivityType >( bulkGeneric::simpleKineticTestRateParams, + bulkGeneric::simpleIdentityActivityTestParams, + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); + } + { + using ActivityType = Bdot< double, int, IonicStrengthType >; + 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, + ActivityType >( bulkGeneric::simpleKineticTestRateParams, + bulkGeneric::simpleActivityTestParams, + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); + computeReactionRatesTest< double, + true, + ActivityType >( bulkGeneric::simpleKineticTestRateParams, + bulkGeneric::simpleActivityTestParams, + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); + } + + } TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams ) { + using IonicStrengthType = bulkGeneric::simpleIonicStrengthType; + using ActivityType = Identity< double, int, IonicStrengthType >; + double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; double const expectedSpeciesRates[5] = { -2.0, 1.0, 0.75, -0.25, 0.5 }; double const expectedSpeciesRatesDerivatives[5][5] = { { -4.0, 1.0, 0.0, 0.0, 0.0 }, @@ -52,12 +98,18 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams { 0.0, 0.0, -0.5, -0.25, 0.0 }, { 0.0, 0.0, 1.0, 0.5, 0.0 } }; - computeSpeciesRatesTest< double, false >( bulkGeneric::simpleKineticTestRateParams, + computeSpeciesRatesTest< double, + false, + ActivityType >( bulkGeneric::simpleKineticTestRateParams, + bulkGeneric::simpleIdentityActivityTestParams, initialSpeciesConcentration, expectedSpeciesRates, expectedSpeciesRatesDerivatives ); - computeSpeciesRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams, + computeSpeciesRatesTest< double, + true, + ActivityType >( bulkGeneric::simpleKineticTestRateParams, + bulkGeneric::simpleIdentityActivityTestParams, initialSpeciesConcentration, expectedSpeciesRates, expectedSpeciesRatesDerivatives ); @@ -65,24 +117,18 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams } -TEST( testKineticReactions, testTimeStep ) -{ - double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; - double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01, 7.02014627734060e-01, 5.95970744531880e-01 }; - - timeStepTest< double, false >( bulkGeneric::simpleKineticTestRateParams, - 2.0, - 10, - initialSpeciesConcentration, - expectedSpeciesConcentrations ); - - // ln(c) as the primary variable results in a singular system. - // timeStepTest< double, true >( simpleKineticTestRateParams, - // 2.0, - // 10, - // initialSpeciesConcentration, - // expectedSpeciesConcentrations ); -} +// TEST( testKineticReactions, testTimeStep ) +// { +// double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; +// double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01, 7.02014627734060e-01, 5.95970744531880e-01 }; + +// timeStepTest< double, false >( bulkGeneric::simpleKineticTestRateParams, +// 2.0, +// 10, +// initialSpeciesConcentration, +// expectedSpeciesConcentrations ); + +// } int main( int argc, char * * argv ) { diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp index 4c6a3d7..b53aad1 100644 --- a/src/reactions/reactionsSystems/KineticReactions.hpp +++ b/src/reactions/reactionsSystems/KineticReactions.hpp @@ -63,12 +63,14 @@ class KineticReactions static HPCREACT_HOST_DEVICE inline void computeReactionRates( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & reactionRates, ARRAY_2D & reactionRatesDerivatives ) { computeReactionRates_impl< PARAMS_DATA, true >( temperature, params, + activityParams, speciesConcentration, reactionRates, reactionRatesDerivatives ); @@ -90,12 +92,14 @@ class KineticReactions static HPCREACT_HOST_DEVICE inline void computeReactionRates( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & reactionRates ) { REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions()][PARAMS_DATA::numSpecies()] = { {0.0} }; computeReactionRates_impl< PARAMS_DATA, false >( temperature, params, + activityParams, speciesConcentration, reactionRates, reactionRatesDerivatives ); @@ -128,6 +132,7 @@ class KineticReactions static HPCREACT_HOST_DEVICE inline void computeReactionRates( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D_SA const & surfaceArea, ARRAY_1D & reactionRates, @@ -137,6 +142,7 @@ class KineticReactions { computeReactionRates_impl< PARAMS_DATA, true >( temperature, params, + activityParams, speciesConcentration, reactionRates, reactionRatesDerivatives ); @@ -145,6 +151,7 @@ class KineticReactions { computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature, params, + activityParams, speciesConcentration, surfaceArea, reactionRates, @@ -164,12 +171,14 @@ class KineticReactions static HPCREACT_HOST_DEVICE inline void computeSpeciesRates( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & speciesRates, ARRAY_2D & speciesRatesDerivatives ) { computeSpeciesRates_impl< PARAMS_DATA, true >( temperature, params, + activityParams, speciesConcentration, speciesRates, speciesRatesDerivatives ); @@ -191,12 +200,14 @@ class KineticReactions static HPCREACT_HOST_DEVICE inline void computeSpeciesRates( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & speciesRates ) { char speciesRatesDerivatives; computeSpeciesRates_impl< PARAMS_DATA, false >( temperature, params, + activityParams, speciesConcentration, speciesRates, speciesRatesDerivatives ); @@ -271,6 +282,7 @@ class KineticReactions static HPCREACT_HOST_DEVICE void computeReactionRates_impl( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & reactionRates, ARRAY_2D & reactionRatesDerivatives ); @@ -300,6 +312,7 @@ class KineticReactions static HPCREACT_HOST_DEVICE void computeReactionRatesQuotient_impl( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D_SA const & surfaceArea, ARRAY_1D & reactionRates, @@ -326,6 +339,7 @@ class KineticReactions static HPCREACT_HOST_DEVICE void computeSpeciesRates_impl( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & speciesRates, ARRAY_2D & speciesRatesDerivatives ); diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp index 8307036..d8242f3 100644 --- a/src/reactions/reactionsSystems/KineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/KineticReactions_impl.hpp @@ -49,6 +49,7 @@ KineticReactions< REAL_TYPE, LOGE_CONCENTRATION >::computeReactionRates_impl( RealType const &, //temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & reactionRates, ARRAY_2D & reactionRatesDerivatives ) @@ -68,11 +69,16 @@ KineticReactions< REAL_TYPE, { expSpeciesConcentration[i] = exp( speciesConcentration[i] ); } - ACTIVITY_MODEL::calculateActivities( typename ACTIVITY_MODEL::Params(), expSpeciesConcentration, activities ); + ACTIVITY_MODEL::calculateActivities( activityParams, expSpeciesConcentration, activities ); + for( INT_TYPE i=0; i::computeReactionRatesQuotient_impl( RealType const &, //temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D_SA const & surfaceArea, ARRAY_1D & reactionRates, @@ -251,11 +258,11 @@ KineticReactions< REAL_TYPE, { expSpeciesConcentration[i] = exp( speciesConcentration[i] ); } - ACTIVITY_MODEL::calculateActivities( typename ACTIVITY_MODEL::Params(), expSpeciesConcentration, activities ); + ACTIVITY_MODEL::calculateActivities( activityParams, expSpeciesConcentration, activities ); } else { - ACTIVITY_MODEL::calculateActivities( typename ACTIVITY_MODEL::Params(), speciesConcentration, activities ); + ACTIVITY_MODEL::calculateActivities( activityParams, speciesConcentration, activities ); } // loop over each reaction @@ -348,6 +355,7 @@ KineticReactions< REAL_TYPE, LOGE_CONCENTRATION >::computeSpeciesRates_impl( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & speciesRates, ARRAY_2D & speciesRatesDerivatives ) @@ -360,9 +368,12 @@ KineticReactions< REAL_TYPE, HPCREACT_UNUSED_VAR( speciesRatesDerivatives ); } - - - computeReactionRates< PARAMS_DATA >( temperature, params, speciesConcentration, reactionRates, reactionRatesDerivatives ); + computeReactionRates< PARAMS_DATA >( temperature, + params, + activityParams, + speciesConcentration, + reactionRates, + reactionRatesDerivatives ); for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) { diff --git a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp index 46addd3..74697e1 100644 --- a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp @@ -15,8 +15,6 @@ #include "common/pmpl.hpp" #include "common/printers.hpp" #include "common/CArrayWrapper.hpp" -#include "constitutive/activity/Bdot.hpp" -#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" #include "reactions/reactionsSystems/KineticReactions.hpp" #include @@ -58,21 +56,23 @@ struct ComputeReactionRatesTestData template< typename REAL_TYPE, bool LOGE_CONCENTRATION, - 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()] ) + typename ACTIVITY_MODEL, + typename KINETIC_PARAMS_DATA > +void computeReactionRatesTest( KINETIC_PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, + REAL_TYPE const (&initialSpeciesConcentration)[KINETIC_PARAMS_DATA::numSpecies()], + REAL_TYPE const (&surfaceArea)[KINETIC_PARAMS_DATA::numReactions()], + REAL_TYPE const (&expectedReactionRates)[KINETIC_PARAMS_DATA::numReactions()], + REAL_TYPE const (&expectedReactionRatesDerivatives)[KINETIC_PARAMS_DATA::numReactions()][KINETIC_PARAMS_DATA::numSpecies()] ) { using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE, int, int, - Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >, + ACTIVITY_MODEL, LOGE_CONCENTRATION >; - static constexpr int numSpecies = PARAMS_DATA::numSpecies(); - static constexpr int numReactions = PARAMS_DATA::numReactions(); + static constexpr int numSpecies = KINETIC_PARAMS_DATA::numSpecies(); + static constexpr int numReactions = KINETIC_PARAMS_DATA::numReactions(); double const temperature = 298.15; ComputeReactionRatesTestData< numReactions, numSpecies > data; @@ -105,10 +105,11 @@ void computeReactionRatesTest( PARAMS_DATA const & params, } - pmpl::genericKernelWrapper( 1, &data, [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy ) + pmpl::genericKernelWrapper( 1, &data, [params, temperature, activityParams] HPCREACT_DEVICE ( auto * const dataCopy ) { KineticReactionsType::computeReactionRates( temperature, params, + activityParams, dataCopy->speciesConcentration, dataCopy->surfaceArea, dataCopy->reactionRates, @@ -157,20 +158,22 @@ struct ComputeSpeciesRatesTestData template< typename REAL_TYPE, bool LOGE_CONCENTRATION, - typename PARAMS_DATA > -void computeSpeciesRatesTest( PARAMS_DATA const & params, - REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], - REAL_TYPE const (&expectedSpeciesRates)[PARAMS_DATA::numSpecies()], - REAL_TYPE const (&expectedSpeciesRatesDerivatives)[PARAMS_DATA::numSpecies()][PARAMS_DATA::numSpecies()] ) + typename ACTIVITY_MODEL, + typename KINETIC_PARAMS_DATA > +void computeSpeciesRatesTest( KINETIC_PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, + REAL_TYPE const (&initialSpeciesConcentration)[KINETIC_PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedSpeciesRates)[KINETIC_PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedSpeciesRatesDerivatives)[KINETIC_PARAMS_DATA::numSpecies()][KINETIC_PARAMS_DATA::numSpecies()] ) { using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE, int, int, - Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >, + ACTIVITY_MODEL, LOGE_CONCENTRATION >; - static constexpr int numSpecies = PARAMS_DATA::numSpecies(); + static constexpr int numSpecies = KINETIC_PARAMS_DATA::numSpecies(); double const temperature = 298.15; ComputeSpeciesRatesTestData< numSpecies > data; @@ -190,10 +193,11 @@ void computeSpeciesRatesTest( PARAMS_DATA const & params, } } - pmpl::genericKernelWrapper( 1, &data, [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy ) + pmpl::genericKernelWrapper( 1, &data, [params, temperature, activityParams] HPCREACT_DEVICE ( auto * const dataCopy ) { KineticReactionsType::computeSpeciesRates( temperature, params, + activityParams, dataCopy->speciesConcentration, dataCopy->speciesRates, dataCopy->speciesRatesDerivatives ); @@ -236,20 +240,21 @@ struct TimeStepTestData template< typename REAL_TYPE, bool LOGE_CONCENTRATION, - typename PARAMS_DATA > -void timeStepTest( PARAMS_DATA const & params, + typename ACTIVITY_MODEL, + typename KINETIC_PARAMS_DATA > +void timeStepTest( KINETIC_PARAMS_DATA const & params, REAL_TYPE const dt, int const numSteps, - REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], - REAL_TYPE const (&expectedSpeciesConcentrations)[PARAMS_DATA::numSpecies()] ) + REAL_TYPE const (&initialSpeciesConcentration)[KINETIC_PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedSpeciesConcentrations)[KINETIC_PARAMS_DATA::numSpecies()] ) { using KineticReactionsType = reactionsSystems::KineticReactions< REAL_TYPE, int, int, - Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >, + ACTIVITY_MODEL, LOGE_CONCENTRATION >; - static constexpr int numSpecies = PARAMS_DATA::numSpecies(); + static constexpr int numSpecies = KINETIC_PARAMS_DATA::numSpecies(); double const temperature = 298.15; TimeStepTestData< numSpecies > data; From 60498c3bd4b20af69922f6de8ff92c1b787250e5 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Mon, 12 Jan 2026 14:13:14 -0800 Subject: [PATCH 4/5] wip --- src/constitutive/activity/Bdot.hpp | 5 +- src/reactions/exampleSystems/BulkGeneric.hpp | 1 + src/reactions/exampleSystems/ChainGeneric.hpp | 22 ++++ .../unitTests/testGenericChainReactions.cpp | 13 +- src/reactions/geochemistry/Carbonate.hpp | 85 ++++++++++++ .../testGeochemicalKineticReactions.cpp | 124 +++++++++++------- 6 files changed, 195 insertions(+), 55 deletions(-) diff --git a/src/constitutive/activity/Bdot.hpp b/src/constitutive/activity/Bdot.hpp index 57c37da..2ac61e2 100644 --- a/src/constitutive/activity/Bdot.hpp +++ b/src/constitutive/activity/Bdot.hpp @@ -11,6 +11,7 @@ #pragma once #include "DebyeHuckel.hpp" +#include "common/CArrayWrapper.hpp" namespace hpcReact { @@ -28,8 +29,8 @@ class Bdot struct Params : public IONIC_STRENGTH_TYPE::Params { - RealType m_ionSizeParameter[IONIC_STRENGTH_TYPE::Params::numSpecies()]; - RealType m_bdotParameter[IONIC_STRENGTH_TYPE::Params::numSpecies()]; + CArrayWrapper m_ionSizeParameter; + CArrayWrapper m_bdotParameter; }; diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp index 7b157ab..3ea4d13 100644 --- a/src/reactions/exampleSystems/BulkGeneric.hpp +++ b/src/reactions/exampleSystems/BulkGeneric.hpp @@ -72,6 +72,7 @@ simpleKineticTestRateParams = using simpleIonicStrengthType = SpeciatedIonicStrength< double, int, 5 >; using simpleActivityParamsType = Bdot< double, int, simpleIonicStrengthType >::Params; + constexpr simpleActivityParamsType simpleActivityTestParams = diff --git a/src/reactions/exampleSystems/ChainGeneric.hpp b/src/reactions/exampleSystems/ChainGeneric.hpp index 8b0acd8..bf082c0 100644 --- a/src/reactions/exampleSystems/ChainGeneric.hpp +++ b/src/reactions/exampleSystems/ChainGeneric.hpp @@ -12,6 +12,9 @@ #pragma once #include "../reactionsSystems/Parameters.hpp" +#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" +#include "constitutive/activity/Bdot.hpp" +#include "constitutive/activity/Identity.hpp" namespace hpcReact { @@ -63,6 +66,25 @@ namespace ChainGeneric 0 // Use the forward and reverse to calculate the kinetic reaction rates }; + using serialAllKineticIonicStrengthType = SpeciatedIonicStrength< double, int, 3 >; + + using serialAllKineticActivityParamsType = Bdot< double, int, serialAllKineticIonicStrengthType >::Params; + + constexpr + serialAllKineticActivityParamsType + serialAllKineticActivityParams = + { + // species charge + {{ 0.0, 0.0, 0.0 }}, + // ion size parameter + { 3.5, 3.5, 3.5 }, + // bdot parameter + { 0.0, 0.0, 0.0 } + }; + + Identity< double, int, serialAllKineticIonicStrengthType >::Params serialAllKineticIdentityActivityParams = {}; + + // *****UNCRUSTIFY-ON****** } // namespace ChainGeneric } // namespace hpcReact diff --git a/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp index 6e0f9a5..5fded7d 100644 --- a/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp +++ b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp @@ -21,6 +21,9 @@ TEST( testChainGenericKineticReactions, computeReactionRatesTest_chainReactionPa { using namespace hpcReact::ChainGeneric; + using IonicStrengthType = ChainGeneric::serialAllKineticIonicStrengthType; + using ActivityType = Identity< double, int, IonicStrengthType >; + double const initialSpeciesConcentration[] = { 1.0, // C1 @@ -44,12 +47,18 @@ TEST( testChainGenericKineticReactions, computeReactionRatesTest_chainReactionPa { { 0.05, 0.0, 0.0 }, { 0.0, 0.03, 0.0 }, { 0.0, 0.0, 0.02 } }; - computeReactionRatesTest< double, false >( serialAllKineticParams.kineticReactionsParameters(), + computeReactionRatesTest< double, + false, + ActivityType >( serialAllKineticParams.kineticReactionsParameters(), + ChainGeneric::serialAllKineticIdentityActivityParams, initialSpeciesConcentration, surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); - computeReactionRatesTest< double, true >( serialAllKineticParams.kineticReactionsParameters(), + computeReactionRatesTest< double, + true, + ActivityType >( serialAllKineticParams.kineticReactionsParameters(), + ChainGeneric::serialAllKineticIdentityActivityParams, initialSpeciesConcentration, surfaceArea, // No use. Just to pass something here expectedReactionRates, diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp index 52ae3ab..9b4b9f6 100644 --- a/src/reactions/geochemistry/Carbonate.hpp +++ b/src/reactions/geochemistry/Carbonate.hpp @@ -12,6 +12,9 @@ #pragma once #include "../reactionsSystems/Parameters.hpp" +#include "constitutive/ionicStrength/SpeciatedIonicStrength.hpp" +#include "constitutive/activity/Bdot.hpp" +#include "constitutive/activity/Identity.hpp" namespace hpcReact { @@ -107,15 +110,97 @@ constexpr CArrayWrapper mobileSpeciesFlag = 1 // CaCO3 + H+ = Ca+2 + HCO3- }; + + + + +constexpr CArrayWrapper speciesCharge = + { -1.0, // OH- + 0.0, // CO2(aq) + -2.0, // CO3-2 + 1.0, // CaHCO3+ + 0.0, // CaSO4(aq) + 1.0, // CaCl+ + 0.0, // CaCl2(aq) + 0.0, // MgSO4(aq) + -1.0, // NaSO4- + 0.0, // CaCO3(aq) + 1.0, // H+ + -1.0, // HCO3- + 2.0, // Ca+2 + -2.0, // SO4-2 + -1.0, // Cl- + 2.0, // Mg+2 + 1.0 // Na+ + }; + + constexpr CArrayWrapper ionSize = + { + 3.5, // OH- (from H2O = OH- + H+, -gamma 3.5 0.0) + 0.0, // CO2(aq) (neutral, no -gamma; typically gamma ≈ 1) + 5.4, // CO3-2 + 5.4, // CaHCO3+ + 0.0, // CaSO4(aq) (neutral) + 0.0, // CaCl+ (no -gamma in phreeqc.dat) + 0.0, // CaCl2(aq) (neutral) + 0.0, // MgSO4(aq) (neutral) + 0.0, // NaSO4- (no -gamma in phreeqc.dat) + 0.0, // CaCO3(aq) (neutral) + 9.0, // H+ + 5.4, // HCO3- (from CO3-2 + H+ = HCO3-, -gamma 5.4 0.0) + 5.0, // Ca+2 + 5.0, // SO4-2 + 3.5, // Cl- + 5.5, // Mg+2 + 4.0 // Na+ + }; + + constexpr CArrayWrapper bdotParameters = + { + 0.0, // OH- + 0.0, // CO2(aq) + 0.0, // CO3-2 + 0.0, // CaHCO3+ + 0.0, // CaSO4(aq) + 0.0, // CaCl+ + 0.0, // CaCl2(aq) + 0.0, // MgSO4(aq) + 0.0, // NaSO4- + 0.0, // CaCO3(aq) + 0.0, // H+ + 0.0, // HCO3- + 0.165, // Ca+2 + -0.040, // SO4-2 + 0.015, // Cl- + 0.200, // Mg+2 + 0.075 // Na+ + }; + + } using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 10, 0 >; using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 10, 10 >; using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 16, 10, 9 >; +using carbonateIonicStrengthType = SpeciatedIonicStrength< double, int, 17 >; +using carbonateActivityType = Bdot< double, int, carbonateIonicStrengthType >; +using carbonateIdentityActivityType = Identity< double, int, carbonateIonicStrengthType >; + + 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 ); +constexpr carbonateActivityType::Params carbonateActivityParams = +{ + {carbonate::speciesCharge}, + carbonate::ionSize, + carbonate::bdotParameters +}; + + + +Identity< double, int, carbonateIonicStrengthType >::Params carbonateIdentityActivityParams = {}; // *****UNCRUSTIFY-ON****** } // namespace geochemistry diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp index 8c4e116..8235611 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp @@ -81,12 +81,21 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic ) }; - computeReactionRatesTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(), + using ActivityType = carbonateIdentityActivityType; + + + computeReactionRatesTest< double, + false, + ActivityType >( carbonateSystemAllKinetic.kineticReactionsParameters(), + typename ActivityType::Params(), initialSpeciesConcentration, surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); - computeReactionRatesTest< double, true >( carbonateSystemAllKinetic.kineticReactionsParameters(), + computeReactionRatesTest< double, + true, + ActivityType >( carbonateSystemAllKinetic.kineticReactionsParameters(), + typename ActivityType::Params(), initialSpeciesConcentration, surfaceArea, // No use. Just to pass something here expectedReactionRates, @@ -124,12 +133,20 @@ TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem ) { 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(), + using ActivityType = carbonateIdentityActivityType; + + computeReactionRatesTest< double, + false, + ActivityType >( carbonateSystem.kineticReactionsParameters(), + typename ActivityType::Params(), initialSpeciesConcentration, surfaceArea, expectedReactionRates, expectedReactionRatesDerivatives ); - computeReactionRatesTest< double, true >( carbonateSystem.kineticReactionsParameters(), + computeReactionRatesTest< double, + true, + ActivityType >( carbonateSystem.kineticReactionsParameters(), + typename ActivityType::Params(), initialSpeciesConcentration, surfaceArea, expectedReactionRates, @@ -175,54 +192,59 @@ TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem ) // } -TEST( testKineticReactions, testTimeStep_carbonateSystemAllKinetic ) -{ - double const initialSpeciesConcentration[17] = - { - 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- - 1.0e-16, // CaCO3 - 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 - }; +// TEST( testKineticReactions, testTimeStep_carbonateSystemAllKinetic ) +// { +// double const initialSpeciesConcentration[17] = +// { +// 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- +// 1.0e-16, // CaCO3 +// 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 expectedSpeciesConcentrations[17] = - { 2.327841695586879e-11, // OH- - 0.37555955033916549, // CO2 - 3.956656978189456e-11, // CO3-2 - 6.739226982791492e-05, // CaHCO3+ - 5.298329882666738e-03, // CaSO4 - 5.844517547638333e-03, // CaCl+ - 1.277319392670652e-02, // CaCl2 - 6.618125707964991e-03, // MgSO4 - 1.769217213462983e-02, // NaSO4- - 1.065032288527957e-09, // CaCO3 - 4.396954721488358e-04, // H+ - 3.723009698453808e-04, // HCO3- - 1.471656530812871e-02, // Ca+2 - 2.491372274738741e-03, // SO4-2 - 1.858609094598949e+00, // Cl- - 9.881874292035110e-03, // Mg+2 - 1.072307827865370e+00 // Na+1 - }; +// double const expectedSpeciesConcentrations[17] = +// { 2.327841695586879e-11, // OH- +// 0.37555955033916549, // CO2 +// 3.956656978189456e-11, // CO3-2 +// 6.739226982791492e-05, // CaHCO3+ +// 5.298329882666738e-03, // CaSO4 +// 5.844517547638333e-03, // CaCl+ +// 1.277319392670652e-02, // CaCl2 +// 6.618125707964991e-03, // MgSO4 +// 1.769217213462983e-02, // NaSO4- +// 1.065032288527957e-09, // CaCO3 +// 4.396954721488358e-04, // H+ +// 3.723009698453808e-04, // HCO3- +// 1.471656530812871e-02, // Ca+2 +// 2.491372274738741e-03, // SO4-2 +// 1.858609094598949e+00, // Cl- +// 9.881874292035110e-03, // Mg+2 +// 1.072307827865370e+00 // Na+1 +// }; - timeStepTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(), - 10.0, - 10000, - initialSpeciesConcentration, - expectedSpeciesConcentrations ); +// using ActivityType = carbonateIdentityActivityType; + +// timeStepTest< double, +// false, +// ActivityType >( carbonateSystemAllKinetic.kineticReactionsParameters(), +// ActivityType::Params(), +// 10.0, +// 10000, +// initialSpeciesConcentration, +// expectedSpeciesConcentrations ); // ln(c) as the primary variable results in a singular system. // timeStepTest< double, true >( simpleKineticTestRateParams, @@ -230,7 +252,7 @@ TEST( testKineticReactions, testTimeStep_carbonateSystemAllKinetic ) // 10, // initialSpeciesConcentration, // expectedSpeciesConcentrations ); -} +//} int main( int argc, char * * argv ) { From 73861a1560bd8ecb316cd5f0a2bee3714a1de7a8 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Tue, 13 Jan 2026 09:01:18 -0800 Subject: [PATCH 5/5] change parameter specification approach --- src/reactions/exampleSystems/BulkGeneric.hpp | 89 +++++++----- src/reactions/exampleSystems/ChainGeneric.hpp | 91 ++++++------ .../exampleSystems/MoMasBenchmark.hpp | 132 ++++++++++-------- .../MixedEquilibriumKineticReactions.hpp | 9 ++ 4 files changed, 190 insertions(+), 131 deletions(-) diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp index 3ea4d13..2550b1d 100644 --- a/src/reactions/exampleSystems/BulkGeneric.hpp +++ b/src/reactions/exampleSystems/BulkGeneric.hpp @@ -50,39 +50,52 @@ namespace bulkGeneric using simpleKineticParamsType = reactionsSystems::KineticReactionsParameters< double, int, int, 5, 2 >; -constexpr -simpleKineticParamsType -simpleKineticTestRateParams = -{ +constexpr CArrayWrapper< double, 2, 5 > simpleKineticStoichMatrix = +{ // stoichiometric matrix { { -2, 1, 1, 0, 0 }, { 0, 0, -1, -1, 2 } - }, - // forward rate constants - { 1.0, 0.5 }, - // reverse rate constants - { 1.0, 0.5 }, - // equilibrium constants - { 1.0, 1.0 }, - // Use the forward and reverse to calculate the kinetic reaction rates - 0 + } }; +constexpr CArrayWrapper< double, 2 > simpleKineticForwardRates = +{ 1.0, 0.5 }; + +constexpr CArrayWrapper< double, 2 > simpleKineticReverseRates = +{ 1.0, 0.5 }; + +constexpr CArrayWrapper< double, 2 > simpleKineticEquilibriumConstants = +{ 1.0, 1.0 }; + +constexpr simpleKineticParamsType simpleKineticTestRateParams( + simpleKineticStoichMatrix, + simpleKineticForwardRates, + simpleKineticReverseRates, + simpleKineticEquilibriumConstants, + 0 ); + using simpleIonicStrengthType = SpeciatedIonicStrength< double, int, 5 >; using simpleActivityParamsType = Bdot< double, int, simpleIonicStrengthType >::Params; -constexpr -simpleActivityParamsType -simpleActivityTestParams = +constexpr CArrayWrapper< double, 5 > simpleSpeciesCharge = +{ 2.0, -1.0, 1.0, 0.0, -1.0 }; + +constexpr CArrayWrapper< double, 5 > simpleIonSize = +{ 4.0, 3.5, 3.5, 3.5, 3.5 }; + +constexpr CArrayWrapper< double, 5 > simpleBdotParameters = +{ 0.1, 0.1, 0.1, 0.0, 0.1 }; + +constexpr simpleActivityParamsType simpleActivityTestParams = { // species charge - {{ 2.0, -1.0, 1.0, 0.0, -1.0 }}, + {{ simpleSpeciesCharge }}, // ion size parameter - { 4.0, 3.5, 3.5, 3.5, 3.5 }, + simpleIonSize, // bdot parameter - { 0.1, 0.1, 0.1, 0.0, 0.1 } + simpleBdotParameters }; Identity< double, int, simpleIonicStrengthType >::Params simpleIdentityActivityTestParams = {}; @@ -90,27 +103,35 @@ Identity< double, int, simpleIonicStrengthType >::Params simpleIdentityActivityT using simpleMixedParamsType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >; -constexpr -simpleMixedParamsType -simpleTestRateParams = -{ +constexpr CArrayWrapper< double, 2, 5 > simpleMixedStoichMatrix = +{ // stoichiometric matrix { { -2, 1, 1, 0, 0 }, { 0, 0, -1, -1, 2 } - }, - // equilibrium constants - { 1.0, 1.0 }, - // forward rate constants - { 1.0, 0.5 }, - // reverse rate constants - { 1.0, 0.5 }, - // flag of mobile secondary species - { 1, 1 }, - // Use the forward and reverse to calculate the kinetic reaction rates - 0 + } }; +constexpr CArrayWrapper< double, 2 > simpleMixedEquilibriumConstants = +{ 1.0, 1.0 }; + +constexpr CArrayWrapper< double, 2 > simpleMixedForwardRates = +{ 1.0, 0.5 }; + +constexpr CArrayWrapper< double, 2 > simpleMixedReverseRates = +{ 1.0, 0.5 }; + +constexpr CArrayWrapper< int, 2 > simpleMixedMobileSpeciesFlag = +{ 1, 1 }; + +constexpr simpleMixedParamsType simpleTestRateParams( + simpleMixedStoichMatrix, + simpleMixedEquilibriumConstants, + simpleMixedForwardRates, + simpleMixedReverseRates, + simpleMixedMobileSpeciesFlag, + 0 ); + // *****UNCRUSTIFY-ON****** } // namespace bulkGeneric } // namespace hpcReact diff --git a/src/reactions/exampleSystems/ChainGeneric.hpp b/src/reactions/exampleSystems/ChainGeneric.hpp index bf082c0..773ec56 100644 --- a/src/reactions/exampleSystems/ChainGeneric.hpp +++ b/src/reactions/exampleSystems/ChainGeneric.hpp @@ -24,62 +24,75 @@ namespace ChainGeneric using serialAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 3, 3, 0 >; - constexpr serialAllKineticType serialAllKineticParams = + constexpr CArrayWrapper< double, 3, 3 > serialAllKineticStoichMatrix = { // Stoichiometric matrix [3 rows × 3 columns] // Columns 0–3 are primary species: {C1, C2, C3 } { - // C1 C2 C3 + // C1 C2 C3 { -1, 1, 0 }, // C1 = C2 { 0, -1, 1 }, // C2 = C3 - { 0, 0, -1 }, // 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 - }, - - 0 // Use the forward and reverse to calculate the kinetic reaction rates + constexpr CArrayWrapper< double, 3 > serialAllKineticEquilibriumConstants = + { + 0, // C1 = C2 + 0, // C2 = C3 + 0 // C3 + }; + + constexpr CArrayWrapper< double, 3 > serialAllKineticForwardRates = + { + 0.05, // C1 = C2 + 0.03, // C2 = C3 + 0.02 // C3 + }; + + constexpr CArrayWrapper< double, 3 > serialAllKineticReverseRates = + { + 0.0, // C1 = C2 + 0.0, // C2 = C3 + 0.0 // C3 + }; + + constexpr CArrayWrapper< int, 3 > serialAllKineticMobileSpeciesFlag = + { + 1, // C1 = C2 + 1, // C2 = C3 + 1 // C3 }; + constexpr serialAllKineticType serialAllKineticParams( + serialAllKineticStoichMatrix, + serialAllKineticEquilibriumConstants, + serialAllKineticForwardRates, + serialAllKineticReverseRates, + serialAllKineticMobileSpeciesFlag, + 0 ); + using serialAllKineticIonicStrengthType = SpeciatedIonicStrength< double, int, 3 >; using serialAllKineticActivityParamsType = Bdot< double, int, serialAllKineticIonicStrengthType >::Params; - constexpr - serialAllKineticActivityParamsType - serialAllKineticActivityParams = + constexpr CArrayWrapper< double, 3 > serialAllKineticSpeciesCharge = + { 0.0, 0.0, 0.0 }; + + constexpr CArrayWrapper< double, 3 > serialAllKineticIonSize = + { 3.5, 3.5, 3.5 }; + + constexpr CArrayWrapper< double, 3 > serialAllKineticBdotParameters = + { 0.0, 0.0, 0.0 }; + + constexpr serialAllKineticActivityParamsType serialAllKineticActivityParams = { // species charge - {{ 0.0, 0.0, 0.0 }}, + {{ serialAllKineticSpeciesCharge }}, // ion size parameter - { 3.5, 3.5, 3.5 }, + serialAllKineticIonSize, // bdot parameter - { 0.0, 0.0, 0.0 } + serialAllKineticBdotParameters }; Identity< double, int, serialAllKineticIonicStrengthType >::Params serialAllKineticIdentityActivityParams = {}; diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index c0c8f06..cef63e2 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -22,23 +22,24 @@ namespace MoMasBenchmark using easyCaseType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >; using mediumCaseType = reactionsSystems::MixedReactionsParameters< double, int, int, 14, 10, 9 >; - constexpr easyCaseType easyCaseParams = -{ - // Stoichiometric matrix [7 rows × 12 columns] - // Columns 0–6 are secondary species (must be -1 on the diagonal) - // Columns 7–11 are primary species: {X1, X2, X3, X4, S} + constexpr CArrayWrapper< double, 7, 12 > easyCaseStoichMatrix = { - // 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, 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 - { 0, 0, 0, 0, 0, 0, -1, 0, -3, 0, 1, 2 } // CS2 = -3X2 + X4 + 2S - }, - - // Equilibrium constants K + // Stoichiometric matrix [7 rows × 12 columns] + // Columns 0–6 are secondary species (must be -1 on the diagonal) + // Columns 7–11 are primary species: {X1, X2, X3, X4, S} + { + // 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, 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 + { 0, 0, 0, 0, 0, 0, -1, 0, -3, 0, 1, 2 } // CS2 = -3X2 + X4 + 2S + } + }; + + constexpr CArrayWrapper< double, 7 > easyCaseEquilibriumConstants = { 1.0e12, // C1 + X2 = inf 1.0, // C2 = X2 + X3 @@ -47,10 +48,10 @@ namespace MoMasBenchmark 1.0e-35, // C5 = 4X2 + 3X3 + X4 1.0e-6, // CS1 = 3X2 + X3 + S 1.0e1 // CS2 + 3X2 = + X4 + 2S - }, + }; - // Forward rate constants - { + constexpr CArrayWrapper< double, 7 > easyCaseForwardRates = + { 0.0, // C1 = -X2 0.0, // C2 = X2 + X3 0.0, // C3 = -X2 + X4 @@ -58,10 +59,10 @@ namespace MoMasBenchmark 0.0, // C5 = 4X2 + 3X3 + X4 0.0, // CS1 = 3X2 + X3 + S 0.0 // CS2 = -3X2 + X4 + 2S - }, + }; - // Reverse rate constants - { + constexpr CArrayWrapper< double, 7 > easyCaseReverseRates = + { 0.0, // C1 = -X2 0.0, // C2 = X2 + X3 0.0, // C3 = -X2 + X4 @@ -69,39 +70,47 @@ namespace MoMasBenchmark 0.0, // C5 = 4X2 + 3X3 + X4 0.0, // CS1 = 3X2 + X3 + S 0.0 // CS2 = -3X2 + X4 + 2S - }, + }; - // Flag of mobile secondary species - { 1, // C1 = -X2 + constexpr CArrayWrapper< int, 7 > easyCaseMobileSpeciesFlag = + { + 1, // C1 = -X2 1, // C2 = X2 + X3 1, // C3 = -X2 + X4 1, // C4 = -4X2 + X3 + 3X4 1, // C5 = 4X2 + 3X3 + X4 0, // CS1 = 3X2 + X3 + S 0 // CS2 = -3X2 + X4 + 2S - } -}; + }; -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} + constexpr easyCaseType easyCaseParams( + easyCaseStoichMatrix, + easyCaseEquilibriumConstants, + easyCaseForwardRates, + easyCaseReverseRates, + easyCaseMobileSpeciesFlag ); + + constexpr CArrayWrapper< double, 10, 14 > mediumCaseStoichMatrix = { - // 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 + // 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) + } + }; + + constexpr CArrayWrapper< double, 10 > mediumCaseEquilibriumConstants = { 1.0e12, // C1 + X2 = inf 1.0, // C2 = X2 + X3 @@ -113,10 +122,10 @@ constexpr mediumCaseType mediumCaseParams = 1.0e-6, // CS1 = 3X2 + X3 + S 1.0e1, // CS2 + 3X2 = X4 + 2S 5 // Cc + 3X2 = X4 (kinetic) - }, + }; - // Forward rate constants - { + constexpr CArrayWrapper< double, 10 > mediumCaseForwardRates = + { 0.0, // C1 = -X2 0.0, // C2 = X2 + X3 0.0, // C3 = -X2 + X4 @@ -126,11 +135,11 @@ constexpr mediumCaseType mediumCaseParams = 0.0, // C7 = -8X2 + 2X4 0.0, // CS1 = 3X2 + X3 + S 0.0, // CS2 = -3X2 + X4 + 2S - 10.0 // Cc = -3X2 + X4 (kinetic) - }, + 10.0 // Cc = -3X2 + X4 (kinetic) + }; - // Reverse rate constants - { + constexpr CArrayWrapper< double, 10 > mediumCaseReverseRates = + { 0.0, // C1 = -X2 0.0, // C2 = X2 + X3 0.0, // C3 = -X2 + X4 @@ -141,10 +150,11 @@ constexpr mediumCaseType mediumCaseParams = 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 + constexpr CArrayWrapper< int, 10 > mediumCaseMobileSpeciesFlag = + { + 1, // C1 = -X2 1, // C2 = X2 + X3 1, // C3 = -X2 + X4 1, // C4 = -4X2 + X3 + 3X4 @@ -154,8 +164,14 @@ constexpr mediumCaseType mediumCaseParams = 0, // CS1 = 3X2 + X3 + S 0, // CS2 = -3X2 + X4 + 2S 1 // Cc = -3X2 + X4 (kinetic) - } -}; + }; + + constexpr mediumCaseType mediumCaseParams( + mediumCaseStoichMatrix, + mediumCaseEquilibriumConstants, + mediumCaseForwardRates, + mediumCaseReverseRates, + mediumCaseMobileSpeciesFlag ); // *****UNCRUSTIFY-ON****** } // namespace MoMasBenchmark diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp index 0950cb6..58f18c0 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp @@ -94,6 +94,7 @@ class MixedEquilibriumKineticReactions static HPCREACT_HOST_DEVICE inline void updateMixedSystem( RealType const & temperature, PARAMS_DATA const & params, + ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, ARRAY_1D_TO_CONST_KINETIC const & surfaceArea, ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, @@ -108,6 +109,7 @@ class MixedEquilibriumKineticReactions { updateMixedSystem_impl( temperature, params, + activityParams, logPrimarySpeciesConcentrations, surfaceArea, logSecondarySpeciesConcentrations, @@ -147,6 +149,7 @@ class MixedEquilibriumKineticReactions static HPCREACT_HOST_DEVICE inline void computeReactionRates( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, ARRAY_1D_TO_CONST2 const & logSecondarySpeciesConcentrations, ARRAY_1D_TO_CONST_KINETIC const & surfaceArea, @@ -156,6 +159,7 @@ class MixedEquilibriumKineticReactions { computeReactionRates_impl( temperature, params, + activityParams, logPrimarySpeciesConcentrations, logSecondarySpeciesConcentrations, surfaceArea, @@ -188,6 +192,7 @@ class MixedEquilibriumKineticReactions typename ARRAY_2D > static HPCREACT_HOST_DEVICE inline void computeAggregateSpeciesRates( PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D_TO_CONST2 const & reactionRates, ARRAY_2D_TO_CONST const & reactionRatesDerivatives, @@ -201,6 +206,7 @@ class MixedEquilibriumKineticReactions ARRAY_1D, ARRAY_2D, true >( params, + activityParams, speciesConcentration, reactionRates, reactionRatesDerivatives, @@ -248,6 +254,7 @@ class MixedEquilibriumKineticReactions static HPCREACT_HOST_DEVICE void updateMixedSystem_impl( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, ARRAY_1D_TO_CONST_KINETIC const & surfaceArea, ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, @@ -279,6 +286,7 @@ class MixedEquilibriumKineticReactions static HPCREACT_HOST_DEVICE void computeReactionRates_impl( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, ARRAY_1D_TO_CONST2 const & logSecondarySpeciesConcentrations, ARRAY_1D_TO_CONST_KINETIC const & surfaceArea, @@ -312,6 +320,7 @@ class MixedEquilibriumKineticReactions bool CALCULATE_DERIVATIVES > static HPCREACT_HOST_DEVICE void computeAggregateSpeciesRates_impl( PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D_TO_CONST2 const & reactionRates, ARRAY_2D_TO_CONST const & reactionRatesDerivatives,