Skip to content

Commit 2ffbbc4

Browse files
rrsettgastrrsettgast
committed
rename utilities to SpeciesUtilities
1 parent 668db6a commit 2ffbbc4

File tree

6 files changed

+549
-205
lines changed

6 files changed

+549
-205
lines changed

‎src/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ set( hpcReact_headers
77
reactions/bulkGeneric/KineticReactions_impl.hpp
88
reactions/bulkGeneric/Parameters.hpp
99
reactions/bulkGeneric/ParametersPredefined.hpp
10+
reactions/bulkGeneric/SpeciesUtilities.hpp
1011
)
1112

1213
set( hpcReact_sources

‎src/reactions/bulkGeneric/SpeciesUtilities.hpp‎

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
#pragma once
2+
3+
#include "common/macros.hpp"
4+
#include <math.h>
5+
#include <functional>
6+
7+
namespace hpcReact
8+
{
9+
namespace bulkGeneric
10+
{
11+
12+
namespace utilities_impl
13+
{
14+
15+
template< typename REAL_TYPE,
16+
typename INT_TYPE,
17+
typename INDEX_TYPE,
18+
typename PARAMS_DATA,
19+
typename ARRAY_1D_TO_CONST,
20+
typename ARRAY_1D,
21+
typename FUNC >
22+
HPCREACT_HOST_DEVICE
23+
inline
24+
void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params,
25+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
26+
ARRAY_1D & logSecondarySpeciesConcentrations,
27+
FUNC && derivativeFunc )
28+
{
29+
constexpr int numSpecies = PARAMS_DATA::numSpecies;
30+
constexpr int numSecondarySpecies = PARAMS_DATA::numReactions;
31+
constexpr int numPrimarySpecies = numSpecies - numSecondarySpecies;
32+
33+
for( int j=0; j<numSecondarySpecies; ++j )
34+
{
35+
REAL_TYPE const gamma = 1;
36+
logSecondarySpeciesConcentrations[j] = -log( params.equilibriumConstant( j ) ) - log( gamma );
37+
for( int k=0; k<numPrimarySpecies; ++k )
38+
{
39+
logSecondarySpeciesConcentrations[j] += params.stoichiometricMatrix( j, k+numSecondarySpecies ) * ( logPrimarySpeciesConcentrations[k] + log( gamma ) );
40+
derivativeFunc( j, k, params.stoichiometricMatrix( j, k+numSecondarySpecies ) );
41+
}
42+
}
43+
}
44+
45+
}
46+
47+
48+
template< typename REAL_TYPE,
49+
typename INT_TYPE,
50+
typename INDEX_TYPE,
51+
typename PARAMS_DATA,
52+
typename ARRAY_1D_TO_CONST,
53+
typename ARRAY_1D >
54+
HPCREACT_HOST_DEVICE
55+
inline
56+
void calculateLogSecondarySpeciesConcentration( PARAMS_DATA const & params,
57+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
58+
ARRAY_1D & logSecondarySpeciesConcentrations )
59+
{
60+
utilities_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE,
61+
INT_TYPE,
62+
INDEX_TYPE >( params,
63+
logPrimarySpeciesConcentrations,
64+
logSecondarySpeciesConcentrations,
65+
[](INDEX_TYPE, INDEX_TYPE, REAL_TYPE ){} );
66+
}
67+
68+
69+
template< typename REAL_TYPE,
70+
typename INT_TYPE,
71+
typename INDEX_TYPE,
72+
typename PARAMS_DATA,
73+
typename ARRAY_1D_TO_CONST,
74+
typename ARRAY_1D,
75+
typename ARRAY_2D >
76+
HPCREACT_HOST_DEVICE
77+
inline
78+
void calculateLogSecondarySpeciesConcentrationWrtLogC( PARAMS_DATA const & params,
79+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
80+
ARRAY_1D & logSecondarySpeciesConcentrations,
81+
ARRAY_2D & dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations )
82+
{
83+
utilities_impl::calculateLogSecondarySpeciesConcentration< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params,
84+
logPrimarySpeciesConcentrations,
85+
logSecondarySpeciesConcentrations,
86+
[&]( const int j, const int k, REAL_TYPE const value )
87+
{
88+
dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[j][k] = value;
89+
} );
90+
}
91+
92+
template< typename REAL_TYPE,
93+
typename INT_TYPE,
94+
typename INDEX_TYPE,
95+
typename PARAMS_DATA,
96+
typename ARRAY_1D_TO_CONST,
97+
typename ARRAY_1D,
98+
typename ARRAY_2D >
99+
HPCREACT_HOST_DEVICE
100+
inline
101+
void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
102+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
103+
ARRAY_1D & aggregatePrimarySpeciesConcentrations,
104+
ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations )
105+
{
106+
constexpr int numSpecies = PARAMS_DATA::numSpecies;
107+
constexpr int numSecondarySpecies = PARAMS_DATA::numReactions;
108+
constexpr int numPrimarySpecies = numSpecies - numSecondarySpecies;
109+
110+
REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0};
111+
112+
calculateLogSecondarySpeciesConcentration< REAL_TYPE,
113+
INT_TYPE,
114+
INDEX_TYPE >( params,
115+
logPrimarySpeciesConcentrations,
116+
logSecondarySpeciesConcentrations );
117+
118+
for( int i = 0; i < numPrimarySpecies; ++i )
119+
{
120+
REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] );
121+
aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
122+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
123+
for( int j = 0; j < numSecondarySpecies; ++j )
124+
{
125+
REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] );
126+
aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
127+
for( int k=0; k<numPrimarySpecies; ++k )
128+
{
129+
REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = params.stoichiometricMatrix( j, k+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
130+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j,
131+
i+numSecondarySpecies ) *
132+
dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration;
133+
}
134+
}
135+
}
136+
}
137+
138+
// template< typename REAL_TYPE,
139+
// typename INT_TYPE,
140+
// typename INDEX_TYPE,
141+
// typename PARAMS_DATA,
142+
// typename ARRAY_1D_TO_CONST,
143+
// typename ARRAY_1D,
144+
// typename ARRAY_2D >
145+
// void calculateAggregateBasedResidualAndJacobianWrtLogC( )
146+
147+
} // namespace bulkGeneric
148+
} // namespace hpcReact

