Skip to content

Commit 24c115a

Browse files
rrsettgastrrsettgast
committed
move definitions from EquilibriumReactions.hpp to EquilibriumReactions_impl.hpp
1 parent bf6a258 commit 24c115a

File tree

8 files changed

+94
-204
lines changed

8 files changed

+94
-204
lines changed

‎src/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,9 @@ set( hpcReact_headers
33
common/macros.hpp
44
common/CArrayWrapper.hpp
55
reactions/bulkGeneric/EquilibriumReactions.hpp
6+
reactions/bulkGeneric/EquilibriumReactions_impl.hpp
67
reactions/bulkGeneric/KineticReactions.hpp
8+
reactions/bulkGeneric/KineticReactions_impl.hpp
79
reactions/bulkGeneric/Parameters.hpp
810
reactions/bulkGeneric/ParametersPredefined.hpp
911
)

‎src/common/macros.hpp‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,3 +19,26 @@
1919
/// This macro is used to ignore warnings that that a variable is
2020
/// unused.
2121
#define HPCREACT_UNUSED_VAR( ... ) (void)( __VA_ARGS__ )
22+
23+
24+
#if defined( __clang__ )
25+
#define HPCREACT_NO_MISSING_BRACES( ... ) \
26+
_Pragma("clang diagnostic push") \
27+
_Pragma("clang diagnostic ignored \"-Wmissing-braces\"") \
28+
__VA_ARGS__ \
29+
_Pragma("clang diagnostic pop")
30+
#elif defined(__GNUC__)
31+
#define HPCREACT_NO_MISSING_BRACES( ... ) \
32+
_Pragma("GCC diagnostic push") \
33+
_Pragma("GCC diagnostic ignored \"-Wmissing-braces\"") \
34+
__VA_ARGS__ \
35+
_Pragma("GCC diagnostic pop")
36+
#elif defined(_MSC_VER)
37+
#define HPCREACT_NO_MISSING_BRACES( ... ) \
38+
__pragma(warning(push)) \
39+
__pragma(warning(disable : 4351)) \
40+
__VA_ARGS__ \
41+
__pragma(warning(pop))
42+
#else
43+
#define HPCREACT_NO_MISSING_BRACES( ... ) __VA_ARGS__ // No-op for unknown compilers
44+
#endif

‎src/reactions/bulkGeneric/EquilibriumReactions.hpp‎

Lines changed: 5 additions & 118 deletions
Original file line numberDiff line numberDiff line change
@@ -36,75 +36,7 @@ class EquilibriumReactions
3636
ARRAY_1D_TO_CONST const & speciesConcentration0,
3737
ARRAY_1D_TO_CONST2 const & xi,
3838
ARRAY_1D & residual,
39-
ARRAY_2D & jacobian )
40-
{
41-
HPCREACT_UNUSED_VAR( temperature );
42-
constexpr int numSpecies = PARAMS_DATA::numSpecies;
43-
constexpr int numReactions = PARAMS_DATA::numReactions;
44-
45-
// initialize the species concentration
46-
RealType speciesConcentration[numSpecies];
47-
for( IndexType i=0; i<numSpecies; ++i )
48-
{
49-
speciesConcentration[i] = speciesConcentration0[i];
50-
for( IndexType r=0; r<numReactions; ++r )
51-
{
52-
speciesConcentration[i] += params.stoichiometricMatrix( r, i ) * xi[r];
53-
}
54-
}
55-
56-
// loop over reactions
57-
for( IndexType a=0; a<numReactions; ++a )
58-
{
59-
// get the equilibrium constant for this reaction
60-
RealType const Keq = params.equilibriumConstant( a );
61-
62-
RealType forwardProduct = 1.0;
63-
RealType reverseProduct = 1.0;
64-
RealType dForwardProduct_dxi[numReactions] = {0.0};
65-
RealType dReverseProduct_dxi[numReactions] = {0.0};
66-
// loop over species
67-
for( IndexType i=0; i<numSpecies; ++i )
68-
{
69-
RealType const s_ai = params.stoichiometricMatrix( a, i );
70-
71-
if( s_ai < 0.0 )
72-
{
73-
// forward reaction
74-
forwardProduct *= pow( speciesConcentration[i], -s_ai );
75-
// derivative of forward product with respect to xi
76-
for( IndexType b=0; b<numReactions; ++b )
77-
{
78-
dForwardProduct_dxi[b] += -s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
79-
}
80-
}
81-
else if( s_ai > 0.0 )
82-
{
83-
// reverse reaction
84-
reverseProduct *= pow( speciesConcentration[i], s_ai );
85-
// derivative of reverse product with respect to xi
86-
for( IndexType b=0; b<numReactions; ++b )
87-
{
88-
dReverseProduct_dxi[b] += s_ai / speciesConcentration[i] * params.stoichiometricMatrix( b, i );
89-
}
90-
}
91-
}
92-
// compute the residual for this reaction
93-
residual[a] = forwardProduct - Keq * reverseProduct;
94-
// printf( "residual[%d] = %g - %g * %g = %g\n", a, forwardProduct, Keq, reverseProduct, residual[a] );
95-
96-
// Finish the derivatives of the product terms with respect to xi
97-
for( IndexType b=0; b<numReactions; ++b )
98-
{
99-
dForwardProduct_dxi[b] *= forwardProduct;
100-
dReverseProduct_dxi[b] *= reverseProduct;
101-
// printf( "(%d) dForwardProduct_dxi[%d] = %f\n", a, b, dForwardProduct_dxi[b] );
102-
// printf( "(%d) dReverseProduct_dxi[%d] = %f\n", a, b, dReverseProduct_dxi[b] );
103-
// printf( "Keq = %f\n", Keq );
104-
jacobian( a, b ) = dForwardProduct_dxi[b] - Keq * dReverseProduct_dxi[b];
105-
}
106-
}
107-
}
39+
ARRAY_2D & jacobian );
10840

