Flow123d  master-f44eb46
darcy_flow_lmh.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_lmh.hh
15  * @brief Lumped 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_LMH_HH
34 #define DARCY_FLOW_LMH_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 #include "la/vector_mpi.hh" // for VectorMPI
58 
59 
60 class Balance;
61 class DarcyFlowMHOutput;
62 class Element;
63 class Intersection;
64 class LinSys;
65 // class LinSys_BDDC;
66 class LocalSystem;
67 class LocalConstraint;
68 namespace Input {
69  class AbstractRecord;
70  class Record;
71  namespace Type {
72  class Record;
73  class Selection;
74  }
75 }
76 template<unsigned int dim> class ReadInitCondAssemblyLMH;
78 template< template<IntDim...> class DimAssembly> class GenericAssembly;
79 
80 
81 /**
82  * @brief Mixed-hybrid model of linear Darcy flow, possibly unsteady.
83  *
84  * Abstract class for various implementations of Darcy flow. In future there should be
85  * one further level of abstraction for general time dependent problem.
86  *
87  * maybe TODO:
88  * split compute_one_step to :
89  * 1) prepare_next_timestep
90  * 2) actualize_solution - this is for iterative nonlinear solvers
91  *
92  */
93 
94 
95 /**
96  * @brief Mixed-hybrid of steady Darcy flow with sources and variable density.
97  *
98  * solve equations:
99  * @f[
100  * q= -{\mathbf{K}} \nabla h -{\mathbf{K}} R \nabla z
101  * @f]
102  * @f[
103  * \mathrm{div} q = f
104  * @f]
105  *
106  * where
107  * - @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.
108  * - @f$ \mathbf{K} @f$ is hydraulic tensor ( its orientation for 2d, 1d case is questionable )
109  * - @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.
110  * Assumes gravity force acts counter to the direction of the @f$ z @f$ axis.
111  * - @f$ R @f$ is destity or gravity variability coefficient. For density driven flow it should be
112  * @f[
113  * R = \frac{\rho}{\rho_0} -1 = \rho_0^{-1}\sum_{i=1}^s c_i
114  * @f]
115  * where @f$ c_i @f$ is concentration in @f$ kg m^{-3} @f$.
116  *
117  * The time key is optional, when not specified the equation is forced to steady regime. Using Steady TimeGovernor which have no dt constraints.
118  *
119  *
120  * TODO:
121  * Make solution regular field (need FeSeystem and parallel DofHandler for edge pressures), then remove get_solution_vector from
122  * Equation interface.
123  */
124 /**
125  * Model for transition coefficients due to Martin, Jaffre, Roberts (see manual for full reference)
126  *
127  * TODO:
128  * - how we can reuse field values computed during assembly
129  *
130  */
131 
133 {
134 public:
135  TYPEDEF_ERR_INFO( EI_Reason, string);
136  DECLARE_EXCEPTION(ExcSolverDiverge,
137  << "Diverged nonlinear solver. Reason: " << EI_Reason::val
138  );
139  DECLARE_INPUT_EXCEPTION(ExcMissingTimeGovernor,
140  << "Missing the key 'time', obligatory for the transient problems.");
141 
142  /** Class with all fields used in the equation DarcyFlow.
143  * This is common to all implementations since this provides interface
144  * to this equation for possible coupling.
145  *
146  * This class uses the common output class DarcyFlowMHOutput.
147  * It is base class of RichardsLMH::EqFields.
148  * */
149  class EqFields : public FieldSet {
150  public:
151  /**
152  * For compatibility with old BCD file we have to assign integer codes starting from 1.
153  */
154  enum BC_Type {
155  none=0,
159  river=6
160  };
161 
162  /// Return a Selection corresponding to enum BC_Type.
164 
165  /// Creation of all fields.
166  EqFields();
167 
168  /// Return coords field
170  return this->X_;
171  }
172 
173 
179 
180  BCField<3, FieldValue<3>::Enum > bc_type; // Discrete need Selection for initialization
185 
188  Field<3, FieldValue<3>::Scalar > extra_storativity; /// Externally added storativity.
189  Field<3, FieldValue<3>::Scalar > extra_source; /// Externally added water source.
190 
196 
197  Field<3, FieldValue<3>::VectorFixed > gravity_field; /// Holds gravity vector acceptable in FieldModel
198  BCField<3, FieldValue<3>::VectorFixed > bc_gravity; /// Same as previous but used in boundary fields
202 
203  Field<3, FieldValue<3>::Scalar> ref_pressure; /// Precompute l2 difference outputs
206  };
207 
208  class EqData {
209  public:
210 
211  EqData();
212 
213  void init(); ///< Initialize vectors, ...
214  void reset(); ///< Reset data members
215 
216  /**
217  * Gravity vector and constant shift of pressure potential. Used to convert piezometric head
218  * to pressure head and vice versa.
219  */
220  arma::vec4 gravity_;
222 
223  // Mirroring the following members of DarcyLMH:
225  std::shared_ptr<DOFHandlerMultiDim> dh_; ///< full DOF handler represents DOFs of sides, elements and edges
226  std::shared_ptr<SubDOFHandlerMultiDim> dh_cr_; ///< DOF handler represents DOFs of edges
227  std::shared_ptr<DOFHandlerMultiDim> dh_cr_disc_; ///< DOF handler represents DOFs of sides
228  std::shared_ptr<SubDOFHandlerMultiDim> dh_p_; ///< DOF handler represents DOFs of element pressure
229 
230 
232 
233  // TODO: check usage of MortarMethod in LMH
235 
236  int is_linear; ///< Hack fo BDDC solver.
237  bool force_no_neumann_bc; ///< auxiliary flag for switchting Dirichlet like BC
238 
239  /// Idicator of dirichlet or neumann type of switch boundary conditions.
241 
242  VectorMPI full_solution; //< full solution [vel,press,lambda] from 2. Schur complement
243 
244  // Propagate test for the time term to the assembly.
245  // This flag is necessary for switching BC to avoid setting zero neumann on the whole boundary in the steady case.
247 
248  // for time term assembly
249  double time_step_;
250 
251  std::shared_ptr<LinSys> lin_sys_schur; //< Linear system of the 2. Schur complement.
252  VectorMPI p_edge_solution; //< 2. Schur complement solution
253  VectorMPI p_edge_solution_previous; //< 2. Schur complement previous solution (iterative)
254  VectorMPI p_edge_solution_previous_time; //< 2. Schur complement previous solution (time)
255 
257 
258  /// Shared Balance object
259  std::shared_ptr<Balance> balance_;
260 
261  unsigned int nonlinear_iteration_; //< Actual number of completed nonlinear iterations, need to pass this information into assembly.
262 
263  /// Following data members are stored in vectors, one item for every cell.
264 // /** TODO: Investigate why the hell do we need this flag.
265 // * If removed, it does not break any of the integration tests,
266 // * however it must influence the Dirichlet rows in matrix.
267 // */
268 // std::vector<unsigned int> dirichlet_edge;
269 
273  std::array<std::vector<unsigned int>, 3> loc_side_dofs;
274  std::array<std::vector<unsigned int>, 3> loc_edge_dofs;
275  std::array<unsigned int, 3> loc_ele_dof;
276 
277 // // std::shared_ptr<MortarAssemblyBase> mortar_assembly;
278 
279  std::vector<bool> save_local_system_; ///< Flag for saving the local system. Currently used only in case of seepage BC.
280  std::vector<bool> bc_fluxes_reconstruted; ///< Flag indicating whether the fluxes for seepage BC has been reconstructed already.
281  std::array<unsigned int, 3> schur_offset_; ///< Index offset in the local system for the Schur complement (of dim = 1,2,3).
282  };
283 
284  /// Selection for enum MortarMethod.
286 
287 
288 
289 
290 
291  DarcyLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm = nullptr);
292 
294  static const Input::Type::Record & get_input_type();
295 
296  void init_eq_data();
297  void initialize() override;
298  virtual void initialize_specific();
299  void zero_time_step() override;
300  void update_solution() override;
301 
302  /// Solve the problem without moving to next time and without output.
303  void solve_time_step(bool output = true);
304 
305  /// postprocess velocity field (add sources)
306  virtual void accept_time_step();
307  virtual void postprocess();
308  virtual void output_data() override;
309 
310  virtual double solved_time() override;
311 
312  inline EqFields &eq_fields() { return *eq_fields_; }
313  inline EqData &eq_data() { return *eq_data_; }
314 
315  /// Sets external storarivity field (coupling with other equation).
317  { eq_fields_->extra_storativity = extra_stor; }
318 
319  /// Sets external source field (coupling with other equation).
320  void set_extra_source(const Field<3, FieldValue<3>::Scalar> &extra_src)
321  { eq_fields_->extra_source = extra_src; }
322 
323  virtual ~DarcyLMH() override;
324 
325 
326 protected:
327  /**
328  * Returns true is the fields involved in the time term have values that makes the time term zero.
329  * For time_global==true, it returns true if there are no field descriptors in the input list, so the
330  * fields )of the time ter) have their default values for whole simulation.
331  * If time_global==false (default), only the actual values are considered.
332  */
333  virtual bool zero_time_term(bool time_global=false);
334 
335  /// Solve method common to zero_time_step and update solution.
336  void solve_nonlinear();
337 
338  /**
339  * Create and preallocate MH linear system (including matrix, rhs and solution vectors)
340  */
342 
343  /**
344  * Read initial condition into solution vector.
345  * Must be called after create_linear_system.
346  *
347  * For the LMH scheme we have to be able to save edge pressures in order to
348  * restart simulation or use results of one simulation as initial condition for other one.
349  */
350 // void read_initial_condition();
351 
352  /**
353  * In some circumstances, the intial condition must be processed.
354  * It is called at the end of @p read_initial_condition().
355  * This is used in Richards equation due the update of water content.
356  */
357 // virtual void initial_condition_postprocess();
358 
359  /**
360  * Allocates linear system matrix for MH.
361  * TODO:
362  * - use general preallocation methods in DofHandler
363  */
364  void allocate_mh_matrix();
365 
366  /**
367  * Assembles linear system matrix for MH.
368  * Element by element assembly is done using dim-template assembly class.
369  * Assembles only steady part of the equation.
370  * TODO:
371  * - include time term - DONE
372  * - add support for Robin type sources
373  * - support for nonlinear solvers - assembly either residual vector, matrix, or both (using FADBAD++)
374  */
375 // void assembly_mh_matrix(MultidimAssembly& assembler);
376 
377 // void reconstruct_solution_from_schur(MultidimAssembly& assembler);
378 
379  /**
380  * Assembly or update whole linear system.
381  */
382  virtual void assembly_linear_system();
383 
384 // void set_mesh_data_for_bddc(LinSys_BDDC * bddc_ls);
385  /**
386  * Return a norm of residual vector.
387  * TODO: Introduce Equation::compute_residual() updating
388  * residual field, standard part of EqData.
389  */
390 // virtual double solution_precision() const;
391 
392  /// Print darcy flow matrix in matlab format into a file.
393  void print_matlab_matrix(string matlab_file);
394 
395  /// Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
396  std::vector<int> get_component_indices_vec(unsigned int component) const;
397 
398  /// Getter for the linear system of the 2. Schur complement.
400  { return *(eq_data_->lin_sys_schur); }
401 
402  /// Create and initialize assembly objects
403  virtual void initialize_asm();
404 
405  /// Call assemble of read_init_cond_assembly_
406  virtual void read_init_cond_asm();
407 
408  std::shared_ptr<Balance> balance_;
409 
411 
412  int size; // global size of MH matrix
413 
415 
416  // Setting of the nonlinear solver. TODO: Move to the solver class later on.
417  double tolerance_;
418  unsigned int min_n_it_;
419  unsigned int max_n_it_;
420 
421  std::shared_ptr<EqFields> eq_fields_;
422  std::shared_ptr<EqData> eq_data_;
423 
424 
425  friend class DarcyFlowMHOutput;
426  //friend class P0_CouplingAssembler;
427  //friend class P1_CouplingAssembler;
428 
429  /// general assembly objects, hold assembly objects of appropriate dimension
433 private:
434  /// Registrar of class to factory
435  static const int registrar;
436 
437  static std::string equation_name()
438  { return "Flow_Darcy_LMH";}
439 };
440 
441 #endif //DARCY_FLOW_LMH_HH
442 //-----------------------------------------------------------------------------
443 // vim: set cindent:
444 
MortarMethod
Type of experimental Mortar-like method for non-compatible 1d-2d interaction.
std::array< unsigned int, 3 > schur_offset_
Index offset in the local system for the Schur complement (of dim = 1,2,3).
unsigned int nonlinear_iteration_
arma::vec3 gravity_vec_
std::vector< bool > bc_fluxes_reconstruted
Flag indicating whether the fluxes for seepage BC has been reconstructed already.
std::map< LongIdx, LocalSystem > seepage_bc_systems
int is_linear
Hack fo BDDC solver.
std::shared_ptr< SubDOFHandlerMultiDim > dh_p_
DOF handler represents DOFs of element pressure.
VectorMPI full_solution
std::vector< bool > save_local_system_
Flag for saving the local system. Currently used only in case of seepage BC.
void reset()
Reset data members.
std::vector< arma::vec > postprocess_solution_
MortarMethod mortar_method_
std::shared_ptr< SubDOFHandlerMultiDim > dh_cr_
DOF handler represents DOFs of edges.
bool force_no_neumann_bc
auxiliary flag for switchting Dirichlet like BC
VectorMPI p_edge_solution_previous_time
std::shared_ptr< DOFHandlerMultiDim > dh_cr_disc_
DOF handler represents DOFs of sides.
VectorMPI p_edge_solution
VectorMPI p_edge_solution_previous
std::shared_ptr< DOFHandlerMultiDim > dh_
full DOF handler represents DOFs of sides, elements and edges
std::vector< LocalConstraint > loc_constraint_
std::vector< LocalSystem > loc_system_
Following data members are stored in vectors, one item for every cell.
std::vector< char > bc_switch_dirichlet
Idicator of dirichlet or neumann type of switch boundary conditions.
std::shared_ptr< LinSys > lin_sys_schur
std::array< unsigned int, 3 > loc_ele_dof
std::array< std::vector< unsigned int >, 3 > loc_edge_dofs
void init()
Initialize vectors, ...
std::shared_ptr< Balance > balance_
Shared Balance object.
std::array< std::vector< unsigned int >, 3 > loc_side_dofs
Field< 3, FieldValue< 3 >::Scalar > water_source_density
Field< 3, FieldValue< 3 >::Scalar > ref_divergence
Field< 3, FieldValue< 3 >::Scalar > extra_storativity
Field< 3, FieldValue< 3 >::Scalar > field_ele_pressure
Externally added water source.
BCField< 3, FieldValue< 3 >::Enum > bc_type
Field< 3, FieldValue< 3 >::VectorFixed > gravity_field
Field< 3, FieldValue< 3 >::Scalar > sigma
Field< 3, FieldValue< 3 >::Scalar > init_pressure
Field< 3, FieldValue< 3 >::Scalar > storativity
BCField< 3, FieldValue< 3 >::Scalar > bc_pressure
BCField< 3, FieldValue< 3 >::Scalar > bc_flux
Field< 3, FieldValue< 3 >::TensorFixed > anisotropy
Field< 3, FieldValue< 3 >::Scalar > cross_section
FieldCoords & X()
Return coords field.
Field< 3, FieldValue< 3 >::Scalar > conductivity
Field< 3, FieldValue< 3 >::VectorFixed > field_ele_velocity
Field< 3, FieldValue< 3 >::Scalar > field_edge_pressure
BCField< 3, FieldValue< 3 >::Scalar > bc_switch_piezo_head
Field< 3, FieldValue< 3 >::Scalar > field_ele_piezo_head
Field< 3, FieldValue< 3 >::Scalar > init_piezo_head
Same as previous but used in boundary fields.
BCField< 3, FieldValue< 3 >::Scalar > bc_robin_sigma
static const Input::Type::Selection & get_bc_type_selection()
Return a Selection corresponding to enum BC_Type.
Field< 3, FieldValue< 3 >::Scalar > extra_source
Externally added storativity.
Field< 3, FieldValue< 3 >::VectorFixed > ref_velocity
Precompute l2 difference outputs.
Field< 3, FieldValue< 3 >::VectorFixed > flux
EqFields()
Creation of all fields.
BCField< 3, FieldValue< 3 >::Scalar > bc_switch_pressure
BCField< 3, FieldValue< 3 >::VectorFixed > bc_gravity
Holds gravity vector acceptable in FieldModel.
BCField< 3, FieldValue< 3 >::Scalar > bc_piezo_head
Field< 3, FieldValue< 3 >::Scalar > ref_pressure
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
GenericAssemblyBase * mh_matrix_assembly_
void initialize() override
void print_matlab_matrix(string matlab_file)
Print darcy flow matrix in matlab format into a file.
virtual double solved_time() override
void zero_time_step() override
virtual void postprocess()
EqFields & eq_fields()
void solve_nonlinear()
Solve method common to zero_time_step and update solution.
static std::string equation_name()
void set_extra_source(const Field< 3, FieldValue< 3 >::Scalar > &extra_src)
Sets external source field (coupling with other equation).
virtual void initialize_specific()
static const int registrar
Registrar of class to factory.
DarcyFlowMHOutput * output_object
std::shared_ptr< EqData > eq_data_
unsigned int max_n_it_
DECLARE_INPUT_EXCEPTION(ExcMissingTimeGovernor,<< "Missing the key 'time', obligatory for the transient problems.")
GenericAssembly< ReadInitCondAssemblyLMH > * read_init_cond_assembly_
general assembly objects, hold assembly objects of appropriate dimension
bool data_changed_
GenericAssemblyBase * reconstruct_schur_assembly_
double tolerance_
virtual void initialize_asm()
Create and initialize assembly objects.
void set_extra_storativity(const Field< 3, FieldValue< 3 >::Scalar > &extra_stor)
Sets external storarivity field (coupling with other equation).
void update_solution() override
void allocate_mh_matrix()
unsigned int min_n_it_
void solve_time_step(bool output=true)
Solve the problem without moving to next time and without output.
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)
static const Input::Type::Record & type_field_descriptor()
EqData & eq_data()
void create_linear_system(Input::AbstractRecord rec)
static const Input::Type::Record & get_input_type()
DECLARE_EXCEPTION(ExcSolverDiverge,<< "Diverged nonlinear solver. Reason: "<< EI_Reason::val)
DarcyLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL.
virtual bool zero_time_term(bool time_global=false)
std::shared_ptr< EqFields > eq_fields_
static const Input::Type::Selection & get_mh_mortar_selection()
Selection for enum MortarMethod.
std::shared_ptr< Balance > balance_
void init_eq_data()
virtual void accept_time_step()
postprocess velocity field (add sources)
virtual ~DarcyLMH() override
virtual void output_data() override
Write computed fields.
LinSys & lin_sys_schur()
Getter for the linear system of the 2. Schur complement.
TYPEDEF_ERR_INFO(EI_Reason, string)
virtual void assembly_linear_system()
virtual void read_init_cond_asm()
Call assemble of read_init_cond_assembly_.
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
Generic class of assemblation.
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).
unsigned int uint
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Abstract linear system class.
Definition: balance.hh:40
Basic time management class.