Flow123d  master-f44eb46
sorption_base.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file sorption_base.cc
15  * @brief
16  */
17 
23 #include "reaction/isotherm.hh"
25 
26 #include "system/system.hh"
27 #include "system/sys_profiler.hh"
28 
29 #include "la/distribution.hh"
30 #include "mesh/mesh.h"
31 #include "mesh/elements.h"
32 #include "mesh/accessors.hh"
33 #include "input/type_selection.hh"
34 
35 #include "fields/field_set.hh"
36 #include "fields/field_fe.hh"
37 
38 using namespace Input::Type;
39 
41  return Selection("SorptionType")
42  .add_value(Isotherm::none,"none", "No sorption considered.")
43  .add_value(Isotherm::linear, "linear",
44  "Linear isotherm runs the concentration exchange between liquid and solid.")
45  .add_value(Isotherm::langmuir, "langmuir",
46  "Langmuir isotherm runs the concentration exchange between liquid and solid.")
47  .add_value(Isotherm::freundlich, "freundlich",
48  "Freundlich isotherm runs the concentration exchange between liquid and solid.")
49  .close();
50 }
51 
52 
53 
55  return Record("SorptionAux", "AUXILIARY RECORD. Should not be directly part of the input tree.")
56  .declare_key("substances", Array(String(),1), Default::obligatory(),
57  "Names of the substances that take part in the sorption model.")
58  .declare_key("solvent_density", Double(0.0), Default("1.0"),
59  "Density of the solvent.")
60  .declare_key("substeps", Integer(1), Default("1000"),
61  "Number of equidistant substeps, molar mass and isotherm intersections")
62  .declare_key("solubility", Array(Double(0.0)), Default::optional(), //("-1.0"), //
63  "Specifies solubility limits of all the sorbing species.")
64  .declare_key("table_limits", Array(Double(-1.0)), Default::optional(), //("-1.0"), //
65  "Specifies the highest aqueous concentration in the isotherm function interpolation table. "
66  "Use any negative value for an automatic choice according to current maximal concentration (default and recommended). "
67  "Use '0' to always evaluate isotherm function directly (can be very slow). "
68  "Use a positive value to set the interpolation table limit manually "
69  "(if aqueous concentration is higher, then the isotherm function is evaluated directly).")
70  .declare_key("input_fields", Array(EqFields("","").input_field_set_.make_field_descriptor_type("Sorption")), Default::obligatory(), //
71  "Containes region specific data necessary to construct isotherms.")//;
72  .declare_key("reaction_liquid", ReactionTerm::it_abstract_reaction(), Default::optional(), "Reaction model following the sorption in the liquid.")
73  .declare_key("reaction_solid", ReactionTerm::it_abstract_reaction(), Default::optional(), "Reaction model following the sorption in the solid.")
74  .close();
75 }
76 
77 
78 SorptionBase::EqFields::EqFields(const string &output_field_name, const string &output_field_desc)
80 {
81  *this += rock_density.name("rock_density")
82  .description("Rock matrix density.")
83  .input_default("0.0")
84  .units( UnitSI().kg().m(-3) );
85 
86  *this += sorption_type.name("sorption_type")
87  .description("Considered sorption is described by selected isotherm.\n"
88  "If porosity on an element is equal to 1.0 (or even higher), meaning no sorbing surface, then type 'none' will be selected automatically.")
90  .input_default("\"none\"")
92 
93  *this += distribution_coefficient.name("distribution_coefficient")
94  .description("Distribution coefficient (( $k_l, k_F, k_L $)) of linear, Freundlich or Langmuir isotherm respectively.")
95  .input_default("1.0")
96  .units( UnitSI().m(3).kg(-1) );
97 
98  *this += isotherm_other.name("isotherm_other")
99  .description("Additional parameter (($ \\alpha $)) of nonlinear isotherms.")
100  .input_default("1.0")
102 
103  *this += init_conc_solid.name("init_conc_solid")
104  .description("Initial solid concentration of substances. It is a vector: one value for every substance.")
105  .input_default("0")
106  .units( UnitSI().dimensionless() );
107 
108  input_field_set_ += *this;
109 
110  // porosity field is set from governing equation (transport) later
111  // hence we do not add it to the input_field_set_
112  *this += porosity
113  .name("porosity")
116  .set_limits(0.0);
117 
118  output_fields += *this;
119  output_fields += conc_solid.name(output_field_name)
120  .description(output_field_desc)
121  .units( UnitSI().dimensionless() )
123 
124  *this += scale_aqua.name("scale_aqua")
125  .description("Scale aqua computed by model.")
126  .input_default("0.0")
128 
129  *this += scale_sorbed.name("scale_sorbed")
130  .description("Scale sorbed computed by model.")
131  .input_default("0.0")
132  .units( UnitSI().kg().m(-3) );
133 
134  *this += no_sorbing_surface_cond.name("no_sorbing_surface_cond")
135  .description("No sorbing surface condition computed by model.")
136  .input_default("0.0")
138 
139 }
140 
141 
143 : ReactionTerm::EqData() {}
144 
145 
147  : ReactionTerm(init_mesh, in_rec),
148  init_condition_assembly_(nullptr),
149  reaction_assembly_(nullptr)
150 {
151  eq_data_ = std::make_shared<EqData>();
152  this->eq_data_base_ = std::static_pointer_cast<ReactionTerm::EqData>(eq_data_);
153 
154  // creating reaction from input and setting their parameters
155  make_reactions();
156 }
157 
158 
160 {
162  if (reaction_assembly_!=nullptr) delete reaction_assembly_;
163 }
164 
166 {
168 
169  reactions_it = input_record_.find<Input::AbstractRecord>("reaction_liquid");
170  if ( reactions_it )
171  {
172  // TODO: allowed instances in this case are only
173  // FirstOrderReaction, RadioactiveDecay
174  reaction_liquid = (*reactions_it).factory< ReactionTerm, Mesh &, Input::Record >(*mesh_, *reactions_it);
175  } else
176  {
177  reaction_liquid = nullptr;
178  }
179 
180  reactions_it = input_record_.find<Input::AbstractRecord>("reaction_solid");
181  if ( reactions_it )
182  {
183  // TODO: allowed instances in this case are only
184  // FirstOrderReaction, RadioactiveDecay
185  reaction_solid = (*reactions_it).factory< ReactionTerm, Mesh &, Input::Record >(*mesh_, *reactions_it);
186  } else
187  {
188  reaction_solid = nullptr;
189  }
190 }
191 
193 {
194  ASSERT_PERMANENT_PTR(time_).error("Time governor has not been set yet.\n");
195  ASSERT_PERMANENT(output_stream_).error("Null output stream.\n");
196  ASSERT_PERMANENT_LT(0, eq_data_->substances_.size());
197 
198  initialize_substance_ids(); //computes present substances and sets indices
199  initialize_from_input(); //reads non-field data from input
200 
201  //isotherms array resized bellow
202  unsigned int nr_of_regions = mesh_->region_db().bulk_size();
203  eq_data_->isotherms.resize(nr_of_regions);
204  eq_data_->max_conc.resize(nr_of_regions);
205  for(unsigned int i_reg = 0; i_reg < nr_of_regions; i_reg++)
206  {
207  eq_data_->isotherms[i_reg].resize(eq_data_->n_substances_);
208  eq_data_->max_conc[i_reg].resize(eq_data_->n_substances_, 0.0);
209  for(unsigned int i_subst = 0; i_subst < eq_data_->n_substances_; i_subst++)
210  {
211  eq_data_->isotherms[i_reg][i_subst] = Isotherm();
212  }
213  }
214 
217 
218  if(reaction_liquid)
219  {
220  reaction_liquid->substances(eq_data_->substances_)
221  .concentration_fields(eq_fields_->conc_mobile_fe)
222  .set_time_governor(*time_);
223  reaction_liquid->initialize();
224  }
225  if(reaction_solid)
226  {
227  reaction_solid->substances(eq_data_->substances_)
228  .concentration_fields(eq_fields_->conc_solid_fe)
229  .set_time_governor(*time_);
230  reaction_solid->initialize();
231  }
232 
235 }
236 
237 
239 {
240  Input::Array substances_array = input_record_.val<Input::Array>("substances");
241  unsigned int k, global_idx, i_subst = 0;
242  bool found;
243  Input::Iterator<string> spec_iter = substances_array.begin<string>();
244 
245  for(; spec_iter != substances_array.end(); ++spec_iter, i_subst++)
246  {
247  //finding the name of a substance in the global array of names
248  found = false;
249  for(k = 0; k < eq_data_->substances_.size(); k++)
250  {
251  if (*spec_iter == eq_data_->substances_[k].name())
252  {
253  global_idx = k;
254  found = true;
255  break;
256  }
257  }
258 
259  if(!found)
260  THROW(ReactionTerm::ExcUnknownSubstance()
261  << ReactionTerm::EI_Substance(*spec_iter)
262  << substances_array.ei_address());
263 
264  //finding the global index of substance in the local array
265  found = false;
266  for(k = 0; k < eq_data_->substance_global_idx_.size(); k++)
267  {
268  if(eq_data_->substance_global_idx_[k] == global_idx)
269  {
270  found = true;
271  break;
272  }
273  }
274 
275  if(!found)
276  {
277  eq_data_->substance_global_idx_.push_back(global_idx);
278  }
279 
280  }
281  eq_data_->n_substances_ = eq_data_->substance_global_idx_.size();
282 }
283 
285 {
286  // read number of interpolation steps - value checked by the record definition
287  n_interpolation_steps_ = input_record_.val<int>("substeps");
288 
289  // read the density of solvent - value checked by the record definition
290  eq_data_->solvent_density_ = input_record_.val<double>("solvent_density");
291 
292  // read the solubility vector
293  Input::Iterator<Input::Array> solub_iter = input_record_.find<Input::Array>("solubility");
294  if( solub_iter )
295  {
296  if (solub_iter->Array::size() != eq_data_->n_substances_)
297  {
298  THROW(SorptionBase::ExcSubstanceCountMatch()
299  << SorptionBase::EI_ArrayName("solubility")
301  // there is no way to get ei_address from 'solub_iter', only as a string
302  }
303 
304  else solub_iter->copy_to(eq_data_->solubility_vec_);
305  }
306  else{
307  // fill solubility_vec_ with zeros
308  eq_data_->solubility_vec_.clear();
309  eq_data_->solubility_vec_.resize(eq_data_->n_substances_,0.0);
310  }
311 
312  // read the interpolation table limits
313  Input::Iterator<Input::Array> interp_table_limits = input_record_.find<Input::Array>("table_limits");
314  if( interp_table_limits )
315  {
316  if (interp_table_limits->Array::size() != eq_data_->n_substances_)
317  {
318  THROW(SorptionBase::ExcSubstanceCountMatch()
319  << SorptionBase::EI_ArrayName("table_limits")
321  // there is no way to get ei_address from 'interp_table_limits', only as a string
322  }
323 
324  else interp_table_limits->copy_to(eq_data_->table_limit_);
325  }
326  else{
327  // fill table_limit_ with negative values -> automatic limit
328  eq_data_->table_limit_.clear();
329  eq_data_->table_limit_.resize(eq_data_->n_substances_,-1.0);
330  }
331 }
332 
334 {
335  ASSERT_GT(eq_data_->n_substances_, 0).error("Number of substances is wrong, they might have not been set yet.\n");
336 
337  // create vector of substances that are involved in sorption
338  // and initialize eq_fields_ with their names
339  std::vector<std::string> substances_sorption;
340  for (unsigned int i : eq_data_->substance_global_idx_)
341  substances_sorption.push_back(eq_data_->substances_[i].name());
342  eq_fields_->set_components(substances_sorption);
343 
344  // read fields from input file
345  eq_fields_->input_field_set_.set_input_list(input_record_.val<Input::Array>("input_fields"), *time_);
346 
347  eq_fields_->set_mesh(*mesh_);
348 
349  //initialization of output
350  //output_array = input_record_.val<Input::Array>("output_fields");
351  eq_fields_->conc_solid.set_components(eq_data_->substances_.names());
352  eq_fields_->output_fields.set_mesh(*mesh_);
353  eq_fields_->output_fields.output_type(OutputTime::ELEM_DATA);
354  eq_fields_->conc_solid.setup_components();
355 
356  //creating field fe and output multifield for sorbed concentrations
357  eq_fields_->conc_solid_fe.resize(eq_data_->substances_.size());
358  for (unsigned int sbi = 0; sbi < eq_data_->substances_.size(); sbi++)
359  {
360  // create shared pointer to a FieldFE and push this Field to output_field on all regions
361  eq_fields_->conc_solid_fe[sbi] = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dof_handler_);
362  eq_fields_->conc_solid[sbi].set(eq_fields_->conc_solid_fe[sbi], 0);
363  }
364  //output_stream_->add_admissible_field_names(output_array);
365  eq_fields_->output_fields.initialize(output_stream_, mesh_, input_record_.val<Input::Record>("output"), time());
366 }
367 
368 
370 {
371  ASSERT_PTR(time_).error("Time governor has not been set yet.\n");
372  ASSERT(output_stream_).error("Null output stream.\n");
373  ASSERT_LT(0, eq_data_->substances_.size());
374 
375  eq_fields_->set_time(time_->step(), LimitSide::right);
376  std::stringstream ss; // print warning message with table of uninitialized fields
377  if ( FieldCommon::print_message_table(ss, "sorption") ) {
378  WarningOut() << ss.str();
379  }
381 
382  update_max_conc();
383  make_tables();
384 
385  // write initial condition
386  //eq_fields_->output_fields.set_time(time_->step(), LimitSide::right);
387  //eq_fields_->output_fields.output(output_stream_);
388 
389  if(reaction_liquid) reaction_liquid->zero_time_step();
390  if(reaction_solid) reaction_solid->zero_time_step();
391 
392  output_data();
393 }
394 
395 
397 {
398  eq_fields_->set_time(time_->step(), LimitSide::right); // set to the last computed time
399 
400  // if parameters changed during last time step, reinit isotherms and eventualy
401  // update interpolation tables in the case of constant rock matrix parameters
402  make_tables();
403  clear_max_conc();
404 
405  START_TIMER("Sorption");
406  try{
407  reaction_assembly_->assemble(eq_data_->dof_handler_);
408  }
409  catch(ExceptionBase const &e)
410  {
411  e << input_record_.ei_address();
412  throw;
413  }
414 // for ( DHCellAccessor dh_cell : eq_data_->dof_handler_->own_range() )
415 // {
416 // compute_reaction(dh_cell);
417 // }
418  END_TIMER("Sorption");
419 
420  if(reaction_liquid) reaction_liquid->update_solution();
421  if(reaction_solid) reaction_solid->update_solution();
422 }
423 
424 //void SorptionBase::isotherm_reinit(unsigned int i_subst, const ElementAccessor<3> &elem)
425 //{
426 // START_TIMER("SorptionBase::isotherm_reinit");
427 //
428 // double mult_coef = eq_fields_->distribution_coefficient[i_subst].value(elem.centre(),elem);
429 // double second_coef = eq_fields_->isotherm_other[i_subst].value(elem.centre(),elem);
430 //
431 // int reg_idx = elem.region().bulk_idx();
432 // Isotherm & isotherm = eq_data_->isotherms[reg_idx][i_subst];
433 //
434 // bool limited_solubility_on = eq_data_->solubility_vec_[i_subst] > 0.0;
435 //
436 // double scale_aqua = eq_fields_->scale_aqua.value(elem.centre(), elem);
437 // double scale_sorbed = eq_fields_->scale_sorbed.value(elem.centre(), elem);
438 // double no_sorbing_surface_cond = eq_fields_->no_sorbing_surface_cond.value(elem.centre(), elem);
439 //
440 // // in case of no sorbing surface, set isotherm type None
441 // if( no_sorbing_surface_cond <= std::numeric_limits<double>::epsilon())
442 // {
443 // isotherm.reinit(Isotherm::none, false, eq_data_->solvent_density_,
444 // scale_aqua, scale_sorbed,
445 // 0,0,0);
446 // return;
447 // }
448 //
449 // if ( scale_sorbed <= 0.0)
450 // THROW( ExcNotPositiveScaling() << EI_Subst(i_subst) );
451 //
452 // isotherm.reinit(Isotherm::SorptionType(eq_fields_->sorption_type[i_subst].value(elem.centre(),elem)),
453 // limited_solubility_on, eq_data_->solvent_density_,
454 // scale_aqua, scale_sorbed,
455 // eq_data_->solubility_vec_[i_subst], mult_coef, second_coef);
456 //}
457 //
458 //void SorptionBase::isotherm_reinit_all(const ElementAccessor<3> &elem)
459 //{
460 // for(unsigned int i_subst = 0; i_subst < eq_data_->n_substances_; i_subst++)
461 // {
462 // isotherm_reinit(i_subst, elem);
463 // }
464 //}
465 
467 {
468  unsigned int reg_idx, i_subst;
469 
470  // clear max concetrations array
471  unsigned int nr_of_regions = mesh_->region_db().bulk_size();
472  for(reg_idx = 0; reg_idx < nr_of_regions; reg_idx++)
473  for(i_subst = 0; i_subst < eq_data_->n_substances_; i_subst++)
474  eq_data_->max_conc[reg_idx][i_subst] = 0.0;
475 }
476 
478 {
479  unsigned int reg_idx, i_subst, subst_id;
480 
481  clear_max_conc();
482 
483  for ( DHCellAccessor dh_cell : eq_data_->dof_handler_->own_range() ) {
484  IntIdx dof_p0 = dh_cell.get_loc_dof_indices()[0];
485  reg_idx = dh_cell.elm().region().bulk_idx();
486  for(i_subst = 0; i_subst < eq_data_->n_substances_; i_subst++){
487  subst_id = eq_data_->substance_global_idx_[i_subst];
488  eq_data_->max_conc[reg_idx][i_subst] = std::max(eq_data_->max_conc[reg_idx][i_subst], eq_fields_->conc_mobile_fe[subst_id]->vec().get(dof_p0));
489  }
490  }
491 }
492 
494 {
495  START_TIMER("SorptionBase::make_tables");
496  try
497  {
498  ElementAccessor<3> elm;
499  for(const Region &reg_iter: this->mesh_->region_db().get_region_set("BULK"))
500  {
501  int reg_idx = reg_iter.bulk_idx();
502  // true if data has been changed and are constant on the region
503  bool call_reinit = eq_fields_->changed() && eq_fields_->is_constant(reg_iter);
504 
505 // if(call_reinit)
506 // {
507 // ElementAccessor<3> elm(this->mesh_, reg_iter);
508 //// DebugOut().fmt("isotherm reinit\n");
509 // isotherm_reinit_all(elm);
510 // }
511 
512  // find table limit and create interpolation table for every substance
513  for(unsigned int i_subst = 0; i_subst < eq_data_->n_substances_; i_subst++){
514 
515  // clear interpolation tables, if not spacially constant OR switched off
516  if(! eq_fields_->is_constant(reg_iter) || eq_data_->table_limit_[i_subst] == 0.0){
517  eq_data_->isotherms[reg_idx][i_subst].clear_table();
518 // DebugOut().fmt("limit: 0.0 -> clear table\n");
519  continue;
520  }
521 
522  // if true then make_table will be called at the end
523  bool call_make_table = call_reinit;
524  // initialy try to keep the current table limit (it is zero at zero time step)
525  double subst_table_limit = eq_data_->isotherms[reg_idx][i_subst].table_limit();
526 
527  // if automatic, possibly remake tables with doubled range when table maximum was reached
528  if(eq_data_->table_limit_[i_subst] < 0.0)
529  {
530  if(subst_table_limit < eq_data_->max_conc[reg_idx][i_subst])
531  {
532  call_make_table = true;
533  subst_table_limit = 2*eq_data_->max_conc[reg_idx][i_subst];
534 // DebugOut().fmt("limit: max conc\n");
535  }
536  }
537  // if not automatic, set given table limit
538  else
539  {
540  subst_table_limit = eq_data_->table_limit_[i_subst];
541  }
542 
543  if(call_make_table){
544  eq_data_->isotherms[reg_idx][i_subst].make_table(n_interpolation_steps_, subst_table_limit);
545 // DebugOut().fmt("reg: {} i_subst {}: table_limit = {}\n", reg_idx, i_subst, eq_data_->isotherms[reg_idx][i_subst].table_limit());
546  }
547  }
548  }
549  }
550  catch(ExceptionBase const &e)
551  {
552  e << input_record_.ei_address();
553  throw;
554  }
555 }
556 
558 {
559 // const ElementAccessor<3> ele = dh_cell.elm();
560 // int reg_idx = ele.region().bulk_idx();
561 // IntIdx dof_p0 = dh_cell.get_loc_dof_indices()[0];
562 // unsigned int i_subst, subst_id;
563 // // for checking, that the common element data are computed once at maximum
564 //
565 // try{
566 // for(i_subst = 0; i_subst < eq_data_->n_substances_; i_subst++)
567 // {
568 // subst_id = eq_data_->substance_global_idx_[i_subst];
569 // Isotherm & isotherm = eq_data_->isotherms[reg_idx][i_subst];
570 // if (isotherm.is_precomputed()){
571 //// DebugOut().fmt("isotherms precomputed - interpolate, subst[{}]\n", i_subst);
572 // double c_aqua = eq_fields_->conc_mobile_fe[subst_id]->vec().get(dof_p0);
573 // double c_sorbed = eq_fields_->conc_solid_fe[subst_id]->vec().get(dof_p0);
574 // isotherm.interpolate(c_aqua, c_sorbed);
575 // eq_fields_->conc_mobile_fe[subst_id]->vec().set(dof_p0, c_aqua);
576 // eq_fields_->conc_solid_fe[subst_id]->vec().set(dof_p0, c_sorbed);
577 // }
578 // else{
579 // DebugOut().fmt("isotherms reinit - compute , subst[{}]\n", i_subst);
580 //
581 // isotherm_reinit(i_subst, ele);
582 // double c_aqua = eq_fields_->conc_mobile_fe[subst_id]->vec().get(dof_p0);
583 // double c_sorbed = eq_fields_->conc_solid_fe[subst_id]->vec().get(dof_p0);
584 // isotherm.compute(c_aqua, c_sorbed);
585 // eq_fields_->conc_mobile_fe[subst_id]->vec().set(dof_p0, c_aqua);
586 // eq_fields_->conc_solid_fe[subst_id]->vec().set(dof_p0, c_sorbed);
587 // }
588 //
589 // // update maximal concentration per region (optimization for interpolation)
590 // if(eq_data_->table_limit_[i_subst] < 0)
591 // eq_data_->max_conc[reg_idx][i_subst] = std::max(eq_data_->max_conc[reg_idx][i_subst],
592 // eq_fields_->conc_mobile_fe[subst_id]->vec().get(dof_p0));
593 // }
594 // }
595 // catch(ExceptionBase const &e)
596 // {
597 // e << input_record_.ei_address();
598 // throw;
599 // }
600 }
601 
602 
603 /**************************************** OUTPUT ***************************************************/
604 
606 {
607  eq_fields_->output_fields.set_time(time().step(), LimitSide::right);
608  // Register fresh output data
609  eq_fields_->output_fields.output(time().step());
610 }
#define ASSERT_PERMANENT_LT(a, b)
Definition of comparative assert macro (Less Than)
Definition: asserts.hh:297
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR)
Definition: asserts.hh:337
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
Cell accessor allow iterate over DOF handler cells.
Input::Record input_record_
Definition: equation.hh:242
TimeGovernor * time_
Definition: equation.hh:241
Mesh * mesh_
Definition: equation.hh:240
TimeGovernor & time()
Definition: equation.hh:151
Base of exceptions used in Flow123d.
Definition: exceptions.hh:76
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:95
FieldCommon & input_selection(Input::Type::Selection element_selection)
FieldCommon & description(const string &description)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
FieldCommon & set_limits(double min, double max=std::numeric_limits< double >::max())
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
FieldCommon & input_default(const string &input_default)
static constexpr Mask input_copy
Definition: field_flag.hh:44
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
Input::Type::Record make_field_descriptor_type(const std::string &equation_name) const
Definition: field_set.cc:95
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
EI_Address ei_address() const
Definition: accessors.cc:314
Iterator< ValueType > begin() const
IteratorBase end() const
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
EI_Address ei_address() const
Definition: accessors.cc:178
const Ret val(const string &key) const
Iterator< Ret > find(const string &key) const
Class for declaration of inputs sequences.
Definition: type_base.hh:339
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
static Default optional()
The factory function to make an empty default value which is optional.
Definition: type_record.hh:124
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Class for declaration of the integral input data.
Definition: type_base.hh:483
Record type proxy class.
Definition: type_record.hh:182
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Template for classes storing finite set of named values.
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
const Selection & close() const
Close the Selection, no more values can be added.
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
@ freundlich
Definition: isotherm.hh:177
@ langmuir
Definition: isotherm.hh:178
const RegionDB & region_db() const
Definition: mesh.h:175
Definition: mesh.h:362
EqFields()
Constructor.
std::shared_ptr< OutputTime > output_stream_
Pointer to a transport output stream.
std::shared_ptr< EqData > eq_data_base_
Equation data - all data needs in assembly class.
static Input::Type::Abstract & it_abstract_reaction()
ReactionTerm(Mesh &init_mesh, Input::Record in_rec)
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
unsigned int bulk_size() const
Definition: region.cc:268
DualPorosity data.
EqData()
Collect all fields.
MultiField< 3, FieldValue< 3 >::Scalar > init_conc_solid
Initial sorbed concentrations.
MultiField< 3, FieldValue< 3 >::Scalar > isotherm_other
Langmuir sorption coeficients alpha (in fraction c_s = omega * (alpha*c_a)/(1- alpha*c_a)).
FieldSet input_field_set_
Input data set - fields in this set are read from the input file.
Field< 3, FieldValue< 3 >::Scalar > scale_sorbed
Field< 3, FieldValue< 3 >::Scalar > rock_density
Rock matrix density.
MultiField< 3, FieldValue< 3 >::Scalar > conc_solid
Calculated sorbed concentrations, for output only.
MultiField< 3, FieldValue< 3 >::Enum > sorption_type
Discrete need Selection for initialization.
Field< 3, FieldValue< 3 >::Scalar > scale_aqua
Instances of FieldModel used in assembly methods.
static const Input::Type::Selection & get_sorption_type_selection()
Field< 3, FieldValue< 3 >::Scalar > porosity
Porosity field copied from transport.
Field< 3, FieldValue< 3 >::Scalar > no_sorbing_surface_cond
MultiField< 3, FieldValue< 3 >::Scalar > distribution_coefficient
Multiplication coefficients (k, omega) for all types of isotherms.
EquationOutput output_fields
Fields indended for output, i.e. all input fields plus those representing solution.
void compute_reaction(const DHCellAccessor &dh_cell) override
Compute reaction on a single element.
void zero_time_step() override
void output_data(void) override
Output method.
static const Input::Type::Record & get_input_type()
void make_reactions()
std::shared_ptr< ReactionTerm > reaction_liquid
GenericAssembly< ReactionAssemblySorp > * reaction_assembly_
void initialize() override
Prepares the object to usage.
virtual ~SorptionBase(void)
std::shared_ptr< ReactionTerm > reaction_solid
std::shared_ptr< EqData > eq_data_
Equation data.
void initialize_from_input()
Initializes private members of sorption from the input record.
unsigned int n_interpolation_steps_
void clear_max_conc()
Sets max conc to zeros on all regins.
void initialize_substance_ids()
Reads names of substances from input and creates indexing to global vector of substance.
GenericAssembly< InitConditionAssemblySorp > * init_condition_assembly_
general assembly objects, hold assembly objects of appropriate dimension
void update_max_conc()
Computes maximal aqueous concentration at the current step.
void initialize_fields()
Initializes field sets.
std::shared_ptr< EqFields > eq_fields_
Pointer to equation fields. The object is constructed in descendants.
virtual void init_field_models()=0
Initialize FieldModels, method is implemented in descendants.
void make_tables(void)
void update_solution(void) override
Updates the solution.
const TimeStep & step(int index=-1) const
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
Support classes for parallel programing.
Class Dual_por_exchange implements the model of dual porosity.
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
int IntIdx
Definition: index_types.hh:25
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define FMT_UNUSED
Definition: posix.h:75
Class ReactionTerm is an abstract class representing reaction term in transport.
Class SorptionBase is abstract class representing model of sorption in transport.
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.