10941

11042
template< typename PARAMS_DATA,
@@ -115,58 +47,13 @@ class EquilibriumReactions
11547
enforceEquilibrium( RealType const & temperature,
11648
PARAMS_DATA const & params,
11749
ARRAY_1D_TO_CONST const & speciesConcentration0,
118-
ARRAY_1D & speciesConcentration )
119-
{
120-
HPCREACT_UNUSED_VAR( temperature );
121-
constexpr int numSpecies = PARAMS_DATA::numSpecies;
122-
constexpr int numReactions = PARAMS_DATA::numReactions;
123-
double residual[numReactions] = { 0.0 };
124-
double xi[numReactions] = { 0.0 };
125-
double dxi[numReactions] = { 0.0 };
126-
CArrayWrapper< double, numReactions, numReactions > jacobian;
127-
128-
REAL_TYPE residualNorm = 0.0;
129-
for( int k=0; k<10; ++k )
130-
{
131-
computeResidualAndJacobian( temperature,
132-
params,
133-
speciesConcentration0,
134-
xi,
135-
residual,
136-
jacobian );
137-
138-
residualNorm = 0.0;
139-
for( int j = 0; j < numReactions; ++j )
140-
{
141-
residualNorm += residual[j] * residual[j];
142-
}
143-
residualNorm = sqrt( residualNorm );
144-
if( residualNorm < 1.0e-14 )
145-
{
146-
break;
147-
}
148-
149-
solveNxN_pivoted< double, numReactions >( jacobian.data, residual, dxi );
150-
151-
// solve for the change in xi
152-
for( IndexType r=0; r<numReactions; ++r )
153-
{
154-
xi[r] -= dxi[r];
155-
}
156-
}
157-
158-
for( IndexType i=0; i<numSpecies; ++i )
159-
{
160-
speciesConcentration[i] = speciesConcentration0[i];
161-
for( IndexType r=0; r<numReactions; ++r )
162-
{
163-
speciesConcentration[i] += params.stoichiometricMatrix( r, i ) * xi[r];
164-
}
165-
}
166-
}
50+
ARRAY_1D & speciesConcentration );
16751

16852
};
16953

17054

17155
} // namespace bulkGeneric
17256
} // namespace hpcReact
57+
58+
#include "EquilibriumReactions_impl.hpp"
59+
#include "common/macrosCleanup.hpp"

‎src/reactions/bulkGeneric/KineticReactions.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -175,4 +175,5 @@ class KineticReactions
175175
} // namespace bulkGeneric
176176
} // namespace hpcReact
177177

178+
#include "KineticReactions_impl.hpp"
178179
#include "common/macrosCleanup.hpp"

‎src/reactions/bulkGeneric/KineticReactions_impl.hpp‎

