Skip to content

Commit 76fceb7

Browse files
CusiniMCusiniM
committed
mixed reactions test working.
1 parent 8dbbaac commit 76fceb7

File tree

4 files changed

+87
-37
lines changed

4 files changed

+87
-37
lines changed

‎src/common/nonlinearSolvers.hpp‎

Lines changed: 46 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,39 @@ void scale( double (&x)[N], double const value )
3939

4040
}
4141

42+
namespace utils
43+
{
44+
template< int N >
45+
HPCREACT_HOST_DEVICE
46+
void print(double const (&J)[N][N], double const (&r)[N], double const (&dx)[N])
47+
{
48+
printf("=======================\nJacobian matrix:\n=======================\n");
49+
printf(" RowID ColID Value\n"); for ( int i = 0; i < N; ++i )
50+
{
51+
for ( int j = 0; j < N; ++j )
52+
{
53+
printf("%10d%16d%27.16e\n", i, j, J[i][j]);
54+
}
55+
}
56+
57+
printf("\n=======================\nSystem right-hand side:\n=======================\n");
58+
printf(" RowID Value\n");
59+
60+
for ( int i = 0; i < N; ++i )
61+
{
62+
printf("%10d%27.16e\n", i, r[i]);
63+
}
64+
65+
printf("\n=======================\nDelta update vector:\n=======================\n");
66+
printf(" RowID Value\n");
67+
68+
for ( int i = 0; i < N; ++i )
69+
{
70+
printf("%10d%27.16e\n", i, dx[i]);
71+
}
72+
}
73+
}
74+
4275
template< int N,
4376
typename REAL_TYPE,
4477
typename ResidualFunc,
@@ -74,12 +107,13 @@ template< int N,
74107
HPCREACT_HOST_DEVICE
75108
bool newtonRaphson( REAL_TYPE (&x)[N],
76109
FUNCTION_TYPE computeResidualAndJacobian,
77-
int maxIters = 25,
78-
double tol = 1e-10 )
110+
int maxIters = 12,
111+
double tol = 1e-10,
112+
bool const do_print = false )
79113
{
80-
REAL_TYPE residual[N];
81-
REAL_TYPE dx[N];
82-
REAL_TYPE jacobian[N][N];
114+
REAL_TYPE residual[N]{};
115+
REAL_TYPE dx[N]{};
116+
REAL_TYPE jacobian[N][N]{};
83117

84118
for ( int iter = 0; iter < maxIters; ++iter )
85119
{
@@ -94,8 +128,15 @@ bool newtonRaphson( REAL_TYPE (&x)[N],
94128
return true;
95129
}
96130
internal::scale<N>( residual, -1.0);
131+
132+
if( do_print )
133+
{
134+
utils::print( jacobian, residual, dx );
135+
}
136+
97137
solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx );
98138
internal::add< N >( x, dx );
139+
99140
}
100141

101142
printf( "--Newton solver error: Max iterations reached without convergence.\n" );

‎src/reactions/bulkGeneric/MixedEquilibriumKineticReactions_impl.hpp‎

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -101,13 +101,21 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
101101
constexpr IntType numKineticReactions = PARAMS_DATA::numKineticReactions();
102102

103103
RealType logSpeciesConcentration[numSpecies] {};
104-
for ( IntType i = 0; i < numSecondarySpecies; ++i )
104+
for ( INDEX_TYPE i = 0; i < numSecondarySpecies; ++i )
105105
{
106106
logSpeciesConcentration[i] = logSecondarySpeciesConcentrations[i];
107107
}
108-
for ( IntType i = numSecondarySpecies; i < numSpecies; ++i )
108+
for ( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i )
109109
{
110-
logSpeciesConcentration[i] = logPrimarySpeciesConcentrations[i];
110+
logSpeciesConcentration[i+numSecondarySpecies] = logPrimarySpeciesConcentrations[i];
111+
}
112+
113+
for( INDEX_TYPE i = 0; i < numKineticReactions; ++i )
114+
{
115+
for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j )
116+
{
117+
dReactionRates_dLogPrimarySpeciesConcentrations[i][j] = 0.0;
118+
}
111119
}
112120

113121
CArrayWrapper< RealType, numKineticReactions, numSpecies > reactionRatesDerivatives;
@@ -165,12 +173,12 @@ MixedEquilibriumKineticReactions< REAL_TYPE,
165173
constexpr IntType numKineticReactions = PARAMS_DATA::numKineticReactions();
166174
constexpr IntType numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
167175

