Skip to content

Commit 7bc1f24

Browse files
rrsettgastrrsettgast
committed
added timestep integrator in kinetic reactions
1 parent a217b59 commit 7bc1f24

File tree

4 files changed

+222
-85
lines changed

4 files changed

+222
-85
lines changed

‎src/common/DirectSystemSolve.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1+
#pragma once
12
#include "macros.hpp"
2-
33
#include<cmath>
44

55
namespace hpcReact

‎src/reactions/bulkGeneric/KineticReactions.hpp‎

Lines changed: 30 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -28,13 +28,13 @@ class KineticReactions
2828
static HPCREACT_HOST_DEVICE inline void
2929
computeReactionRates( RealType const & temperature,
3030
PARAMS_DATA const & params,
31-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
31+
ARRAY_1D_TO_CONST const & speciesConcentration,
3232
ARRAY_1D & reactionRates,
3333
ARRAY_2D & reactionRatesDerivatives )
3434
{
3535
computeReactionRates_impl< PARAMS_DATA, true >( temperature,
3636
params,
37-
primarySpeciesConcentration,
37+
speciesConcentration,
3838
reactionRates,
3939
reactionRatesDerivatives );
4040
}
@@ -45,13 +45,13 @@ class KineticReactions
4545
static HPCREACT_HOST_DEVICE inline void
4646
computeReactionRates( RealType const & temperature,
4747
PARAMS_DATA const & params,
48-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
48+
ARRAY_1D_TO_CONST const & speciesConcentration,
4949
ARRAY_1D & reactionRates)
5050
{
5151
REAL_TYPE reactionRatesDerivatives[PARAMS_DATA::numReactions][PARAMS_DATA::numSpecies] = { 0.0 };
5252
computeReactionRates_impl< PARAMS_DATA, false >( temperature,
5353
params,
54-
primarySpeciesConcentration,
54+
speciesConcentration,
5555
reactionRates,
5656
reactionRatesDerivatives );
5757
}
@@ -64,15 +64,15 @@ class KineticReactions
6464
static HPCREACT_HOST_DEVICE inline void
6565
computeSpeciesRates( RealType const & temperature,
6666
PARAMS_DATA const & params,
67-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
67+
ARRAY_1D_TO_CONST const & speciesConcentration,
6868
ARRAY_1D & speciesRates,
69-
ARRAY_2D & primarySpeciesRatesDerivatives )
69+
ARRAY_2D & speciesRatesDerivatives )
7070
{
7171
computeSpeciesRates_impl< PARAMS_DATA, true >( temperature,
7272
params,
73-
primarySpeciesConcentration,
73+
speciesConcentration,
7474
speciesRates,
75-
primarySpeciesRatesDerivatives );
75+
speciesRatesDerivatives );
7676
}
7777

7878
template< typename PARAMS_DATA,
@@ -81,28 +81,41 @@ class KineticReactions
8181
static HPCREACT_HOST_DEVICE inline void
8282
computeSpeciesRates( RealType const & temperature,
8383
PARAMS_DATA const & params,
84-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
84+
ARRAY_1D_TO_CONST const & speciesConcentration,
8585
ARRAY_1D & speciesRates )
8686
{
87-
double primarySpeciesRatesDerivatives;
87+
double speciesRatesDerivatives;
8888
computeSpeciesRates_impl< PARAMS_DATA, false >( temperature,
8989
params,
90-
primarySpeciesConcentration,
90+
speciesConcentration,
9191
speciesRates,
92-
primarySpeciesRatesDerivatives );
92+
speciesRatesDerivatives );
9393
}
9494

95+
template< typename PARAMS_DATA,
96+
typename ARRAY_1D,
97+
typename ARRAY_1D_TO_CONST,
98+
typename ARRAY_2D >
99+
static HPCREACT_HOST_DEVICE void
100+
timeStep( RealType const dt,
101+
RealType const & temperature,
102+
PARAMS_DATA const & params,
103+
ARRAY_1D_TO_CONST const & speciesConcentration_n,
104+
ARRAY_1D & speciesConcentration,
105+
ARRAY_1D & speciesRates,
106+
ARRAY_2D & speciesRatesDerivatives );
107+
95108