Lines changed: 23 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,7 @@
11
#pragma once
22

3-
#include "KineticReactions.hpp"
43
#include "common/constants.hpp"
54
#include "common/CArrayWrapper.hpp"
6-
#include "common/macros.hpp"
75
#include "common/DirectSystemSolve.hpp"
86

97
#include <cmath>
@@ -329,26 +327,29 @@ KineticReactions< REAL_TYPE,
329327
}
330328

331329

332-
printf( "residual = { ");
333-
for( int i = 0; i < numSpecies; ++i )
334-
{
335-
printf( " %g, ", residual[i] );
336-
}
337-
printf( "}\n" );
338-
339-
printf( "Jacobian = { \n" );
340-
for( int i = 0; i < numSpecies; ++i )
341-
{
342-
printf( " { ");
343-
for( int j = 0; j < numSpecies; ++j )
344-
{
345-
printf( " %g ", speciesRatesDerivatives( i, j ) );
346-
// printf( " %g ", speciesRatesDerivatives( i, j ) / exp(speciesConcentration[j]) );
347-
if( j < numSpecies-1 ) { printf( ", " ); }
348-
}
349-
printf( "}, \n" );
350-
}
351-
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" );
352353

353354
solveNxN_pivoted< double, numSpecies >( speciesRatesDerivatives.data, residual, deltaPrimarySpeciesConcentration );
354355

‎src/reactions/bulkGeneric/Parameters.hpp‎

Lines changed: 38 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,6 @@ namespace bulkGeneric
1414

1515

1616

17-
1817
template< typename REAL_TYPE,
1918
typename INT_TYPE,
2019
typename INDEX_TYPE,
@@ -32,10 +31,10 @@ struct EquilibriumReactionsParameters
3231
constexpr
3332
EquilibriumReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies],
3433
RealType const (&equilibriumConstant)[numReactions] ):
35-
EquilibriumReactionsParameters( stoichiometricMatrix,
36-
equilibriumConstant,
37-
std::make_index_sequence<NUM_REACTIONS>(),
38-
std::make_index_sequence<NUM_REACTIONS*NUM_SPECIES>() )
34+
EquilibriumReactionsParameters( stoichiometricMatrix,
35+
equilibriumConstant,
36+
std::make_index_sequence< NUM_REACTIONS >(),
37+
std::make_index_sequence< NUM_REACTIONS *NUM_SPECIES >() )
3938
{}
4039

4140

@@ -45,30 +44,18 @@ struct EquilibriumReactionsParameters
4544
RealType m_stoichiometricMatrix[numReactions][numSpecies];
4645
RealType m_equilibriumConstant[numReactions];
4746

48-
private:
49-
#if defined(__GNUC__) || defined(__clang__)
50-
#pragma GCC diagnostic push
51-
#pragma GCC diagnostic ignored "-Wmissing-braces"
52-
#elif defined(_MSC_VER)
53-
#pragma warning(push)
54-
#pragma warning(disable : 4351) // MSVC equivalent of missing-braces warning
55-
#endif
56-
57-
template< std::size_t ... R, std::size_t ... RxS >
58-
constexpr
59-
EquilibriumReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies],
60-
RealType const (&equilibriumConstant)[numReactions],
61-
std::index_sequence<R...>,
62-
std::index_sequence<RxS...> ):
63-
m_stoichiometricMatrix{ stoichiometricMatrix[RxS/numSpecies][RxS%numSpecies]... },
64-
m_equilibriumConstant{ equilibriumConstant[R]... }
65-
{}
66-
67-
#if defined(__GNUC__) || defined(__clang__)
68-
#pragma GCC diagnostic pop
69-
#elif defined(_MSC_VER)
70-
#pragma warning(pop)
71-
#endif
47+
private:
48+
HPCREACT_NO_MISSING_BRACES (
49+
template< std::size_t ... R, std::size_t ... RxS >
50+
constexpr
51+
EquilibriumReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies],
52+
RealType const (&equilibriumConstant)[numReactions],
53+
std::index_sequence< R... >,
54+
std::index_sequence< RxS... > ) :
55+
m_stoichiometricMatrix{ stoichiometricMatrix[RxS/numSpecies][RxS%numSpecies] ... },
56+
m_equilibriumConstant{ equilibriumConstant[R] ... }
57+
{}
58+
);
7259
};
7360

