Skip to content

Commit 4fd890f

Browse files
rrsettgastrrsettgast
committed
break up equilbrium reactions to impl file:
1 parent 24c115a commit 4fd890f

File tree

4 files changed

+159
-46
lines changed

4 files changed

+159
-46
lines changed

‎src/reactions/bulkGeneric/EquilibriumReactions.hpp‎

Lines changed: 118 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,75 @@ class EquilibriumReactions
3636
ARRAY_1D_TO_CONST const & speciesConcentration0,
3737
ARRAY_1D_TO_CONST2 const & xi,
3838
ARRAY_1D & residual,
39-
ARRAY_2D & jacobian );
39+
ARRAY_2D & jacobian )
40+
{
41+
HPCREACT_UNUSED_VAR( temperature );
42+
constexpr int numSpecies = PARAMS_DATA::numSpecies;
43+
constexpr int numReactions = PARAMS_DATA::numReactions;
44+
45+
// initialize the species concentration
46+
RealType speciesConcentration[numSpecies];
47+
for( IndexType i=0; i<numSpecies; ++i )
48+
{
49+
speciesConcentration[i] = speciesConcentration0[i];
50+
for( IndexType r=0; r<numReactions; ++r )
51+
{
52+
speciesConcentration[i] += params.stoichiometricMatrix( r, i ) * xi[r];
53+
}
54+
}
55+
56+
// loop over reactions
57+
for( IndexType a=0; a<numReactions; ++a )
58+
{
59+
// get the equilibrium constant for this reaction
60+
RealType const Keq = params.equilibriumConstant( a );
61+
62+
RealType forwardProduct = 1.0;
63+
RealType reverseProduct = 1.0;
64+
RealType dForwardProduct_dxi[numReactions] = {0.0};
65+
RealType dReverseProduct_dxi[numReactions] = {0.0};
66+
// loop over species
67+
for( IndexType i=0; i<numSpecies; ++i )
68+
{
69+
RealType const s_ai = params.stoichiometricMatrix( a, i );
70+
71+
if( s_ai < 0.0 )
72+
{
73+
// forward reaction
74+
forwardProduct *= pow( speciesConcentration[i], -s_ai );
75+
// derivative of forward product with respect to xi
76+
for( IndexType b=0; b<numReactions; ++b )
77+
{
78+
dForwardProduct_dxi[b] += -s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
79+
}
80+
}
81+
else if( s_ai > 0.0 )
82+
{
83+
// reverse reaction
84+
reverseProduct *= pow( speciesConcentration[i], s_ai );
85+
// derivative of reverse product with respect to xi
86+
for( IndexType b=0; b<numReactions; ++b )
87+
{
88+
dReverseProduct_dxi[b] += s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
89+
}
90+
}
91+
}
92+
// compute the residual for this reaction
93+
residual[a] = forwardProduct - Keq * reverseProduct;
94+
// printf( "residual[%d] = %g - %g * %g = %g\n", a, forwardProduct, Keq, reverseProduct, residual[a] );
95+
96+
// Finish the derivatives of the product terms with respect to xi
97+
for( IndexType b=0; b<numReactions; ++b )
98+
{
99+
dForwardProduct_dxi[b] *= forwardProduct;
100+
dReverseProduct_dxi[b] *= reverseProduct;
101+
// printf( "(%d) dForwardProduct_dxi[%d] = %f\n", a, b, dForwardProduct_dxi[b] );
102+
// printf( "(%d) dReverseProduct_dxi[%d] = %f\n", a, b, dReverseProduct_dxi[b] );
103+
// printf( "Keq = %f\n", Keq );
104+
jacobian( a, b ) = dForwardProduct_dxi[b] - Keq * dReverseProduct_dxi[b];
105+
}
106+
}
107+
}
40108

41109

42110
template< typename PARAMS_DATA,
@@ -47,13 +115,58 @@ class EquilibriumReactions
47115
enforceEquilibrium( RealType const & temperature,
48116
PARAMS_DATA const & params,
49117
ARRAY_1D_TO_CONST const & speciesConcentration0,
50-
ARRAY_1D & speciesConcentration );
118+
ARRAY_1D & speciesConcentration )
119+
{
120+
HPCREACT_UNUSED_VAR( temperature );
121+
constexpr int numSpecies = PARAMS_DATA::numSpecies;
122+
constexpr int numReactions = PARAMS_DATA::numReactions;
123+
double residual[numReactions] = { 0.0 };
124+
double xi[numReactions] = { 0.0 };
125+
double dxi[numReactions] = { 0.0 };
126+
CArrayWrapper< double, numReactions, numReactions > jacobian;
127+
128+
REAL_TYPE residualNorm = 0.0;
129+
for( int k=0; k<10; ++k )
130+
{
131+
computeResidualAndJacobian( temperature,
132+
params,
133+
speciesConcentration0,
134+
xi,
135+
residual,
136+
jacobian );
137+
138+
residualNorm = 0.0;
139+
for( int j = 0; j < numReactions; ++j )
140+
{
141+
residualNorm += residual[j] * residual[j];
142+
}
143+
residualNorm = sqrt( residualNorm );
144+
if( residualNorm < 1.0e-14 )
145+
{
146+
break;
147+
}
148+
149+
solveNxN_pivoted< double, numReactions >( jacobian.data, residual, dxi );
150+
151+
// solve for the change in xi
152+
for( IndexType r=0; r<numReactions; ++r )
153+
{
154+
xi[r] -= dxi[r];
155+
}
156+
}
157+
158+
for( IndexType i=0; i<numSpecies; ++i )
159+
{
160+
speciesConcentration[i] = speciesConcentration0[i];
161+
for( IndexType r=0; r<numReactions; ++r )
162+
{
163+
speciesConcentration[i] += params.stoichiometricMatrix( r, i ) * xi[r];
164+
}
165+
}
166+
}
51167

