Flow123d  release_2.2.0-914-gf1a3a4f
dual_porosity.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 dual_porosity.cc
15  * @brief
16  */
17 
18 #include <iostream>
19 #include <stdlib.h>
20 #include <math.h>
21 
24 #include "system/system.hh"
25 #include "system/sys_profiler.hh"
26 
27 #include "la/distribution.hh"
28 #include "mesh/mesh.h"
29 #include "mesh/elements.h"
30 #include "mesh/region.hh"
32 
33 #include "reaction/sorption.hh"
37 #include "input/factory.hh"
38 
39 FLOW123D_FORCE_LINK_IN_CHILD(dualPorosity)
40 
41 
42 using namespace Input::Type;
43 
44 
45 
46 const Record & DualPorosity::get_input_type() {
47  return Record("DualPorosity",
48  "Dual porosity model in transport problems.\n"
49  "Provides computing the concentration of substances in mobile and immobile zone.\n"
50  )
52  .declare_key("input_fields", Array(DualPorosity::EqData().make_field_descriptor_type("DualPorosity")), Default::obligatory(),
53  "Containes region specific data necessary to construct dual porosity model.")
54  .declare_key("scheme_tolerance", Double(0.0), Default("1e-3"),
55  "Tolerance according to which the explicit Euler scheme is used or not."
56  "Set 0.0 to use analytic formula only (can be slower).")
57 
58  .declare_key("reaction_mobile", ReactionTerm::it_abstract_mobile_term(), Default::optional(), "Reaction model in mobile zone.")
59  .declare_key("reaction_immobile", ReactionTerm::it_abstract_immobile_term(), Default::optional(), "Reaction model in immobile zone.")
60  .declare_key("output",
61  EqData().output_fields.make_output_type("DualPorosity", ""),
62  IT::Default("{ \"fields\": [ \"conc_immobile\" ] }"),
63  "Setting of the fields output.")
64  .close();
65 }
66 
67 const int DualPorosity::registrar =
68  Input::register_class< DualPorosity, Mesh &, Input::Record >("DualPorosity") +
70 
72 {
74  .name("diffusion_rate_immobile")
75  .description("Diffusion coefficient of non-equilibrium linear exchange between mobile and immobile zone.")
76  .input_default("0")
77  .units( UnitSI().s(-1) );
78 
79  *this += porosity_immobile
80  .name("porosity_immobile")
81  .description("Porosity of the immobile zone.")
82  .input_default("0")
84  .set_limits(0.0);
85 
86  *this += init_conc_immobile
87  .name("init_conc_immobile")
88  .description("Initial concentration of substances in the immobile zone.")
89  .units( UnitSI().kg().m(-3) );
90 
91  //creating field for porosity that is set later from the governing equation (transport)
92  *this +=porosity
93  .name("porosity")
96  .set_limits(0.0);
97 
98  *this += conc_immobile
99  .name("conc_immobile")
100  .units( UnitSI().kg().m(-3) )
102 
103  output_fields += *this;
104 
105 }
106 
108  : ReactionTerm(init_mesh, in_rec)
109 {
110  //set pointer to equation data fieldset
111  this->eq_data_ = &data_;
112 
113  //reads input and creates possibly other reaction terms
114  make_reactions();
115  //read value from input
116  scheme_tolerance_ = input_record_.val<double>("scheme_tolerance");
117 }
118 
120 {
121  VecScatterDestroy(&(vconc_out_scatter));
122 
123 
124  for (unsigned int sbi = 0; sbi < substances_.size(); sbi++)
125  {
126 
127  //no mpi vectors
128  VecDestroy(&(vconc_immobile[sbi]));
129  delete [] conc_immobile[sbi];
130  }
131 
132  delete [] vconc_immobile;
133  delete [] conc_immobile;
134 }
135 
136 
139  if ( reactions_it )
140  {
141  // TODO: allowed instances in this case are only
142  // FirstOrderReaction, RadioactiveDecay and SorptionMob
143  reaction_mobile = (*reactions_it).factory< ReactionTerm, Mesh &, Input::Record >(*mesh_, *reactions_it);
144  } else
145  {
146  reaction_mobile = nullptr;
147  }
148 
149  reactions_it = input_record_.find<Input::AbstractRecord>("reaction_immobile");
150  if ( reactions_it )
151  {
152  // TODO: allowed instances in this case are only
153  // FirstOrderReaction, RadioactiveDecay and SorptionImmob
154  reaction_immobile = (*reactions_it).factory< ReactionTerm, Mesh &, Input::Record >(*mesh_, *reactions_it);
155  } else
156  {
157  reaction_immobile = nullptr;
158  }
159 
160 }
161 
163 {
164  OLD_ASSERT(distribution_ != nullptr, "Distribution has not been set yet.\n");
165  OLD_ASSERT(time_ != nullptr, "Time governor has not been set yet.\n");
166  OLD_ASSERT(output_stream_,"Null output stream.");
168 
169  //allocating memory for immobile concentration matrix
170  conc_immobile = new double* [substances_.size()];
171  conc_immobile_out.clear();
172  conc_immobile_out.resize( substances_.size() );
173  for (unsigned int sbi = 0; sbi < substances_.size(); sbi++)
174  {
175  conc_immobile[sbi] = new double [distribution_->lsize()];
176  conc_immobile_out[sbi].resize( distribution_->size() );
177  }
179 
181 
182  if(reaction_mobile)
183  {
184  reaction_mobile->substances(substances_)
185  .output_stream(output_stream_)
186  .concentration_matrix(concentration_matrix_, distribution_, el_4_loc_, row_4_el_)
187  .set_time_governor(*time_);
188  reaction_mobile->initialize();
189  }
190 
192  {
193  reaction_immobile->substances(substances_)
194  .output_stream(output_stream_)
195  .concentration_matrix(conc_immobile, distribution_, el_4_loc_, row_4_el_)
196  .set_time_governor(*time_);
197  reaction_immobile->initialize();
198  }
199 
200 }
201 
203 {
205  //setting fields that are set from input file
208 
209  //setting fields in data
210  data_.set_mesh(*mesh_);
211 
212  //initialization of output
213  //output_array = input_record_.val<Input::Array>("output_fields");
218  for (unsigned int sbi=0; sbi<substances_.size(); sbi++)
219  {
220  // create shared pointer to a FieldElementwise and push this Field to output_field on all regions
221  auto output_field_ptr = conc_immobile_out[sbi].create_field<3, FieldValue<3>::Scalar>(substances_.size());
222  data_.conc_immobile[sbi].set_field(mesh_->region_db().get_region_set("ALL"), output_field_ptr, 0);
223  }
224  //output_stream_->add_admissible_field_names(output_array);
226 }
227 
228 
230 {
231  OLD_ASSERT(distribution_ != nullptr, "Distribution has not been set yet.\n");
232  OLD_ASSERT(time_ != nullptr, "Time governor has not been set yet.\n");
233  OLD_ASSERT(output_stream_,"Null output stream.");
235 
236  //coupling - passing fields
237  if(reaction_mobile)
238  if (typeid(*reaction_mobile) == typeid(SorptionMob))
239  {
240  reaction_mobile->data().set_field("porosity", data_["porosity"]);
241  reaction_mobile->data().set_field("porosity_immobile", data_["porosity_immobile"]);
242  }
244  if (typeid(*reaction_immobile) == typeid(SorptionImmob))
245  {
246  reaction_immobile->data().set_field("porosity", data_["porosity"]);
247  reaction_immobile->data().set_field("porosity_immobile", data_["porosity_immobile"]);
248  }
249 
251  std::stringstream ss; // print warning message with table of uninitialized fields
252  if ( FieldCommon::print_message_table(ss, "dual porosity") ) {
253  WarningOut() << ss.str();
254  }
256 
257  // write initial condition
258  //output_vector_gather();
259  //data_.output_fields.set_time(time_->step(0), LimitSide::right);
260  //data_.output_fields.output(output_stream_);
261 
262  output_data();
263 
264  if(reaction_mobile)
265  reaction_mobile->zero_time_step();
266 
268  reaction_immobile->zero_time_step();
269 
270 }
271 
273 {
274  //setting initial condition for immobile concentration matrix
275  for (unsigned int loc_el = 0; loc_el < distribution_->lsize(); loc_el++)
276  { // Optimize: SWAP LOOPS
277  unsigned int index = el_4_loc_[loc_el];
278  ElementAccessor<3> ele_acc = mesh_->element_accessor(index);
279 
280  for (unsigned int sbi=0; sbi < substances_.size(); sbi++)
281  {
282  conc_immobile[sbi][loc_el] = data_.init_conc_immobile[sbi].value(ele_acc.centre(), ele_acc);
283  }
284  }
285 }
286 
288 {
290 
291  START_TIMER("dual_por_exchange_step");
292  for (unsigned int loc_el = 0; loc_el < distribution_->lsize(); loc_el++)
293  {
295  }
296  END_TIMER("dual_por_exchange_step");
297 
298  if(reaction_mobile) reaction_mobile->update_solution();
299  if(reaction_immobile) reaction_immobile->update_solution();
300 }
301 
302 
303 double **DualPorosity::compute_reaction(double **concentrations, int loc_el)
304 {
305  unsigned int sbi;
306  double conc_average, // weighted (by porosity) average of concentration
307  conc_mob, conc_immob, // new mobile and immobile concentration
308  previous_conc_mob, previous_conc_immob, // mobile and immobile concentration in previous time step
309  conc_max, //difference between concentration and average concentration
310  por_mob, por_immob; // mobile and immobile porosity
311 
312  // get data from fields
313  ElementFullIter ele = mesh_->element(el_4_loc_[loc_el]);
314  por_mob = data_.porosity.value(ele->centre(),ele->element_accessor());
315  por_immob = data_.porosity_immobile.value(ele->centre(),ele->element_accessor());
316  arma::Col<double> diff_vec(substances_.size());
317  for (sbi=0; sbi<substances_.size(); sbi++) // Optimize: SWAP LOOPS
318  diff_vec[sbi] = data_.diffusion_rate_immobile[sbi].value(ele->centre(), ele->element_accessor());
319 
320  // if porosity_immobile == 0 then mobile concentration stays the same
321  // and immobile concentration cannot change
322  if (por_immob == 0.0) return conc_immobile;
323 
324  double exponent,
325  temp_exponent = (por_mob + por_immob) / (por_mob * por_immob) * time_->dt();
326 
327  for (sbi = 0; sbi < substances_.size(); sbi++) //over all substances
328  {
329  exponent = diff_vec[sbi] * temp_exponent;
330  //previous values
331  previous_conc_mob = concentration_matrix_[sbi][loc_el];
332  previous_conc_immob = conc_immobile[sbi][loc_el];
333 
334  // ---compute average concentration------------------------------------------
335  conc_average = ((por_mob * previous_conc_mob) + (por_immob * previous_conc_immob))
336  / (por_mob + por_immob);
337 
338  conc_max = std::max(previous_conc_mob-conc_average, previous_conc_immob-conc_average);
339 
340  // the following 2 conditions guarantee:
341  // 1) stability of forward Euler's method
342  // 2) the error of forward Euler's method will not be large
343  if(time_->dt() <= por_mob*por_immob/(max(diff_vec)*(por_mob+por_immob)) &&
344  conc_max <= (2*scheme_tolerance_/(exponent*exponent)*conc_average)) // forward euler
345  {
346  double temp = diff_vec[sbi]*(previous_conc_immob - previous_conc_mob) * time_->dt();
347  // ---compute concentration in mobile area
348  conc_mob = temp / por_mob + previous_conc_mob;
349 
350  // ---compute concentration in immobile area
351  conc_immob = -temp / por_immob + previous_conc_immob;
352  }
353  else //analytic solution
354  {
355  double temp = exp(-exponent);
356  // ---compute concentration in mobile area
357  conc_mob = (previous_conc_mob - conc_average) * temp + conc_average;
358 
359  // ---compute concentration in immobile area
360  conc_immob = (previous_conc_immob - conc_average) * temp + conc_average;
361  }
362 
363  concentration_matrix_[sbi][loc_el] = conc_mob;
364  conc_immobile[sbi][loc_el] = conc_immob;
365  }
366 
367  return conc_immobile;
368 }
369 
370 
372 {
373  int sbi, n_subst;
374  n_subst = substances_.size();
375 
376  vconc_immobile = new Vec [n_subst];
377 
378 
379  for (sbi = 0; sbi < n_subst; sbi++) {
380  VecCreateMPIWithArray(PETSC_COMM_WORLD,1, distribution_->lsize(), mesh_->n_elements(), conc_immobile[sbi],
381  &vconc_immobile[sbi]);
382  VecZeroEntries(vconc_immobile[sbi]);
383 
384  // if(rank == 0)
385  VecZeroEntries(conc_immobile_out[sbi].get_data_petsc());
386  }
387 
388  // create output vector scatter
389  IS is;
390  ISCreateGeneral(PETSC_COMM_SELF, mesh_->n_elements(), row_4_el_, PETSC_COPY_VALUES, &is); //WithArray
391  VecScatterCreate(vconc_immobile[0], is, conc_immobile_out[0].get_data_petsc(), PETSC_NULL, &vconc_out_scatter);
392  ISDestroy(&(is));
393 }
394 
395 
397 {
398  unsigned int sbi;
399 
400  for (sbi = 0; sbi < substances_.size(); sbi++) {
401  VecScatterBegin(vconc_out_scatter, vconc_immobile[sbi], conc_immobile_out[sbi].get_data_petsc(), INSERT_VALUES, SCATTER_FORWARD);
402  VecScatterEnd(vconc_out_scatter, vconc_immobile[sbi], conc_immobile_out[sbi].get_data_petsc(), INSERT_VALUES, SCATTER_FORWARD);
403  }
404 }
405 
406 
408 {
412  }
413 
414  // Register fresh output data
416 
417  if (time_->tlevel() !=0) {
418  // zero_time_step call zero_time_Step of subreactions which performs its own output
419  if (reaction_mobile) reaction_mobile->output_data();
420  if (reaction_immobile) reaction_immobile->output_data();
421  }
422 }
423 
424 
425 bool DualPorosity::evaluate_time_constraint(double &time_constraint)
426 {
427  bool cfl_changed = false;
428  if (reaction_mobile)
429  {
430  if (reaction_mobile->evaluate_time_constraint(time_constraint))
431  cfl_changed = true;
432  }
433  if (reaction_immobile)
434  {
435  double cfl_immobile;
436  if (reaction_immobile->evaluate_time_constraint(cfl_immobile))
437  {
438  time_constraint = std::min(time_constraint, cfl_immobile);
439  cfl_changed = true;
440  }
441  }
442 
443  return cfl_changed;
444 }
TimeGovernor & time()
Definition: equation.hh:148
void output_type(OutputTime::DiscreteSpace rt)
Definition: field_set.hh:201
unsigned int size() const
get global size
void set_input_list(Input::Array input_list)
Definition: field_set.hh:180
Sorption model in immobile zone in case dual porosity is considered.
Definition: sorption.hh:123
FieldSet * eq_data_
Definition: equation.hh:232
void set_initial_condition()
Sets initial condition from input.
VecScatter vconc_out_scatter
Output vector scatter.
std::shared_ptr< ReactionTerm > reaction_mobile
Reaction running in mobile zone.
Accessor to input data conforming to declared Array.
Definition: accessors.hh:567
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:598
int tlevel() const
void output_data(void) override
Main output routine.
double ** concentration_matrix_
const std::vector< std::string > & names()
Definition: substance.hh:85
MultiField< 3, FieldValue< 3 >::Scalar > conc_immobile
Calculated concentrations in the immobile zone.
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
void output(TimeStep step)
Abstract linear system class.
Definition: equation.hh:37
SubstanceList substances_
Names belonging to substances.
EquationOutput output_fields
Fields indended for output, i.e. all input fields plus those representing solution.
RegionSet get_region_set(const string &set_name) const
Definition: region.cc:328
void zero_time_step() override
Definition: mesh.h:99
Field< 3, FieldValue< 3 >::Scalar > porosity_immobile
Immobile porosity field.
Iterator< Ret > find(const string &key) const
~DualPorosity(void)
Destructor.
std::vector< VectorSeqDouble > conc_immobile_out
concentration array output for immobile phase (gathered - sequential)
const RegionDB & region_db() const
Definition: mesh.h:170
IdxInt * el_4_loc_
Indices of elements belonging to local dofs.
ElementAccessor< 3 > element_accessor(unsigned int idx, bool boundary=false)
Definition: mesh.cc:668
void allocate_output_mpi(void)
Allocates petsc vectors, prepares them for output and creates output vector scatter.
const TimeStep & step(int index=-1) const
static Input::Type::Abstract & it_abstract_mobile_term()
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
unsigned int size() const
Definition: substance.hh:87
Class for declaration of inputs sequences.
Definition: type_base.hh:345
std::shared_ptr< OutputTime > output_stream_
Pointer to a transport output stream.
void update_solution(void) override
ReactionTerm(Mesh &init_mesh, Input::Record in_rec)
Constructor.
Class ReactionTerm is an abstract class representing reaction term in transport.
Class Dual_por_exchange implements the model of dual porosity.
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:195
double scheme_tolerance_
Dual porosity computational scheme tolerance.
#define OLD_ASSERT(...)
Definition: global_defs.h:131
void setup_components()
std::shared_ptr< ReactionTerm > reaction_immobile
Reaction running in immobile zone.
void make_reactions()
Resolves construction of following reactions.
unsigned int n_elements() const
Definition: mesh.h:156
void initialize() override
Prepares the object to usage.
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:540
Distribution * distribution_
Pointer to reference to distribution of elements between processors.
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
static constexpr Mask input_copy
Definition: field_flag.hh:44
FieldCommon & input_default(const string &input_default)
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
const Ret val(const string &key) const
#define START_TIMER(tag)
Starts a timer with specified tag.
#define OLD_ASSERT_LESS(a, b)
Definition: global_defs.h:134
bool evaluate_time_constraint(double &time_constraint) override
Computes a constraint for time step.
IdxInt * row_4_el_
Indices of rows belonging to elements.
void initialize_fields()
Initializes field sets.
static Input::Type::Abstract & it_abstract_immobile_term()
Mesh * mesh_
Definition: equation.hh:223
virtual Value::return_type const & value(const Point &p, const ElementAccessor< spacedim > &elm) const
Definition: field.hh:368
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:490
EqData()
Collect all fields.
Class representing dual porosity model in transport.
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:459
double ** conc_immobile
static Input::Type::Abstract & it_abstract_term()
Support classes for parallel programing.
Sorption model in mobile zone in case dual porosity is considered.
Definition: sorption.hh:95
void output_vector_gather(void) override
Gathers all the parallel vectors to enable them to be output.
FieldCommon & description(const string &description)
void initialize(std::shared_ptr< OutputTime > stream, Mesh *mesh, Input::Record in_rec, const TimeGovernor &tg)
Input::Record input_record_
Definition: equation.hh:225
void set_components(const std::vector< string > &names)
Definition: field_set.hh:167
double dt() const
double ** compute_reaction(double **concentrations, int loc_el) override
bool set_time(const TimeStep &time, LimitSide limit_side)
Definition: field_set.cc:157
#define WarningOut()
Macro defining &#39;warning&#39; record of log.
Definition: logger.hh:246
FieldCommon & name(const string &name)
Definition: field_common.hh:97
#define END_TIMER(tag)
Ends a timer with specified tag.
MultiField< 3, FieldValue< 3 >::Scalar > diffusion_rate_immobile
Mass transfer coefficients between mobile and immobile pores.
arma::vec::fixed< spacedim > centre() const
Definition: accessors.hh:91
static const int registrar
Registrar of class to factory.
void set_mesh(const Mesh &mesh)
Definition: field_set.hh:173
MultiField< 3, FieldValue< 3 >::Scalar > init_conc_immobile
Initial concentrations in the immobile zone.
Record type proxy class.
Definition: type_record.hh:182
bool is_field_output_time(const FieldCommon &field, TimeStep step) const
FieldCommon & set_limits(double min, double max=std::numeric_limits< double >::max())
Field< 3, FieldValue< 3 >::Scalar > porosity
Porosity field.
FieldCommon & flags(FieldFlag::Flags::Mask mask)
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static const Input::Type::Record & get_input_type()
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:91
DualPorosity data.
Vec * vconc_immobile
PETSC concentration vector for immobile phase (parallel).
This file contains classes representing sorption model. Sorption model can be computed both in case t...
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:180
TimeGovernor * time_
Definition: equation.hh:224
FieldSet input_data_set_
ElementVector element
Vector of elements of the mesh.
Definition: mesh.h:260
unsigned int lsize(int proc) const
get local size