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 1266 : GeochemistryActivityCoefficientsDebyeHuckel::GeochemistryActivityCoefficientsDebyeHuckel(
14 1266 : const GeochemistryIonicStrength & is_calculator, const GeochemicalDatabaseReader & db)
15 : : GeochemistryActivityCoefficients(),
16 2532 : _numT(db.getTemperatures().size()),
17 1266 : _database_dh_params(db.getDebyeHuckel()),
18 2529 : _database_dh_water((db.getNeutralSpeciesActivity().count("h2o") == 1)
19 2529 : ? db.getNeutralSpeciesActivity().at("h2o")
20 : : GeochemistryNeutralSpeciesActivity()),
21 2529 : _database_dh_neutral((db.getNeutralSpeciesActivity().count("co2") == 1)
22 2529 : ? db.getNeutralSpeciesActivity().at("co2")
23 : : GeochemistryNeutralSpeciesActivity()),
24 1265 : _is_calculator(is_calculator),
25 1265 : _ionic_strength(1.0),
26 1265 : _sqrt_ionic_strength(1.0),
27 1265 : _stoichiometric_ionic_strength(1.0),
28 1265 : _num_basis(0),
29 1265 : _num_eqm(0),
30 : _dh(),
31 1265 : _interp_A(db.getTemperatures(),
32 1265 : (_database_dh_params.adh.size() == _numT) ? _database_dh_params.adh
33 1265 : : std::vector<Real>(_numT, 0.0),
34 1265 : db.getLogKModel()),
35 1265 : _interp_B(db.getTemperatures(),
36 1265 : (_database_dh_params.bdh.size() == _numT) ? _database_dh_params.bdh
37 1265 : : std::vector<Real>(_numT, 0.0),
38 1265 : db.getLogKModel()),
39 1265 : _interp_Bdot(db.getTemperatures(),
40 1265 : (_database_dh_params.bdot.size() == _numT) ? _database_dh_params.bdot
41 1265 : : std::vector<Real>(_numT, 0.0),
42 1265 : db.getLogKModel()),
43 1265 : _interp_a_water(db.getTemperatures(),
44 1265 : (_database_dh_water.a.size() == _numT) ? _database_dh_water.a
45 1266 : : std::vector<Real>(_numT, 0.0),
46 1265 : db.getLogKModel()),
47 1265 : _interp_b_water(db.getTemperatures(),
48 1265 : (_database_dh_water.b.size() == _numT) ? _database_dh_water.b
49 1266 : : std::vector<Real>(_numT, 0.0),
50 1265 : db.getLogKModel()),
51 1265 : _interp_c_water(db.getTemperatures(),
52 1265 : (_database_dh_water.c.size() == _numT) ? _database_dh_water.c
53 1266 : : std::vector<Real>(_numT, 0.0),
54 1265 : db.getLogKModel()),
55 1265 : _interp_d_water(db.getTemperatures(),
56 1265 : (_database_dh_water.d.size() == _numT) ? _database_dh_water.d
57 1266 : : std::vector<Real>(_numT, 0.0),
58 1265 : db.getLogKModel()),
59 1265 : _interp_a_neutral(db.getTemperatures(),
60 1265 : (_database_dh_neutral.a.size() == _numT) ? _database_dh_neutral.a
61 1266 : : std::vector<Real>(_numT, 0.0),
62 1265 : db.getLogKModel()),
63 1265 : _interp_b_neutral(db.getTemperatures(),
64 1265 : (_database_dh_neutral.b.size() == _numT) ? _database_dh_neutral.b
65 1266 : : std::vector<Real>(_numT, 0.0),
66 1265 : db.getLogKModel()),
67 1265 : _interp_c_neutral(db.getTemperatures(),
68 1265 : (_database_dh_neutral.c.size() == _numT) ? _database_dh_neutral.c
69 1266 : : std::vector<Real>(_numT, 0.0),
70 1265 : db.getLogKModel()),
71 1265 : _interp_d_neutral(db.getTemperatures(),
72 1265 : (_database_dh_neutral.d.size() == _numT) ? _database_dh_neutral.d
73 1597 : : std::vector<Real>(_numT, 0.0),
74 2531 : db.getLogKModel())
75 : {
76 1265 : _interp_A.generate();
77 1265 : _interp_B.generate();
78 1265 : _interp_Bdot.generate();
79 1265 : _interp_a_water.generate();
80 1265 : _interp_b_water.generate();
81 1265 : _interp_c_water.generate();
82 1265 : _interp_d_water.generate();
83 1265 : _interp_a_neutral.generate();
84 1265 : _interp_b_neutral.generate();
85 1265 : _interp_c_neutral.generate();
86 1265 : _interp_d_neutral.generate();
87 1265 : }
88 :
89 : void
90 54953 : 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 54953 : _ionic_strength = _is_calculator.ionicStrength(
98 : mgd, basis_species_molality, eqm_species_molality, kin_species_molality);
99 54953 : _sqrt_ionic_strength = std::sqrt(_ionic_strength);
100 54953 : _stoichiometric_ionic_strength = _is_calculator.stoichiometricIonicStrength(
101 : mgd, basis_species_molality, eqm_species_molality, kin_species_molality);
102 54953 : _num_basis = mgd.basis_species_index.size();
103 54953 : _num_eqm = mgd.eqm_species_index.size();
104 : // Debye-Huckel base parameters
105 54953 : _dh.A = _interp_A.sample(temperature);
106 54953 : _dh.B = _interp_B.sample(temperature);
107 54953 : _dh.Bdot = _interp_Bdot.sample(temperature);
108 : // water Debye-Huckel
109 54953 : _dh.a_water = _interp_a_water.sample(temperature);
110 54953 : _dh.b_water = _interp_b_water.sample(temperature);
111 54953 : _dh.c_water = _interp_c_water.sample(temperature);
112 54953 : _dh.d_water = _interp_d_water.sample(temperature);
113 : // neutral species Debye-Huckel
114 54953 : _dh.a_neutral = _interp_a_neutral.sample(temperature);
115 54953 : _dh.b_neutral = _interp_b_neutral.sample(temperature);
116 54953 : _dh.c_neutral = _interp_c_neutral.sample(temperature);
117 54953 : _dh.d_neutral = _interp_d_neutral.sample(temperature);
118 54953 : }
119 :
120 : const DebyeHuckelParameters &
121 4 : GeochemistryActivityCoefficientsDebyeHuckel::getDebyeHuckel() const
122 : {
123 4 : return _dh;
124 : }
125 :
126 : Real
127 54030 : GeochemistryActivityCoefficientsDebyeHuckel::waterActivity() const
128 : {
129 108060 : return std::exp(GeochemistryActivityCalculators::lnActivityDHBdotWater(
130 54030 : _stoichiometric_ionic_strength, _dh.A, _dh.a_water, _dh.b_water, _dh.c_water, _dh.d_water));
131 : }
132 :
133 : void
134 54386 : GeochemistryActivityCoefficientsDebyeHuckel::buildActivityCoefficients(
135 : const ModelGeochemicalDatabase & mgd,
136 : std::vector<Real> & basis_activity_coef,
137 : std::vector<Real> & eqm_activity_coef) const
138 : {
139 54386 : basis_activity_coef.resize(_num_basis);
140 54386 : eqm_activity_coef.resize(_num_eqm);
141 :
142 362282 : for (unsigned basis_i = 1; basis_i < _num_basis; ++basis_i) // don't loop over water
143 307896 : if (mgd.basis_species_mineral[basis_i] || mgd.basis_species_gas[basis_i])
144 53654 : continue; // these should never be used
145 254242 : else if (mgd.basis_species_radius[basis_i] == -0.5)
146 : {
147 15396 : basis_activity_coef[basis_i] = std::pow(
148 : 10.0,
149 : GeochemistryActivityCalculators::log10ActCoeffDHBdotNeutral(
150 15396 : _ionic_strength, _dh.a_neutral, _dh.b_neutral, _dh.c_neutral, _dh.d_neutral));
151 : }
152 238846 : 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 238845 : else if (mgd.basis_species_radius[basis_i] == -1.5)
160 : {
161 193 : basis_activity_coef[basis_i] = 1.0;
162 : }
163 238652 : else if (mgd.basis_species_charge[basis_i] == 0.0)
164 : {
165 8646 : basis_activity_coef[basis_i] = 1.0;
166 : }
167 : else
168 : {
169 230006 : 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 230006 : _sqrt_ionic_strength,
174 230006 : _dh.A,
175 230006 : _dh.B,
176 230006 : _dh.Bdot));
177 : }
178 :
179 1638061 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
180 1583675 : if (mgd.eqm_species_mineral[eqm_j] || mgd.eqm_species_gas[eqm_j])
181 77467 : continue;
182 1506208 : else if (mgd.eqm_species_radius[eqm_j] == -0.5)
183 : {
184 29437 : eqm_activity_coef[eqm_j] = std::pow(
185 : 10.0,
186 : GeochemistryActivityCalculators::log10ActCoeffDHBdotNeutral(
187 29437 : _ionic_strength, _dh.a_neutral, _dh.b_neutral, _dh.c_neutral, _dh.d_neutral));
188 : }
189 1476771 : 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 1476770 : else if (mgd.eqm_species_radius[eqm_j] == -1.5)
197 : {
198 20627 : eqm_activity_coef[eqm_j] = 1.0;
199 : }
200 1456143 : else if (mgd.eqm_species_charge[eqm_j] == 0.0)
201 : {
202 488371 : eqm_activity_coef[eqm_j] = 1.0;
203 : }
204 : else
205 : {
206 967772 : 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 967772 : _sqrt_ionic_strength,
211 967772 : _dh.A,
212 967772 : _dh.B,
213 967772 : _dh.Bdot));
214 : }
215 54386 : }
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 : }
|
