Skip to content

Commit 56211c5

Browse files
CusiniMCusiniM
committed
added pmpl
1 parent 76fceb7 commit 56211c5

File tree

1 file changed

+96
-78
lines changed

1 file changed

+96
-78
lines changed

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

Lines changed: 96 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include "common/macros.hpp"
66
#include "common/printers.hpp"
77
#include "common/nonlinearSolvers.hpp"
8+
#include "common/pmpl.hpp"
89

910
#include <gtest/gtest.h>
1011

@@ -30,97 +31,114 @@ void timeStepTest( PARAMS_DATA const & params,
3031
{
3132
HPCREACT_UNUSED_VAR( expectedSpeciesConcentrations );
3233

33-
using MixedReactionsType = MixedEquilibriumKineticReactions< REAL_TYPE,
34+
CArrayWrapper< REAL_TYPE, PARAMS_DATA::numPrimarySpecies() > primarySpeciesConcentration;
35+
36+
for( int i = 0; i < PARAMS_DATA::numPrimarySpecies(); ++i )
37+
{
38+
primarySpeciesConcentration[i] = initialSpeciesConcentration[i];
39+
}
40+
41+
pmpl::genericKernelWrapper( PARAMS_DATA::numPrimarySpecies(),
42+
primarySpeciesConcentration.data,
43+
[=] HPCREACT_HOST_DEVICE ( auto * speciesConcentration )
44+
{
45+
using MixedReactionsType = MixedEquilibriumKineticReactions< REAL_TYPE,
3446
int,
3547
int,
3648
LOGE_CONCENTRATION >;
37-
using EquilibriumReactionsType = EquilibriumReactions< REAL_TYPE,
38-
int,
39-
int >;
40-
41-
42-
// constexpr int numSpecies = PARAMS_DATA::numSpecies();
43-
constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
44-
constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
45-
constexpr int numKineticReactions = PARAMS_DATA::numKineticReactions();
49+
using EquilibriumReactionsType = EquilibriumReactions< REAL_TYPE,
50+
int,
51+
int >;
52+
53+
// constexpr int numSpecies = PARAMS_DATA::numSpecies();
54+
constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
55+
constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
56+
constexpr int numKineticReactions = PARAMS_DATA::numKineticReactions();
4657

47-
// define variables
48-
double const temperature = 298.15;
49-
REAL_TYPE logPrimarySpeciesConcentration[numPrimarySpecies];
50-
// must use CArrayWrapper to ensure correct capture in the lambda functions
51-
CArrayWrapper< REAL_TYPE, numSecondarySpecies > logSecondarySpeciesConcentration;
52-
CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration;
53-
CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration_n;
54-
CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration;
55-
CArrayWrapper< REAL_TYPE, numKineticReactions > reactionRates;
56-
CArrayWrapper< REAL_TYPE, numKineticReactions, numPrimarySpecies > dReactionRates_dlogPrimarySpeciesConcentration;
57-
CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregateSpeciesRates;
58-
CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dAggregateSpeciesRates_dlogPrimarySpeciesConcentration;
59-
60-
61-
// Initialize species concentrations
62-
for( int i = 0; i < numPrimarySpecies; ++i )
63-
{
64-
logPrimarySpeciesConcentration[i] = ::log( initialSpeciesConcentration[i] );
65-
aggregatePrimarySpeciesConcentration[i] = initialSpeciesConcentration[i];
66-
}
67-
68-
EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature,
69-
carbonateSystem.equilibriumReactionsParameters(),
70-
logPrimarySpeciesConcentration,
71-
logPrimarySpeciesConcentration );
72-
58+
// define variables
59+
double const temperature = 298.15;
60+
REAL_TYPE logPrimarySpeciesConcentration[numPrimarySpecies];
61+
// must use CArrayWrapper to ensure correct capture in the lambda functions
62+
CArrayWrapper< REAL_TYPE, numSecondarySpecies > logSecondarySpeciesConcentration;
63+
CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration;
64+
CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregatePrimarySpeciesConcentration_n;
65+
CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration;
66+
CArrayWrapper< REAL_TYPE, numKineticReactions > reactionRates;
67+
CArrayWrapper< REAL_TYPE, numKineticReactions, numPrimarySpecies > dReactionRates_dlogPrimarySpeciesConcentration;
68+
CArrayWrapper< REAL_TYPE, numPrimarySpecies > aggregateSpeciesRates;
69+
CArrayWrapper< REAL_TYPE, numPrimarySpecies, numPrimarySpecies > dAggregateSpeciesRates_dlogPrimarySpeciesConcentration;
7370

74-
/// Time step loop
75-
double time = 0.0;
76-
for( int t = 0; t < numSteps; ++t )
77-
{
78-
printf( "Timestep %d, Time = %.6e\n", t, time );
79-
80-
for( int i=0; i < numPrimarySpecies; ++i )
71+
// Initialize species concentrations
72+
for( int i = 0; i < numPrimarySpecies; ++i )
8173
{
82-
aggregatePrimarySpeciesConcentration_n[i] = aggregatePrimarySpeciesConcentration[i];
74+
logPrimarySpeciesConcentration[i] = ::log( speciesConcentration[i] );
75+
aggregatePrimarySpeciesConcentration[i] = speciesConcentration[i];
8376
}
77+
78+
EquilibriumReactionsType::enforceEquilibrium_Aggregate( temperature,
79+
carbonateSystem.equilibriumReactionsParameters(),
80+
logPrimarySpeciesConcentration,
81+
logPrimarySpeciesConcentration );
82+
83+
/// Time step loop
84+
double time = 0.0;
85+
for( int t = 0; t < numSteps; ++t )
86+
{
87+
printf( "Timestep %d, Time = %.6e\n", t, time );
8488

85-
auto computeResidualAndJacobian = [&] HPCREACT_HOST_DEVICE ( REAL_TYPE const (&logPrimarySpeciesConcentration)[numPrimarySpecies],
89+
for( int i=0; i < numPrimarySpecies; ++i )
90+
{
91+
aggregatePrimarySpeciesConcentration_n[i] = aggregatePrimarySpeciesConcentration[i];
92+
}
93+
94+
auto computeResidualAndJacobian = [&] HPCREACT_HOST_DEVICE ( REAL_TYPE const (&logPrimarySpeciesConcentration)[numPrimarySpecies],
8695
REAL_TYPE (&r)[numPrimarySpecies],
8796
REAL_TYPE (&J)[numPrimarySpecies][numPrimarySpecies] )
88-
{
89-
MixedReactionsType::updateMixedSystem( temperature,
90-
params,
91-
logPrimarySpeciesConcentration,
92-
logSecondarySpeciesConcentration,
93-
aggregatePrimarySpeciesConcentration,
94-
dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration,
95-
reactionRates,
96-
dReactionRates_dlogPrimarySpeciesConcentration,
97-
aggregateSpeciesRates,
98-
dAggregateSpeciesRates_dlogPrimarySpeciesConcentration );
97+
{
98+
MixedReactionsType::updateMixedSystem( temperature,
99+
params,
100+
logPrimarySpeciesConcentration,
101+
logSecondarySpeciesConcentration,
102+
aggregatePrimarySpeciesConcentration,
103+
dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration,
104+
reactionRates,
105+
dReactionRates_dlogPrimarySpeciesConcentration,
106+
aggregateSpeciesRates,
107+
dAggregateSpeciesRates_dlogPrimarySpeciesConcentration );
99108

100109

101-
for ( int i = 0; i < numPrimarySpecies; ++i )
102-
{
103-
r[i] = ( aggregatePrimarySpeciesConcentration[i] - aggregatePrimarySpeciesConcentration_n[i] ) - aggregateSpeciesRates[i] * dt;
104-
for ( int j = 0; j < numPrimarySpecies; ++j )
110+
for ( int i = 0; i < numPrimarySpecies; ++i )
105111
{
106-
J[i][j] = dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration[i][j] - dAggregateSpeciesRates_dlogPrimarySpeciesConcentration[i][j] * dt;
107-
}
108-
}
109-
};
112+
r[i] = ( aggregatePrimarySpeciesConcentration[i] - aggregatePrimarySpeciesConcentration_n[i] ) - aggregateSpeciesRates[i] * dt;
113+
for ( int j = 0; j < numPrimarySpecies; ++j )
114+
{
115+
J[i][j] = dAggregatePrimarySpeciesConcentrations_dlogPrimarySpeciesConcentration[i][j] - dAggregateSpeciesRates_dlogPrimarySpeciesConcentration[i][j] * dt;
116+
}
117+
}
118+
};
110119

111-
nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian );
120+
nonlinearSolvers::newtonRaphson< numPrimarySpecies >( logPrimarySpeciesConcentration, computeResidualAndJacobian );
112121

113-
time += dt;
122+
time += dt;
123+
}
124+
for(int i = 0; i < numPrimarySpecies; ++i )
125+
{
126+
speciesConcentration[i] = exp( logPrimarySpeciesConcentration[i] );
127+
}
128+
} );
129+
130+
// Check results
131+
for(int i = 0; i < PARAMS_DATA::numPrimarySpecies(); ++i )
132+
{
133+
EXPECT_NEAR( primarySpeciesConcentration[ i ], expectedSpeciesConcentrations[ i ], 1.0e-8 * expectedSpeciesConcentrations[ i ]);
114134
}
115-
116-
EXPECT_NEAR( reactionRatesDerivatives( r, i ), expectedReactionRatesDerivatives[r][i], std::max( magScale, fabs( expectedReactionRatesDerivatives[r][i] ) ) * 1.0e-8 );
117135
}
118136

