From 877df19cc3cb7378e3bdbdc60820f647ec8bb297 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Mon, 14 Jul 2025 11:23:23 -0700 Subject: [PATCH 1/8] 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, From 952ba2ec79965db03ea6156ccf8f951148883774 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Fri, 25 Jul 2025 13:15:12 -0700 Subject: [PATCH 2/8] revised notations in parameters --- .../exampleSystems/MoMasBenchmark.hpp | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 8325882..7ddc72d 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -28,9 +28,9 @@ namespace MomMasBenchmark // 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 + { -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 @@ -69,13 +69,13 @@ namespace MomMasBenchmark }, // Flag of mobile secondary species - { 1, - 1, - 1, - 1, - 1, - 0, - 0 + { 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 } }; From b172c9677489ed3abd1f50e21eb73ecfb1406736 Mon Sep 17 00:00:00 2001 From: jkgolla Date: Mon, 28 Jul 2025 14:13:50 -0700 Subject: [PATCH 3/8] Implement unit test for MoMaS benchmark easy case --- src/reactions/geochemistry/Momas.hpp | 93 ++++++++++++++++++ .../geochemistry/unitTests/CMakeLists.txt | 1 + .../unitTests/testMomasEasyCase.cpp | 95 +++++++++++++++++++ 3 files changed, 189 insertions(+) create mode 100644 src/reactions/geochemistry/Momas.hpp create mode 100644 src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp diff --git a/src/reactions/geochemistry/Momas.hpp b/src/reactions/geochemistry/Momas.hpp new file mode 100644 index 0000000..252bb9e --- /dev/null +++ b/src/reactions/geochemistry/Momas.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 "../reactionsSystems/Parameters.hpp" + +namespace hpcReact +{ + +namespace geochemistry +{ +// turn off uncrustify to allow for better readability of the parameters +// *****UNCRUSTIFY-OFF****** + +namespace momas +{ + +constexpr CArrayWrapper stoichMatrix = + { + // 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 equilibriumConstants = + { + 1.0e-12, // C1 = -X2 + 1.0, // C2 = X2 + X3 + 1.0, // C3 = X2 + X4 + 0.1, // C4 = -4X2 + X3 + 3X4 + 1.0e35, // C5 = 4X2 + 3X3 + X4 + 1.0e6, // CS1 = 3X2 + X3 + S + 1.0e-1 // CS2 = -3X2 + X4 + 2S + }; + +constexpr CArrayWrapper forwardRates = + { + 0.0, // C1 = -X2 + 0.0, // C2 = X2 + X3 + 0.0, // C3 = X2 + X4 + 0.0, // C4 = -4X2 + X3 + 3X4 + 0.0, // C5 = 4X2 + 3X3 + X4 + 0.0, // CS1 = 3X2 + X3 + S + 0.0 // CS2 = -3X2 + X4 + 2S + }; + +constexpr CArrayWrapper reverseRates = + { + 0.0, // C1 = -X2 + 0.0, // C2 = X2 + X3 + 0.0, // C3 = X2 + X4 + 0.0, // C4 = -4X2 + X3 + 3X4 + 0.0, // C5 = 4X2 + 3X3 + X4 + 0.0, // CS1 = 3X2 + X3 + S + 0.0 // CS2 = -3X2 + X4 + 2S + }; + +constexpr CArrayWrapper mobileSpeciesFlag = + { + 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 + }; + } + using momasSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 0 >; + using momasSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >; + using momasSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 12, 7, 7 >; + + constexpr momasSystemAllKineticType momasSystemAllKinetic( momas::stoichMatrix, momas::equilibriumConstants, momas::forwardRates, momas::reverseRates, momas::mobileSpeciesFlag ); + constexpr momasSystemAllEquilibriumType momasSystemAllEquilibrium( momas::stoichMatrix, momas::equilibriumConstants, momas::forwardRates, momas::reverseRates, momas::mobileSpeciesFlag ); + constexpr momasSystemType carbonateSystem( momas::stoichMatrix, momas::equilibriumConstants, momas::forwardRates, momas::reverseRates, momas::mobileSpeciesFlag ); + +// *****UNCRUSTIFY-ON****** +} // namespace geochemistry +} // namespace hpcReact diff --git a/src/reactions/geochemistry/unitTests/CMakeLists.txt b/src/reactions/geochemistry/unitTests/CMakeLists.txt index a75b366..0827a05 100644 --- a/src/reactions/geochemistry/unitTests/CMakeLists.txt +++ b/src/reactions/geochemistry/unitTests/CMakeLists.txt @@ -3,6 +3,7 @@ set( testSourceFiles testGeochemicalEquilibriumReactions.cpp testGeochemicalKineticReactions.cpp testGeochemicalMixedReactions.cpp + testMomasEasyCase.cpp ) set( dependencyList hpcReact gtest ) diff --git a/src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp b/src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp new file mode 100644 index 0000000..d8bf5c4 --- /dev/null +++ b/src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp @@ -0,0 +1,95 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: (BSD-3-Clause) + * + * Copyright (c) 2025- Lawrence Livermore National Security LLC + * All rights reserved + * + * See top level LICENSE files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +#include "reactions/unitTestUtilities/equilibriumReactionsTestUtilities.hpp" +#include "../Momas.hpp" + +using namespace hpcReact; +using namespace hpcReact::geochemistry; +using namespace hpcReact::unitTest_utilities; + + +// //****************************************************************************** +// TEST( testEquilibriumReactions, testEnforceEquilibrium ) +// { +// double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; +// double const expectedSpeciesConcentrations[5] = { 3.92138294e-01, 3.03930853e-01, 5.05945481e-01, 7.02014628e-01, 5.95970745e-01 }; + + +// std::cout<<" RESIDUAL_FORM 2:"<( simpleTestRateParams.equilibriumReactionsParameters(), +// initialSpeciesConcentration, +// expectedSpeciesConcentrations ); + +// } + + +//****************************************************************************** + +TEST( testEquilibriumReactions, testMoMasAllEquilibrium ) +{ + + using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, + int, + int >; + + constexpr int numPrimarySpecies = hpcReact::geochemistry::momasSystemAllEquilibrium.numPrimarySpecies(); + + double const initialPrimarySpeciesConcentration[numPrimarySpecies] = + { + 1.00e-20, // X1 + 2.00, // X2 + 1.00e-20, // X3 + 2.00, // X4 + 1.00 // S + }; + + + + double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] = + { + log( initialPrimarySpeciesConcentration[0] ), + log( initialPrimarySpeciesConcentration[1] ), + log( initialPrimarySpeciesConcentration[2] ), + log( initialPrimarySpeciesConcentration[3] ), + log( initialPrimarySpeciesConcentration[4] ) + }; + + double logPrimarySpeciesConcentration[numPrimarySpecies]; + EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0, + hpcReact::geochemistry::momasSystemAllEquilibrium.equilibriumReactionsParameters(), + logInitialPrimarySpeciesConcentration, + logPrimarySpeciesConcentration ); + + double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] = + { + 9.7051090442170804E-21, // X1 + 5.0023298955833342E-12, // X2 + 1.9327372426296357E-33, // X3 + 7.3929274619958745E-12, // X4 + 9.8708294125907346E-13 // S + }; + + for( int r=0; r Date: Tue, 29 Jul 2025 16:15:14 -0700 Subject: [PATCH 4/8] cleanup of massActions...and change definition of K in the momas example --- src/reactions/geochemistry/Momas.hpp | 14 +++---- src/reactions/massActions/MassActions.hpp | 39 ++++--------------- ...ionsAggregatePrimaryConcentration_impl.hpp | 19 ++++++++- 3 files changed, 32 insertions(+), 40 deletions(-) diff --git a/src/reactions/geochemistry/Momas.hpp b/src/reactions/geochemistry/Momas.hpp index 252bb9e..d57cb0e 100644 --- a/src/reactions/geochemistry/Momas.hpp +++ b/src/reactions/geochemistry/Momas.hpp @@ -38,13 +38,13 @@ constexpr CArrayWrapper stoichMatrix = constexpr CArrayWrapper equilibriumConstants = { - 1.0e-12, // C1 = -X2 - 1.0, // C2 = X2 + X3 - 1.0, // C3 = X2 + X4 - 0.1, // C4 = -4X2 + X3 + 3X4 - 1.0e35, // C5 = 4X2 + 3X3 + X4 - 1.0e6, // CS1 = 3X2 + X3 + S - 1.0e-1 // CS2 = -3X2 + X4 + 2S + 1.0e12, // C1 + X2 = ?? + 1.0, // C2 = X2 + X3 + 1.0, // C3 = X2 + X4 + 1.0e1, // C4 + 4X2 = X3 + 3X4 + 1.0e-35, // C5 = 4X2 + 3X3 + X4 + 1.0e-6, // CS1 = 3X2 + X3 + S + 1.0e1 // CS2 + 3X2 = + X4 + 2S }; constexpr CArrayWrapper forwardRates = diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index 28d1d4f..61a39cb 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -172,42 +172,17 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ) { constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); - constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0}; - for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) - { - for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) - { - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0; - } - } - - calculateLogSecondarySpeciesConcentration< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations ); + calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE>( params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); - 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 Date: Wed, 30 Jul 2025 05:49:23 -0700 Subject: [PATCH 5/8] add new interface function to specify target aggregate concentrations for equilibrium. Fix momas test --- .../unitTests/testMomasEasyCase.cpp | 20 ++++++--- .../reactionsSystems/EquilibriumReactions.hpp | 11 +++++ ...ionsAggregatePrimaryConcentration_impl.hpp | 42 +++++++++++++++++-- 3 files changed, 63 insertions(+), 10 deletions(-) diff --git a/src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp b/src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp index d8bf5c4..bbf419e 100644 --- a/src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp +++ b/src/reactions/geochemistry/unitTests/testMomasEasyCase.cpp @@ -44,17 +44,24 @@ TEST( testEquilibriumReactions, testMoMasAllEquilibrium ) constexpr int numPrimarySpecies = hpcReact::geochemistry::momasSystemAllEquilibrium.numPrimarySpecies(); + double const targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] = + { + 1.0e-20, // X1 + -2.0, // X2 + 1.0e-20, // X3 + 2.0, // X4 + 1.0 // S + }; + double const initialPrimarySpeciesConcentration[numPrimarySpecies] = { - 1.00e-20, // X1 - 2.00, // X2 - 1.00e-20, // X3 - 2.00, // X4 + 1.0e-20, // X1 + 0.02, // X2 + 1.0e-20, // X3 + 1.0, // X4 1.00 // S }; - - double const logInitialPrimarySpeciesConcentration[numPrimarySpecies] = { log( initialPrimarySpeciesConcentration[0] ), @@ -67,6 +74,7 @@ TEST( testEquilibriumReactions, testMoMasAllEquilibrium ) double logPrimarySpeciesConcentration[numPrimarySpecies]; EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0, hpcReact::geochemistry::momasSystemAllEquilibrium.equilibriumReactionsParameters(), + targetAggregatePrimarySpeciesConcentration, logInitialPrimarySpeciesConcentration, logPrimarySpeciesConcentration ); diff --git a/src/reactions/reactionsSystems/EquilibriumReactions.hpp b/src/reactions/reactionsSystems/EquilibriumReactions.hpp index b1a7eb7..c8d52c5 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactions.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactions.hpp @@ -97,6 +97,17 @@ class EquilibriumReactions ARRAY_1D_TO_CONST const & speciesConcentration0, ARRAY_1D & speciesConcentration ); + template< typename PARAMS_DATA, + typename ARRAY_1D, + typename ARRAY_1D_TO_CONST > + static HPCREACT_HOST_DEVICE + void + enforceEquilibrium_Aggregate( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & targetAggregatePrimarySpeciesConcentration, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, + ARRAY_1D & speciesConcentration ); + /** * @brief This method computes the residual and jacobian when using reaction extents to solve * for the equilibrium of a given set of species. diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index 3c2eb7f..a58cc28 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -79,17 +79,50 @@ EquilibriumReactions< REAL_TYPE, { HPCREACT_UNUSED_VAR( temperature ); constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + double targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 }; + + + + for( int i=0; i +template< typename PARAMS_DATA, + typename ARRAY_1D, + typename ARRAY_1D_TO_CONST > +HPCREACT_HOST_DEVICE inline +void +EquilibriumReactions< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >::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 ) +{ + HPCREACT_UNUSED_VAR( temperature ); + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); double residual[numPrimarySpecies] = { 0.0 }; // double aggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 }; - double targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 }; double dLogCp[numPrimarySpecies] = { 0.0 }; CArrayWrapper< double, numPrimarySpecies, numPrimarySpecies > jacobian; - for( int i=0; i Date: Tue, 5 Aug 2025 16:06:25 -0700 Subject: [PATCH 6/8] a few changes to make GEOS work --- src/reactions/exampleSystems/MoMasBenchmark.hpp | 14 +++++++------- .../testGeochemicalEquilibriumReactions.cpp | 8 ++++---- .../reactionsSystems/EquilibriumReactions.hpp | 8 ++++---- ...ReactionsAggregatePrimaryConcentration_impl.hpp | 10 +++++----- .../mixedReactionsTestUtilities.hpp | 8 ++++---- 5 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/reactions/exampleSystems/MoMasBenchmark.hpp b/src/reactions/exampleSystems/MoMasBenchmark.hpp index 7ddc72d..f39410c 100644 --- a/src/reactions/exampleSystems/MoMasBenchmark.hpp +++ b/src/reactions/exampleSystems/MoMasBenchmark.hpp @@ -39,13 +39,13 @@ namespace MomMasBenchmark // Equilibrium constants K { - 1.0e-12, // C1 - 1.0, // C2 - 1.0, // C3 - 0.1, // C4 - 1.0e35, // C5 - 1.0e6, // CS1 - 1.0e-1 // CS2 + 1.0e12, // C1 + X2 = ?? + 1.0, // C2 = X2 + X3 + 1.0, // C3 = X2 + X4 + 1.0e1, // C4 + 4X2 = X3 + 3X4 + 1.0e-35, // C5 = 4X2 + 3X3 + X4 + 1.0e-6, // CS1 = 3X2 + X3 + S + 1.0e1 // CS2 + 3X2 = + X4 + 2S }, // Forward rate constants diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp index b3a6eeb..1343c7f 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp @@ -133,10 +133,10 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 ) }; double logPrimarySpeciesConcentration[numPrimarySpecies]; - EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0, - hpcReact::geochemistry::carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), - logInitialPrimarySpeciesConcentration, - logPrimarySpeciesConcentration ); + EquilibriumReactionsType::enforceEquilibrium_LogAggregate( 0, + hpcReact::geochemistry::carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), + logInitialPrimarySpeciesConcentration, + logPrimarySpeciesConcentration ); double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] = { diff --git a/src/reactions/reactionsSystems/EquilibriumReactions.hpp b/src/reactions/reactionsSystems/EquilibriumReactions.hpp index c8d52c5..da24f44 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactions.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactions.hpp @@ -92,10 +92,10 @@ class EquilibriumReactions typename ARRAY_1D_TO_CONST > static HPCREACT_HOST_DEVICE void - enforceEquilibrium_Aggregate( RealType const & temperature, - PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & speciesConcentration0, - ARRAY_1D & speciesConcentration ); + enforceEquilibrium_LogAggregate( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & speciesConcentration0, + ARRAY_1D & speciesConcentration ); template< typename PARAMS_DATA, typename ARRAY_1D, diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index a58cc28..a90e038 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -72,10 +72,10 @@ HPCREACT_HOST_DEVICE inline void EquilibriumReactions< REAL_TYPE, INT_TYPE, - INDEX_TYPE >::enforceEquilibrium_Aggregate( REAL_TYPE const & temperature, - PARAMS_DATA const & params, - ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, - ARRAY_1D & logPrimarySpeciesConcentration ) + INDEX_TYPE >::enforceEquilibrium_LogAggregate( REAL_TYPE const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, + ARRAY_1D & logPrimarySpeciesConcentration ) { HPCREACT_UNUSED_VAR( temperature ); constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); @@ -131,7 +131,7 @@ EquilibriumReactions< REAL_TYPE, // 0: 1e-20 -0 2 -2.5e+11 1e-20 7 2 1.8 1 5 printf( "iter X1 R0 X2 R1 X3 R2 X4 R3 S R4\n"); printf( "---- --------------- --------------- --------------- --------------- ---------------\n"); - for( int k=0; k<40; ++k ) + for( int k=0; k<80; ++k ) { computeResidualAndJacobianAggregatePrimaryConcentrations( temperature, params, diff --git a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp index 00e211d..c9575ba 100644 --- a/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/mixedReactionsTestUtilities.hpp @@ -85,10 +85,10 @@ void timeStepTest( PARAMS_DATA const & params, aggregatePrimarySpeciesConcentration[i] = speciesConcentration[i]; } - EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature, - params.equilibriumReactionsParameters(), - logPrimarySpeciesConcentration, - logPrimarySpeciesConcentration ); + EquilibriumReactionsType::enforceEquilibrium_LogAggregate( temperature, + params.equilibriumReactionsParameters(), + logPrimarySpeciesConcentration, + logPrimarySpeciesConcentration ); /// Time step loop double time = 0.0; From 3a687e52f6489ac322246c83a3488a76edbed683 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Wed, 13 Aug 2025 13:07:55 -0700 Subject: [PATCH 7/8] more iterations to solve equilibrium --- .../EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index a90e038..c47f309 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -131,7 +131,7 @@ EquilibriumReactions< REAL_TYPE, // 0: 1e-20 -0 2 -2.5e+11 1e-20 7 2 1.8 1 5 printf( "iter X1 R0 X2 R1 X3 R2 X4 R3 S R4\n"); printf( "---- --------------- --------------- --------------- --------------- ---------------\n"); - for( int k=0; k<80; ++k ) + for( int k=0; k<150; ++k ) { computeResidualAndJacobianAggregatePrimaryConcentrations( temperature, params, From 6d0282aa944e8fb2efbe0abe8fd1acf5e2e645a1 Mon Sep 17 00:00:00 2001 From: frankfeifan Date: Thu, 14 Aug 2025 17:28:29 -0700 Subject: [PATCH 8/8] fix gcc error --- src/reactions/geochemistry/Ultramafics.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactions/geochemistry/Ultramafics.hpp b/src/reactions/geochemistry/Ultramafics.hpp index e4d6994..d69d517 100644 --- a/src/reactions/geochemistry/Ultramafics.hpp +++ b/src/reactions/geochemistry/Ultramafics.hpp @@ -151,7 +151,7 @@ constexpr CArrayWrapper mobileSpeciesFlag = 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 >;