Flow123d  JS_before_hm-62-ge56f9d5
richards_lmh.cc
Go to the documentation of this file.
1 /*
2  * richards_lmh.cc
3  *
4  * Created on: Sep 16, 2015
5  * Author: jb
6  */
7 
8 
9 #include "system/global_defs.h"
10 #include "system/sys_profiler.hh"
11 #include "system/asserts.hh"
12 
13 
14 #include "input/input_type.hh"
15 #include "input/factory.hh"
16 #include "flow/richards_lmh.hh"
17 #include "flow/assembly_lmh.hh"
19 #include "tools/time_governor.hh"
20 
21 #include "petscmat.h"
22 #include "petscviewer.h"
23 #include "petscerror.h"
24 #include <armadillo>
25 
26 #include "la/schur.hh"
27 
28 #include "coupling/balance.hh"
29 
30 #include "la/vector_mpi.hh"
31 
32 
33 #include "tools/include_fadbad.hh" // for "fadbad.h", "badiff.h", "fadiff.h"
34 
35 
36 FLOW123D_FORCE_LINK_IN_CHILD(richards_lmh)
37 
38 
39 namespace it=Input::Type;
40 
41 
43 {
44  *this += water_content_saturated.name("water_content_saturated")
45  .description(R"(Saturated water content (($ \theta_s $)).
46  Relative volume of water in a reference volume of a saturated porous media.)")
47  .input_default("0.0")
48  .units( UnitSI::dimensionless() );
49 
50  *this += water_content_residual.name("water_content_residual")
51  .description(R"(Residual water content (($ \theta_r $)).
52  Relative volume of water in a reference volume of an ideally dry porous media.)")
53  .input_default("0.0")
54  .units( UnitSI::dimensionless() );
55 
56  *this += genuchten_p_head_scale.name("genuchten_p_head_scale")
57  .description(R"(The van Genuchten pressure head scaling parameter (($ \alpha $)).
58  It is related to the inverse of the air entry pressure, i.e. the pressure
59  where the relative water content starts to decrease below 1.)")
60  .input_default("0.0")
61  .units( UnitSI().m(-1) );
62 
63  *this += genuchten_n_exponent.name("genuchten_n_exponent")
64  .description("The van Genuchten exponent parameter (($ n $)).")
65  .input_default("2.0")
66  .units( UnitSI::dimensionless() );
67 }
68 
69 
71  it::Record field_descriptor = it::Record("RichardsLMH_Data",FieldCommon::field_descriptor_record_description("RichardsLMH_Data"))
73  .copy_keys( RichardsLMH::EqData().make_field_descriptor_type("RichardsLMH_Data_aux") )
74  .close();
75 
76  auto model_selection = it::Selection("Soil_Model_Type", "")
77  .add_value(SoilModelBase::van_genuchten, "van_genuchten", "Van Genuchten soil model with cutting near zero.")
78  .add_value(SoilModelBase::irmay, "irmay", "Irmay model for conductivity, Van Genuchten model for the water content. Suitable for bentonite.")
79  .close();
80 
81  auto soil_rec = it::Record("SoilModel", "Soil model settings.")
82  .allow_auto_conversion("model_type")
83  .declare_key("model_type", model_selection, it::Default("\"van_genuchten\""),
84  "Selection of the globally applied soil model. In future we replace this key by a field for selection of the model."
85  "That will allow usage of different soil model in a single simulation.")
86  .declare_key("cut_fraction", it::Double(0.0,1.0), it::Default("0.999"),
87  "Fraction of the water content where we cut and rescale the curve.")
88  .close();
89 
90  return it::Record("Flow_Richards_LMH", "Lumped Mixed-Hybrid solver for unsteady unsaturated Darcy flow.")
93  .declare_key("input_fields", it::Array( field_descriptor ), it::Default::obligatory(),
94  "Input data for Darcy flow model.")
95  .declare_key("soil_model", soil_rec, it::Default("\"van_genuchten\""),
96  "Soil model settings.")
97  .close();
98 }
99 
100 
101 const int RichardsLMH::registrar =
102  Input::register_class< RichardsLMH, Mesh &, const Input::Record >("Flow_Richards_LMH") +
104 
105 
106 
108  : DarcyMH(mesh_in, in_rec, tm)
109 {
110  data_ = make_shared<EqData>();
113  //data_->edge_new_local_4_mesh_idx_ = &(this->edge_new_local_4_mesh_idx_);
114 }
115 
117 
118  auto model_rec = input_record_.val<Input::Record>("soil_model");
119  auto model_type = model_rec.val<SoilModelBase::SoilModelType>("model_type");
120  double fraction= model_rec.val<double>("cut_fraction");
121  if (model_type == SoilModelBase::van_genuchten)
122  data_->soil_model_ = std::make_shared<SoilModel_VanGenuchten>(fraction);
123  else if (model_type == SoilModelBase::irmay)
124  data_->soil_model_ = std::make_shared<SoilModel_Irmay>(fraction);
125  else
126  ASSERT(false);
127 
128  // create edge vectors
129  unsigned int n_local_edges = mh_dh.edge_new_local_4_mesh_idx_.size();
130  unsigned int n_local_sides = mh_dh.side_ds->lsize();
131  data_->phead_edge_.resize( n_local_edges);
132  data_->water_content_previous_it.resize(n_local_sides);
133  data_->water_content_previous_time.resize(n_local_sides);
134  data_->capacity.resize(n_local_sides);
135 
136  Distribution ds_split_edges(n_local_edges, PETSC_COMM_WORLD);
137  vector<int> local_edge_rows(n_local_edges);
138 
139  IS is_loc;
140  for(auto item : mh_dh.edge_new_local_4_mesh_idx_) {
141  local_edge_rows[item.second]=mh_dh.row_4_edge[item.first];
142  }
143  ISCreateGeneral(PETSC_COMM_SELF, local_edge_rows.size(),
144  &(local_edge_rows[0]), PETSC_COPY_VALUES, &(is_loc));
145 
146  VecScatterCreate(schur0->get_solution(), is_loc,
147  data_->phead_edge_.petsc_vec(), PETSC_NULL, &solution_2_edge_scatter_);
148  chkerr(ISDestroy(&is_loc));
149 
150 }
151 
152 
154 {
155 
156 }
157 
158 
160 {
161  // apply initial condition
162  // cycle over local element rows
163  double init_value;
164 
165  for (unsigned int i_loc_el = 0; i_loc_el < mh_dh.el_ds->lsize(); i_loc_el++) {
166  auto ele_ac = mh_dh.accessor(i_loc_el);
167 
168  init_value = data_->init_pressure.value(ele_ac.centre(), ele_ac.element_accessor());
169 
170  for (unsigned int i=0; i<ele_ac.element_accessor()->n_sides(); i++) {
171  int edge_row = ele_ac.edge_row(i);
172  uint n_sides_of_edge = ele_ac.element_accessor().side(i)->edge().n_sides();
173  VecSetValue(schur0->get_solution(),edge_row, init_value/n_sides_of_edge, ADD_VALUES);
174  }
175  VecSetValue(schur0->get_solution(),ele_ac.ele_row(), init_value,ADD_VALUES);
176  }
177  VecAssemblyBegin(schur0->get_solution());
178  VecAssemblyEnd(schur0->get_solution());
179 
180  // set water_content
181  // pretty ugly since postprocess change fluxes, which cause bad balance, so we must set them back
182  VecCopy(schur0->get_solution(), previous_solution); // store solution vector
183  postprocess();
184  VecSwap(schur0->get_solution(), previous_solution); // restore solution vector
185 
186  //DebugOut() << "init sol:\n";
187  //VecView( schur0->get_solution(), PETSC_VIEWER_STDOUT_WORLD);
188  //DebugOut() << "init water content:\n";
189  //VecView( data_->water_content_previous_it.petsc_vec(), PETSC_VIEWER_STDOUT_WORLD);
190 
192 }
193 
194 
196 {
198  data_->water_content_previous_time.copy(data_->water_content_previous_it);
199  //VecCopy(schur0->get_solution(), previous_solution);
200 }
201 
202 bool RichardsLMH::zero_time_term(bool time_global) {
203  if (time_global) {
204  return (data_->storativity.input_list_size() == 0)
205  && (data_->water_content_saturated.input_list_size() == 0);
206 
207  } else {
208  return (data_->storativity.field_result(mesh_->region_db().get_region_set("BULK"))
209  == result_zeros)
210  && (data_->water_content_saturated.field_result(mesh_->region_db().get_region_set("BULK"))
211  == result_zeros);
212  }
213 }
214 
215 
217 {
218 
219  START_TIMER("RicharsLMH::assembly_linear_system");
220 
221  VecScatterBegin(solution_2_edge_scatter_, schur0->get_solution(), data_->phead_edge_.petsc_vec() , INSERT_VALUES, SCATTER_FORWARD);
222  VecScatterEnd(solution_2_edge_scatter_, schur0->get_solution(), data_->phead_edge_.petsc_vec() , INSERT_VALUES, SCATTER_FORWARD);
223 
224  data_->is_linear = data_->genuchten_p_head_scale.field_result(mesh_->region_db().get_region_set("BULK")) == result_zeros;
225 
226  bool is_steady = zero_time_term();
227  //DebugOut() << "Assembly linear system\n";
228  START_TIMER("full assembly");
229  if (typeid(*schur0) != typeid(LinSys_BDDC)) {
230  schur0->start_add_assembly(); // finish allocation and create matrix
231  }
232  data_->time_step_ = time_->dt();
233  auto multidim_assembler = AssemblyBase::create< AssemblyLMH >(data_);
234 
235 
238 
239  balance_->start_source_assembly(data_->water_balance_idx);
240  balance_->start_mass_assembly(data_->water_balance_idx);
241 
242  assembly_mh_matrix( multidim_assembler ); // fill matrix
243 
244  balance_->finish_source_assembly(data_->water_balance_idx);
245  balance_->finish_mass_assembly(data_->water_balance_idx);
246  //MatView( *const_cast<Mat*>(schur0->get_matrix()), PETSC_VIEWER_STDOUT_WORLD );
247  //VecView( *const_cast<Vec*>(schur0->get_rhs()), PETSC_VIEWER_STDOUT_WORLD);
248 
251 
252 
253  if (! is_steady) {
254  START_TIMER("fix time term");
255  //DebugOut() << "setup time term\n";
256  // assembly time term and rhs
258  }
259 }
260 
261 
262 
264 {
265  FEAL_ASSERT(false).error("Shold not be called.");
266 }
267 
268 
269 
270 
271 
273 
274  // update structures for balance of water volume
276 
277 
278 
279  int side_rows[4];
280  double values[4];
281 
282 
283  VecScatterBegin(solution_2_edge_scatter_, schur0->get_solution(), data_->phead_edge_.petsc_vec() , INSERT_VALUES, SCATTER_FORWARD);
284  VecScatterEnd(solution_2_edge_scatter_, schur0->get_solution(), data_->phead_edge_.petsc_vec() , INSERT_VALUES, SCATTER_FORWARD);
285 
286 
287  // modify side fluxes in parallel
288  // for every local edge take time term on digonal and add it to the corresponding flux
289  //PetscScalar *loc_prev_sol;
290  auto multidim_assembler = AssemblyBase::create< AssemblyLMH >(data_);
291 
292  //VecGetArray(previous_solution, &loc_prev_sol);
293  for (unsigned int i_loc = 0; i_loc < mh_dh.el_ds->lsize(); i_loc++) {
294  auto ele_ac = mh_dh.accessor(i_loc);
295  multidim_assembler[ele_ac.dim()-1]->update_water_content(ele_ac);
296 
297  double ele_scale = ele_ac.measure() *
298  data_->cross_section.value(ele_ac.centre(), ele_ac.element_accessor()) / ele_ac.n_sides();
299  double ele_source = data_->water_source_density.value(ele_ac.centre(), ele_ac.element_accessor())
300  + data_->extra_source.value(ele_ac.centre(), ele_ac.element_accessor());
301  //double storativity = data_->storativity.value(ele_ac.centre(), ele_ac.element_accessor());
302 
303  for (unsigned int i=0; i<ele_ac.element_accessor()->n_sides(); i++) {
304  //unsigned int loc_edge_row = ele_ac.edge_local_row(i);
305  side_rows[i] = ele_ac.side_row(i);
306  double water_content = data_->water_content_previous_it[ele_ac.side_local_idx(i)];
307  double water_content_previous_time = data_->water_content_previous_time[ele_ac.side_local_idx(i)];
308 
309  values[i] = ele_scale * ele_source - ele_scale * (water_content - water_content_previous_time) / time_->dt();
310  }
311  VecSetValues(schur0->get_solution(), ele_ac.n_sides(), side_rows, values, ADD_VALUES);
312  }
313 
314 
315  VecAssemblyBegin(schur0->get_solution());
316  //VecRestoreArray(previous_solution, &loc_prev_sol);
317  VecAssemblyEnd(schur0->get_solution());
318 
319 }
void initialize_specific() override
FieldSet * eq_data_
Definition: equation.hh:232
Distribution * side_ds
Output class for darcy_flow_mh model.
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:329
void assembly_mh_matrix(MultidimAssembly &assembler)
unsigned int uint
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:598
void postprocess() override
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
virtual void start_add_assembly()
Definition: linsys.hh:341
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:264
Abstract linear system class.
Definition: balance.hh:37
bool zero_time_term(bool time_global=false) override
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
Definitions of ASSERTS.
Definition: mesh.h:76
static const Input::Type::Record & type_field_descriptor()
Distribution * el_ds
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:148
virtual void finish_assembly()=0
const RegionDB & region_db() const
Definition: mesh.h:141
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:347
static const std::string field_descriptor_record_description(const string &record_name)
Definition: field_common.cc:73
Basic time management functionality for unsteady (and steady) solvers (class Equation).
bool solution_changed_for_scatter
std::shared_ptr< EqData > data_
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Basic time management class.
Class for declaration of inputs sequences.
Definition: type_base.hh:346
std::unordered_map< unsigned int, unsigned int > edge_new_local_4_mesh_idx_
void prepare_new_time_step() override
postprocess velocity field (add sources)
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
VecScatter solution_2_edge_scatter_
PETSC scatter from the solution vector to the parallel edge vector with ghost values.
Global macros to enhance readability and debugging, general constants.
virtual Record & allow_auto_conversion(const string &from_key)
Allows shorter input of the Record providing only value of the from_key given as the parameter...
Definition: type_record.cc:133
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:541
void setup_time_term() override
std::shared_ptr< Balance > balance_
Edge lumped mixed-hybrid solution of unsteady Darcy flow.
Definition: richards_lmh.hh:62
const Vec & get_solution()
Definition: linsys.hh:282
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
const Ret val(const string &key) const
#define START_TIMER(tag)
Starts a timer with specified tag.
Mesh * mesh_
Definition: equation.hh:223
Selection & add_value(const int value, const std::string &key, const std::string &description="", TypeBase::attribute_map attributes=TypeBase::attribute_map())
Adds one new value with name given by key to the Selection.
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
static Input::Type::Abstract & get_input_type()
LinSys * schur0
virtual PetscErrorCode rhs_zero_entries()
Definition: linsys.hh:273
std::shared_ptr< EqData > data_
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
MH_DofHandler mh_dh
LongIdx * row_4_edge
const Selection & close() const
Close the Selection, no more values can be added.
Input::Record input_record_
Definition: equation.hh:225
static const Input::Type::Record & get_input_type()
Definition: richards_lmh.cc:70
double dt() const
RichardsLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
LocalElementAccessorBase< 3 > accessor(uint local_ele_idx)
void read_initial_condition() override
void assembly_linear_system() override
void set_matrix_changed()
Definition: linsys.hh:212
Record type proxy class.
Definition: type_record.hh:182
static const int registrar
Registrar of class to factory.
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
Vec previous_solution
Template for classes storing finite set of named values.
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:180
static const Input::Type::Record & get_input_type()
TimeGovernor * time_
Definition: equation.hh:224
void assembly_source_term() override
Source term is implemented differently in LMH version.
#define FEAL_ASSERT(expr)
Definition of assert for debug and release mode.
Definition: asserts.hh:280
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
unsigned int lsize(int proc) const
get local size