Skip to content

Commit 8dbbaac

Browse files
CusiniMCusiniM
committed
wip: add mixedReactions unit test.
1 parent 462beae commit 8dbbaac

File tree

6 files changed

+336
-42
lines changed

6 files changed

+336
-42
lines changed

‎src/common/nonlinearSolvers.hpp‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
#pragma once
2+
3+
#include "macros.hpp"
4+
#include "DirectSystemSolve.hpp"
5+
#include <cmath>
6+
7+
namespace hpcReact
8+
{
9+
10+
namespace nonlinearSolvers
11+
{
12+
13+
namespace internal
14+
{
15+
16+
template< int N >
17+
HPCREACT_HOST_DEVICE
18+
double norm( double const (&r)[N])
19+
{
20+
double sum = 0.0;
21+
for ( int i = 0; i < N; ++i ) sum += r[i] * r[i];
22+
return ::sqrt( sum );
23+
}
24+
25+
26+
template< int N >
27+
HPCREACT_HOST_DEVICE
28+
void add(double (&x)[N], double const (&dx)[N])
29+
{
30+
for (int i = 0; i < N; ++i) x[i] += dx[i];
31+
}
32+
33+
template< int N >
34+
HPCREACT_HOST_DEVICE
35+
void scale( double (&x)[N], double const value )
36+
{
37+
for (int i = 0; i < N; ++i) x[i] *= value;
38+
}
39+
40+
}
41+
42+
template< int N,
43+
typename REAL_TYPE,
44+
typename ResidualFunc,
45+
typename JacobianFunc >
46+
HPCREACT_HOST_DEVICE
47+
void newtonRaphson( REAL_TYPE (&x)[N],
48+
ResidualFunc computeResidual,
49+
JacobianFunc computeJacobian,
50+
int maxIters = 25,
51+
double tol = 1e-10 )
52+
{
53+
REAL_TYPE residual[N];
54+
REAL_TYPE dx[N];
55+
REAL_TYPE jacobian[N][N];
56+
57+
for ( int iter = 0; iter < maxIters; ++iter )
58+
{
59+
computeResidual( x, residual );
60+
61+
if ( internal::norm< N >( residual ) < tol ) return;
62+
63+
computeJacobian( x, jacobian );
64+
65+
solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx );
66+
67+
internal::add< N >( x, dx );
68+
}
69+
}
70+
71+
template< int N,
72+
typename REAL_TYPE,
73+
typename FUNCTION_TYPE >
74+
HPCREACT_HOST_DEVICE
75+
bool newtonRaphson( REAL_TYPE (&x)[N],
76+
FUNCTION_TYPE computeResidualAndJacobian,
77+
int maxIters = 25,
78+
double tol = 1e-10 )
79+
{
80+
REAL_TYPE residual[N];
81+
REAL_TYPE dx[N];
82+
REAL_TYPE jacobian[N][N];
83+
84+
for ( int iter = 0; iter < maxIters; ++iter )
85+
{
86+
computeResidualAndJacobian( x, residual, jacobian );
87+
88+
double const norm = internal::norm< N >( residual );
89+
90+
printf( "--Iter %d: Residual norm = %.12e\n", iter, norm );
91+
92+
if ( norm < tol ) {
93+
printf( "--Converged.\n" );
94+
return true;
95+
}
96+
internal::scale<N>( residual, -1.0);
97+
solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx );
98+
internal::add< N >( x, dx );
99+
}
100+
101+
printf( "--Newton solver error: Max iterations reached without convergence.\n" );
102+
103+
return false;
104+
}
105+
106+
}
107+
}

‎src/reactions/bulkGeneric/MixedEquilibriumKineticReactions.hpp‎

Lines changed: 24 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -44,19 +44,22 @@ class MixedEquilibriumKineticReactions
4444

