Flow123d  master-f44eb46
hm_iterative.hh
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.hh
15  * @brief
16  * @author Jan Stebel
17  */
18 
19 #ifndef HM_ITERATIVE_HH_
20 #define HM_ITERATIVE_HH_
21 
22 #include <memory>
23 #include <string>
24 #include <vector>
27 #include "coupling/equation.hh"
29 #include "mechanics/elasticity.hh"
30 #include "system/exceptions.hh"
31 
32 class Mesh;
33 class FieldCommon;
34 class DarcyLMH;
35 
36 template<unsigned int dim> class FlowPotentialAssemblyHM;
37 template<unsigned int dim> class ResidualAssemblyHM;
38 
39 namespace it = Input::Type;
40 
41 
43 public:
44  TYPEDEF_ERR_INFO( EI_Reason, string);
45  DECLARE_EXCEPTION(ExcSolverDiverge,
46  << "Nonlinear solver did not converge. Reason: " << EI_Reason::val
47  );
48 
50  return it::Record("Coupling_Iterative_AUX",
51  "Record with data for iterative coupling.\n")
52  .declare_key( "max_it", it::Integer(0), it::Default("100"),
53  "Maximal count of HM iterations." )
54  .declare_key( "min_it", it::Integer(0), it::Default("1"),
55  "Minimal count of HM iterations." )
56  .declare_key( "a_tol", it::Double(0), it::Default("0"),
57  "Absolute tolerance for difference in HM iteration." )
58  .declare_key( "r_tol", it::Double(0), it::Default("1e-7"),
59  "Relative tolerance for difference in HM iteration." )
60  .close();
61  }
62 
64  : it(0)
65  {
66  min_it_ = in_record.val<unsigned int>("min_it");
67  max_it_ = in_record.val<unsigned int>("max_it");
68  a_tol_ = in_record.val<double>("a_tol");
69  r_tol_ = in_record.val<double>("r_tol");
70  }
71 
72  void solve_step()
73  {
74  it = 0;
75  double abs_error = std::numeric_limits<double>::max();
76  double rel_error = std::numeric_limits<double>::max();
77 
78  while ( it < min_it_ || (abs_error > a_tol_ && rel_error > r_tol_ && it < max_it_) )
79  {
80  it++;
82  compute_iteration_error(abs_error, rel_error);
84  }
86  }
87 
88  unsigned int iteration()
89  { return it; }
90 
91 protected:
92 
93  /// Solve equations and update data (fields).
94  virtual void solve_iteration() = 0;
95 
96  /// Compute absolute and relative error in the solution.
97  virtual void compute_iteration_error(double &abs_error, double &rel_error) = 0;
98 
99  /// Save data (e.g. solution fields) for the next iteration.
100  virtual void update_after_iteration() = 0;
101 
102  /// Save data after iterations have finished.
103  virtual void update_after_converged() = 0;
104 
105 
106  /// Minimal number of iterations to perform.
107  unsigned int min_it_;
108 
109  /// Maximal number of iterations.
110  unsigned int max_it_;
111 
112  /// Absolute tolerance for difference between two succeeding iterations.
113  double a_tol_;
114 
115  /// Relative tolerance for difference between two succeeding iterations.
116  double r_tol_;
117 
118 private:
119 
120  /// Iteration index.
121  unsigned int it;
122 
123 };
124 
125 
126 /**
127  * @brief Class for solution of fully coupled flow and mechanics using fixed-stress iterative splitting.
128  *
129  * Flow and mechanics are solved separately and within each iteration the coupling terms are updated.
130  * Here we use the fixed-stress splitting [see Mikelic&Wheeler, Comput. Geosci. 17(3), 2013] which uses
131  * a tuning parameter "beta" to speed up the convergence.
132  */
134 public:
135 
136  class EqData
137  {
138  public:
139  /// steady or unsteady water flow simulator based on MH scheme
140  std::shared_ptr<DarcyLMH> flow_;
141 
142  /// solute transport with chemistry through operator splitting
143  std::shared_ptr<Elasticity> mechanics_;
144 
145  double p_dif2; ///< Squared norm of pressure difference in two subsequent iterations.
146  double p_norm2; ///< Squared pressure norm in the last iteration.
147  };
148 
149 
150  class EqFields : public FieldSet
151  {
152  public:
153  EqFields();
154 
155  void initialize(Mesh &mesh, HM_Iterative::EqData &eq_data, const TimeGovernor *time_, double beta_);
156 
157  Field<3, FieldValue<3>::Scalar> alpha; ///< Biot coefficient.
158  Field<3, FieldValue<3>::Scalar> density; ///< Density of fluid.
159  Field<3, FieldValue<3>::Scalar> gravity; ///< Standard gravity.
161 
162  /// Potential -alpha*pressure whose gradient is passed to mechanics as additional load.
164  Field<3, FieldValue<3>::Scalar> ref_pressure_potential; ///< Potential of reference (prescribed) pressure from flow b.c. TODO: Swith to BCField when possible.
169 
170  /// FieldFE for pressure_potential field.
171  std::shared_ptr<FieldFE<3, FieldValue<3>::Scalar> > ref_potential_ptr_;
172  std::shared_ptr<FieldFE<3, FieldValue<3>::Scalar> > old_iter_pressure_ptr_;
173  std::shared_ptr<FieldFE<3, FieldValue<3>::Scalar> > old_div_u_ptr_;
174  };
175 
176  /// Define input record.
177  static const Input::Type::Record & get_input_type();
178 
179  HM_Iterative(Mesh &mesh, Input::Record in_record);
180  void initialize() override;
181  void zero_time_step() override;
182  void update_solution() override;
183  ~HM_Iterative();
184 
185 private:
186 
187  void update_potential();
188 
189  void update_flow_fields();
190 
191  void solve_iteration() override;
192 
193  void update_after_iteration() override;
194 
195  void update_after_converged() override;
196 
197  void compute_iteration_error(double &abs_error, double &rel_error) override;
198 
199  static const int registrar;
200 
203 
204  std::shared_ptr<EqFields> eq_fields_;
205 
206  std::shared_ptr<EqData> eq_data_;
207 
208 };
209 
210 #endif /* HC_EXPLICIT_SEQUENTIAL_HH_ */
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
Mesh & mesh()
Definition: equation.hh:181
TimeGovernor * time_
Definition: equation.hh:241
Common abstract parent of all Field<...> classes.
Definition: field_common.hh:77
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
const Mesh * mesh() const
Returns pointer to mesh.
Definition: field_set.hh:393
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
double p_dif2
Squared norm of pressure difference in two subsequent iterations.
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
double p_norm2
Squared pressure norm in the last iteration.
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > old_iter_pressure_ptr_
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.
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > ref_potential_ptr_
FieldFE for pressure_potential field.
Field< 3, FieldValue< 3 >::Scalar > alpha
Biot coefficient.
Field< 3, FieldValue< 3 >::Scalar > old_div_u
std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > old_div_u_ptr_
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.
Class for solution of fully coupled flow and mechanics using fixed-stress iterative splitting.
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 the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
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
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
unsigned int it
Iteration index.
double r_tol_
Relative tolerance for difference between two succeeding iterations.
virtual void solve_iteration()=0
Solve equations and update data (fields).
double a_tol_
Absolute tolerance for difference between two succeeding iterations.
DECLARE_EXCEPTION(ExcSolverDiverge,<< "Nonlinear solver did not converge. Reason: "<< EI_Reason::val)
unsigned int min_it_
Minimal number of iterations to perform.
virtual void compute_iteration_error(double &abs_error, double &rel_error)=0
Compute absolute and relative error in the solution.
unsigned int max_it_
Maximal number of iterations.
unsigned int iteration()
Definition: hm_iterative.hh:88
virtual void update_after_iteration()=0
Save data (e.g. solution fields) for the next iteration.
IterativeCoupling(Input::Record in_record)
Definition: hm_iterative.hh:63
TYPEDEF_ERR_INFO(EI_Reason, string)
static const Input::Type::Record & record_template()
Definition: hm_iterative.hh:49
virtual void update_after_converged()=0
Save data after iterations have finished.
Definition: mesh.h:362
Basic time management functionality for unsteady (and steady) solvers (class Equation).
FEM for linear elasticity.
Abstract base class for equation clasess.