Skip to content

Commit 49013a9

Browse files
rrsettgastrrsettgast
committed
added log(C) option
1 parent e06483a commit 49013a9

File tree

3 files changed

+185
-123
lines changed

3 files changed

+185
-123
lines changed

‎src/reactions/bulkGeneric/KineticReactions.hpp‎

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -15,15 +15,13 @@ namespace bulkGeneric
1515

1616
template< typename REAL_TYPE,
1717
typename INT_TYPE,
18-
typename INDEX_TYPE >
18+
typename INDEX_TYPE,
19+
bool LOGE_CONCENTRATION >
1920
class KineticReactions
2021
{
2122
public:
2223

2324
using RealType = REAL_TYPE;
24-
// using RealDataArrayView1d = REAL_DATA_ARRAY_1D_VIEW_TYPE;
25-
// using RealConstDataArrayView1d = REAL_CONST_DATA_ARRAY_1D_VIEW_TYPE;
26-
// using RealDataArrayView2d = REAL_DATA_ARRAY_2D_VIEW_TYPE;
2725
using IntType = INT_TYPE;
2826
using IndexType = INDEX_TYPE;
2927

@@ -57,7 +55,7 @@ class KineticReactions
5755
ARRAY_1D_TO_CONST const & speciesConcentration,
5856
ARRAY_1D & reactionRates)
5957
{
60-
REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions][PARAMS_DATA::numSpecies] = { 0.0 };
58+
REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions][PARAMS_DATA::numSpecies] = { {0.0} };
6159
computeReactionRates_impl< PARAMS_DATA, false >( temperature,
6260
params,
6361
speciesConcentration,
@@ -93,7 +91,7 @@ class KineticReactions
9391
ARRAY_1D_TO_CONST const & speciesConcentration,
9492
ARRAY_1D & speciesRates )
9593
{
96-
double speciesRatesDerivatives;
94+
char speciesRatesDerivatives;
9795
computeSpeciesRates_impl< PARAMS_DATA, false >( temperature,
9896
params,
9997
speciesConcentration,

‎src/reactions/bulkGeneric/KineticReactions_impl.hpp‎

Lines changed: 102 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,8 @@ namespace bulkGeneric
2222

2323
template< typename REAL_TYPE,
2424
typename INT_TYPE,
25-
typename INDEX_TYPE >
25+
typename INDEX_TYPE,
26+
bool LOGE_CONCENTRATION >
2627
template< typename PARAMS_DATA,
2728
bool CALCULATE_DERIVATIVES,
2829
typename ARRAY_1D_TO_CONST,
@@ -31,7 +32,8 @@ template< typename PARAMS_DATA,
3132
HPCREACT_HOST_DEVICE inline void
3233
KineticReactions< REAL_TYPE,
3334
INT_TYPE,
34-
INDEX_TYPE
35+
INDEX_TYPE,
36+
LOGE_CONCENTRATION
3537
>::computeReactionRates_impl( RealType const & ,//temperature,
3638
PARAMS_DATA const & params,
3739
ARRAY_1D_TO_CONST const & speciesConcentration,
@@ -47,93 +49,122 @@ KineticReactions< REAL_TYPE,
4749
// loop over each reaction
4850
for( IntType r=0; r<PARAMS_DATA::numReactions; ++r )
4951
{
52+
// set reaction rate to zero
5053
reactionRates[r] = 0.0;
54+
// get/calculate the forward and reverse rate constants for this reaction
5155
RealType const forwardRateConstant = params.rateConstantForward(r) ;//* exp( -params.m_activationEnergy[r] / ( constants::R * temperature ) );
5256
RealType const reverseRateConstant = params.rateConstantReverse(r) ;
5357

54-
RealType productConcForward = 1.0;
55-
RealType productConcReverse = 1.0;
56-
57-
RealType dProductConcForward_dC[PARAMS_DATA::numSpecies];
58-
RealType dProductConcReverse_dC[PARAMS_DATA::numSpecies];
59-
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
60-
{
61-
dProductConcForward_dC[i] = 1.0;
62-
dProductConcReverse_dC[i] = 1.0;
63-
}
64-
65-
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
58+
if constexpr ( LOGE_CONCENTRATION )
6659
{
60+
RealType productConcForward = 0.0;
61+
RealType productConcReverse = 0.0;
62+
// build the products for the forward and reverse reaction rates
63+
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
64+
{
65+
RealType const s_ri = params.stoichiometricMatrix( r, i );
66+
67+
if( s_ri < 0.0 )
68+
{
69+
productConcForward += (-s_ri) * speciesConcentration[i];
70+
}
71+
else if( s_ri > 0.0 )
72+
{
73+
productConcReverse += s_ri * speciesConcentration[i];
74+
}
75+
}
76+
reactionRates[r] = forwardRateConstant * exp( productConcForward )
77+
- reverseRateConstant * exp( productConcReverse );
6778

68-
RealType const s_ri = params.stoichiometricMatrix( r, i );
69-
RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], abs(s_ri) ) : 0.0;
70-
71-
if( s_ri < 0.0 )
79+
if constexpr ( CALCULATE_DERIVATIVES )
7280
{
73-
productConcForward *= productTerm_i;
81+
for( IntType j=0; j<PARAMS_DATA::numSpecies; ++j )
82+
{
83+
reactionRatesDerivatives( r, j ) = 0.0;
84+
}
7485
}
75-
else if( s_ri > 0.0 )
86+
}
87+
else
88+
{
89+
// variables used to build the product terms for the forward and reverse reaction rates
90+
RealType productConcForward = 1.0;
91+
RealType productConcReverse = 1.0;
92+
93+
RealType dProductConcForward_dC[PARAMS_DATA::numSpecies];
94+
RealType dProductConcReverse_dC[PARAMS_DATA::numSpecies];
95+
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
7696
{
77-
productConcReverse *= productTerm_i;
97+
dProductConcForward_dC[i] = 1.0;
98+
dProductConcReverse_dC[i] = 1.0;
7899
}
79-
80-
if constexpr ( CALCULATE_DERIVATIVES )
100+
101+
// build the products for the forward and reverse reaction rates
102+
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
81103
{
82104

105+
RealType const s_ri = params.stoichiometricMatrix( r, i );
106+
RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], abs(s_ri) ) : 0.0;
107+
83108
if( s_ri < 0.0 )
84109
{
85-
for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j )
86-
{
87-
if( i==j )
88-
{
89-
dProductConcForward_dC[j] *= -s_ri * pow( speciesConcentration[i], -s_ri-1 );
90-
dProductConcReverse_dC[j] = 0.0;
91-
}
92-
else
93-
{
94-
dProductConcForward_dC[j] *= productTerm_i;
95-
}
96-
}
110+
productConcForward *= productTerm_i;
97111
}
98112
else if( s_ri > 0.0 )
99113
{
100-
for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j )
101-
{
102-
if( i==j )
114+
productConcReverse *= productTerm_i;
115+
}
116+
117+
if constexpr ( CALCULATE_DERIVATIVES )
118+
{
119+
120+
if( s_ri < 0.0 )
121+
{
122+
for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j )
103123
{
104-
dProductConcReverse_dC[j] *= s_ri * pow( speciesConcentration[i], s_ri-1 );
105-
dProductConcForward_dC[j] = 0.0;
124+
if( i==j )
125+
{
126+
dProductConcForward_dC[j] *= -s_ri * pow( speciesConcentration[i], -s_ri-1 );
127+
dProductConcReverse_dC[j] = 0.0;
128+
}
129+
else
130+
{
131+
dProductConcForward_dC[j] *= productTerm_i;
132+
}
106133
}
107-
else
108-
{
109-
dProductConcReverse_dC[j] *= productTerm_i;
134+
}
135+
else if( s_ri > 0.0 )
136+
{
137+
for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j )
138+
{
139+
if( i==j )
140+
{
141+
dProductConcReverse_dC[j] *= s_ri * pow( speciesConcentration[i], s_ri-1 );
142+
dProductConcForward_dC[j] = 0.0;
143+
}
144+
else
145+
{
146+
dProductConcReverse_dC[j] *= productTerm_i;
147+
}
110148
}
111149
}
112-
}
113-
else
114-
{
115-
dProductConcForward_dC[i] = 0.0;
116-
dProductConcReverse_dC[i] = 0.0;
150+
else
151+
{
152+
dProductConcForward_dC[i] = 0.0;
153+
dProductConcReverse_dC[i] = 0.0;
154+
}
117155
}
118156
}
119-
}
120-
121-
122-
reactionRates[r] = forwardRateConstant * productConcForward - reverseRateConstant * productConcReverse;
123-
// printf( " kf, kr, reactionRates[%2d] = % 8.4e, % 8.4e, % 8.4e\n", r, forwardRateConstant, reverseRateConstant, reactionRates[r] );
157+
reactionRates[r] = forwardRateConstant * productConcForward - reverseRateConstant * productConcReverse;
124158