4545
template< typename PARAMS_DATA,
4646
typename ARRAY_1D_TO_CONST,
47-
typename ARRAY_1D,
48-
typename ARRAY_2D >
47+
typename ARRAY_1D_PRIMARY,
48+
typename ARRAY_1D_SECONDARY,
49+
typename ARRAY_1D_KINETIC,
50+
typename ARRAY_2D_PRIMARY,
51+
typename ARRAY_2D_KINETIC >
4952
static HPCREACT_HOST_DEVICE inline void
5053
updateMixedSystem( RealType const & temperature,
5154
PARAMS_DATA const & params,
5255
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
53-
ARRAY_1D & logSecondarySpeciesConcentrations,
54-
ARRAY_1D & aggregatePrimarySpeciesConcentrations,
55-
ARRAY_2D & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
56-
ARRAY_1D & reactionRates,
57-
ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations,
58-
ARRAY_1D & aggregateSpeciesRates,
59-
ARRAY_2D & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations )
56+
ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations,
57+
ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations,
58+
ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
59+
ARRAY_1D_KINETIC & reactionRates,
60+
ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations,
61+
ARRAY_1D_PRIMARY & aggregateSpeciesRates,
62+
ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations )
6063
{
6164
updateMixedSystem_impl( temperature,
6265
params,
@@ -124,19 +127,22 @@ class MixedEquilibriumKineticReactions
124127

125128
template< typename PARAMS_DATA,
126129
typename ARRAY_1D_TO_CONST,
127-
typename ARRAY_1D,
128-
typename ARRAY_2D >
130+
typename ARRAY_1D_PRIMARY,
131+
typename ARRAY_1D_SECONDARY,
132+
typename ARRAY_1D_KINETIC,
133+
typename ARRAY_2D_PRIMARY,
134+
typename ARRAY_2D_KINETIC >
129135
static HPCREACT_HOST_DEVICE void
130136
updateMixedSystem_impl( RealType const & temperature,
131137
PARAMS_DATA const & params,
132138
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
133-
ARRAY_1D & logSecondarySpeciesConcentrations,
134-
ARRAY_1D & aggregatePrimarySpeciesConcentrations,
135-
ARRAY_2D & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
136-
ARRAY_1D & reactionRates,
137-
ARRAY_2D & dReactionRates_dLogPrimarySpeciesConcentrations,
138-
ARRAY_1D & aggregateSpeciesRates,
139-
ARRAY_2D & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations );
139+
ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations,
140+
ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations,
141+
ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
142+
ARRAY_1D_KINETIC & reactionRates,
143+
ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations,
144+
ARRAY_1D_PRIMARY & aggregateSpeciesRates,
145+
ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations );
140146

141147
template< typename PARAMS_DATA,
142148
typename ARRAY_1D_TO_CONST,

‎src/reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp‎

Lines changed: 20 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -21,33 +21,33 @@ template< typename REAL_TYPE,
2121
typename INDEX_TYPE,
2222
bool LOGE_CONCENTRATION >
2323
template< typename PARAMS_DATA,
24-
typename ARRAY_1D_TO_CONST,
25-
typename ARRAY_1D,
26-
typename ARRAY_2D >
24+
typename ARRAY_1D_TO_CONST,
25+
typename ARRAY_1D_PRIMARY,
26+
typename ARRAY_1D_SECONDARY,
27+
typename ARRAY_1D_KINETIC,
28+
typename ARRAY_2D_PRIMARY,
29+
typename ARRAY_2D_KINETIC >
2730
HPCREACT_HOST_DEVICE inline void
2831
MixedEquilibriumKineticReactions< REAL_TYPE,
2932
INT_TYPE,
3033
INDEX_TYPE,
3134
LOGE_CONCENTRATION
3235
>::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 )
36+
PARAMS_DATA const & params,
37+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
38+
ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations,
39+
ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations,
40+
ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
41+
ARRAY_1D_KINETIC & reactionRates,
42+
ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations,
43+
ARRAY_1D_PRIMARY & aggregateSpeciesRates,
44+
ARRAY_2D_PRIMARY & dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations )
4245
{
43-
// constexpr IntType numSpecies = PARAMS_DATA::numSpecies();
44-
// constexpr IntType numSecondarySpecies = PARAMS_DATA::numReactions();
45-
// constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
4646

4747
// 1. Compute new aggregate species from primary species
4848
calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE,
49-
INT_TYPE,
50-
INDEX_TYPE >( params,
49+
INT_TYPE,
50+
INDEX_TYPE >( params.equilibriumReactionsParameters(),
5151
logPrimarySpeciesConcentrations,
5252
logSecondarySpeciesConcentrations,
5353
aggregatePrimarySpeciesConcentrations,
@@ -105,9 +105,9 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
105105
{
106106
logSpeciesConcentration[i] = logSecondarySpeciesConcentrations[i];
107107
}
108-
for ( IntType i = 0; i < numPrimarySpecies; ++i )
108+
for ( IntType i = numSecondarySpecies; i < numSpecies; ++i )
109109
{
110-
logSpeciesConcentration[i+numSecondarySpecies] = logPrimarySpeciesConcentrations[i];
110+
logSpeciesConcentration[i] = logPrimarySpeciesConcentrations[i];
111111
}
112112

113113
CArrayWrapper< RealType, numKineticReactions, numSpecies > reactionRatesDerivatives;
@@ -159,7 +159,7 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
159159
ARRAY_1D & aggregatesRates,
160160
ARRAY_2D & aggregatesRatesDerivatives )
161161
{
162-
GEOS_UNUSED_VAR( speciesConcentration );
162+
HPCREACT_UNUSED_VAR( speciesConcentration );
163163
// constexpr IntType numSpecies = PARAMS_DATA::numSpecies();
164164
constexpr IntType numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
165165
constexpr IntType numKineticReactions = PARAMS_DATA::numKineticReactions();

‎src/reactions/bulkGeneric/SpeciesUtilities.hpp‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -94,14 +94,15 @@ template< typename REAL_TYPE,
9494
typename INDEX_TYPE,
9595
typename PARAMS_DATA,
9696
typename ARRAY_1D_TO_CONST,
97-
typename ARRAY_1D,
97+
typename ARRAY_1D_PRIMARY,
98+
typename ARRAY_1D_SECONDARY,
9899
typename ARRAY_2D >
99100
HPCREACT_HOST_DEVICE
100101
inline
101-
void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA & params,
102+
void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
102103
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
103-
ARRAY_1D & logSecondarySpeciesConcentrations,
104-
ARRAY_1D & aggregatePrimarySpeciesConcentrations,
104+
ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations,
105+
ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations,
105106
ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations )
106107
{
107108
constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();

‎src/reactions/bulkGeneric/unitTests/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ set( testSourceFiles
33
testEquilibriumReactions.cpp
44
testKineticReactions.cpp
55
testSpeciesUtilities.cpp
6+
testMixedReactions.cpp
67
)
78

89
set( dependencyList hpcReact gtest )

0 commit comments

Comments
 (0)