168-
for( IntType i = 0; i < numPrimarySpecies; ++i )
176+
for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i )
169177
{
170178
aggregatesRates[i] = 0.0;
171179
if constexpr( CALCULATE_DERIVATIVES )
172180
{
173-
for( IntType j = 0; j < numPrimarySpecies; ++j )
181+
for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j )
174182
{
175183
aggregatesRatesDerivatives( i, j ) = 0.0;
176184
}

‎src/reactions/bulkGeneric/SpeciesUtilities.hpp‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,11 @@ void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params,
3030
constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
3131
constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
3232

33+
for (INDEX_TYPE i = 0; i < numSecondarySpecies; ++i)
34+
{
35+
logSecondarySpeciesConcentrations[i] = 0.0;
36+
}
37+
3338
for( int j=0; j<numSecondarySpecies; ++j )
3439
{
3540
REAL_TYPE const gamma = 1;
@@ -114,6 +119,13 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
114119
INDEX_TYPE >( params,
115120
logPrimarySpeciesConcentrations,
116121
logSecondarySpeciesConcentrations );
122+
for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i )
123+
{
124+
for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j )
125+
{
126+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0;
127+
}
128+
}
117129

118130
for( int i = 0; i < numPrimarySpecies; ++i )
119131
{
@@ -154,6 +166,14 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
154166

155167
REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0};
156168

169+
for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i )
170+
{
171+
for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j )
172+
{
173+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0;
174+
}
175+
}
176+
157177
calculateLogSecondarySpeciesConcentration< REAL_TYPE,
158178
INT_TYPE,
159179
INDEX_TYPE >( params,

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

Lines changed: 8 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -65,20 +65,11 @@ void timeStepTest( PARAMS_DATA const & params,
6565
aggregatePrimarySpeciesConcentration[i] = initialSpeciesConcentration[i];
6666
}
6767

68-
printf( "Initial equilibrium solve\n");
6968
EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature,
7069
carbonateSystem.equilibriumReactionsParameters(),
7170
logPrimarySpeciesConcentration,
7271
logPrimarySpeciesConcentration );
73-
74-
for( int i = 0; i < numPrimarySpecies; ++i )
75-
{
76-
std::cout << "aggregatePrimarySpeciesConcentration[" << i << "] = "
77-
<< aggregatePrimarySpeciesConcentration[i] << ", ";
78-
79-
std::cout << "logPrimarySpeciesConcentration[" << i << "] = "
80-
<< logPrimarySpeciesConcentration[i] << std::endl;
81-
}
72+
8273

8374
/// Time step loop
8475
double time = 0.0;
@@ -106,33 +97,23 @@ void timeStepTest( PARAMS_DATA const & params,
10697
aggregateSpeciesRates,
10798
dAggregateSpeciesRates_dlogPrimarySpeciesConcentration );
10899

109-
for( int i = 0; i < numPrimarySpecies; ++i )
110-
{
111-
std::cout << "aggregateSpeciesRates[" << i << "] = "
112-
<< aggregateSpeciesRates[i] << ", ";
113-
std::cout << std::endl;
114-
for ( int j = 0; j < numPrimarySpecies; ++j )
115-
{
116-
std::cout << "dAggregateSpeciesRates_dlogPrimarySpeciesConcentration[" << i << "][" << j << "] = "
117-
<< dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration[i][j] << ", ";
118-
}
119-
std::cout << std::endl;
120-
}
121-
100+
122101
for ( int i = 0; i < numPrimarySpecies; ++i )
123102
{
124103
r[i] = ( aggregatePrimarySpeciesConcentration[i] - aggregatePrimarySpeciesConcentration_n[i] ) - aggregateSpeciesRates[i] * dt;
125104
for ( int j = 0; j < numPrimarySpecies; ++j )
126105
{
127106
J[i][j] = dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration[i][j] - dAggregateSpeciesRates_dlogPrimarySpeciesConcentration[i][j] * dt;
128107
}
129-
}
108+
}
130109
};
131110

132111
nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian );
133112

134113
time += dt;
135114
}
115+
116+
EXPECT_NEAR( reactionRatesDerivatives( r, i ), expectedReactionRatesDerivatives[r][i], std::max( magScale, fabs( expectedReactionRatesDerivatives[r][i] ) ) * 1.0e-8 );
136117
}
137118

138119
TEST( testMixedReactions, testTimeStep_carbonateSystem )
@@ -141,7 +122,7 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem )
141122

142123
double const initialSpeciesConcentration[numPrimarySpecies] =
143124
{
144-
3.76e-3, //
125+
3.76e-3, // CaCO3
145126
3.76e-1, // H+
146127
3.76e-1, // HCO3-
147128
3.87e-2, // Ca+2
@@ -153,7 +134,7 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem )
153134

154135
double const expectedSpeciesConcentrations[numPrimarySpecies] =
155136
{
156-
1.0e-3, //
137+
1.0e-3, // CaCO3
157138
3.76e-1, // H+
158139
3.76e-1, // HCO3-
159140
3.87e-2, // Ca+2
@@ -165,7 +146,7 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem )
165146

166147
timeStepTest< double, true >( carbonateSystem,
167148
0.2,
168-
1,
149+
10,
169150
initialSpeciesConcentration,
170151
expectedSpeciesConcentrations );
171152

0 commit comments

Comments
 (0)