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