52168
};
53169

54170

55171
} // namespace bulkGeneric
56172
} // namespace hpcReact
57-
58-
#include "EquilibriumReactions_impl.hpp"
59-
#include "common/macrosCleanup.hpp"

‎src/reactions/bulkGeneric/KineticReactions_impl.hpp‎

Lines changed: 23 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -327,29 +327,29 @@ KineticReactions< REAL_TYPE,
327327
}
328328

329329

330-
// printf( "residual = { " );
331-
// for( int i = 0; i < numSpecies; ++i )
332-
// {
333-
// printf( " %g, ", residual[i] );
334-
// }
335-
// printf( "}\n" );
336-
337-
// printf( "Jacobian = { \n" );
338-
// for( int i = 0; i < numSpecies; ++i )
339-
// {
340-
// printf( " { " );
341-
// for( int j = 0; j < numSpecies; ++j )
342-
// {
343-
// printf( " %g ", speciesRatesDerivatives( i, j ) );
344-
// // printf( " %g ", speciesRatesDerivatives( i, j ) / exp(speciesConcentration[j]) );
345-
// if( j < numSpecies-1 )
346-
// {
347-
// printf( ", " );
348-
// }
349-
// }
350-
// printf( "}, \n" );
351-
// }
352-
// printf( "}\n" );
330+
printf( "residual = { " );
331+
for( int i = 0; i < numSpecies; ++i )
332+
{
333+
printf( " %g, ", residual[i] );
334+
}
335+
printf( "}\n" );
336+
337+
printf( "Jacobian = { \n" );
338+
for( int i = 0; i < numSpecies; ++i )
339+
{
340+
printf( " { " );
341+
for( int j = 0; j < numSpecies; ++j )
342+
{
343+
printf( " %g ", speciesRatesDerivatives( i, j ) );
344+
// printf( " %g ", speciesRatesDerivatives( i, j ) / exp(speciesConcentration[j]) );
345+
if( j < numSpecies-1 )
346+
{
347+
printf( ", " );
348+
}
349+
}
350+
printf( "}, \n" );
351+
}
352+
printf( "}\n" );
353353

354354
solveNxN_pivoted< double, numSpecies >( speciesRatesDerivatives.data, residual, deltaPrimarySpeciesConcentration );
355355

‎src/reactions/bulkGeneric/unitTests/testKineticReactions.cpp‎

Lines changed: 17 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -264,7 +264,7 @@ void timeStepTest( PARAMS_DATA const & params,
264264

265265
TEST( testKineticReactions, testTimeStep )
266266
{
267-
double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
267+
double const initialSpeciesConcentration[5] = { 0.9, 1.0e-16, 0.5, 0.9, 1.0e-16 };
268268
double const expectedSpeciesConcentrations[5] = { 3.92138294e-01, 3.03930853e-01, 5.05945481e-01, 7.02014628e-01, 5.95970745e-01 };
269269
double const expectedSpeciesRates[5] = { -2.0, 1.0, 0.75, -0.25, 0.5 };
270270
double const expectedSpeciesRatesDerivatives[5][5] = { { -4.0, 1.0, 0.0, 0.0, 0.0 },
@@ -273,22 +273,22 @@ TEST( testKineticReactions, testTimeStep )
273273
{ 0.0, 0.0, -0.5, -0.25, 0.0 },
274274
{ 0.0, 0.0, 1.0, 0.5, 0.0 } };
275275

276-
timeStepTest< double, false >( simpleTestRateParams,
277-
2.0,
278-
10,
279-
initialSpeciesConcentration,
280-
expectedSpeciesConcentrations,
281-
expectedSpeciesRates,
282-
expectedSpeciesRatesDerivatives );
283-
284-
// ln(c) as the primary variable results in a singular system.
285-
// timeStepTest< double, true >( simpleTestRateParams,
286-
// 2.0,
287-
// 10,
288-
// initialSpeciesConcentration,
289-
// expectedSpeciesConcentrations,
290-
// expectedSpeciesRates,
291-
// expectedSpeciesRatesDerivatives );
276+
// timeStepTest< double, false >( simpleTestRateParams,
277+
// 2.0,
278+
// 10,
279+
// initialSpeciesConcentration,
280+
// expectedSpeciesConcentrations,
281+
// expectedSpeciesRates,
282+
// expectedSpeciesRatesDerivatives );
283+
284+
// ln(c) as the primary variable results in a singular system.
285+
timeStepTest< double, true >( simpleTestRateParams,
286+
2.0,
287+
10,
288+
initialSpeciesConcentration,
289+
expectedSpeciesConcentrations,
290+
expectedSpeciesRates,
291+
expectedSpeciesRatesDerivatives );
292292
}
293293

294294
int main( int argc, char * * argv )

‎src/uncrustify.cfg‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2147,4 +2147,4 @@ warn_level_tabs_found_in_verbatim_string_literals = 2 # unsigned number
21472147
## option(s) with 'not default' value: 46
21482148
#
21492149

2150-
set EXECUTION_CONTEXT __host__ __device__ GEOS_DEVICE GEOS_HPCREACT_HOST_DEVICE GEOS_HYPRE_DEVICE LVARRY_DEVICE LVARRY_HPCREACT_HOST_DEVICE HPCREACT_NO_MISSING_BRACES()
2150+
set EXECUTION_CONTEXT __host__ __device__ GEOS_DEVICE GEOS_HPCREACT_HOST_DEVICE GEOS_HYPRE_DEVICE LVARRY_DEVICE LVARRY_HPCREACT_HOST_DEVICE

0 commit comments

Comments
 (0)