Flow123d  JS_before_hm-1013-g06f2edc
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 #include "coupling/balance.hh"
14 
15 #include "input/input_type.hh"
16 #include "input/factory.hh"
17 #include "flow/richards_lmh.hh"
18 #include "flow/soil_models.hh"
21 #include "tools/time_governor.hh"
22 
23 #include "petscmat.h"
24 #include "petscviewer.h"
25 #include "petscerror.h"
26 #include <armadillo>
27 
28 #include "la/schur.hh"
29 #include "la/vector_mpi.hh"
30 
31 
32 #include "tools/include_fadbad.hh" // for "fadbad.h", "badiff.h", "fadiff.h"
33 
34 
35 FLOW123D_FORCE_LINK_IN_CHILD(richards_lmh)
36 
37 
38 namespace it=Input::Type;
39 
40 
42 {
43  *this += water_content.name("water_content")
44  .units(UnitSI::dimensionless())
46  .description(R"(Water content.
47  It is a fraction of water volume to the whole volume.)");
48  *this += conductivity_richards.name("conductivity_richards")
49  .units( UnitSI().m().s(-1) )
51  .description("Computed isotropic scalar conductivity by the soil model.");
52 
53  *this += water_content_saturated.name("water_content_saturated")
54  .description(R"(Saturated water content (($ \theta_s $)).
55  Relative volume of water in a reference volume of a saturated porous media.)")
56  .input_default("0.0")
57  .units( UnitSI::dimensionless() );
58 
59  *this += water_content_residual.name("water_content_residual")
60  .description(R"(Residual water content (($ \theta_r $)).
61  Relative volume of water in a reference volume of an ideally dry porous media.)")
62  .input_default("0.0")
63  .units( UnitSI::dimensionless() );
64 
65  *this += genuchten_p_head_scale.name("genuchten_p_head_scale")
66  .description(R"(The van Genuchten pressure head scaling parameter (($ \alpha $)).
67  It is related to the inverse of the air entry pressure, i.e. the pressure
68  where the relative water content starts to decrease below 1.)")
69  .input_default("0.0")
70  .units( UnitSI().m(-1) );
71 
72  *this += genuchten_n_exponent.name("genuchten_n_exponent")
73  .description("The van Genuchten exponent parameter (($ n $)).")
74  .input_default("2.0")
75  .units( UnitSI::dimensionless() );
76 }
77 
78 
80  it::Record field_descriptor = it::Record("RichardsLMH_Data",FieldCommon::field_descriptor_record_description("RichardsLMH_Data"))
82  .copy_keys( RichardsLMH::EqData().make_field_descriptor_type("RichardsLMH_Data_aux") )
83  .close();
84 
85  auto model_selection = it::Selection("Soil_Model_Type", "")
86  .add_value(SoilModelBase::van_genuchten, "van_genuchten", "Van Genuchten soil model with cutting near zero.")
87  .add_value(SoilModelBase::irmay, "irmay", "Irmay model for conductivity, Van Genuchten model for the water content. Suitable for bentonite.")
88  .close();
89 
90  auto soil_rec = it::Record("SoilModel", "Soil model settings.")
91  .allow_auto_conversion("model_type")
92  .declare_key("model_type", model_selection, it::Default("\"van_genuchten\""),
93  "Selection of the globally applied soil model. In future we replace this key by a field for selection of the model."
94  "That will allow usage of different soil model in a single simulation.")
95  .declare_key("cut_fraction", it::Double(0.0,1.0), it::Default("0.999"),
96  "Fraction of the water content where we cut and rescale the curve.")
97  .close();
98 
99  RichardsLMH::EqData eq_data;
100 
101  return it::Record("Flow_Richards_LMH", "Lumped Mixed-Hybrid solver for unsteady unsaturated Darcy flow.")
104  .declare_key("input_fields", it::Array( field_descriptor ), it::Default::obligatory(),
105  "Input data for Darcy flow model.")
106  .declare_key("output", DarcyFlowMHOutput::get_input_type(eq_data, "Flow_Richards_LMH"),
107  IT::Default("{ \"fields\": [ \"pressure_p0\", \"velocity_p0\" ] }"),
108  "Specification of output fields and output times.")
109  .declare_key("soil_model", soil_rec, it::Default("\"van_genuchten\""),
110  "Soil model settings.")
111  .close();
112 }
113 
114 
115 const int RichardsLMH::registrar =
116  Input::register_class< RichardsLMH, Mesh &, const Input::Record >("Flow_Richards_LMH") +
118 
119 
120 
122  : DarcyLMH(mesh_in, in_rec, tm)
123 {
124  data_ = make_shared<EqData>();
127  //data_->edge_new_local_4_mesh_idx_ = &(this->edge_new_local_4_mesh_idx_);
128 }
129 
130 
132 
133  auto model_rec = input_record_.val<Input::Record>("soil_model");
134  auto model_type = model_rec.val<SoilModelBase::SoilModelType>("model_type");
135  double fraction= model_rec.val<double>("cut_fraction");
136  if (model_type == SoilModelBase::van_genuchten)
137  data_->soil_model_ = std::make_shared<SoilModel_VanGenuchten>(fraction);
138  else if (model_type == SoilModelBase::irmay)
139  data_->soil_model_ = std::make_shared<SoilModel_Irmay>(fraction);
140  else
141  ASSERT(false);
142 
143  // create edge vectors
144  data_->water_content_previous_time = data_->dh_cr_disc_->create_vector();
145  data_->capacity = data_->dh_cr_disc_->create_vector();
146 
147  ASSERT_PTR(mesh_);
148  data_->mesh = mesh_;
149  data_->set_mesh(*mesh_);
150 
151  data_->water_content_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(data_->dh_cr_disc_);
152  data_->water_content.set_field(mesh_->region_db().get_region_set("ALL"), data_->water_content_ptr);
153 
154  data_->conductivity_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(data_->dh_p_);
155  data_->conductivity_richards.set_field(mesh_->region_db().get_region_set("ALL"), data_->conductivity_ptr);
156 
157 
158  data_->multidim_assembler = AssemblyBase::create< AssemblyRichards >(data_);
159 }
160 
161 
163 {
164  // modify side fluxes in parallel
165  // for every local edge take time term on diagonal and add it to the corresponding flux
166 
167  for ( DHCellAccessor dh_cell : data_->dh_->own_range() ) {
168  data_->multidim_assembler[dh_cell.elm().dim()-1]->update_water_content(dh_cell);
169  }
170 }
171 
172 
174 {
175  data_->p_edge_solution_previous_time.copy_from(data_->p_edge_solution);
176  VectorMPI water_content_vec = data_->water_content_ptr->vec();
177  data_->water_content_previous_time.copy_from(water_content_vec);
178 
179  data_->p_edge_solution_previous_time.local_to_ghost_begin();
180  data_->p_edge_solution_previous_time.local_to_ghost_end();
181 }
182 
183 bool RichardsLMH::zero_time_term(bool time_global) {
184  if (time_global) {
185  return (data_->storativity.input_list_size() == 0)
186  && (data_->water_content_saturated.input_list_size() == 0);
187 
188  } else {
189  return (data_->storativity.field_result(mesh_->region_db().get_region_set("BULK"))
190  == result_zeros)
191  && (data_->water_content_saturated.field_result(mesh_->region_db().get_region_set("BULK"))
192  == result_zeros);
193  }
194 }
195 
196 
198 {
199  START_TIMER("RicharsLMH::assembly_linear_system");
200 
201  data_->p_edge_solution.local_to_ghost_begin();
202  data_->p_edge_solution.local_to_ghost_end();
203 
204  data_->is_linear = data_->genuchten_p_head_scale.field_result(mesh_->region_db().get_region_set("BULK")) == result_zeros;
205 
206  //DebugOut() << "Assembly linear system\n";
207  START_TIMER("full assembly");
208 // if (typeid(*schur0) != typeid(LinSys_BDDC)) {
209 // schur0->start_add_assembly(); // finish allocation and create matrix
210 // schur_compl->start_add_assembly();
211 // }
212 
214 
215  data_->time_step_ = time_->dt();
216 
219 
220  assembly_mh_matrix( data_->multidim_assembler ); // fill matrix
221 
224 }
void initialize_specific() override
FieldSet * eq_data_
Definition: equation.hh:227
Output class for darcy_flow_mh model.
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
void initial_condition_postprocess() override
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static const Input::Type::Record & get_input_type()
virtual void start_add_assembly()
Definition: linsys.hh:341
void assembly_mh_matrix(MultidimAssembly &assembler)
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:264
Abstract linear system class.
Definition: balance.hh:40
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:78
Cell accessor allow iterate over DOF handler cells.
virtual void finish_assembly()=0
static const Input::Type::Instance & get_input_type(FieldSet &eq_data, const std::string &equation_name)
const RegionDB & region_db() const
Definition: mesh.h:143
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:347
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
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).
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:339
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
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:534
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
Edge lumped mixed-hybrid solution of unsteady Darcy flow.
Definition: richards_lmh.hh:62
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
#define START_TIMER(tag)
Starts a timer with specified tag.
Mesh * mesh_
Definition: equation.hh:218
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()
std::shared_ptr< EqData > data_
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
LinSys & lin_sys_schur()
Getter for the linear system of the 2. Schur complement.
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR)
Definition: asserts.hh:336
const Selection & close() const
Close the Selection, no more values can be added.
Input::Record input_record_
Definition: equation.hh:220
static const Input::Type::Record & get_input_type()
Definition: richards_lmh.cc:79
double dt() const
RichardsLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
void assembly_linear_system() override
void set_matrix_changed()
Definition: linsys.hh:212
static const Input::Type::Record & type_field_descriptor()
Record type proxy class.
Definition: type_record.hh:182
void accept_time_step() override
postprocess velocity field (add sources)
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
Template for classes storing finite set of named values.
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:180
TimeGovernor * time_
Definition: equation.hh:219