Flow123d  master-1fea4ce
transport_operator_splitting.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 transport_operator_splitting.cc
15  * @brief
16  */
17 
18 #include <iostream>
19 #include <iomanip>
20 
21 #include "system/system.hh"
22 #include "system/sys_profiler.hh"
23 #include "system/index_types.hh"
24 
26 #include <petscmat.h>
27 
28 #include "io/output_time.hh"
29 #include "tools/time_governor.hh"
30 #include "coupling/equation.hh"
31 #include "coupling/balance.hh"
32 #include "transport/transport.h"
33 #include "mesh/mesh.h"
34 
38 #include "reaction/sorption.hh"
40 
41 #include "la/distribution.hh"
42 #include "input/input_type.hh"
43 #include "input/accessors.hh"
44 #include "input/factory.hh"
46 #include "fem/fe_p.hh"
47 
48 FLOW123D_FORCE_LINK_IN_CHILD(transportOperatorSplitting)
49 
50 FLOW123D_FORCE_LINK_IN_PARENT(firstOrderReaction)
51 FLOW123D_FORCE_LINK_IN_PARENT(radioactiveDecay)
53 FLOW123D_FORCE_LINK_IN_PARENT(sorptionMobile)
54 FLOW123D_FORCE_LINK_IN_PARENT(sorptionImmobile)
56 
57 
58 using namespace Input::Type;
59 
60 
61 
63  return Abstract("Solute",
64  "Transport of soluted substances.")
65  .close();
66 }
67 
68 
70  return Record("Coupling_OperatorSplitting",
71  "Transport by convection and/or diffusion\n"
72  "coupled with reaction and adsorption model (ODE per element)\n"
73  " via operator splitting.")
75  .add_attribute( FlowAttribute::subfields_address(), "\"/problem/solute_equation/substances/*/name\"")
77  .declare_key("balance", Balance::get_input_type(), Default("{}"),
78  "Settings for computing mass balance.")
79  .declare_key("output_stream", OutputTime::get_input_type(), Default("{}"),
80  "Output stream settings.\n Specify file format, precision etc.")
82  "Specification of transported substances.")
83  // input data
85  "Type of the numerical method for the transport equation.")
87  "Reaction model involved in transport.")
88 /*
89  .declare_key("output_fields", Array(ConvectionTransport::get_output_selection()),
90  Default("\"conc\""),
91  "List of fields to write to output file.")*/
92  .close();
93 }
94 
95 
97  Input::register_class< TransportOperatorSplitting, Mesh &, const Input::Record>("Coupling_OperatorSplitting") +
99 
100 
101 
102 
103 
104 
105 
107 {
108  *this += porosity.name("porosity")
109  .description("Porosity of the mobile phase.")
110  .input_default("1.0")
111  .units( UnitSI::dimensionless() )
112  .flags_add(in_main_matrix & in_rhs);
113 
114  *this += water_content.name("water_content")
115  .description("INTERNAL. Water content passed from unsaturated Darcy flow model.")
116  .units( UnitSI::dimensionless() )
117  .flags_add(input_copy & in_time_term & in_main_matrix & in_rhs);
118 
119  *this += cross_section.name("cross_section")
120  .flags( FieldFlag::input_copy )
121  .flags_add(in_time_term & in_main_matrix & in_rhs);
122 
123  *this += flow_flux.name("flow_flux")
124  .flags( FieldFlag::input_copy )
125  .flags_add(in_time_term & in_main_matrix & in_rhs);
126 
127  *this += sources_density.name("sources_density")
128  .description("Density of concentration sources.")
129  .input_default("0.0")
130  .units( UnitSI().kg().m(-3).s(-1) )
131  .flags_add(in_rhs);
132 
133  *this += sources_sigma.name("sources_sigma")
134  .description("Concentration flux.")
135  .input_default("0.0")
136  .units( UnitSI().s(-1) )
137  .flags_add(in_main_matrix & in_rhs);
138 
139  *this += sources_conc.name("sources_conc")
140  .description("Concentration sources threshold.")
141  .input_default("0.0")
142  .units( UnitSI().kg().m(-3) )
143  .flags_add(in_rhs);
144 }
145 
146 
147 
148 
149 
150 
151 
152 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
153 
154 
155 
157 : AdvectionProcessBase(init_mesh, in_rec),
158  cfl_convection(numeric_limits<double>::max())
159 {
160  START_TIMER("TransportOperatorSpliting");
161 
162  Input::AbstractRecord trans = in_rec.val<Input::AbstractRecord>("transport");
163  convection = trans.factory< ConcentrationTransportBase, Mesh &, const Input::Record >(init_mesh, trans);
164 
165  time_ = new TimeGovernor(in_rec.val<Input::Record>("time"), TimeMark::none_type, false);
166  ASSERT( time_->is_default() == false ).error("Missing key 'time' in Coupling_OperatorSplitting.");
167  convection->set_time_governor(time());
168 
169  // Initialize list of substances.
170  convection->substances().initialize(in_rec.val<Input::Array>("substances"));
171 
172  // Initialize output stream.
173  std::shared_ptr<OutputTime> stream = OutputTime::create_output_stream("solute", in_rec.val<Input::Record>("output_stream"), time().get_unit_conversion());
174  convection->set_output_stream(stream);
175 
176 
177  // initialization of balance object
178 
179  balance_ = make_shared<Balance>("mass", mesh_);
180  balance_->init_from_input(in_rec.val<Input::Record>("balance"), this->time());
181 
182  balance_->units(UnitSI().kg(1));
183  convection->set_balance_object(balance_);
184 
185  convection->initialize(); //
186 
187 
188 
189 
190  time_ = new TimeGovernor(in_rec.val<Input::Record>("time"), convection->mark_type());
191 
192  this->eq_fieldset_ = convection->eq_fieldset_ptr();
193 
194  Input::Iterator<Input::AbstractRecord> reactions_it = in_rec.find<Input::AbstractRecord>("reaction_term");
195  if ( reactions_it ) {
196  // TODO: allowed instances in this case are only
197  // FirstOrderReaction, RadioactiveDecay, SorptionSimple and DualPorosity
198  reaction = (*reactions_it).factory< ReactionTerm, Mesh &, Input::Record >(init_mesh, *reactions_it);
199 
200  //initialization of DOF handler
201  static MixedPtr<FE_P_disc> fe(0);
202  shared_ptr<DOFHandlerMultiDim> dof_handler = make_shared<DOFHandlerMultiDim>(*mesh_);
203  shared_ptr<DiscreteSpace> ds = make_shared<EqualOrderDiscreteSpace>( mesh_, fe);
204  dof_handler->distribute_dofs(ds);
205 
206  reaction->substances(convection->substances())
207  .concentration_fields(convection->get_p0_interpolation())
208  .output_stream(stream)
209  .set_time_governor((TimeGovernor &)convection->time());
210 
211  reaction->initialize();
212 
213  } else {
214  reaction = nullptr;
215  //Semchem_reactions = nullptr;
216  }
217 }
218 
220 {
221  //delete field_output;
222  //if (Semchem_reactions) delete Semchem_reactions;
223  delete time_;
224 }
225 
226 
228 {
229  //coupling - passing fields
230  if(reaction)
231  if( typeid(*reaction) == typeid(SorptionSimple) ||
232  typeid(*reaction) == typeid(DualPorosity)
233  )
234  {
235  reaction->eq_fieldset().set_field("porosity", convection->eq_fieldset()["porosity"]);
236  }
237 }
238 
239 
240 
241 
243 
244 
245  START_TIMER("TOS-output data");
246 
247 
248  convection->output_data();
249  if(reaction) reaction->output_data(); // do not perform write_time_frame
250 
251  END_TIMER("TOS-output data");
252 }
253 
254 
256 {
257  //DebugOut() << "tos ZERO TIME STEP.\n";
258  convection->zero_time_step();
259  convection->compute_p0_interpolation(); // due to reading of init_conc in reactions
260  if(reaction)
261  {
262  reaction->zero_time_step();
263  reaction->output_data(); // do not perform write_time_frame
264  }
265 
266 }
267 
268 
269 
271 
272  vector<double> source(convection->n_substances()), region_mass(mesh_->region_db().bulk_size());
273 
274  time_->next_time();
275  time_->view("TOS"); //show time governor
276 
277  convection->set_target_time(time_->t());
278 
279  START_TIMER("TOS-one step");
280  int steps=0;
281  while ( convection->time().step().lt(time_->t()) )
282  {
283  steps++;
284  // one internal step
285  bool cfl_changed = convection->evaluate_time_constraint(cfl_convection);
286 
287  if (steps == 1 || cfl_changed)
288  {
289  convection->time().set_upper_constraint(cfl_convection, "Time step constrained by transport CFL condition (including both flow and sources).");
290 
291  // fix step with new constraint
292  convection->time().fix_dt_until_mark();
293 
294  convection->time().view("Convection"); // write TG only once on change
295  }
296 
297  convection->update_solution();
298 
299 
300  if (balance_->cumulative())
301  {
302  START_TIMER("TOS-balance");
303 
304  // save mass after transport step
305  for (unsigned int sbi=0; sbi<convection->n_substances(); sbi++)
306  {
307  balance_->calculate_mass(convection->get_subst_idx()[sbi], convection->get_component_vec(sbi), region_mass);
308  source[sbi] = 0;
309  for (unsigned int ri=0; ri<mesh_->region_db().bulk_size(); ri++)
310  source[sbi] -= region_mass[ri];
311  }
312 
313  END_TIMER("TOS-balance");
314  }
315 
316  if(reaction) {
317  convection->compute_p0_interpolation();
318  reaction->update_solution();
319  convection->update_after_reactions(true);
320  }
321  else
322  convection->update_after_reactions(false);
323 
324  //if(Semchem_reactions) Semchem_reactions->update_solution();
325 
326 
327 
328  if (balance_->cumulative())
329  {
330  START_TIMER("TOS-balance");
331 
332  for (unsigned int sbi=0; sbi<convection->n_substances(); sbi++)
333  {
334  // compute mass difference due to reactions
335  balance_->calculate_mass(convection->get_subst_idx()[sbi], convection->get_component_vec(sbi), region_mass);
336  for (unsigned int ri=0; ri<mesh_->region_db().bulk_size(); ri++)
337  source[sbi] += region_mass[ri];
338 
339  // update balance of sources due to reactions
340  balance_->add_cumulative_source(sbi, source[sbi]);
341  }
342 
343  END_TIMER("TOS-balance");
344  }
345  }
346 
347  LogOut().fmt("CONVECTION: steps: {}\n", steps);
348 }
#define ASSERT(expr)
Definition: asserts.hh:351
static Input::Type::Abstract & get_input_type()
Common specification of the input record for secondary equations.
static const Input::Type::Record & get_input_type()
Main balance input record type.
Definition: balance.cc:53
static Input::Type::Abstract & get_input_type()
Common specification of the input record for secondary equations.
Class representing dual porosity model in transport.
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
Definition: equation.cc:39
std::shared_ptr< FieldSet > eq_fieldset_
Definition: equation.hh:249
TimeGovernor * time_
Definition: equation.hh:241
std::shared_ptr< Balance > balance_
object for calculation and writing the mass balance to file.
Definition: equation.hh:252
Mesh * mesh_
Definition: equation.hh:240
TimeGovernor & time()
Definition: equation.hh:151
static constexpr Mask input_copy
Definition: field_flag.hh:44
static string subfields_address()
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
const std::shared_ptr< Type > factory(Arguments... arguments) const
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 polymorphic Record.
Abstract & close()
Close the Abstract and add its to type repository (see TypeRepository::add_type).
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
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
Record & add_attribute(std::string key, TypeBase::json_string value)
Add TYPE key as obligatory.
Definition: type_record.cc:354
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 & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
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
const RegionDB & region_db() const
Definition: mesh.h:175
Definition: mesh.h:362
static const Input::Type::Record & get_input_type()
The specification of output stream.
Definition: output_time.cc:38
static std::shared_ptr< OutputTime > create_output_stream(const std::string &equation_name, const Input::Record &in_rec, const std::shared_ptr< TimeUnitConversion > &time_unit_conv)
This method delete all object instances of class OutputTime stored in output_streams vector.
Definition: output_time.cc:187
static Input::Type::Abstract & it_abstract_term()
unsigned int bulk_size() const
Definition: region.cc:268
Simple sorption model without dual porosity.
Definition: sorption.hh:48
static const Input::Type::Record & get_input_type()
Input type for a substance.
Definition: substance.cc:29
Basic time management functionality for unsteady (and steady) solvers (class Equation).
double t() const
void view(const char *name="") const
void next_time()
Proceed to the next time according to current estimated time step.
static const Type none_type
Mark Type with all bits unset.
Definition: time_marks.hh:95
void output_data() override
Write computed fields.
virtual ~TransportOperatorSplitting()
Destructor.
double cfl_convection
Time restriction due to transport.
static const int registrar
Registrar of class to factory.
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportOperatorSplittiong.
std::shared_ptr< ConcentrationTransportBase > convection
TransportOperatorSplitting(Mesh &init_mesh, const Input::Record in_rec)
Constructor.
std::shared_ptr< ReactionTerm > reaction
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.
Abstract base class for equation clasess.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
#define FLOW123D_FORCE_LINK_IN_PARENT(x)
Definition: global_defs.h:107
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
#define LogOut()
Macro defining 'log' record of log.
Definition: logger.hh:281
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.
Basic time management class.