Flow123d  master-f44eb46
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"
31 #include "mesh/accessors.hh"
32 #include "fields/field_fe.hh"
33 
34 #include "reaction/sorption.hh"
38 #include "input/factory.hh"
39 
40 FLOW123D_FORCE_LINK_IN_CHILD(dualPorosity)
41 
42 
43 using namespace Input::Type;
44 
45 
46 
48  return Record("DualPorosity",
49  "Dual porosity model in transport problems.\n"
50  "Provides computing the concentration of substances in mobile and immobile zone.\n"
51  )
53  .declare_key("input_fields", Array(DualPorosity::EqFields().make_field_descriptor_type("DualPorosity")), Default::obligatory(),
54  "Containes region specific data necessary to construct dual porosity model.")
55  .declare_key("scheme_tolerance", Double(0.0), Default("1e-3"),
56  "Tolerance according to which the explicit Euler scheme is used or not."
57  "Set 0.0 to use analytic formula only (can be slower).")
58 
59  .declare_key("reaction_mobile", ReactionTerm::it_abstract_mobile_term(), Default::optional(), "Reaction model in mobile zone.")
60  .declare_key("reaction_immobile", ReactionTerm::it_abstract_immobile_term(), Default::optional(), "Reaction model in immobile zone.")
61  .declare_key("output",
62  EqFields().output_fields.make_output_type("DualPorosity", ""),
63  IT::Default("{ \"fields\": [ \"conc_immobile\" ] }"),
64  "Setting of the fields output.")
65  .close();
66 }
67 
68 const int DualPorosity::registrar =
69  Input::register_class< DualPorosity, Mesh &, Input::Record >("DualPorosity") +
71 
74 {
76  .name("diffusion_rate_immobile")
77  .description("Diffusion coefficient of non-equilibrium linear exchange between mobile and immobile zone.")
78  .input_default("0")
79  .units( UnitSI().s(-1) );
80 
81  *this += porosity_immobile
82  .name("porosity_immobile")
83  .description("Porosity of the immobile zone.")
84  .input_default("0")
86  .set_limits(0.0);
87 
88  *this += init_conc_immobile
89  .name("init_conc_immobile")
90  .description("Initial concentration of substances in the immobile zone.")
91  .input_default("0")
92  .units( UnitSI().kg().m(-3) );
93 
94  //creating field for porosity that is set later from the governing equation (transport)
95  *this +=porosity
96  .name("porosity")
97  .description("Concentration solution in the mobile phase.")
100  .set_limits(0.0);
101 
102  *this += conc_immobile
103  .name("conc_immobile")
104  .units( UnitSI().kg().m(-3) )
106 
107  output_fields += *this;
108 
109 }
110 
112 : ReactionTerm::EqData() {}
113 
115  : ReactionTerm(init_mesh, in_rec),
116  init_condition_assembly_(nullptr),
117  reaction_assembly_(nullptr)
118 {
119  eq_fields_ = std::make_shared<EqFields>();
120  eq_fields_->add_coords_field();
121  //set pointer to equation data fieldset
122  this->eq_fieldset_ = eq_fields_;
123  this->eq_fields_base_ = std::static_pointer_cast<ReactionTerm::EqFields>(eq_fields_);
124 
125  eq_data_ = std::make_shared<EqData>();
126  this->eq_data_base_ = std::static_pointer_cast<ReactionTerm::EqData>(eq_data_);
127 
128  //reads input and creates possibly other reaction terms
129  make_reactions();
130  //read value from input
131  eq_data_->scheme_tolerance_ = input_record_.val<double>("scheme_tolerance");
132 }
133 
135 {
137  if (reaction_assembly_!=nullptr) delete reaction_assembly_;
138 }
139 
140 
143  if ( reactions_it )
144  {
145  // TODO: allowed instances in this case are only
146  // FirstOrderReaction, RadioactiveDecay and SorptionMob
147  reaction_mobile = (*reactions_it).factory< ReactionTerm, Mesh &, Input::Record >(*mesh_, *reactions_it);
148  } else
149  {
150  reaction_mobile = nullptr;
151  }
152 
153  reactions_it = input_record_.find<Input::AbstractRecord>("reaction_immobile");
154  if ( reactions_it )
155  {
156  // TODO: allowed instances in this case are only
157  // FirstOrderReaction, RadioactiveDecay and SorptionImmob
158  reaction_immobile = (*reactions_it).factory< ReactionTerm, Mesh &, Input::Record >(*mesh_, *reactions_it);
159  } else
160  {
161  reaction_immobile = nullptr;
162  }
163 
164 }
165 
167 {
168  ASSERT_PERMANENT(time_ != nullptr).error("Time governor has not been set yet.\n");
169  ASSERT_PERMANENT_LT(0, eq_data_->substances_.size()).error("No substances for rection term.\n");
170  ASSERT_PERMANENT(output_stream_ != nullptr).error("Null output stream.\n");
171 
172  eq_data_->time_ = this->time_;
173 
175 
176  if(reaction_mobile)
177  {
178  reaction_mobile->substances(eq_data_->substances_)
179  .output_stream(output_stream_)
180  .concentration_fields(eq_fields_->conc_mobile_fe)
181  .set_time_governor(*time_);
182  reaction_mobile->initialize();
183  }
184 
186  {
187  reaction_immobile->substances(eq_data_->substances_)
188  .output_stream(output_stream_)
189  .concentration_fields(eq_fields_->conc_immobile_fe)
190  .set_time_governor(*time_);
191  reaction_immobile->initialize();
192  }
193 
196 
197 }
198 
200 {
201  eq_fields_->set_components(eq_data_->substances_.names());
202  //setting fields that are set from input file
205 
206  //setting fields in data
207  eq_fields_->set_mesh(*mesh_);
208 
209  //initialization of output
210  eq_fields_->output_fields.set_components(eq_data_->substances_.names());
211  eq_fields_->output_fields.set_mesh(*mesh_);
212  eq_fields_->output_fields.output_type(OutputTime::ELEM_DATA);
213  eq_fields_->conc_immobile.setup_components();
214 
215 
216  //creating field fe and output multifield for sorbed concentrations
217  eq_fields_->conc_immobile_fe.resize(eq_data_->substances_.size());
218  for (unsigned int sbi = 0; sbi < eq_data_->substances_.size(); sbi++)
219  {
220  // create shared pointer to a FieldFE and push this Field to output_field on all regions
221  eq_fields_->conc_immobile_fe[sbi] = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dof_handler_);
222  eq_fields_->conc_immobile[sbi].set(eq_fields_->conc_immobile_fe[sbi], 0);
223  }
224 
225  eq_fields_->output_fields.initialize(output_stream_, mesh_, input_record_.val<Input::Record>("output"),time());
226 }
227 
228 
230 {
231  ASSERT(time_ != nullptr).error("Time governor has not been set yet.\n");
232  ASSERT_LT(0, eq_data_->substances_.size()).error("No substances for rection term.\n");
233  ASSERT(output_stream_ != nullptr).error("Null output stream.\n");
234 
235  //coupling - passing fields
236  if(reaction_mobile)
237  if (typeid(*reaction_mobile) == typeid(SorptionMob))
238  {
239  reaction_mobile->eq_fieldset().set_field("porosity", (*eq_fields_)["porosity"]);
240  reaction_mobile->eq_fieldset().set_field("porosity_immobile", (*eq_fields_)["porosity_immobile"]);
241  }
243  if (typeid(*reaction_immobile) == typeid(SorptionImmob))
244  {
245  reaction_immobile->eq_fieldset().set_field("porosity", (*eq_fields_)["porosity"]);
246  reaction_immobile->eq_fieldset().set_field("porosity_immobile", (*eq_fields_)["porosity_immobile"]);
247  }
248 
249  eq_fields_->set_time(time_->step(0), LimitSide::right);
250  std::stringstream ss; // print warning message with table of uninitialized fields
251  if ( FieldCommon::print_message_table(ss, "dual porosity") ) {
252  WarningOut() << ss.str();
253  }
255 
256  output_data();
257 
258  if(reaction_mobile)
259  reaction_mobile->zero_time_step();
260 
262  reaction_immobile->zero_time_step();
263 
264 }
265 
267 {
268  eq_fields_->set_time(time_->step(-2), LimitSide::right);
269 
270  START_TIMER("dual_por_exchange_step");
271  reaction_assembly_->assemble(eq_data_->dof_handler_);
272  END_TIMER("dual_por_exchange_step");
273 
274  if(reaction_mobile) reaction_mobile->update_solution();
275  if(reaction_immobile) reaction_immobile->update_solution();
276 }
277 
278 
280 {
281 }
282 
283 
285 {
286  eq_fields_->output_fields.set_time(time_->step(), LimitSide::right);
287 
288  // Register fresh output data
289  eq_fields_->output_fields.output(time_->step());
290 
291  if (time_->tlevel() !=0) {
292  // zero_time_step call zero_time_Step of subreactions which performs its own output
293  if (reaction_mobile) reaction_mobile->output_data();
294  if (reaction_immobile) reaction_immobile->output_data();
295  }
296 }
#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(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
Cell accessor allow iterate over DOF handler cells.
DualPorosity data.
EqData()
Constructor.
DualPorosity fields.
EquationOutput output_fields
Fields indended for output, i.e. all input fields plus those representing solution.
MultiField< 3, FieldValue< 3 >::Scalar > conc_immobile
Calculated concentrations in the immobile zone.
EqFields()
Collect all fields.
MultiField< 3, FieldValue< 3 >::Scalar > init_conc_immobile
Initial concentrations in the immobile zone.
MultiField< 3, FieldValue< 3 >::Scalar > diffusion_rate_immobile
Mass transfer coefficients between mobile and immobile pores.
Field< 3, FieldValue< 3 >::Scalar > porosity_immobile
Immobile porosity field.
Field< 3, FieldValue< 3 >::Scalar > porosity
Porosity field.
GenericAssembly< InitConditionAssemblyDp > * init_condition_assembly_
general assembly objects, hold assembly objects of appropriate dimension
void update_solution(void) override
GenericAssembly< ReactionAssemblyDp > * reaction_assembly_
~DualPorosity(void)
Destructor.
std::shared_ptr< ReactionTerm > reaction_mobile
Reaction running in mobile zone.
void zero_time_step() override
static const Input::Type::Record & get_input_type()
void compute_reaction(const DHCellAccessor &dh_cell) override
Compute reaction on a single element.
void initialize_fields()
Initializes field sets.
std::shared_ptr< EqData > eq_data_
Equation data.
void initialize() override
Prepares the object to usage.
static const int registrar
Registrar of class to factory.
std::shared_ptr< EqFields > eq_fields_
Equation fields - all fields are in this set.
void output_data(void) override
Main output routine.
std::shared_ptr< ReactionTerm > reaction_immobile
Reaction running in immobile zone.
void make_reactions()
Resolves construction of following reactions.
FieldSet input_field_set_
Input::Record input_record_
Definition: equation.hh:242
std::shared_ptr< FieldSet > eq_fieldset_
Definition: equation.hh:249
TimeGovernor * time_
Definition: equation.hh:241
Mesh * mesh_
Definition: equation.hh:240
TimeGovernor & time()
Definition: equation.hh:151
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:95
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
void set_input_list(Input::Array input_list, const TimeGovernor &tg)
Definition: field_set.hh:305
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
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
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
Record type proxy class.
Definition: type_record.hh:182
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
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
Definition: mesh.h:362
std::shared_ptr< OutputTime > output_stream_
Pointer to a transport output stream.
static Input::Type::Abstract & it_abstract_immobile_term()
std::shared_ptr< EqFields > eq_fields_base_
Equation data - all fields needs in assembly class.
std::shared_ptr< EqData > eq_data_base_
Equation data - all data needs in assembly class.
static Input::Type::Abstract & it_abstract_term()
static Input::Type::Abstract & it_abstract_mobile_term()
ReactionTerm(Mesh &init_mesh, Input::Record in_rec)
Sorption model in immobile zone in case dual porosity is considered.
Definition: sorption.hh:132
Sorption model in mobile zone in case dual porosity is considered.
Definition: sorption.hh:106
const TimeStep & step(int index=-1) const
int tlevel() 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 FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
#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.
This file contains classes representing sorption model. Sorption model can be computed both in case t...
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.