Flow123d  release_3.0.0-1212-g8801db3
transport_dg.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 transport_dg.hh
15  * @brief Discontinuous Galerkin method for equation of transport with dispersion.
16  * @author Jan Stebel
17  */
18 
19 #ifndef TRANSPORT_DG_HH_
20 #define TRANSPORT_DG_HH_
21 
22 #include <math.h> // for fabs
23 #include <string.h> // for memcpy
24 #include <algorithm> // for max
25 #include <boost/exception/info.hpp> // for operator<<, error_info...
26 #include <string> // for operator<<
27 #include <vector> // for vector
28 #include <armadillo>
29 #include "fem/update_flags.hh" // for operator|
30 #include "fem/mapping_p1.hh"
31 #include "fields/field_values.hh" // for FieldValue<>::Scalar
32 #include "fields/field.hh"
33 #include "fields/multi_field.hh"
34 #include "la/vector_mpi.hh"
36 #include "la/linsys.hh"
37 #include "input/accessors.hh" // for ExcAccessorForNullStorage
38 #include "input/accessors_impl.hh" // for Record::val
39 #include "input/storage.hh" // for ExcStorageTypeMismatch
40 #include "input/type_base.hh" // for Array
41 #include "input/type_generic.hh" // for Instance
42 #include "input/type_record.hh" // for Record::ExcRecordKeyNo...
43 #include "mesh/long_idx.hh" // for LongIdx
44 #include "mesh/accessors.hh" // for ElementAccessor
45 #include "mesh/elements.h" // for Element::dim, Element:...
46 #include "mesh/neighbours.h" // for Neighbour::element
47 #include "mesh/side_impl.hh" // for Side::cond, Side::cond...
48 #include "mesh/sides.h" // for SideIter
49 #include "mpi.h" // for MPI_Comm_rank
50 #include "petscmat.h" // for Mat, MatDestroy
51 #include "petscvec.h" // for Vec, VecDestroy, VecSc...
52 #include "transport/concentration_model.hh" // for ConcentrationTransport...
53 #include "transport/heat_model.hh" // for HeatTransferModel, Hea...
54 
55 class DiscreteSpace;
56 class Distribution;
57 class OutputTime;
58 class DOFHandlerMultiDim;
59 template<unsigned int dim, unsigned int spacedim> class FEValuesBase;
60 template<unsigned int dim> class FiniteElement;
61 template<unsigned int dim, unsigned int spacedim> class Mapping;
62 class Quadrature;
63 namespace Input { namespace Type { class Selection; } }
64 
65 
66 
67 /**
68  * Auxiliary container class for Finite element and related objects of all dimensions.
69  * Its purpose is to provide templated access to these objects, applicable in
70  * the assembling methods.
71  */
72 class FEObjects {
73 public:
74 
75  FEObjects(Mesh *mesh_, unsigned int fe_order);
76  ~FEObjects();
77 
78  template<unsigned int dim>
79  inline FiniteElement<dim> *fe();
80 
81  template<unsigned int dim>
82  inline FiniteElement<dim> *fe_rt();
83 
84  template<unsigned int dim>
85  inline Quadrature *q() { return q_[dim]; }
86 
87  template<unsigned int dim>
88  inline MappingP1<dim,3> *mapping();
89 
90  inline std::shared_ptr<DOFHandlerMultiDim> dh();
91 
92 private:
93 
94  /// Finite elements for the solution of the advection-diffusion equation.
99 
100  /// Finite elements for the water velocity field.
104 
105  /// Quadratures used in assembling methods.
106  Quadrature *q_[4];
107 
108  /// Auxiliary mappings of reference elements.
112 
113  std::shared_ptr<DiscreteSpace> ds_;
114 
115  /// Object for distribution of dofs.
116  std::shared_ptr<DOFHandlerMultiDim> dh_;
117 };
118 
119 
120 
121 /**
122  * @brief Transport with dispersion implemented using discontinuous Galerkin method.
123  *
124  * TransportDG implements the discontinuous Galerkin method for the transport and diffusion of substances.
125  * The concentration @f$ c_i ~[kg/m^3]@f$ of the i-th substance is governed by the advection-diffusion equation
126  * @f[
127  * \partial_t c_i + \mathbf v\cdot\nabla c_i - \mathrm{div}(D\nabla c_i) = F \mbox{ in }\Omega^d,
128  * @f]
129  * where @f$\mathbf v@f$ is the fluid velocity and @f$\Omega^d@f$ the @f$d@f$-dimensional domain, respectively.
130  * The hydrodynamic dispersivity tensor @f$\mathbf D ~[m^2/s]@f$ is given by:
131  * @f[
132  * \mathbf D = D_m\mathbf I + |\mathbf v|\left(\alpha_T\mathbf I + (\alpha_L-\alpha_T)\frac{\mathbf v\otimes\mathbf v}{|\mathbf v|^2}\right).
133  * @f]
134  * The molecular dispersivity @f$D_m~[m^2/s]@f$, as well as the longitudal and transversal dispersivity @f$\alpha_L,~\alpha_T~[m]@f$ are input parameters of the model.
135  *
136  * For lower dimensions @f$d=1,2@f$ the advection-diffusion equation is multiplied by the fracture cross-cut @f$\delta^d~[m^{3-d}]@f$.
137  *
138  * The boundary @f$\partial\Omega^d@f$ is divided into three disjoint parts @f$\Gamma^d_D\cup\Gamma^d_N\cup\Gamma^d_F@f$.
139  * We prescribe the following boundary conditions:
140  * @f{eqnarray*}{
141  * c_i^d &= c_{iD}^d &\mbox{ on }\Gamma^d_D \mbox{ (Dirichlet)},\\
142  * \mathbf D^d\nabla c_i^d\cdot\mathbf n &= 0 &\mbox{ on }\Gamma^d_N \mbox{ (Neumann)},
143  * @f}
144  * The transfer of mass through fractures is described by the transmission conditions on @f$\Gamma^d_F@f$:
145  * @f[
146  * -\mathbf D^d\nabla c_i^d\cdot\mathbf n = \sigma(c_i^d-c_i^{d-1}) + \left\{\begin{array}{cl}0 &\mbox{ if }\mathbf v^d\cdot\mathbf n\ge 0\\\mathbf v^d\cdot\mathbf n(c_i^{d-1}-c_i^d) & \mbox{ if }\mathbf v^d\cdot\mathbf n<0\end{array}\right.,\qquad
147  * F^{d-1} = (\sigma + |\mathbf v^d\cdot\mathbf n|)(c_i^d-c_i^{d-1}).
148  * @f]
149  * Here @f$\mathbf n@f$ stands for the unit outward normal vector to @f$\partial\Omega^d@f$.
150  * The coefficient @f$\sigma@f$ determines the transfer of mass through fractures due to diffusion.
151  *
152  * @ingroup transport_mod
153  *
154  */
155 template<class Model>
156 class TransportDG : public Model
157 {
158 public:
159 
160  class EqData : public Model::ModelEqData {
161  public:
162 
163  EqData();
164 
165  MultiField<3, FieldValue<3>::Scalar> fracture_sigma; ///< Transition parameter for diffusive transfer on fractures (for each substance).
166  MultiField<3, FieldValue<3>::Scalar> dg_penalty; ///< Penalty enforcing inter-element continuity of solution (for each substance).
169 
171 
172  };
173 
174 
175 
176  enum DGVariant {
177  // Non-symmetric weighted interior penalty DG
178  non_symmetric = -1,
179 
180  // Incomplete weighted interior penalty DG
181  incomplete = 0,
182 
183  // Symmetric weighted interior penalty DG
184  symmetric = 1
185  };
186 
187  /**
188  * @brief Constructor.
189  * @param init_mesh computational mesh
190  * @param in_rec input record
191  */
192  TransportDG(Mesh &init_mesh, const Input::Record in_rec);
193  /**
194 
195  * @brief Declare input record type for the equation TransportDG.
196  */
197  static const Input::Type::Record & get_input_type();
198 
199  /**
200  * @brief Input type for the DG variant selection.
201  */
202  static const Input::Type::Selection & get_dg_variant_selection_input_type();
203 
204  /**
205  * @brief Initialize solution in the zero time.
206  */
207  void zero_time_step() override;
208 
209  bool evaluate_time_constraint(double &time_constraint)
210  { return false; }
211 
212  /**
213  * @brief Computes the solution in one time instant.
214  */
215  void update_solution() override;
216 
217  /**
218  * @brief Postprocesses the solution and writes to output file.
219  */
220  void output_data();
221 
222  /**
223  * @brief Destructor.
224  */
225  ~TransportDG();
226 
227  void initialize() override;
228 
229  void calculate_cumulative_balance();
230 
231  const Vec &get_solution(unsigned int sbi)
232  { return ls[sbi]->get_solution(); }
233 
235  { return solution_elem_; }
236 
237  void calculate_concentration_matrix();
238 
239  void update_after_reactions(bool solution_changed);
240 
241  void get_par_info(LongIdx * &el_4_loc, Distribution * &el_ds);
242 
243  LongIdx *get_row_4_el();
244 
245 
246 
247 
248 private:
249  /// Registrar of class to factory
250  static const int registrar;
251 
252  inline typename Model::ModelEqData &data() { return data_; }
253 
254  void preallocate();
255 
256  /**
257  * @brief Assembles the mass matrix.
258  *
259  * The routine just calls templated method assemble_mass_matrix() for each
260  * space dimension.
261  */
262  void assemble_mass_matrix();
263 
264  /**
265  * @brief Assembles the mass matrix for the given dimension.
266  */
267  template<unsigned int dim>
268  void assemble_mass_matrix();
269 
270  /**
271  * @brief Assembles the stiffness matrix.
272  *
273  * This routine just calls assemble_volume_integrals(), assemble_fluxes_boundary(),
274  * assemble_fluxes_element_element() and assemble_fluxes_element_side() for each
275  * space dimension.
276  */
277  void assemble_stiffness_matrix();
278 
279  /**
280  * @brief Assembles the volume integrals into the stiffness matrix.
281  */
282  template<unsigned int dim>
283  void assemble_volume_integrals();
284 
285  /**
286  * @brief Assembles the right hand side due to volume sources.
287  *
288  * This method just calls set_sources() for each space dimension.
289  */
290  void set_sources();
291 
292  /**
293  * @brief Assembles the right hand side vector due to volume sources.
294  */
295  template<unsigned int dim>
296  void set_sources();
297 
298  /**
299  * @brief Assembles the fluxes on the boundary.
300  */
301  template<unsigned int dim>
302  void assemble_fluxes_boundary();
303 
304  /**
305  * @brief Assembles the fluxes between elements of the same dimension.
306  */
307  template<unsigned int dim>
308  void assemble_fluxes_element_element();
309 
310  /**
311  * @brief Assembles the fluxes between elements of different dimensions.
312  */
313  template<unsigned int dim>
314  void assemble_fluxes_element_side();
315 
316 
317  /**
318  * @brief Assembles the r.h.s. components corresponding to the Dirichlet boundary conditions.
319  *
320  * The routine just calls templated method set_boundary_condition() for each space dimension.
321  */
322  void set_boundary_conditions();
323 
324  /**
325  * @brief Assembles the r.h.s. components corresponding to the Dirichlet boundary conditions
326  * for a given space dimension.
327  */
328  template<unsigned int dim>
329  void set_boundary_conditions();
330 
331  /**
332  * @brief Calculates the velocity field on a given @p dim dimensional cell.
333  *
334  * @param cell The cell.
335  * @param velocity The computed velocity field (at quadrature points).
336  * @param fv The FEValues class providing the quadrature points
337  * and the shape functions for velocity.
338  */
339  template<unsigned int dim>
340  void calculate_velocity(const ElementAccessor<3> &cell, std::vector<arma::vec3> &velocity, FEValuesBase<dim,3> &fv);
341 
342  /**
343  * @brief Calculates the dispersivity (diffusivity) tensor from the velocity field.
344  *
345  * @param K The computed dispersivity tensor.
346  * @param velocity The velocity field (at quadrature points).
347  * @param Dm Molecular diffusivities.
348  * @param alphaL Longitudal dispersivities.
349  * @param alphaT Transversal dispersivities.
350  * @param porosity Porosities.
351  * @param cross_cut Cross-cuts of higher dimension.
352  */
353 // void calculate_dispersivity_tensor(arma::mat33 &K, const arma::vec3 &velocity,
354 // double Dm, double alphaL, double alphaT, double porosity,
355 // double cross_cut);
356 
357  /**
358  * @brief Sets up parameters of the DG method on a given boundary edge.
359  *
360  * Assumption is that the edge consists of only 1 side.
361  * @param side The boundary side.
362  * @param K_size Size of vector of tensors K.
363  * @param K Dispersivity tensor.
364  * @param ad_vector Advection vector.
365  * @param normal_vector Normal vector (assumed constant along the edge).
366  * @param alpha Penalty parameter that influences the continuity
367  * of the solution (large value=more continuity).
368  * @param gamma Computed penalty parameters.
369  */
370  void set_DG_parameters_boundary(Side side,
371  const int K_size,
372  const std::vector<arma::mat33> &K,
373  const double flux,
374  const arma::vec3 &normal_vector,
375  const double alpha,
376  double &gamma);
377 
378 
379  /**
380  * @brief Sets the initial condition.
381  */
382  void set_initial_condition();
383 
384  /**
385  * @brief Assembles the auxiliary linear system to calculate the initial solution
386  * as L^2-projection of the prescribed initial condition.
387  */
388  template<unsigned int dim>
389  void prepare_initial_condition();
390 
391 
392 
393  /// @name Physical parameters
394  // @{
395 
396  /// Field data for model parameters.
398 
399  // @}
400 
401 
402  /// @name Parameters of the numerical method
403  // @{
404 
405  /// Finite element objects
407 
408  /// Penalty parameters.
410 
411  /// DG variant ((non-)symmetric/incomplete
413 
414  /// Polynomial order of finite elements.
415  unsigned int dg_order;
416 
417  // @}
418 
419 
420 
421  /// @name Solution of algebraic system
422  // @{
423 
424  /// Vector of right hand side.
426 
427  /// The stiffness matrix.
429 
430  /// The mass matrix.
432 
433  /// Mass from previous time instant (necessary when coefficients of mass matrix change in time).
435 
436  /// Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
438 
439  /// Linear algebra system for the transport equation.
441 
442  /// Linear algebra system for the time derivative (actually it is used only for handling the matrix structures).
444 
445  /// Element averages of solution (the array is passed to reactions in operator splitting).
446  double **solution_elem_;
447 
448  // @}
449 
450 
451  /// @name Output to file
452  // @{
453 
454  /// Array for storing the output solution data.
455  //vector<double*> output_solution;
456 
457  /// Vector of solution data.
459 
460  /// Record with input specification.
462 
463 
464  // @}
465 
466 
467  /// @name Auxiliary fields used during assembly
468  // @{
469 
470  /// Mass matrix coefficients.
472  /// Retardation coefficient due to sorption.
474  /// Temporary values of increments due to retardation (e.g. sorption)
476  /// Advection coefficients.
478  /// Diffusion coefficients.
480  /// Advection coefficients on edges.
482  /// Diffusion coefficients on edges.
484 
485  // @}
486 
487 
488 
489 
490  /// @name Other
491  // @{
492 
493  /// Indicates whether matrices have been preallocated.
495 
496  // @}
497 };
498 
499 
500 
501 
502 
503 
504 #endif /* TRANSPORT_DG_HH_ */
Definition: sides.h:39
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: long_idx.hh:22
Input::Record input_rec
Record with input specification.
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell...
FiniteElement< 2 > * fe_rt2_
vector< double > mm_coef
Mass matrix coefficients.
MappingP1< 3, 3 > * map3_
Transport with dispersion implemented using discontinuous Galerkin method.
Abstract linear system class.
Definition: balance.hh:37
Wrappers for linear systems based on MPIAIJ and MATIS format.
Class template representing a field with values dependent on: point, element, and region...
Definition: field.hh:83
vector< VectorMPI > output_vec
Array for storing the output solution data.
int dg_variant
DG variant ((non-)symmetric/incomplete.
Definition: mesh.h:76
vector< vector< arma::vec3 > > ad_coef
Advection coefficients.
vector< vector< vector< arma::mat33 > > > dif_coef_edg
Diffusion coefficients on edges.
Discontinuous Galerkin method for equation of transport with dispersion.
const Vec & get_solution(unsigned int sbi)
vector< double > ret_sources_prev
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell...
EquationOutput output_fields
std::vector< Mat > mass_matrix
The mass matrix.
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
vector< vector< double > > ret_coef
Retardation coefficient due to sorption.
vector< vector< vector< arma::vec3 > > > ad_coef_edg
Advection coefficients on edges.
double ** get_concentration_matrix()
FiniteElement< 0 > * fe0_
Finite elements for the solution of the advection-diffusion equation.
Definition: transport_dg.hh:95
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
vector< vector< arma::mat33 > > dif_coef
Diffusion coefficients.
FEObjects * feo
Finite element objects.
static const int registrar
Registrar of class to factory.
FiniteElement< 2 > * fe2_
Definition: transport_dg.hh:97
Abstract class for the mapping between reference and actual cell.
Definition: fe_values.hh:38
Quadrature * q()
Definition: transport_dg.hh:85
Field< 3, FieldValue< 3 >::Scalar > subdomain
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
Provides the numbering of the finite element degrees of freedom on the computational mesh...
Definition: dofhandler.hh:152
std::shared_ptr< DiscreteSpace > ds_
std::vector< Vec > rhs
Vector of right hand side.
bool evaluate_time_constraint(double &time_constraint)
Field< 3, FieldValue< 3 >::Scalar > region_id
unsigned int dg_order
Polynomial order of finite elements.
The class for outputting data during time.
Definition: output_time.hh:50
FiniteElement< 1 > * fe_rt1_
Finite elements for the water velocity field.
MappingP1< 2, 3 > * map2_
double ** solution_elem_
Element averages of solution (the array is passed to reactions in operator splitting).
Model::ModelEqData & data()
std::vector< Mat > stiffness_matrix
The stiffness matrix.
EqData data_
Field data for model parameters.
bool allocation_done
Indicates whether matrices have been preallocated.
std::vector< std::vector< double > > gamma
Penalty parameters.
FiniteElement< 1 > * fe1_
Definition: transport_dg.hh:96
std::vector< Vec > ret_vec
Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
FiniteElement< 3 > * fe3_
Definition: transport_dg.hh:98
FiniteElement< 3 > * fe_rt3_
Discontinuous Galerkin method for equation of transport with dispersion.
std::vector< Vec > mass_vec
Mass from previous time instant (necessary when coefficients of mass matrix change in time)...
Record type proxy class.
Definition: type_record.hh:182
Class for representation of a vector of fields of the same physical quantity.
Definition: multi_field.hh:87
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
Base class for FEValues and FESideValues.
Definition: fe_values.hh:37
MappingP1< 1, 3 > * map1_
Auxiliary mappings of reference elements.
LinSys ** ls
Linear algebra system for the transport equation.
LinSys ** ls_dt
Linear algebra system for the time derivative (actually it is used only for handling the matrix struc...
Template for classes storing finite set of named values.
std::shared_ptr< DOFHandlerMultiDim > dh_
Object for distribution of dofs.