Flow123d  JS_before_hm-984-g3a19f2f
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 
26 #include <string> // for operator<<
27 #include <vector> // for vector
28 #include <armadillo>
29 #include "fem/update_flags.hh" // for operator|
30 #include "fields/field_values.hh" // for FieldValue<>::Scalar
31 #include "fields/field.hh"
32 #include "fields/field_fe.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 "system/index_types.hh" // for LongIdx
44 #include "mesh/accessors.hh" // for ElementAccessor, SIdeIter
45 #include "mesh/elements.h" // for Element::dim, Element:...
46 #include "mesh/neighbours.h" // for Neighbour::element
47 #include "mpi.h" // for MPI_Comm_rank
48 #include "petscmat.h" // for Mat, MatDestroy
49 #include "petscvec.h" // for Vec, VecDestroy, VecSc...
50 #include "transport/concentration_model.hh" // for ConcentrationTransport...
51 #include "transport/heat_model.hh" // for HeatTransferModel, Hea...
52 #include "tools/mixed.hh"
53 class DiscreteSpace;
54 class Distribution;
55 class OutputTime;
56 class DOFHandlerMultiDim;
57 template<unsigned int dim, class Model> class AssemblyDG;
58 template<unsigned int dim, class Model> class MassAssemblyDG;
59 template<unsigned int dim, class Model> class StiffnessAssemblyDG;
60 template<unsigned int dim, class Model> class SourcesAssemblyDG;
61 template<unsigned int dim, class Model> class BdrConditionAssemblyDG;
62 template<unsigned int dim, class Model> class InitConditionAssemblyDG;
63 template< template<IntDim...> class DimAssembly> class GenericAssembly;
64 template<unsigned int dim, unsigned int spacedim> class FEValuesBase;
65 template<unsigned int dim> class FiniteElement;
66 template<unsigned int dim, unsigned int spacedim> class Mapping;
67 class Quadrature;
68 namespace Input { namespace Type { class Selection; } }
69 class ElementCacheMap;
70 
71 /*class FEObjects {
72 public:
73 
74  inline FEObjects(Mesh *mesh_, unsigned int fe_order)
75  {
76  fe = MixedPtr<FE_P_disc>(fe_order);
77  fe_rt = MixedPtr<FE_RT0>();
78  q = MixedPtr<QGauss>(2*fe_order);
79 
80  auto ds = std::make_shared<EqualOrderDiscreteSpace>(mesh_, fe);
81  dh = std::make_shared<DOFHandlerMultiDim>(*mesh_);
82  dh->distribute_dofs(ds_);
83  }
84 
85  ~FEObjects();
86 
87  MixedPtr<FiniteElement> fe;
88  MixedPtr<FiniteElement> fe_rt;
89  MixedPtr<Quadrature> q;
90 
91  /// Object for distribution of dofs.
92  std::shared_ptr<DOFHandlerMultiDim> dh;
93 };*/
94 
95 
96 
97 /**
98  * @brief Transport with dispersion implemented using discontinuous Galerkin method.
99  *
100  * TransportDG implements the discontinuous Galerkin method for the transport and diffusion of substances.
101  * The concentration @f$ c_i ~[kg/m^3]@f$ of the i-th substance is governed by the advection-diffusion equation
102  * @f[
103  * \partial_t c_i + \mathbf v\cdot\nabla c_i - \mathrm{div}(D\nabla c_i) = F \mbox{ in }\Omega^d,
104  * @f]
105  * where @f$\mathbf v@f$ is the fluid velocity and @f$\Omega^d@f$ the @f$d@f$-dimensional domain, respectively.
106  * The hydrodynamic dispersivity tensor @f$\mathbf D ~[m^2/s]@f$ is given by:
107  * @f[
108  * \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).
109  * @f]
110  * 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.
111  *
112  * 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$.
113  *
114  * 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$.
115  * We prescribe the following boundary conditions:
116  * @f{eqnarray*}{
117  * c_i^d &= c_{iD}^d &\mbox{ on }\Gamma^d_D \mbox{ (Dirichlet)},\\
118  * \mathbf D^d\nabla c_i^d\cdot\mathbf n &= 0 &\mbox{ on }\Gamma^d_N \mbox{ (Neumann)},
119  * @f}
120  * The transfer of mass through fractures is described by the transmission conditions on @f$\Gamma^d_F@f$:
121  * @f[
122  * -\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
123  * F^{d-1} = (\sigma + |\mathbf v^d\cdot\mathbf n|)(c_i^d-c_i^{d-1}).
124  * @f]
125  * Here @f$\mathbf n@f$ stands for the unit outward normal vector to @f$\partial\Omega^d@f$.
126  * The coefficient @f$\sigma@f$ determines the transfer of mass through fractures due to diffusion.
127  *
128  * @ingroup transport_mod
129  *
130  */
131 template<class Model>
132 class TransportDG : public Model
133 {
134 public:
135 
136  template<unsigned int dim> using MassAssemblyDim = MassAssemblyDG<dim, Model>;
137  template<unsigned int dim> using StiffnessAssemblyDim = StiffnessAssemblyDG<dim, Model>;
138  template<unsigned int dim> using SourcesAssemblyDim = SourcesAssemblyDG<dim, Model>;
141 
143 
144  class EqData : public Model::ModelEqData {
145  public:
146 
147  EqData();
148 
149  /**
150  * @brief Sets up parameters of the DG method on a given boundary edge.
151  *
152  * Assumption is that the edge consists of only 1 side.
153  * @param side The boundary side.
154  * @param K_size Size of vector of tensors K.
155  * @param K Dispersivity tensor.
156  * @param ad_vector Advection vector.
157  * @param normal_vector Normal vector (assumed constant along the edge).
158  * @param alpha Penalty parameter that influences the continuity
159  * of the solution (large value=more continuity).
160  * @param gamma Computed penalty parameters.
161  */
162  void set_DG_parameters_boundary(Side side,
163  const int K_size,
164  const std::vector<arma::mat33> &K,
165  const double flux,
166  const arma::vec3 &normal_vector,
167  const double alpha,
168  double &gamma);
169 
170 
171  /// Compute and return anisotropy of given element
172  double elem_anisotropy(ElementAccessor<3> e) const;
173 
174 
175 
176  MultiField<3, FieldValue<3>::Scalar> fracture_sigma; ///< Transition parameter for diffusive transfer on fractures (for each substance).
177  MultiField<3, FieldValue<3>::Scalar> dg_penalty; ///< Penalty enforcing inter-element continuity of solution (for each substance).
180 
182 
183 
184  /// @name Parameters of the numerical method
185  // @{
186 
187  /// Penalty parameters.
189 
190  /// DG variant ((non-)symmetric/incomplete
192 
193  /// Polynomial order of finite elements.
194  unsigned int dg_order;
195 
196  // @}
197 
198 
199  /// Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
201 
202  /// Linear algebra system for the transport equation.
204 
205  /// Linear algebra system for the time derivative (actually it is used only for handling the matrix structures).
207 
208  /// @name Auxiliary fields used during assembly
209  // @{
210 
211  /// Advection coefficients.
213  /// Diffusion coefficients.
215  /// Advection coefficients on edges.
217  /// Diffusion coefficients on edges.
219 
220  // @}
221 
222  /// Object for distribution of dofs.
223  std::shared_ptr<DOFHandlerMultiDim> dh_;
224 
225  /// general assembly objects, hold assembly objects of appropriate dimension
231 
232  FieldFEScalarVec conc_fe;
233  std::shared_ptr<DOFHandlerMultiDim> dh_p0;
234  };
235 
236  enum DGVariant {
237  // Non-symmetric weighted interior penalty DG
238  non_symmetric = -1,
239 
240  // Incomplete weighted interior penalty DG
241  incomplete = 0,
242 
243  // Symmetric weighted interior penalty DG
244  symmetric = 1
245  };
246 
247  /**
248  * @brief Constructor.
249  * @param init_mesh computational mesh
250  * @param in_rec input record
251  */
252  TransportDG(Mesh &init_mesh, const Input::Record in_rec);
253  /**
254 
255  * @brief Declare input record type for the equation TransportDG.
256  */
257  static const Input::Type::Record & get_input_type();
258 
259  /**
260  * @brief Input type for the DG variant selection.
261  */
262  static const Input::Type::Selection & get_dg_variant_selection_input_type();
263 
264  /**
265  * @brief Initialize solution in the zero time.
266  */
267  void zero_time_step() override;
268 
270  { return false; }
271 
272  /**
273  * @brief Computes the solution in one time instant.
274  */
275  void update_solution() override;
276 
277  /**
278  * @brief Postprocesses the solution and writes to output file.
279  */
280  void output_data();
281 
282  /**
283  * @brief Destructor.
284  */
285  ~TransportDG() override;
286 
287  void initialize() override;
288 
289  void calculate_cumulative_balance();
290 
291  /// Return PETSc vector with solution for sbi-th component.
292  Vec get_component_vec(unsigned int sbi)
293  { return data_->ls[sbi]->get_solution(); }
294 
295  /// Getter for P0 interpolation by FieldFE.
296  FieldFEScalarVec& get_p0_interpolation()
297  { return data_->conc_fe;}
298 
299  /// Compute P0 interpolation of the solution (used in reaction term).
300  void compute_p0_interpolation();
301 
302  void update_after_reactions(bool solution_changed);
303 
304  void get_par_info(LongIdx * &el_4_loc, Distribution * &el_ds);
305 
306  LongIdx *get_row_4_el();
307 
308  /// Access to balance object of Model
309  inline std::shared_ptr<Balance> balance() const {
310  return Model::balance_;
311  }
312 
313  /// Return vector of substances indices
314  inline const vector<unsigned int> subst_idx() const {
315  return Model::subst_idx;
316  }
317 
318 
319  inline typename Model::ModelEqData &data() { return *data_; }
320 
321 private:
322  /// Registrar of class to factory
323  static const int registrar;
324 
325  void preallocate();
326 
327  /**
328  * @brief Calculates the dispersivity (diffusivity) tensor from the velocity field.
329  *
330  * @param K The computed dispersivity tensor.
331  * @param velocity The velocity field (at quadrature points).
332  * @param Dm Molecular diffusivities.
333  * @param alphaL Longitudal dispersivities.
334  * @param alphaT Transversal dispersivities.
335  * @param porosity Porosities.
336  * @param cross_cut Cross-cuts of higher dimension.
337  */
338 // void calculate_dispersivity_tensor(arma::mat33 &K, const arma::vec3 &velocity,
339 // double Dm, double alphaL, double alphaT, double porosity,
340 // double cross_cut);
341 
342  /**
343  * @brief Sets the initial condition.
344  *
345  * This method just calls AssemblyDG::prepare_initial_condition() for each elements.
346  */
347  void set_initial_condition();
348 
349 
350 
351  void output_region_statistics();
352 
353 
354 
355  /// @name Physical parameters
356  // @{
357 
358  /// Field data for model parameters.
359  std::shared_ptr<EqData> data_;
360 
361  // @}
362 
363 
364 
365  /// @name Solution of algebraic system
366  // @{
367 
368  /// Vector of right hand side.
370 
371  /// The stiffness matrix.
373 
374  /// The mass matrix.
376 
377  /// Mass from previous time instant (necessary when coefficients of mass matrix change in time).
379  // @}
380 
381 
382  /// @name Output to file
383  // @{
384 
385  /// Array for storing the output solution data.
386  //vector<double*> output_solution;
387 
388  /// Vector of solution data.
390 
391  /// Record with input specification.
393 
394  ofstream reg_stat_stream;
395 
396 
397  // @}
398 
399 
400  /// @name Auxiliary fields used during assembly
401  // @{
402 
403  /// Temporary values of increments due to retardation (e.g. sorption)
405 
406  // @}
407 
408 
409 
410 
411  /// @name Other
412  // @{
413 
414  /// Indicates whether matrices have been preallocated.
416 
417  // @}
418 
419 };
420 
421 
422 
423 
424 
425 
426 #endif /* TRANSPORT_DG_HH_ */
GenericAssembly< MassAssemblyDim > * mass_assembly_
general assembly objects, hold assembly objects of appropriate dimension
Input::Record input_rec
Record with input specification.
std::vector< std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > > FieldFEScalarVec
vector< vector< arma::vec3 > > ad_coef
Advection coefficients.
std::vector< std::vector< double > > gamma
Penalty parameters.
Transport with dispersion implemented using discontinuous Galerkin method.
unsigned int dg_order
Polynomial order of finite elements.
const vector< unsigned int > subst_idx() const
Return vector of substances indices.
Abstract linear system class.
Definition: balance.hh:40
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:92
vector< VectorMPI > output_vec
Array for storing the output solution data.
Vec get_component_vec(unsigned int sbi)
Return PETSc vector with solution for sbi-th component.
FieldFEScalarVec & get_p0_interpolation()
Getter for P0 interpolation by FieldFE.
FieldFEScalarVec conc_fe
Directing class of FieldValueCache.
Definition: mesh.h:78
int dg_variant
DG variant ((non-)symmetric/incomplete.
Discontinuous Galerkin method for equation of transport with dispersion.
vector< double > ret_sources_prev
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
std::shared_ptr< DOFHandlerMultiDim > dh_p0
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.
vector< vector< vector< arma::mat33 > > > dif_coef_edg
Diffusion coefficients on edges.
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
bool evaluate_time_constraint(double &)
GenericAssembly< SourcesAssemblyDim > * sources_assembly_
vector< vector< arma::mat33 > > dif_coef
Diffusion coefficients.
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
GenericAssembly< BdrConditionAssemblyDim > * bdr_cond_assembly_
static const int registrar
Registrar of class to factory.
std::vector< Vec > ret_vec
Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
LinSys ** ls
Linear algebra system for the transport equation.
Field< 3, FieldValue< 3 >::Scalar > subdomain
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
ofstream reg_stat_stream
Provides the numbering of the finite element degrees of freedom on the computational mesh...
Definition: dofhandler.hh:151
std::vector< Vec > rhs
Vector of right hand side.
Field< 3, FieldValue< 3 >::Scalar > region_id
The class for outputting data during time.
Definition: output_time.hh:50
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
vector< vector< vector< arma::vec3 > > > ad_coef_edg
Advection coefficients on edges.
Model::ModelEqData & data()
std::vector< Mat > stiffness_matrix
The stiffness matrix.
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
bool allocation_done
Indicates whether matrices have been preallocated.
GenericAssembly< StiffnessAssemblyDim > * stiffness_assembly_
Generic class of assemblation.
Definition: assembly_dg.hh:64
std::shared_ptr< DOFHandlerMultiDim > dh_
Object for distribution of dofs.
std::shared_ptr< Balance > balance() const
Access to balance object of Model.
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...
GenericAssembly< InitConditionAssemblyDim > * init_cond_assembly_
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< EqData > data_
Field data for model parameters.