Skip to content

Commit 1667909

Browse files
rrsettgastrrsettgast
committed
some tweaks
1 parent 50587d6 commit 1667909

File tree

5 files changed

+155
-130
lines changed

5 files changed

+155
-130
lines changed

‎src/reactions/bulkGeneric/EquilibriumReactions_impl.hpp‎

Lines changed: 25 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1+
#include "common/printers.hpp"
2+
13
#include <iostream>
24

35

@@ -245,20 +247,24 @@ EquilibriumReactions< REAL_TYPE,
245247

246248

247249
// solve for the change in xi
248-
solveNxN_Cholesky< double, numReactions >( jacobian.data, residual, dxi );
250+
for( int r=0; r<numReactions; ++r )
251+
{
252+
dxi[r] = 0.0;
253+
residual[r] = -residual[r];
254+
}
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+
}
249263

250264
if constexpr( debugPrinting )
251265
{
252-
printf( "dxi = { " );
253-
for( int i=0; i<numReactions; ++i )
254-
{
255-
printf( "%18.12g", dxi[i] );
256-
if( i < numReactions-1 )
257-
{
258-
printf( ", " );
259-
}
260-
}
261-
printf( " }\n" );
266+
print( xi, "xi", 12 );
267+
print( dxi, "dxi", 12 );
262268
}
263269

264270
// scaling
@@ -271,25 +277,21 @@ EquilibriumReactions< REAL_TYPE,
271277
{
272278
dc += params.stoichiometricMatrix( r, i ) * dxi[r];
273279
cn += params.stoichiometricMatrix( r, i ) * xi[r];
274-
if( cn-dc < 1.0e-30 )
280+
}
281+
if( cn+dc < 1.0e-30 )
282+
{
283+
REAL_TYPE const fscale = ( 1.0e-30 - cn ) / (dc);
284+
if( fscale < scale )
275285
{
276-
REAL_TYPE const fscale = ( 1.0e-30 - cn ) / (-dc);
277-
if( fscale < scale )
278-
{
279-
scale = 0.99*fscale;
280-
}
286+
scale = 0.9*fscale;
287+
//printf( "i, cn, dc, scale = %d, %18.12g %18.12g %18.12g\n", i, cn, dc, scale );
281288
}
282289
}
283290
}
284291

285-
if constexpr( debugPrinting )
286-
{
287-
printf( "scale = %18.12g\n", scale );
288-
}
289-
290292
for( IndexType r=0; r<numReactions; ++r )
291293
{
292-
xi[r] -= scale * dxi[r];
294+
xi[r] += scale * dxi[r];
293295
}
294296

295297

‎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/ParametersPredefined.hpp‎

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -83,11 +83,22 @@ carbonateSystem =
8383
3.72E-03, // MgSO4 = Mg+2 + SO4-2
8484
1.51E-01 }, // NaSO4- = Na+ + SO4-2
8585
// forward rate constants
86-
{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
86+
{ 9.77E+13,
87+
4.37E-07,
88+
2.14E+10,
89+
1.70E-04,
90+
8.13E-02,
91+
1.17E+07,
92+
6.92E-03,
93+
4.68E+00,
94+
3.98E+00,
95+
3.72E-03,
96+
1.51E-01 },
8797
// reverse rate constants
88-
{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
98+
{ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 }
8999
};
90100

101+
91102
// *****UNCRUSTIFY-ON******
92103
} // namespace bulkGeneric
93104
} // namespace hpcReact

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

Lines changed: 76 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -54,16 +54,16 @@ void computeResidualAndJacobianTest( PARAMS_DATA const & params,
5454
residual,
5555
jacobian );
5656

57-
printf( "R = | %8.4g %8.4g |\n", residual[0], residual[1] );
57+
// printf( "R = { %8.4g, %8.4g }\n", residual[0], residual[1] );
5858
for( int r=0; r<numReactions; ++r )
5959
{
6060
EXPECT_NEAR( residual[r], expectedResidual[r], 1.0e-8 );
6161
}
6262