96109
private:
97110
template< typename PARAMS_DATA,
98111
bool CALCULATE_DERIVATIVES,
99112
typename ARRAY_1D_TO_CONST,
100113
typename ARRAY_1D,
101114
typename ARRAY_2D >
102-
static HPCREACT_HOST_DEVICE inline void
115+
static HPCREACT_HOST_DEVICE void
103116
computeReactionRates_impl( RealType const & temperature,
104117
PARAMS_DATA const & params,
105-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
118+
ARRAY_1D_TO_CONST const & speciesConcentration,
106119
ARRAY_1D & reactionRates,
107120
ARRAY_2D & reactionRatesDerivatives );
108121

@@ -111,12 +124,12 @@ class KineticReactions
111124
typename ARRAY_1D_TO_CONST,
112125
typename ARRAY_1D,
113126
typename ARRAY_2D >
114-
static HPCREACT_HOST_DEVICE inline void
127+
static HPCREACT_HOST_DEVICE void
115128
computeSpeciesRates_impl( RealType const & temperature,
116129
PARAMS_DATA const & params,
117-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
130+
ARRAY_1D_TO_CONST const & speciesConcentration,
118131
ARRAY_1D & speciesRates,
119-
ARRAY_2D & primarySpeciesRatesDerivatives );
132+
ARRAY_2D & speciesRatesDerivatives );
120133

121134
};
122135

‎src/reactions/bulkGeneric/KineticReactions_impl.hpp‎

Lines changed: 83 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#include "common/constants.hpp"
33
#include "common/CArrayWrapper.hpp"
44
#include "common/macros.hpp"
5+
#include "common/DirectSystemSolve.hpp"
56

