Flow123d  master-f44eb46
darcy_flow_mh.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 darcy_flow_mh.hh
15  * @brief mixed-hybrid model of linear Darcy flow, possibly unsteady.
16  * @author Jan Brezina
17  *
18  * Main object for mixed-hybrid discretization of the linear elliptic PDE (Laplace)
19  * on a multidimensional domain. Discretization of saturated Darcy flow using
20  * RT0 approximation for the velocity
21  */
22 
23 /*
24  * list of files dependent on this one:
25  *
26  * posprocess.cc
27  * problem.cc
28  * main.hh
29  * transport.cc
30  */
31 
32 
33 #ifndef DARCY_FLOW_MH_HH
34 #define DARCY_FLOW_MH_HH
35 
36 #include <petscmat.h> // for Mat
37 #include <string.h> // for memcpy
38 #include <algorithm> // for swap
39 
40 #include <memory> // for shared_ptr, allocator...
41 #include <new> // for operator new[]
42 #include <string> // for string, operator<<
43 #include <vector> // for vector, vector<>::con...
44 #include <armadillo>
45 #include "fields/bc_field.hh" // for BCField
46 #include "fields/field.hh" // for Field
47 #include "fields/field_set.hh" // for FieldSet
48 #include "fields/field_values.hh" // for FieldValue<>::Scalar
49 #include "flow/darcy_flow_interface.hh" // for DarcyFlowInterface
50 #include "input/input_exception.hh" // for DECLARE_INPUT_EXCEPTION
51 #include "input/type_base.hh" // for Array
52 #include "input/type_generic.hh" // for Instance
53 #include "mesh/mesh.h" // for Mesh
54 #include "petscvec.h" // for Vec, _p_Vec, VecScatter
55 #include "system/exceptions.hh" // for ExcStream, operator<<
56 #include "tools/time_governor.hh" // for TimeGovernor
57 
58 class AssemblyFlowBase;
59 class Balance;
60 class DarcyFlowMHOutput;
61 class Element;
62 class Intersection;
63 class LinSys;
64 class LinSys_BDDC;
65 namespace Input {
66  class AbstractRecord;
67  class Record;
68  namespace Type {
69  class Record;
70  class Selection;
71  }
72 }
74 
75 /**
76  * @brief Mixed-hybrid model of linear Darcy flow, possibly unsteady.
77  *
78  * Abstract class for various implementations of Darcy flow. In future there should be
79  * one further level of abstraction for general time dependent problem.
80  *
81  * maybe TODO:
82  * split compute_one_step to :
83  * 1) prepare_next_timestep
84  * 2) actualize_solution - this is for iterative nonlinear solvers
85  *
86  */
87 
88 
89 /**
90  * @brief Mixed-hybrid of steady Darcy flow with sources and variable density.
91  *
92  * solve equations:
93  * @f[
94  * q= -{\mathbf{K}} \nabla h -{\mathbf{K}} R \nabla z
95  * @f]
96  * @f[
97  * \mathrm{div} q = f
98  * @f]
99  *
100  * where
101  * - @f$ q @f$ is flux @f$[ms^{-1}]@f$ for 3d, @f$[m^2s^{-1}]@f$ for 2d and @f$[m^3s^{-1}]@f$ for 1d.
102  * - @f$ \mathbf{K} @f$ is hydraulic tensor ( its orientation for 2d, 1d case is questionable )
103  * - @f$ h = \frac{\pi}{\rho_0 g}+z @f$ is pressure head, @f$ \pi, \rho_0, g @f$ are the pressure, water density, and acceleration of gravity , respectively.
104  * Assumes gravity force acts counter to the direction of the @f$ z @f$ axis.
105  * - @f$ R @f$ is destity or gravity variability coefficient. For density driven flow it should be
106  * @f[
107  * R = \frac{\rho}{\rho_0} -1 = \rho_0^{-1}\sum_{i=1}^s c_i
108  * @f]
109  * where @f$ c_i @f$ is concentration in @f$ kg m^{-3} @f$.
110  *
111  * The time key is optional, when not specified the equation is forced to steady regime. Using Steady TimeGovernor which have no dt constraints.
112  *
113  */
114 /**
115  * Model for transition coefficients due to Martin, Jaffre, Roberts (see manual for full reference)
116  *
117  * TODO:
118  * - how we can reuse field values computed during assembly
119  *
120  * TODO: Remove in future. It is supposed not to be improved anymore,
121  * however it is kept functioning aside of the LMH lumped version until
122  * the LMH version is stable and optimized.
123  */
124 
126 {
127 public:
128  TYPEDEF_ERR_INFO( EI_Reason, string);
129  DECLARE_EXCEPTION(ExcSolverDiverge,
130  << "Diverged nonlinear solver. Reason: " << EI_Reason::val
131  );
132  DECLARE_INPUT_EXCEPTION(ExcMissingTimeGovernor,
133  << "Missing the key 'time', obligatory for the transient problems.");
134 
135  /** Class with all fields used in the equation DarcyFlow.
136  * This is common to all implementations since this provides interface
137  * to this equation for possible coupling.
138  *
139  * This class is the base class for equation data also in DarcyLMH and RichardsLMH classes
140  * especially due to the common output class DarcyFlowMHOutput.
141  * This is the only dependence between DarcyMH and DarcyLMH classes.
142  */
143  class EqFields : public FieldSet {
144  public:
145 
146  /**
147  * For compatibility with old BCD file we have to assign integer codes starting from 1.
148  */
149  enum BC_Type {
150  none=0,
154  river=6
155  };
156 
157  /// Return a Selection corresponding to enum BC_Type.
159 
160  /// Creation of all fields.
161  EqFields();
162 
163  /// Return coords field
165  return this->X_;
166  }
167 
168 
174 
175  BCField<3, FieldValue<3>::Enum > bc_type; // Discrete need Selection for initialization
180 
183  Field<3, FieldValue<3>::Scalar > extra_storativity; /// Externally added storativity.
184  Field<3, FieldValue<3>::Scalar > extra_source; /// Externally added water source.
185 
191 
192  Field<3, FieldValue<3>::VectorFixed > gravity_field; /// Holds gravity vector acceptable in FieldModel
193  BCField<3, FieldValue<3>::VectorFixed > bc_gravity; /// Same as previous but used in boundary fields
197  };
198 
199  /**
200  * Class with all data used in the equation DarcyFlow.
201  * This is common to all implementations since this provides interface
202  * to this equation for possible coupling.
203  */
204  class EqData {
205  public:
206  /// Constructor.
207  EqData();
208 
209  /**
210  * Gravity vector and constant shift of pressure potential. Used to convert piezometric head
211  * to pressure head and vice versa.
212  */
213  arma::vec4 gravity_;
215 
216  // Mirroring the following members of DarcyMH:
218  std::shared_ptr<DOFHandlerMultiDim> dh_; ///< full DOF handler represents DOFs of sides, elements and edges
219  std::shared_ptr<SubDOFHandlerMultiDim> dh_cr_; ///< DOF handler represents DOFs of edges
220  std::shared_ptr<DOFHandlerMultiDim> dh_cr_disc_; ///< DOF handler represents DOFs of sides
221 
222 
224 
226 
227  std::shared_ptr<Balance> balance;
229 
230  unsigned int n_schur_compls;
231  int is_linear; ///< Hack fo BDDC solver.
232  bool force_no_neumann_bc; ///< auxiliary flag for switchting Dirichlet like BC
233 
234  /// Idicator of dirichlet or neumann type of switch boundary conditions.
236 
237  VectorMPI full_solution; //< full solution [vel,press,lambda] from 2. Schur complement
238 
240  };
241 
242  /// Selection for enum MortarMethod.
244 
245 
246 
247 
248 
249  DarcyMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm = nullptr);
250 
252  static const Input::Type::Record & get_input_type();
253 
254  void init_eq_data();
255  void initialize() override;
256  virtual void initialize_specific();
257  void zero_time_step() override;
258  void update_solution() override;
259  /// Solve the problem without moving to next time and without output.
260  void solve_time_step(bool output = true);
261 
262  /// postprocess velocity field (add sources)
263  virtual void postprocess();
264  virtual void output_data() override;
265 
266  virtual double solved_time() override;
267 
268  inline EqFields &eq_fields() { return *eq_fields_; }
269  inline EqData &eq_data() { return *eq_data_; }
270 
272  { eq_fields_->extra_storativity = extra_stor; }
273 
274  void set_extra_source(const Field<3, FieldValue<3>::Scalar> &extra_src)
275  { eq_fields_->extra_source = extra_src; }
276 
277  virtual ~DarcyMH() override;
278 
279 
280 protected:
281  /**
282  * Returns true is the fields involved in the time term have values that makes the time term zero.
283  * For time_global==true, it returns true if there are no field descriptors in the input list, so the
284  * fields )of the time ter) have their default values for whole simulation.
285  * If time_global==false (default), only the actual values are considered.
286  */
287  virtual bool zero_time_term(bool time_global=false);
288 
289  /// Solve method common to zero_time_step and update solution.
290  void solve_nonlinear();
291  void modify_system();
292  virtual void setup_time_term();
293  void prepare_new_time_step();
294 
295 
296  //void prepare_parallel();
297  //void make_row_numberings();
298 
299  /**
300  * Create and preallocate MH linear system (including matrix, rhs and solution vectors)
301  */
303 
304  /**
305  * Read initial condition into solution vector.
306  * Must be called after create_linear_system.
307  *
308  * For the LMH scheme we have to be able to save edge pressures in order to
309  * restart simulation or use results of one simulation as initial condition for other one.
310  */
311  virtual void read_initial_condition();
312 
314 
315  /**
316  * Part of per element assembly that is specific for MH and LMH respectively.
317  *
318  * This implemnets MH case:
319  * - compute conductivity scaling
320  * - assembly source term
321  * - no time term, managed by diagonal extraction etc.
322  */
323  //virtual void local_assembly_specific(AssemblyData &local_data);
324 
325  /**
326  * Allocates linear system matrix for MH.
327  * TODO:
328  * - use general preallocation methods in DofHandler
329  */
330  void allocate_mh_matrix();
331 
332  /**
333  * Assembles linear system matrix for MH.
334  * Element by element assembly is done using dim-template assembly class.
335  * Assembles only steady part of the equation.
336  * TODO:
337  * - include time term
338  * - add support for Robin type sources
339  * - support for nonlinear solvers - assembly either residual vector, matrix, or both (using FADBAD++)
340  */
341  void assembly_mh_matrix(MultidimAssembly& assembler);
342 
343  /// Source term is implemented differently in LMH version.
344  virtual void assembly_source_term();
345 
346  /**
347  * Assembly or update whole linear system.
348  */
349  virtual void assembly_linear_system();
350 
351  void set_mesh_data_for_bddc(LinSys_BDDC * bddc_ls);
352  /**
353  * Return a norm of residual vector.
354  * TODO: Introduce Equation::compute_residual() updating
355  * residual field, standard part of EqData.
356  */
357  virtual double solution_precision() const;
358 
359  /// Print darcy flow matrix in matlab format into a file.
360  void print_matlab_matrix(string matlab_file);
361 
362  /// Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
363  std::vector<int> get_component_indices_vec(unsigned int component) const;
364 
365  std::shared_ptr<Balance> balance_;
366 
368 
369  int size; // global size of MH matrix
370  int n_schur_compls; // number of shur complements to make
371 
372  // Propagate test for the time term to the assembly.
373  // This flag is necessary for switching BC to avoid setting zero neumann on the whole boundary in the steady case.
376 
377  // Setting of the nonlinear solver. TODO: Move to the solver class later on.
378  double tolerance_;
379  unsigned int min_n_it_;
380  unsigned int max_n_it_;
381  unsigned int nonlinear_iteration_; //< Actual number of completed nonlinear iterations, need to pass this information into assembly.
382 
383 
384  LinSys *schur0; //< whole MH Linear System
385 
390 
391  std::shared_ptr<EqFields> eq_fields_;
392  std::shared_ptr<EqData> eq_data_;
393 
394  friend class DarcyFlowMHOutput;
395  //friend class P0_CouplingAssembler;
396  //friend class P1_CouplingAssembler;
397 
398 private:
399  /// Registrar of class to factory
400  static const int registrar;
401 
402  static std::string equation_name()
403  { return "Flow_Darcy_MH";}
404 };
405 
406 
407 
408 void mat_count_off_proc_values(Mat m, Vec v);
409 /// Helper method fills range (min and max) of given component
410 void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component);
411 
412 
413 #endif //DARCY_FLOW_MH_HH
414 //-----------------------------------------------------------------------------
415 // vim: set cindent:
416 
MortarMethod
Type of experimental Mortar-like method for non-compatible 1d-2d interaction.
arma::vec3 gravity_vec_
std::shared_ptr< SubDOFHandlerMultiDim > dh_cr_
DOF handler represents DOFs of edges.
VectorMPI full_solution
MortarMethod mortar_method_
int is_linear
Hack fo BDDC solver.
MultidimAssembly multidim_assembler
EqData()
Constructor.
arma::vec4 gravity_
unsigned int n_schur_compls
std::shared_ptr< DOFHandlerMultiDim > dh_
full DOF handler represents DOFs of sides, elements and edges
std::vector< char > bc_switch_dirichlet
Idicator of dirichlet or neumann type of switch boundary conditions.
std::shared_ptr< DOFHandlerMultiDim > dh_cr_disc_
DOF handler represents DOFs of sides.
std::shared_ptr< Balance > balance
bool force_no_neumann_bc
auxiliary flag for switchting Dirichlet like BC
Field< 3, FieldValue< 3 >::Scalar > init_piezo_head
Same as previous but used in boundary fields.
BCField< 3, FieldValue< 3 >::Scalar > bc_piezo_head
BCField< 3, FieldValue< 3 >::Enum > bc_type
BCField< 3, FieldValue< 3 >::Scalar > bc_switch_piezo_head
BCField< 3, FieldValue< 3 >::VectorFixed > bc_gravity
Holds gravity vector acceptable in FieldModel.
Field< 3, FieldValue< 3 >::Scalar > extra_source
Externally added storativity.
Field< 3, FieldValue< 3 >::Scalar > field_ele_piezo_head
Field< 3, FieldValue< 3 >::TensorFixed > anisotropy
Field< 3, FieldValue< 3 >::Scalar > conductivity
Field< 3, FieldValue< 3 >::Scalar > field_ele_pressure
Externally added water source.
EqFields()
Creation of all fields.
Field< 3, FieldValue< 3 >::VectorFixed > gravity_field
Field< 3, FieldValue< 3 >::Scalar > field_edge_pressure
Field< 3, FieldValue< 3 >::Scalar > sigma
Field< 3, FieldValue< 3 >::Scalar > init_pressure
Field< 3, FieldValue< 3 >::Scalar > cross_section
Field< 3, FieldValue< 3 >::Scalar > water_source_density
FieldCoords & X()
Return coords field.
BCField< 3, FieldValue< 3 >::Scalar > bc_pressure
BCField< 3, FieldValue< 3 >::Scalar > bc_switch_pressure
BCField< 3, FieldValue< 3 >::Scalar > bc_robin_sigma
BCField< 3, FieldValue< 3 >::Scalar > bc_flux
static const Input::Type::Selection & get_bc_type_selection()
Return a Selection corresponding to enum BC_Type.
Field< 3, FieldValue< 3 >::VectorFixed > field_ele_velocity
Field< 3, FieldValue< 3 >::VectorFixed > flux
Field< 3, FieldValue< 3 >::Scalar > extra_storativity
Field< 3, FieldValue< 3 >::Scalar > storativity
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
DarcyFlowMHOutput * output_object
unsigned int max_n_it_
void set_mesh_data_for_bddc(LinSys_BDDC *bddc_ls)
Vec previous_solution
virtual double solved_time() override
void prepare_new_time_step()
int n_schur_compls
DECLARE_INPUT_EXCEPTION(ExcMissingTimeGovernor,<< "Missing the key 'time', obligatory for the transient problems.")
void assembly_mh_matrix(MultidimAssembly &assembler)
static const Input::Type::Record & type_field_descriptor()
virtual void read_initial_condition()
void print_matlab_matrix(string matlab_file)
Print darcy flow matrix in matlab format into a file.
virtual void initialize_specific()
std::shared_ptr< Balance > balance_
double tolerance_
EqData & eq_data()
std::vector< int > get_component_indices_vec(unsigned int component) const
Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
void allocate_mh_matrix()
Vec new_diagonal
void initialize() override
LinSys * schur0
static std::string equation_name()
void solve_nonlinear()
Solve method common to zero_time_step and update solution.
void init_eq_data()
virtual ~DarcyMH() override
void set_extra_source(const Field< 3, FieldValue< 3 >::Scalar > &extra_src)
std::shared_ptr< EqData > eq_data_
bool use_steady_assembly_
void zero_time_step() override
Vec steady_rhs
std::shared_ptr< EqFields > eq_fields_
static const int registrar
Registrar of class to factory.
unsigned int min_n_it_
virtual double solution_precision() const
unsigned int nonlinear_iteration_
void set_extra_storativity(const Field< 3, FieldValue< 3 >::Scalar > &extra_stor)
virtual void output_data() override
Write computed fields.
virtual void assembly_linear_system()
static const Input::Type::Selection & get_mh_mortar_selection()
Selection for enum MortarMethod.
void reconstruct_solution_from_schur(MultidimAssembly &assembler)
DECLARE_EXCEPTION(ExcSolverDiverge,<< "Diverged nonlinear solver. Reason: "<< EI_Reason::val)
virtual bool zero_time_term(bool time_global=false)
Vec steady_diagonal
void modify_system()
EqFields & eq_fields()
DarcyMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL.
virtual void postprocess()
postprocess velocity field (add sources)
TYPEDEF_ERR_INFO(EI_Reason, string)
static const Input::Type::Record & get_input_type()
void update_solution() override
void create_linear_system(Input::AbstractRecord rec)
bool data_changed_
void solve_time_step(bool output=true)
Solve the problem without moving to next time and without output.
virtual void setup_time_term()
virtual void assembly_source_term()
Source term is implemented differently in LMH version.
Mesh & mesh()
Definition: equation.hh:181
Container for various descendants of FieldCommonBase.
Definition: field_set.hh:159
FieldCoords X_
Field holds coordinates for computing of FieldFormulas.
Definition: field_set.hh:436
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Record type proxy class.
Definition: type_record.hh:182
Template for classes storing finite set of named values.
Definition: mesh.h:362
Basic time management functionality for unsteady (and steady) solvers (class Equation).
std::vector< std::shared_ptr< AssemblyFlowBase > > MultidimAssembly
void dofs_range(unsigned int n_dofs, unsigned int &min, unsigned int &max, unsigned int component)
Helper method fills range (min and max) of given component.
void mat_count_off_proc_values(Mat m, Vec v)
unsigned int uint
Abstract linear system class.
Definition: balance.hh:40
Basic time management class.