63-
HPCREACT_UNUSED_VAR( expectedJacobian );
64-
printf( "J = | %8.4g %8.4g |\n | %8.4g %8.4g |\n\n",
65-
jacobian( 0, 0 ), jacobian( 0, 1 ),
66-
jacobian( 1, 0 ), jacobian( 1, 1 ) );
63+
// HPCREACT_UNUSED_VAR( expectedJacobian );
64+
// printf( "J = { { %8.4g, %8.4g },\n { %8.4g %8.4g } }\n\n",
65+
// jacobian( 0, 0 ), jacobian( 0, 1 ),
66+
// jacobian( 1, 0 ), jacobian( 1, 1 ) );
6767
for( int a = 0; a < numReactions; ++a )
6868
{
6969
for( int b = 0; b < numReactions; ++b )
@@ -75,47 +75,47 @@ void computeResidualAndJacobianTest( PARAMS_DATA const & params,
7575

7676

7777
//******************************************************************************
78-
// TEST( testEquilibriumReactions, computeResidualAndJacobianTest )
79-
// {
80-
// double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
81-
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-
// }
105-
106-
// {
107-
// std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
108-
// double const expectedResiduals[] = { -37.534508668465, -72.989575795250 };
109-
// double const expectedJacobian[2][2] =
110-
// { { 1.0e16,-2.0 },
111-
// { -2.0,4.0e16 } };
112-
// computeResidualAndJacobianTest< double, 2 >( simpleTestRateParams,
113-
// initialSpeciesConcentration,
114-
// expectedResiduals,
115-
// expectedJacobian );
116-
// }
117-
118-
// }
78+
TEST( testEquilibriumReactions, computeResidualAndJacobianTest )
79+
{
80+
double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
81+
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+
}
105+
106+
{
107+
std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
108+
double const expectedResiduals[] = { -37.534508668465, -72.989575795250 };
109+
double const expectedJacobian[2][2] =
110+
{ { 1.0e16,-2.0 },
111+
{ -2.0,4.0e16 } };
112+
computeResidualAndJacobianTest< double, 2 >( simpleTestRateParams,
113+
initialSpeciesConcentration,
114+
expectedResiduals,
115+
expectedJacobian );
116+
}
117+
118+
}
119119

120120

121121
//******************************************************************************
@@ -150,35 +150,36 @@ void testEnforceEquilibrium( PARAMS_DATA const & params,
150150

151151
for( int r=0; r<numSpecies; ++r )
152152
{
153+
// printf( "c[%d] = %22.14e\n", r, speciesConcentration[r] );
153154
EXPECT_NEAR( speciesConcentration[r], expectedSpeciesConcentrations[r], 1.0e-8 * expectedSpeciesConcentrations[r] );
154155
}
155156

156157
}
157158

158159

159-
// //******************************************************************************
160-
// TEST( testEquilibriumReactions, testEnforceEquilibrium )
161-
// {
162-
// double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
163-
// double const expectedSpeciesConcentrations[5] = { 3.92138294e-01, 3.03930853e-01, 5.05945481e-01, 7.02014628e-01, 5.95970745e-01 };
160+
//******************************************************************************
161+
TEST( testEquilibriumReactions, testEnforceEquilibrium )
162+
{
163+
double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
164+
double const expectedSpeciesConcentrations[5] = { 3.92138294e-01, 3.03930853e-01, 5.05945481e-01, 7.02014628e-01, 5.95970745e-01 };
164165

165-
// std::cout<<" RESIDUAL_FORM 0:"<<std::endl;
166-
// testEnforceEquilibrium< double, 0 >( simpleTestRateParams.equilibriumReactionsParameters(),
167-
// initialSpeciesConcentration,
168-
// expectedSpeciesConcentrations );
166+
std::cout<<" RESIDUAL_FORM 0:"<<std::endl;
167+
testEnforceEquilibrium< double, 0 >( simpleTestRateParams.equilibriumReactionsParameters(),
168+
initialSpeciesConcentration,
169+
expectedSpeciesConcentrations );
169170

170-
// // std::cout<<std::endl;
171-
// // std::cout<<" RESIDUAL_FORM 1:"<<std::endl;
172-
// // testEnforceEquilibrium< double, 1 >( simpleTestRateParams.equilibriumReactionsParameters(),
173-
// // initialSpeciesConcentration,
174-
// // expectedSpeciesConcentrations );
171+
// std::cout<<std::endl;
172+
// std::cout<<" RESIDUAL_FORM 1:"<<std::endl;
173+
// testEnforceEquilibrium< double, 1 >( simpleTestRateParams.equilibriumReactionsParameters(),
174+
// initialSpeciesConcentration,
175+
// expectedSpeciesConcentrations );
175176

176-
// std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
177-
// testEnforceEquilibrium< double, 2 >( simpleTestRateParams.equilibriumReactionsParameters(),
178-
// initialSpeciesConcentration,
179-
// expectedSpeciesConcentrations );
177+
std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
178+
testEnforceEquilibrium< double, 2 >( simpleTestRateParams.equilibriumReactionsParameters(),
179+
initialSpeciesConcentration,
180+
expectedSpeciesConcentrations );
180181

181-
// }
182+
}
182183

183184

184185
//******************************************************************************
@@ -227,6 +228,16 @@ TEST( testEquilibriumReactions, testCarbonateSystem )
227228
1.072307827865370e+00 // Na+1
228229
};
229230