‎src/reactions/bulkGeneric/unitTests/CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
set( testSourceFiles
33
testEquilibriumReactions.cpp
44
testKineticReactions.cpp
5-
testUtilities.cpp
5+
testSpeciesUtilities.cpp
66
)
77

88
set( dependencyList hpcReact gtest )

src/reactions/bulkGeneric/unitTests/testUtilities.cpp renamed to src/reactions/bulkGeneric/unitTests/testSpeciesUtilities.cpp

Lines changed: 89 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
#include "../utilities.hpp"
1+
#include "../SpeciesUtilities.hpp"
22
#include "../ParametersPredefined.hpp"
33
#include "common/printers.hpp"
44

@@ -8,6 +8,94 @@ using namespace hpcReact;
88
using namespace hpcReact::bulkGeneric;
99

1010

11+
12+
TEST( testUtilities, test_calculateLogSecondarySpeciesConcentration )
13+
{
14+
constexpr int numReactions = carbonateSystem.numReactions;
15+
constexpr int numSpecies = carbonateSystem.numSpecies;
16+
constexpr int numPrimarySpecies = numSpecies - numReactions;
17+
constexpr int numSecondarySpecies = numReactions;
18+
19+
double const logPrimarySpeciesSolution[numPrimarySpecies] =
20+
{
21+
log( 0.00043969547214915125 ),
22+
log( 0.00037230096984514874 ),
23+
log( 0.014716565308128551 ),
24+
log( 0.0024913722747387217 ),
25+
log( 1.8586090945989489 ),
26+
log( 0.009881874292035079 ),
27+
log( 1.0723078278653704 )
28+
};
29+
30+
double logSecondarySpeciesConcentrations[numSecondarySpecies] = {0};
31+
32+
calculateLogSecondarySpeciesConcentration< double,
33+
int,
34+
int >( carbonateSystem,
35+
logPrimarySpeciesSolution,
36+
logSecondarySpeciesConcentrations );
37+
38+
double expectedSecondarySpeciesConcentrations[numSecondarySpecies] =
39+
{
40+
2.3278416955869804e-11,
41+
0.3745973700632361,
42+
3.956656978189425e-11,
43+
0.0009629355924566718,
44+
0.00006739226982791149,
45+
1.065032288527949e-9,
46+
0.005298329882666744,
47+
0.005844517547638335,
48+
0.012773193926706526,
49+
0.006618125707964999,
50+
0.01769217213462983
51+
};
52+
53+
for( int j=0; j<numSecondarySpecies; ++j )
54+
{
55+
EXPECT_NEAR( exp( logSecondarySpeciesConcentrations[j] ),
56+
expectedSecondarySpeciesConcentrations[j],
57+
1.0e-8 );
58+
}
59+
60+
61+
double dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[numSecondarySpecies][numPrimarySpecies] = {{0}};
62+
63+
calculateLogSecondarySpeciesConcentrationWrtLogC< double,
64+
int,
65+
int >( carbonateSystem,
66+
logPrimarySpeciesSolution,
67+
logSecondarySpeciesConcentrations,
68+
dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations );
69+
70+
71+
72+
double expected_dLogSdLogC[numSecondarySpecies][numPrimarySpecies] =
73+
{
74+
{ -1, 0, 0, 0, 0, 0, 0 },
75+
{ 1, 1, 0, 0, 0, 0, 0 },
76+
{ -1, 1, 0, 0, 0, 0, 0 },
77+
{ 1, 1, 0, 0, 0, 0, 0 },
78+
{ 0, 1, 1, 0, 0, 0, 0 },
79+
{ -1, 1, 1, 0, 0, 0, 0 },
80+
{ 0, 0, 1, 1, 0, 0, 0 },
81+
{ 0, 0, 1, 0, 1, 0, 0 },
82+
{ 0, 0, 1, 0, 2, 0, 0 },
83+
{ 0, 0, 0, 1, 0, 1, 0 },
84+
{ 0, 0, 0, 1, 0, 0, 1 }
85+
};
86+
87+
for( int i=0; i<numSecondarySpecies; ++i )
88+
{
89+
for( int j=0; j<numPrimarySpecies; ++j )
90+
{
91+
EXPECT_NEAR( dLogSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentrations[i][j],
92+
expected_dLogSdLogC[i][j],
93+
1.0e-8 );
94+
}
95+
}
96+
}
97+
98+
1199
TEST( testUtilities, testcalculateAggregatePrimaryConcentrationsWrtLogC )
12100
{
13101
constexpr int numReactions = carbonateSystem.numReactions;
@@ -72,9 +160,6 @@ TEST( testUtilities, testcalculateAggregatePrimaryConcentrationsWrtLogC )
72160
};
73161

74162

75-
76-
print( dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations, "dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations", 4 );
77-
78163
for( int i=0; i<numPrimarySpecies; ++i )
79164
{
80165
for( int j=0; j<numPrimarySpecies; ++j )

0 commit comments

Comments
 (0)