Flow123d  master-f44eb46
hm_iterative.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 hm_iterative.cc
15  * @brief
16  * @author Jan Stebel
17  */
18 
19 #include "hm_iterative.hh"
20 #include "system/sys_profiler.hh"
21 #include "input/input_type.hh"
22 #include "flow/darcy_flow_lmh.hh"
23 #include "fields/field_fe.hh" // for create_field_fe()
24 #include "fields/field_model.hh" // for Model
25 #include "assembly_hm.hh"
26 
27 
28 FLOW123D_FORCE_LINK_IN_CHILD(coupling_iterative)
29 
30 
31 namespace it = Input::Type;
32 
33 
35  inline double operator() (double alpha, double density, double gravity, double pressure)
36  {
37  return -alpha*density*gravity*pressure;
38  }
39 };
40 
41 
43 
44  fn_hm_coupling_beta(double beta_f) : beta_factor(beta_f) {}
45 
46 
47  inline double operator() (double alpha, double lame_mu, double lame_lambda, double density, double gravity)
48  {
49  return beta_factor*0.5*alpha*alpha/(2*lame_mu/3 + lame_lambda)*density*gravity;
50  }
51 
52 private:
53 
54  const double beta_factor; ///< User-defined factor for iteration parameter.
55 
56 };
57 
58 
60 
61  fn_fluid_source(const TimeGovernor *time) : time_(time) {}
62 
63  inline double operator() (double alpha, double beta, double pressure, double old_pressure, double div_u, double old_div_u)
64  {
65  return (beta*(pressure-old_pressure) + alpha*(old_div_u - div_u)) / time_->dt();
66  }
67 
68 private:
69 
71 };
72 
73 
75  std::string equation_name = "Coupling_Iterative";
76  return it::Record(equation_name,
77  "Record with data for iterative coupling of flow and mechanics.\n")
82  .declare_key("flow_equation", DarcyLMH::get_input_type(),
84  "Flow equation, provides the velocity field as a result.")
85  .declare_key("mechanics_equation", Elasticity::get_input_type(),
87  "Mechanics, provides the displacement field.")
88  .declare_key("input_fields", it::Array(
90  .make_field_descriptor_type(equation_name)),
92  "Input fields of the HM coupling.")
93  .declare_key( "iteration_parameter", it::Double(), it::Default("1"),
94  "Tuning parameter for iterative splitting. Its default value "
95  "corresponds to a theoretically optimal value with fastest convergence." )
96  .declare_key( "max_it", it::Integer(0), it::Default("100"),
97  "Maximal count of HM iterations." )
98  .declare_key( "min_it", it::Integer(0), it::Default("1"),
99  "Minimal count of HM iterations." )
100  .declare_key( "a_tol", it::Double(0), it::Default("0"),
101  "Absolute tolerance for difference in HM iteration." )
102  .declare_key( "r_tol", it::Double(0), it::Default("1e-7"),
103  "Relative tolerance for difference in HM iteration." )
104  .close();
105 }
106 
107 
108 const int HM_Iterative::registrar = Input::register_class< HM_Iterative, Mesh &, const Input::Record >("Coupling_Iterative")
110 
111 
113 {
114  *this += alpha.name("biot_alpha")
115  .description("Biot poroelastic coefficient.")
116  .units(UnitSI().dimensionless())
117  .input_default("0.0")
119 
120  *this += density.name("fluid_density")
121  .description("Volumetric mass density of the fluid.")
122  .units(UnitSI().kg().m(-3))
123  .input_default("0.0")
125 
126  *this += gravity.name("gravity")
127  .description("Gravitational acceleration constant.")
128  .units(UnitSI().m().s(-2))
129  .input_default("9.81")
131 
132  *this += beta.name("relaxation_beta")
133  .description("Parameter of numerical method for iterative solution of hydro-mechanical coupling.")
134  .units(UnitSI().dimensionless())
136 
137  *this += pressure_potential.name("pressure_potential")
138  .description("Coupling term entering the mechanics equation.")
139  .units(UnitSI().m())
141 
142  *this += old_pressure.name("old_pressure")
143  .description("Pressure from last computed time.")
144  .units(UnitSI().m())
146 
147  *this += old_iter_pressure.name("old_iter_pressure")
148  .description("Pressure from last computed iteration.")
149  .units(UnitSI().m())
151 
152  *this += old_div_u.name("old_displacement_divergence")
153  .description("Displacement divergence from last computed time.")
154  .units(UnitSI().dimensionless())
156 
157  *this += ref_pressure_potential.name("ref_pressure_potential")
158  .description("Pressure potential on boundary (taking into account the flow boundary condition.")
159  .units(UnitSI().m())
161 
162  *this += flow_source.name("extra_flow_source")
163  .description("Coupling term entering the flow equation.")
164  .units(UnitSI().s(-1))
166 
167  this->set_default_fieldset();
168 }
169 
170 
172 {
173  // initialize coupling fields with FieldFE
174  set_mesh(mesh);
175 
176  pressure_potential.set(Model<3, FieldValue<3>::Scalar>::create(
178  alpha,
179  density,
180  gravity,
181  eq_data.flow_->eq_fields().field_edge_pressure
182  ), 0.0);
183 
184  ref_potential_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(mesh, MixedPtr<FE_CR>());
185  ref_pressure_potential.set(ref_potential_ptr_, 0.0);
186 
187  beta.set(Model<3, FieldValue<3>::Scalar>::create(
188  fn_hm_coupling_beta(beta_),
189  alpha,
190  eq_data.mechanics_->eq_fields().lame_mu,
191  eq_data.mechanics_->eq_fields().lame_lambda,
192  density,
193  gravity
194  ), 0.0);
195 
196  auto old_pressure_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(eq_data.flow_->eq_data().dh_cr_, &eq_data.flow_->eq_data().p_edge_solution_previous_time);
197  old_pressure.set(old_pressure_ptr_, 0.0);
198 
199  // TODO: implement a method of FieldFE to easily duplicate the field "field_ele_pressure" to "old_iter_pressure".
200  old_iter_pressure_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(eq_data.flow_->eq_fields().field_ele_pressure.get_field_fe()->get_dofhandler(), nullptr, 1);
201  old_iter_pressure.set(old_iter_pressure_ptr_, 0.0);
202 
203  old_div_u_ptr_ = create_field_fe<3, FieldValue<3>::Scalar>(eq_data.mechanics_->eq_fields().output_div_ptr->get_dofhandler());
204  old_div_u.set(old_div_u_ptr_, 0.0);
205 
206  flow_source.set(Model<3, FieldValue<3>::Scalar>::create(
208  alpha,
209  beta,
210  eq_data.flow_->eq_fields().field_edge_pressure,
211  old_pressure,
212  eq_data.mechanics_->eq_fields().output_divergence,
213  old_div_u
214  ), 0.0);
215 }
216 
217 
218 
220 : DarcyFlowInterface(mesh, in_record),
221  IterativeCoupling(in_record),
222  flow_potential_assembly_(nullptr),
223  residual_assembly_(nullptr)
224 {
225  START_TIMER("HM constructor");
226  using namespace Input;
227 
228  eq_fields_ = std::make_shared<EqFields>();
229  eq_data_ = std::make_shared<EqData>();
230 
231  time_ = new TimeGovernor(in_record.val<Record>("time"));
232  ASSERT( time_->is_default() == false ).error("Missing key 'time' in Coupling_Iterative.");
233 
234  // setup flow equation
235  Record flow_rec = in_record.val<Record>("flow_equation");
236  // Need explicit template types here, since reference is used (automatically passing by value)
237  eq_data_->flow_ = std::make_shared<DarcyLMH>(*mesh_, flow_rec, time_);
238 
239  // setup mechanics
240  Record mech_rec = in_record.val<Record>("mechanics_equation");
241  eq_data_->mechanics_ = std::make_shared<Elasticity>(*mesh_, mech_rec, this->time_);
242  eq_data_->mechanics_->initialize();
243 
244  // setup coupling fields and finish initialization of flow
245  eq_data_->mechanics_->eq_fields()["cross_section"].copy_from(eq_data_->flow_->eq_fields()["cross_section"]);
246  eq_data_->flow_->eq_fields() += eq_data_->mechanics_->eq_fields()["cross_section_updated"];
247  eq_data_->flow_->eq_fields() += eq_data_->mechanics_->eq_fields()["stress"];
248  eq_data_->flow_->eq_fields() += eq_data_->mechanics_->eq_fields()["von_mises_stress"];
249  eq_data_->flow_->eq_fields() += eq_data_->mechanics_->eq_fields()["mean_stress"];
250  eq_data_->flow_->initialize();
251  std::stringstream ss; // print warning message with table of uninitialized fields
252  if ( FieldCommon::print_message_table(ss, "flow") )
253  WarningOut() << ss.str();
254 
255  *eq_fields_ += *eq_data_->flow_->eq_fields().field("pressure_edge");
256 
257  this->eq_fieldset_ = eq_fields_;
258 
259  // setup input fields
260  eq_fields_->set_input_list( in_record.val<Input::Array>("input_fields"), time() );
261 
262  eq_fields_->initialize(*mesh_, *eq_data_, time_, input_record_.val<double>("iteration_parameter"));
263  eq_data_->mechanics_->set_potential_load(eq_fields_->pressure_potential, eq_fields_->ref_pressure_potential);
264 
265  eq_fields_->add_coords_field();
266 }
267 
268 
270 {
273 
274  Input::Array user_fields_arr;
275  if (input_record_.opt_val("user_fields", user_fields_arr)) {
276  FieldSet sham_eq_output; // only for correct call of init_user_fields method
277  this->init_user_fields(user_fields_arr, sham_eq_output);
278  }
279 }
280 
281 
282 template<int dim, class Value>
283 void copy_field(const FieldFE<dim, Value> &from_field, FieldFE<dim, Value> &to_field)
284 {
285  ASSERT( from_field.get_dofhandler() == to_field.get_dofhandler() );
286  to_field.vec().copy_from( from_field.vec() );
287 }
288 
289 
290 
292 {
293  eq_fields_->set_time(time_->step(), LimitSide::right);
294  std::stringstream ss;
295  if ( FieldCommon::print_message_table(ss, "coupling_iterative") )
296  WarningOut() << ss.str();
297 
298  eq_data_->mechanics_->update_output_fields(); // init field values for use in flow
299  eq_data_->flow_->zero_time_step();
301  eq_data_->mechanics_->zero_time_step();
302 
303  copy_field(*eq_data_->flow_->eq_fields().field_ele_pressure.get_field_fe(), *eq_fields_->old_iter_pressure_ptr_);
304  copy_field(*eq_data_->mechanics_->eq_fields().output_divergence.get_field_fe(), *eq_fields_->old_div_u_ptr_);
305  eq_fields_->old_iter_pressure.set_time_result_changed();
306 }
307 
308 
310 {
311  time_->next_time();
312  time_->view("HM");
313  eq_fields_->set_time(time_->step(), LimitSide::right);
314 
315  solve_step();
316 }
317 
319 {
320  // pass displacement (divergence) to flow
321  // and solve flow problem
323  eq_data_->flow_->solve_time_step(false);
324 
325  // pass pressure to mechanics and solve mechanics
327  eq_data_->mechanics_->solve_linear_system();
328 }
329 
330 
332 {
333  eq_data_->mechanics_->update_output_fields();
334  copy_field(*eq_data_->flow_->eq_fields().field_ele_pressure.get_field_fe(), *eq_fields_->old_iter_pressure_ptr_);
335  eq_fields_->old_iter_pressure.set_time_result_changed();
336 }
337 
338 
340 {
341  eq_data_->flow_->accept_time_step();
342  eq_data_->flow_->output_data();
343  eq_data_->mechanics_->output_data();
344 
345  copy_field(*eq_data_->mechanics_->eq_fields().output_divergence.get_field_fe(), *eq_fields_->old_div_u_ptr_);
346 }
347 
348 
350 {
351  auto ref_potential_vec_ = eq_fields_->ref_potential_ptr_->vec();
352  auto dh = eq_fields_->ref_potential_ptr_->get_dofhandler();
353 
354  ref_potential_vec_.zero_entries();
356 
357  ref_potential_vec_.local_to_ghost_begin();
358  ref_potential_vec_.local_to_ghost_end();
359  eq_fields_->pressure_potential.set_time_result_changed();
360  eq_fields_->ref_pressure_potential.set_time_result_changed();
361  eq_data_->mechanics_->set_potential_load(eq_fields_->pressure_potential, eq_fields_->ref_pressure_potential);
362 }
363 
364 
366 {
367  eq_fields_->beta.set_time_result_changed();
368  eq_fields_->flow_source.set_time_result_changed();
369  eq_data_->flow_->set_extra_storativity(eq_fields_->beta);
370  eq_data_->flow_->set_extra_source(eq_fields_->flow_source);
371 }
372 
373 
374 void HM_Iterative::compute_iteration_error(double& abs_error, double& rel_error)
375 {
376  auto dh = eq_fields_->old_iter_pressure_ptr_->get_dofhandler();
377  eq_data_->p_dif2 = 0;
378  eq_data_->p_norm2 = 0;
379 
381 
382  double send_data[] = { eq_data_->p_dif2, eq_data_->p_norm2 };
383  double recv_data[2];
384  MPI_Allreduce(&send_data, &recv_data, 2, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
385  abs_error = sqrt(recv_data[0]);
386  rel_error = abs_error / (sqrt(recv_data[1]) + std::numeric_limits<double>::min());
387 
388  MessageOut().fmt("HM Iteration {} abs. difference: {} rel. difference: {}\n"
389  "--------------------------------------------------------",
390  iteration(), abs_error, rel_error);
391 
392  if(iteration() >= max_it_ && (abs_error > a_tol_ || rel_error > r_tol_))
393  THROW(ExcSolverDiverge() << EI_Reason("Reached max_it."));
394 }
395 
396 
397 
399  eq_data_->flow_.reset();
400  eq_data_->mechanics_.reset();
401  if (flow_potential_assembly_ != nullptr) delete flow_potential_assembly_;
402  if (residual_assembly_ != nullptr) delete residual_assembly_;
403 }
404 
405 
406 
#define ASSERT(expr)
Definition: asserts.hh:351
static Input::Type::Abstract & get_input_type()
static const Input::Type::Record & get_input_type()
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
Definition: elasticity.cc:48
Input::Record input_record_
Definition: equation.hh:242
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
Mesh & mesh()
Definition: equation.hh:181
void init_user_fields(Input::Array user_fields, FieldSet &output_fields)
Definition: equation.cc:84
TimeGovernor * time_
Definition: equation.hh:241
static Input::Type::Record & user_fields_template(std::string equation_name)
Template Record with common key user_fields for derived equations.
Definition: equation.cc:46
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_add(FieldFlag::Flags::Mask mask)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
FieldCommon & input_default(const string &input_default)
std::shared_ptr< DOFHandlerMultiDim > get_dofhandler() const
Definition: field_fe.hh:176
VectorMPI & vec()
Definition: field_fe.hh:180
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 equation_external_output
Match an output field, that can be also copy of other field.
Definition: field_flag.hh:58
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Definition: field_flag.hh:51
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
void set_default_fieldset()
Definition: field_set.hh:408
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
std::shared_ptr< DarcyLMH > flow_
steady or unsteady water flow simulator based on MH scheme
std::shared_ptr< Elasticity > mechanics_
solute transport with chemistry through operator splitting
Field< 3, FieldValue< 3 >::Scalar > flow_source
Field< 3, FieldValue< 3 >::Scalar > density
Density of fluid.
Field< 3, FieldValue< 3 >::Scalar > old_pressure
Field< 3, FieldValue< 3 >::Scalar > old_iter_pressure
Field< 3, FieldValue< 3 >::Scalar > beta
Field< 3, FieldValue< 3 >::Scalar > ref_pressure_potential
Potential of reference (prescribed) pressure from flow b.c. TODO: Swith to BCField when possible.
Field< 3, FieldValue< 3 >::Scalar > alpha
Biot coefficient.
Field< 3, FieldValue< 3 >::Scalar > old_div_u
void initialize(Mesh &mesh, HM_Iterative::EqData &eq_data, const TimeGovernor *time_, double beta_)
Field< 3, FieldValue< 3 >::Scalar > pressure_potential
Potential -alpha*pressure whose gradient is passed to mechanics as additional load.
Field< 3, FieldValue< 3 >::Scalar > gravity
Standard gravity.
void solve_iteration() override
Solve equations and update data (fields).
void update_potential()
GenericAssembly< FlowPotentialAssemblyHM > * flow_potential_assembly_
GenericAssembly< ResidualAssemblyHM > * residual_assembly_
void initialize() override
static const Input::Type::Record & get_input_type()
Define input record.
Definition: hm_iterative.cc:74
static const int registrar
void update_after_iteration() override
Save data (e.g. solution fields) for the next iteration.
std::shared_ptr< EqData > eq_data_
void update_solution() override
void compute_iteration_error(double &abs_error, double &rel_error) override
Compute absolute and relative error in the solution.
void zero_time_step() override
void update_flow_fields()
HM_Iterative(Mesh &mesh, Input::Record in_record)
void update_after_converged() override
Save data after iterations have finished.
std::shared_ptr< EqFields > eq_fields_
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
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
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
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
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
double r_tol_
Relative tolerance for difference between two succeeding iterations.
double a_tol_
Absolute tolerance for difference between two succeeding iterations.
unsigned int max_it_
Maximal number of iterations.
unsigned int iteration()
Definition: hm_iterative.hh:88
static const Input::Type::Record & record_template()
Definition: hm_iterative.hh:49
Definition: mesh.h:362
Basic time management functionality for unsteady (and steady) solvers (class Equation).
double dt() const
void view(const char *name="") const
const TimeStep & step(int index=-1) const
void next_time()
Proceed to the next time according to current estimated time step.
Class for representation SI units of Fields.
Definition: unit_si.hh:40
void copy_from(const VectorMPI &other)
Definition: vector_mpi.cc:103
Lumped mixed-hybrid model of linear Darcy flow, possibly unsteady.
double lame_mu(double young, double poisson)
Definition: elasticity.cc:80
double lame_lambda(double young, double poisson)
Definition: elasticity.cc:86
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
void copy_field(const FieldFE< dim, Value > &from_field, FieldFE< dim, Value > &to_field)
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define MessageOut()
Macro defining 'message' record of log.
Definition: logger.hh:275
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
Definition: mpi.h:612
#define MPI_SUM
Definition: mpi.h:196
#define MPI_DOUBLE
Definition: mpi.h:156
Abstract linear system class.
Definition: balance.hh:40
fn_fluid_source(const TimeGovernor *time)
Definition: hm_iterative.cc:61
double operator()(double alpha, double beta, double pressure, double old_pressure, double div_u, double old_div_u)
Definition: hm_iterative.cc:63
const TimeGovernor * time_
Definition: hm_iterative.cc:70
fn_hm_coupling_beta(double beta_f)
Definition: hm_iterative.cc:44
double operator()(double alpha, double lame_mu, double lame_lambda, double density, double gravity)
Definition: hm_iterative.cc:47
const double beta_factor
User-defined factor for iteration parameter.
Definition: hm_iterative.cc:54
double operator()(double alpha, double density, double gravity, double pressure)
Definition: hm_iterative.cc:35
#define START_TIMER(tag)
Starts a timer with specified tag.