231+
std::cout<<" RESIDUAL_FORM 0:"<<std::endl;
232+
testEnforceEquilibrium< double, 0 >( carbonateSystem.equilibriumReactionsParameters(),
233+
initialSpeciesConcentration,
234+
expectedSpeciesConcentrations );
235+
236+
// std::cout<<" RESIDUAL_FORM 1:"<<std::endl;
237+
// testEnforceEquilibrium< double, 1 >( carbonateSystem.equilibriumReactionsParameters(),
238+
// initialSpeciesConcentration,
239+
// expectedSpeciesConcentrations );
240+
230241
std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
231242
testEnforceEquilibrium< double, 2 >( carbonateSystem.equilibriumReactionsParameters(),
232243
initialSpeciesConcentration,

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

Lines changed: 18 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -80,7 +80,8 @@ TEST( testKineticReactions, computeReactionRatesTest )
8080
{
8181
double const initialSpeciesConcentration[] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
8282
double const expectedReactionRates[] = { 1.0, 0.25 };
83-
double const expectedReactionRatesDerivatives[][5] = { { 2.0, -0.5, 0.0, 0.0, 0.0 },
83+
double const expectedReactionRatesDerivatives[][5] =
84+
{ { 2.0, -0.5, 0.0, 0.0, 0.0 },
8485
{ 0.0, 0.0, 0.5, 0.25, 0.0 } };
8586
computeReactionRatesTest< double, false >( simpleTestRateParams.kineticReactionsParameters(),
8687
initialSpeciesConcentration,
@@ -264,31 +265,31 @@ void timeStepTest( PARAMS_DATA const & params,
264265

265266
TEST( testKineticReactions, testTimeStep )
266267
{
267-
double const initialSpeciesConcentration[5] = { 0.9, 1.0e-16, 0.5, 0.9, 1.0e-16 };
268-
double const expectedSpeciesConcentrations[5] = { 3.92138294e-01, 3.03930853e-01, 5.05945481e-01, 7.02014628e-01, 5.95970745e-01 };
268+
double const initialSpeciesConcentration[5] = { 1.0, 1.0e-16, 0.5, 1.0, 1.0e-16 };
269+
double const expectedSpeciesConcentrations[5] = { 3.92138293924124e-01, 3.03930853037938e-01, 5.05945480771998e-01, 7.02014627734060e-01, 5.95970744531880e-01 };
269270
double const expectedSpeciesRates[5] = { -2.0, 1.0, 0.75, -0.25, 0.5 };
270271
double const expectedSpeciesRatesDerivatives[5][5] = { { -4.0, 1.0, 0.0, 0.0, 0.0 },
271272
{ 2.0, -0.5, 0.0, 0.0, 0.0 },
272273
{ 2.0, -0.5, -0.5, -0.25, 0.0 },
273274
{ 0.0, 0.0, -0.5, -0.25, 0.0 },
274275
{ 0.0, 0.0, 1.0, 0.5, 0.0 } };
275276

276-
// timeStepTest< double, false >( simpleTestRateParams,
277-
// 2.0,
278-
// 10,
279-
// initialSpeciesConcentration,
280-
// expectedSpeciesConcentrations,
281-
// expectedSpeciesRates,
282-
// expectedSpeciesRatesDerivatives );
277+
timeStepTest< double, false >( simpleTestRateParams,
278+
2.0,
279+
10,
280+
initialSpeciesConcentration,
281+
expectedSpeciesConcentrations,
282+
expectedSpeciesRates,
283+
expectedSpeciesRatesDerivatives );
283284

284285
// 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 );
286+
// timeStepTest< double, true >( simpleTestRateParams,
287+
// 2.0,
288+
// 10,
289+
// initialSpeciesConcentration,
290+
// expectedSpeciesConcentrations,
291+
// expectedSpeciesRates,
292+
// expectedSpeciesRatesDerivatives );
292293
}
293294

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

0 commit comments

Comments
 (0)