Flow123d
transport_dg.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2007 Technical University of Liberec. All rights reserved.
4  *
5  * Please make a following refer to Flow123d on your project site if you use the program for any purpose,
6  * especially for academic research:
7  * Flow123d, Research Centre: Advanced Remedial Technologies, Technical University of Liberec, Czech Republic
8  *
9  * This program is free software; you can redistribute it and/or modify it under the terms
10  * of the GNU General Public License version 3 as published by the Free Software Foundation.
11  *
12  * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
13  * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14  * See the GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License along with this program; if not,
17  * write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 021110-1307, USA.
18  *
19  *
20  * $Id$
21  * $Revision$
22  * $LastChangedBy$
23  * $LastChangedDate$
24  *
25  * @file
26  * @brief Discontinuous Galerkin method for equation of transport with dispersion.
27  * @author Jan Stebel
28  */
29 
30 #ifndef TRANSPORT_DG_HH_
31 #define TRANSPORT_DG_HH_
32 
34 #include "la/linsys.hh"
35 #include "flow/mh_dofhandler.hh"
36 
37 class Distribution;
38 class OutputTime;
39 class DOFHandlerBase;
40 class DOFHandlerMultiDim;
41 template<unsigned int dim, unsigned int spacedim> class FEValuesBase;
42 template<unsigned int dim, unsigned int spacedim> class FiniteElement;
43 template<unsigned int dim, unsigned int spacedim> class Mapping;
44 template<unsigned int dim> class Quadrature;
45 
46 
47 
48 /**
49  * Auxiliary container class for Finite element and related objects of all dimensions.
50  * Its purpose is to provide templated access to these objects, applicable in
51  * the assembling methods.
52  */
53 class FEObjects {
54 public:
55 
56  FEObjects(Mesh *mesh_, unsigned int fe_order);
57  ~FEObjects();
58 
59  template<unsigned int dim>
60  inline FiniteElement<dim,3> *fe();
61 
62  template<unsigned int dim>
63  inline FiniteElement<dim,3> *fe_rt();
64 
65  template<unsigned int dim>
66  inline Quadrature<dim> *q();
67 
68  template<unsigned int dim>
69  inline Mapping<dim,3> *mapping();
70 
71  inline DOFHandlerMultiDim *dh();
72 
73 private:
74 
75  /// Finite elements for the solution of the advection-diffusion equation.
79 
80  /// Finite elements for the water velocity field.
84 
85  /// Quadratures used in assembling methods.
90 
91  /// Auxiliary mappings of reference elements.
96 
97  /// Object for distribution of dofs.
99 };
100 
101 
102 
103 /**
104  * @brief Transport with dispersion implemented using discontinuous Galerkin method.
105  *
106  * TransportDG implements the discontinuous Galerkin method for the transport and diffusion of substances.
107  * The concentration @f$ c_i ~[kg/m^3]@f$ of the i-th substance is governed by the advection-diffusion equation
108  * @f[
109  * \partial_t c_i + \mathbf v\cdot\nabla c_i - \mathrm{div}(D\nabla c_i) = F \mbox{ in }\Omega^d,
110  * @f]
111  * where @f$\mathbf v@f$ is the fluid velocity and @f$\Omega^d@f$ the @f$d@f$-dimensional domain, respectively.
112  * The hydrodynamic dispersivity tensor @f$\mathbf D ~[m^2/s]@f$ is given by:
113  * @f[
114  * \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).
115  * @f]
116  * 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.
117  *
118  * 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$.
119  *
120  * 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$.
121  * We prescribe the following boundary conditions:
122  * @f{eqnarray*}{
123  * c_i^d &= c_{iD}^d &\mbox{ on }\Gamma^d_D \mbox{ (Dirichlet)},\\
124  * \mathbf D^d\nabla c_i^d\cdot\mathbf n &= 0 &\mbox{ on }\Gamma^d_N \mbox{ (Neumann)},
125  * @f}
126  * The transfer of mass through fractures is described by the transmission conditions on @f$\Gamma^d_F@f$:
127  * @f[
128  * -\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
129  * F^{d-1} = (\sigma + |\mathbf v^d\cdot\mathbf n|)(c_i^d-c_i^{d-1}).
130  * @f]
131  * Here @f$\mathbf n@f$ stands for the unit outward normal vector to @f$\partial\Omega^d@f$.
132  * The coefficient @f$\sigma@f$ determines the transfer of mass through fractures due to diffusion.
133  *
134  * @ingroup transport_mod
135  *
136  */
137 template<class Model>
138 class TransportDG : public TransportBase, public Model
139 {
140 public:
141 
142  class EqData : public Model::ModelEqData {
143  public:
144 
145  enum BC_Type {
150  };
152 
154 
155  EqData();
156 
157  Field<3, FieldValue<3>::Vector> fracture_sigma; ///< Transition parameter for diffusive transfer on fractures (for each substance).
158  Field<3, FieldValue<3>::Vector> dg_penalty; ///< Penalty enforcing inter-element continuity of solution (for each substance).
159 
163 
164  };
165 
166 
167 
168  enum DGVariant {
169  // Non-symmetric weighted interior penalty DG
171 
172  // Incomplete weighted interior penalty DG
174 
175  // Symmetric weighted interior penalty DG
177  };
178 
179  /**
180  * @brief Constructor.
181  * @param init_mesh computational mesh
182  * @param in_rec input record
183  */
184  TransportDG(Mesh &init_mesh, const Input::Record &in_rec);
185  /**
186 
187  * @brief Declare input record type for the equation TransportDG.
188  */
190 
191  /**
192  * @brief Input type for the DG variant selection.
193  */
195 
196  /**
197  * @brief Initialize solution in the zero time.
198  */
199  void zero_time_step() override;
200 
201  /**
202  * @brief Computes the solution in one time instant.
203  */
204  void update_solution() override;
205 
206  /**
207  * @brief Returns the serialized solution array.
208  * @param vector the solution
209  * @param size size of the array
210  */
211  void get_solution_vector(double * &vector, unsigned int &size);
212 
213  /**
214  * @brief Returns the (possibly) parallel solution vector.
215  * @param vector
216  */
218 
219  /**
220  * @brief Updates the velocity field which determines some coefficients of the transport equation.
221  *
222  * @param dh mixed hybrid dof handler
223  *
224  * (So far it does not work since the flow module returns a vector of zeros.)
225  * @param velocity_vector Input array of velocity values.
226  */
227  virtual void set_velocity_field(const MH_DofHandler &dh);
228 
229  /**
230  * @brief Postprocesses the solution and writes to output file.
231  */
232  void output_data();
233 
234  /**
235  * @brief Sets pointer to data of other equations.
236  * TODO: there should be also passed the sigma parameter between dimensions
237  * @param cross_section is pointer to cross_section data of Darcy flow equation
238  */
239  //void set_cross_section_field(const Field< 3, FieldValue<3>::Scalar > &cross_section) { Model::set_cross_section_field(cross_section); }
240 
241  /**
242  * @brief Getter for field data.
243  */
244  virtual EqData *get_data() { return &data_; }
245 
247 
248  /**
249  * @brief Destructor.
250  */
251  ~TransportDG();
252 
253 private:
254 
255  typename Model::ModelEqData &data() { return data_; }
256 
257  void output_vector_gather();
258 
260  return Model::stiffness_matrix_changed() ||
263  }
264 
265  bool mass_matrix_changed() { return Model::mass_matrix_changed(); }
266 
267  bool rhs_changed() {
268  return Model::rhs_changed() ||
270  }
271 
272  /**
273  * @brief Assembles the mass matrix.
274  *
275  * The routine just calls templated method assemble_mass_matrix() for each
276  * space dimension.
277  */
278  void assemble_mass_matrix();
279 
280  /**
281  * @brief Assembles the mass matrix for the given dimension.
282  */
283  template<unsigned int dim>
284  void assemble_mass_matrix();
285 
286  /**
287  * @brief Assembles the stiffness matrix.
288  *
289  * This routine just calls assemble_volume_integrals(), assemble_fluxes_boundary(),
290  * assemble_fluxes_element_element() and assemble_fluxes_element_side() for each
291  * space dimension.
292  */
294 
295  /**
296  * @brief Assembles the volume integrals into the stiffness matrix.
297  */
298  template<unsigned int dim>
300 
301  /**
302  * @brief Assembles the right hand side due to volume sources.
303  *
304  * This method just calls set_sources() for each space dimension.
305  */
306  void set_sources();
307 
308  /**
309  * @brief Assembles the right hand side vector due to volume sources.
310  */
311  template<unsigned int dim>
312  void set_sources();
313 
314  /**
315  * @brief Assembles the fluxes on the boundary.
316  */
317  template<unsigned int dim>
319 
320  /**
321  * @brief Assembles the fluxes between elements of the same dimension.
322  */
323  template<unsigned int dim>
325 
326  /**
327  * @brief Assembles the fluxes between elements of different dimensions.
328  */
329  template<unsigned int dim>
331 
332 
333  /**
334  * @brief Assembles the r.h.s. components corresponding to the Dirichlet boundary conditions.
335  *
336  * The routine just calls templated method set_boundary_condition() for each space dimension.
337  */
339 
340  /**
341  * @brief Assembles the r.h.s. components corresponding to the Dirichlet boundary conditions
342  * for a given space dimension.
343  */
344  template<unsigned int dim>
346 
347  /**
348  * @brief Calculates the velocity field on a given @p dim dimensional cell.
349  *
350  * @param cell The cell.
351  * @param velocity The computed velocity field (at quadrature points).
352  * @param fv The FEValues class providing the quadrature points
353  * and the shape functions for velocity.
354  */
355  template<unsigned int dim>
357 
358  /**
359  * @brief Calculates the dispersivity (diffusivity) tensor from the velocity field.
360  *
361  * @param K The computed dispersivity tensor.
362  * @param velocity The velocity field (at quadrature points).
363  * @param Dm Molecular diffusivities.
364  * @param alphaL Longitudal dispersivities.
365  * @param alphaT Transversal dispersivities.
366  * @param porosity Porosities.
367  * @param cross_cut Cross-cuts of higher dimension.
368  */
369  void calculate_dispersivity_tensor(arma::mat33 &K, const arma::vec3 &velocity,
370  double Dm, double alphaL, double alphaT, double porosity,
371  double cross_cut);
372 
373  /**
374  * @brief Sets up some parameters of the DG method for two sides of an edge.
375  *
376  * @param edg The edge.
377  * @param s1 Side 1.
378  * @param s2 Side 2.
379  * @param K_size Size of vector of tensors K.
380  * @param K1 Dispersivity tensors on side s1 (in quadrature points).
381  * @param K2 Dispersivity tensors on side s2 (in quadrature points).
382  * @param normal_vector Normal vector to side 0 of the neighbour
383  * (assumed constant along the side).
384  * @param alpha1, alpha2 Penalty parameter that influences the continuity
385  * of the solution (large value=more continuity).
386  * @param gamma Computed penalty parameters.
387  * @param omega Computed weights.
388  * @param transport_flux Computed flux from side s1 to side s2.
389  */
390  void set_DG_parameters_edge(const Edge &edg,
391  const int s1,
392  const int s2,
393  const int K_size,
394  const std::vector<arma::mat33> &K1,
395  const std::vector<arma::mat33> &K2,
396  const std::vector<double> &fluxes,
397  const arma::vec3 &normal_vector,
398  const double alpha1,
399  const double alpha2,
400  double &gamma,
401  double *omega,
402  double &transport_flux);
403 
404  /**
405  * @brief Sets up parameters of the DG method on a given boundary edge.
406  *
407  * Assumption is that the edge consists of only 1 side.
408  * @param side The boundary side.
409  * @param K_size Size of vector of tensors K.
410  * @param K Dispersivity tensor.
411  * @param ad_vector Advection vector.
412  * @param normal_vector Normal vector (assumed constant along the edge).
413  * @param alpha Penalty parameter that influences the continuity
414  * of the solution (large value=more continuity).
415  * @param gamma Computed penalty parameters.
416  */
417  void set_DG_parameters_boundary(const SideIter side,
418  const int K_size,
419  const std::vector<arma::mat33> &K,
420  const double flux,
421  const arma::vec3 &normal_vector,
422  const double alpha,
423  double &gamma);
424 
425 
426  /**
427  * @brief Sets the initial condition.
428  */
429  void set_initial_condition();
430 
431  /**
432  * @brief Assembles the auxiliary linear system to calculate the initial solution
433  * as L^2-projection of the prescribed initial condition.
434  */
435  template<unsigned int dim>
437 
438  /**
439  * @brief Calculates flux through boundary of each region.
440  *
441  * This actually calls calc_fluxes<dim>() for each space dimension.
442  * @param bcd_balance Total fluxes.
443  * @param bcd_plus_balance Incoming fluxes.
444  * @param bcd_minus_balance Outgoing fluxes.
445  */
446  void calc_fluxes(vector<vector<double> > &bcd_balance, vector<vector<double> > &bcd_plus_balance, vector<vector<double> > &bcd_minus_balance);
447 
448  /**
449  * @brief Calculates flux through boundary of each region of specific dimension.
450  * @param bcd_balance Total fluxes.
451  * @param bcd_plus_balance Incoming fluxes.
452  * @param bcd_minus_balance Outgoing fluxes.
453  */
454  template<unsigned int dim>
455  void calc_fluxes(vector<vector<double> > &bcd_balance, vector<vector<double> > &bcd_plus_balance, vector<vector<double> > &bcd_minus_balance);
456 
457  /**
458  * @brief Calculates volume sources for each region.
459  *
460  * This method actually calls calc_elem_sources<dim>() for each space dimension.
461  * @param mass Vector of substance mass per region.
462  * @param src_balance Vector of sources per region.
463  */
464  void calc_elem_sources(vector<vector<double> > &mass, vector< vector<double> > &src_balance);
465 
466  /**
467  * @brief Calculates volume sources for each region of specific dimension.
468  * @param mass Vector of substance mass per region.
469  * @param src_balance Vector of sources per region.
470  */
471  template<unsigned int dim>
472  void calc_elem_sources(vector<vector<double> > &mass, vector< vector<double> > &src_balance);
473 
474 
475 
476  /// @name Physical parameters
477  // @{
478 
479  /// Field data for model parameters.
481 
482  // @}
483 
484 
485  /// @name Parameters of the numerical method
486  // @{
487 
488  /// Finite element objects
490 
491  /// Penalty parameters.
493 
494  /// DG variant ((non-)symmetric/incomplete
496 
497  /// Polynomial order of finite elements.
498  unsigned int dg_order;
499 
500  // @}
501 
502 
503 
504  /// @name Solution of algebraic system
505  // @{
506 
507  /// Vector of right hand side.
508  Vec *rhs;
509 
510  /// The stiffness matrix.
512 
513  /// The mass matrix.
515 
516  /// Linear algebra system for the transport equation.
518 
519  /// Linear algebra system for the time derivative (actually it is used only for handling the matrix structures).
521 
522  // @}
523 
524 
525  /// @name Output to file
526  // @{
527 
528  /// Array for storing the output solution data.
530 
531  /// Vector of solution data.
533 
534  /// Record with output specification.
536 
538 
539 
540  // @}
541 
542 
543  /// @name Auxiliary fields used during assembly
544  // @{
545 
546  /// Mass matrix coefficients.
548  /// Advection coefficients.
550  /// Diffusion coefficients.
552  /// Advection coefficients on edges.
554  /// Diffusion coefficients on edges.
556 
557  // @}
558 
559 
560 
561 
562  /// @name Other
563  // @{
564 
565  /// Indicates whether matrices have been preallocated.
567 
568  // @}
569 };
570 
571 
572 
573 
574 
575 
576 #endif /* TRANSPORT_DG_HH_ */