67
#include <cmath>
78
#include <string>
@@ -26,7 +27,7 @@ KineticReactions< REAL_TYPE,
2627
INDEX_TYPE
2728
>::computeReactionRates_impl( RealType const & ,//temperature,
2829
PARAMS_DATA const & params,
29-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
30+
ARRAY_1D_TO_CONST const & speciesConcentration,
3031
ARRAY_1D & reactionRates,
3132
ARRAY_2D & reactionRatesDerivatives )
3233
{
@@ -58,7 +59,7 @@ KineticReactions< REAL_TYPE,
5859
{
5960

6061
RealType const s_ri = params.stoichiometricMatrix( r, i );
61-
RealType const productTerm_i = primarySpeciesConcentration[i] > 1e-100 ? pow( primarySpeciesConcentration[i], abs(s_ri) ) : 0.0;
62+
RealType const productTerm_i = speciesConcentration[i] > 1e-100 ? pow( speciesConcentration[i], abs(s_ri) ) : 0.0;
6263

6364
if( s_ri < 0.0 )
6465
{
@@ -78,7 +79,7 @@ KineticReactions< REAL_TYPE,
7879
{
7980
if( i==j )
8081
{
81-
dProductConcForward_dC[j] *= -s_ri * pow( primarySpeciesConcentration[i], -s_ri-1 );
82+
dProductConcForward_dC[j] *= -s_ri * pow( speciesConcentration[i], -s_ri-1 );
8283
dProductConcReverse_dC[j] = 0.0;
8384
}
8485
else
@@ -93,7 +94,7 @@ KineticReactions< REAL_TYPE,
9394
{
9495
if( i==j )
9596
{
96-
dProductConcReverse_dC[j] *= s_ri * pow( primarySpeciesConcentration[i], s_ri-1 );
97+
dProductConcReverse_dC[j] *= s_ri * pow( speciesConcentration[i], s_ri-1 );
9798
dProductConcForward_dC[j] = 0.0;
9899
}
99100
else
@@ -150,45 +151,114 @@ KineticReactions< REAL_TYPE,
150151
INDEX_TYPE
151152
>::computeSpeciesRates_impl( RealType const & temperature,
152153
PARAMS_DATA const & params,
153-
ARRAY_1D_TO_CONST const & primarySpeciesConcentration,
154-
ARRAY_1D & primarySpeciesRates,
155-
ARRAY_2D & primarySpeciesRatesDerivatives )
154+
ARRAY_1D_TO_CONST const & speciesConcentration,
155+
ARRAY_1D & speciesRates,
156+
ARRAY_2D & speciesRatesDerivatives )
156157
{
157158
RealType reactionRates[PARAMS_DATA::numReactions] = { 0.0 };
158159
CArrayWrapper< double, PARAMS_DATA::numReactions, PARAMS_DATA::numSpecies > reactionRatesDerivatives;
159160

160161
if constexpr ( !CALCULATE_DERIVATIVES )
161162
{
162-
HPCREACT_UNUSED_VAR(primarySpeciesRatesDerivatives);
163+
HPCREACT_UNUSED_VAR(speciesRatesDerivatives);
163164
}
164165

165-
computeReactionRates< PARAMS_DATA >( temperature, params, primarySpeciesConcentration, reactionRates, reactionRatesDerivatives );
166+
computeReactionRates< PARAMS_DATA >( temperature, params, speciesConcentration, reactionRates, reactionRatesDerivatives );
166167

167168
for( IntType i = 0; i < PARAMS_DATA::numSpecies; ++i )
168169
{
169-
primarySpeciesRates[i] = 0.0;
170+
speciesRates[i] = 0.0;
170171
if constexpr ( CALCULATE_DERIVATIVES )
171172
{
172173
for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j )
173174
{
174-
primarySpeciesRatesDerivatives( i, j ) = 0.0;
175+
speciesRatesDerivatives( i, j ) = 0.0;
175176
}
176177
}
177178
for( IntType r=0; r<PARAMS_DATA::numReactions; ++r )
178179
{
179180
RealType const s_ir = params.stoichiometricMatrix( r, i );
180-
primarySpeciesRates[i] += s_ir * reactionRates[r];
181+
speciesRates[i] += s_ir * reactionRates[r];
181182
if constexpr ( CALCULATE_DERIVATIVES )
182183
{
183184
for( IntType j = 0; j < PARAMS_DATA::numSpecies; ++j )
184185
{
185-
primarySpeciesRatesDerivatives( i, j ) += s_ir * reactionRatesDerivatives( r, j );
186+
speciesRatesDerivatives( i, j ) += s_ir * reactionRatesDerivatives( r, j );
186187
}
187188
}
188189
}
189190
}
190191
}
191192

193+
template< typename REAL_TYPE,
194+
typename INT_TYPE,
195+
typename INDEX_TYPE >
196+
template< typename PARAMS_DATA,
197+
typename ARRAY_1D,
198+
typename ARRAY_1D_TO_CONST,
199+
typename ARRAY_2D >
200+
HPCREACT_HOST_DEVICE inline void
201+
KineticReactions< REAL_TYPE,
202+
INT_TYPE,
203+
INDEX_TYPE
204+
>::timeStep(
205+
RealType const dt,
206+
RealType const & temperature,
207+
PARAMS_DATA const & params,
208+
ARRAY_1D_TO_CONST const & speciesConcentration_n,
209+
ARRAY_1D & speciesConcentration,
210+
ARRAY_1D & speciesRates,
211+
ARRAY_2D & speciesRatesDerivatives )
212+
{
213+
// constexpr int numReactions = PARAMS_DATA::numReactions;
214+
constexpr int numSpecies = PARAMS_DATA::numSpecies;
215+
216+
217+
218+
219+
REAL_TYPE residualNorm = 0.0;
220+
for( int k=0; k<10; ++k )
221+
{
222+
// printf( "iteration %2d: \n", k );
223+
224+
computeSpeciesRates( temperature,
225+
params,
226+
speciesConcentration,
227+
speciesRates,
228+
speciesRatesDerivatives );
229+
230+
double residual[numSpecies] = { 0.0 };
231+
double deltaPrimarySpeciesConcentration[numSpecies] = { 0.0 };
232+
for( int i = 0; i < numSpecies; ++i )
233+
{
234+
residual[i] = -(speciesConcentration[i] - speciesConcentration_n[i] - dt * speciesRates[i]);
235+
for( int j = 0; j < numSpecies; ++j )
236+
{
237+
speciesRatesDerivatives( i, j ) = - dt * speciesRatesDerivatives( i, j );
238+
}
239+
speciesRatesDerivatives( i, i ) += 1.0;
240+
}
241+
242+
residualNorm = 0.0;
243+
for( int j = 0; j < numSpecies; ++j )
244+
{
245+
residualNorm += residual[j] * residual[j];
246+
}
247+
residualNorm = sqrt( residualNorm );
248+
if( residualNorm < 1.0e-8 )
249+
{
250+
break;
251+
}
252+
253+
solveNxN_pivoted<double,numSpecies>( speciesRatesDerivatives.data, residual, deltaPrimarySpeciesConcentration );
254+
255+
for( int i = 0; i < numSpecies; ++i )
256+
{
257+
speciesConcentration[i] += deltaPrimarySpeciesConcentration[i];
258+
}
259+
260+
}
261+
}
192262
} // namespace bulkGeneric
193263
} // namespace hpcReact
194264

0 commit comments

Comments
 (0)