From 877df19cc3cb7378e3bdbdc60820f647ec8bb297 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Mon, 14 Jul 2025 11:23:23 -0700 Subject: [PATCH] initial commit --- src/reactions/exampleSystems/BulkGeneric.hpp | 8 ++- .../exampleSystems/MoMasBenchmark.hpp | 10 +++ src/reactions/geochemistry/Carbonate.hpp | 26 ++++++-- src/reactions/geochemistry/Forge.hpp | 25 +++++++- src/reactions/geochemistry/Ultramafics.hpp | 31 ++++++++- src/reactions/massActions/MassActions.hpp | 63 +++++++++++++++++++ .../MixedEquilibriumKineticReactions.hpp | 10 +++ .../MixedEquilibriumKineticReactions_impl.hpp | 19 +++--- src/reactions/reactionsSystems/Parameters.hpp | 19 ++++-- .../mixedReactionsTestUtilities.hpp | 4 ++ 10 files changed, 190 insertions(+), 25 deletions(-) diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp index 99288ac..036a6cd 100644 --- a/src/reactions/exampleSystems/BulkGeneric.hpp +++ b/src/reactions/exampleSystems/BulkGeneric.hpp @@ -61,7 +61,9 @@ simpleKineticTestRateParams = // forward rate constants { 1.0, 0.5 }, // reverse rate constants - { 1.0, 0.5 } + { 1.0, 0.5 }, + // flag of mobile secondary species + { 1, 1 } }; using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >; @@ -80,7 +82,9 @@ simpleTestRateParams = // forward rate constants { 1.0, 0.5 }, // reverse rate constants - { 1.0, 0.5 } + { 1.0, 0.5 }, + // flag of mobile secondary species + { 1, 1 } }; // *****UNCRUSTIFY-ON****** diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 46053e0..8325882 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -67,6 +67,16 @@ namespace MomMasBenchmark 0.0, 0.0 }, + + // Flag of mobile secondary species + { 1, + 1, + 1, + 1, + 1, + 0, + 0 + } }; // *****UNCRUSTIFY-ON****** diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp index 137903c..1662677 100644 --- a/src/reactions/geochemistry/Carbonate.hpp +++ b/src/reactions/geochemistry/Carbonate.hpp @@ -69,7 +69,7 @@ constexpr CArrayWrapper equilibriumConstants = 3.98E+00, // CaCl2 = Ca+2 + 2Cl- 3.72E-03, // MgSO4 = Mg+2 + SO4-2 1.51E-01, // NaSO4- = Na+ + SO4-2 - 1.17E+07 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + 70.55 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) // 1 }; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) @@ -85,7 +85,7 @@ constexpr CArrayWrapper forwardRates = 1.0e7, // CaCl2 = Ca+2 + 2Cl- 1.0e5, // MgSO4 = Mg+2 + SO4-2 1.0e7, // NaSO4- = Na+ + SO4-2 - 1.0e5 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + 1.55e-6 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) // 1 }; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) @@ -101,17 +101,31 @@ constexpr CArrayWrapper reverseRates = 2.51E+06, // CaCl2 = Ca+2 + 2Cl- 2.69E+07, // MgSO4 = Mg+2 + SO4-2 6.62E+07, // NaSO4- = Na+ + SO4-2 - 8.55E-03 // CaCO3 + H+ = Ca+2 + HCO3- + 2.197e-8 // CaCO3 + H+ = Ca+2 + HCO3- // 1 // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) }; + +constexpr CArrayWrapper mobileSpeciesFlag = + { 1, // OH- + H+ = H2O + 1, // CO2 + H2O = H+ + HCO3- + 1, // CO3-2 + H+ = HCO3- + 1, // H2CO3 = H+ + HCO3- + 1, // CaHCO3+ = Ca+2 + HCO3- + 1, // CaSO4 = Ca+2 + SO4-2 + 1, // CaCl+ = Ca+2 + Cl- + 1, // CaCl2 = Ca+2 + 2Cl- + 1, // MgSO4 = Mg+2 + SO4-2 + 1, // NaSO4- = Na+ + SO4-2 + 1 // CaCO3 + H+ = Ca+2 + HCO3- + }; } using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 0 >; using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 11 >; using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 11, 10 >; - constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates ); - constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates ); - constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates ); + constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag ); + 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 ); // *****UNCRUSTIFY-ON****** } // namespace geochemistry diff --git a/src/reactions/geochemistry/Forge.hpp b/src/reactions/geochemistry/Forge.hpp index 618a6fb..276577c 100644 --- a/src/reactions/geochemistry/Forge.hpp +++ b/src/reactions/geochemistry/Forge.hpp @@ -119,12 +119,35 @@ constexpr CArrayWrapper< double, 19 > reverseRateConstant = 1.0 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq) }; +constexpr CArrayWrapper< int, 19 > mobileSpeciesFlag = +{ + 1, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻ + 1, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻ + 1, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻ + 1, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻ + 1, // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source) + 1, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻ + 1, // MgCO₃(aq) + H⁺ ⇌ Mg²⁺ + HCO₃⁻ + 1, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻ + 1, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻ + 1, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻ + 1, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻ + 1, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq) + 1, // NaHSilO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq) + 1, // NaCl ⇌ Na⁺ + Cl⁻ + 1, // KCl ⇌ K⁺ + Cl⁻ + 1, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻ + 1, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻ + 1, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq) + 1 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq) +}; + } using forgeSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 26, 19, 16 >; -constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant ); +constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant, forge::mobileSpeciesFlag ); // *****UNCRUSTIFY-ON****** diff --git a/src/reactions/geochemistry/Ultramafics.hpp b/src/reactions/geochemistry/Ultramafics.hpp index 15e9c63..e4d6994 100644 --- a/src/reactions/geochemistry/Ultramafics.hpp +++ b/src/reactions/geochemistry/Ultramafics.hpp @@ -126,15 +126,40 @@ constexpr CArrayWrapper reverseRates = 2.83E-44, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O 2.10E-25 // Mg(OH)2 + 2H+ = Mg++ + 2H2O }; + +constexpr CArrayWrapper mobileSpeciesFlag = + { + 1, // OH- + H+ = H2O + 1, // CO2(aq) + H2O = HCO3- + H+ + 1, // CO3-- + H+ = HCO3- + 1, // Mg2OH+++ + H+ = 2Mg++ + H2O + 1, // Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O + 1, // MgOH+ + H+ = Mg++ + H2O + 1, // Mg2CO3++ + H+ = 2Mg++ + HCO3- + 1, // MgCO3 + H+ = Mg++ + HCO3- + 1, // MgHCO3+ = Mg++ + HCO3- + 1, // Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O + 1, // MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O + 1, // MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O + 1, // H2SiO4-- + 2H+ = SiO2(aq) + 2H2O + 1, // H3SiO4- + H+ = SiO2(aq) + 2H2O + 1, // H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O + 1, // H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O + 1, // Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O + 1, // MgCO3 + H+ = Mg++ + HCO3- + 1, // SiO2 = SiO2(aq) + 1, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O + 1 // Mg(OH)2 + 2H+ = Mg++ + 2H2O + }; }; using ultramaficSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 0 >; using ultramaficSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 21 >; using ultramaficSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 16 >; - constexpr ultramaficSystemAllKineticType ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates ); - constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates ); - constexpr ultramaficSystemType ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates ); + constexpr ultramaficSystemAllKineticType ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag ); + constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag ); + constexpr ultramaficSystemType ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag ); // *****UNCRUSTIFY-ON****** } // namespace geochemistry diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index 87c2819..28d1d4f 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -210,5 +210,68 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( 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_PRIMARY, + typename ARRAY_1D_SECONDARY, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE +inline +void calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_1D_PRIMARY & mobileAggregatePrimarySpeciesConcentrations, + ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations, + ARRAY_2D & dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) +{ + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + + calculateLogSecondarySpeciesConcentration< 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 ) + { + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; + dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.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; + + 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( params.equilibriumReactionsParameters(), - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations, - aggregatePrimarySpeciesConcentrations, - dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + massActions::calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params.equilibriumReactionsParameters(), + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + mobileAggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, + dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); if constexpr( PARAMS_DATA::numKineticReactions() > 0 ) { diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index ecfd827..d2ca7ba 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -54,17 +54,21 @@ struct EquilibriumReactionsParameters constexpr EquilibriumReactionsParameters( CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > const & stoichiometricMatrix, - CArrayWrapper< RealType, NUM_REACTIONS > equilibriumConstant ): + CArrayWrapper< RealType, NUM_REACTIONS > equilibriumConstant, + CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag ): m_stoichiometricMatrix( stoichiometricMatrix ), - m_equilibriumConstant( equilibriumConstant ) + m_equilibriumConstant( equilibriumConstant ), + m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ) {} RealType stoichiometricMatrix( IndexType const r, int const i ) const { return m_stoichiometricMatrix[r][i]; } RealType equilibriumConstant( IndexType const r ) const { return m_equilibriumConstant[r]; } + IntType mobileSecondarySpeciesFlag( IndexType const r ) const { return m_mobileSecondarySpeciesFlag[r]; } CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; CArrayWrapper< RealType, NUM_REACTIONS > m_equilibriumConstant; + CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag; }; template< typename REAL_TYPE, @@ -124,11 +128,13 @@ struct MixedReactionsParameters constexpr MixedReactionsParameters( CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > const & stoichiometricMatrix, CArrayWrapper< RealType, NUM_REACTIONS > const & equilibriumConstant, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, - CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse ): + CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse, + CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag ): m_stoichiometricMatrix( stoichiometricMatrix ), m_equilibriumConstant( equilibriumConstant ), m_rateConstantForward( rateConstantForward ), - m_rateConstantReverse( rateConstantReverse ) + m_rateConstantReverse( rateConstantReverse ), + m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ) {} static constexpr IndexType numReactions() { return NUM_REACTIONS; } @@ -149,6 +155,7 @@ struct MixedReactionsParameters { CArrayWrapper< RealType, numEquilibriumReactions(), numSpecies() > eqMatrix{}; CArrayWrapper< RealType, numEquilibriumReactions() > eqConstants{}; + CArrayWrapper< IntType, numEquilibriumReactions() > mobileSpeciesFlags{}; for( IntType i = 0; i < numEquilibriumReactions(); ++i ) { @@ -157,9 +164,10 @@ struct MixedReactionsParameters eqMatrix( i, j ) = m_stoichiometricMatrix( i, j ); } eqConstants( i ) = m_equilibriumConstant( i ); + mobileSpeciesFlags( i ) = m_mobileSecondarySpeciesFlag( i ); } - return { eqMatrix, eqConstants }; + return { eqMatrix, eqConstants, mobileSpeciesFlags }; } constexpr @@ -229,6 +237,7 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > m_equilibriumConstant; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; + CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag; }; diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp index 0888c7f..00e211d 100644 --- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp @@ -70,7 +70,9 @@ void timeStepTest( PARAMS_DATA const & params, CArrayWrapper< REAL_TYPE, numSecondarySpecies > logSecondarySpeciesConcentration; CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration; CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration_n; + CArrayWrapper< REAL_TYPE, numPrimarySpecies > mobileAggregatePrimarySpeciesConcentration; CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration; + CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dMobileAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration; CArrayWrapper< REAL_TYPE, numKineticReactions > reactionRates; CArrayWrapper< REAL_TYPE, numKineticReactions, numPrimarySpecies > dReactionRates_dlogPrimarySpeciesConcentration; CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregateSpeciesRates; @@ -109,7 +111,9 @@ void timeStepTest( PARAMS_DATA const & params, surfaceArea, logSecondarySpeciesConcentration, aggregatePrimarySpeciesConcentration, + mobileAggregatePrimarySpeciesConcentration, dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration, + dMobileAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration, reactionRates, dReactionRates_dlogPrimarySpeciesConcentration, aggregateSpeciesRates,