Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "GeochemistryActivityCoefficientsDebyeHuckel.h"
11 : #include "GeochemistryActivityCalculators.h"
12 :
13 1210 : GeochemistryActivityCoefficientsDebyeHuckel::GeochemistryActivityCoefficientsDebyeHuckel(
14 1210 : const GeochemistryIonicStrength & is_calculator, const GeochemicalDatabaseReader & db)
15 : : GeochemistryActivityCoefficients(),
16 2420 : _numT(db.getTemperatures().size()),
17 1210 : _database_dh_params(db.getDebyeHuckel()),
18 2417 : _database_dh_water((db.getNeutralSpeciesActivity().count("h2o") == 1)
19 2417 : ? db.getNeutralSpeciesActivity().at("h2o")
20 : : GeochemistryNeutralSpeciesActivity()),
21 2417 : _database_dh_neutral((db.getNeutralSpeciesActivity().count("co2") == 1)
22 2417 : ? db.getNeutralSpeciesActivity().at("co2")
23 : : GeochemistryNeutralSpeciesActivity()),
24 1209 : _is_calculator(is_calculator),
25 1209 : _ionic_strength(1.0),
26 1209 : _sqrt_ionic_strength(1.0),
27 1209 : _stoichiometric_ionic_strength(1.0),
28 1209 : _num_basis(0),
29 1209 : _num_eqm(0),
30 : _dh(),
31 1209 : _interp_A(db.getTemperatures(),
32 1209 : (_database_dh_params.adh.size() == _numT) ? _database_dh_params.adh
33 1209 : : std::vector<Real>(_numT, 0.0),
34 1209 : db.getLogKModel()),
35 1209 : _interp_B(db.getTemperatures(),
36 1209 : (_database_dh_params.bdh.size() == _numT) ? _database_dh_params.bdh
37 1209 : : std::vector<Real>(_numT, 0.0),
38 1209 : db.getLogKModel()),
39 1209 : _interp_Bdot(db.getTemperatures(),
40 1209 : (_database_dh_params.bdot.size() == _numT) ? _database_dh_params.bdot
41 1209 : : std::vector<Real>(_numT, 0.0),
42 1209 : db.getLogKModel()),
43 1209 : _interp_a_water(db.getTemperatures(),
44 1209 : (_database_dh_water.a.size() == _numT) ? _database_dh_water.a
45 1210 : : std::vector<Real>(_numT, 0.0),
46 1209 : db.getLogKModel()),
47 1209 : _interp_b_water(db.getTemperatures(),
48 1209 : (_database_dh_water.b.size() == _numT) ? _database_dh_water.b
49 1210 : : std::vector<Real>(_numT, 0.0),
50 1209 : db.getLogKModel()),
51 1209 : _interp_c_water(db.getTemperatures(),
52 1209 : (_database_dh_water.c.size() == _numT) ? _database_dh_water.c
53 1210 : : std::vector<Real>(_numT, 0.0),
54 1209 : db.getLogKModel()),
55 1209 : _interp_d_water(db.getTemperatures(),
56 1209 : (_database_dh_water.d.size() == _numT) ? _database_dh_water.d
57 1210 : : std::vector<Real>(_numT, 0.0),
58 1209 : db.getLogKModel()),
59 1209 : _interp_a_neutral(db.getTemperatures(),
60 1209 : (_database_dh_neutral.a.size() == _numT) ? _database_dh_neutral.a
61 1210 : : std::vector<Real>(_numT, 0.0),
62 1209 : db.getLogKModel()),
63 1209 : _interp_b_neutral(db.getTemperatures(),
64 1209 : (_database_dh_neutral.b.size() == _numT) ? _database_dh_neutral.b
65 1210 : : std::vector<Real>(_numT, 0.0),
66 1209 : db.getLogKModel()),
67 1209 : _interp_c_neutral(db.getTemperatures(),
68 1209 : (_database_dh_neutral.c.size() == _numT) ? _database_dh_neutral.c
69 1210 : : std::vector<Real>(_numT, 0.0),
70 1209 : db.getLogKModel()),
71 1209 : _interp_d_neutral(db.getTemperatures(),
72 1209 : (_database_dh_neutral.d.size() == _numT) ? _database_dh_neutral.d
73 1541 : : std::vector<Real>(_numT, 0.0),
74 2419 : db.getLogKModel())
75 : {
76 1209 : _interp_A.generate();
77 1209 : _interp_B.generate();
78 1209 : _interp_Bdot.generate();
79 1209 : _interp_a_water.generate();
80 1209 : _interp_b_water.generate();
81 1209 : _interp_c_water.generate();
82 1209 : _interp_d_water.generate();
83 1209 : _interp_a_neutral.generate();
84 1209 : _interp_b_neutral.generate();
85 1209 : _interp_c_neutral.generate();
86 1209 : _interp_d_neutral.generate();
87 1209 : }
88 :
89 : void
90 48457 : GeochemistryActivityCoefficientsDebyeHuckel::setInternalParameters(
91 : Real temperature,
92 : const ModelGeochemicalDatabase & mgd,
93 : const std::vector<Real> & basis_species_molality,
94 : const std::vector<Real> & eqm_species_molality,
95 : const std::vector<Real> & kin_species_molality)
96 : {
97 48457 : _ionic_strength = _is_calculator.ionicStrength(
98 : mgd, basis_species_molality, eqm_species_molality, kin_species_molality);
99 48457 : _sqrt_ionic_strength = std::sqrt(_ionic_strength);
100 48457 : _stoichiometric_ionic_strength = _is_calculator.stoichiometricIonicStrength(
101 : mgd, basis_species_molality, eqm_species_molality, kin_species_molality);
102 48457 : _num_basis = mgd.basis_species_index.size();
103 48457 : _num_eqm = mgd.eqm_species_index.size();
104 : // Debye-Huckel base parameters
105 48457 : _dh.A = _interp_A.sample(temperature);
106 48457 : _dh.B = _interp_B.sample(temperature);
107 48457 : _dh.Bdot = _interp_Bdot.sample(temperature);
108 : // water Debye-Huckel
109 48457 : _dh.a_water = _interp_a_water.sample(temperature);
110 48457 : _dh.b_water = _interp_b_water.sample(temperature);
111 48457 : _dh.c_water = _interp_c_water.sample(temperature);
112 48457 : _dh.d_water = _interp_d_water.sample(temperature);
113 : // neutral species Debye-Huckel
114 48457 : _dh.a_neutral = _interp_a_neutral.sample(temperature);
115 48457 : _dh.b_neutral = _interp_b_neutral.sample(temperature);
116 48457 : _dh.c_neutral = _interp_c_neutral.sample(temperature);
117 48457 : _dh.d_neutral = _interp_d_neutral.sample(temperature);
118 48457 : }
119 :
120 : const DebyeHuckelParameters &
121 4 : GeochemistryActivityCoefficientsDebyeHuckel::getDebyeHuckel() const
122 : {
123 4 : return _dh;
124 : }
125 :
126 : Real
127 47627 : GeochemistryActivityCoefficientsDebyeHuckel::waterActivity() const
128 : {
129 95254 : return std::exp(GeochemistryActivityCalculators::lnActivityDHBdotWater(
130 47627 : _stoichiometric_ionic_strength, _dh.A, _dh.a_water, _dh.b_water, _dh.c_water, _dh.d_water));
131 : }
132 :
133 : void
134 47983 : GeochemistryActivityCoefficientsDebyeHuckel::buildActivityCoefficients(
135 : const ModelGeochemicalDatabase & mgd,
136 : std::vector<Real> & basis_activity_coef,
137 : std::vector<Real> & eqm_activity_coef) const
138 : {
139 47983 : basis_activity_coef.resize(_num_basis);
140 47983 : eqm_activity_coef.resize(_num_eqm);
141 :
142 323208 : for (unsigned basis_i = 1; basis_i < _num_basis; ++basis_i) // don't loop over water
143 275225 : if (mgd.basis_species_mineral[basis_i] || mgd.basis_species_gas[basis_i])
144 45717 : continue; // these should never be used
145 229508 : else if (mgd.basis_species_radius[basis_i] == -0.5)
146 : {
147 14401 : basis_activity_coef[basis_i] = std::pow(
148 : 10.0,
149 : GeochemistryActivityCalculators::log10ActCoeffDHBdotNeutral(
150 14401 : _ionic_strength, _dh.a_neutral, _dh.b_neutral, _dh.c_neutral, _dh.d_neutral));
151 : }
152 215107 : else if (mgd.basis_species_radius[basis_i] == -1.0)
153 : {
154 1 : basis_activity_coef[basis_i] =
155 1 : std::pow(10.0,
156 1 : GeochemistryActivityCalculators::log10ActCoeffDHBdotAlternative(_ionic_strength,
157 1 : _dh.Bdot));
158 : }
159 215106 : else if (mgd.basis_species_radius[basis_i] == -1.5)
160 : {
161 168 : basis_activity_coef[basis_i] = 1.0;
162 : }
163 214938 : else if (mgd.basis_species_charge[basis_i] == 0.0)
164 : {
165 7618 : basis_activity_coef[basis_i] = 1.0;
166 : }
167 : else
168 : {
169 207320 : basis_activity_coef[basis_i] = std::pow(
170 : 10.0,
171 : GeochemistryActivityCalculators::log10ActCoeffDHBdot(mgd.basis_species_charge[basis_i],
172 : mgd.basis_species_radius[basis_i],
173 207320 : _sqrt_ionic_strength,
174 207320 : _dh.A,
175 207320 : _dh.B,
176 207320 : _dh.Bdot));
177 : }
178 :
179 1488487 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
180 1440504 : if (mgd.eqm_species_mineral[eqm_j] || mgd.eqm_species_gas[eqm_j])
181 70973 : continue;
182 1369531 : else if (mgd.eqm_species_radius[eqm_j] == -0.5)
183 : {
184 26721 : eqm_activity_coef[eqm_j] = std::pow(
185 : 10.0,
186 : GeochemistryActivityCalculators::log10ActCoeffDHBdotNeutral(
187 26721 : _ionic_strength, _dh.a_neutral, _dh.b_neutral, _dh.c_neutral, _dh.d_neutral));
188 : }
189 1342810 : else if (mgd.eqm_species_radius[eqm_j] == -1.0)
190 : {
191 1 : eqm_activity_coef[eqm_j] =
192 1 : std::pow(10.0,
193 1 : GeochemistryActivityCalculators::log10ActCoeffDHBdotAlternative(_ionic_strength,
194 1 : _dh.Bdot));
195 : }
196 1342809 : else if (mgd.eqm_species_radius[eqm_j] == -1.5)
197 : {
198 18739 : eqm_activity_coef[eqm_j] = 1.0;
199 : }
200 1324070 : else if (mgd.eqm_species_charge[eqm_j] == 0.0)
201 : {
202 442949 : eqm_activity_coef[eqm_j] = 1.0;
203 : }
204 : else
205 : {
206 881121 : eqm_activity_coef[eqm_j] = std::pow(
207 : 10.0,
208 : GeochemistryActivityCalculators::log10ActCoeffDHBdot(mgd.eqm_species_charge[eqm_j],
209 : mgd.eqm_species_radius[eqm_j],
210 881121 : _sqrt_ionic_strength,
211 881121 : _dh.A,
212 881121 : _dh.B,
213 881121 : _dh.Bdot));
214 : }
215 47983 : }
216 :
217 : Real
218 2 : GeochemistryActivityCoefficientsDebyeHuckel::getIonicStrength() const
219 : {
220 2 : return _ionic_strength;
221 : }
222 :
223 : Real
224 2 : GeochemistryActivityCoefficientsDebyeHuckel::getStoichiometricIonicStrength() const
225 : {
226 2 : return _stoichiometric_ionic_strength;
227 : }
|