7461

@@ -87,13 +74,13 @@ struct KineticReactionsParameters
8774
static constexpr IndexType numReactions = NUM_REACTIONS;
8875

8976
KineticReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies],
90-
RealType const (&rateConstantForward)[numReactions],
91-
RealType const (&rateConstantReverse)[numReactions] ):
92-
KineticReactionsParameters( stoichiometricMatrix,
93-
rateConstantForward,
77+
RealType const (&rateConstantForward)[numReactions],
78+
RealType const (&rateConstantReverse)[numReactions] ):
79+
KineticReactionsParameters( stoichiometricMatrix,
80+
rateConstantForward,
9481
rateConstantReverse,
95-
std::make_index_sequence<NUM_REACTIONS>(),
96-
std::make_index_sequence<NUM_REACTIONS*NUM_SPECIES>() )
82+
std::make_index_sequence< NUM_REACTIONS >(),
83+
std::make_index_sequence< NUM_REACTIONS *NUM_SPECIES >() )
9784
{}
9885

9986

@@ -106,29 +93,20 @@ struct KineticReactionsParameters
10693
RealType m_rateConstantForward[numReactions];
10794
RealType m_rateConstantReverse[numReactions];
10895

109-
private:
110-
#if defined(__GNUC__) || defined(__clang__)
111-
#pragma GCC diagnostic push
112-
#pragma GCC diagnostic ignored "-Wmissing-braces"
113-
#elif defined(_MSC_VER)
114-
#pragma warning(push)
115-
#pragma warning(disable : 4351) // MSVC equivalent of missing-braces warning
116-
#endif
117-
template< std::size_t ... R, std::size_t ... RxS >
118-
KineticReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies],
119-
RealType const (&rateConstantForward)[numReactions],
120-
RealType const (&rateConstantReverse)[numReactions],
121-
std::index_sequence<R...>,
122-
std::index_sequence<RxS...> ):
123-
m_stoichiometricMatrix{ stoichiometricMatrix[RxS/numSpecies][RxS%numSpecies]... },
124-
m_rateConstantForward{ rateConstantForward[R]... },
125-
m_rateConstantReverse{ rateConstantReverse[R]... }
126-
{}
127-
#if defined(__GNUC__) || defined(__clang__)
128-
#pragma GCC diagnostic pop
129-
#elif defined(_MSC_VER)
130-
#pragma warning(pop)
131-
#endif
96+
private:
97+
HPCREACT_NO_MISSING_BRACES(
98+
template< std::size_t ... R, std::size_t ... RxS >
99+
KineticReactionsParameters( RealType const (&stoichiometricMatrix)[numReactions][numSpecies],
100+
RealType const (&rateConstantForward)[numReactions],
101+
RealType const (&rateConstantReverse)[numReactions],
102+
std::index_sequence< R... >,
103+
std::index_sequence< RxS... > ) :
104+
m_stoichiometricMatrix{ stoichiometricMatrix[RxS/numSpecies][RxS%numSpecies] ... },
105+
m_rateConstantForward{ rateConstantForward[R] ... },
106+
m_rateConstantReverse{ rateConstantReverse[R] ... }
107+
{}
108+
)
109+
132110
};
133111

134112

@@ -146,14 +124,14 @@ struct MixedReactionsParameters
146124
static constexpr IndexType numReactions = NUM_REACTIONS;
147125

148126
constexpr
149-
EquilibriumReactionsParameters< RealType, IntType, IndexType, numSpecies, numReactions >
127+
EquilibriumReactionsParameters< RealType, IntType, IndexType, numSpecies, numReactions >
150128
equilibriumReactionsParameters() const
151129
{
152130
return {m_stoichiometricMatrix, m_equilibriumConstant};
153131
}
154132

155133
constexpr
156-
KineticReactionsParameters< RealType, IntType, IndexType, numSpecies, numReactions >
134+
KineticReactionsParameters< RealType, IntType, IndexType, numSpecies, numReactions >
157135
kineticReactionsParameters() const
158136
{
159137
return {m_stoichiometricMatrix, m_rateConstantForward, m_rateConstantReverse};
@@ -207,7 +185,5 @@ struct MixedReactionsParameters
207185

208186

209187

210-
211-
212188
} // namespace solidStateBattery
213189
} // namespace hpcReact

0 commit comments

Comments
 (0)