Skip to content

Commit d6cfe4f

Browse files
rrsettgastrrsettgast
committed
../src/common/printers.hpp
1 parent a32b9b8 commit d6cfe4f

File tree

3 files changed

+15
-179
lines changed

3 files changed

+15
-179
lines changed

‎src/reactions/bulkGeneric/EquilibriumReactions.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,13 +13,14 @@ namespace bulkGeneric
1313
template< typename REAL_TYPE,
1414
typename INT_TYPE,
1515
typename INDEX_TYPE,
16-
int RESIDUAL_FORM >
16+
int FORMULATION >
1717
class EquilibriumReactions
1818
{
1919
public:
2020
using RealType = REAL_TYPE;
2121
using IntType = INT_TYPE;
2222
using IndexType = INDEX_TYPE;
23+
// constexpr static int formulation = FORMULATION;
2324

2425
template< typename PARAMS_DATA,
2526
typename ARRAY_1D,

‎src/reactions/bulkGeneric/EquilibriumReactions_impl.hpp‎

Lines changed: 13 additions & 145 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ constexpr bool debugPrinting = false;
1313
template< typename REAL_TYPE,
1414
typename INT_TYPE,
1515
typename INDEX_TYPE,
16-
int RESIDUAL_FORM >
16+
int FORMULATION >
1717
template< typename PARAMS_DATA,
1818
typename ARRAY_1D,
1919
typename ARRAY_1D_TO_CONST,
@@ -24,23 +24,14 @@ void
2424
EquilibriumReactions< REAL_TYPE,
2525
INT_TYPE,
2626
INDEX_TYPE,
27-
RESIDUAL_FORM >::computeResidualAndJacobian( RealType const & temperature,
27+
FORMULATION >::computeResidualAndJacobian( REAL_TYPE const & temperature,
2828
PARAMS_DATA const & params,
2929
ARRAY_1D_TO_CONST const & speciesConcentration0,
3030
ARRAY_1D_TO_CONST2 const & xi,
3131
ARRAY_1D & residual,
3232
ARRAY_2D & jacobian )
3333
{
3434

35-
// printf( "speciesConcentration0 = | %g %g %g %g %g |\n",
36-
// speciesConcentration0[0],
37-
// speciesConcentration0[1],
38-
// speciesConcentration0[2],
39-
// speciesConcentration0[3],
40-
// speciesConcentration0[4] );
41-
// printf( "xi = | %g %g |\n", xi[0], xi[1] );
42-
43-
4435
HPCREACT_UNUSED_VAR( temperature );
4536
constexpr int numSpecies = PARAMS_DATA::numSpecies;
4637
constexpr int numReactions = PARAMS_DATA::numReactions;
@@ -53,10 +44,6 @@ EquilibriumReactions< REAL_TYPE,
5344
for( IndexType r=0; r<numReactions; ++r )
5445
{
5546
speciesConcentration[i] += params.stoichiometricMatrix( r, i ) * xi[r];
56-
// if( speciesConcentration[i] < 1.0e-17 )
57-
// {
58-
// speciesConcentration[i] = 1.0e-17;
59-
// }
6047
}
6148
}
6249

@@ -73,8 +60,8 @@ EquilibriumReactions< REAL_TYPE,
7360

7461
// these actually only hold the derivatives of the single concentration term for the product...not the full product.
7562
// it will have to be multiplied by the product itself to get the derivative of the product.
76-
RealType dForwardProduct_dxi[numReactions] = {0.0};
77-
RealType dReverseProduct_dxi[numReactions] = {0.0};
63+
RealType dForwardProduct_dxi_divProduct[numReactions] = {0.0};
64+
RealType dReverseProduct_dxi_divProduct[numReactions] = {0.0};
7865
// loop over species
7966
for( IndexType i=0; i<numSpecies; ++i )
8067
{
@@ -86,7 +73,7 @@ EquilibriumReactions< REAL_TYPE,
8673
// derivative of forward product with respect to xi
8774
for( IndexType b=0; b<numReactions; ++b )
8875
{
89-
dForwardProduct_dxi[b] += -s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
76+
dForwardProduct_dxi_divProduct[b] += -s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
9077
}
9178
}
9279
else if( s_ai > 0.0 )
@@ -96,60 +83,25 @@ EquilibriumReactions< REAL_TYPE,
9683
// derivative of reverse product with respect to xi
9784
for( IndexType b=0; b<numReactions; ++b )
9885
{
99-
dReverseProduct_dxi[b] += s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
86+
dReverseProduct_dxi_divProduct[b] += s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
10087
}
10188
}
10289
}
10390
// compute the residual for this reaction
104-
if constexpr( RESIDUAL_FORM == 0 )
105-
{
106-
residual[a] = Keq * forwardProduct - reverseProduct;
107-
}
108-
else if constexpr( RESIDUAL_FORM == 1 )
109-
{
110-
residual[a] = 1.0 - reverseProduct / ( forwardProduct * Keq );
111-
}
112-
else if constexpr( RESIDUAL_FORM == 2 )
113-
{
114-
residual[a] = log( reverseProduct / ( forwardProduct * Keq ) );
115-
}
116-
// printf( "residual[%d] = %g - %g * %g = %g\n", a, forwardProduct, Keq, reverseProduct, residual[a] );
91+
residual[a] = log( reverseProduct / ( forwardProduct * Keq ) );
11792

