37 UG_LOG (
">>> Start SmartChem <<<\n");
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.");
64 r = sscanf(nuclides,
"%s %d %n", name, ind_nuc + no_nuc, &n);
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.");
72 if (no_nuc == 0)
UG_THROW(
"SmartChem: List of nuclides and indices malformed.");
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++)
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]);
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]);
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);
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]);
164 double R = 1.0 + (1.0-porosity) / porosity * Kd * rock_density;
165 DoFRef(*spGF, ind[0]) = R * c;
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]);
230 UG_THROW (
"!!! Fatal error in SmartEquilibrate!");
234 DoFRef(*spGF2, ind_Cc [0]) = 1e3*Cc_T;
235 DoFRef(*spGF, ind_Ca [0]) = 1e3*Ca_T;
236 DoFRef(*spGF, ind_DIC[0]) = 1e3*DIC_T;
239 DoFRef(*spGF, ind_Al [0]) = 1e3*Al_T;
240 DoFRef(*spGF, ind_H [0]) = 1e3*H_T;
241 DoFRef(*spGF, ind_OH[0]) = 1e3*OH_T;
245 for (
size_t i = 0; i < no_nuc; i++)
247 std::vector<DoFIndex> ind(1);
248 spGF->inner_dof_indices (pVertex, fct_nuc[i], ind);
249 double c =
DoFRef(*spGF, ind[0]);
266 double R = 1.0 + (1.0-porosity) / porosity * Kd * rock_density;
267 DoFRef(*spGF, ind[0]) = c / R;
288 UG_LOG (
"<<< End SmartChem >>>\n");
301 UG_LOG (
">>> Start SmartChemEquilibrate <<<\n");
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.");
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++)
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]);
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]);
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]);
423 UG_THROW (
"!!! Fatal error in SmartEquilibrate!");
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;
432 DoFRef(*spGF, ind_Al [0]) = 1e3*Al_T;
433 DoFRef(*spGF, ind_H [0]) = 1e3*H_T;
434 DoFRef(*spGF, ind_OH[0]) = 1e3*OH_T;
436 UG_LOG (
"<<< End SmartChemEquilibrate >>>\n");