125-
126-
if constexpr( CALCULATE_DERIVATIVES )
127-
{
128-
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
159+
if constexpr( CALCULATE_DERIVATIVES )
129160
{
130-
reactionRatesDerivatives( r, i ) = forwardRateConstant * dProductConcForward_dC[i] - reverseRateConstant * dProductConcReverse_dC[i];
131-
// printf( " dR%ddC%d = % 3.1f * %3.2f - %3.1f * % 3.2f = % 3.2f \n", r, i, forwardRateConstant, dProductConcForward_dC[i], reverseRateConstant, dProductConcReverse_dC[i], reactionRatesDerivatives( r, i ) );
161+
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
162+
{
163+
reactionRatesDerivatives( r, i ) = forwardRateConstant * dProductConcForward_dC[i] - reverseRateConstant * dProductConcReverse_dC[i];
164+
}
132165
}
133-
}
134-
135-
136-
}
166+
} // end of if constexpr ( LOGE_CONCENTRATION )
167+
} // end of loop over reactions
137168
}
138169

139170

@@ -146,7 +177,8 @@ KineticReactions< REAL_TYPE,
146177
// function to the reaction rate. Includes impact of temperature, concentration, surface area, volume fraction and porosity
147178
template< typename REAL_TYPE,
148179
typename INT_TYPE,
149-
typename INDEX_TYPE >
180+
typename INDEX_TYPE,
181+
bool LOGE_CONCENTRATION >
150182
template< typename PARAMS_DATA,
151183
bool CALCULATE_DERIVATIVES,
152184
typename ARRAY_1D_TO_CONST,
@@ -155,7 +187,8 @@ template< typename PARAMS_DATA,
155187
HPCREACT_HOST_DEVICE inline void
156188
KineticReactions< REAL_TYPE,
157189
INT_TYPE,
158-
INDEX_TYPE
190+
INDEX_TYPE,
191+
LOGE_CONCENTRATION
159192
>::computeSpeciesRates_impl( RealType const & temperature,
160193
PARAMS_DATA const & params,
161194
ARRAY_1D_TO_CONST const & speciesConcentration,
@@ -199,15 +232,17 @@ KineticReactions< REAL_TYPE,
199232

200233
template< typename REAL_TYPE,
201234
typename INT_TYPE,
202-
typename INDEX_TYPE >
235+
typename INDEX_TYPE,
236+
bool LOGE_CONCENTRATION >
203237
template< typename PARAMS_DATA,
204238
typename ARRAY_1D,
205239
typename ARRAY_1D_TO_CONST,
206240
typename ARRAY_2D >
207241
HPCREACT_HOST_DEVICE inline void
208242
KineticReactions< REAL_TYPE,
209243
INT_TYPE,
210-
INDEX_TYPE
244+
INDEX_TYPE,
245+
LOGE_CONCENTRATION
211246
>::timeStep(
212247
RealType const dt,
213248
RealType const & temperature,

0 commit comments

Comments
 (0)