118-
// Finish the derivatives of the product terms with respect to xi
93+
// compute the jacobian
11994
for( IndexType b=0; b<numReactions; ++b )
12095
{
121-
// printf( "(%d) dForwardProduct_dxi[%d] = %f\n", a, b, dForwardProduct_dxi[b] );
122-
// printf( "(%d) dReverseProduct_dxi[%d] = %f\n", a, b, dReverseProduct_dxi[b] );
123-
// printf( "Keq = %f\n", Keq );
124-
125-
126-
if constexpr( RESIDUAL_FORM == 0 )
127-
{
128-
dForwardProduct_dxi[b] *= forwardProduct;
129-
dReverseProduct_dxi[b] *= reverseProduct;
130-
jacobian( a, b ) = Keq * dForwardProduct_dxi[b] - dReverseProduct_dxi[b];
131-
}
132-
else if constexpr( RESIDUAL_FORM == 1 )
133-
{
134-
dForwardProduct_dxi[b] *= forwardProduct;
135-
dReverseProduct_dxi[b] *= reverseProduct;
136-
jacobian( a, b ) = -1/Keq * ( dReverseProduct_dxi[b] / forwardProduct - dForwardProduct_dxi[b] * reverseProduct / ( forwardProduct * forwardProduct ) );
137-
}
138-
else if constexpr( RESIDUAL_FORM == 2 )
139-
{
140-
// printf( "(%d) dForwardProduct_dxi[%d] = %f\n", a, b, dForwardProduct_dxi[b] );
141-
// printf( "(%d) dReverseProduct_dxi[%d] = %f\n", a, b, dReverseProduct_dxi[b] );
142-
jacobian( a, b ) = -dForwardProduct_dxi[b] + dReverseProduct_dxi[b];
143-
// printf( " jacobian( a, b ) = %g + %g = %g\n", -dForwardProduct_dxi[b], dReverseProduct_dxi[b], jacobian( a, b ) );
144-
}
96+
jacobian( a, b ) = -dForwardProduct_dxi_divProduct[b] + dReverseProduct_dxi_divProduct[b];
14597
}
14698
}
14799
}
148100

149101
template< typename REAL_TYPE,
150102
typename INT_TYPE,
151103
typename INDEX_TYPE,
152-
int RESIDUAL_FORM >
104+
int FORMULATION >
153105
template< typename PARAMS_DATA,
154106
typename ARRAY_1D,
155107
typename ARRAY_1D_TO_CONST >
@@ -158,7 +110,7 @@ void
158110
EquilibriumReactions< REAL_TYPE,
159111
INT_TYPE,
160112
INDEX_TYPE,
161-
RESIDUAL_FORM >::enforceEquilibrium( RealType const & temperature,
113+
FORMULATION >::computeResidualAndJacobianReactionExtents( REAL_TYPE const & temperature,
162114
PARAMS_DATA const & params,
163115
ARRAY_1D_TO_CONST const & speciesConcentration0,
164116
ARRAY_1D & speciesConcentration )
@@ -194,78 +146,15 @@ EquilibriumReactions< REAL_TYPE,
194146
break;
195147
}
196148

