From 6c89cc570fa725328ae802ccc585771feab6f015 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Wed, 26 Nov 2025 17:17:29 -0800 Subject: [PATCH 1/3] 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/3] 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/3] 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;