diff --git a/.gitignore b/.gitignore index 7583b4b..980f50f 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ build* -cmake/blt* \ No newline at end of file +install* +cmake/blt* +.vscode* \ No newline at end of file diff --git a/hostconfigs/apple/macOS_base.cmake b/hostconfigs/apple/macOS_base.cmake index 9a056a4..358672c 100644 --- a/hostconfigs/apple/macOS_base.cmake +++ b/hostconfigs/apple/macOS_base.cmake @@ -1,5 +1,5 @@ message( "this hostconfig assumes you are using homebrew") -message( "brew install bison cmake gcc git-lfs open-mpi openblas python3") +message( "brew install bison cmake gcc git-lfs open-mpi openblas python3 llvm cppcheck lcov") message( "CMAKE_SYSTEM_PROCESSOR = ${CMAKE_SYSTEM_PROCESSOR}" ) message("CONFIG_NAME = ${CONFIG_NAME}") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 69f694e..01c7ed5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,13 +7,14 @@ set( hpcReact_headers reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp reactions/bulkGeneric/KineticReactions.hpp reactions/bulkGeneric/KineticReactions_impl.hpp + reactions/bulkGeneric/MixedEquilibriumKineticReactions.hpp + reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp reactions/bulkGeneric/Parameters.hpp reactions/bulkGeneric/ParametersPredefined.hpp reactions/bulkGeneric/SpeciesUtilities.hpp ) -set( hpcReact_sources - ) +set( hpcReact_sources) find_package(LAPACK REQUIRED) find_package(BLAS REQUIRED) @@ -35,11 +36,13 @@ blt_add_library( NAME hpcReact target_include_directories( hpcReact INTERFACE - $ - $ + $ + $ $ ) blt_print_target_properties( TARGET hpcReact ) +message(STATUS "HPCReact/src CMAKE_CURRENT_SOURCE_DIR: ${CMAKE_CURRENT_SOURCE_DIR}") + # install( FILES ${hpcReact_headers} # DESTINATION include ) diff --git a/src/common/CArrayWrapper.hpp b/src/common/CArrayWrapper.hpp index 397c58b..c5a5c7c 100644 --- a/src/common/CArrayWrapper.hpp +++ b/src/common/CArrayWrapper.hpp @@ -1,5 +1,6 @@ #pragma once #include "macros.hpp" +#include /** @@ -26,37 +27,72 @@ struct CArrayWrapper; template< typename T, int DIM0 > struct CArrayWrapper< T, DIM0 > { + // default constructor + constexpr CArrayWrapper() = default; + + /** + * @brief Construct a CArrayWrapper from an initializer list. + * + * Allows brace-initialization with a list of values: + * @code + * CArrayWrapper< double, 3 > arr = {1.0, 2.0, 3.0}; + * @endcode + * + * @param init An initializer list with exactly DIM0 elements. + * + * @note No runtime bounds checking is performed on the initializer size. + */ + constexpr CArrayWrapper( std::initializer_list< T > init ) + { + // static_assert(init.size() == DIM0, "Size mismatch"); // needs c++20 + int i = 0; + for( auto const & val : init ) + { + data[i++] = val; + } + } + /** + * @brief Copy constructor. + * @param src The source CArrayWrapper to copy from. + */ + constexpr CArrayWrapper( CArrayWrapper const & src ) + { + for( std::size_t i = 0; i < DIM0; i++ ) + { + data[i] = src.data[i]; + } + } /** * @brief Read/write access to an element by index. * @param dim The index (must be in range [0, DIM0)). * @return Reference to the element at the specified index. */ - HPCREACT_HOST_DEVICE inline T & operator()( int const dim ) { return data[dim]; } + constexpr HPCREACT_HOST_DEVICE inline T & operator()( int const dim ) { return data[dim]; } /** * @brief Read-only access to an element by index (const overload). * @param dim The index (must be in range [0, DIM0)). * @return Const reference to the element at the specified index. */ - HPCREACT_HOST_DEVICE inline T const & operator()( int const dim ) const { return data[dim]; } + constexpr HPCREACT_HOST_DEVICE inline T const & operator()( int const dim ) const { return data[dim]; } /** * @brief Subscript operator for read/write access. * @param dim The index (must be in range [0, DIM0)). * @return Reference to the element at the specified index. */ - HPCREACT_HOST_DEVICE inline T & operator[]( int const dim ) { return data[dim]; } + constexpr HPCREACT_HOST_DEVICE inline T & operator[]( int const dim ) { return data[dim]; } /** * @brief Subscript operator for read-only access (const overload). * @param dim The index (must be in range [0, DIM0)). * @return Const reference to the element at the specified index. */ - HPCREACT_HOST_DEVICE inline T const & operator[]( int const dim ) const { return data[dim]; } + constexpr HPCREACT_HOST_DEVICE inline T const & operator[]( int const dim ) const { return data[dim]; } /// The underlying 1D C-style array. - T data[DIM0]; + T data[DIM0]{}; }; /** @@ -72,13 +108,60 @@ struct CArrayWrapper< T, DIM0 > template< typename T, int DIM0, int DIM1 > struct CArrayWrapper< T, DIM0, DIM1 > { + // default constructor + constexpr CArrayWrapper() = default; + + /** + * @brief Copy constructor. + * @param src The source CArrayWrapper to copy from. + */ + constexpr CArrayWrapper( CArrayWrapper const & src ) + { + for( std::size_t i = 0; i < DIM0; i++ ) + { + for( std::size_t j = 0; j < DIM1; j++ ) + data[i][j] = src.data[i][j]; + } + } + + /** + * @brief Construct a 2D CArrayWrapper from nested initializer lists. + * + * Allows brace-initialization with a matrix-like structure: + * @code + * CArrayWrapper mat = { + * {1.0, 2.0, 3.0}, + * {4.0, 5.0, 6.0} + * }; + * @endcode + * + * @param init A nested initializer list with exactly D0 rows and D1 elements per row. + * + * @note No runtime bounds checking is performed on the initializer dimensions. + */ + constexpr CArrayWrapper( std::initializer_list< std::initializer_list< T > > init ) + { + // static_assert(init.size() == DIM0, "Size mismatch"); // needs c++20 + int i = 0; + for( auto const & row : init ) + { + // static_assert(row.size() == DIM1, "Size mismatch"); // needs c++20 + int j = 0; + for( auto const & val : row ) + { + data[i][j++] = val; + } + ++i; + } + } + /** * @brief Read/write access to an element by 2D indices. * @param dim0 Index in the first dimension (range [0, DIM0)). * @param dim1 Index in the second dimension (range [0, DIM1)). * @return Reference to the element at the specified 2D location. */ - HPCREACT_HOST_DEVICE inline T & operator()( int const dim0, int const dim1 ) + constexpr HPCREACT_HOST_DEVICE inline T & operator()( int const dim0, int const dim1 ) { return data[dim0][dim1]; } @@ -89,7 +172,7 @@ struct CArrayWrapper< T, DIM0, DIM1 > * @param dim1 Index in the second dimension (range [0, DIM1)). * @return Const reference to the element at the specified 2D location. */ - HPCREACT_HOST_DEVICE inline T const & operator()( int const dim0, int const dim1 ) const + constexpr HPCREACT_HOST_DEVICE inline T const & operator()( int const dim0, int const dim1 ) const { return data[dim0][dim1]; } @@ -101,7 +184,7 @@ struct CArrayWrapper< T, DIM0, DIM1 > * * This allows usage like `obj[dim0][dim1]`. */ - HPCREACT_HOST_DEVICE inline T ( & operator[]( int const dim0 ))[DIM1] + constexpr HPCREACT_HOST_DEVICE inline T ( & operator[]( int const dim0 ))[DIM1] { return data[dim0]; } @@ -111,13 +194,13 @@ struct CArrayWrapper< T, DIM0, DIM1 > * @param dim0 The row index (range [0, DIM0)). * @return Const reference to an array of type T[DIM1]. */ - HPCREACT_HOST_DEVICE inline T const (&operator[]( int const dim0 ) const)[DIM1] + constexpr HPCREACT_HOST_DEVICE inline T const (&operator[]( int const dim0 ) const)[DIM1] { return data[dim0]; } /// The underlying 2D C-style array of size DIM0 x DIM1. - T data[DIM0][DIM1]; + T data[DIM0][DIM1]{}; }; /** @@ -134,6 +217,54 @@ struct CArrayWrapper< T, DIM0, DIM1 > template< typename T, int DIM0, int DIM1, int DIM2 > struct CArrayWrapper< T, DIM0, DIM1, DIM2 > { + // default constructor + constexpr CArrayWrapper() = default; + + /** + * @brief Construct a 3D CArrayWrapper from nested initializer lists. + * + * Enables tensor-like initialization using triple-nested braces: + * @code + * CArrayWrapper cube = { + * { + * {1.0, 2.0}, + * {3.0, 4.0} + * }, + * { + * {5.0, 6.0}, + * {7.0, 8.0} + * } + * }; + * @endcode + * + * @param init A three-level nested initializer list with D0 planes, D1 rows per plane, + * and D2 elements per row. + * + * @note This constructor does not perform size validation. Incorrect initializer sizes + * may lead to undefined behavior. + */ + constexpr CArrayWrapper( std::initializer_list< std::initializer_list< std::initializer_list< T > > > init ) + { + // static_assert(init.size() == DIM0, "Size mismatch"); // needs c++20 + int i = 0; + for( auto const & plane : init ) + { + // static_assert(plane.size() == DIM1, "Size mismatch"); // needs c++20 + int j = 0; + for( auto const & row : plane ) + { + // static_assert(row.size() == DIM2, "Size mismatch"); // needs c++20 + int k = 0; + for( auto const & val : row ) + { + data[i][j][k++] = val; + } + ++j; + } + ++i; + } + } + /** * @brief Read/write access to an element by 3D indices. * @param dim0 Index in the first dimension (range [0, DIM0)). @@ -144,7 +275,7 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 > * @note Currently, this function incorrectly indexes data[dim0][dim1], missing dim2. * It should be `data[dim0][dim1][dim2]`. Please correct if intended. */ - HPCREACT_HOST_DEVICE inline T & operator()( int const dim0, int const dim1, int const dim2 ) + constexpr HPCREACT_HOST_DEVICE inline T & operator()( int const dim0, int const dim1, int const dim2 ) { // NOTE: This looks like a bug in your original code. Should be data[dim0][dim1][dim2]. return data[dim0][dim1][dim2]; @@ -157,7 +288,7 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 > * @param dim2 Index in the third dimension (range [0, DIM2)). * @return Const reference to the element at the specified 3D location. */ - HPCREACT_HOST_DEVICE inline T const & operator()( int const dim0, int const dim1, int const dim2 ) const + constexpr HPCREACT_HOST_DEVICE inline T const & operator()( int const dim0, int const dim1, int const dim2 ) const { // NOTE: Same potential bug as above. Should be data[dim0][dim1][dim2]. return data[dim0][dim1][dim2]; @@ -170,7 +301,7 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 > * * This allows usage like `obj[dim0][dim1][dim2]`. */ - HPCREACT_HOST_DEVICE inline T ( & operator[]( int const dim0 ))[DIM1][DIM2] + constexpr HPCREACT_HOST_DEVICE inline T ( & operator[]( int const dim0 ))[DIM1][DIM2] { return data[dim0]; } @@ -180,11 +311,11 @@ struct CArrayWrapper< T, DIM0, DIM1, DIM2 > * @param dim0 The slice index (range [0, DIM0)). * @return Const reference to an array of type T[DIM1][DIM2]. */ - HPCREACT_HOST_DEVICE inline T const (&operator[]( int const dim0 ) const)[DIM1][DIM2] + constexpr HPCREACT_HOST_DEVICE inline T const (&operator[]( int const dim0 ) const)[DIM1][DIM2] { return data[dim0]; } /// The underlying 3D C-style array of size DIM0 x DIM1 x DIM2. - T data[DIM0][DIM1][DIM2]; + T data[DIM0][DIM1][DIM2]{}; }; diff --git a/src/common/nonlinearSolvers.hpp b/src/common/nonlinearSolvers.hpp new file mode 100644 index 0000000..dd2dd19 --- /dev/null +++ b/src/common/nonlinearSolvers.hpp @@ -0,0 +1,191 @@ +#pragma once + +#include "macros.hpp" +#include "DirectSystemSolve.hpp" +#include + +namespace hpcReact +{ + +namespace nonlinearSolvers +{ + +namespace internal +{ + +/** + * Computes the norm of a vector. + * This function calculates the Euclidean norm (L2 norm) of a vector. + * @tparam N The size of the vector. + * @param r The input vector. + * @return The Euclidean norm of the vector. + */ +template< int N > +HPCREACT_HOST_DEVICE +double norm( double const (&r)[N] ) +{ + double sum = 0.0; + for( int i = 0; i < N; ++i ) + sum += r[i] * r[i]; + return ::sqrt( sum ); +} + +/** + * Adds two vectors element-wise. + * This function adds the elements of the second vector to the first vector. + * @tparam N The size of the vectors. + * @param x The first vector, which will be modified. + * @param dx The second vector, which will be added to the first vector. + */ +template< int N > +HPCREACT_HOST_DEVICE +void add( double (& x)[N], double const (&dx)[N] ) +{ + for( int i = 0; i < N; ++i ) + x[i] += dx[i]; +} + +/** + * Scales a vector by a constant value. + * This function multiplies each element of the vector by a given value. + * @tparam N The size of the vector. + * @param x The vector to be scaled. + * @param value The scaling factor. + */ +template< int N > +HPCREACT_HOST_DEVICE +void scale( double (& x)[N], double const value ) +{ + for( int i = 0; i < N; ++i ) + x[i] *= value; +} + +} + +namespace utils +{ +// LCOV_EXCL_START + +/** + * Prints the Jacobian matrix, residual vector, and delta update vector. + * This function is useful for debugging and understanding the + * state of the Newton-Raphson method at each iteration. + * @tparam N The size of the vectors and matrix. + * @param J The Jacobian matrix. + * @param r The residual vector. + * @param dx The delta update vector. + */ +template< int N > +HPCREACT_HOST_DEVICE +void print( double const (&J)[N][N], double const (&r)[N], double const (&dx)[N] ) +{ + printf( "=======================\nJacobian matrix:\n=======================\n" ); + printf( " RowID ColID Value\n" ); for( int i = 0; i < N; ++i ) + { + for( int j = 0; j < N; ++j ) + { + printf( "%10d%16d%27.16e\n", i, j, J[i][j] ); + } + } + + printf( "\n=======================\nSystem right-hand side:\n=======================\n" ); + printf( " RowID Value\n" ); + + for( int i = 0; i < N; ++i ) + { + printf( "%10d%27.16e\n", i, r[i] ); + } + + printf( "\n=======================\nDelta update vector:\n=======================\n" ); + printf( " RowID Value\n" ); + + for( int i = 0; i < N; ++i ) + { + printf( "%10d%27.16e\n", i, dx[i] ); + } +} +// LCOV_EXCL_STOP + +} + +template< int N, + typename REAL_TYPE, + typename ResidualFunc, + typename JacobianFunc > +HPCREACT_HOST_DEVICE +void newtonRaphson( REAL_TYPE (& x)[N], + ResidualFunc computeResidual, + JacobianFunc computeJacobian, + int maxIters = 25, + double tol = 1e-10 ) +{ + REAL_TYPE residual[N]; + REAL_TYPE dx[N]; + REAL_TYPE jacobian[N][N]; + + for( int iter = 0; iter < maxIters; ++iter ) + { + computeResidual( x, residual ); + + if( internal::norm< N >( residual ) < tol ) + return; + + computeJacobian( x, jacobian ); + + solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx ); + + internal::add< N >( x, dx ); + } +} + +template< int N, + typename REAL_TYPE, + typename FUNCTION_TYPE > +HPCREACT_HOST_DEVICE +bool newtonRaphson( REAL_TYPE (& x)[N], + FUNCTION_TYPE computeResidualAndJacobian, + int maxIters = 12, + double tol = 1e-10, + bool const do_print = false ) +{ + REAL_TYPE residual[N]{}; + REAL_TYPE dx[N]{}; + REAL_TYPE jacobian[N][N]{}; + bool isConverged = false; + + for( int iter = 0; iter < maxIters; ++iter ) + { + computeResidualAndJacobian( x, residual, jacobian ); + + double const norm = internal::norm< N >( residual ); + + printf( "--Iter %d: Residual norm = %.12e\n", iter, norm ); + + if( norm < tol ) + { + printf( "--Converged.\n" ); + isConverged = true; + break; + } + internal::scale< N >( residual, -1.0 ); + + if( do_print ) + { + utils::print( jacobian, residual, dx ); // LCOV_EXCL_LINE + } + + solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx ); + internal::add< N >( x, dx ); + + } + + if( !isConverged ) + { + printf( "--Newton solver error: Max iterations reached without convergence.\n" ); // LCOV_EXCL_LINE + } + + return isConverged; +} + +} +} diff --git a/src/reactions/bulkGeneric/EquilibriumReactions.hpp b/src/reactions/bulkGeneric/EquilibriumReactions.hpp index f6568ec..fca5a2c 100644 --- a/src/reactions/bulkGeneric/EquilibriumReactions.hpp +++ b/src/reactions/bulkGeneric/EquilibriumReactions.hpp @@ -132,12 +132,13 @@ class EquilibriumReactions template< typename PARAMS_DATA, typename ARRAY_1D, typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, typename ARRAY_2D > static HPCREACT_HOST_DEVICE void computeResidualAndJacobianAggregatePrimaryConcentrations( RealType const & temperature, PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & targetAggregatePrimaryConcentrations, - ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration, + ARRAY_1D_TO_CONST2 const & logPrimarySpeciesConcentration, ARRAY_1D & residual, ARRAY_2D & jacobian ); }; diff --git a/src/reactions/bulkGeneric/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp b/src/reactions/bulkGeneric/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp index 3454a12..1b8e534 100644 --- a/src/reactions/bulkGeneric/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp +++ b/src/reactions/bulkGeneric/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp @@ -16,6 +16,7 @@ template< typename REAL_TYPE, template< typename PARAMS_DATA, typename ARRAY_1D, typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, typename ARRAY_2D > HPCREACT_HOST_DEVICE inline @@ -25,14 +26,12 @@ EquilibriumReactions< REAL_TYPE, INDEX_TYPE >::computeResidualAndJacobianAggregatePrimaryConcentrations( RealType const & temperature, PARAMS_DATA const & params, ARRAY_1D_TO_CONST const & targetAggregatePrimaryConcentrations, - ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration, + ARRAY_1D_TO_CONST2 const & logPrimarySpeciesConcentration, ARRAY_1D & residual, ARRAY_2D & jacobian ) { HPCREACT_UNUSED_VAR( temperature ); - constexpr int numSpecies = PARAMS_DATA::numSpecies; - constexpr int numSecondarySpecies = PARAMS_DATA::numReactions; - constexpr int numPrimarySpecies = numSpecies - numSecondarySpecies; + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); RealType aggregatePrimaryConcentrations[numPrimarySpecies] = {0.0}; ARRAY_2D dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations = {{{0.0}}}; @@ -68,9 +67,7 @@ EquilibriumReactions< REAL_TYPE, ARRAY_1D & logPrimarySpeciesConcentration ) { HPCREACT_UNUSED_VAR( temperature ); - constexpr int numSpecies = PARAMS_DATA::numSpecies; - constexpr int numReactions = PARAMS_DATA::numReactions; - constexpr int numPrimarySpecies = numSpecies - numReactions; + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); double residual[numPrimarySpecies] = { 0.0 }; // double aggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 }; diff --git a/src/reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp b/src/reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp index 923b0c2..1078138 100644 --- a/src/reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp +++ b/src/reactions/bulkGeneric/EquilibriumReactionsReactionExtents_impl.hpp @@ -32,8 +32,8 @@ EquilibriumReactions< REAL_TYPE, { HPCREACT_UNUSED_VAR( temperature ); - constexpr int numSpecies = PARAMS_DATA::numSpecies; - constexpr int numReactions = PARAMS_DATA::numReactions; + constexpr int numSpecies = PARAMS_DATA::numSpecies(); + constexpr int numReactions = PARAMS_DATA::numReactions(); // initialize the species concentration RealType speciesConcentration[numSpecies]; @@ -115,8 +115,8 @@ EquilibriumReactions< REAL_TYPE, ARRAY_1D & speciesConcentration ) { HPCREACT_UNUSED_VAR( temperature ); - constexpr int numSpecies = PARAMS_DATA::numSpecies; - constexpr int numReactions = PARAMS_DATA::numReactions; + constexpr int numSpecies = PARAMS_DATA::numSpecies(); + constexpr int numReactions = PARAMS_DATA::numReactions(); double residual[numReactions] = { 0.0 }; double xi[numReactions] = { 0.0 }; double dxi[numReactions] = { 0.0 }; diff --git a/src/reactions/bulkGeneric/KineticReactions.hpp b/src/reactions/bulkGeneric/KineticReactions.hpp index 48b79f2..a094608 100644 --- a/src/reactions/bulkGeneric/KineticReactions.hpp +++ b/src/reactions/bulkGeneric/KineticReactions.hpp @@ -79,7 +79,7 @@ class KineticReactions ARRAY_1D_TO_CONST const & speciesConcentration, ARRAY_1D & reactionRates ) { - REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions][PARAMS_DATA::numSpecies] = { {0.0} }; + REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions()][PARAMS_DATA::numSpecies()] = { {0.0} }; computeReactionRates_impl< PARAMS_DATA, false >( temperature, params, speciesConcentration, diff --git a/src/reactions/bulkGeneric/KineticReactions_impl.hpp b/src/reactions/bulkGeneric/KineticReactions_impl.hpp index 225f239..322261e 100644 --- a/src/reactions/bulkGeneric/KineticReactions_impl.hpp +++ b/src/reactions/bulkGeneric/KineticReactions_impl.hpp @@ -47,7 +47,7 @@ KineticReactions< REAL_TYPE, } // loop over each reaction - for( IntType r=0; r 0.0 ) { - for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j ) + for( IntType j = 0; j < PARAMS_DATA::numSpecies(); ++j ) { if( i==j ) { @@ -173,7 +174,7 @@ KineticReactions< REAL_TYPE, if constexpr( CALCULATE_DERIVATIVES ) { - for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i ) + for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) { reactionRatesDerivatives( r, i ) = forwardRateConstant * dProductConcForward_dC[i] - reverseRateConstant * dProductConcReverse_dC[i]; } @@ -205,8 +206,8 @@ KineticReactions< REAL_TYPE, ARRAY_1D & speciesRates, ARRAY_2D & speciesRatesDerivatives ) { - RealType reactionRates[PARAMS_DATA::numReactions] = { 0.0 }; - CArrayWrapper< double, PARAMS_DATA::numReactions, PARAMS_DATA::numSpecies > reactionRatesDerivatives; + RealType reactionRates[PARAMS_DATA::numReactions()] = { 0.0 }; + CArrayWrapper< double, PARAMS_DATA::numReactions(), PARAMS_DATA::numSpecies() > reactionRatesDerivatives; if constexpr( !CALCULATE_DERIVATIVES ) { @@ -215,23 +216,23 @@ KineticReactions< REAL_TYPE, computeReactionRates< PARAMS_DATA >( temperature, params, speciesConcentration, reactionRates, reactionRatesDerivatives ); - for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i ) + for( IntType i = 0; i < PARAMS_DATA::numSpecies(); ++i ) { speciesRates[i] = 0.0; if constexpr( CALCULATE_DERIVATIVES ) { - for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j ) + for( IntType j = 0; j < PARAMS_DATA::numSpecies(); ++j ) { speciesRatesDerivatives( i, j ) = 0.0; } } - for( IntType r=0; r +class MixedEquilibriumKineticReactions +{ +public: + + /// Type alias for the real type used in the class. + using RealType = REAL_TYPE; + + /// Type alias for the integer type used in the class. + using IntType = INT_TYPE; + + /// Type alias for the index type used in the class. + using IndexType = INDEX_TYPE; + + /// Type alias for the Kinetic reactions type used in the class. + using kineticReactions = KineticReactions< REAL_TYPE, INT_TYPE, INDEX_TYPE, LOGE_CONCENTRATION >; + + /** + * @brief Update a mixed chemical system by computing secondary species concentrations, + * aggregate primary species concentrations, and reaction rates. + * + * @tparam PARAMS_DATA Struct providing all parameter access (stoichiometry, rate constants, etc.) + * @tparam ARRAY_1D_TO_CONST Read-only 1D array type for primary log-concentrations + * @tparam ARRAY_1D_PRIMARY Mutable 1D array type for primary species outputs + * @tparam ARRAY_1D_SECONDARY Mutable 1D array type for secondary log-concentrations + * @tparam ARRAY_1D_KINETIC Mutable 1D array type for reaction rates + * @tparam ARRAY_2D_PRIMARY Mutable 2D array type for primary derivatives + * @tparam ARRAY_2D_KINETIC Mutable 2D array type for reaction rate derivatives + * + * @param temperature Temperature of the system (in Kelvin) + * @param params Parameter object for stoichiometry, rates, etc. + * @param logPrimarySpeciesConcentrations Log of primary species concentrations + * @param logSecondarySpeciesConcentrations Output log concentrations for secondary species + * @param aggregatePrimarySpeciesConcentrations Output aggregate concentrations (per primary) + * @param dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations Derivatives of aggregate concentrations w.r.t. log + * primary + * @param reactionRates Output vector of kinetic reaction rates + * @param dReactionRates_dLogPrimarySpeciesConcentrations Derivatives of reaction rates w.r.t. log primary species + * @param aggregateSpeciesRates Output net source/sink for each primary species + * @param dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations Derivatives of aggregate source terms + */ + template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_PRIMARY, + typename ARRAY_1D_SECONDARY, + typename ARRAY_1D_KINETIC, + typename ARRAY_2D_PRIMARY, + typename ARRAY_2D_KINETIC > + static HPCREACT_HOST_DEVICE inline void + updateMixedSystem( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, + ARRAY_1D_KINETIC & reactionRates, + ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregateSpeciesRates, + ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ) + { + updateMixedSystem_impl( temperature, + params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, + reactionRates, + dReactionRates_dLogPrimarySpeciesConcentrations, + aggregateSpeciesRates, + dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); + } + + /** + * @brief Compute reaction rates and their derivatives. + * + * @tparam PARAMS_DATA Struct providing reaction parameters + * @tparam ARRAY_1D_TO_CONST Read-only array of primary species (log-space) + * @tparam ARRAY_1D_TO_CONST2 Read-only array of secondary species (log-space) + * @tparam ARRAY_1D Output array type for reaction rates + * @tparam ARRAY_2D Output array type for reaction rate derivatives + * + * @param temperature Temperature in Kelvin + * @param params Parameter data for the reaction system + * @param logPrimarySpeciesConcentrations Log concentrations of primary species + * @param logSecondarySpeciesConcentrations Log concentrations of secondary species + * @param reactionRates Output reaction rates for each kinetic reaction + * @param dReactionRates_dLogPrimarySpeciesConcentrations Derivatives of reaction rates w.r.t. log primary species + */ + template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_1D, + typename ARRAY_2D > + static HPCREACT_HOST_DEVICE inline void + computeReactionRates( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST2 const & logSecondarySpeciesConcentrations, + ARRAY_1D & reactionRates, + ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations ) + + { + computeReactionRates_impl( temperature, + params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + reactionRates, + dReactionRates_dLogPrimarySpeciesConcentrations ); + } + + /** + * @brief Compute net reaction rate for each primary species by aggregating contributions from all reactions. + * + * @tparam PARAMS_DATA Struct with stoichiometry and mappings + * @tparam ARRAY_1D_TO_CONST Array type for primary species concentrations + * @tparam ARRAY_1D_TO_CONST2 Array type for reaction rates + * @tparam ARRAY_2D_TO_CONST Array type for reaction rate derivatives + * @tparam ARRAY_1D Output type for net species rates + * @tparam ARRAY_2D Output type for rate derivatives + * + * @param params Reaction parameters + * @param speciesConcentration Current concentrations of primary species + * @param reactionRates Computed reaction rates + * @param reactionRatesDerivatives Derivatives of reaction rates w.r.t. log concentrations + * @param aggregatesRates Output: net rate for each primary species + * @param aggregatesRatesDerivatives Output: derivative of net rates w.r.t. log concentrations + */ + template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_2D_TO_CONST, + typename ARRAY_1D, + typename ARRAY_2D > + static HPCREACT_HOST_DEVICE inline void + computeAggregateSpeciesRates( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D_TO_CONST2 const & reactionRates, + ARRAY_2D_TO_CONST const & reactionRatesDerivatives, + ARRAY_1D & aggregatesRates, + ARRAY_2D & aggregatesRatesDerivatives ) + { + computeAggregateSpeciesRates_impl< PARAMS_DATA, + ARRAY_1D_TO_CONST, + ARRAY_1D_TO_CONST2, + ARRAY_2D_TO_CONST, + ARRAY_1D, + ARRAY_2D, + true >( params, + speciesConcentration, + reactionRates, + reactionRatesDerivatives, + aggregatesRates, + aggregatesRatesDerivatives ); + } + +private: + /** + * @brief Internal implementation of updateMixedSystem with template-dispatched logic. + * + * @details Called by the public `updateMixedSystem` function. Handles the complete chain: + * secondary speciation, aggregation, reaction rate evaluation, and net source terms. + * @tparam PARAMS_DATA Struct providing all parameter access (stoichiometry, rate constants, etc.) + * @tparam ARRAY_1D_TO_CONST Read-only 1D array type for primary log-concentrations + * @tparam ARRAY_1D_PRIMARY Mutable 1D array type for primary species outputs + * @tparam ARRAY_1D_SECONDARY Mutable 1D array type for secondary log-concentrations + * @tparam ARRAY_1D_KINETIC Mutable 1D array type for reaction rates + * @tparam ARRAY_2D_PRIMARY Mutable 2D array type for primary derivatives + * @tparam ARRAY_2D_KINETIC Mutable 2D array type for reaction rate derivatives + * + * @param temperature Temperature of the system (in Kelvin) + * @param params Parameter object for stoichiometry, rates, etc. + * @param logPrimarySpeciesConcentrations Log of primary species concentrations + * @param logSecondarySpeciesConcentrations Output log concentrations for secondary species + * @param aggregatePrimarySpeciesConcentrations Output aggregate concentrations (per primary) + * @param dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations Derivatives of aggregate concentrations w.r.t. log + * primary + * @param reactionRates Output vector of kinetic reaction rates + * @param dReactionRates_dLogPrimarySpeciesConcentrations Derivatives of reaction rates w.r.t. log primary species + * @param aggregateSpeciesRates Output net source/sink for each primary species + * @param dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations Derivatives of aggregate source terms + */ + template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_PRIMARY, + typename ARRAY_1D_SECONDARY, + typename ARRAY_1D_KINETIC, + typename ARRAY_2D_PRIMARY, + typename ARRAY_2D_KINETIC > + static HPCREACT_HOST_DEVICE void + updateMixedSystem_impl( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, + ARRAY_1D_KINETIC & reactionRates, + ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregateSpeciesRates, + ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); + /** + * @brief Internal implementation of computeReactionRates. + * + * @details Handles kinetic rate law evaluation for forward and reverse reactions. + * @param temperature Temperature in Kelvin + * @param params Parameter data for the reaction system + * @param logPrimarySpeciesConcentrations Log concentrations of primary species + * @param logSecondarySpeciesConcentrations Log concentrations of secondary species + * @param reactionRates Output reaction rates for each kinetic reaction + * @param dReactionRates_dLogPrimarySpeciesConcentrations Derivatives of reaction rates w.r.t. log primary species + */ + template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_1D, + typename ARRAY_2D > + static HPCREACT_HOST_DEVICE void + computeReactionRates_impl( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST2 const & logSecondarySpeciesConcentrations, + ARRAY_1D & reactionRates, + ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations ); + + + /** + * @brief Internal implementation of computeAggregateSpeciesRates. + * + * @tparam CALCULATE_DERIVATIVES Whether to compute Jacobian derivatives + * @tparam ARRAY_1D_TO_CONST Array type for primary species concentrations + * @tparam ARRAY_1D_TO_CONST2 Array type for reaction rates + * @tparam ARRAY_2D_TO_CONST Array type for reaction rate derivatives + * @tparam ARRAY_1D Output type for net species rates + * @tparam ARRAY_2D Output type for rate derivatives + * + * @param params Reaction parameters + * @param speciesConcentration Current concentrations of primary species + * @param reactionRates Computed reaction rates + * @param reactionRatesDerivatives Derivatives of reaction rates w.r.t. log concentrations + * @param aggregatesRates Output: net rate for each primary species + * @param aggregatesRatesDerivatives Output: derivative of net rates w.r.t. log concentrations + */ + template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_2D_TO_CONST, + typename ARRAY_1D, + typename ARRAY_2D, + bool CALCULATE_DERIVATIVES > + static HPCREACT_HOST_DEVICE void + computeAggregateSpeciesRates_impl( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D_TO_CONST2 const & reactionRates, + ARRAY_2D_TO_CONST const & reactionRatesDerivatives, + ARRAY_1D & aggregatesRates, + ARRAY_2D & aggregatesRatesDerivatives ); + + +}; + +} // namespace bulkGeneric +} // namespace hpcReact + +#include "MixedEquilibriumKineticReactions_impl.hpp" +#include "common/macrosCleanup.hpp" diff --git a/src/reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp b/src/reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp new file mode 100644 index 0000000..627527e --- /dev/null +++ b/src/reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp @@ -0,0 +1,214 @@ +#pragma once + +#include "common/constants.hpp" +#include "common/CArrayWrapper.hpp" +#include "SpeciesUtilities.hpp" + + +/** @file MixedEquilibriumKineticReactions_impl.hpp + * @brief Header file for the MixedEquilibriumKineticReactions implementation. + * @author HPC-REACT Team + * @date 2025 + */ + +namespace hpcReact +{ +namespace bulkGeneric +{ + +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + bool LOGE_CONCENTRATION > +template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_PRIMARY, + typename ARRAY_1D_SECONDARY, + typename ARRAY_1D_KINETIC, + typename ARRAY_2D_PRIMARY, + typename ARRAY_2D_KINETIC > +HPCREACT_HOST_DEVICE inline void +MixedEquilibriumKineticReactions< REAL_TYPE, + INT_TYPE, + INDEX_TYPE, + LOGE_CONCENTRATION + >::updateMixedSystem_impl( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations, + ARRAY_1D_KINETIC & reactionRates, + ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregateSpeciesRates, + ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ) +{ + + // 1. Compute new aggregate species from primary species + calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, + INT_TYPE, + INDEX_TYPE >( params.equilibriumReactionsParameters(), + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + aggregatePrimarySpeciesConcentrations, + dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); + + if constexpr( PARAMS_DATA::numKineticReactions() > 0 ) + { + + + // 2. Compute the reaction rates for all kinetic reactions + computeReactionRates( temperature, + params, + logPrimarySpeciesConcentrations, + logSecondarySpeciesConcentrations, + reactionRates, + dReactionRates_dLogPrimarySpeciesConcentrations ); + + + // 3. Compute aggregate species rates + computeAggregateSpeciesRates( params, + logPrimarySpeciesConcentrations, + reactionRates, + dReactionRates_dLogPrimarySpeciesConcentrations, + aggregateSpeciesRates, + dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); + } + else + { + GEOS_UNUSED_VAR( reactionRates, dReactionRates_dLogPrimarySpeciesConcentrations, aggregateSpeciesRates, dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations ); + } + +} + +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + bool LOGE_CONCENTRATION > +template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_1D, + typename ARRAY_2D > +HPCREACT_HOST_DEVICE inline void +MixedEquilibriumKineticReactions< REAL_TYPE, + INT_TYPE, + INDEX_TYPE, + LOGE_CONCENTRATION + >::computeReactionRates_impl( RealType const & temperature, + PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_TO_CONST2 const & logSecondarySpeciesConcentrations, + ARRAY_1D & reactionRates, + ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations ) + +{ + constexpr IntType numSpecies = PARAMS_DATA::numSpecies(); + constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + constexpr IntType numKineticReactions = PARAMS_DATA::numKineticReactions(); + + RealType logSpeciesConcentration[numSpecies] {}; + for( INDEX_TYPE i = 0; i < numSecondarySpecies; ++i ) + { + logSpeciesConcentration[i] = logSecondarySpeciesConcentrations[i]; + } + for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) + { + logSpeciesConcentration[i+numSecondarySpecies] = logPrimarySpeciesConcentrations[i]; + } + + for( INDEX_TYPE i = 0; i < numKineticReactions; ++i ) + { + for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) + { + dReactionRates_dLogPrimarySpeciesConcentrations[i][j] = 0.0; + } + } + + CArrayWrapper< RealType, numKineticReactions, numSpecies > reactionRatesDerivatives; + + kineticReactions::computeReactionRates( temperature, + params.kineticReactionsParameters(), + logSpeciesConcentration, + reactionRates, + reactionRatesDerivatives ); + + // Compute the reaction rates derivatives w.r.t. log primary species concentrations + for( IntType i = 0; i < numKineticReactions; ++i ) + { + for( IntType j = 0; j < numPrimarySpecies; ++j ) + { + dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) = reactionRatesDerivatives( i, j + numSecondarySpecies ); + + for( IntType k = 0; k < numSecondarySpecies; ++k ) + { + RealType const dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations = params.stoichiometricMatrix( k, j + numSecondarySpecies ); + + dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) += + reactionRatesDerivatives( i, k ) * dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations; + } + } + } +} + +template< typename REAL_TYPE, + typename INT_TYPE, + typename INDEX_TYPE, + bool LOGE_CONCENTRATION > +template< typename PARAMS_DATA, + typename ARRAY_1D_TO_CONST, + typename ARRAY_1D_TO_CONST2, + typename ARRAY_2D_TO_CONST, + typename ARRAY_1D, + typename ARRAY_2D, + bool CALCULATE_DERIVATIVES > +HPCREACT_HOST_DEVICE inline void +MixedEquilibriumKineticReactions< REAL_TYPE, + INT_TYPE, + INDEX_TYPE, + LOGE_CONCENTRATION + >::computeAggregateSpeciesRates_impl( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & speciesConcentration, + ARRAY_1D_TO_CONST2 const & reactionRates, + ARRAY_2D_TO_CONST const & reactionRatesDerivatives, + ARRAY_1D & aggregatesRates, + ARRAY_2D & aggregatesRatesDerivatives ) +{ + HPCREACT_UNUSED_VAR( speciesConcentration ); + // constexpr IntType numSpecies = PARAMS_DATA::numSpecies(); + constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + constexpr IntType numKineticReactions = PARAMS_DATA::numKineticReactions(); + constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + + for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i ) + { + aggregatesRates[i] = 0.0; + if constexpr( CALCULATE_DERIVATIVES ) + { + for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j ) + { + aggregatesRatesDerivatives( i, j ) = 0.0; + } + } + for( IntType r=0; r(), - std::make_index_sequence< NUM_REACTIONS *NUM_SPECIES >() ) + EquilibriumReactionsParameters( CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > const & stoichiometricMatrix, + CArrayWrapper< RealType, NUM_REACTIONS > equilibriumConstant ): + m_stoichiometricMatrix( stoichiometricMatrix ), + m_equilibriumConstant( equilibriumConstant ) {} RealType stoichiometricMatrix( IndexType const r, int const i ) const { return m_stoichiometricMatrix[r][i]; } RealType equilibriumConstant( IndexType const r ) const { return m_equilibriumConstant[r]; } - RealType m_stoichiometricMatrix[numReactions][numSpecies]; - RealType m_equilibriumConstant[numReactions]; - -private: - HPCREACT_NO_MISSING_BRACES_OPEN - template< std::size_t ... R, std::size_t ... RxS > - constexpr - EquilibriumReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies], - RealType const (&equilibriumConstant)[numReactions], - std::index_sequence< R... >, - std::index_sequence< RxS... > ): - m_stoichiometricMatrix{ stoichiometricMatrix[RxS/numSpecies][RxS%numSpecies] ... }, - m_equilibriumConstant{ equilibriumConstant[R] ... } - {} - HPCREACT_NO_MISSING_BRACES_CLOSE + CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; + CArrayWrapper< RealType, NUM_REACTIONS > m_equilibriumConstant; }; @@ -72,17 +62,16 @@ struct KineticReactionsParameters using IntType = INT_TYPE; using IndexType = INDEX_TYPE; - static constexpr IndexType numSpecies = NUM_SPECIES; - static constexpr IndexType numReactions = NUM_REACTIONS; - - KineticReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies], - RealType const (&rateConstantForward)[numReactions], - RealType const (&rateConstantReverse)[numReactions] ): - KineticReactionsParameters( stoichiometricMatrix, - rateConstantForward, - rateConstantReverse, - std::make_index_sequence< NUM_REACTIONS >(), - std::make_index_sequence< NUM_REACTIONS *NUM_SPECIES >() ) + static constexpr IndexType numSpecies() { return NUM_SPECIES; } + + static constexpr IndexType numReactions() { return NUM_REACTIONS; } + + constexpr KineticReactionsParameters( CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > const & stoichiometricMatrix, + CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantForward, + CArrayWrapper< RealType, NUM_REACTIONS > const & rateConstantReverse ): + m_stoichiometricMatrix( stoichiometricMatrix ), + m_rateConstantForward( rateConstantForward ), + m_rateConstantReverse( rateConstantReverse ) {} @@ -91,24 +80,9 @@ struct KineticReactionsParameters RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; } - RealType m_stoichiometricMatrix[numReactions][numSpecies]; - RealType m_rateConstantForward[numReactions]; - RealType m_rateConstantReverse[numReactions]; - -private: - HPCREACT_NO_MISSING_BRACES( - template< std::size_t ... R, std::size_t ... RxS > - KineticReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies], - RealType const (&rateConstantForward)[numReactions], - RealType const (&rateConstantReverse)[numReactions], - std::index_sequence< R... >, - std::index_sequence< RxS... > ) : - m_stoichiometricMatrix{ stoichiometricMatrix[RxS/numSpecies][RxS%numSpecies] ... }, - m_rateConstantForward{ rateConstantForward[R] ... }, - m_rateConstantReverse{ rateConstantReverse[R] ... } - {} - ) - + CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; + CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward; + CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; }; @@ -116,33 +90,83 @@ template< typename REAL_TYPE, typename INT_TYPE, typename INDEX_TYPE, int NUM_SPECIES, - int NUM_REACTIONS > + int NUM_REACTIONS, + int NUM_EQ_REACTIONS > struct MixedReactionsParameters { + using RealType = REAL_TYPE; using IntType = INT_TYPE; using IndexType = INDEX_TYPE; - static constexpr IndexType numSpecies = NUM_SPECIES; - static constexpr IndexType numReactions = NUM_REACTIONS; + + constexpr MixedReactionsParameters() = default; + + 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 ): + m_stoichiometricMatrix( stoichiometricMatrix ), + m_equilibriumConstant( equilibriumConstant ), + m_rateConstantForward( rateConstantForward ), + m_rateConstantReverse( rateConstantReverse ) + {} + + static constexpr IndexType numReactions() { return NUM_REACTIONS; } + + static constexpr IndexType numKineticReactions() { return NUM_REACTIONS - NUM_EQ_REACTIONS; } + + static constexpr IndexType numEquilibriumReactions() { return NUM_EQ_REACTIONS; } + + static constexpr IndexType numSpecies() { return NUM_SPECIES; } + + static constexpr IndexType numPrimarySpecies() { return NUM_SPECIES - NUM_EQ_REACTIONS; } + + static constexpr IndexType numSecondarySpecies() { return NUM_EQ_REACTIONS; } constexpr - EquilibriumReactionsParameters< RealType, IntType, IndexType, numSpecies, numReactions > + EquilibriumReactionsParameters< RealType, IntType, IndexType, numSpecies(), numEquilibriumReactions() > equilibriumReactionsParameters() const { - return {m_stoichiometricMatrix, m_equilibriumConstant}; + CArrayWrapper< RealType, numEquilibriumReactions(), numSpecies() > eqMatrix{}; + CArrayWrapper< RealType, numEquilibriumReactions() > eqConstants{}; + + for( IntType i = 0; i < numEquilibriumReactions(); ++i ) + { + for( IntType j = 0; j < numSpecies(); ++j ) + { + eqMatrix( i, j ) = m_stoichiometricMatrix( i, j ); + } + eqConstants( i ) = m_equilibriumConstant( i ); + } + + return { eqMatrix, eqConstants }; } constexpr - KineticReactionsParameters< RealType, IntType, IndexType, numSpecies, numReactions > + KineticReactionsParameters< RealType, IntType, IndexType, numSpecies(), numKineticReactions() > kineticReactionsParameters() const { - return {m_stoichiometricMatrix, m_rateConstantForward, m_rateConstantReverse}; + CArrayWrapper< RealType, numKineticReactions(), numSpecies() > kineticMatrix{}; + CArrayWrapper< RealType, numKineticReactions() > rateConstantForward{}; + CArrayWrapper< RealType, numKineticReactions() > rateConstantReverse{}; + + for( IndexType i = 0; i < numKineticReactions(); ++i ) + { + for( IndexType j = 0; j < numSpecies(); ++j ) + { + kineticMatrix( i, j ) = m_stoichiometricMatrix( numEquilibriumReactions() + i, j ); + } + rateConstantForward( i ) = m_rateConstantForward( numEquilibriumReactions() + i ); + rateConstantReverse( i ) = m_rateConstantReverse( numEquilibriumReactions() + i ); + } + + return { kineticMatrix, rateConstantForward, rateConstantReverse }; } void verifyParameterConsistency() { static constexpr int num_digits = 12; - for( int i = 0; i < numReactions; ++i ) + for( int i = 0; i < numReactions(); ++i ) { RealType & K = m_equilibriumConstant[i]; RealType & kf = m_rateConstantForward[i]; @@ -179,10 +203,10 @@ struct MixedReactionsParameters RealType rateConstantForward( IndexType const r ) const { return m_rateConstantForward[r]; } RealType rateConstantReverse( IndexType const r ) const { return m_rateConstantReverse[r]; } - RealType m_stoichiometricMatrix[numReactions][numSpecies]; - RealType m_equilibriumConstant[numReactions]; - RealType m_rateConstantForward[numReactions]; - RealType m_rateConstantReverse[numReactions]; + CArrayWrapper< RealType, NUM_REACTIONS, NUM_SPECIES > m_stoichiometricMatrix; + CArrayWrapper< RealType, NUM_REACTIONS > m_equilibriumConstant; + CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantForward; + CArrayWrapper< RealType, NUM_REACTIONS > m_rateConstantReverse; }; diff --git a/src/reactions/bulkGeneric/ParametersPredefined.hpp b/src/reactions/bulkGeneric/ParametersPredefined.hpp index be1b0bc..2f6edc3 100644 --- a/src/reactions/bulkGeneric/ParametersPredefined.hpp +++ b/src/reactions/bulkGeneric/ParametersPredefined.hpp @@ -34,14 +34,16 @@ namespace bulkGeneric // um1Constants }; +using simpleKineticTestType = MixedReactionsParameters< double, int, int, 5, 2, 0 >; + constexpr -MixedReactionsParameters< double, int, int, 5, 2 > -simpleTestRateParams = +simpleKineticTestType +simpleKineticTestRateParams = { // stoichiometric matrix { { -2, 1, 1, 0, 0 }, - { 0, 0, -1, -1, 2 } + { 0, 0, -1, -1, 2 } }, // equilibrium constants { 1.0, 1.0 }, @@ -51,65 +53,215 @@ simpleTestRateParams = { 1.0, 0.5 } }; +using simpleTestType = MixedReactionsParameters< double, int, int, 5, 2, 2 >; constexpr -MixedReactionsParameters< double, int, int, 18, 11 > -carbonateSystem = +simpleTestType +simpleTestRateParams = { // stoichiometric matrix - {// OH- CO2 CO3-2 H2CO3 CaHCO3+ CaCO3 CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+ + { + { -2, 1, 1, 0, 0 }, + { 0, 0, -1, -1, 2 } + }, + // equilibrium constants + { 1.0, 1.0 }, + // forward rate constants + { 1.0, 0.5 }, + // reverse rate constants + { 1.0, 0.5 } +}; + +namespace carbonate +{ + +constexpr CArrayWrapper stoichMatrix = + { // OH- CO2 CO3-2 H2CO3 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- CaCO3 H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+ { -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O { 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3- { 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0 }, // CO3-2 + H+ = HCO3- { 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // H2CO3 = H+ + HCO3- { 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0 }, // CaHCO3+ = Ca+2 + HCO3- - { 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 1, 1, 0, 0, 0, 0 }, // CaCO3 + H+ = Ca+2 + HCO3- - { 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2 - { 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl- - { 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl- - { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2 - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 1 } // NaSO4- = Na+ + SO4-2 - }, - // equilibrium constants - { 9.77E+13, // OH- + H+ = H2O + { 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2 + { 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl- + { 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl- + { 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2 + { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 1 }, // NaSO4- = Na+ + SO4-2 + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0 } // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + // { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 1, 0, 0, 0, 0 } // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) + }; + +constexpr CArrayWrapper equilibriumConstants = + { + 9.77E+13, // OH- + H+ = H2O 4.37E-07, // CO2 + H2O = H+ + HCO3- 2.14E+10, // CO3-2 + H+ = HCO3- 1.70E-04, // H2CO3 = H+ + HCO3- 8.13E-02, // CaHCO3+ = Ca+2 + HCO3- - 1.17E+07, // CaCO3 + H+ = Ca+2 + HCO3- 6.92E-03, // CaSO4 = Ca+2 + SO4-2 4.68E+00, // CaCl+ = Ca+2 + Cl- 3.98E+00, // CaCl2 = Ca+2 + 2Cl- 3.72E-03, // MgSO4 = Mg+2 + SO4-2 - 1.51E-01 }, // NaSO4- = Na+ + SO4-2 - // forward rate constants - { 1.4e11, // OH- + H+ = H2O - 0.039, // CO2 + H2O = H+ + HCO3- - 1.0e10, // CO3-2 + H+ = HCO3- - 0.57, // H2CO3 = H+ + HCO3- - 1.5e6, // CaHCO3+ = Ca+2 + HCO3- - 1.0e5, // CaCO3 + H+ = Ca+2 + HCO3- - 1.0e5, // CaSO4 = Ca+2 + SO4-2 - 1.0e8, // CaCl+ = Ca+2 + Cl- - 1.0e7, // CaCl2 = Ca+2 + 2Cl- - 1.0e5, // MgSO4 = Mg+2 + SO4-2 - 1.0e7 // NaSO4- = Na+ + SO4-2 - }, - // reverse rate constants - { 1.43E-03, - 8.92E+04, - 4.67E-01, - 3.35E+03, - 1.85E+07, - 8.55E-03, - 1.45E+07, - 2.14E+07, - 2.51E+06, - 2.69E+07, - 6.62E+07 } + 1.51E-01, // NaSO4- = Na+ + SO4-2 + 1.17E+07 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + // 1 + }; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) + +constexpr CArrayWrapper forwardRates = + { + 1.4e11, // OH- + H+ = H2O + 0.039, // CO2 + H2O = H+ + HCO3- + 1.0e10, // CO3-2 + H+ = HCO3- + 0.57, // H2CO3 = H+ + HCO3- + 1.5e6, // CaHCO3+ = Ca+2 + HCO3- + 1.0e5, // CaSO4 = Ca+2 + SO4-2 + 1.0e8, // CaCl+ = Ca+2 + Cl- + 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 + }; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) + +constexpr CArrayWrapper reverseRates = + { 1.43E-03, // OH- + H+ = H2O + 8.92E+04, // CO2 + H2O = H+ + HCO3- + 4.67E-01, // CO3-2 + H+ = HCO3- + 3.35E+03, // H2CO3 = H+ + HCO3- + 1.85E+07, // CaHCO3+ = Ca+2 + HCO3- + 1.45E+07, // CaSO4 = Ca+2 + SO4-2 + 2.14E+07, // CaCl+ = Ca+2 + Cl- + 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- + // 1 // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) + }; + } + using carbonateSystemAllKineticType = MixedReactionsParameters< double, int, int, 18, 11, 0 >; + using carbonateSystemAllEquilibriumType = MixedReactionsParameters< double, int, int, 18, 11, 11 >; + using carbonateSystemType = MixedReactionsParameters< double, int, int, 18, 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::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates ); + +// ################################## Ultramafic rxn set ################################## +namespace ultramafics +{ + +constexpr CArrayWrapper stoichMatrix = +{ // OH- CO2(aq) CO3-- Mg2OH+++ Mg4(OH)++++ MgOH+ Mg2CO3++ MgCO3(aq) MgHCO3+ Mg(H3SiO4)2 MgH2SiO4 MgH3SiO4+ H2SiO4-- H3SiO4- H4(H2SiO4)---- H6(H2SiO4)-- Mg2SiO4 MgCO3 SiO2 Mg3Si2O5(OH)4 Mg(OH)2 H+ HCO3- Mg++ SiO2(aq) + { -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0 }, // OH- + H+ = H2O + { 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0 }, // CO2(aq) + H2O = HCO3- + H+ + { 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0 }, // CO3-- + H+ = HCO3- + { 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 2, 0 }, // Mg2OH+++ + H+ = 2Mg++ + H2O + { 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 0, 4, 0 }, // Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O + { 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0 }, // MgOH+ + H+ = Mg++ + H2O + { 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 2, 0 }, // Mg2CO3++ + H+ = 2Mg++ + HCO3- + { 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 1, 0 }, // MgCO3 + H+ = Mg++ + HCO3- + { 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0 }, // MgHCO3+ = Mg++ + HCO3- + { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 1, 1 }, // Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 1, 1 }, // MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 1 }, // MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 1 }, // H2SiO4-- + 2H+ = SiO2(aq) + 2H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1 }, // H3SiO4- + H+ = SiO2(aq) + 2H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -4, 0, 0, 4 }, // H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -2, 0, 0, 4 }, // H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -4, 0, 2, 1 }, // Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 1, 1, 0 }, // MgCO3 + H+ = Mg++ + HCO3- + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1 }, // SiO2 = SiO2(aq) + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -6, 0, 3, 2 }, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -2, 0, 1, 0 } // Mg(OH)2 + 2H+ = Mg++ + 2H2O + }; + +// 2Mg2SiO4 + 3H2O → Mg3Si2O5(OH)4 + Mg(OH)2 Serpentinization reaction + +constexpr CArrayWrapper equilibriumConstants = + { + 9.89E+13, // OH- + H+ = H2O + 4.42E-07, // CO2(aq) + H2O = HCO3- + H+ + 2.23E+10, // CO3-- + H+ = HCO3- + 2.32E+13, // Mg2OH+++ + H+ = 2Mg++ + H2O + 4.47E+39, // Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O + 6.18E+11, // MgOH+ + H+ = Mg++ + H2O + 7.66E+06, // Mg2CO3++ + H+ = 2Mg++ + HCO3- + 2.67E+07, // MgCO3 + H+ = Mg++ + HCO3- + 9.77E-02, // MgHCO3+ = Mg++ + HCO3- + 3.45E+14, // Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O + 9.49E+16, // MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O + 1.96E+08, // MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O + 8.08E+22, // H2SiO4-- + 2H+ = SiO2(aq) + 2H2O + 6.44E+09, // H3SiO4- + H+ = SiO2(aq) + 2H2O + 5.39E+35, // H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O + 2.72E+13, // H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O + 1.40E+28, // Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O + 2.73E+02, // MgCO3 + H+ = Mg++ + HCO3- + 1.93E-03, // SiO2 = SiO2(aq) + 3.54E+31, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O + 2.75E+16 // Mg(OH)2 + 2H+ = Mg++ + 2H2O + }; + +constexpr CArrayWrapper forwardRates = + { + 1.00E+10, // OH- + H+ = H2O + 1.00E+10, // CO2(aq) + H2O = HCO3- + H+ + 1.00E+10, // CO3-- + H+ = HCO3- + 1.00E+10, // Mg2OH+++ + H+ = 2Mg++ + H2O + 1.00E+10, // Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O + 1.00E+10, // MgOH+ + H+ = Mg++ + H2O + 1.00E+10, // Mg2CO3++ + H+ = 2Mg++ + HCO3- + 1.00E+10, // MgCO3 + H+ = Mg++ + HCO3- + 1.00E+10, // MgHCO3+ = Mg++ + HCO3- + 1.00E+10, // Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O + 1.00E+10, // MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O + 1.00E+10, // MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O + 1.00E+10, // H2SiO4-- + 2H+ = SiO2(aq) + 2H2O + 1.00E+10, // H3SiO4- + H+ = SiO2(aq) + 2H2O + 1.00E+10, // H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O + 1.00E+10, // H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O + 2.29E-11, // Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O + 4.57E-10, // MgCO3 + H+ = Mg++ + HCO3- + 1.70E-13, // SiO2 = SiO2(aq) + 1.00E-12, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O + 5.75E-09 // Mg(OH)2 + 2H+ = Mg++ + 2H2O + }; + +constexpr CArrayWrapper reverseRates = + { + 1.00E+10, // OH- + H+ = H2O + 1.00E+10, // CO2(aq) + H2O = HCO3- + H+ + 1.00E+10, // CO3-- + H+ = HCO3- + 1.00E+10, // Mg2OH+++ + H+ = 2Mg++ + H2O + 1.00E+10, // Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O + 1.00E+10, // MgOH+ + H+ = Mg++ + H2O + 1.00E+10, // Mg2CO3++ + H+ = 2Mg++ + HCO3- + 1.00E+10, // MgCO3 + H+ = Mg++ + HCO3- + 1.00E+10, // MgHCO3+ = Mg++ + HCO3- + 1.00E+10, // Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O + 1.00E+10, // MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O + 1.00E+10, // MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O + 1.00E+10, // H2SiO4-- + 2H+ = SiO2(aq) + 2H2O + 1.00E+10, // H3SiO4- + H+ = SiO2(aq) + 2H2O + 1.00E+10, // H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O + 1.00E+10, // H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O + 1.65E-39, // Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O + 1.67E-12, // MgCO3 + H+ = Mg++ + HCO3- + 8.78E-11, // SiO2 = SiO2(aq) + 2.83E-44, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O + 2.10E-25 // Mg(OH)2 + 2H+ = Mg++ + 2H2O + }; }; + using ultramaficSystemAllKineticType = MixedReactionsParameters< double, int, int, 25, 21, 0 >; + using ultramaficSystemAllEquilibriumType = MixedReactionsParameters< double, int, int, 25, 21, 21 >; + using ultramaficSystemType = 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 ); -// *****UNCRUSTIFY-ON****** +// UNCRUSTIFY-ON } // namespace bulkGeneric } // namespace hpcReact diff --git a/src/reactions/bulkGeneric/SpeciesUtilities.hpp b/src/reactions/bulkGeneric/SpeciesUtilities.hpp index 0f99ebd..be205c4 100644 --- a/src/reactions/bulkGeneric/SpeciesUtilities.hpp +++ b/src/reactions/bulkGeneric/SpeciesUtilities.hpp @@ -3,6 +3,7 @@ #include "common/macros.hpp" #include #include +#include namespace hpcReact { @@ -26,9 +27,13 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params, ARRAY_1D & logSecondarySpeciesConcentrations, FUNC && derivativeFunc ) { - constexpr int numSpecies = PARAMS_DATA::numSpecies; - constexpr int numSecondarySpecies = PARAMS_DATA::numReactions; - constexpr int numPrimarySpecies = numSpecies - numSecondarySpecies; + constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + + for( INDEX_TYPE i = 0; i < numSecondarySpecies; ++i ) + { + logSecondarySpeciesConcentrations[i] = 0.0; + } for( int j=0; j +HPCREACT_HOST_DEVICE +inline +void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params, + ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations, + ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations, + ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations, + ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_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; + } + } + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] ); + aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i; + dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i; + for( int j = 0; j < numSecondarySpecies; ++j ) + { + REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] ); + aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j; + for( int k=0; k( params, @@ -135,14 +200,5 @@ 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, -// typename ARRAY_2D > -// void calculateAggregateBasedResidualAndJacobianWrtLogC( ) - } // namespace bulkGeneric } // namespace hpcReact diff --git a/src/reactions/bulkGeneric/unitTests/CMakeLists.txt b/src/reactions/bulkGeneric/unitTests/CMakeLists.txt index 6a30651..eabb67e 100644 --- a/src/reactions/bulkGeneric/unitTests/CMakeLists.txt +++ b/src/reactions/bulkGeneric/unitTests/CMakeLists.txt @@ -3,6 +3,7 @@ set( testSourceFiles testEquilibriumReactions.cpp testKineticReactions.cpp testSpeciesUtilities.cpp + testMixedReactions.cpp ) set( dependencyList hpcReact gtest ) diff --git a/src/reactions/bulkGeneric/unitTests/testEquilibriumReactions.cpp b/src/reactions/bulkGeneric/unitTests/testEquilibriumReactions.cpp index 3cadb3b..e83ecf9 100644 --- a/src/reactions/bulkGeneric/unitTests/testEquilibriumReactions.cpp +++ b/src/reactions/bulkGeneric/unitTests/testEquilibriumReactions.cpp @@ -21,16 +21,16 @@ template< typename REAL_TYPE, int RESIDUAL_FORM, typename PARAMS_DATA > void computeResidualAndJacobianTest( PARAMS_DATA const & params, - REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies], - REAL_TYPE const (&expectedResidual)[PARAMS_DATA::numReactions], - REAL_TYPE const (&expectedJacobian)[PARAMS_DATA::numReactions][PARAMS_DATA::numReactions] ) + REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedResidual)[PARAMS_DATA::numReactions()], + REAL_TYPE const (&expectedJacobian)[PARAMS_DATA::numReactions()][PARAMS_DATA::numReactions()] ) { using EquilibriumReactionsType = EquilibriumReactions< REAL_TYPE, int, int >; - constexpr int numSpecies = PARAMS_DATA::numSpecies; - constexpr int numReactions = PARAMS_DATA::numReactions; + constexpr int numSpecies = PARAMS_DATA::numSpecies(); + constexpr int numReactions = PARAMS_DATA::numReactions(); double const temperature = 298.15; double speciesConcentration[numSpecies]; @@ -99,14 +99,14 @@ template< typename REAL_TYPE, int RESIDUAL_FORM, typename PARAMS_DATA > void testEnforceEquilibrium( PARAMS_DATA const & params, - REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies], - REAL_TYPE const (&expectedSpeciesConcentrations)[PARAMS_DATA::numSpecies] ) + REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedSpeciesConcentrations)[PARAMS_DATA::numSpecies()] ) { using EquilibriumReactionsType = EquilibriumReactions< REAL_TYPE, int, int >; - constexpr int numSpecies = PARAMS_DATA::numSpecies; + constexpr int numSpecies = PARAMS_DATA::numSpecies(); double const temperature = 298.15; double speciesConcentration0[numSpecies]; @@ -148,7 +148,7 @@ TEST( testEquilibriumReactions, testEnforceEquilibrium ) //****************************************************************************** -TEST( testEquilibriumReactions, testCarbonateSystem ) +TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium ) { double const initialSpeciesConcentration[18] = { @@ -157,12 +157,12 @@ TEST( testEquilibriumReactions, testCarbonateSystem ) 1.0e-16, // CO3-2 1.0e-16, // H2CO3 1.0e-16, // CaHCO3+ - 1.0e-16, // CaCO3 1.0e-16, // CaSO4 1.0e-16, // CaCl+ 1.0e-16, // CaCl2 1.0e-16, // MgSO4 1.0e-16, // NaSO4- + 1.0e-16, // CaCO3 3.76e-1, // H+ 3.76e-1, // HCO3- 3.87e-2, // Ca+2 @@ -178,12 +178,12 @@ TEST( testEquilibriumReactions, testCarbonateSystem ) 3.956656978189456e-11, // CO3-2 9.629355924567627e-04, // H2CO3 6.739226982791492e-05, // CaHCO3+ - 1.065032288527957e-09, // CaCO3 5.298329882666738e-03, // CaSO4 5.844517547638333e-03, // CaCl+ 1.277319392670652e-02, // CaCl2 6.618125707964991e-03, // MgSO4 1.769217213462983e-02, // NaSO4- + 1.065032288527957e-09, // CaCO3 4.396954721488358e-04, // H+ 3.723009698453808e-04, // HCO3- 1.471656530812871e-02, // Ca+2 @@ -194,33 +194,30 @@ TEST( testEquilibriumReactions, testCarbonateSystem ) }; std::cout<<" RESIDUAL_FORM 0:"<( carbonateSystem.equilibriumReactionsParameters(), + testEnforceEquilibrium< double, 0 >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), initialSpeciesConcentration, expectedSpeciesConcentrations ); // std::cout<<" RESIDUAL_FORM 1:"<( carbonateSystem.equilibriumReactionsParameters(), + // testEnforceEquilibrium< double, 1 >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), // initialSpeciesConcentration, // expectedSpeciesConcentrations ); std::cout<<" RESIDUAL_FORM 2:"<( carbonateSystem.equilibriumReactionsParameters(), + testEnforceEquilibrium< double, 2 >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), initialSpeciesConcentration, expectedSpeciesConcentrations ); } -TEST( testEquilibriumReactions, testCarbonateSystem2 ) +TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 ) { using EquilibriumReactionsType = EquilibriumReactions< double, int, int >; - constexpr int numSpecies = carbonateSystem.numSpecies; - constexpr int numReactions = carbonateSystem.numReactions; - constexpr int numPrimarySpecies = numSpecies - numReactions; -// constexpr int numSecondarySpecies = numReactions; + constexpr int numPrimarySpecies = carbonateSystemAllEquilibrium.numPrimarySpecies(); double const initialPrimarySpeciesConcentration[numPrimarySpecies] = { @@ -248,7 +245,7 @@ TEST( testEquilibriumReactions, testCarbonateSystem2 ) double logPrimarySpeciesConcentration[numPrimarySpecies]; EquilibriumReactionsType::enforceEquilibrium_Aggregate( 0, - carbonateSystem, + carbonateSystemAllEquilibrium, logInitialPrimarySpeciesConcentration, logPrimarySpeciesConcentration ); diff --git a/src/reactions/bulkGeneric/unitTests/testKineticReactions.cpp b/src/reactions/bulkGeneric/unitTests/testKineticReactions.cpp index 3945511..ae089a4 100644 --- a/src/reactions/bulkGeneric/unitTests/testKineticReactions.cpp +++ b/src/reactions/bulkGeneric/unitTests/testKineticReactions.cpp @@ -22,17 +22,17 @@ template< typename REAL_TYPE, bool LOGE_CONCENTRATION, typename PARAMS_DATA > void computeReactionRatesTest( PARAMS_DATA const & params, - REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies], - REAL_TYPE const (&expectedReactionRates)[PARAMS_DATA::numReactions], - REAL_TYPE const (&expectedReactionRatesDerivatives)[PARAMS_DATA::numReactions][PARAMS_DATA::numSpecies] ) + REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedReactionRates)[PARAMS_DATA::numReactions()], + REAL_TYPE const (&expectedReactionRatesDerivatives)[PARAMS_DATA::numReactions()][PARAMS_DATA::numSpecies()] ) { using KineticReactionsType = KineticReactions< REAL_TYPE, int, int, LOGE_CONCENTRATION >; - constexpr int numSpecies = PARAMS_DATA::numSpecies; - constexpr int numReactions = PARAMS_DATA::numReactions; + constexpr int numSpecies = PARAMS_DATA::numSpecies(); + constexpr int numReactions = PARAMS_DATA::numReactions(); double const temperature = 298.15; double speciesConcentration[numSpecies]; @@ -93,24 +93,24 @@ void computeReactionRatesTest( PARAMS_DATA const & params, //****************************************************************************** -TEST( testKineticReactions, computeReactionRatesTest_simpleTestRateParams ) +TEST( testKineticReactions, computeReactionRatesTest_simpleKineticTestRateParams ) { double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; 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 >( simpleTestRateParams.kineticReactionsParameters(), + computeReactionRatesTest< double, false >( simpleKineticTestRateParams.kineticReactionsParameters(), initialSpeciesConcentration, expectedReactionRates, expectedReactionRatesDerivatives ); - computeReactionRatesTest< double, true >( simpleTestRateParams, + computeReactionRatesTest< double, true >( simpleKineticTestRateParams.kineticReactionsParameters(), initialSpeciesConcentration, expectedReactionRates, expectedReactionRatesDerivatives ); } -TEST( testKineticReactions, computeReactionRatesTest_carbonateSystem ) +TEST( testKineticReactions, computeReactionRatesTest_carbonateSystemAllKinetic ) { double const initialSpeciesConcentration[18] = { @@ -119,12 +119,12 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystem ) 1.0e-16, // CO3-2 1.0e-16, // H2CO3 1.0e-16, // CaHCO3+ - 1.0e-16, // CaCO3 1.0e-16, // CaSO4 1.0e-16, // CaCl+ 1.0e-16, // CaCl2 1.0e-16, // MgSO4 1.0e-16, // NaSO4- + 1.0e-16, // CaCO3 3.76e-1, // H+ 3.76e-1, // HCO3- 3.87e-2, // Ca+2 @@ -134,8 +134,18 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystem ) 1.09 // Na+1 }; - double const expectedReactionRates[11] = { -0.001424736, -12610.7392, -0.175591624, -473.6096, -269197.19999999984, -0.00012441275624000003, -18012.914999999986, -1.56526019999999e6, - -346983.07769999903, -14247.58499999999, -2.316271799999999e6 }; + double const expectedReactionRates[11] = { -0.001424736, // OH- + H+ = H2O + -12610.7392, // CO2 + H2O = H+ + HCO3- + -0.175591624, // CO3-2 + H+ = HCO3- + -473.6096, // H2CO3 = H+ + HCO3- + -269197.19999999984, // CaHCO3+ = Ca+2 + HCO3- + -18012.914999999986, // CaSO4 = Ca+2 + SO4-2 + -1.56526019999999e6, // CaCl+ = Ca+2 + Cl- + -346983.07769999903, // CaCl2 = Ca+2 + 2Cl- + -14247.58499999999, // MgSO4 = Mg+2 + SO4-2 + -2.316271799999999e6, // NaSO4- = Na+ + SO4-2 + -0.00012441275624000003 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + }; double const expectedReactionRatesDerivatives[11][18] = { { 5.264e10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.000014, 0, 0, 0, 0, 0, 0 }, @@ -143,19 +153,20 @@ TEST( testKineticReactions, computeReactionRatesTest_carbonateSystem ) { 0, 0, 3.76e9, 0, 0, 0, 0, 0, 0, 0, 0, 1.e-6, -0.467, 0, 0, 0, 0, 0 }, { 0, 0, 0, 0.57, 0, 0, 0, 0, 0, 0, 0, -1259.6, -1259.6, 0, 0, 0, 0, 0 }, { 0, 0, 0, 0, 1.5e6, 0, 0, 0, 0, 0, 0, 0, -715950., -6.956e6, 0, 0, 0, 0 }, - { 0, 0, 0, 0, 0, 37600., 0, 0, 0, 0, 0, 1.e-11, -0.000330885, -0.0032148000000000003, 0, 0, 0, 0 }, - { 0, 0, 0, 0, 0, 0, 100000., 0, 0, 0, 0, 0, 0, -465449.99999999994, -561150., 0, 0, 0 }, - { 0, 0, 0, 0, 0, 0, 0, 1.e8, 0, 0, 0, 0, 0, -4.0446e7, 0, -828180., 0, 0 }, - { 0, 0, 0, 0, 0, 0, 0, 0, 1.e7, 0, 0, 0, 0, -8.965971e6, 0, -367177.86, 0, 0 }, - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 100000., 0, 0, 0, 0, -443850., 0, -863489.9999999999, 0 }, - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.e7, 0, 0, 0, -7.2158e7, 0, 0, -2.12502e6 } + { 0, 0, 0, 0, 0, 100000., 0, 0, 0, 0, 0, 0, 0, -465449.99999999994, -561150., 0, 0, 0 }, + { 0, 0, 0, 0, 0, 0, 1.e8, 0, 0, 0, 0, 0, 0, -4.0446e7, 0, -828180., 0, 0 }, + { 0, 0, 0, 0, 0, 0, 0, 1.e7, 0, 0, 0, 0, 0, -8.965971e6, 0, -367177.86, 0, 0 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 100000., 0, 0, 0, 0, 0, -443850., 0, -863489.9999999999, 0 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.e7, 0, 0, 0, 0, -7.2158e7, 0, 0, -2.12502e6 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 37600., 1.e-11, -0.000330885, -0.0032148000000000003, 0, 0, 0, 0 } + }; - computeReactionRatesTest< double, false >( carbonateSystem, + computeReactionRatesTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(), initialSpeciesConcentration, expectedReactionRates, expectedReactionRatesDerivatives ); - computeReactionRatesTest< double, true >( carbonateSystem, + computeReactionRatesTest< double, true >( carbonateSystemAllKinetic.kineticReactionsParameters(), initialSpeciesConcentration, expectedReactionRates, expectedReactionRatesDerivatives ); @@ -168,9 +179,9 @@ template< typename REAL_TYPE, bool LOGE_CONCENTRATION, typename PARAMS_DATA > void computeSpeciesRatesTest( PARAMS_DATA const & params, - REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies], - REAL_TYPE const (&expectedSpeciesRates)[PARAMS_DATA::numSpecies], - REAL_TYPE const (&expectedSpeciesRatesDerivatives)[PARAMS_DATA::numSpecies][PARAMS_DATA::numSpecies] ) + REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedSpeciesRates)[PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedSpeciesRatesDerivatives)[PARAMS_DATA::numSpecies()][PARAMS_DATA::numSpecies()] ) { using KineticReactionsType = KineticReactions< REAL_TYPE, @@ -178,7 +189,7 @@ void computeSpeciesRatesTest( PARAMS_DATA const & params, int, LOGE_CONCENTRATION >; - constexpr int numSpecies = PARAMS_DATA::numSpecies; + constexpr int numSpecies = PARAMS_DATA::numSpecies(); double const temperature = 298.15; double speciesConcentration[numSpecies]; @@ -225,7 +236,7 @@ void computeSpeciesRatesTest( PARAMS_DATA const & params, } } -TEST( testKineticReactions, computeSpeciesRatesTest_simpleTestRateParams ) +TEST( testKineticReactions, computeSpeciesRatesTest_simpleKineticTestRateParams ) { double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; double const expectedSpeciesRates[5] = { -2.0, 1.0, 0.75, -0.25, 0.5 }; @@ -235,19 +246,19 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleTestRateParams ) { 0.0, 0.0, -0.5, -0.25, 0.0 }, { 0.0, 0.0, 1.0, 0.5, 0.0 } }; - computeSpeciesRatesTest< double, false >( simpleTestRateParams, + computeSpeciesRatesTest< double, false >( simpleKineticTestRateParams.kineticReactionsParameters(), initialSpeciesConcentration, expectedSpeciesRates, expectedSpeciesRatesDerivatives ); - computeSpeciesRatesTest< double, true >( simpleTestRateParams, + computeSpeciesRatesTest< double, true >( simpleKineticTestRateParams.kineticReactionsParameters(), initialSpeciesConcentration, expectedSpeciesRates, expectedSpeciesRatesDerivatives ); } -// TEST( testKineticReactions, computeSpeciesRatesTest_carbonateSystem ) +// TEST( testKineticReactions, computeSpeciesRatesTest_carbonateSystemAllKinetic ) // { // double const initialSpeciesConcentration[18] = // { @@ -274,7 +285,7 @@ TEST( testKineticReactions, computeSpeciesRatesTest_simpleTestRateParams ) // double const expectedSpeciesRates[18] = { 0 }; // double const expectedSpeciesRatesDerivatives[18][18] = {{ 0}}; -// computeSpeciesRatesTest< double, false >( carbonateSystem, +// computeSpeciesRatesTest< double, false >( carbonateSystemAllKinetic, // initialSpeciesConcentration, // expectedSpeciesRates, // expectedSpeciesRatesDerivatives ); @@ -289,15 +300,15 @@ template< typename REAL_TYPE, void timeStepTest( PARAMS_DATA const & params, REAL_TYPE const dt, int const numSteps, - REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies], - REAL_TYPE const (&expectedSpeciesConcentrations)[PARAMS_DATA::numSpecies] ) + REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numSpecies()], + REAL_TYPE const (&expectedSpeciesConcentrations)[PARAMS_DATA::numSpecies()] ) { using KineticReactionsType = KineticReactions< double, int, int, LOGE_CONCENTRATION >; - constexpr int numSpecies = PARAMS_DATA::numSpecies; + constexpr int numSpecies = PARAMS_DATA::numSpecies(); double const temperature = 298.15; double speciesConcentration[numSpecies]; @@ -357,14 +368,14 @@ TEST( testKineticReactions, testTimeStep ) double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 }; double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01, 7.02014627734060e-01, 5.95970744531880e-01 }; - timeStepTest< double, false >( simpleTestRateParams, + timeStepTest< double, false >( simpleKineticTestRateParams.kineticReactionsParameters(), 2.0, 10, initialSpeciesConcentration, expectedSpeciesConcentrations ); // ln(c) as the primary variable results in a singular system. - // timeStepTest< double, true >( simpleTestRateParams, + // timeStepTest< double, true >( simpleKineticTestRateParams, // 2.0, // 10, // initialSpeciesConcentration, @@ -372,7 +383,7 @@ TEST( testKineticReactions, testTimeStep ) } -TEST( testKineticReactions, testTimeStep_carbonateSystem ) +TEST( testKineticReactions, testTimeStep_carbonateSystemAllKinetic ) { double const initialSpeciesConcentration[18] = { @@ -381,12 +392,12 @@ TEST( testKineticReactions, testTimeStep_carbonateSystem ) 1.0e-16, // CO3-2 1.0e-16, // H2CO3 1.0e-16, // CaHCO3+ - 1.0e-16, // CaCO3 1.0e-16, // CaSO4 1.0e-16, // CaCl+ 1.0e-16, // CaCl2 1.0e-16, // MgSO4 1.0e-16, // NaSO4- + 1.0e-16, // CaCO3 3.76e-1, // H+ 3.76e-1, // HCO3- 3.87e-2, // Ca+2 @@ -402,12 +413,12 @@ TEST( testKineticReactions, testTimeStep_carbonateSystem ) 3.956656978189456e-11, // CO3-2 9.629355924567627e-04, // H2CO3 6.739226982791492e-05, // CaHCO3+ - 1.065032288527957e-09, // CaCO3 5.298329882666738e-03, // CaSO4 5.844517547638333e-03, // CaCl+ 1.277319392670652e-02, // CaCl2 6.618125707964991e-03, // MgSO4 1.769217213462983e-02, // NaSO4- + 1.065032288527957e-09, // CaCO3 4.396954721488358e-04, // H+ 3.723009698453808e-04, // HCO3- 1.471656530812871e-02, // Ca+2 @@ -417,14 +428,14 @@ TEST( testKineticReactions, testTimeStep_carbonateSystem ) 1.072307827865370e+00 // Na+1 }; - timeStepTest< double, false >( carbonateSystem, + timeStepTest< double, false >( carbonateSystemAllKinetic.kineticReactionsParameters(), 2.0, 100000, initialSpeciesConcentration, expectedSpeciesConcentrations ); // ln(c) as the primary variable results in a singular system. - // timeStepTest< double, true >( simpleTestRateParams, + // timeStepTest< double, true >( simpleKineticTestRateParams, // 2.0, // 10, // initialSpeciesConcentration, diff --git a/src/reactions/bulkGeneric/unitTests/testMixedReactions.cpp b/src/reactions/bulkGeneric/unitTests/testMixedReactions.cpp new file mode 100644 index 0000000..6ae37d3 --- /dev/null +++ b/src/reactions/bulkGeneric/unitTests/testMixedReactions.cpp @@ -0,0 +1,178 @@ + +#include "../MixedEquilibriumKineticReactions.hpp" +#include "../EquilibriumReactions.hpp" +#include "../ParametersPredefined.hpp" +#include "common/macros.hpp" +#include "common/printers.hpp" +#include "common/nonlinearSolvers.hpp" +#include "common/pmpl.hpp" + +#include + + +using namespace hpcReact; +using namespace hpcReact::bulkGeneric; + +template< typename REAL_TYPE > +REAL_TYPE tolerance( REAL_TYPE const a, REAL_TYPE const b, REAL_TYPE const ndigits ) +{ + return std::numeric_limits< double >::epsilon() * std::max( fabs( a ), fabs( b ) ) * pow( 10, ndigits ); +} + +//****************************************************************************** +template< typename REAL_TYPE, + bool LOGE_CONCENTRATION, + typename PARAMS_DATA > +void timeStepTest( PARAMS_DATA const & params, + REAL_TYPE const dt, + int const numSteps, + REAL_TYPE const (&initialSpeciesConcentration)[PARAMS_DATA::numPrimarySpecies()], + REAL_TYPE const (&expectedSpeciesConcentrations)[PARAMS_DATA::numPrimarySpecies()] ) +{ + HPCREACT_UNUSED_VAR( expectedSpeciesConcentrations ); + + CArrayWrapper< REAL_TYPE, PARAMS_DATA::numPrimarySpecies() > primarySpeciesConcentration; + + for( int i = 0; i < PARAMS_DATA::numPrimarySpecies(); ++i ) + { + primarySpeciesConcentration[i] = initialSpeciesConcentration[i]; + } + + pmpl::genericKernelWrapper( PARAMS_DATA::numPrimarySpecies(), + primarySpeciesConcentration.data, + [=] HPCREACT_HOST_DEVICE ( auto * speciesConcentration ) + { + using MixedReactionsType = MixedEquilibriumKineticReactions< REAL_TYPE, + int, + int, + LOGE_CONCENTRATION >; + using EquilibriumReactionsType = EquilibriumReactions< REAL_TYPE, + int, + int >; + + // constexpr int numSpecies = PARAMS_DATA::numSpecies(); + constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies(); + constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies(); + constexpr int numKineticReactions = PARAMS_DATA::numKineticReactions(); + + // define variables + double const temperature = 298.15; + REAL_TYPE logPrimarySpeciesConcentration[numPrimarySpecies]; + // must use CArrayWrapper to ensure correct capture in the lambda functions + CArrayWrapper< REAL_TYPE, numSecondarySpecies > logSecondarySpeciesConcentration; + CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration; + CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration_n; + CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration; + CArrayWrapper< REAL_TYPE, numKineticReactions > reactionRates; + CArrayWrapper< REAL_TYPE, numKineticReactions, numPrimarySpecies > dReactionRates_dlogPrimarySpeciesConcentration; + CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregateSpeciesRates; + CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dAggregateSpeciesRates_dlogPrimarySpeciesConcentration; + + // Initialize species concentrations + for( int i = 0; i < numPrimarySpecies; ++i ) + { + logPrimarySpeciesConcentration[i] = ::log( speciesConcentration[i] ); + aggregatePrimarySpeciesConcentration[i] = speciesConcentration[i]; + } + + EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature, + carbonateSystem.equilibriumReactionsParameters(), + logPrimarySpeciesConcentration, + logPrimarySpeciesConcentration ); + + /// Time step loop + double time = 0.0; + for( int t = 0; t < numSteps; ++t ) + { + printf( "Timestep %d, Time = %.6e\n", t, time ); + + for( int i=0; i < numPrimarySpecies; ++i ) + { + aggregatePrimarySpeciesConcentration_n[i] = aggregatePrimarySpeciesConcentration[i]; + } + + auto computeResidualAndJacobian = [&] HPCREACT_HOST_DEVICE ( REAL_TYPE const (&X)[numPrimarySpecies], + REAL_TYPE ( &r )[numPrimarySpecies], + REAL_TYPE ( &J )[numPrimarySpecies][numPrimarySpecies] ) + { + MixedReactionsType::updateMixedSystem( temperature, + params, + X, + logSecondarySpeciesConcentration, + aggregatePrimarySpeciesConcentration, + dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration, + reactionRates, + dReactionRates_dlogPrimarySpeciesConcentration, + aggregateSpeciesRates, + dAggregateSpeciesRates_dlogPrimarySpeciesConcentration ); + + + for( int i = 0; i < numPrimarySpecies; ++i ) + { + r[i] = ( aggregatePrimarySpeciesConcentration[i] - aggregatePrimarySpeciesConcentration_n[i] ) - aggregateSpeciesRates[i] * dt; + for( int j = 0; j < numPrimarySpecies; ++j ) + { + J[i][j] = dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration[i][j] - dAggregateSpeciesRates_dlogPrimarySpeciesConcentration[i][j] * dt; + } + } + }; + + nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian ); + + time += dt; + } + for( int i = 0; i < numPrimarySpecies; ++i ) + { + speciesConcentration[i] = exp( logPrimarySpeciesConcentration[i] ); + } + } ); + + // Check results + for( int i = 0; i < PARAMS_DATA::numPrimarySpecies(); ++i ) + { + EXPECT_NEAR( primarySpeciesConcentration[ i ], expectedSpeciesConcentrations[ i ], 1.0e-8 * expectedSpeciesConcentrations[ i ] ); + } +} + +TEST( testMixedReactions, testTimeStep_carbonateSystem ) +{ + constexpr int numPrimarySpecies = carbonateSystemType::numPrimarySpecies(); + + double const initialAggregateSpeciesConcentration[numPrimarySpecies] = + { + 3.76e-3, // CaCO3 + 3.76e-1, // H+ + 3.76e-1, // HCO3- + 3.87e-2, // Ca+2 + 3.21e-2, // SO4-2 + 1.89, // Cl- + 1.65e-2, // Mg+2 + 1.09 // Na+1 + }; + + double const expectedSpeciesConcentrations[numPrimarySpecies] = + { + 3.3318075516669661e-05, // CaCO3 + 2.5894448848121536e-05, // H+ + 0.0062660162912796741, // HCO3- + 0.015741214773567921, // Ca+2 + 0.0024602709470074127, // SO4-2 + 1.8564927944498291, // Cl- + 0.0099316034080546619, // Mg+2 + 1.0725251492409775 // Na+1 + }; + + timeStepTest< double, true >( carbonateSystem, + 0.2, + 10, + initialAggregateSpeciesConcentration, + expectedSpeciesConcentrations ); + +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + int const result = RUN_ALL_TESTS(); + return result; +} diff --git a/src/reactions/bulkGeneric/unitTests/testSpeciesUtilities.cpp b/src/reactions/bulkGeneric/unitTests/testSpeciesUtilities.cpp index 889888c..1378453 100644 --- a/src/reactions/bulkGeneric/unitTests/testSpeciesUtilities.cpp +++ b/src/reactions/bulkGeneric/unitTests/testSpeciesUtilities.cpp @@ -11,10 +11,8 @@ using namespace hpcReact::bulkGeneric; TEST( testUtilities, test_calculateLogSecondarySpeciesConcentration ) { - constexpr int numReactions = carbonateSystem.numReactions; - constexpr int numSpecies = carbonateSystem.numSpecies; - constexpr int numPrimarySpecies = numSpecies - numReactions; - constexpr int numSecondarySpecies = numReactions; + constexpr int numPrimarySpecies = carbonateSystemAllEquilibrium.numPrimarySpecies(); + constexpr int numSecondarySpecies = carbonateSystemAllEquilibrium.numSecondarySpecies(); double const logPrimarySpeciesSolution[numPrimarySpecies] = { @@ -31,7 +29,7 @@ TEST( testUtilities, test_calculateLogSecondarySpeciesConcentration ) calculateLogSecondarySpeciesConcentration< double, int, - int >( carbonateSystem, + int >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), logPrimarySpeciesSolution, logSecondarySpeciesConcentrations ); @@ -42,12 +40,12 @@ TEST( testUtilities, test_calculateLogSecondarySpeciesConcentration ) 3.956656978189425e-11, 0.0009629355924566718, 0.00006739226982791149, - 1.065032288527949e-9, 0.005298329882666744, 0.005844517547638335, 0.012773193926706526, 0.006618125707964999, - 0.01769217213462983 + 0.01769217213462983, + 1.065032288527949e-9 }; for( int j=0; j( carbonateSystem, + int >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), logPrimarySpeciesSolution, logSecondarySpeciesConcentrations, dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations ); @@ -76,12 +74,13 @@ TEST( testUtilities, test_calculateLogSecondarySpeciesConcentration ) { -1, 1, 0, 0, 0, 0, 0 }, { 1, 1, 0, 0, 0, 0, 0 }, { 0, 1, 1, 0, 0, 0, 0 }, - { -1, 1, 1, 0, 0, 0, 0 }, { 0, 0, 1, 1, 0, 0, 0 }, { 0, 0, 1, 0, 1, 0, 0 }, { 0, 0, 1, 0, 2, 0, 0 }, { 0, 0, 0, 1, 0, 1, 0 }, - { 0, 0, 0, 1, 0, 0, 1 } + { 0, 0, 0, 1, 0, 0, 1 }, + { -1, 1, 1, 0, 0, 0, 0 } + }; for( int i=0; i dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations = {{{0.0}}}; + CArrayWrapper< double, numPrimarySpecies, numPrimarySpecies > dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations; + for( int i = 0; i < numPrimarySpecies; ++i ) + { + for( int k=0; k( carbonateSystem, + calculateAggregatePrimaryConcentrationsWrtLogC< double, int, int >( carbonateSystemAllEquilibrium.equilibriumReactionsParameters(), primarySpeciesSolution, aggregatePrimarySpeciesConcentration, dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations );