Skip to content

Commit 877df19

Browse files
frankfeifanfrankfeifan
committed
initial commit
1 parent 057bb05 commit 877df19

File tree

10 files changed

+190
-25
lines changed

10 files changed

+190
-25
lines changed

‎src/reactions/exampleSystems/BulkGeneric.hpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -61,7 +61,9 @@ simpleKineticTestRateParams =
6161
// forward rate constants
6262
{ 1.0, 0.5 },
6363
// reverse rate constants
64-
{ 1.0, 0.5 }
64+
{ 1.0, 0.5 },
65+
// flag of mobile secondary species
66+
{ 1, 1 }
6567
};
6668

6769
using simpleTestType = reactionsSystems::MixedReactionsParameters< double, int, int, 5, 2, 2 >;
@@ -80,7 +82,9 @@ simpleTestRateParams =
8082
// forward rate constants
8183
{ 1.0, 0.5 },
8284
// reverse rate constants
83-
{ 1.0, 0.5 }
85+
{ 1.0, 0.5 },
86+
// flag of mobile secondary species
87+
{ 1, 1 }
8488
};
8589

8690
// *****UNCRUSTIFY-ON******

‎src/reactions/exampleSystems/MoMasBenchmark.hpp‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,16 @@ namespace MomMasBenchmark
6767
0.0,
6868
0.0
6969
},
70+
71+
// Flag of mobile secondary species
72+
{ 1,
73+
1,
74+
1,
75+
1,
76+
1,
77+
0,
78+
0
79+
}
7080
};
7181

7282
// *****UNCRUSTIFY-ON******

‎src/reactions/geochemistry/Carbonate.hpp‎

Lines changed: 20 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,7 @@ constexpr CArrayWrapper<double, 11> equilibriumConstants =
6969
3.98E+00, // CaCl2 = Ca+2 + 2Cl-
7070
3.72E-03, // MgSO4 = Mg+2 + SO4-2
7171
1.51E-01, // NaSO4- = Na+ + SO4-2
72-
1.17E+07 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
72+
70.55 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
7373
// 1
7474
}; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic)
7575

@@ -85,7 +85,7 @@ constexpr CArrayWrapper<double, 11> forwardRates =
8585
1.0e7, // CaCl2 = Ca+2 + 2Cl-
8686
1.0e5, // MgSO4 = Mg+2 + SO4-2
8787
1.0e7, // NaSO4- = Na+ + SO4-2
88-
1.0e5 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
88+
1.55e-6 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
8989

9090
// 1
9191
}; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic)
@@ -101,17 +101,31 @@ constexpr CArrayWrapper<double, 11> reverseRates =
101101
2.51E+06, // CaCl2 = Ca+2 + 2Cl-
102102
2.69E+07, // MgSO4 = Mg+2 + SO4-2
103103
6.62E+07, // NaSO4- = Na+ + SO4-2
104-
8.55E-03 // CaCO3 + H+ = Ca+2 + HCO3-
104+
2.197e-8 // CaCO3 + H+ = Ca+2 + HCO3-
105105
// 1 // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic)
106106
};
107+
108+
constexpr CArrayWrapper<int, 11> mobileSpeciesFlag =
109+
{ 1, // OH- + H+ = H2O
110+
1, // CO2 + H2O = H+ + HCO3-
111+
1, // CO3-2 + H+ = HCO3-
112+
1, // H2CO3 = H+ + HCO3-
113+
1, // CaHCO3+ = Ca+2 + HCO3-
114+
1, // CaSO4 = Ca+2 + SO4-2
115+
1, // CaCl+ = Ca+2 + Cl-
116+
1, // CaCl2 = Ca+2 + 2Cl-
117+
1, // MgSO4 = Mg+2 + SO4-2
118+
1, // NaSO4- = Na+ + SO4-2
119+
1 // CaCO3 + H+ = Ca+2 + HCO3-
120+
};
107121
}
108122
using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 0 >;
109123
using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 11 >;
110124
using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 11, 10 >;
111125

112-
constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
113-
constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
114-
constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );
126+
constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
127+
constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
128+
constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
115129

116130
// *****UNCRUSTIFY-ON******
117131
} // namespace geochemistry

‎src/reactions/geochemistry/Forge.hpp‎

Lines changed: 24 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -119,12 +119,35 @@ constexpr CArrayWrapper< double, 19 > reverseRateConstant =
119119
1.0 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq)
120120
};
121121

122+
constexpr CArrayWrapper< int, 19 > mobileSpeciesFlag =
123+
{
124+
1, // CaCO₃(aq) + H⁺ ⇌ Ca²⁺ + HCO₃⁻
125+
1, // CaHCO₃⁺ ⇌ Ca²⁺ + HCO₃⁻
126+
1, // CaSO₄ ⇌ Ca²⁺ + SO₄²⁻
127+
1, // CaCl⁺ ⇌ Ca²⁺ + Cl⁻
128+
1, // CaCl₂ ⇌ Ca²⁺ + 2Cl⁻ (approximate, same source)
129+
1, // MgHCO₃⁺ ⇌ Mg²⁺ + HCO₃⁻
130+
1, // MgCO₃(aq) + H⁺ ⇌ Mg²⁺ + HCO₃⁻
131+
1, // MgCl⁺ ⇌ Mg²⁺ + Cl⁻
132+
1, // CO₂(aq) + H₂O ⇌ H⁺ + HCO₃⁻
133+
1, // HSO₄⁻ ⇌ H⁺ + SO₄²⁻
134+
1, // KHSO₄ ⇌ H⁺ + K⁺ + SO₄²⁻
135+
1, // HSiO₃⁻ ⇌ H⁺ + SiO₂(aq)
136+
1, // NaHSilO₃ ⇌ H⁺ + Na⁺ + SiO₂(aq)
137+
1, // NaCl ⇌ Na⁺ + Cl⁻
138+
1, // KCl ⇌ K⁺ + Cl⁻
139+
1, // KSO₄⁻ ⇌ K⁺ + SO₄²⁻
140+
1, // Dolomite: CaMg(CO₃)₂(s) + 2H⁺ ⇌ Ca²⁺ + Mg²⁺ + 2HCO₃⁻
141+
1, // Microcline: KAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + K⁺ + 3SiO₂(aq)
142+
1 // Albite: NaAlSi₃O₈(s) + 4H⁺ ⇌ Al³⁺ + Na⁺ + 3SiO₂(aq)
143+
};
144+
122145
}
123146

124147
using forgeSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 26, 19, 16 >;
125148

126149

127-
constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant );
150+
constexpr forgeSystemType forgeSystem( forge::soichMatrix, forge::equilibriumConstants, forge::fwRateConstant, forge::reverseRateConstant, forge::mobileSpeciesFlag );
128151

129152

130153
// *****UNCRUSTIFY-ON******

‎src/reactions/geochemistry/Ultramafics.hpp‎

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -126,15 +126,40 @@ constexpr CArrayWrapper<double, 21> reverseRates =
126126
2.83E-44, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O
127127
2.10E-25 // Mg(OH)2 + 2H+ = Mg++ + 2H2O
128128
};
129+
130+
constexpr CArrayWrapper<int, 21> mobileSpeciesFlag =
131+
{
132+
1, // OH- + H+ = H2O
133+
1, // CO2(aq) + H2O = HCO3- + H+
134+
1, // CO3-- + H+ = HCO3-
135+
1, // Mg2OH+++ + H+ = 2Mg++ + H2O
136+
1, // Mg4(OH)++++ + 4H+ = 4Mg++ + 4H2O
137+
1, // MgOH+ + H+ = Mg++ + H2O
138+
1, // Mg2CO3++ + H+ = 2Mg++ + HCO3-
139+
1, // MgCO3 + H+ = Mg++ + HCO3-
140+
1, // MgHCO3+ = Mg++ + HCO3-
141+
1, // Mg(H3SiO4)2 + 2H+ = Mg++ + SiO2(aq) + 4H2O
142+
1, // MgH2SiO4 + 2H+ = Mg++ + SiO2(aq) + 2H2O
143+
1, // MgH3SiO4+ + H+ = Mg++ + SiO2(aq) + 2H2O
144+
1, // H2SiO4-- + 2H+ = SiO2(aq) + 2H2O
145+
1, // H3SiO4- + H+ = SiO2(aq) + 2H2O
146+
1, // H4(H2SiO4)---- + 4H+ = 4SiO2(aq) + 8H2O
147+
1, // H6(H2SiO4)-- + 2H+ = 4SiO2 + 8H2O
148+
1, // Mg2SiO4 + 4H+ = 2Mg++ + SiO2(aq) + 2H2O
149+
1, // MgCO3 + H+ = Mg++ + HCO3-
150+
1, // SiO2 = SiO2(aq)
151+
1, // Mg3Si2O5(OH)4 + 6H+ = 3Mg++ + 2SiO2(aq) + 5H2O
152+
1 // Mg(OH)2 + 2H+ = Mg++ + 2H2O
153+
};
129154
};
130155

131156
using ultramaficSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 0 >;
132157
using ultramaficSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 21 >;
133158
using ultramaficSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 25, 21, 16 >;
134159

135-
constexpr ultramaficSystemAllKineticType ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
136-
constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
137-
constexpr ultramaficSystemType ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates );
160+
constexpr ultramaficSystemAllKineticType ultramaficSystemAllKinetic( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
161+
constexpr ultramaficSystemAllEquilibriumType ultramaficSystemAllEquilibrium( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
162+
constexpr ultramaficSystemType ultramaficSystem( ultramafics::stoichMatrix, ultramafics::equilibriumConstants, ultramafics::forwardRates, ultramafics::reverseRates, ultramafics::mobileSpeciesFlag );
138163

139164
// *****UNCRUSTIFY-ON******
140165
} // namespace geochemistry

‎src/reactions/massActions/MassActions.hpp‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -210,5 +210,68 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
210210
}
211211
}
212212

213+
template< typename REAL_TYPE,
214+
typename INT_TYPE,
215+
typename INDEX_TYPE,
216+
typename PARAMS_DATA,
217+
typename ARRAY_1D_TO_CONST,
218+
typename ARRAY_1D_PRIMARY,
219+
typename ARRAY_1D_SECONDARY,
220+
typename ARRAY_2D >
221+
HPCREACT_HOST_DEVICE
222+
inline
223+
void calculateTotalAndMobileAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
224+
ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentrations,
225+
ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations,
226+
ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations,
227+
ARRAY_1D_PRIMARY & mobileAggregatePrimarySpeciesConcentrations,
228+
ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations,
229+
ARRAY_2D & dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations )
230+
{
231+
constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
232+
constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
233+
234+
calculateLogSecondarySpeciesConcentration< REAL_TYPE,
235+
INT_TYPE,
236+
INDEX_TYPE >( params,
237+
logPrimarySpeciesConcentrations,
238+
logSecondarySpeciesConcentrations );
239+
for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i )
240+
{
241+
for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j )
242+
{
243+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0;
244+
dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0;
245+
}
246+
}
247+
248+
for( int i = 0; i < numPrimarySpecies; ++i )
249+
{
250+
REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] );
251+
aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
252+
mobileAggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
253+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
254+
dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
255+
256+
for( int j = 0; j < numSecondarySpecies; ++j )
257+
{
258+
REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] );
259+
aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
260+
mobileAggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j * params.mobileSecondarySpeciesFlag( j );
261+
for( int k=0; k<numPrimarySpecies; ++k )
262+
{
263+
REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = params.stoichiometricMatrix( j, k+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
264+
dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j,
265+
i+numSecondarySpecies ) *
266+
dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration;
267+
268+
dMobileAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j,
269+
i+numSecondarySpecies ) *
270+
dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration * params.mobileSecondarySpeciesFlag( j );
271+
}
272+
}
273+
}
274+
}
275+
213276
} // namespace massActions
214277
} // namespace hpcReact

‎src/reactions/reactionsSystems/MixedEquilibriumKineticReactions.hpp‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,7 +72,9 @@ class MixedEquilibriumKineticReactions
7272
* @param surfaceArea surface Aread for kinetic reactions
7373
* @param logSecondarySpeciesConcentrations Output log concentrations for secondary species
7474
* @param aggregatePrimarySpeciesConcentrations Output aggregate concentrations (per primary)
75+
* @param mobileAggregatePrimarySpeciesConcentrations Output mobile aggregate concentrations (per primary)
7576
* @param dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations Derivatives of aggregate concentrations w.r.t. log
77+
* @param dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations Derivatives of mobile aggregate concentrations w.r.t. log
7678
* primary
7779
* @param reactionRates Output vector of kinetic reaction rates
7880
* @param dReactionRates_dLogPrimarySpeciesConcentrations Derivatives of reaction rates w.r.t. log primary species
@@ -94,7 +96,9 @@ class MixedEquilibriumKineticReactions
9496
ARRAY_1D_TO_CONST_KINETIC const & surfaceArea,
9597
ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations,
9698
ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations,
99+
ARRAY_1D_PRIMARY & mobileAggregatePrimarySpeciesConcentrations,
97100
ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
101+
ARRAY_2D_PRIMARY & dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
98102
ARRAY_1D_KINETIC & reactionRates,
99103
ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations,
100104
ARRAY_1D_PRIMARY & aggregateSpeciesRates,
@@ -106,7 +110,9 @@ class MixedEquilibriumKineticReactions
106110
surfaceArea,
107111
logSecondarySpeciesConcentrations,
108112
aggregatePrimarySpeciesConcentrations,
113+
mobileAggregatePrimarySpeciesConcentrations,
109114
dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
115+
dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
110116
reactionRates,
111117
dReactionRates_dLogPrimarySpeciesConcentrations,
112118
aggregateSpeciesRates,
@@ -219,7 +225,9 @@ class MixedEquilibriumKineticReactions
219225
* @param logPrimarySpeciesConcentrations Log of primary species concentrations
220226
* @param logSecondarySpeciesConcentrations Output log concentrations for secondary species
221227
* @param aggregatePrimarySpeciesConcentrations Output aggregate concentrations (per primary)
228+
* @param mobileAggregatePrimarySpeciesConcentrations Output mobile aggregate concentrations (per primary)
222229
* @param dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations Derivatives of aggregate concentrations w.r.t. log
230+
* @param dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations Derivatives of mobile aggregate concentrations w.r.t. log
223231
* primary
224232
* @param reactionRates Output vector of kinetic reaction rates
225233
* @param dReactionRates_dLogPrimarySpeciesConcentrations Derivatives of reaction rates w.r.t. log primary species
@@ -241,7 +249,9 @@ class MixedEquilibriumKineticReactions
241249
ARRAY_1D_TO_CONST_KINETIC const & surfaceArea,
242250
ARRAY_1D_SECONDARY & logSecondarySpeciesConcentrations,
243251
ARRAY_1D_PRIMARY & aggregatePrimarySpeciesConcentrations,
252+
ARRAY_1D_PRIMARY & mobileAggregatePrimarySpeciesConcentrations,
244253
ARRAY_2D_PRIMARY & dAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
254+
ARRAY_2D_PRIMARY & dMobileAggregatePrimarySpeciesConcentrations_dLogPrimarySpeciesConcentrations,
245255
ARRAY_1D_KINETIC & reactionRates,
246256
ARRAY_2D_KINETIC & dReactionRates_dLogPrimarySpeciesConcentrations,
247257
ARRAY_1D_PRIMARY & aggregateSpeciesRates,

0 commit comments

Comments
 (0)