Flow123d  master-f44eb46
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 
40 {
41  *this += water_content.name("water_content")
42  .units(UnitSI::dimensionless())
44  .description(R"(Water content.
45  It is a fraction of water volume to the whole volume.)");
46  *this += conductivity_richards.name("conductivity_richards")
47  .units( UnitSI().m().s(-1) )
49  .description("Computed isotropic scalar conductivity by the soil model.");
50 
51  *this += water_content_saturated.name("water_content_saturated")
52  .description(R"(Saturated water content (($ \theta_s $)).
53  Relative volume of water in a reference volume of a saturated porous media.)")
54  .input_default("0.0")
55  .units( UnitSI::dimensionless() );
56 
57  *this += water_content_residual.name("water_content_residual")
58  .description(R"(Residual water content (($ \theta_r $)).
59  Relative volume of water in a reference volume of an ideally dry porous media.)")
60  .input_default("0.0")
61  .units( UnitSI::dimensionless() );
62 
63  *this += genuchten_p_head_scale.name("genuchten_p_head_scale")
64  .description(R"(The van Genuchten pressure head scaling parameter (($ \alpha $)).
65  It is related to the inverse of the air entry pressure, i.e. the pressure
66  where the relative water content starts to decrease below 1.)")
67  .input_default("0.0")
68  .units( UnitSI().m(-1) );
69 
70  *this += genuchten_n_exponent.name("genuchten_n_exponent")
71  .description("The van Genuchten exponent parameter (($ n $)).")
72  .input_default("2.0")
73  .units( UnitSI::dimensionless() );
74 
75  this->set_default_fieldset();
76 }
77 
78 
80 : DarcyLMH::EqData::EqData() {}
81 
82 
84  namespace it=Input::Type;
85 
86  it::Record field_descriptor = it::Record(equation_name() + "_Data",
89  .copy_keys( RichardsLMH::EqFields().make_field_descriptor_type(equation_name() + "_Data_aux") )
90  .close();
91 
92  auto model_selection = it::Selection("Soil_Model_Type", "")
93  .add_value(SoilModelBase::van_genuchten, "van_genuchten", "Van Genuchten soil model with cutting near zero.")
94  .add_value(SoilModelBase::irmay, "irmay", "Irmay model for conductivity, Van Genuchten model for the water content. Suitable for bentonite.")
95  .close();
96 
97  auto soil_rec = it::Record("SoilModel", "Soil model settings.")
98  .allow_auto_conversion("model_type")
99  .declare_key("model_type", model_selection, it::Default("\"van_genuchten\""),
100  "Selection of the globally applied soil model. In future we replace this key by a field for selection of the model."
101  "That will allow usage of different soil model in a single simulation.")
102  .declare_key("cut_fraction", it::Double(0.0,1.0), it::Default("0.999"),
103  "Fraction of the water content where we cut and rescale the curve.")
104  .close();
105 
107 
108  return it::Record(equation_name(), "Lumped Mixed-Hybrid solver for unsteady unsaturated Darcy flow.")
111  .declare_key("input_fields", it::Array( field_descriptor ), it::Default::obligatory(),
112  "Input data for Darcy flow model.")
114  IT::Default("{ \"fields\": [ \"pressure_p0\", \"velocity_p0\" ] }"),
115  "Specification of output fields and output times.")
116  .declare_key("soil_model", soil_rec, it::Default("\"van_genuchten\""),
117  "Soil model settings.")
118  .close();
119 }
120 
121 
122 const int RichardsLMH::registrar =
123  Input::register_class< RichardsLMH, Mesh &, const Input::Record >(equation_name()) +
125 
126 
127 
129  : DarcyLMH(mesh_in, in_rec, tm),
131 {
132  eq_fields_ = make_shared<EqFields>();
133  eq_data_ = make_shared<EqData>();
136  this->eq_fieldset_ = eq_fields_;
137  //eq_data_->edge_new_local_4_mesh_idx_ = &(this->edge_new_local_4_mesh_idx_);
138 
139  eq_fields_->set_mesh(*mesh_);
140 }
141 
142 
144 
145  auto model_rec = input_record_.val<Input::Record>("soil_model");
146  auto model_type = model_rec.val<SoilModelBase::SoilModelType>("model_type");
147  double fraction= model_rec.val<double>("cut_fraction");
148  if (model_type == SoilModelBase::van_genuchten)
149  eq_data_->soil_model_ = std::make_shared<SoilModel_VanGenuchten>(fraction);
150  else if (model_type == SoilModelBase::irmay)
151  eq_data_->soil_model_ = std::make_shared<SoilModel_Irmay>(fraction);
152  else
153  ASSERT_PERMANENT(false);
154 
155  // create edge vectors
156  eq_data_->water_content_previous_time = eq_data_->dh_cr_disc_->create_vector();
157  eq_data_->capacity = eq_data_->dh_cr_disc_->create_vector();
158 
159  ASSERT_PTR(mesh_);
160  eq_data_->mesh = mesh_;
161 
162  eq_fields_->water_content_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_cr_disc_);
163  eq_fields_->water_content.set(eq_fields_->water_content_ptr, 0.0);
164 
165  eq_fields_->conductivity_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_p_);
166  eq_fields_->conductivity_richards.set(eq_fields_->conductivity_ptr, 0.0);
167 
168 
169  //eq_data_->multidim_assembler = AssemblyFlowBase::create< AssemblyRichards >(eq_fields_, eq_data_);
170 }
171 
172 
173 //void RichardsLMH::initial_condition_postprocess()
174 //{
175 // // modify side fluxes in parallel
176 // // for every local edge take time term on diagonal and add it to the corresponding flux
177 //
178 // for ( DHCellAccessor dh_cell : eq_data_->dh_->own_range() ) {
179 // eq_data_->multidim_assembler[dh_cell.elm().dim()-1]->update_water_content(dh_cell);
180 // }
181 //}
182 
183 
185 {
186  eq_data_->p_edge_solution_previous_time.copy_from(eq_data_->p_edge_solution);
187  VectorMPI water_content_vec = eq_fields_->water_content_ptr->vec();
188  eq_data_->water_content_previous_time.copy_from(water_content_vec);
189 
190  eq_data_->p_edge_solution_previous_time.local_to_ghost_begin();
191  eq_data_->p_edge_solution_previous_time.local_to_ghost_end();
192 }
193 
194 bool RichardsLMH::zero_time_term(bool time_global) {
195  if (time_global) {
196  return (eq_fields_->storativity.input_list_size() == 0)
197  && (eq_fields_->water_content_saturated.input_list_size() == 0);
198 
199  } else {
200  return (eq_fields_->storativity.field_result(mesh_->region_db().get_region_set("BULK"))
201  == result_zeros)
202  && (eq_fields_->water_content_saturated.field_result(mesh_->region_db().get_region_set("BULK"))
203  == result_zeros);
204  }
205 }
206 
207 
209 {
210  START_TIMER("RicharsLMH::assembly_linear_system");
211 
212  eq_data_->p_edge_solution.local_to_ghost_begin();
213  eq_data_->p_edge_solution.local_to_ghost_end();
214 
215  eq_data_->is_linear = eq_fields_->genuchten_p_head_scale.field_result(mesh_->region_db().get_region_set("BULK")) == result_zeros;
216 
217  //DebugOut() << "Assembly linear system\n";
218  START_TIMER("full assembly");
219 // if (typeid(*schur0) != typeid(LinSys_BDDC)) {
220 // schur0->start_add_assembly(); // finish allocation and create matrix
221 // schur_compl->start_add_assembly();
222 // }
223 
225 
226  eq_data_->time_step_ = time_->dt();
227 
230 
231 // assembly_mh_matrix( eq_data_->multidim_assembler ); // fill matrix
232  START_TIMER("RichardsLMH::assembly_steady_mh_matrix");
233  this->mh_matrix_assembly_->assemble(eq_data_->dh_); // fill matrix
234  END_TIMER("RichardsLMH::assembly_steady_mh_matrix");
235 
238 }
239 
240 
244  this->mh_matrix_assembly_ = new GenericAssembly< MHMatrixAssemblyRichards >(this->eq_fields_.get(), this->eq_data_.get());
246 }
247 
248 
252 }
253 
254 
256  if (init_cond_postprocess_assembly_!=nullptr) {
259  }
260 }
Definitions of ASSERTS.
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
static Input::Type::Abstract & get_input_type()
static const Input::Type::Instance & get_input_type(FieldSet &eq_data, const std::string &equation_name)
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
GenericAssemblyBase * mh_matrix_assembly_
EqFields & eq_fields()
std::shared_ptr< EqData > eq_data_
GenericAssembly< ReadInitCondAssemblyLMH > * read_init_cond_assembly_
general assembly objects, hold assembly objects of appropriate dimension
GenericAssemblyBase * reconstruct_schur_assembly_
static const Input::Type::Record & type_field_descriptor()
static const Input::Type::Record & get_input_type()
std::shared_ptr< EqFields > eq_fields_
LinSys & lin_sys_schur()
Getter for the linear system of the 2. Schur complement.
Input::Record input_record_
Definition: equation.hh:242
std::shared_ptr< FieldSet > eq_fieldset_
Definition: equation.hh:249
TimeGovernor * time_
Definition: equation.hh:241
Mesh * mesh_
Definition: equation.hh:240
static const std::string field_descriptor_record_description(const string &record_name)
Definition: field_common.cc:72
static constexpr Mask equation_result
Match result fields. These are never given by input or copy of input.
Definition: field_flag.hh:55
virtual void assemble(std::shared_ptr< DOFHandlerMultiDim > dh)=0
void assemble(std::shared_ptr< DOFHandlerMultiDim > dh) override
General assemble methods.
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
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
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 & 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
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
Template for classes storing finite set of named values.
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.
const Selection & close() const
Close the Selection, no more values can be added.
void set_matrix_changed()
Definition: linsys.hh:212
virtual void start_add_assembly()
Definition: linsys.hh:341
virtual void finish_assembly()=0
virtual PetscErrorCode rhs_zero_entries()
Definition: linsys.hh:273
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:264
const RegionDB & region_db() const
Definition: mesh.h:175
Definition: mesh.h:362
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:328
EqData()
Constructor.
Definition: richards_lmh.cc:79
Edge lumped mixed-hybrid solution of unsteady Darcy flow.
Definition: richards_lmh.hh:66
bool zero_time_term(bool time_global=false) override
static const int registrar
Registrar of class to factory.
void read_init_cond_asm() override
Call assemble of read_init_cond_assembly_ and init_cond_postprocess_assembly_.
virtual ~RichardsLMH() override
void accept_time_step() override
postprocess velocity field (add sources)
static const Input::Type::Record & get_input_type()
Definition: richards_lmh.cc:83
GenericAssembly< InitCondPostprocessAssembly > * init_cond_postprocess_assembly_
general assembly object, hold assembly objects of appropriate dimension
void initialize_specific() override
void assembly_linear_system() override
RichardsLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
std::shared_ptr< EqFields > eq_fields_
void initialize_asm() override
Create and initialize assembly objects.
std::shared_ptr< EqData > eq_data_
static std::string equation_name()
Basic time management functionality for unsteady (and steady) solvers (class Equation).
double dt() const
Class for representation SI units of Fields.
Definition: unit_si.hh:40
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
Output class for darcy_flow_mh model.
@ result_zeros
Global macros to enhance readability and debugging, general constants.
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.
Basic time management class.