diff --git a/cmake/Macros.cmake b/cmake/Macros.cmake index d94bd85..faa8256 100644 --- a/cmake/Macros.cmake +++ b/cmake/Macros.cmake @@ -42,6 +42,7 @@ macro(hpcReact_add_code_checks) --suppress=unusedFunction --suppress=constStatement --suppress=unusedStructMember + -DHPCREACT_HOST_DEVICE= -I../hpcReact/src ) if( UNCRUSTIFY_FOUND ) 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..4eae0d3 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 @@ -39,34 +39,46 @@ 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 - { 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0 + { + 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 }, // Reverse rate constants - { 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0 + { + 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 }, + + // Flag of mobile secondary species + { 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 + } }; // *****UNCRUSTIFY-ON****** diff --git a/src/reactions/exampleSystems/unitTests/CMakeLists.txt b/src/reactions/exampleSystems/unitTests/CMakeLists.txt index 99b0388..a5174b1 100644 --- a/src/reactions/exampleSystems/unitTests/CMakeLists.txt +++ b/src/reactions/exampleSystems/unitTests/CMakeLists.txt @@ -2,6 +2,7 @@ set( testSourceFiles testEquilibriumReactions.cpp testKineticReactions.cpp + testMomasEasyCase.cpp ) set( dependencyList hpcReact gtest ) diff --git a/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp new file mode 100644 index 0000000..2f77886 --- /dev/null +++ b/src/reactions/exampleSystems/unitTests/testMomasEasyCase.cpp @@ -0,0 +1,87 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "../MoMasBenchmark.hpp" + +using namespace hpcReact; +using namespace hpcReact::MomMasBenchmark; +using namespace hpcReact::unitTest_utilities; + +//****************************************************************************** + +TEST( testEquilibriumReactions, testMoMasAllEquilibrium ) +{ + + using EquilibriumReactionsType = reactionsSystems::EquilibriumReactions< double, + int, + int >; + + constexpr int numPrimarySpecies = hpcReact::MomMasBenchmark::simpleSystemParams.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.0e-20, // X1 + 0.02, // X2 + 1.0e-20, // X3 + 1.0, // 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::MomMasBenchmark::simpleSystemParams.equilibriumReactionsParameters(), + targetAggregatePrimarySpeciesConcentration, + logInitialPrimarySpeciesConcentration, + logPrimarySpeciesConcentration ); + + double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] = + { + 9.9999999999999919e-21, // X1 + 0.25971841330881928, // X2 + 1.4603613417111526e-24, // X3 + 0.3495378685828045, // X4 + 0.39074371811222675 // S + }; + + for( int r=0; r 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- + 8.55E-03 // 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..d69d517 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/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/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index 87c2819..f8d6b96 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -172,39 +172,78 @@ 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}; + calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + +} + +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; } } - calculateLogSecondarySpeciesConcentration< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations ); - 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 static HPCREACT_HOST_DEVICE void + enforceEquilibrium_LogAggregate( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & speciesConcentration0, + ARRAY_1D & speciesConcentration ); + + + /** + * @brief This method enforces equilibrium for a given set of species using + * log of aggregate primary concentrations. + * @tparam PARAMS_DATA The type of the parameters data. + * @tparam ARRAY_1D The type of the array of species concentrations. + * @tparam ARRAY_1D_TO_CONST The type of the array of species concentrations. + * @param temperature The temperature of the system. + * @param params The parameters for the equilibrium reactions. + * @param targetAggregatePrimarySpeciesConcentration The target aggregate + * primary species concentration. + * @param logPrimarySpeciesConcentration0 The initial value of the log of + * the primary species concentrations. + * @param speciesConcentration The species concentrations to be updated. + * @details This method uses the log of aggregate primary concentrations to enforce + * equilibrium for a given set of species. It uses the + * computeResidualAndJacobianLogAggregate method to compute the residual and + * jacobian for the system and then uses a direct solver to solve the system. + * The solution is then used to update the species concentrations. + */ + 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 & speciesConcentration0, + ARRAY_1D_TO_CONST const & targetAggregatePrimarySpeciesConcentration, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, ARRAY_1D & speciesConcentration ); /** diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index 975dbd7..07e99f3 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -62,6 +62,41 @@ EquilibriumReactions< REAL_TYPE, } } +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE > +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_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(); + double targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 }; + + + + for( int i=0; i @@ -74,6 +109,7 @@ 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 ) { @@ -82,20 +118,21 @@ EquilibriumReactions< REAL_TYPE, 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( 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..c9575ba 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; @@ -83,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; @@ -109,7 +111,9 @@ void timeStepTest( PARAMS_DATA const & params, surfaceArea, logSecondarySpeciesConcentration, aggregatePrimarySpeciesConcentration, + mobileAggregatePrimarySpeciesConcentration, dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration, + dMobileAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration, reactionRates, dReactionRates_dlogPrimarySpeciesConcentration, aggregateSpeciesRates,