Skip to content

Commit c320232

Browse files
CusiniMCusiniM
committed
add MixedEquilibriumKineticReactions_impl.hpp
1 parent 9581b6a commit c320232

File tree

1 file changed

+195
-0
lines changed

1 file changed

+195
-0
lines changed

‎src/reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp‎

Lines changed: 195 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,195 @@
1+
#pragma once
2+
3+
#include "common/constants.hpp"
4+
#include "common/CArrayWrapper.hpp"
5+
#include "SpeciesUtilities.hpp"
6+
7+
8+
/** @file MixedEquilibriumKineticReactions_impl.hpp
9+
* @brief Header file for the MixedEquilibriumKineticReactions implementation.
10+
* @author HPC-REACT Team
11+
* @date 2025
12+
*/
13+
14+
namespace hpcReact
15+
{
16+
namespace bulkGeneric
17+
{
18+
19+
template< typename REAL_TYPE,
20+
typename INT_TYPE,
21+
typename INDEX_TYPE,
22+
bool LOGE_CONCENTRATION >
23+
template< typename PARAMS_DATA,
24+
typename ARRAY_1D_TO_CONST,
25+
typename ARRAY_1D,
26+
typename ARRAY_2D >
27+
HPCREACT_HOST_DEVICE inline void
28+
MixedEquilibriumKineticReactions< REAL_TYPE,
29+
INT_TYPE,
30+
INDEX_TYPE,
31+
LOGE_CONCENTRATION
32+
>::updateMixedSystem_impl( RealType const & temperature,
33+
PARAMS_DATA const & params,
34+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
35+
ARRAY_1D & logSecondarySpeciesConcentrations,
36+
ARRAY_1D & aggregatePrimarySpeciesConcentrations,
37+
ARRAY_2D & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
38+
ARRAY_1D & reactionRates,
39+
ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations,
40+
ARRAY_1D & aggregateSpeciesRates,
41+
ARRAY_2D & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations )
42+
{
43+
constexpr IntType numSpecies = PARAMS_DATA::numSpecies();
44+
constexpr IntType numSecondarySpecies = PARAMS_DATA::numReactions();
45+
constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
46+
47+
// 1. Compute new aggregate species from primary species
48+
calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE,
49+
INT_TYPE,
50+
INDEX_TYPE >( params,
51+
logPrimarySpeciesConcentrations,
52+
logSecondarySpeciesConcentrations,
53+
aggregatePrimarySpeciesConcentrations,
54+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations );
55+
56+
// 2. Compute the reaction rates for all kinetic reactions
57+
computeReactionRates( temperature,
58+
params,
59+
logPrimarySpeciesConcentrations,
60+
logSecondarySpeciesConcentrations,
61+
aggregatePrimarySpeciesConcentrations,
62+
dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
63+
reactionRates,
64+
dReactionRates_dLogPrimarySpeciesConcentrations );
65+
66+
67+
// 3. Compute aggregate species rates
68+
computeAggregateSpeciesRates( params,
69+
logPrimarySpeciesConcentrations,
70+
reactionRates,
71+
dReactionRates_dLogPrimarySpeciesConcentrations,
72+
aggregateSpeciesRates,
73+
dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations );
74+
75+
76+
}
77+
78+
template< typename REAL_TYPE,
79+
typename INT_TYPE,
80+
typename INDEX_TYPE,
81+
bool LOGE_CONCENTRATION >
82+
template< typename PARAMS_DATA,
83+
typename ARRAY_1D_TO_CONST,
84+
typename ARRAY_1D,
85+
typename ARRAY_2D >
86+
HPCREACT_HOST_DEVICE inline void
87+
MixedEquilibriumKineticReactions< REAL_TYPE,
88+
INT_TYPE,
89+
INDEX_TYPE,
90+
LOGE_CONCENTRATION
91+
>::computeReactionRates_impl( RealType const & temperature,
92+
PARAMS_DATA const & params,
93+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
94+
ARRAY_1D_TO_CONST const & logSecondarySpeciesConcentrations,
95+
ARRAY_1D & reactionRates,
96+
ARRAY_2D & dReactionRates_dPrimarySpeciesConcentrations )
97+
98+
{
99+
constexpr IntType numSpecies = PARAMS_DATA::numSpecies();
100+
constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
101+
constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
102+
constexpr IntType numKineticReactions = PARAMS_DATA::numKineticReactions();
103+
104+
RealType logSpeciesConcentration[numSpecies] {};
105+
for ( IntType i = 0; i < numPrimarySpecies; ++i )
106+
{
107+
logSpeciesConcentration[i] = logPrimarySpeciesConcentrations[i];
108+
}
109+
for ( IntType i = 0; i < numSecondarySpecies; ++i )
110+
{
111+
logSpeciesConcentration[i+numPrimarySpecies] = logSecondarySpeciesConcentrations[i];
112+
}
113+
114+
CArrayWrapper< RealType, numKineticReactions, numSpecies > reactionRatesDerivatives;
115+
kineticReactions::computeReactionRates( temperature,
116+
params,
117+
logSpeciesConcentration,
118+
reactionRates,
119+
reactionRatesDerivatives );
120+
121+
// Compute the reaction rates derivatives w.r.t. log primary species concentrations
122+
for( IntType i = 0; i < numKineticReactions; ++i )
123+
{
124+
for( IntType j = 0; j < numPrimarySpecies; ++j )
125+
{
126+
dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) = reactionRatesDerivatives( i, j );
127+
for( IntType k = 0; k < numSecondarySpecies; ++k )
128+
{
129+
RealType const dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = params.stoichiometricMatrix( k, j+numSecondarySpecies );
130+
131+
dReactionRates_dLogPrimarySpeciesConcentrations( i, j ) +=
132+
reactionRatesDerivatives( i, numPrimarySpecies + k ) * dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations( k, j );
133+
}
134+
}
135+
}
136+
}
137+
138+
template< typename REAL_TYPE,
139+
typename INT_TYPE,
140+
typename INDEX_TYPE,
141+
bool LOGE_CONCENTRATION >
142+
template< typename PARAMS_DATA,
143+
typename ARRAY_1D_TO_CONST,
144+
typename ARRAY_2D_TO_CONST,
145+
typename ARRAY_1D,
146+
typename ARRAY_2D,
147+
bool CALCULATE_DERIVATIVES >
148+
HPCREACT_HOST_DEVICE inline void
149+
MixedEquilibriumKineticReactions< REAL_TYPE,
150+
INT_TYPE,
151+
INDEX_TYPE,
152+
LOGE_CONCENTRATION
153+
>::computeAggregateSpeciesRates_impl( PARAMS_DATA const & params,
154+
ARRAY_1D_TO_CONST const & speciesConcentration,
155+
ARRAY_1D_TO_CONST const & reactionRates,
156+
ARRAY_2D_TO_CONST const & reactionRatesDerivatives,
157+
ARRAY_1D & aggregatesRates,
158+
ARRAY_2D & aggregatesRatesDerivatives )
159+
{
160+
constexpr IntType numSpecies = PARAMS_DATA::numSpecies();
161+
constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
162+
constexpr IntType numKineticReactions = PARAMS_DATA::numKineticReactions();
163+
constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
164+
165+
for( IntType i = 0; i < numPrimarySpecies; ++i )
166+
{
167+
speciesRates[i] = 0.0;
168+
if constexpr( CALCULATE_DERIVATIVES )
169+
{
170+
for( IntType j = 0; j < numPrimarySpecies; ++j )
171+
{
172+
speciesRatesDerivatives( i, j ) = 0.0;
173+
}
174+
}
175+
for( IntType r=0; r<numKineticReactions; ++r )
176+
{
177+
RealType const s_ir = params.stoichiometricMatrix( r, i );
178+
speciesRates[i] += s_ir * reactionRates[r];
179+
if constexpr( CALCULATE_DERIVATIVES )
180+
{
181+
for( IntType j = 0; j < numPrimarySpecies; ++j )
182+
{
183+
speciesRatesDerivatives( i, j ) += s_ir * reactionRatesDerivatives( r, j );
184+
}
185+
}
186+
}
187+
}
188+
189+
}
190+
191+
} // namespace bulkGeneric
192+
193+
} // namespace hpcReact
194+
195+
#include "common/macrosCleanup.hpp"

0 commit comments

Comments
 (0)