197-
if constexpr( debugPrinting )
198-
{
199-
printf( "************************************************** \n" );
200-
printf( "R = { " );
201-
for( int i=0; i<numReactions; ++i )
202-
{
203-
printf( "%18.12g", residual[i] );
204-
if( i < numReactions-1 )
205-
{
206-
printf( ", " );
207-
}
208-
}
209-
printf( " }\n" );
210-
211-
printf( "J = { \n" );
212-
for( int i=0; i<numReactions; ++i )
213-
{
214-
printf( " { " );
215-
for( int j=0; j<numReactions; ++j )
216-
{
217-
printf( "%22.14g", jacobian( i, j ) );
218-
if( j < numReactions-1 )
219-
{
220-
printf( ", " );
221-
}
222-
}
223-
printf( "},\n" );
224-
}
225-
printf( " }\n" );
226-
227-
int numNonSymmetric = 0;
228-
for( int i=0; i<numReactions; ++i )
229-
{
230-
for( int j=i; j<numReactions; ++j )
231-
{
232-
if( fabs( jacobian( i, j ) - jacobian( j, i ) ) > 0.5 * ( jacobian( i, j ) + jacobian( j, i ) )*1.0e-14 )
233-
{
234-
++numNonSymmetric;
235-
printf( "jacobian not symmetric: i = %d, j = %d, jacobian(i,j) = %g, jacobian(j,i) = %g\n", i, j, jacobian( i, j ), jacobian( j, i ) );
236-
}
237-
}
238-
}
239-
if( numNonSymmetric > 0 )
240-
{
241-
printf( "numNonSymmetric = %d\n", numNonSymmetric );
242-
}
243-
244-
printf( " is J PD = %d\n", isPositiveDefinite< double, numReactions >( jacobian.data ) );
245-
246-
}
247-
248-
249149
// solve for the change in xi
250150
for( int r=0; r<numReactions; ++r )
251151
{
252152
dxi[r] = 0.0;
253153
residual[r] = -residual[r];
254154
}
255-
if constexpr( RESIDUAL_FORM == 2 )
256-
{
257-
solveNxN_Cholesky< double, numReactions >( jacobian.data, residual, dxi );
258-
}
259-
else
260-
{
261-
solveNxN_pivoted< double, numReactions >( jacobian.data, residual, dxi );
262-
}
263155

264-
if constexpr( debugPrinting )
265-
{
266-
print( xi, "xi", 12 );
267-
print( dxi, "dxi", 12 );
268-
}
156+
solveNxN_Cholesky< double, numReactions >( jacobian.data, residual, dxi );
157+
269158

270159
// scaling
271160
REAL_TYPE scale = 1.0;
@@ -294,27 +183,6 @@ EquilibriumReactions< REAL_TYPE,
294183
xi[r] += scale * dxi[r];
295184
}
296185

297-
298-
if constexpr( debugPrinting )
299-
{
300-
printf( "c = { " );
301-
for( IndexType i=0; i<numSpecies; ++i )
302-
{
303-
REAL_TYPE c = speciesConcentration0[i];
304-
for( IndexType r=0; r<numReactions; ++r )
305-
{
306-
c += params.stoichiometricMatrix( r, i ) * xi[r];
307-
}
308-
printf( "%18.12g", c );
309-
if( i < numSpecies-1 )
310-
{
311-
printf( ", " );
312-
}
313-
}
314-
printf( " }\n" );
315-
}
316-
317-
318186
}
319187

320188
for( IndexType i=0; i<numSpecies; ++i )

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

Lines changed: 0 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -79,29 +79,6 @@ TEST( testEquilibriumReactions, computeResidualAndJacobianTest )
7979
{
8080
double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
8181

82-
std::cout<<" RESIDUAL_FORM 0:"<<std::endl;
83-
{
84-
double const expectedResiduals[] = { 1.0, 0.5 };
85-
double const expectedJacobian[2][2] =
86-
{ { -4.5, 1.0e-16 },
87-
{ 1.0, -1.5 } };
88-
computeResidualAndJacobianTest< double, 0 >( simpleTestRateParams,
89-
initialSpeciesConcentration,
90-
expectedResiduals,
91-
expectedJacobian );
92-
}
93-
94-
{
95-
std::cout<<" RESIDUAL_FORM 1:"<<std::endl;
96-
double const expectedResiduals[] = { 1.0, 1.0 };
97-
double const expectedJacobian[2][2] =
98-
{ { -0.5,1.0e-16 },
99-
{ 4.0e-32,-8.0e-16 } };
100-
computeResidualAndJacobianTest< double, 1 >( simpleTestRateParams,
101-
initialSpeciesConcentration,
102-
expectedResiduals,
103-
expectedJacobian );
104-
}
10582

10683
{
10784
std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
@@ -163,16 +140,6 @@ TEST( testEquilibriumReactions, testEnforceEquilibrium )
163140
double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
164141
double const expectedSpeciesConcentrations[5] = { 3.92138294e-01, 3.03930853e-01, 5.05945481e-01, 7.02014628e-01, 5.95970745e-01 };
165142

166-
std::cout<<" RESIDUAL_FORM 0:"<<std::endl;
167-
testEnforceEquilibrium< double, 0 >( simpleTestRateParams.equilibriumReactionsParameters(),
168-
initialSpeciesConcentration,
169-
expectedSpeciesConcentrations );
170-
171-
// std::cout<<std::endl;
172-
// std::cout<<" RESIDUAL_FORM 1:"<<std::endl;
173-
// testEnforceEquilibrium< double, 1 >( simpleTestRateParams.equilibriumReactionsParameters(),
174-
// initialSpeciesConcentration,
175-
// expectedSpeciesConcentrations );
176143

177144
std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
178145
testEnforceEquilibrium< double, 2 >( simpleTestRateParams.equilibriumReactionsParameters(),

0 commit comments

Comments
 (0)