Plugins
Loading...
Searching...
No Matches
smile_chem.hpp
Go to the documentation of this file.
1/*
2 * smile_chem.hpp
3 */
4
5#ifndef _SMILE_CHEM_HPP_
6#define _SMILE_CHEM_HPP_
7
8#include <cmath>
9#include <vector>
10
11#include "common/common.h"
15
16#include "smile_chem_eq.hpp"
17#include "smile_cloud.hpp"
18
19namespace ug {
20namespace smile {
21
26template <typename TGridFunction>
27static void SmartChem
28(
31 const char* nuclides,
32 Cloud *cloud,
33 number porosity,
34 number rock_density
35)
36{
37 UG_LOG (">>> Start SmartChem <<<\n");
38
39 // get function ids by name (agents, fixed names)
40 const size_t fct_H = spGF->fct_id_by_name("H");
41 if (fct_H > spGF->num_fct()) UG_THROW("SmileChem: Component 'H' not found.");
42 const size_t fct_OH = spGF->fct_id_by_name("OH");
43 if (fct_OH > spGF->num_fct()) UG_THROW("SmileChem: Component 'OH' not found.");
44 const size_t fct_S = spGF->fct_id_by_name("c");
45 if (fct_S > spGF->num_fct()) UG_THROW("SmileChem: Component 'c' not found.");
46 const size_t fct_Ca = spGF->fct_id_by_name("Ca");
47 if (fct_Ca > spGF->num_fct()) UG_THROW("SmileChem: Component 'Ca' not found.");
48 const size_t fct_DIC = spGF->fct_id_by_name("DIC");
49 if (fct_DIC > spGF->num_fct()) UG_THROW("SmileChem: Component 'DIC' not found.");
50 const size_t fct_SO4 = spGF->fct_id_by_name("SO4");
51 if (fct_SO4 > spGF->num_fct()) UG_THROW("SmileChem: Component 'SO4' not found.");
52 const size_t fct_Al = spGF->fct_id_by_name("Al");
53 if (fct_Al > spGF->num_fct()) UG_THROW("SmileChem: Component 'Al' not found.");
54 const size_t fct_Cc = spGF2->fct_id_by_name("c");
55 if (fct_Cc > spGF2->num_fct()) UG_THROW("SmileChem: Component 'Calcite' not found.");
56
57 // get function ids by name (nuclides)
58 size_t no_nuc = 0;
59 size_t fct_nuc[CLOUD_MAX_NUCLIDES];
60 int ind_nuc[CLOUD_MAX_NUCLIDES];
61 while (1) {
62 char name[32];
63 int r, n;
64 r = sscanf(nuclides, "%s %d %n", name, ind_nuc + no_nuc, &n);
65 if (r <= 0) break;
66 if (r != 2) UG_THROW("SmartChem: List of nuclides and indices malformed.");
67 fct_nuc[no_nuc] = spGF->fct_id_by_name(name);
68 if (fct_nuc[no_nuc] > spGF->num_fct()) UG_THROW("SmartChem: Component '" << name << "' not found.");
69 no_nuc++;
70 nuclides +=n;
71 }
72 if (no_nuc == 0) UG_THROW("SmartChem: List of nuclides and indices malformed.");
73
74 // Loop the vertices
75 typedef typename TGridFunction::template traits<Vertex>::const_iterator t_vert_iterator;
76 t_vert_iterator iter = spGF->template begin<Vertex> ();
77 t_vert_iterator iterEnd = spGF->template end<Vertex> ();
78 for (; iter != iterEnd; iter++)
79 {
80 Vertex * pVertex = *iter;
81 /*-- For debugging only: --*
82 bool is_my_vertex = false;
83 typedef typename TGridFunction::domain_type domain_type;
84 static const int dim = domain_type::dim;
85 typedef typename TGridFunction::domain_type::position_type coord_type; // actually, MathVector<dim>
86 std::vector<coord_type> vrt_coord(1);
87 CollectCornerCoordinates (vrt_coord, *pVertex, * (spGF->domain ()));
88 {
89 MathVector<dim> bb_low, bb_high;
90
91 // interval # 1
92 bb_low[0] = 0.15; bb_low[1] = 0.75; // bb_low[2] = 0;
93 bb_high[0] = 0.25; bb_high[1] = 0.85; // bb_high[2] = 0;
94 is_my_vertex |= VecIsInBB (vrt_coord[0], bb_low, bb_high);
95
96 // interval # 2
97 //bb_low[0] = 2; bb_low[1] = 0.4; // bb_low[2] = 0;
98 //bb_high[0] = 2.2; bb_high[1] = 0.4; // bb_high[2] = 0;
99 //is_my_vertex |= VecIsInBB (vrt_coord[0], bb_low, bb_high);
100 }
101 *--*/
102
103 // get transported concentrations of agents
104 std::vector<DoFIndex> ind_H(1), ind_OH(1), ind_S(1), ind_Ca(1),
105 ind_DIC(1), ind_Cc(1), ind_SO4(1), ind_Al(1);
106 spGF ->inner_dof_indices (pVertex, fct_H, ind_H);
107 spGF ->inner_dof_indices (pVertex, fct_OH, ind_OH);
108 spGF ->inner_dof_indices (pVertex, fct_S, ind_S);
109 spGF ->inner_dof_indices (pVertex, fct_Ca, ind_Ca);
110 spGF ->inner_dof_indices (pVertex, fct_DIC, ind_DIC);
111 spGF ->inner_dof_indices (pVertex, fct_SO4, ind_SO4);
112 spGF ->inner_dof_indices (pVertex, fct_Al, ind_Al);
113 spGF2->inner_dof_indices (pVertex, fct_Cc, ind_Cc);
114 double H_T = DoFRef(*spGF, ind_H [0]); // all in mol/m^3
115 double OH_T = DoFRef(*spGF, ind_OH [0]);
116 double s = DoFRef(*spGF, ind_S [0]);
117 double Ca_T = DoFRef(*spGF, ind_Ca [0]);
118 double DIC_T = DoFRef(*spGF, ind_DIC[0]);
119 double SO4_T = DoFRef(*spGF, ind_SO4[0]);
120 double Al_T = DoFRef(*spGF, ind_Al [0]);
121 double Cc_T = DoFRef(*spGF2,ind_Cc [0]);
122
123
124 /*-- For debugging only --*
125 if (is_my_vertex)
126 {
127 UG_LOG ("\n\n>>>Vertex @ (" << vrt_coord[0] << ")\n");
128 UG_LOG (" Input values: \n");
129 UG_LOG ("\n H_T = " << H_T << " OH_T = " << OH_T << " DIC_T = " << DIC_T << " \n");
130 UG_LOG ("\n s = " << s << " Ca_T = " << Ca_T << " Al_T = " << Al_T << " \n");
131 UG_LOG ("\n SO4_T = " << SO4_T << " \n");
132 }
133 *--*/
134
135 /* chlorid, ion strength & uncorrected pH */
136 double I, D;
137 t_gamma gamma; // this is a local copy
138 s *= 5428.78; // convert brine mass fraction to concentration
139 Compute_I_D(1e-3*s, 1e-3*Ca_T, 1e-3*Al_T, 1e-3*DIC_T, 1e-3*SO4_T, &I, &D, &gamma);
140 double pH_T = pH_of_sHD(1e-3*s, 1e-3*H_T, D);
141
142 /*-- For debugging only --*
143 if (is_my_vertex)
144 {
145 UG_LOG ("\n\n Rescaled s = " << s << " \n");
146 UG_LOG (" After 1st Compute_I_D: \n");
147 UG_LOG ("\n I = " << I << " D = " << D << " pH_T = " << pH_T << " \n");
148 }
149 *--*/
150
151 /* scale all nuclide concentrations with R_old */
152 for (size_t i = 0; i < no_nuc; i++) {
153 std::vector<DoFIndex> ind(1);
154 spGF->inner_dof_indices (pVertex, fct_nuc[i], ind);
155 double c = DoFRef(*spGF, ind[0]);
156 double x[6];
157 x[0] = pH_T;
158 x[1] = DIC_T;
159 x[2] = I;
160 x[3] = Ca_T;
161 x[4] = Al_T;
162 x[5] = SO4_T;
163 double Kd = cloud->interpolate(x, ind_nuc[i]);
164 double R = 1.0 + (1.0-porosity) / porosity * Kd * rock_density;
165 DoFRef(*spGF, ind[0]) = R * c;
166
167 /*-- For debugging only --*
168 if (is_my_vertex && i == 3)
169 {
170 UG_LOG ("\n\n Before equilibration: \n");
171 UG_LOG ("\n c = " << c << " Kd = " << Kd << " R = " << R << " \n");
172 }
173 *--*/
174 }
175
176 H_T *= 1e-3; // scale down to mol/l
177 OH_T *= 1e-3; // scale down to mol/l
178 Al_T *= 1e-3; // scale down to mol/l
179 Ca_T *= 1e-3; // scale down to mol/l
180 DIC_T *= 1e-3; // scale down to mol/l
181 Cc_T *= 1e-3; // scale down to mol/l
182 SO4_T *= 1e-3; // scale down to mol/l
183 s *= 1e-3; // scale down to mol/l
184
185 /*-- For debugging only --*
186 if (is_my_vertex)
187 {
188 UG_LOG ("\n\n SmartEquilibrate args before the call: \n");
189 UG_LOG (" H_T = " << H_T << ',');
190 UG_LOG (" OH_T = " << OH_T << ',');
191 UG_LOG (" s = " << s << ',');
192 UG_LOG (" Ca_T = " << Ca_T << ',');
193 UG_LOG (" DIC_T = " << DIC_T << ',');
194 UG_LOG (" SO4_T = " << SO4_T << ',');
195 UG_LOG (" Al_T = " << Al_T << ',');
196 UG_LOG (" Cc_T = " << Cc_T << ',');
197 UG_LOG ("\n");
198 }
199 *--*/
200
201 /* Compute the chemical equilibrium */
202 int err_code = SmartEquilibrate (H_T, OH_T, s, Ca_T, DIC_T, SO4_T, Al_T, Cc_T);
203
204 /*-- For debugging only --*
205 if (is_my_vertex)
206 {
207 UG_LOG ("\n\n SmartEquilibrate args after the call: \n");
208 UG_LOG (" H_T = " << H_T << ',');
209 UG_LOG (" OH_T = " << OH_T << ',');
210 UG_LOG (" s = " << s << ',');
211 UG_LOG (" Ca_T = " << Ca_T << ',');
212 UG_LOG (" DIC_T = " << DIC_T << ',');
213 UG_LOG (" SO4_T = " << SO4_T << ',');
214 UG_LOG (" Al_T = " << Al_T << ',');
215 UG_LOG (" Cc_T = " << Cc_T << ',');
216 UG_LOG ("\n");
217 }
218 *--*/
219
220 if (err_code != 0)
221 {
222 typedef typename TGridFunction::domain_type domain_type;
223 const typename domain_type::position_accessor_type & aaPos
224 = spGF->domain()->position_accessor();
225 typename domain_type::position_type pos = aaPos [pVertex];
226 UG_LOG ("!!! Error in SmartEquilibrate at vertex grid data idx" << pVertex->grid_data_index () << " # (" << pos[0]);
227 for (size_t i = 1; i < domain_type::dim; i++) UG_LOG (", " << pos[i]);
228 UG_LOG (")\n");
229 if (err_code > 0)
230 UG_THROW ("!!! Fatal error in SmartEquilibrate!");
231 }
232
233 /* set Ca, DIC, Calcite */
234 DoFRef(*spGF2, ind_Cc [0]) = 1e3*Cc_T; // scale up to mol/m^3
235 DoFRef(*spGF, ind_Ca [0]) = 1e3*Ca_T; // scale up to mol/m^3
236 DoFRef(*spGF, ind_DIC[0]) = 1e3*DIC_T;// scale up to mol/m^3
237
238 /* set H, OH, Al */
239 DoFRef(*spGF, ind_Al [0]) = 1e3*Al_T; // scale up to mol/m^3
240 DoFRef(*spGF, ind_H [0]) = 1e3*H_T; // scale up to mol/m^3
241 DoFRef(*spGF, ind_OH[0]) = 1e3*OH_T; // scale up to mol/m^3
242
243 /* scale all nuclide concentrations with 1/R_new */
244 Compute_I_pH (s, H_T, Ca_T, Al_T, DIC_T, SO4_T, I, pH_T);
245 for (size_t i = 0; i < no_nuc; i++)
246 {
247 std::vector<DoFIndex> ind(1);
248 spGF->inner_dof_indices (pVertex, fct_nuc[i], ind);
249 double c = DoFRef(*spGF, ind[0]);
250 double x[6];
251 x[0] = pH_T;
252 x[1] = 1e3*DIC_T; // scale up to mol/m^3
253 x[2] = I;
254 x[3] = 1e3*Ca_T; // scale up to mol/m^3
255 x[4] = 1e3*Al_T; // scale up to mol/m^3
256 x[5] = 1e3*SO4_T;// scale up to mol/m^3
257 /*-- For debugging only --*
258 if (is_my_vertex && i == 3)
259 {
260 UG_LOG ("\n\n Preparation of the rescaling: \n");
261 UG_LOG ("\n pH_T = " << x[0] << " DIC_T = " << x[1] << " I = " << x[2] << " \n");
262 UG_LOG ("\n Ca_T = " << x[3] << " Al_T = " << x[4] << " SO4_T = " << x[5] << " \n");
263 }
264 *--*/
265 double Kd = cloud->interpolate(x, ind_nuc[i]);
266 double R = 1.0 + (1.0-porosity) / porosity * Kd * rock_density;
267 DoFRef(*spGF, ind[0]) = c / R;
268
269 /*-- For debugging only --*
270 if (is_my_vertex && i == 3)
271 {
272 UG_LOG ("\n\n After equilibration: \n");
273 UG_LOG ("\n c = " << c << " Kd = " << Kd << " R = " << R << " \n");
274 }
275 *--*/
276 }
277
278 /*-- For debugging only --*
279 if (is_my_vertex)
280 {
281 UG_LOG ("\n\n Output values: \n");
282 UG_LOG ("\n H_T = " << H_T << " OH_T = " << OH_T << " DIC_T = " << DIC_T << " \n");
283 UG_LOG ("\n s = " << s << " Ca_T = " << Ca_T << " Al_T = " << Al_T << " \n");
284 UG_LOG ("\n SO4_T = " << SO4_T << " \n");
285 }
286 *--*/
287 }
288 UG_LOG ("<<< End SmartChem >>>\n");
289}
290
294template <typename TGridFunction>
296(
299)
300{
301 UG_LOG (">>> Start SmartChemEquilibrate <<<\n");
302
303 // get function ids by name (agents, fixed names)
304 const size_t fct_H = spGF->fct_id_by_name("H");
305 if (fct_H > spGF->num_fct()) UG_THROW("SmileChem: Component 'H' not found.");
306 const size_t fct_OH = spGF->fct_id_by_name("OH");
307 if (fct_OH > spGF->num_fct()) UG_THROW("SmileChem: Component 'OH' not found.");
308 const size_t fct_S = spGF->fct_id_by_name("c");
309 if (fct_S > spGF->num_fct()) UG_THROW("SmileChem: Component 'c' not found.");
310 const size_t fct_Ca = spGF->fct_id_by_name("Ca");
311 if (fct_Ca > spGF->num_fct()) UG_THROW("SmileChem: Component 'Ca' not found.");
312 const size_t fct_DIC = spGF->fct_id_by_name("DIC");
313 if (fct_DIC > spGF->num_fct()) UG_THROW("SmileChem: Component 'DIC' not found.");
314 const size_t fct_SO4 = spGF->fct_id_by_name("SO4");
315 if (fct_SO4 > spGF->num_fct()) UG_THROW("SmileChem: Component 'SO4' not found.");
316 const size_t fct_Al = spGF->fct_id_by_name("Al");
317 if (fct_Al > spGF->num_fct()) UG_THROW("SmileChem: Component 'Al' not found.");
318 const size_t fct_Cc = spGF2->fct_id_by_name("c");
319 if (fct_Cc > spGF2->num_fct()) UG_THROW("SmileChem: Component 'Calcite' not found.");
320
321 // Loop the vertices
322 typedef typename TGridFunction::template traits<Vertex>::const_iterator t_vert_iterator;
323 t_vert_iterator iter = spGF->template begin<Vertex> ();
324 t_vert_iterator iterEnd = spGF->template end<Vertex> ();
325 for (; iter != iterEnd; iter++)
326 {
327 Vertex * pVertex = *iter;
328 /*-- For debugging only: --*
329 bool is_my_vertex = false;
330 typedef typename TGridFunction::domain_type domain_type;
331 static const int dim = domain_type::dim;
332 typedef typename TGridFunction::domain_type::position_type coord_type; // actually, MathVector<dim>
333 std::vector<coord_type> vrt_coord(1);
334 CollectCornerCoordinates (vrt_coord, *pVertex, * (spGF->domain ()));
335 {
336 MathVector<dim> bb_low, bb_high;
337
338 // interval # 1
339 bb_low[0] = 0.15; bb_low[1] = 0.75; // bb_low[2] = 0;
340 bb_high[0] = 0.25; bb_high[1] = 0.85; // bb_high[2] = 0;
341 is_my_vertex |= VecIsInBB (vrt_coord[0], bb_low, bb_high);
342
343 // interval # 2
344 //bb_low[0] = 2; bb_low[1] = 0.4; // bb_low[2] = 0;
345 //bb_high[0] = 2.2; bb_high[1] = 0.4; // bb_high[2] = 0;
346 //is_my_vertex |= VecIsInBB (vrt_coord[0], bb_low, bb_high);
347 }
348 *--*/
349
350 // get transported concentrations of agents
351 std::vector<DoFIndex> ind_H(1), ind_OH(1), ind_S(1), ind_Ca(1),
352 ind_DIC(1), ind_Cc(1), ind_SO4(1), ind_Al(1);
353 spGF ->inner_dof_indices (pVertex, fct_H, ind_H);
354 spGF ->inner_dof_indices (pVertex, fct_OH, ind_OH);
355 spGF ->inner_dof_indices (pVertex, fct_S, ind_S);
356 spGF ->inner_dof_indices (pVertex, fct_Ca, ind_Ca);
357 spGF ->inner_dof_indices (pVertex, fct_DIC, ind_DIC);
358 spGF ->inner_dof_indices (pVertex, fct_SO4, ind_SO4);
359 spGF ->inner_dof_indices (pVertex, fct_Al, ind_Al);
360 spGF2->inner_dof_indices (pVertex, fct_Cc, ind_Cc);
361 double H_T = DoFRef(*spGF, ind_H [0]); // all in mol/m^3
362 double OH_T = DoFRef(*spGF, ind_OH [0]);
363 double s = DoFRef(*spGF, ind_S [0]);
364 double Ca_T = DoFRef(*spGF, ind_Ca [0]);
365 double DIC_T = DoFRef(*spGF, ind_DIC[0]);
366 double SO4_T = DoFRef(*spGF, ind_SO4[0]);
367 double Al_T = DoFRef(*spGF, ind_Al [0]);
368 double Cc_T = DoFRef(*spGF2,ind_Cc [0]);
369
370 H_T *= 1e-3; // scale down to mol/l
371 Al_T *= 1e-3; // scale down to mol/l
372 Ca_T *= 1e-3; // scale down to mol/l
373 DIC_T *= 1e-3; // scale down to mol/l
374 Cc_T *= 1e-3; // scale down to mol/l
375 SO4_T *= 1e-3; // scale down to mol/l
376 s *= 1e-3 * NaCl_sat; // scale the mass fraction to molar concentraton in mol/l
377
378 /*-- For debugging only --*
379 if (is_my_vertex)
380 {
381 UG_LOG ("\n\n SmartEquilibrate args before the call: \n");
382 UG_LOG (" H_T = " << H_T << ',');
383 UG_LOG (" OH_T = " << OH_T << ',');
384 UG_LOG (" s = " << s << ',');
385 UG_LOG (" Ca_T = " << Ca_T << ',');
386 UG_LOG (" DIC_T = " << DIC_T << ',');
387 UG_LOG (" SO4_T = " << SO4_T << ',');
388 UG_LOG (" Al_T = " << Al_T << ',');
389 UG_LOG (" Cc_T = " << Cc_T << ',');
390 UG_LOG ("\n");
391 }
392 *--*/
393
394 /* Compute the chemical equilibrium */
395 int err_code = SmartEquilibrate (H_T, OH_T, s, Ca_T, DIC_T, SO4_T, Al_T, Cc_T);
396
397 /*-- For debugging only --*
398 if (is_my_vertex)
399 {
400 UG_LOG ("\n\n SmartEquilibrate args after the call: \n");
401 UG_LOG (" H_T = " << H_T << ',');
402 UG_LOG (" OH_T = " << OH_T << ',');
403 UG_LOG (" s = " << s << ',');
404 UG_LOG (" Ca_T = " << Ca_T << ',');
405 UG_LOG (" DIC_T = " << DIC_T << ',');
406 UG_LOG (" SO4_T = " << SO4_T << ',');
407 UG_LOG (" Al_T = " << Al_T << ',');
408 UG_LOG (" Cc_T = " << Cc_T << ',');
409 UG_LOG ("\n");
410 }
411 *--*/
412
413 if (err_code != 0)
414 {
415 typedef typename TGridFunction::domain_type domain_type;
416 const typename domain_type::position_accessor_type & aaPos
417 = spGF->domain()->position_accessor();
418 typename domain_type::position_type pos = aaPos [pVertex];
419 UG_LOG ("!!! Error in SmartEquilibrate at vertex grid data idx" << pVertex->grid_data_index () << " # (" << pos[0]);
420 for (size_t i = 1; i < domain_type::dim; i++) UG_LOG (", " << pos[i]);
421 UG_LOG (")\n");
422 if (err_code > 0)
423 UG_THROW ("!!! Fatal error in SmartEquilibrate!");
424 }
425
426 /* set Ca, DIC, Calcite */
427 DoFRef(*spGF2, ind_Cc [0]) = 1e3*Cc_T;
428 DoFRef(*spGF, ind_Ca [0]) = 1e3*Ca_T;
429 DoFRef(*spGF, ind_DIC[0]) = 1e3*DIC_T;
430
431 /* set H, OH, Al */
432 DoFRef(*spGF, ind_Al [0]) = 1e3*Al_T; // scale up to mol/m^3
433 DoFRef(*spGF, ind_H [0]) = 1e3*H_T;
434 DoFRef(*spGF, ind_OH[0]) = 1e3*OH_T;
435 }
436 UG_LOG ("<<< End SmartChemEquilibrate >>>\n");
437}
438
439} // end namespace smile
440} // end namespace ug
441
442#endif
parameterString s
Definition Biogas.lua:2
uint grid_data_index() const
Definition smile_cloud.hpp:15
double interpolate(double *x, int k)
Definition smile_cloud.cpp:68
#define UG_THROW(msg)
#define UG_LOG(msg)
double number
void Compute_I_D(double s, double Ca, double Al, double DIC, double SO4, double *I, double *D, t_gamma *gamma)
Definition smile_chem_eq.cpp:36
void Compute_I_pH(number s, number H, number Ca, number Al, number DIC, number SO4, number &I, number &pH)
Definition smile_chem_eq.cpp:24
static void SmartChem(SmartPtr< TGridFunction > spGF, SmartPtr< TGridFunction > spGF2, const char *nuclides, Cloud *cloud, number porosity, number rock_density)
Definition smile_chem.hpp:28
static void SmartChemEquilibrate(SmartPtr< TGridFunction > spGF, SmartPtr< TGridFunction > spGF2)
Definition smile_chem.hpp:296
number pH_of_sHD(number s, number H, number D)
Definition smile_chem_eq.cpp:17
int SmartEquilibrate(double &H_T, double &OH_T, double &s, double &Ca_T, double &DIC_T, double &SO4_T, double &Al_T, double &Cc_T)
Definition smile_chem_eq.cpp:139
const double NaCl_sat
Definition smile_chem_eq.hpp:55
const number & DoFRef(const TMatrix &mat, const DoFIndex &iInd, const DoFIndex &jInd)
#define CLOUD_MAX_NUCLIDES
Definition smile_cloud.hpp:9
Definition smile_chem_eq.hpp:67