Flow123d  master-f44eb46
hc_explicit_sequential.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 hc_explicit_sequential.cc
15  * @brief
16  * @author Jan Brezina
17  */
18 
21 // TODO:
22 // After having general default values:
23 // make TransportNoting default for AdvectionProcessBase abstract
24 // use default "{}" for secondary equation.
25 // Then we can remove following include.
26 
28 #include "fields/field_common.hh"
29 #include "transport/heat_model.hh"
30 
31 #include "fields/field_set.hh"
32 #include "mesh/mesh.h"
33 #include "io/msh_gmshreader.h"
34 #include "system/sys_profiler.hh"
35 #include "input/input_type.hh"
36 #include "input/accessors.hh"
37 
38 
39 FLOW123D_FORCE_LINK_IN_PARENT(transportOperatorSplitting)
40 FLOW123D_FORCE_LINK_IN_PARENT(concentrationTransportModel)
41 FLOW123D_FORCE_LINK_IN_PARENT(convectionTransport)
43 
44 FLOW123D_FORCE_LINK_IN_PARENT(darcy_flow_lmh)
46 FLOW123D_FORCE_LINK_IN_PARENT(coupling_iterative)
47 
48 
49 namespace it = Input::Type;
50 
52  return it::Abstract("Coupling_Base", "The root record of description of particular the problem to solve.")
53  .close();
54 }
55 
56 
58  return it::Record("Coupling_Sequential",
59  "Record with data for a general sequential coupling.\n")
62  .declare_key("description",it::String(),
63  "Short description of the solved problem.\n"
64  "Is displayed in the main log, and possibly in other text output files.")
66  "Computational mesh common to all equations.")
69  "Flow equation, provides the velocity field as a result.")
71  "Transport of soluted substances, depends on the velocity field from a Flow equation.")
73  "Heat transfer, depends on the velocity field from a Flow equation.")
74  .close();
75 }
76 
77 
79 
80 
81 
82 std::shared_ptr<AdvectionProcessBase> HC_ExplicitSequential::make_advection_process(string process_key)
83 {
84  using namespace Input;
85  // setup heat object
87 
88  if (it) {
89  auto process = (*it).factory< AdvectionProcessBase, Mesh &, const Input::Record >(*mesh, *it);
90 
91  // setup fields
92  process->eq_fieldset()["cross_section"]
93  .copy_from(water->eq_fieldset()["cross_section"]);
94  /*
95  if (water_content_saturated_) // only for unsteady Richards water model
96  process->data()["porosity"].copy_from(*water_content_saturated_);
97 
98  if (water_content_p0_)
99  process->data()["water_content"].copy_from(*water_content_p0_);
100  else {
101 
102  }*/
103 
104  FieldCommon *porosity = process->eq_fieldset().field("porosity");
105  process->eq_fieldset()["water_content"].copy_from( *porosity );
106 
107 
108  process->initialize();
109  return process;
110  } else {
111  return std::make_shared<TransportNothing>(*mesh);
112  }
113 }
114 
115 /**
116  * FUNCTION "MAIN" FOR COMPUTING MIXED-HYBRID PROBLEM FOR UNSTEADY SATURATED FLOW
117  */
119 : in_record_(in_record)
120 {
121  START_TIMER("HC constructor");
122  using namespace Input;
123 
124  // Read mesh
125  {
126  START_TIMER("HC read mesh");
127 
128  mesh = BaseMeshReader::mesh_factory( in_record.val<Record>("mesh") );
129 
130  //getting description for the Profiler
131  string description;
132  in_record.opt_val<string>("description", description);
133 
135  description,
136  //"Description has to be set in main. by different method.",
137  mesh->n_elements());
138  }
139 
140  // setup primary equation - water flow object
141  AbstractRecord prim_eq = in_record.val<AbstractRecord>("flow_equation");
142  // Need explicit template types here, since reference is used (automatically passing by value)
143  water = prim_eq.factory< DarcyFlowInterface, Mesh &, const Input::Record>(*mesh, prim_eq);
144  water->initialize();
145  std::stringstream ss; // print warning message with table of uninitialized fields
146  if ( FieldCommon::print_message_table(ss, "HC explicit sequential") ) {
147  WarningOut() << ss.str();
148  }
149 
150  RegionSet bulk_set = mesh->region_db().get_region_set("BULK");
151  water_content_saturated_ = water->eq_fieldset().field("water_content_saturated");
153  water_content_saturated_ = nullptr;
154 
155  water_content_p0_ = water->eq_fieldset().field("water_content_p0");
157  water_content_p0_ = nullptr;
158 
159  processes_.push_back(AdvectionData(make_advection_process("solute_equation")));
160  processes_.push_back(AdvectionData(make_advection_process("heat_equation")));
161 }
162 
164 {
165  if (pdata.process->time().is_end()) return;
166 
167  is_end_all_=false;
168  if ( pdata.process->time().step().le(pdata.velocity_time) ) {
169  // having information about velocity field we can perform transport step
170 
171  // here should be interpolation of the velocity at least if the interpolation time
172  // is not close to the solved_time of the water module
173  // for simplicity we use only last velocity field
174  if (pdata.velocity_changed) {
175  // DebugOut() << "flow_flux update\n";
176  auto& flux = pdata.process->eq_fieldset()["flow_flux"];
177  flux.copy_from(water->eq_fieldset()["flux"]);
178  flux.set_time_result_changed();
179  pdata.velocity_changed = false;
180  }
181  if (pdata.process->time().tlevel() == 0) pdata.process->zero_time_step();
182 
183  pdata.process->update_solution();
184  pdata.process->output_data();
185  }
186 }
187 
188 void HC_ExplicitSequential::flow_step(double requested_time)
189 {
190  // DebugOut() << "water->solved_time() = " << water->solved_time() << " requested_time = " << requested_time << "\n";
191  if (water->solved_time() < requested_time) {
192  // solve water over the requested_time (nearest transport interval)
193  water->update_solution();
194 
195  // here possibly save solution from water for interpolation in time
196 
197  //water->time().view("WATER"); //show water time governor
198  water->choose_next_time();
199 
200  for(auto &pdata : processes_)
201  pdata.velocity_changed = true;
202  }
203 }
204 
205 /**
206  * TODO:
207  * - have support for steady problems in TimeGovernor, make Noting problems steady
208  * - apply splitting of compute_one_step to particular models
209  * - how to set output time marks for steady problems (we need solved time == infinity) but
210  * add no time marks
211  * - allow create steady time governor without time marks (at least in nothing models)
212  * - pass refference to time marks in time governor constructor?
213  */
214 
216 {
217  START_TIMER("HC run simulation");
218  // following should be specified in constructor:
219  // value for velocity interpolation :
220  // theta = 0 velocity from beginning of transport interval (fully explicit method)
221  // theta = 0.5 velocity from center of transport interval ( mimic Crank-Nicholson)
222  // theta = 1.0 velocity from end of transport interval (partialy explicit scheme)
223  const double theta=0.5;
224 
225  {
226  START_TIMER("HC water zero time step");
227  water->zero_time_step();
228  for(auto &process : processes_)
229  process.velocity_changed = true;
230  }
231 
232 
233  // following cycle is designed to support independent time stepping of
234  // both processes. The question is which value of the water field use to compute a transport step.
235  // Meaningful cases are
236  // 1) beginning (fully explicit method)
237  // 2) center ( mimic Crank-Nicholson)
238  // 3) end of the interval (partialy explicit scheme)
239  // However with current implementation of the explicit transport on have to assembly transport matrix for
240  // every new value of the velocity field. So we have to keep same velocity field over some time interval t_dt
241  // which is further split into shorter time intervals ts_dt dictated by the CFL condition.
242  // One can consider t_dt as the transport time step and apply one of the previous three cases.
243  //
244  // The question is how to choose intervals t_dt. That should depend on variability of the velocity field in time.
245  // Currently we simply use t_dt == w_dt.
246 
247  is_end_all_=false;
248  while (! is_end_all_) {
249  is_end_all_ = true;
250 
251  // if any of the advection processes reached the current time of flow model,
252  // take a new flow step before setting a new time constraint for adv. process
253  for(auto &pdata : processes_){
254  if(! pdata.process->time().is_end())
255  flow_step(pdata.process->solved_time());
256  }
257 
258  double water_dt=water->time().estimate_dt();
259  if (water->time().is_end()) water_dt = TimeGovernor::inf_time;
260 
261  // in future here could be re-estimation of transport planed time according to
262  // evolution of the velocity field. Consider the case w_dt << t_dt and velocity almost constant in time
263  // which suddenly rise in time 3*w_dt. First we the planed transport time step t_dt could be quite big, but
264  // in time 3*w_dt we can reconsider value of t_dt to better capture changing velocity.
266  for(auto &pdata : processes_) {
267  if(! pdata.process->time().is_end()){
268  pdata.process->set_time_upper_constraint(water_dt, "Flow time step");
269  pdata.velocity_time = theta * pdata.process->planned_time() + (1-theta) * pdata.process->solved_time();
270  min_velocity_time = min(min_velocity_time, pdata.velocity_time);
271  }
272  }
273  // DebugOut() << "min_velocity_time = " << min_velocity_time << "\n";
274 
275  // printing water and transport times every step
276  //MessageOut().fmt("HC_EXPL_SEQ: velocity_interpolation_time: {}, water_time: {} transport time: {}\n",
277  // velocity_interpolation_time, water->time().t(), transport_reaction->time().t());
278 
279  // if transport is off, transport should return infinity solved and planned times so that
280  // only water branch takes the place
281 
282  if (! water->time().is_end() ) {
283  is_end_all_=false;
285  }
286  advection_process_step(processes_[0]); // solute
288  }
289  //MessageOut().fmt("End of simulation at time: {}\n", max(solute->solved_time(), heat->solved_time()));
290 }
291 
292 
294  water.reset();
295  for(auto &pdata : processes_) pdata.process.reset();
296  delete mesh;
297 }
298 
299 
300 
301 
302 //-----------------------------------------------------------------------------
303 // vim: set cindent:
304 //-----------------------------------------------------------------------------
static Input::Type::Abstract & get_input_type()
Common specification of the input record for secondary equations.
static Mesh * mesh_factory(const Input::Record &input_mesh_rec)
static Input::Type::Abstract & get_input_type()
static Input::Type::Abstract & get_input_type()
static Input::Type::Record & record_template()
Template Record with common keys for derived equations.
Definition: equation.cc:39
Common abstract parent of all Field<...> classes.
Definition: field_common.hh:77
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:95
virtual FieldResult field_result(RegionSet region_set) const =0
Indicates special field states.
virtual void copy_from(const FieldCommon &other)=0
void advection_process_step(AdvectionData &pdata)
std::vector< AdvectionData > processes_
solute transport with chemistry through operator splitting
static const Input::Type::Record & get_input_type()
Mesh * mesh
mesh common to darcy flow and transport
AdvectionPtr make_advection_process(std::string process_key)
HC_ExplicitSequential(Input::Record in_record)
void flow_step(double requested_time)
std::shared_ptr< DarcyFlowInterface > water
steady or unsteady water flow simulator based on MH scheme
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 the data with type Type::Record.
Definition: accessors.hh:291
bool opt_val(const string &key, Ret &value) const
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).
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
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 & 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
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
const RegionDB & region_db() const
Definition: mesh.h:175
unsigned int n_elements() const
Definition: mesh.h:111
Definition: mesh.h:362
static const Input::Type::Record & get_input_type()
Definition: mesh.cc:213
static Profiler * instance(bool clear=false)
void set_task_info(string, int)
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
static const double inf_time
Infinity time used for steady case.
@ result_zeros
#define FLOW123D_FORCE_LINK_IN_PARENT(x)
Definition: global_defs.h:107
Discontinuous Galerkin method for equation of transport with dispersion.
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
Abstract linear system class.
Definition: balance.hh:40
#define START_TIMER(tag)
Starts a timer with specified tag.