119137
TEST( testMixedReactions, testTimeStep_carbonateSystem )
120138
{
121139
constexpr int numPrimarySpecies = carbonateSystemType::numPrimarySpecies();
122140

123-
double const initialSpeciesConcentration[numPrimarySpecies] =
141+
double const initialAggregateSpeciesConcentration[numPrimarySpecies] =
124142
{
125143
3.76e-3, // CaCO3
126144
3.76e-1, // H+
@@ -134,20 +152,20 @@ TEST( testMixedReactions, testTimeStep_carbonateSystem )
134152

135153
double const expectedSpeciesConcentrations[numPrimarySpecies] =
136154
{
137-
1.0e-3, // CaCO3
138-
3.76e-1, // H+
139-
3.76e-1, // HCO3-
140-
3.87e-2, // Ca+2
141-
3.21e-2, // SO4-2
142-
1.89, // Cl-
143-
1.65e-2, // Mg+2
144-
1.09 // Na+1
155+
3.3318075516669661e-05, // CaCO3
156+
2.5894448848121536e-05, // H+
157+
0.0062660162912796741, // HCO3-
158+
0.015741214773567921, // Ca+2
159+
0.0024602709470074127, // SO4-2
160+
1.8564927944498291, // Cl-
161+
0.0099316034080546619, // Mg+2
162+
1.0725251492409775 // Na+1
145163
};
146164

147165
timeStepTest< double, true >( carbonateSystem,
148166
0.2,
149167
10,
150-
initialSpeciesConcentration,
168+
initialAggregateSpeciesConcentration,
151169
expectedSpeciesConcentrations );
152170

153171
}

0 commit comments

Comments
 (0)