Flow123d  JS_before_hm-2206-gca2a8366f
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 #include "flow/darcy_flow_mh.hh" // for DarcyMH::EqData
60 
61 class Balance;
62 class DarcyFlowMHOutput;
63 class Element;
64 class Intersection;
65 class LinSys;
66 // class LinSys_BDDC;
67 class LocalSystem;
68 namespace Input {
69  class AbstractRecord;
70  class Record;
71  namespace Type {
72  class Record;
73  class Selection;
74  }
75 }
76 
77 /**
78  * @brief Mixed-hybrid model of linear Darcy flow, possibly unsteady.
79  *
80  * Abstract class for various implementations of Darcy flow. In future there should be
81  * one further level of abstraction for general time dependent problem.
82  *
83  * maybe TODO:
84  * split compute_one_step to :
85  * 1) prepare_next_timestep
86  * 2) actualize_solution - this is for iterative nonlinear solvers
87  *
88  */
89 
90 
91 /**
92  * @brief Mixed-hybrid of steady Darcy flow with sources and variable density.
93  *
94  * solve equations:
95  * @f[
96  * q= -{\mathbf{K}} \nabla h -{\mathbf{K}} R \nabla z
97  * @f]
98  * @f[
99  * \mathrm{div} q = f
100  * @f]
101  *
102  * where
103  * - @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.
104  * - @f$ \mathbf{K} @f$ is hydraulic tensor ( its orientation for 2d, 1d case is questionable )
105  * - @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.
106  * Assumes gravity force acts counter to the direction of the @f$ z @f$ axis.
107  * - @f$ R @f$ is destity or gravity variability coefficient. For density driven flow it should be
108  * @f[
109  * R = \frac{\rho}{\rho_0} -1 = \rho_0^{-1}\sum_{i=1}^s c_i
110  * @f]
111  * where @f$ c_i @f$ is concentration in @f$ kg m^{-3} @f$.
112  *
113  * The time key is optional, when not specified the equation is forced to steady regime. Using Steady TimeGovernor which have no dt constraints.
114  *
115  *
116  * TODO:
117  * Make solution regular field (need FeSeystem and parallel DofHandler for edge pressures), then remove get_solution_vector from
118  * Equation interface.
119  */
120 /**
121  * Model for transition coefficients due to Martin, Jaffre, Roberts (see manual for full reference)
122  *
123  * TODO:
124  * - how we can reuse field values computed during assembly
125  *
126  */
127 
129 {
130 public:
131  TYPEDEF_ERR_INFO( EI_Reason, string);
132  DECLARE_EXCEPTION(ExcSolverDiverge,
133  << "Diverged nonlinear solver. Reason: " << EI_Reason::val
134  );
135  DECLARE_INPUT_EXCEPTION(ExcMissingTimeGovernor,
136  << "Missing the key 'time', obligatory for the transient problems.");
137 
138  /** Class with all fields used in the equation DarcyFlow.
139  * This is common to all implementations since this provides interface
140  * to this equation for possible coupling.
141  *
142  * This class is derived from DarcyMH::EqData especially due to the common output class DarcyFlowMHOutput.
143  * This is the only dependence between DarcyMH and DarcyLMH classes.
144  * It is also base class of RichardsLMH::EqData.
145  * */
146  class EqFields : public DarcyMH::EqFields {
147  public:
148 
149  EqFields();
150  };
151 
152  class EqData : public DarcyMH::EqData {
153  public:
154 
155  EqData();
156 
157  std::shared_ptr<SubDOFHandlerMultiDim> dh_p_; ///< DOF handler represents DOFs of element pressure
158 
159  // Propagate test for the time term to the assembly.
160  // This flag is necessary for switching BC to avoid setting zero neumann on the whole boundary in the steady case.
162 
163  // for time term assembly
164  double time_step_;
165 
166  std::shared_ptr<LinSys> lin_sys_schur; //< Linear system of the 2. Schur complement.
167  VectorMPI p_edge_solution; //< 2. Schur complement solution
168  VectorMPI p_edge_solution_previous; //< 2. Schur complement previous solution (iterative)
169  VectorMPI p_edge_solution_previous_time; //< 2. Schur complement previous solution (time)
170 
172  };
173 
174  /// Selection for enum MortarMethod.
176 
177 
178 
179 
180 
181  DarcyLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm = nullptr);
182 
184  static const Input::Type::Record & get_input_type();
185 
186  void init_eq_data();
187  void initialize() override;
188  virtual void initialize_specific();
189  void zero_time_step() override;
190  void update_solution() override;
191 
192  /// Solve the problem without moving to next time and without output.
193  void solve_time_step(bool output = true);
194 
195  /// postprocess velocity field (add sources)
196  virtual void accept_time_step();
197  virtual void postprocess();
198  virtual void output_data() override;
199 
200  virtual double solved_time() override;
201 
202  inline EqFields &eq_fields() { return *eq_fields_; }
203  inline EqData &eq_data() { return *eq_data_; }
204 
205  /// Sets external storarivity field (coupling with other equation).
207  { eq_fields_->extra_storativity = extra_stor; }
208 
209  /// Sets external source field (coupling with other equation).
210  void set_extra_source(const Field<3, FieldValue<3>::Scalar> &extra_src)
211  { eq_fields_->extra_source = extra_src; }
212 
213  virtual ~DarcyLMH() override;
214 
215 
216 protected:
217  /**
218  * Returns true is the fields involved in the time term have values that makes the time term zero.
219  * For time_global==true, it returns true if there are no field descriptors in the input list, so the
220  * fields )of the time ter) have their default values for whole simulation.
221  * If time_global==false (default), only the actual values are considered.
222  */
223  virtual bool zero_time_term(bool time_global=false);
224 
225  /// Solve method common to zero_time_step and update solution.
226  void solve_nonlinear();
227 
228  /**
229  * Create and preallocate MH linear system (including matrix, rhs and solution vectors)
230  */
232 
233  /**
234  * Read initial condition into solution vector.
235  * Must be called after create_linear_system.
236  *
237  * For the LMH scheme we have to be able to save edge pressures in order to
238  * restart simulation or use results of one simulation as initial condition for other one.
239  */
240  void read_initial_condition();
241 
242  /**
243  * In some circumstances, the intial condition must be processed.
244  * It is called at the end of @p read_initial_condition().
245  * This is used in Richards equation due the update of water content.
246  */
247  virtual void initial_condition_postprocess();
248 
249  /**
250  * Allocates linear system matrix for MH.
251  * TODO:
252  * - use general preallocation methods in DofHandler
253  */
254  void allocate_mh_matrix();
255 
256  /**
257  * Assembles linear system matrix for MH.
258  * Element by element assembly is done using dim-template assembly class.
259  * Assembles only steady part of the equation.
260  * TODO:
261  * - include time term - DONE
262  * - add support for Robin type sources
263  * - support for nonlinear solvers - assembly either residual vector, matrix, or both (using FADBAD++)
264  */
265  void assembly_mh_matrix(MultidimAssembly& assembler);
266 
268 
269  /**
270  * Assembly or update whole linear system.
271  */
272  virtual void assembly_linear_system();
273 
274 // void set_mesh_data_for_bddc(LinSys_BDDC * bddc_ls);
275  /**
276  * Return a norm of residual vector.
277  * TODO: Introduce Equation::compute_residual() updating
278  * residual field, standard part of EqData.
279  */
280  virtual double solution_precision() const;
281 
282  /// Print darcy flow matrix in matlab format into a file.
283  void print_matlab_matrix(string matlab_file);
284 
285  /// Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
286  std::vector<int> get_component_indices_vec(unsigned int component) const;
287 
288  /// Getter for the linear system of the 2. Schur complement.
290  { return *(eq_data_->lin_sys_schur); }
291 
292  std::shared_ptr<Balance> balance_;
293 
295 
296  int size; // global size of MH matrix
297 
299 
300  // Setting of the nonlinear solver. TODO: Move to the solver class later on.
301  double tolerance_;
302  unsigned int min_n_it_;
303  unsigned int max_n_it_;
304  unsigned int nonlinear_iteration_; //< Actual number of completed nonlinear iterations, need to pass this information into assembly.
305 
306  std::shared_ptr<EqFields> eq_fields_;
307  std::shared_ptr<EqData> eq_data_;
308 
309  friend class DarcyFlowMHOutput;
310  //friend class P0_CouplingAssembler;
311  //friend class P1_CouplingAssembler;
312 
313 private:
314  /// Registrar of class to factory
315  static const int registrar;
316 };
317 
318 #endif //DARCY_FLOW_LMH_HH
319 //-----------------------------------------------------------------------------
320 // vim: set cindent:
321 
bc_field.hh
DarcyLMH::print_matlab_matrix
void print_matlab_matrix(string matlab_file)
Print darcy flow matrix in matlab format into a file.
Definition: darcy_flow_lmh.cc:1062
DarcyLMH::zero_time_term
virtual bool zero_time_term(bool time_global=false)
Definition: darcy_flow_lmh.cc:547
DarcyLMH::reconstruct_solution_from_schur
void reconstruct_solution_from_schur(MultidimAssembly &assembler)
Definition: darcy_flow_lmh.cc:1000
LinSys
Definition: la_linsys_new.hh:169
vector_mpi.hh
DarcyLMH::output_object
DarcyFlowMHOutput * output_object
Definition: darcy_flow_lmh.hh:294
time_governor.hh
Basic time management class.
LocalSystem
Definition: local_system.hh:45
DarcyLMH::EqData
Definition: darcy_flow_lmh.hh:152
DarcyLMH::EqData::p_edge_solution
VectorMPI p_edge_solution
Definition: darcy_flow_lmh.hh:167
Input
Abstract linear system class.
Definition: balance.hh:40
DarcyLMH::accept_time_step
virtual void accept_time_step()
postprocess velocity field (add sources)
Definition: darcy_flow_lmh.cc:647
string.h
Balance
Definition: balance.hh:119
DarcyLMH::create_linear_system
void create_linear_system(Input::AbstractRecord rec)
Definition: darcy_flow_lmh.cc:885
DarcyLMH::EqData::seepage_bc_systems
std::map< LongIdx, LocalSystem > seepage_bc_systems
Definition: darcy_flow_lmh.hh:171
DarcyLMH::TYPEDEF_ERR_INFO
TYPEDEF_ERR_INFO(EI_Reason, string)
DarcyLMH::init_eq_data
void init_eq_data()
Definition: darcy_flow_lmh.cc:254
DarcyLMH::EqFields::EqFields
EqFields()
Definition: darcy_flow_lmh.cc:156
field_set.hh
DarcyLMH::max_n_it_
unsigned int max_n_it_
Definition: darcy_flow_lmh.hh:303
DarcyLMH::nonlinear_iteration_
unsigned int nonlinear_iteration_
Definition: darcy_flow_lmh.hh:304
std::vector< std::shared_ptr< AssemblyFlowBase > >
DarcyLMH::data_changed_
bool data_changed_
Definition: darcy_flow_lmh.hh:298
DarcyLMH::EqData::lin_sys_schur
std::shared_ptr< LinSys > lin_sys_schur
Definition: darcy_flow_lmh.hh:166
DarcyLMH::solve_time_step
void solve_time_step(bool output=true)
Solve the problem without moving to next time and without output.
Definition: darcy_flow_lmh.cc:493
DarcyLMH::update_solution
void update_solution() override
Definition: darcy_flow_lmh.cc:479
DarcyLMH::EqFields
Definition: darcy_flow_lmh.hh:146
DarcyLMH::get_input_type
static const Input::Type::Record & get_input_type()
Definition: darcy_flow_lmh.cc:103
type_base.hh
DarcyLMH::postprocess
virtual void postprocess()
Definition: darcy_flow_lmh.cc:997
DarcyLMH::DECLARE_EXCEPTION
DECLARE_EXCEPTION(ExcSolverDiverge,<< "Diverged nonlinear solver. Reason: "<< EI_Reason::val)
DarcyLMH::assembly_linear_system
virtual void assembly_linear_system()
Definition: darcy_flow_lmh.cc:1025
DarcyLMH::set_extra_source
void set_extra_source(const Field< 3, FieldValue< 3 >::Scalar > &extra_src)
Sets external source field (coupling with other equation).
Definition: darcy_flow_lmh.hh:210
DarcyLMH::min_n_it_
unsigned int min_n_it_
Definition: darcy_flow_lmh.hh:302
exceptions.hh
Element
Definition: elements.h:40
DarcyLMH::type_field_descriptor
static const Input::Type::Record & type_field_descriptor()
Definition: darcy_flow_lmh.cc:88
DarcyMH::EqData
Definition: darcy_flow_mh.hh:204
DarcyLMH::zero_time_step
void zero_time_step() override
Definition: darcy_flow_lmh.cc:439
EquationBase::mesh
Mesh & mesh()
Definition: equation.hh:179
DarcyLMH::assembly_mh_matrix
void assembly_mh_matrix(MultidimAssembly &assembler)
Definition: darcy_flow_lmh.cc:683
DarcyLMH::eq_fields
EqFields & eq_fields()
Definition: darcy_flow_lmh.hh:202
Intersection
Definition: intersectionquadrature.hh:146
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
DarcyLMH::eq_data
EqData & eq_data()
Definition: darcy_flow_lmh.hh:203
type_generic.hh
DarcyLMH::EqData::dh_p_
std::shared_ptr< SubDOFHandlerMultiDim > dh_p_
DOF handler represents DOFs of element pressure.
Definition: darcy_flow_lmh.hh:157
DarcyLMH::solve_nonlinear
void solve_nonlinear()
Solve method common to zero_time_step and update solution.
Definition: darcy_flow_lmh.cc:556
DarcyMH::EqFields
Definition: darcy_flow_mh.hh:143
DarcyLMH::balance_
std::shared_ptr< Balance > balance_
Definition: darcy_flow_lmh.hh:292
DarcyLMH::registrar
static const int registrar
Registrar of class to factory.
Definition: darcy_flow_lmh.hh:315
DarcyLMH::set_extra_storativity
void set_extra_storativity(const Field< 3, FieldValue< 3 >::Scalar > &extra_stor)
Sets external storarivity field (coupling with other equation).
Definition: darcy_flow_lmh.hh:206
TimeGovernor
Basic time management functionality for unsteady (and steady) solvers (class Equation).
Definition: time_governor.hh:317
DarcyLMH::EqData::EqData
EqData()
Definition: darcy_flow_lmh.cc:163
DarcyLMH::solved_time
virtual double solved_time() override
Definition: darcy_flow_lmh.cc:231
field_values.hh
Input::AbstractRecord
Accessor to the polymorphic input data of a type given by an AbstracRecord object.
Definition: accessors.hh:458
input_exception.hh
mesh.h
Input::Type::Selection
Template for classes storing finite set of named values.
Definition: type_selection.hh:65
std::map< LongIdx, LocalSystem >
DarcyLMH::initialize_specific
virtual void initialize_specific()
Definition: darcy_flow_lmh.cc:404
DarcyLMH::solution_precision
virtual double solution_precision() const
Definition: darcy_flow_lmh.cc:671
DarcyLMH::eq_data_
std::shared_ptr< EqData > eq_data_
Definition: darcy_flow_lmh.hh:307
DarcyLMH::~DarcyLMH
virtual ~DarcyLMH() override
Definition: darcy_flow_lmh.cc:1301
DarcyLMH::initialize
void initialize() override
Definition: darcy_flow_lmh.cc:298
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:182
DarcyLMH::get_component_indices_vec
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)
Definition: darcy_flow_lmh.cc:1310
DarcyLMH::EqData::time_step_
double time_step_
Definition: darcy_flow_lmh.hh:164
DarcyLMH::lin_sys_schur
LinSys & lin_sys_schur()
Getter for the linear system of the 2. Schur complement.
Definition: darcy_flow_lmh.hh:289
DarcyLMH::EqData::use_steady_assembly_
bool use_steady_assembly_
Definition: darcy_flow_lmh.hh:161
Mesh
Definition: mesh.h:359
DarcyLMH::tolerance_
double tolerance_
Definition: darcy_flow_lmh.hh:301
DarcyLMH::get_mh_mortar_selection
static const Input::Type::Selection & get_mh_mortar_selection()
Selection for enum MortarMethod.
Definition: darcy_flow_lmh.cc:80
DarcyLMH::allocate_mh_matrix
void allocate_mh_matrix()
Definition: darcy_flow_lmh.cc:710
DarcyLMH
Mixed-hybrid model of linear Darcy flow, possibly unsteady.
Definition: darcy_flow_lmh.hh:128
darcy_flow_mh.hh
mixed-hybrid model of linear Darcy flow, possibly unsteady.
DarcyLMH::initial_condition_postprocess
virtual void initial_condition_postprocess()
Definition: darcy_flow_lmh.cc:436
DarcyLMH::size
int size
Definition: darcy_flow_lmh.hh:296
DarcyLMH::EqData::p_edge_solution_previous
VectorMPI p_edge_solution_previous
Definition: darcy_flow_lmh.hh:168
DarcyLMH::DECLARE_INPUT_EXCEPTION
DECLARE_INPUT_EXCEPTION(ExcMissingTimeGovernor,<< "Missing the key 'time', obligatory for the transient problems.")
VectorMPI
Definition: vector_mpi.hh:43
Field
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
DarcyFlowMHOutput
Definition: darcy_flow_mh_output.hh:74
DarcyLMH::eq_fields_
std::shared_ptr< EqFields > eq_fields_
Definition: darcy_flow_lmh.hh:306
DarcyLMH::output_data
virtual void output_data() override
Write computed fields.
Definition: darcy_flow_lmh.cc:655
DarcyLMH::DarcyLMH
DarcyLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm=nullptr)
CREATE AND FILL GLOBAL MH MATRIX OF THE WATER MODEL.
Definition: darcy_flow_lmh.cc:185
darcy_flow_interface.hh
field.hh
DarcyFlowInterface
Definition: darcy_flow_interface.hh:15
DarcyLMH::read_initial_condition
void read_initial_condition()
Definition: darcy_flow_lmh.cc:409
FieldValue
Definition: field_values.hh:645
DarcyLMH::EqData::p_edge_solution_previous_time
VectorMPI p_edge_solution_previous_time
Definition: darcy_flow_lmh.hh:169