diff --git a/src/reactions/exampleSystems/BulkGeneric.hpp b/src/reactions/exampleSystems/BulkGeneric.hpp index 036a6cd..77547ce 100644 --- a/src/reactions/exampleSystems/BulkGeneric.hpp +++ b/src/reactions/exampleSystems/BulkGeneric.hpp @@ -63,7 +63,9 @@ simpleKineticTestRateParams = // reverse rate constants { 1.0, 0.5 }, // flag of mobile secondary species - { 1, 1 } + { 1, 1 }, + // Use the forward and reverse to calculate the kinetic reaction rates + 0 }; using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >; @@ -84,7 +86,9 @@ simpleTestRateParams = // reverse rate constants { 1.0, 0.5 }, // flag of mobile secondary species - { 1, 1 } + { 1, 1 }, + // Use the forward and reverse to calculate the kinetic reaction rates + 0 }; // *****UNCRUSTIFY-ON****** diff --git a/src/reactions/exampleSystems/ChainGeneric.hpp b/src/reactions/exampleSystems/ChainGeneric.hpp new file mode 100644 index 0000000..8b0acd8 --- /dev/null +++ b/src/reactions/exampleSystems/ChainGeneric.hpp @@ -0,0 +1,68 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ChainGeneric +{ +// *****UNCRUSTIFY-OFF****** + + using serialAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 3, 3, 0 >; + + constexpr serialAllKineticType serialAllKineticParams = + { + // Stoichiometric matrix [3 rows × 3 columns] + // Columns 0–3 are primary species: {C1, C2, C3 } + { + // C1 C2 C3 + { -1, 1, 0 }, // C1 = C2 + { 0, -1, 1 }, // C2 = 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 + }; + +// *****UNCRUSTIFY-ON****** +} // namespace ChainGeneric +} // namespace hpcReact diff --git a/src/reactions/exampleSystems/unitTests/CMakeLists.txt b/src/reactions/exampleSystems/unitTests/CMakeLists.txt index bbc8bab..6892c69 100644 --- a/src/reactions/exampleSystems/unitTests/CMakeLists.txt +++ b/src/reactions/exampleSystems/unitTests/CMakeLists.txt @@ -4,6 +4,7 @@ set( testSourceFiles testKineticReactions.cpp testMomasEasyCase.cpp testMomasMediumCase.cpp + testGenericChainReactions.cpp ) set( dependencyList hpcReact gtest ) diff --git a/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp new file mode 100644 index 0000000..6e0f9a5 --- /dev/null +++ b/src/reactions/exampleSystems/unitTests/testGenericChainReactions.cpp @@ -0,0 +1,64 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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/kineticReactionsTestUtilities.hpp" +#include "../ChainGeneric.hpp" + + +using namespace hpcReact; +using namespace hpcReact::unitTest_utilities; + +//****************************************************************************** +TEST( testChainGenericKineticReactions, computeReactionRatesTest_chainReactionParams ) +{ + using namespace hpcReact::ChainGeneric; + + double const initialSpeciesConcentration[] = + { + 1.0, // C1 + 1e-8, // C2 + 1e-8 // C3 + }; + double const surfaceArea[] = + { + 0.0, // C1 -> C2 + 0.0, // C2 -> C3 + 0.0 // C3 -> + }; + + double const expectedReactionRates[] = + { + 0.05, // C1 -> C2 + 3e-10, // C2 -> C3 + 2e-10 // C3 -> + }; + double const expectedReactionRatesDerivatives[][3] = + { { 0.05, 0.0, 0.0 }, + { 0.0, 0.03, 0.0 }, + { 0.0, 0.0, 0.02 } }; + computeReactionRatesTest< double, false >( serialAllKineticParams.kineticReactionsParameters(), + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); + computeReactionRatesTest< double, true >( serialAllKineticParams.kineticReactionsParameters(), + initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here + expectedReactionRates, + expectedReactionRatesDerivatives ); +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} diff --git a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp index 7baeb0e..1286169 100644 --- a/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp +++ b/src/reactions/exampleSystems/unitTests/testKineticReactions.cpp @@ -24,16 +24,19 @@ using namespace hpcReact::unitTest_utilities; TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams ) { 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.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); computeReactionRatesTest< double, true >( bulkGeneric::simpleKineticTestRateParams.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); } diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp index f36fe0a..52ae3ab 100644 --- a/src/reactions/geochemistry/Carbonate.hpp +++ b/src/reactions/geochemistry/Carbonate.hpp @@ -113,7 +113,7 @@ using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParame using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 10, 10 >; using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 16, 10, 9 >; -constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag ); +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 ); diff --git a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp index 6df4950..8c4e116 100644 --- a/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp +++ b/src/reactions/geochemistry/unitTests/testGeochemicalKineticReactions.cpp @@ -43,6 +43,18 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic ) 1.09 // Na+1 }; + double const surfaceArea[10] = { 0.0, // OH- + H+ = H2O + 0.0, // CO2 + H2O = H+ + HCO3- + 0.0, // CO3-2 + H+ = HCO3- + 0.0, // CaHCO3+ = Ca+2 + HCO3- + 0.0, // CaSO4 = Ca+2 + SO4-2 + 0.0, // CaCl+ = Ca+2 + Cl- + 0.0, // CaCl2 = Ca+2 + 2Cl- + 0.0, // MgSO4 = Mg+2 + SO4-2 + 0.0, // NaSO4- = Na+ + SO4-2 + 0.0, // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + }; + double const expectedReactionRates[10] = { -0.001424736, // OH- + H+ = H2O -12610.7392, // CO2 + H2O = H+ + HCO3- -0.175591624, // CO3-2 + H+ = HCO3- @@ -71,15 +83,58 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic ) computeReactionRatesTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); computeReactionRatesTest< double, true >( carbonateSystemAllKinetic.kineticReactionsParameters(), initialSpeciesConcentration, + surfaceArea, // No use. Just to pass something here expectedReactionRates, expectedReactionRatesDerivatives ); } +TEST( testKineticReactions, computeReactionRatesQuotientTest_carbonateSystem ) +{ + double const initialSpeciesConcentration[16] = + { + 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- + 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 surfaceArea[1] = { 1e6 }; // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + + double const expectedReactionRates[1] = { 1.5488389999999999 }; // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + + double const expectedReactionRatesDerivatives[1][16] = + { + { 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(), + initialSpeciesConcentration, + surfaceArea, + expectedReactionRates, + expectedReactionRatesDerivatives ); + computeReactionRatesTest< double, true >( carbonateSystem.kineticReactionsParameters(), + initialSpeciesConcentration, + surfaceArea, + expectedReactionRates, + expectedReactionRatesDerivatives ); +} //****************************************************************************** diff --git a/src/reactions/massActions/MassActions.hpp b/src/reactions/massActions/MassActions.hpp index d579809..17dbbf8 100644 --- a/src/reactions/massActions/MassActions.hpp +++ b/src/reactions/massActions/MassActions.hpp @@ -72,6 +72,11 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, ARRAY_1D & logSecondarySpeciesConcentrations ) { + if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) + { + return; + } + massActions_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params, @@ -173,15 +178,22 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, { static constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); - REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0}; + if constexpr( numSecondarySpecies > 0 ) + { + REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0}; - calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params, - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations, - aggregatePrimarySpeciesConcentrations, - dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } + else + { + GEOS_UNUSED_VAR( logPrimarySpeciesConcentrations, aggregatePrimarySpeciesConcentrations, dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations ); + } } diff --git a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index 545c7ad..5dad38c 100644 --- a/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -113,6 +113,11 @@ EquilibriumReactions< REAL_TYPE, ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0, ARRAY_1D & logPrimarySpeciesConcentration ) { + if constexpr( PARAMS_DATA::numSecondarySpecies() <= 0 ) + { + return; + } + HPCREACT_UNUSED_VAR( temperature ); static constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); diff --git a/src/reactions/reactionsSystems/KineticReactions.hpp b/src/reactions/reactionsSystems/KineticReactions.hpp index d66f1e9..a5f5575 100644 --- a/src/reactions/reactionsSystems/KineticReactions.hpp +++ b/src/reactions/reactionsSystems/KineticReactions.hpp @@ -13,6 +13,8 @@ #include "common/macros.hpp" +#include + /** @file KineticReactions.hpp * @brief Header file for the KineticReactions class. * @author HPC-REACT Team @@ -130,12 +132,23 @@ class KineticReactions ARRAY_1D & reactionRates, ARRAY_2D & reactionRatesDerivatives ) { - computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature, - params, - speciesConcentration, - surfaceArea, - reactionRates, - reactionRatesDerivatives ); + if( params.reactionRatesUpdateOption() == 0 ) + { + computeReactionRates_impl< PARAMS_DATA, true >( temperature, + params, + speciesConcentration, + reactionRates, + reactionRatesDerivatives ); + } + else if( params.reactionRatesUpdateOption() == 1 ) + { + computeReactionRatesQuotient_impl< PARAMS_DATA, true >( temperature, + params, + speciesConcentration, + surfaceArea, + reactionRates, + reactionRatesDerivatives ); + } } diff --git a/src/reactions/reactionsSystems/KineticReactions_impl.hpp b/src/reactions/reactionsSystems/KineticReactions_impl.hpp index f421f8b..b38a69e 100644 --- a/src/reactions/reactionsSystems/KineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/KineticReactions_impl.hpp @@ -271,7 +271,7 @@ KineticReactions< REAL_TYPE, RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], s_ri ) : 0.0; quotient *= productTerm_i; } - + if constexpr( CALCULATE_DERIVATIVES ) { for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) diff --git a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp index f6dfd18..f24e4d6 100644 --- a/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp +++ b/src/reactions/reactionsSystems/MixedEquilibriumKineticReactions_impl.hpp @@ -58,16 +58,32 @@ MixedEquilibriumKineticReactions< REAL_TYPE, ARRAY_1D_PRIMARY & aggregateSpeciesRates, ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ) { - // 1. Compute new aggregate species from primary species - massActions::calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, - INT_TYPE, - INDEX_TYPE >( params.equilibriumReactionsParameters(), - logPrimarySpeciesConcentrations, - logSecondarySpeciesConcentrations, - aggregatePrimarySpeciesConcentrations, - mobileAggregatePrimarySpeciesConcentrations, - dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, - dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + if constexpr( PARAMS_DATA::numEquilibriumReactions() > 0 ) + { + // 1. Compute new aggregate species from primary species + massActions::calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params.equilibriumReactionsParameters(), + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + mobileAggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, + dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + } + else + { + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); + aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; + mobileAggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; + dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; + dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; + } + } if constexpr( PARAMS_DATA::numKineticReactions() > 0 ) { diff --git a/src/reactions/reactionsSystems/Parameters.hpp b/src/reactions/reactionsSystems/Parameters.hpp index 0b6c1bf..54cddaa 100644 --- a/src/reactions/reactionsSystems/Parameters.hpp +++ b/src/reactions/reactionsSystems/Parameters.hpp @@ -90,11 +90,13 @@ struct KineticReactionsParameters constexpr KineticReactionsParameters( CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > const & stoichiometricMatrix, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse, - CArrayWrapper< RealType, NUM_REACTIONS > const & equilibriumConstant ): + CArrayWrapper< RealType, NUM_REACTIONS > const & equilibriumConstant, + IntType const reactionRatesUpdateOption ): m_stoichiometricMatrix( stoichiometricMatrix ), m_rateConstantForward( rateConstantForward ), m_rateConstantReverse( rateConstantReverse ), - m_equilibiriumConstant( equilibriumConstant ) // Initialize to empty array + m_equilibiriumConstant( equilibriumConstant ), // Initialize to empty array + m_reactionRatesUpdateOption( reactionRatesUpdateOption ) {} @@ -103,11 +105,14 @@ struct KineticReactionsParameters HPCREACT_HOST_DEVICE RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; } HPCREACT_HOST_DEVICE RealType equilibriumConstant( IndexType const r ) const { return m_rateConstantForward[r] / m_rateConstantReverse[r]; } + HPCREACT_HOST_DEVICE IntType reactionRatesUpdateOption() const { return m_reactionRatesUpdateOption; } CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; CArrayWrapper< RealType, NUM_REACTIONS > m_equilibiriumConstant; + + IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. }; @@ -130,12 +135,14 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > const & equilibriumConstant, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse, - CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag ): + CArrayWrapper< IntType, NUM_REACTIONS > mobileSecondarySpeciesFlag, + IntType const reactionRatesUpdateOption = 1 ): m_stoichiometricMatrix( stoichiometricMatrix ), m_equilibriumConstant( equilibriumConstant ), m_rateConstantForward( rateConstantForward ), m_rateConstantReverse( rateConstantReverse ), - m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ) + m_mobileSecondarySpeciesFlag( mobileSecondarySpeciesFlag ), + m_reactionRatesUpdateOption( reactionRatesUpdateOption ) {} HPCREACT_HOST_DEVICE static constexpr IndexType numReactions() { return NUM_REACTIONS; } @@ -193,7 +200,7 @@ struct MixedReactionsParameters equilibriumConstant( i ) = m_equilibriumConstant( numEquilibriumReactions() + i ); } - return { kineticMatrix, rateConstantForward, rateConstantReverse, equilibriumConstant }; + return { kineticMatrix, rateConstantForward, rateConstantReverse, equilibriumConstant, m_reactionRatesUpdateOption }; } HPCREACT_HOST_DEVICE @@ -242,6 +249,8 @@ struct MixedReactionsParameters CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward; CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; CArrayWrapper< IntType, NUM_REACTIONS > m_mobileSecondarySpeciesFlag; + + IntType m_reactionRatesUpdateOption; // 0: forward and reverse rate. 1: quotient form. }; diff --git a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp index 3e0b3b7..3fe18ed 100644 --- a/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp +++ b/src/reactions/unitTestUtilities/kineticReactionsTestUtilities.hpp @@ -49,6 +49,9 @@ struct ComputeReactionRatesTestData /// The reaction rates derivatives CArrayWrapper< double, numReactions, numSpecies > reactionRatesDerivatives; + + /// The surface area + double surfaceArea[numReactions]; }; template< typename REAL_TYPE, @@ -56,6 +59,7 @@ template< typename REAL_TYPE, 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()] ) { @@ -92,12 +96,18 @@ void computeReactionRatesTest( PARAMS_DATA const & params, } } + for( int r = 0; r < numReactions; ++r ) + { + data.surfaceArea[r] = surfaceArea[r]; + } + pmpl::genericKernelWrapper( 1, &data, [params, temperature] HPCREACT_DEVICE ( auto * const dataCopy ) { KineticReactionsType::computeReactionRates( temperature, params, dataCopy->speciesConcentration, + dataCopy->surfaceArea, dataCopy->reactionRates, dataCopy->reactionRatesDerivatives ); });