From 6c89cc570fa725328ae802ccc585771feab6f015 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Wed, 26 Nov 2025 17:17:29 -0800 Subject: [PATCH 1/8] 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/8] 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/8] 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/8] 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/8] 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, From 3d89882708e9d8172fc65794c75b9e87c05d10b3 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Thu, 19 Feb 2026 14:31:29 -0800 Subject: [PATCH 6/8] intermediate commit --- src/constitutive/activity/Bdot.hpp | 25 ++- src/constitutive/activity/DebyeHuckel.hpp | 11 +- src/constitutive/activity/Identity.hpp | 8 +- .../ionicStrength/SpeciatedIonicStrength.hpp | 10 +- src/constitutive/unitTests/testBdot.cpp | 5 +- .../unitTests/testIonicStrength.cpp | 10 +- src/reactions/massActions/MassActions.hpp | 62 +++--- ...ionsAggregatePrimaryConcentration_impl.hpp | 3 + .../reactionsSystems/KineticReactions.hpp | 180 ++++++++++++++---- .../KineticReactions_impl.hpp | 130 +++++++------ .../MixedEquilibriumKineticReactions.hpp | 17 +- .../MixedEquilibriumKineticReactions_impl.hpp | 5 +- .../mixedReactionsTestUtilities.hpp | 4 +- 13 files changed, 297 insertions(+), 173 deletions(-) diff --git a/src/constitutive/activity/Bdot.hpp b/src/constitutive/activity/Bdot.hpp index 2ac61e2..e34addc 100644 --- a/src/constitutive/activity/Bdot.hpp +++ b/src/constitutive/activity/Bdot.hpp @@ -36,15 +36,20 @@ struct Params : public IONIC_STRENGTH_TYPE::Params template< typename ARRAY_1D_TO_CONST, - typename ARRAY_1D > + typename ARRAY_1D, + typename ARRAY_2D > static inline HPCREACT_HOST_DEVICE void calculateActivities( Params const & params, ARRAY_1D_TO_CONST const & speciesConcentrations, - ARRAY_1D & activities ) + ARRAY_1D & activities, + ARRAY_2D & dActivities_dConcentrations ) { - RealType const ionicStrength = IONIC_STRENGTH_TYPE::calculate( params, speciesConcentrations ); + RealType dIonicStrength_dConcentration[ Params::numSpecies() ]; + RealType const ionicStrength = IONIC_STRENGTH_TYPE::calculate( params, + speciesConcentrations, + dIonicStrength_dConcentration ); RealType const sqrtI = sqrt(ionicStrength); RealType const rho_w = 997.0479; // kg/m3 RealType const eps_r = 78.54; // dimensionless @@ -60,12 +65,22 @@ calculateActivities( Params const & params, 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; + B_gamma, + dlog10_gamma_dI ); + RealType const log10gamma = DebyeHuckel_term + b[i] * ionicStrength; + RealType const gamma = pow(10.0, log10gamma); + + activities[i] = speciesConcentrations[i] * gamma; + dActivities_dConcentrations[i][i] = gamma; + for( IndexType j=0; j static inline HPCREACT_HOST_DEVICE void calculateActivities( PARAMS const & , ARRAY_1D_TO_CONST const & speciesConcentrations, - ARRAY_1D & activities ) + ARRAY_1D & activities, + ARRAY_2D & dActivities_dConcentrations ) { constexpr IndexType numSpecies = PARAMS::numSpecies(); for( IndexType i=0; i +template< typename ARRAY_1D_TO_CONST, + typename ARRAY_1D > static inline HPCREACT_HOST_DEVICE REAL_TYPE calculate( Params const & params, - ARRAY_1D_TO_CONST const & speciesConcentration ) + ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D & dIonicStrength_dConcentration ) { REAL_TYPE I = 0.0; auto const & numSpecies = params.numSpecies(); auto const & speciesCharge = params.m_speciesCharge; for( int i=0; i::Params testParams TEST( testBdot, testIonicStrength ) { double speciesConcentration[ testParams.numSpecies() ] = { 0.1, 0.2, 0.3 }; + double dIonicStrength_dConcentration[ testParams.numSpecies() ]; double I = SpeciatedIonicStrength< double, int, 3 >::calculate( testParams, - speciesConcentration ); - + speciesConcentration, + dIonicStrength_dConcentration ); 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 ); diff --git a/src/constitutive/unitTests/testIonicStrength.cpp b/src/constitutive/unitTests/testIonicStrength.cpp index c29ba04..db13320 100644 --- a/src/constitutive/unitTests/testIonicStrength.cpp +++ b/src/constitutive/unitTests/testIonicStrength.cpp @@ -24,17 +24,19 @@ using IonicStrength = SpeciatedIonicStrength; constexpr IonicStrength::Params testParams { // Species charge - { 1.0, -1.0, 2.0 } + { -1.0, -1.0, 2.0 } }; TEST( testBdot, testIonicStrength ) { double speciesConcentration[ testParams.numSpecies() ] = { 0.1, 0.2, 0.3 }; - + double dIonicStrength_dConcentration[ testParams.numSpecies() ]; + double I = IonicStrength::calculate( testParams, - speciesConcentration ); + speciesConcentration, + dIonicStrength_dConcentration ); - double expectedI = 0.5 * ( 0.1 * 1.0 * 1.0 + 0.2 * (-1.0) * (-1.0) + 0.3 * 2.0 * 2.0 ); + 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 ); } diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index 17dbbf8..c77d86b 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -33,9 +33,9 @@ template< typename REAL_TYPE, typename FUNC > HPCREACT_HOST_DEVICE inline -void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, - ARRAY_1D & logSecondarySpeciesConcentrations, +void calculateLogSecondaryActivities( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimaryActivities, + ARRAY_1D & logSecondaryActivities, FUNC && derivativeFunc ) { static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); @@ -43,16 +43,15 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, for( INDEX_TYPE i = 0; i < numSecondarySpecies; ++i ) { - logSecondarySpeciesConcentrations[i] = 0.0; + logSecondaryActivities[i] = 0.0; } for( int j=0; j HPCREACT_HOST_DEVICE inline -void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, - ARRAY_1D & logSecondarySpeciesConcentrations ) +void calculateLogSecondaryActivities( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimaryActivities, + ARRAY_1D & logSecondaryActivities ) { if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) { return; } - massActions_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE, + massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations, + logPrimaryActivities, + logSecondaryActivities, [](INDEX_TYPE, INDEX_TYPE, REAL_TYPE ){} ); } @@ -95,17 +94,17 @@ template< typename REAL_TYPE, typename ARRAY_2D > HPCREACT_HOST_DEVICE inline -void calculateLogSecondarySpeciesConcentrationWrtLogC( PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, - ARRAY_1D & logSecondarySpeciesConcentrations, - ARRAY_2D & dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ) +void calculateLogSecondaryActivitiesWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimaryActivities, + ARRAY_1D & logSecondaryActivities, + ARRAY_2D & dLogSecondaryActivities_dLogPrimaryActivities ) { - massActions_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations, + massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, + logPrimaryActivities, + logSecondaryActivities, [&]( const int j, const int k, REAL_TYPE const value ) { - dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[j][k] = value; + dLogSecondaryActivities_dLogPrimaryActivities[j][k] = value; } ); } @@ -128,12 +127,11 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); - - calculateLogSecondarySpeciesConcentration< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations ); + calculateLogSecondaryActivities< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations ); for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) { for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) @@ -154,8 +152,7 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, for( int k=0; k( params, logPrimarySpeciesConcentration, aggregatePrimaryConcentrations, diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp index b53aad1..1a1b6ea 100644 --- a/src/reactions/reactionsSystems/KineticReactions.hpp +++ b/src/reactions/reactionsSystems/KineticReactions.hpp @@ -12,6 +12,7 @@ #pragma once #include "common/macros.hpp" +#include "constitutive/activity/activity.hpp" #include @@ -54,7 +55,7 @@ class KineticReactions using IndexType = INDEX_TYPE; /** - * @copydoc KineticReactions::computeReactionRates_impl() + * @copydoc KineticReactions::computeReactionRates() */ template< typename PARAMS_DATA, typename ARRAY_1D_TO_CONST, @@ -66,14 +67,36 @@ class KineticReactions typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & reactionRates, - ARRAY_2D & reactionRatesDerivatives ) + ARRAY_2D & dReactionRates_dConcentration ) { - computeReactionRates_impl< PARAMS_DATA, true >( temperature, + + RealType activities[ PARAMS_DATA::numSpecies() ]; + RealType dActivities_dConcentration[ PARAMS_DATA::numSpecies() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; + RealType dReactionRates_dActivities[ PARAMS_DATA::numReactions() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; + + calculateActivities( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); + + computeReactionRates< PARAMS_DATA, true >( temperature, params, - activityParams, - speciesConcentration, + activities, reactionRates, - reactionRatesDerivatives ); + dReactionRates_dActivities ); + + // chain rule to get dReactionRate_dConcentration + for( IntType r=0; r( temperature, + + + RealType activities[ PARAMS_DATA::numSpecies() ]; + RealType dActivities_dConcentration[ PARAMS_DATA::numSpecies() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; + RealType dReactionRates_dActivities[ PARAMS_DATA::numReactions() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; + + calculateActivities( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); + + computeReactionRates< PARAMS_DATA, false >( temperature, params, activityParams, - speciesConcentration, + activities, reactionRates, - reactionRatesDerivatives ); + dReactionRates_dActivities ); + HPCREACT_UNUSED_VAR( dReactionRates_dActivities ); } /** @@ -117,7 +151,7 @@ class KineticReactions * @param speciesConcentration The array of species concentrations. * @param surfaceArea The array of surface area. * @param reactionRates The array of reaction rates. - * @param reactionRatesDerivatives The array of reaction rates derivatives. + * @param dReactionRates_dConcentration The array of reaction rates derivatives. * @details * This function computes the reaction rates for a given set of reactions, * taking into account the surface area of the reactions. If @@ -136,26 +170,51 @@ class KineticReactions ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D_SA const & surfaceArea, ARRAY_1D & reactionRates, - ARRAY_2D & reactionRatesDerivatives ) + ARRAY_2D & dReactionRates_dConcentration ) { + + RealType activities[ PARAMS_DATA::numSpecies() ]; + RealType dActivities_dConcentration[ PARAMS_DATA::numSpecies() ][ PARAMS_DATA::numSpecies() ]{ }; + RealType dReactionRates_dActivities[ PARAMS_DATA::numReactions() ][ PARAMS_DATA::numSpecies() ]{ }; + + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + LOGE_CONCENTRATION >( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); + if( params.reactionRatesUpdateOption() == 0 ) { - computeReactionRates_impl< PARAMS_DATA, true >( temperature, + computeReactionRates< PARAMS_DATA, true >( temperature, params, - activityParams, - speciesConcentration, + activities, reactionRates, - reactionRatesDerivatives ); + dReactionRates_dActivities ); } else if( params.reactionRatesUpdateOption() == 1 ) { computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature, params, - activityParams, - speciesConcentration, + activities, surfaceArea, reactionRates, - reactionRatesDerivatives ); + dReactionRates_dActivities ); + } + + // chain rule to get dReactionRate_dConcentration + for( IntType r=0; r( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); + computeSpeciesRates_impl< PARAMS_DATA, true >( temperature, params, - activityParams, - speciesConcentration, + activities, speciesRates, - speciesRatesDerivatives ); + dSpeciesRates_dActivities ); + // chain rule to get dSpeciesRates_dConcentration + for( IntType i=0; i( temperature, params, - activityParams, - speciesConcentration, + activities, speciesRates, - speciesRatesDerivatives ); + void() ); } /** @@ -243,6 +333,17 @@ class KineticReactions private: + // template< typename ARRAY_1D_TO_CONST, + // typename ARRAY_1D, + // typename ARRAY_2D > + // static HPCREACT_HOST_DEVICE void + // calculateActivities( typename ACTIVITY_MODEL::Params const & params, + // ARRAY_1D_TO_CONST const & speciesConcentration, + // ARRAY_1D & activities, + // ARRAY_2D & dActivities_dConcentration ); + + + /** * @brief Compute the reaction rates for a given set of species concentrations. * @tparam PARAMS_DATA The type of the parameters data. @@ -280,12 +381,11 @@ class KineticReactions typename ARRAY_1D, typename ARRAY_2D > static HPCREACT_HOST_DEVICE void - computeReactionRates_impl( RealType const & temperature, + computeReactionRates( RealType const & temperature, PARAMS_DATA const & params, - typename ACTIVITY_MODEL::Params const & activityParams, - ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D_TO_CONST const & activities, ARRAY_1D & reactionRates, - ARRAY_2D & reactionRatesDerivatives ); + ARRAY_2D & dReactionRates_dActivities ); /** * @brief @@ -297,10 +397,10 @@ class KineticReactions * @tparam ARRAY_2D * @param temperature * @param params - * @param speciesConcentration + * @param activities * @param surfaceArea * @param reactionRates - * @param reactionRatesDerivatives + * @param dReactionRates_dActivities * @return HPCREACT_HOST_DEVICE */ template< typename PARAMS_DATA, @@ -312,11 +412,10 @@ 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_TO_CONST const & activities, ARRAY_1D_SA const & surfaceArea, ARRAY_1D & reactionRates, - ARRAY_2D & reactionRatesDerivatives ); + ARRAY_2D & dReactionRates_dActivities ); /** * @brief Compute the kinetic species rates for a given set of kinetic reactions. @@ -327,9 +426,9 @@ class KineticReactions * @tparam ARRAY_2D The type of the array of species rates derivatives. * @param temperature The temperature of the system. * @param params The parameters data. - * @param speciesConcentration The array of species concentrations. + * @param activities The array of activities. * @param speciesRates The array of species rates. - * @param speciesRatesDerivatives The array of species rates derivatives. + * @param dReactionRates_dActivities The array of species rates derivatives. */ template< typename PARAMS_DATA, bool CALCULATE_DERIVATIVES, @@ -339,10 +438,9 @@ 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_TO_CONST const & activities, ARRAY_1D & speciesRates, - ARRAY_2D & speciesRatesDerivatives ); + ARRAY_2D & dReactionRates_dActivities ); }; diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp index d8242f3..5883563 100644 --- a/src/reactions/reactionsSystems/KineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/KineticReactions_impl.hpp @@ -14,6 +14,7 @@ #include "common/constants.hpp" #include "common/CArrayWrapper.hpp" #include "common/DirectSystemSolve.hpp" +#include "constitutive/activity/activity.hpp" #include #include @@ -31,6 +32,45 @@ namespace reactionsSystems { +// template< typename REAL_TYPE, +// typename INT_TYPE, +// typename INDEX_TYPE, +// typename ACTIVITY_MODEL, +// bool LOGE_CONCENTRATION > +// template< typename ARRAY_1D_TO_CONST, +// typename ARRAY_1D, +// typename ARRAY_2D > +// HPCREACT_HOST_DEVICE inline void +// KineticReactions< REAL_TYPE, +// INT_TYPE, +// INDEX_TYPE, +// ACTIVITY_MODEL, +// LOGE_CONCENTRATION +// >::calculateActivities( typename ACTIVITY_MODEL::Params const & activityParams, +// ARRAY_1D_TO_CONST const & speciesConcentration, +// ARRAY_1D & activities, +// ARRAY_2D & dActivities_dConcentration ) +// { +// if constexpr( LOGE_CONCENTRATION ) +// { +// RealType expSpeciesConcentration[ ACTIVITY_MODEL::Params::numSpecies() ]; +// for( INT_TYPE i=0; i::computeReactionRates_impl( RealType const &, //temperature, + >::computeReactionRates( RealType const &, //temperature, PARAMS_DATA const & params, - typename ACTIVITY_MODEL::Params const & activityParams, - ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D_TO_CONST const & activities, ARRAY_1D & reactionRates, - ARRAY_2D & reactionRatesDerivatives ) + ARRAY_2D & dReactionRate_dActivities ) { if constexpr( !CALCULATE_DERIVATIVES ) { - 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 ) { - reactionRatesDerivatives( r, i ) = -reverseRateConstant * exp( productConcReverse ) * s_ri; + dReactionRate_dActivities[ r ][ i ] = -reverseRateConstant * exp( productConcReverse ) * s_ri; } else { - reactionRatesDerivatives( r, i ) = 0.0; + dReactionRate_dActivities[ r ][ i ] = 0.0; } } } @@ -212,7 +229,7 @@ KineticReactions< REAL_TYPE, { for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) { - reactionRatesDerivatives( r, i ) = forwardRateConstant * dProductConcForward_dC[i] - reverseRateConstant * dProductConcReverse_dC[i]; + dReactionRate_dActivities[ r ][ i ] = forwardRateConstant * dProductConcForward_dC[i] - reverseRateConstant * dProductConcReverse_dC[i]; } } } // end of if constexpr ( LOGE_CONCENTRATION ) @@ -238,31 +255,14 @@ KineticReactions< REAL_TYPE, LOGE_CONCENTRATION >::computeReactionRatesQuotient_impl( RealType const &, //temperature, PARAMS_DATA const & params, - typename ACTIVITY_MODEL::Params const & activityParams, - ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D_TO_CONST const & activities, ARRAY_1D_SA const & surfaceArea, ARRAY_1D & reactionRates, - ARRAY_2D & reactionRatesDerivatives ) + ARRAY_2D & dReactionRates_dActivities ) { if constexpr( !CALCULATE_DERIVATIVES ) { - 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 || s_ri < 0.0 ) { - reactionRatesDerivatives( r, i ) = -rateConstant * surfaceArea[r] * s_ri * quotient / ( equilibriumConstant * activities[i] ); + dReactionRates_dActivities[ r ][ i ] = -rateConstant * surfaceArea[r] * s_ri * quotient / ( equilibriumConstant * activities[i] ); } else { - reactionRatesDerivatives( r, i ) = 0.0; + dReactionRates_dActivities[ r ][ i ] = 0.0; } } } // end of if constexpr ( CALCULATE_DERIVATIVES ) @@ -355,25 +355,23 @@ 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_TO_CONST const & activities, ARRAY_1D & speciesRates, - ARRAY_2D & speciesRatesDerivatives ) + ARRAY_2D & dSpeciesRates_dActivities ) { RealType reactionRates[PARAMS_DATA::numReactions()] = { 0.0 }; - CArrayWrapper< double, PARAMS_DATA::numReactions(), PARAMS_DATA::numSpecies() > reactionRatesDerivatives; + CArrayWrapper< double, PARAMS_DATA::numReactions(), PARAMS_DATA::numSpecies() > dReactionRates_dActivities; if constexpr( !CALCULATE_DERIVATIVES ) { - HPCREACT_UNUSED_VAR( speciesRatesDerivatives ); + HPCREACT_UNUSED_VAR( dSpeciesRates_dActivities ); } - computeReactionRates< PARAMS_DATA >( temperature, + computeReactionRates< PARAMS_DATA, true >( temperature, params, - activityParams, - speciesConcentration, + activities, reactionRates, - reactionRatesDerivatives ); + dReactionRates_dActivities ); for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) { @@ -382,7 +380,7 @@ KineticReactions< REAL_TYPE, { for( IntType j = 0; j < PARAMS_DATA::numSpecies(); ++j ) { - speciesRatesDerivatives( i, j ) = 0.0; + dSpeciesRates_dActivities[ i ][ j ] = 0.0; } } for( IntType r=0; r 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, @@ -206,7 +208,6 @@ class MixedEquilibriumKineticReactions ARRAY_1D, ARRAY_2D, true >( params, - activityParams, speciesConcentration, reactionRates, reactionRatesDerivatives, @@ -254,7 +255,6 @@ 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, @@ -286,7 +286,7 @@ class MixedEquilibriumKineticReactions static HPCREACT_HOST_DEVICE void computeReactionRates_impl( RealType const & temperature, PARAMS_DATA const & params, - typename ACTIVITY_MODEL::Params const & activityParams, + 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, @@ -320,7 +320,6 @@ 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, diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp index 92aabd9..2a23ff8 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp @@ -132,6 +132,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, LOGE_CONCENTRATION >::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, @@ -166,6 +167,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, kineticReactions::computeReactionRates( temperature, params.kineticReactionsParameters(), + activityParams, logSpeciesConcentration, surfaceArea, reactionRates, @@ -182,8 +184,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, { RealType const dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations = params.stoichiometricMatrix( k, j + numSecondarySpecies ); - dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) += - reactionRatesDerivatives( i, k ) * dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations; + dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) += reactionRatesDerivatives( i, k ) * dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations; } } } diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp index 0bbd7d4..06c3b61 100644 --- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp @@ -104,13 +104,13 @@ void timeStepTest( PARAMS_DATA const & params, aggregatePrimarySpeciesConcentration_n[i] = aggregatePrimarySpeciesConcentration[i]; } - auto computeResidualAndJacobian = [&] ( REAL_TYPE const (&X)[numPrimarySpecies], + auto computeResidualAndJacobian = [&] ( REAL_TYPE const (&logPrimarySpeciesConcentration)[numPrimarySpecies], REAL_TYPE ( &r )[numPrimarySpecies], REAL_TYPE ( &J )[numPrimarySpecies][numPrimarySpecies] ) { MixedReactionsType::updateMixedSystem( temperature, params, - X, + logPrimarySpeciesConcentration, surfaceArea, logSecondarySpeciesConcentration, aggregatePrimarySpeciesConcentration, From 0e4e5ceaa96c6c62bd1bbeac2cda1f59b9a621a6 Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Wed, 25 Feb 2026 11:08:01 -0800 Subject: [PATCH 7/8] more changes --- src/CMakeLists.txt | 2 +- .../unitTests/testEquilibriumReactions.cpp | 2 + .../unitTests/testMomasEasyCase.cpp | 5 +- .../unitTests/testMomasMediumCase.cpp | 5 +- src/reactions/geochemistry/Carbonate.hpp | 36 ++ .../testGeochemicalEquilibriumReactions.cpp | 3 + .../testGeochemicalMixedReactions.cpp | 1 + src/reactions/massActions/MassActions.hpp | 486 +++++++++++++++--- .../reactionsSystems/EquilibriumReactions.hpp | 6 + ...ionsAggregatePrimaryConcentration_impl.hpp | 96 +++- ...uilibriumReactionsReactionExtents_impl.hpp | 47 +- .../reactionsSystems/KineticReactions.hpp | 45 +- .../MixedEquilibriumKineticReactions.hpp | 11 +- .../MixedEquilibriumKineticReactions_impl.hpp | 128 ++++- .../equilibriumReactionsTestUtilities.hpp | 8 +- .../mixedReactionsTestUtilities.hpp | 3 + 16 files changed, 745 insertions(+), 139 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0b7dbf8..d3eefb0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,7 @@ set( hpcReact_headers common/macros.hpp common/CArrayWrapper.hpp + constitutive/activity/activity.hpp constitutive/activity/Bdot.hpp reactions/exampleSystems/BulkGeneric.hpp reactions/geochemistry/Carbonate.hpp @@ -78,4 +79,3 @@ add_subdirectory( docs ) hpcReact_add_code_checks( PREFIX hpcReact UNCRUSTIFY_CFG_FILE ${PROJECT_SOURCE_DIR}/src/uncrustify.cfg EXCLUDES "blt/*" ) - diff --git a/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp b/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp index d263fc1..9c9eea5 100644 --- a/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp +++ b/src/reactions/exampleSystems/unitTests/testEquilibriumReactions.cpp @@ -33,6 +33,7 @@ TEST( testEquilibriumReactions, computeResidualAndJacobianTest ) { -2.0, 4.0e16 } }; computeResidualAndJacobianTest< double, 2 >( bulkGeneric::simpleTestRateParams, + bulkGeneric::simpleActivityTestParams, initialSpeciesConcentration, expectedResiduals, expectedJacobian ); @@ -49,6 +50,7 @@ TEST( testEquilibriumReactions, testEnforceEquilibrium ) std::cout<<" RESIDUAL_FORM 2:"<( bulkGeneric::simpleTestRateParams.equilibriumReactionsParameters(), + bulkGeneric::simpleActivityTestParams, initialSpeciesConcentration, expectedSpeciesConcentrations ); diff --git a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp index 6745a96..d5ea6a1 100644 --- a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp +++ b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp @@ -12,6 +12,7 @@ #include "reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp" #include "../MoMasBenchmark.hpp" +#include "constitutive/activity/Identity.hpp" using namespace hpcReact; using namespace hpcReact::MoMasBenchmark; @@ -30,12 +31,13 @@ void testMoMasAllEquilibriumHelper() using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, int, int, - Bdot< double, int, SpeciatedIonicStrength< double, int, numSpecies > >>; + Identity< double, int, SpeciatedIonicStrength< double, int, numSpecies > >>; double logPrimarySpeciesConcentration[numPrimarySpecies]; pmpl::genericKernelWrapper( numPrimarySpecies, logPrimarySpeciesConcentration, [] HPCREACT_DEVICE ( auto * const logPrimarySpeciesConcentrationCopy ) { + Identity< double, int, SpeciatedIonicStrength< double, int, numSpecies > >::Params const activityParams = {}; double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 1.0e-20, // X1 @@ -65,6 +67,7 @@ void testMoMasAllEquilibriumHelper() EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0, hpcReact::MoMasBenchmark::easyCaseParams.equilibriumReactionsParameters(), + activityParams, targetAggregatePrimarySpeciesConcentration, logInitialPrimarySpeciesConcentration, logPrimarySpeciesConcentrationCopy ); diff --git a/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp index ec86b19..3751eb1 100644 --- a/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp +++ b/src/reactions/exampleSystems/unitTests/testMomasMediumCase.cpp @@ -12,6 +12,7 @@ #include "reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp" #include "../MoMasBenchmark.hpp" +#include "constitutive/activity/Identity.hpp" using namespace hpcReact; using namespace hpcReact::MoMasBenchmark; @@ -29,7 +30,7 @@ void testMoMasMediumEquilibriumHelper() using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, int, int, - Bdot< double, int, SpeciatedIonicStrength< double, int, numSpecies > >>; + Identity< double, int, SpeciatedIonicStrength< double, int, numSpecies > >>; @@ -37,6 +38,7 @@ void testMoMasMediumEquilibriumHelper() pmpl::genericKernelWrapper( numPrimarySpecies, logPrimarySpeciesConcentration, [] HPCREACT_DEVICE ( auto * const logPrimarySpeciesConcentrationCopy ) { + Identity< double, int, SpeciatedIonicStrength< double, int, numSpecies > >::Params const activityParams = {}; double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 1.0e-20, // X1 @@ -66,6 +68,7 @@ void testMoMasMediumEquilibriumHelper() EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0, hpcReact::MoMasBenchmark::mediumCaseParams.equilibriumReactionsParameters(), + activityParams, targetAggregatePrimarySpeciesConcentration, logInitialPrimarySpeciesConcentration, logPrimarySpeciesConcentrationCopy ); diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp index 9b4b9f6..2089d4a 100644 --- a/src/reactions/geochemistry/Carbonate.hpp +++ b/src/reactions/geochemistry/Carbonate.hpp @@ -186,11 +186,39 @@ using carbonateSystemType = reactionsSystems::MixedReactionsParame using carbonateIonicStrengthType = SpeciatedIonicStrength< double, int, 17 >; using carbonateActivityType = Bdot< double, int, carbonateIonicStrengthType >; using carbonateIdentityActivityType = Identity< double, int, carbonateIonicStrengthType >; +using carbonateNosolidIonicStrengthType = SpeciatedIonicStrength< double, int, 16 >; +using carbonateNosolidActivityType = Bdot< double, int, carbonateNosolidIonicStrengthType >; +using carbonateNosolidIdentityActivityType = Identity< double, int, carbonateNosolidIonicStrengthType >; 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 CArrayWrapper< double, 16 > carbonateNosolidSpeciesCharge = +{ + carbonate::speciesCharge[0], carbonate::speciesCharge[1], carbonate::speciesCharge[2], carbonate::speciesCharge[3], + carbonate::speciesCharge[4], carbonate::speciesCharge[5], carbonate::speciesCharge[6], carbonate::speciesCharge[7], + carbonate::speciesCharge[8], carbonate::speciesCharge[10], carbonate::speciesCharge[11], carbonate::speciesCharge[12], + carbonate::speciesCharge[13], carbonate::speciesCharge[14], carbonate::speciesCharge[15], carbonate::speciesCharge[16] +}; + +constexpr CArrayWrapper< double, 16 > carbonateNosolidIonSize = +{ + carbonate::ionSize[0], carbonate::ionSize[1], carbonate::ionSize[2], carbonate::ionSize[3], + carbonate::ionSize[4], carbonate::ionSize[5], carbonate::ionSize[6], carbonate::ionSize[7], + carbonate::ionSize[8], carbonate::ionSize[10], carbonate::ionSize[11], carbonate::ionSize[12], + carbonate::ionSize[13], carbonate::ionSize[14], carbonate::ionSize[15], carbonate::ionSize[16] +}; + +constexpr CArrayWrapper< double, 16 > carbonateNosolidBdotParameters = +{ + carbonate::bdotParameters[0], carbonate::bdotParameters[1], carbonate::bdotParameters[2], carbonate::bdotParameters[3], + carbonate::bdotParameters[4], carbonate::bdotParameters[5], carbonate::bdotParameters[6], carbonate::bdotParameters[7], + carbonate::bdotParameters[8], carbonate::bdotParameters[10], carbonate::bdotParameters[11], carbonate::bdotParameters[12], + carbonate::bdotParameters[13], carbonate::bdotParameters[14], carbonate::bdotParameters[15], carbonate::bdotParameters[16] +}; + constexpr carbonateActivityType::Params carbonateActivityParams = { {carbonate::speciesCharge}, @@ -198,9 +226,17 @@ constexpr carbonateActivityType::Params carbonateActivityParams = carbonate::bdotParameters }; +constexpr carbonateNosolidActivityType::Params carbonateNosolidActivityParams = +{ + {carbonateNosolidSpeciesCharge}, + carbonateNosolidIonSize, + carbonateNosolidBdotParameters +}; + Identity< double, int, carbonateIonicStrengthType >::Params carbonateIdentityActivityParams = {}; +Identity< double, int, carbonateNosolidIonicStrengthType >::Params carbonateNosolidIdentityActivityParams = {}; // *****UNCRUSTIFY-ON****** } // namespace geochemistry diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp index 9f9fd3e..dc13788 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp @@ -82,6 +82,7 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium ) std::cout<<" RESIDUAL_FORM 0:"<( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), + hpcReact::geochemistry::carbonateActivityParams, initialSpeciesConcentration, expectedSpeciesConcentrations ); @@ -92,6 +93,7 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium ) std::cout<<" RESIDUAL_FORM 2:"<( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), + hpcReact::geochemistry::carbonateActivityParams, initialSpeciesConcentration, expectedSpeciesConcentrations ); @@ -137,6 +139,7 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 ) double logPrimarySpeciesConcentration[numPrimarySpecies]; EquilibriumReactionsType::enforceEquilibrium_LogAggregate( 0, hpcReact::geochemistry::carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), + hpcReact::geochemistry::carbonateActivityParams, logInitialPrimarySpeciesConcentration, logPrimarySpeciesConcentration ); diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp index 59f90bf..87e2000 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalMixedReactions.cpp @@ -54,6 +54,7 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem ) using ActivityModelType = Bdot< double, int, SpeciatedIonicStrength< double, int, carbonateSystemType::numSpecies() > >; timeStepTest< double, true, ActivityModelType >( carbonateSystem, + carbonateNosolidActivityParams, 1.0, 10, initialAggregateSpeciesConcentration, diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index c77d86b..03d8259 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -34,9 +34,9 @@ template< typename REAL_TYPE, HPCREACT_HOST_DEVICE inline void calculateLogSecondaryActivities( PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & logPrimaryActivities, - ARRAY_1D & logSecondaryActivities, - FUNC && derivativeFunc ) + ARRAY_1D_TO_CONST const & logPrimaryActivities, + ARRAY_1D & logSecondaryActivities, + FUNC && derivativeFunc ) { static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); @@ -57,6 +57,196 @@ void calculateLogSecondaryActivities( PARAMS_DATA const & params, } } +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_2D_TO_CONST, + typename ARRAY_1D_SECONDARY, + typename ARRAY_1D_PRIMARY, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE +inline +void calculateAggregatePrimaryConcentrationsFromActivitiesWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST2 const & logPrimaryActivities, + ARRAY_2D_TO_CONST const & dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondaryActivities, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) +{ + static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + + for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) + { + for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) + { + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; + } + } + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + REAL_TYPE const primarySpeciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); + aggregatePrimarySpeciesConcentrations[i] = primarySpeciesConcentration_i; + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = primarySpeciesConcentration_i; + } + + if constexpr( numSecondarySpecies <= 0 ) + { + return; + } + + massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimaryActivities, + logSecondaryActivities, + []( INDEX_TYPE, INDEX_TYPE, REAL_TYPE ){} ); + + REAL_TYPE dLogSecondaryActivities_dLogPrimaryActivities[numSecondarySpecies][numPrimarySpecies] = {{ 0.0 }}; + massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimaryActivities, + logSecondaryActivities, + [&]( const int j, const int k, REAL_TYPE const value ) + { + dLogSecondaryActivities_dLogPrimaryActivities[j][k] = value; + } ); + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + for( int j = 0; j < numSecondarySpecies; ++j ) + { + REAL_TYPE const secondarySpeciesConcentration_j = exp( logSecondaryActivities[j] ); + REAL_TYPE const nu_ji = params.stoichiometricMatrix( j, i + numSecondarySpecies ); + + aggregatePrimarySpeciesConcentrations[i] += nu_ji * secondarySpeciesConcentration_j; + + for( int k = 0; k < numPrimarySpecies; ++k ) + { + REAL_TYPE dLogSecondaryActivity_j_dLogPrimarySpeciesConcentration_k = 0.0; + for( int l = 0; l < numPrimarySpecies; ++l ) + { + dLogSecondaryActivity_j_dLogPrimarySpeciesConcentration_k += + dLogSecondaryActivities_dLogPrimaryActivities[j][l] * + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[l][k]; + } + + REAL_TYPE const dSecondarySpeciesConcentration_j_dLogPrimarySpeciesConcentration_k = + secondarySpeciesConcentration_j * dLogSecondaryActivity_j_dLogPrimarySpeciesConcentration_k; + + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += + nu_ji * dSecondarySpeciesConcentration_j_dLogPrimarySpeciesConcentration_k; + } + } + } +} + +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_2D_TO_CONST, + typename ARRAY_1D_SECONDARY, + typename ARRAY_1D_PRIMARY, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE +inline +void calculateTotalAndMobileAggregatePrimaryConcentrationsFromActivitiesWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST2 const & logPrimaryActivities, + ARRAY_2D_TO_CONST const & dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondaryActivities, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_1D_PRIMARY & mobileAggregatePrimarySpeciesConcentrations, + ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations, + ARRAY_2D & dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) +{ + static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + + for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) + { + for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) + { + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; + dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; + } + } + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + REAL_TYPE const primarySpeciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); + aggregatePrimarySpeciesConcentrations[i] = primarySpeciesConcentration_i; + mobileAggregatePrimarySpeciesConcentrations[i] = primarySpeciesConcentration_i; + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = primarySpeciesConcentration_i; + dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = primarySpeciesConcentration_i; + } + + if constexpr( numSecondarySpecies <= 0 ) + { + return; + } + + massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimaryActivities, + logSecondaryActivities, + []( INDEX_TYPE, INDEX_TYPE, REAL_TYPE ){} ); + + REAL_TYPE dLogSecondaryActivities_dLogPrimaryActivities[numSecondarySpecies][numPrimarySpecies] = {{ 0.0 }}; + massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimaryActivities, + logSecondaryActivities, + [&]( const int j, const int k, REAL_TYPE const value ) + { + dLogSecondaryActivities_dLogPrimaryActivities[j][k] = value; + } ); + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + for( int j = 0; j < numSecondarySpecies; ++j ) + { + REAL_TYPE const secondarySpeciesConcentration_j = exp( logSecondaryActivities[j] ); + REAL_TYPE const nu_ji = params.stoichiometricMatrix( j, i + numSecondarySpecies ); + REAL_TYPE const mobileFlag = params.mobileSecondarySpeciesFlag( j ); + + aggregatePrimarySpeciesConcentrations[i] += nu_ji * secondarySpeciesConcentration_j; + mobileAggregatePrimarySpeciesConcentrations[i] += nu_ji * secondarySpeciesConcentration_j * mobileFlag; + + for( int k = 0; k < numPrimarySpecies; ++k ) + { + REAL_TYPE dLogSecondaryActivity_j_dLogPrimarySpeciesConcentration_k = 0.0; + for( int l = 0; l < numPrimarySpecies; ++l ) + { + dLogSecondaryActivity_j_dLogPrimarySpeciesConcentration_k += + dLogSecondaryActivities_dLogPrimaryActivities[j][l] * + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[l][k]; + } + + REAL_TYPE const dSecondarySpeciesConcentration_j_dLogPrimarySpeciesConcentration_k = + secondarySpeciesConcentration_j * dLogSecondaryActivity_j_dLogPrimarySpeciesConcentration_k; + + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += + nu_ji * dSecondarySpeciesConcentration_j_dLogPrimarySpeciesConcentration_k; + + dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += + nu_ji * dSecondarySpeciesConcentration_j_dLogPrimarySpeciesConcentration_k * mobileFlag; + } + } + } +} + } // namespace template< typename REAL_TYPE, @@ -68,8 +258,8 @@ template< typename REAL_TYPE, HPCREACT_HOST_DEVICE inline void calculateLogSecondaryActivities( PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & logPrimaryActivities, - ARRAY_1D & logSecondaryActivities ) + ARRAY_1D_TO_CONST const & logPrimaryActivities, + ARRAY_1D & logSecondaryActivities ) { if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) { @@ -77,11 +267,11 @@ void calculateLogSecondaryActivities( PARAMS_DATA const & params, } massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimaryActivities, - logSecondaryActivities, - [](INDEX_TYPE, INDEX_TYPE, REAL_TYPE ){} ); + INT_TYPE, + INDEX_TYPE >( params, + logPrimaryActivities, + logSecondaryActivities, + []( INDEX_TYPE, INDEX_TYPE, REAL_TYPE ){} ); } @@ -100,14 +290,55 @@ void calculateLogSecondaryActivitiesWrtLogC( PARAMS_DATA const & params, ARRAY_2D & dLogSecondaryActivities_dLogPrimaryActivities ) { massActions_impl::calculateLogSecondaryActivities< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, - logPrimaryActivities, - logSecondaryActivities, - [&]( const int j, const int k, REAL_TYPE const value ) + logPrimaryActivities, + logSecondaryActivities, + [&]( const int j, const int k, REAL_TYPE const value ) { dLogSecondaryActivities_dLogPrimaryActivities[j][k] = value; } ); } +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D > +HPCREACT_HOST_DEVICE +inline +void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D & logSecondarySpeciesConcentrations ) +{ + calculateLogSecondaryActivities< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations ); +} + +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE +inline +void calculateLogSecondarySpeciesConcentrationWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D & logSecondarySpeciesConcentrations, + ARRAY_2D & dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ) +{ + calculateLogSecondaryActivitiesWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); +} + template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, @@ -125,38 +356,22 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) { static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); - static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); - calculateLogSecondaryActivities< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations ); + REAL_TYPE dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[numPrimarySpecies][numPrimarySpecies] = {{ 0.0 }}; for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) { - for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) - { - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; - } + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[i][i] = 1.0; } - for( int i = 0; i < numPrimarySpecies; ++i ) - { - REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); - aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; - for( int j = 0; j < numSecondarySpecies; ++j ) - { - REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] ); - aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j; - for( int k=0; k( params, + logPrimarySpeciesConcentrations, + logPrimarySpeciesConcentrations, + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); } template< typename REAL_TYPE, @@ -170,7 +385,6 @@ HPCREACT_HOST_DEVICE inline void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, - ARRAY_1D_TO_CONST const & logPrimaryActivities, ARRAY_1D & aggregatePrimarySpeciesConcentrations, ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) { @@ -178,8 +392,7 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, if constexpr( numSecondarySpecies > 0 ) { - REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0}; - + REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = { 0.0 }; calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, @@ -190,9 +403,100 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, } else { - GEOS_UNUSED_VAR( logPrimarySpeciesConcentrations, aggregatePrimarySpeciesConcentrations, dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( logPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( aggregatePrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); } +} +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE +inline +void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST const & logPrimaryActivities, + ARRAY_1D & aggregatePrimarySpeciesConcentrations, + ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) +{ + static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + + if constexpr( numSecondarySpecies > 0 ) + { + REAL_TYPE dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[numPrimarySpecies][numPrimarySpecies] = {{ 0.0 }}; + for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) + { + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[i][i] = 1.0; + } + + REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0}; + + massActions_impl::calculateAggregatePrimaryConcentrationsFromActivitiesWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logPrimaryActivities, + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } + else + { + HPCREACT_UNUSED_VAR( logPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( logPrimaryActivities ); + HPCREACT_UNUSED_VAR( aggregatePrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } +} + +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_2D_TO_CONST, + typename ARRAY_1D, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE +inline +void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST2 const & logPrimaryActivities, + ARRAY_2D_TO_CONST const & dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + ARRAY_1D & aggregatePrimarySpeciesConcentrations, + ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) +{ + static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + + if constexpr( numSecondarySpecies > 0 ) + { + REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = { 0.0 }; + massActions_impl::calculateAggregatePrimaryConcentrationsFromActivitiesWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logPrimaryActivities, + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } + else + { + HPCREACT_UNUSED_VAR( logPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( logPrimaryActivities ); + HPCREACT_UNUSED_VAR( dLogPrimaryActivities_dLogPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( aggregatePrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } } template< typename REAL_TYPE, @@ -214,46 +518,74 @@ void calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA c ARRAY_2D & dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) { static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); - static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); - calculateLogSecondarySpeciesConcentration< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations ); + REAL_TYPE dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[numPrimarySpecies][numPrimarySpecies] = {{ 0.0 }}; for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) { - for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) - { - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; - dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; - } + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[i][i] = 1.0; } - for( int i = 0; i < numPrimarySpecies; ++i ) - { - REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); - aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; - mobileAggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; - dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; + massActions_impl::calculateTotalAndMobileAggregatePrimaryConcentrationsFromActivitiesWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logPrimarySpeciesConcentrations, + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + mobileAggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations, + dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); +} - for( int j = 0; j < numSecondarySpecies; ++j ) - { - REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] ); - aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j; - mobileAggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j * params.mobileSecondarySpeciesFlag( j ); - for( int k=0; k +HPCREACT_HOST_DEVICE +inline +void calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST2 const & logPrimaryActivities, + ARRAY_2D_TO_CONST const & dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_1D_PRIMARY & mobileAggregatePrimarySpeciesConcentrations, + ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations, + ARRAY_2D & dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) +{ + static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); - dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * - dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration * - params.mobileSecondarySpeciesFlag( j ); - } - } + if constexpr( numSecondarySpecies > 0 ) + { + massActions_impl::calculateTotalAndMobileAggregatePrimaryConcentrationsFromActivitiesWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logPrimaryActivities, + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + mobileAggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations, + dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } + else + { + HPCREACT_UNUSED_VAR( logPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( logPrimaryActivities ); + HPCREACT_UNUSED_VAR( dLogPrimaryActivities_dLogPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( logSecondarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( aggregatePrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( mobileAggregatePrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); } } diff --git a/src/reactions/reactionsSystems/EquilibriumReactions.hpp b/src/reactions/reactionsSystems/EquilibriumReactions.hpp index 8e0b79f..8e0baea 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactions.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactions.hpp @@ -15,6 +15,7 @@ #include "common/CArrayWrapper.hpp" #include "common/DirectSystemSolve.hpp" #include "common/printers.hpp" +#include "constitutive/activity/activity.hpp" #include @@ -71,6 +72,7 @@ class EquilibriumReactions void enforceEquilibrium_Extents( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration0, ARRAY_1D & speciesConcentration ); @@ -95,6 +97,7 @@ class EquilibriumReactions void enforceEquilibrium_LogAggregate( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration0, ARRAY_1D & speciesConcentration ); @@ -125,6 +128,7 @@ class EquilibriumReactions void enforceEquilibrium_Aggregate( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & targetAggregatePrimarySpeciesConcentration, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, ARRAY_1D & speciesConcentration ); @@ -152,6 +156,7 @@ class EquilibriumReactions static HPCREACT_HOST_DEVICE void computeResidualAndJacobianReactionExtents( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & speciesConcentration0, ARRAY_1D_TO_CONST2 const & xi, ARRAY_1D & residual, @@ -180,6 +185,7 @@ class EquilibriumReactions static HPCREACT_HOST_DEVICE void computeResidualAndJacobianAggregatePrimaryConcentrations( RealType const & temperature, PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & targetAggregatePrimaryConcentrations, ARRAY_1D_TO_CONST2 const & logPrimarySpeciesConcentration, ARRAY_1D & residual, diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index c3222f6..380bb08 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -37,24 +37,84 @@ EquilibriumReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, ACTIVITY_MODEL >::computeResidualAndJacobianAggregatePrimaryConcentrations( RealType const & temperature, - PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & targetAggregatePrimaryConcentrations, - ARRAY_1D_TO_CONST2 const & logPrimarySpeciesConcentration, - ARRAY_1D & residual, - ARRAY_2D & jacobian ) + PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, + ARRAY_1D_TO_CONST const & targetAggregatePrimaryConcentrations, + ARRAY_1D_TO_CONST2 const & logPrimarySpeciesConcentration, + ARRAY_1D & residual, + ARRAY_2D & jacobian ) { HPCREACT_UNUSED_VAR( temperature ); + static constexpr int numSpecies = PARAMS_DATA::numSpecies(); + static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + static constexpr int numSecondarySpeciesStorage = numSecondarySpecies > 0 ? numSecondarySpecies : 1; static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); RealType aggregatePrimaryConcentrations[numPrimarySpecies] = {0.0}; + RealType logPrimaryActivities[numPrimarySpecies] = {0.0}; + RealType dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[numPrimarySpecies][numPrimarySpecies] = {{0.0}}; ARRAY_2D dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations = {{{0.0}}}; - + if constexpr( numPrimarySpecies > 0 ) + { + RealType logSpeciesConcentration[numSpecies] = {0.0}; + RealType logActivities[numSpecies] = {0.0}; + RealType dLogActivities_dLogSpeciesConcentrations[numSpecies][numSpecies] = {{0.0}}; + + RealType dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[numSecondarySpeciesStorage][numPrimarySpecies] = {{0.0}}; + if constexpr( numSecondarySpecies > 0 ) + { + RealType logSecondarySpeciesConcentrations[numSecondarySpecies] = {0.0}; + massActions::calculateLogSecondarySpeciesConcentrationWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentration, + logSecondarySpeciesConcentrations, + dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + for( IndexType i = 0; i < numSecondarySpecies; ++i ) + { + logSpeciesConcentration[i] = logSecondarySpeciesConcentrations[i]; + } + } + + for( IndexType i = 0; i < numPrimarySpecies; ++i ) + { + logSpeciesConcentration[i + numSecondarySpecies] = logPrimarySpeciesConcentration[i]; + } + + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + true >( activityParams, + logSpeciesConcentration, + logActivities, + dLogActivities_dLogSpeciesConcentrations ); + + for( IndexType i = 0; i < numPrimarySpecies; ++i ) + { + IndexType const fullRow = i + numSecondarySpecies; + logPrimaryActivities[i] = logActivities[fullRow]; + + for( IndexType j = 0; j < numPrimarySpecies; ++j ) + { + RealType value = dLogActivities_dLogSpeciesConcentrations[fullRow][j + numSecondarySpecies]; + for( IndexType k = 0; k < numSecondarySpecies; ++k ) + { + value += dLogActivities_dLogSpeciesConcentrations[fullRow][k] * + dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[k][j]; + } + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[i][j] = value; + } + } + } massActions::calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, - logPrimarySpeciesConcentration, - aggregatePrimaryConcentrations, - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + logPrimarySpeciesConcentration, + logPrimaryActivities, + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, + aggregatePrimaryConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); for( IndexType i=0; i::enforceEquilibrium_LogAggregate( REAL_TYPE const & temperature, - PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, - ARRAY_1D & logPrimarySpeciesConcentration ) + PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, + ARRAY_1D & logPrimarySpeciesConcentration ) { HPCREACT_UNUSED_VAR( temperature ); static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); @@ -97,6 +158,7 @@ EquilibriumReactions< REAL_TYPE, enforceEquilibrium_Aggregate( temperature, params, + activityParams, targetAggregatePrimarySpeciesConcentration, logPrimarySpeciesConcentration0, logPrimarySpeciesConcentration ); @@ -117,10 +179,11 @@ EquilibriumReactions< REAL_TYPE, INT_TYPE, 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, - ARRAY_1D & logPrimarySpeciesConcentration ) + PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, + ARRAY_1D_TO_CONST const & targetAggregatePrimarySpeciesConcentration, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, + ARRAY_1D & logPrimarySpeciesConcentration ) { if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) { @@ -150,6 +213,7 @@ EquilibriumReactions< REAL_TYPE, { computeResidualAndJacobianAggregatePrimaryConcentrations( temperature, params, + activityParams, targetAggregatePrimarySpeciesConcentration, logPrimarySpeciesConcentration, residual, diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp index d3742d1..c14df63 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsReactionExtents_impl.hpp @@ -37,11 +37,12 @@ EquilibriumReactions< REAL_TYPE, INT_TYPE, 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, - ARRAY_1D & residual, - ARRAY_2D & jacobian ) + PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, + ARRAY_1D_TO_CONST const & speciesConcentration0, + ARRAY_1D_TO_CONST2 const & xi, + ARRAY_1D & residual, + ARRAY_2D & jacobian ) { HPCREACT_UNUSED_VAR( temperature ); @@ -59,6 +60,16 @@ EquilibriumReactions< REAL_TYPE, } } + RealType activities[numSpecies] = { 0.0 }; + RealType dActivities_dConcentration[numSpecies][numSpecies] = {{ 0.0 }}; + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + false >( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); // loop over reactions for( IndexType a=0; a 0.0 ) { // reverse reaction - reverseProduct *= pow( speciesConcentration[i], s_ai ); + reverseProduct *= pow( activities[i], s_ai ); // derivative of reverse product with respect to xi for( IndexType b=0; b::enforceEquilibrium_Extents( REAL_TYPE const & temperature, - PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & speciesConcentration0, - ARRAY_1D & speciesConcentration ) + PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, + ARRAY_1D_TO_CONST const & speciesConcentration0, + ARRAY_1D & speciesConcentration ) { HPCREACT_UNUSED_VAR( temperature ); static constexpr int numSpecies = PARAMS_DATA::numSpecies(); @@ -142,6 +164,7 @@ EquilibriumReactions< REAL_TYPE, { computeResidualAndJacobianReactionExtents( temperature, params, + activityParams, speciesConcentration0, xi, residual, diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp index 1a1b6ea..c609ef1 100644 --- a/src/reactions/reactionsSystems/KineticReactions.hpp +++ b/src/reactions/reactionsSystems/KineticReactions.hpp @@ -74,10 +74,14 @@ class KineticReactions RealType dActivities_dConcentration[ PARAMS_DATA::numSpecies() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; RealType dReactionRates_dActivities[ PARAMS_DATA::numReactions() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; - calculateActivities( activityParams, - speciesConcentration, - activities, - dActivities_dConcentration ); + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + LOGE_CONCENTRATION >( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); computeReactionRates< PARAMS_DATA, true >( temperature, params, @@ -125,17 +129,20 @@ class KineticReactions RealType dActivities_dConcentration[ PARAMS_DATA::numSpecies() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; RealType dReactionRates_dActivities[ PARAMS_DATA::numReactions() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; - calculateActivities( activityParams, - speciesConcentration, - activities, - dActivities_dConcentration ); + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + LOGE_CONCENTRATION >( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); computeReactionRates< PARAMS_DATA, false >( temperature, - params, - activityParams, - activities, - reactionRates, - dReactionRates_dActivities ); + params, + activities, + reactionRates, + dReactionRates_dActivities ); HPCREACT_UNUSED_VAR( dReactionRates_dActivities ); } @@ -291,10 +298,14 @@ class KineticReactions RealType activities[ PARAMS_DATA::numSpecies() ]; RealType dActivities_dConcentration[ PARAMS_DATA::numSpecies() ][ PARAMS_DATA::numSpecies() ] = { 0.0 }; - calculateActivities( activityParams, - speciesConcentration, - activities, - dActivities_dConcentration ); + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + LOGE_CONCENTRATION >( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); computeSpeciesRates_impl< PARAMS_DATA, false >( temperature, params, diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp index 8449757..abe0dbf 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, + typename ACTIVITY_MODEL::Params const & activityParams, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, ARRAY_1D_TO_CONST_KINETIC const & surfaceArea, ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, @@ -106,14 +107,9 @@ class MixedEquilibriumKineticReactions ARRAY_1D_PRIMARY & aggregateSpeciesRates, ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ) { - - RealType logPrimaryActivities[ PARAMS_DATA::numPrimarySpecies() ] = { 0.0 }; - RealType logSecondaryActivities[ PARAMS_DATA::numSecondarySpecies() ] = { 0.0 }; - RealType dLogPrimaryActivites_dLogPrimarySpeciesConcentrations[ PARAMS_DATA::numPrimarySpecies() ][ PARAMS_DATA::numPrimarySpecies() ] = { { 0.0 } }; - RealType dLogSecondaryActivites_dLogPrimarySpeciesConcentrations[ PARAMS_DATA::numSecondarySpecies() ][ PARAMS_DATA::numPrimarySpecies() ] = { { 0.0 } }; - updateMixedSystem_impl( temperature, params, + activityParams, logPrimarySpeciesConcentrations, surfaceArea, logSecondarySpeciesConcentrations, @@ -153,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, @@ -163,6 +160,7 @@ class MixedEquilibriumKineticReactions computeReactionRates_impl( temperature, params, + activityParams, logPrimarySpeciesConcentrations, logSecondarySpeciesConcentrations, surfaceArea, @@ -255,6 +253,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, diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp index 2a23ff8..06d2566 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp @@ -48,6 +48,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, LOGE_CONCENTRATION >::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, @@ -60,13 +61,68 @@ MixedEquilibriumKineticReactions< REAL_TYPE, ARRAY_1D_PRIMARY & aggregateSpeciesRates, ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ) { + constexpr IntType numSpecies = PARAMS_DATA::numSpecies(); + constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + if constexpr( PARAMS_DATA::numEquilibriumReactions() > 0 ) { + RealType logSpeciesConcentration[numSpecies] = { 0.0 }; + RealType logSpeciesActivities[numSpecies] = { 0.0 }; + RealType dLogSpeciesActivities_dLogSpeciesConcentrations[numSpecies][numSpecies] = {{ 0.0 }}; + RealType logPrimaryActivities[numPrimarySpecies] = { 0.0 }; + RealType dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[numPrimarySpecies][numPrimarySpecies] = {{ 0.0 }}; + RealType logSecondarySpeciesConcentrations_guess[numSecondarySpecies] = { 0.0 }; + RealType dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[numSecondarySpecies][numPrimarySpecies] = {{ 0.0 }}; + + massActions::calculateLogSecondarySpeciesConcentrationWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params.equilibriumReactionsParameters(), + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations_guess, + dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + for( IntType i = 0; i < numSecondarySpecies; ++i ) + { + logSpeciesConcentration[i] = logSecondarySpeciesConcentrations_guess[i]; + } + for( IntType i = 0; i < numPrimarySpecies; ++i ) + { + logSpeciesConcentration[i + numSecondarySpecies] = logPrimarySpeciesConcentrations[i]; + } + + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + true >( activityParams, + logSpeciesConcentration, + logSpeciesActivities, + dLogSpeciesActivities_dLogSpeciesConcentrations ); + + for( IntType i = 0; i < numPrimarySpecies; ++i ) + { + IntType const fullRow = i + numSecondarySpecies; + logPrimaryActivities[i] = logSpeciesActivities[fullRow]; + + for( IntType j = 0; j < numPrimarySpecies; ++j ) + { + RealType value = dLogSpeciesActivities_dLogSpeciesConcentrations[fullRow][j + numSecondarySpecies]; + for( IntType k = 0; k < numSecondarySpecies; ++k ) + { + value += dLogSpeciesActivities_dLogSpeciesConcentrations[fullRow][k] * + dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[k][j]; + } + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[i][j] = value; + } + } + // 1. Compute new aggregate species from primary species massActions::calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params.equilibriumReactionsParameters(), logPrimarySpeciesConcentrations, + logPrimaryActivities, + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations, logSecondarySpeciesConcentrations, aggregatePrimarySpeciesConcentrations, mobileAggregatePrimarySpeciesConcentrations, @@ -75,8 +131,6 @@ MixedEquilibriumKineticReactions< REAL_TYPE, } else { - constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); - for( int i = 0; i < numPrimarySpecies; ++i ) { REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); @@ -92,6 +146,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE, // 2. Compute the reaction rates for all kinetic reactions computeReactionRates( temperature, params, + activityParams, logPrimarySpeciesConcentrations, logSecondarySpeciesConcentrations, surfaceArea, @@ -108,7 +163,10 @@ MixedEquilibriumKineticReactions< REAL_TYPE, } else { - GEOS_UNUSED_VAR( reactionRates, dReactionRates_dLogPrimarySpeciesConcentrations, aggregateSpeciesRates, dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( reactionRates ); + HPCREACT_UNUSED_VAR( dReactionRates_dLogPrimarySpeciesConcentrations ); + HPCREACT_UNUSED_VAR( aggregateSpeciesRates ); + HPCREACT_UNUSED_VAR( dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); } } @@ -179,12 +237,70 @@ MixedEquilibriumKineticReactions< REAL_TYPE, for( IntType j = 0; j < numPrimarySpecies; ++j ) { dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) = reactionRatesDerivatives( i, j + numSecondarySpecies ); + } + } + + if constexpr( numSecondarySpecies > 0 ) + { + RealType logSpeciesActivities[numSpecies] = { 0.0 }; + RealType dLogSpeciesActivities_dLogSpeciesConcentrations[numSpecies][numSpecies] = {{ 0.0 }}; + RealType logPrimaryActivities[numPrimarySpecies] = { 0.0 }; + RealType dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[numPrimarySpecies][numPrimarySpecies] = {{ 0.0 }}; + RealType logSecondaryActivities[numSecondarySpecies] = { 0.0 }; + RealType dLogSecondaryActivities_dLogPrimaryActivities[numSecondarySpecies][numPrimarySpecies] = {{ 0.0 }}; + RealType dLogSecondaryActivities_dLogPrimarySpeciesConcentrations[numSecondarySpecies][numPrimarySpecies] = {{ 0.0 }}; - for( IntType k = 0; k < numSecondarySpecies; ++k ) + calculateActivities< RealType, + IntType, + IndexType, + ACTIVITY_MODEL, + true >( activityParams, + logSpeciesConcentration, + logSpeciesActivities, + dLogSpeciesActivities_dLogSpeciesConcentrations ); + + for( IntType i = 0; i < numPrimarySpecies; ++i ) + { + IntType const fullRow = i + numSecondarySpecies; + logPrimaryActivities[i] = logSpeciesActivities[fullRow]; + for( IntType j = 0; j < numPrimarySpecies; ++j ) { - RealType const dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations = params.stoichiometricMatrix( k, j + numSecondarySpecies ); + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[i][j] = + dLogSpeciesActivities_dLogSpeciesConcentrations[fullRow][j + numSecondarySpecies]; + } + } + + massActions::calculateLogSecondaryActivitiesWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params.equilibriumReactionsParameters(), + logPrimaryActivities, + logSecondaryActivities, + dLogSecondaryActivities_dLogPrimaryActivities ); + HPCREACT_UNUSED_VAR( logSecondaryActivities ); - dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) += reactionRatesDerivatives( i, k ) * dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations; + for( IntType k = 0; k < numSecondarySpecies; ++k ) + { + for( IntType j = 0; j < numPrimarySpecies; ++j ) + { + dLogSecondaryActivities_dLogPrimarySpeciesConcentrations[k][j] = 0.0; + for( IntType l = 0; l < numPrimarySpecies; ++l ) + { + dLogSecondaryActivities_dLogPrimarySpeciesConcentrations[k][j] += + dLogSecondaryActivities_dLogPrimaryActivities[k][l] * + dLogPrimaryActivities_dLogPrimarySpeciesConcentrations[l][j]; + } + } + } + + for( IntType i = 0; i < numKineticReactions; ++i ) + { + for( IntType j = 0; j < numPrimarySpecies; ++j ) + { + for( IntType k = 0; k < numSecondarySpecies; ++k ) + { + dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) += + reactionRatesDerivatives( i, k ) * dLogSecondaryActivities_dLogPrimarySpeciesConcentrations[k][j]; + } } } } diff --git a/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp index d259c0e..df1341a 100644 --- a/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp @@ -54,6 +54,7 @@ template< typename REAL_TYPE, int RESIDUAL_FORM, typename PARAMS_DATA > void computeResidualAndJacobianTest( PARAMS_DATA const & params, + typename Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >::Params const & activityParams, REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], REAL_TYPE const (&expectedResidual)[PARAMS_DATA::numReactions()], REAL_TYPE const (&expectedJacobian)[PARAMS_DATA::numReactions()][PARAMS_DATA::numReactions()] ) @@ -74,12 +75,13 @@ void computeResidualAndJacobianTest( PARAMS_DATA const & params, data.speciesConcentration[i] = initialSpeciesConcentration[i]; } - pmpl::genericKernelWrapper( 1, &data, [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy ) + pmpl::genericKernelWrapper( 1, &data, [params, temperature, activityParams] HPCREACT_DEVICE ( auto * const dataCopy ) { double xi[numReactions] = { 0.0 }; EquilibriumReactionsType::computeResidualAndJacobianReactionExtents( temperature, params, + activityParams, dataCopy->speciesConcentration, xi, dataCopy->residual, @@ -126,6 +128,7 @@ template< typename REAL_TYPE, int RESIDUAL_FORM, typename PARAMS_DATA > void testEnforceEquilibrium( PARAMS_DATA const & params, + typename Bdot< REAL_TYPE, int, SpeciatedIonicStrength< REAL_TYPE, int, PARAMS_DATA::numSpecies() > >::Params const & activityParams, REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], REAL_TYPE const (&expectedSpeciesConcentrations)[PARAMS_DATA::numSpecies()] ) { @@ -144,10 +147,11 @@ void testEnforceEquilibrium( PARAMS_DATA const & params, data.speciesConcentration0[i] = initialSpeciesConcentration[i]; } - pmpl::genericKernelWrapper( 1, &data, [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy ) + pmpl::genericKernelWrapper( 1, &data, [params, temperature, activityParams] HPCREACT_DEVICE ( auto * const dataCopy ) { EquilibriumReactionsType::enforceEquilibrium_Extents( temperature, params, + activityParams, dataCopy->speciesConcentration0, dataCopy->speciesConcentration ); }); diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp index 06c3b61..c0b86fb 100644 --- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp @@ -32,6 +32,7 @@ template< typename REAL_TYPE, typename ACTIVITY_MODEL, typename PARAMS_DATA > void timeStepTest( PARAMS_DATA const & params, + typename ACTIVITY_MODEL::Params const & activityParams, REAL_TYPE const dt, int const numSteps, REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numPrimarySpecies()], @@ -90,6 +91,7 @@ void timeStepTest( PARAMS_DATA const & params, EquilibriumReactionsType::enforceEquilibrium_LogAggregate( temperature, params.equilibriumReactionsParameters(), + activityParams, logPrimarySpeciesConcentration, logPrimarySpeciesConcentration ); @@ -110,6 +112,7 @@ void timeStepTest( PARAMS_DATA const & params, { MixedReactionsType::updateMixedSystem( temperature, params, + activityParams, logPrimarySpeciesConcentration, surfaceArea, logSecondarySpeciesConcentration, From 3bd3028562577baab35812cf32d3759fb924ee0d Mon Sep 17 00:00:00 2001 From: Randolph Settgast Date: Wed, 25 Feb 2026 11:24:01 -0800 Subject: [PATCH 8/8] missing file --- src/constitutive/activity/activity.hpp | 93 ++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 src/constitutive/activity/activity.hpp diff --git a/src/constitutive/activity/activity.hpp b/src/constitutive/activity/activity.hpp new file mode 100644 index 0000000..68c7cb7 --- /dev/null +++ b/src/constitutive/activity/activity.hpp @@ -0,0 +1,93 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" + +#include + +namespace hpcReact +{ + +/** + * @brief Convert concentrations to activities and (optionally) map derivatives. + * @details + * When `LOGE_CONCENTRATION == false`, this is a direct pass-through to the activity + * model (`a(c)` and `da/dc`). + * When `LOGE_CONCENTRATION == true`, the input is interpreted as `log(c)` and the + * output as `log(a)`, with derivatives returned as `d log(a) / d log(c)`. + */ +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + typename ACTIVITY_MODEL, + bool LOGE_CONCENTRATION, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE +inline +void calculateActivities( typename ACTIVITY_MODEL::Params const & activityParams, + ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D & activities, + ARRAY_2D & dActivities_dConcentration ) +{ + HPCREACT_UNUSED_VAR( activityParams ); + HPCREACT_UNUSED_VAR( speciesConcentration ); + HPCREACT_UNUSED_VAR( activities ); + HPCREACT_UNUSED_VAR( dActivities_dConcentration ); + HPCREACT_UNUSED_VAR( sizeof( INT_TYPE ) ); + + static constexpr INDEX_TYPE numSpecies = ACTIVITY_MODEL::Params::numSpecies(); + + if constexpr( LOGE_CONCENTRATION ) + { + REAL_TYPE linearConcentration[numSpecies] = { 0.0 }; + REAL_TYPE linearActivities[numSpecies] = { 0.0 }; + REAL_TYPE dLinearActivities_dLinearConcentration[numSpecies][numSpecies] = {{ 0.0 }}; + + for( INDEX_TYPE i = 0; i < numSpecies; ++i ) + { + linearConcentration[i] = exp( speciesConcentration[i] ); + } + + ACTIVITY_MODEL::calculateActivities( activityParams, + linearConcentration, + linearActivities, + dLinearActivities_dLinearConcentration ); + + for( INDEX_TYPE i = 0; i < numSpecies; ++i ) + { + REAL_TYPE const ai = linearActivities[i]; + activities[i] = log( ai ); + + for( INDEX_TYPE j = 0; j < numSpecies; ++j ) + { + dActivities_dConcentration[i][j] = 0.0; + if( ai > 1.0e-300 ) + { + dActivities_dConcentration[i][j] = + dLinearActivities_dLinearConcentration[i][j] * linearConcentration[j] / ai; + } + } + } + } + else + { + ACTIVITY_MODEL::calculateActivities( activityParams, + speciesConcentration, + activities, + dActivities_dConcentration ); + } +} + +} // namespace hpcReact