Flow123d  JS_before_hm-1804-gf2ad740aa
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")
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")
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")
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")
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(data_->water_content_ptr, 0.0);
153 
154  data_->conductivity_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(data_->dh_p_);
155  data_->conductivity_richards.set(data_->conductivity_ptr, 0.0);
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 }
FieldCommon::units
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
Definition: field_common.hh:150
EquationBase::mesh_
Mesh * mesh_
Definition: equation.hh:218
RichardsLMH::accept_time_step
void accept_time_step() override
postprocess velocity field (add sources)
Definition: richards_lmh.cc:173
result_zeros
@ result_zeros
Definition: field_algo_base.hh:74
TimeGovernor::dt
double dt() const
Definition: time_governor.hh:558
UnitSI::dimensionless
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
assembly_richards.hh
factory.hh
vector_mpi.hh
time_governor.hh
Basic time management class.
RichardsLMH::EqData::water_content_saturated
Field< 3, FieldValue< 3 >::Scalar > water_content_saturated
Definition: richards_lmh.hh:72
RichardsLMH::EqData::water_content_residual
Field< 3, FieldValue< 3 >::Scalar > water_content_residual
Definition: richards_lmh.hh:73
ASSERT
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:347
EquationBase::time_
TimeGovernor * time_
Definition: equation.hh:219
RichardsLMH::EqData::conductivity_richards
Field< 3, FieldValue< 3 >::Scalar > conductivity_richards
Definition: richards_lmh.hh:81
include_fadbad.hh
RichardsLMH::EqData::EqData
EqData()
Definition: richards_lmh.cc:41
Input::Record::val
const Ret val(const string &key) const
Definition: accessors_impl.hh:31
asserts.hh
Definitions of ASSERTS.
richards_lmh.hh
Input::Type::Selection::close
const Selection & close() const
Close the Selection, no more values can be added.
Definition: type_selection.cc:65
Input::Type::Double
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
LinSys::rhs_zero_entries
virtual PetscErrorCode rhs_zero_entries()
Definition: linsys.hh:273
RichardsLMH::initial_condition_postprocess
void initial_condition_postprocess() override
Definition: richards_lmh.cc:162
FLOW123D_FORCE_LINK_IN_CHILD
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:157
RichardsLMH::zero_time_term
bool zero_time_term(bool time_global=false) override
Definition: richards_lmh.cc:183
LinSys::mat_zero_entries
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:264
FieldCommon::flags
FieldCommon & flags(FieldFlag::Flags::Mask mask)
Definition: field_common.hh:189
RichardsLMH::EqData::genuchten_p_head_scale
Field< 3, FieldValue< 3 >::Scalar > genuchten_p_head_scale
Definition: richards_lmh.hh:74
DarcyLMH::get_input_type
static const Input::Type::Record & get_input_type()
Definition: darcy_flow_lmh.cc:101
LinSys::set_matrix_changed
void set_matrix_changed()
Definition: linsys.hh:212
RichardsLMH::EqData
Definition: richards_lmh.hh:68
DarcyLMH::type_field_descriptor
static const Input::Type::Record & type_field_descriptor()
Definition: darcy_flow_lmh.cc:86
DarcyFlowMHOutput::get_input_type
static const Input::Type::Instance & get_input_type(FieldSet &eq_data, const std::string &equation_name)
Definition: darcy_flow_mh_output.cc:62
Input::Type::Record::size
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
Input::Type::Default
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
Input::Type::Record::derive_from
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
RichardsLMH::EqData::water_content
Field< 3, FieldValue< 3 >::Scalar > water_content
Definition: richards_lmh.hh:78
DarcyLMH::assembly_mh_matrix
void assembly_mh_matrix(MultidimAssembly &assembler)
Definition: darcy_flow_lmh.cc:681
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
FieldFlag::equation_result
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
sys_profiler.hh
darcy_flow_mh_output.hh
Output class for darcy_flow_mh model.
Input::Type::Record::allow_auto_conversion
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
schur.hh
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
TimeGovernor
Basic time management functionality for unsteady (and steady) solvers (class Equation).
Definition: time_governor.hh:310
RichardsLMH::data_
std::shared_ptr< EqData > data_
Definition: richards_lmh.hh:113
Input::Type::Default::obligatory
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
UnitSI
Class for representation SI units of Fields.
Definition: unit_si.hh:40
Mesh::region_db
const RegionDB & region_db() const
Definition: mesh.h:187
EquationBase::input_record_
Input::Record input_record_
Definition: equation.hh:220
RegionDB::get_region_set
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
FieldCommon::field_descriptor_record_description
static const std::string field_descriptor_record_description(const string &record_name)
Definition: field_common.cc:73
Input::Type::Record::declare_key
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
SoilModelBase::van_genuchten
@ van_genuchten
Definition: soil_models.hh:56
soil_models.hh
Input::Type::Selection
Template for classes storing finite set of named values.
Definition: type_selection.hh:65
DarcyFlowInterface::get_input_type
static Input::Type::Abstract & get_input_type()
Definition: darcy_flow_interface.hh:23
LinSys::start_add_assembly
virtual void start_add_assembly()
Definition: linsys.hh:341
Input::Type::Record::close
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Input::Type
Definition: balance.hh:41
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:182
DarcyLMH::lin_sys_schur
LinSys & lin_sys_schur()
Getter for the linear system of the 2. Schur complement.
Definition: darcy_flow_lmh.hh:284
FieldCommon::input_default
FieldCommon & input_default(const string &input_default)
Definition: field_common.hh:137
input_type.hh
Mesh
Definition: mesh.h:98
Input::Type::Record::copy_keys
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
Input::Type::Array
Class for declaration of inputs sequences.
Definition: type_base.hh:339
RichardsLMH::get_input_type
static const Input::Type::Record & get_input_type()
Definition: richards_lmh.cc:79
DHCellAccessor
Cell accessor allow iterate over DOF handler cells.
Definition: dh_cell_accessor.hh:43
LinSys::finish_assembly
virtual void finish_assembly()=0
global_defs.h
Global macros to enhance readability and debugging, general constants.
RichardsLMH::registrar
static const int registrar
Registrar of class to factory.
Definition: richards_lmh.hh:103
RichardsLMH::EqData::genuchten_n_exponent
Field< 3, FieldValue< 3 >::Scalar > genuchten_n_exponent
Definition: richards_lmh.hh:75
DarcyLMH
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
Definition: darcy_flow_lmh.hh:131
SoilModelBase::SoilModelType
SoilModelType
Definition: soil_models.hh:55
RichardsLMH::RichardsLMH
RichardsLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
Definition: richards_lmh.cc:121
EquationBase::eq_fieldset_
FieldSet * eq_fieldset_
Definition: equation.hh:227
balance.hh
VectorMPI
Definition: vector_mpi.hh:43
FieldCommon::description
FieldCommon & description(const string &description)
Definition: field_common.hh:125
RichardsLMH::assembly_linear_system
void assembly_linear_system() override
Definition: richards_lmh.cc:197
Input::Type::Selection::add_value
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.
Definition: type_selection.cc:50
ASSERT_PTR
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR)
Definition: asserts.hh:336
START_TIMER
#define START_TIMER(tag)
Starts a timer with specified tag.
Definition: sys_profiler.hh:115
DarcyLMH::data_
std::shared_ptr< EqData > data_
Definition: darcy_flow_lmh.hh:301
SoilModelBase::irmay
@ irmay
Definition: soil_models.hh:57
RichardsLMH::initialize_specific
void initialize_specific() override
Definition: richards_lmh.cc:131
FieldCommon::name
FieldCommon & name(const string &name)
Definition: field_common.hh:118