Flow123d  master-f44eb46
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;
58 template<unsigned int dim, class Model> class AssemblyDG;
59 template<unsigned int dim, class Model> class MassAssemblyDG;
60 template<unsigned int dim, class Model> class StiffnessAssemblyDG;
61 template<unsigned int dim, class Model> class SourcesAssemblyDG;
62 template<unsigned int dim, class Model> class BdrConditionAssemblyDG;
63 template<unsigned int dim, class Model> class InitConditionAssemblyDG;
64 template<unsigned int dim, class Model> class InitProjectionAssemblyDG;
65 template< template<IntDim...> class DimAssembly> class GenericAssembly;
66 template<unsigned int dim, unsigned int spacedim> class FEValuesBase;
67 template<unsigned int dim> class FiniteElement;
68 template<unsigned int dim, unsigned int spacedim> class Mapping;
69 class Quadrature;
70 namespace Input { namespace Type { class Selection; } }
71 class ElementCacheMap;
72 
73 /*class FEObjects {
74 public:
75 
76  inline FEObjects(Mesh *mesh_, unsigned int fe_order)
77  {
78  fe = MixedPtr<FE_P_disc>(fe_order);
79  fe_rt = MixedPtr<FE_RT0>();
80  q = MixedPtr<QGauss>(2*fe_order);
81 
82  auto ds = std::make_shared<EqualOrderDiscreteSpace>(mesh_, fe);
83  dh = std::make_shared<DOFHandlerMultiDim>(*mesh_);
84  dh->distribute_dofs(ds_);
85  }
86 
87  ~FEObjects();
88 
89  MixedPtr<FiniteElement> fe;
90  MixedPtr<FiniteElement> fe_rt;
91  MixedPtr<Quadrature> q;
92 
93  /// Object for distribution of dofs.
94  std::shared_ptr<DOFHandlerMultiDim> dh;
95 };*/
96 
97 
98 
99 /**
100  * @brief Transport with dispersion implemented using discontinuous Galerkin method.
101  *
102  * TransportDG implements the discontinuous Galerkin method for the transport and diffusion of substances.
103  * The concentration @f$ c_i ~[kg/m^3]@f$ of the i-th substance is governed by the advection-diffusion equation
104  * @f[
105  * \partial_t c_i + \mathbf v\cdot\nabla c_i - \mathrm{div}(D\nabla c_i) = F \mbox{ in }\Omega^d,
106  * @f]
107  * where @f$\mathbf v@f$ is the fluid velocity and @f$\Omega^d@f$ the @f$d@f$-dimensional domain, respectively.
108  * The hydrodynamic dispersivity tensor @f$\mathbf D ~[m^2/s]@f$ is given by:
109  * @f[
110  * \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).
111  * @f]
112  * 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.
113  *
114  * 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$.
115  *
116  * 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$.
117  * We prescribe the following boundary conditions:
118  * @f{eqnarray*}{
119  * c_i^d &= c_{iD}^d &\mbox{ on }\Gamma^d_D \mbox{ (Dirichlet)},\\
120  * \mathbf D^d\nabla c_i^d\cdot\mathbf n &= 0 &\mbox{ on }\Gamma^d_N \mbox{ (Neumann)},
121  * @f}
122  * The transfer of mass through fractures is described by the transmission conditions on @f$\Gamma^d_F@f$:
123  * @f[
124  * -\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
125  * F^{d-1} = (\sigma + |\mathbf v^d\cdot\mathbf n|)(c_i^d-c_i^{d-1}).
126  * @f]
127  * Here @f$\mathbf n@f$ stands for the unit outward normal vector to @f$\partial\Omega^d@f$.
128  * The coefficient @f$\sigma@f$ determines the transfer of mass through fractures due to diffusion.
129  *
130  * @ingroup transport_mod
131  *
132  */
133 template<class Model>
134 class TransportDG : public Model
135 {
136 public:
137 
138  template<unsigned int dim> using MassAssemblyDim = MassAssemblyDG<dim, Model>;
139  template<unsigned int dim> using StiffnessAssemblyDim = StiffnessAssemblyDG<dim, Model>;
140  template<unsigned int dim> using SourcesAssemblyDim = SourcesAssemblyDG<dim, Model>;
144 
146 
147  class EqFields : public Model::ModelEqFields {
148  public:
149 
150  EqFields();
151 
152  MultiField<3, FieldValue<3>::Scalar> fracture_sigma; ///< Transition parameter for diffusive transfer on fractures (for each substance).
153  MultiField<3, FieldValue<3>::Scalar> dg_penalty; ///< Penalty enforcing inter-element continuity of solution (for each substance).
156 
158  };
159 
160 
161  class EqData : public Model::ModelEqData {
162  public:
163 
164  EqData() : ls(nullptr) {}
165 
166  inline void set_time_governor(TimeGovernor *time) {
167  ASSERT_PTR(time);
168  this->time_ = time;
169  }
170 
171  /// @name Parameters of the numerical method
172  // @{
173 
174  /// DG variant ((non-)symmetric/incomplete
176 
177  /// Polynomial order of finite elements.
178  unsigned int dg_order;
179 
180  // @}
181 
182 
183  /// Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
185 
186  /// Linear algebra system for the transport equation.
188 
189  /// Linear algebra system for the time derivative (actually it is used only for handling the matrix structures).
191 
192  /// @name Auxiliary fields used during assembly
193  // @{
194 
195  /// Maximal number of edge sides (evaluate from dim 1,2,3)
196  unsigned int max_edg_sides;
197 
198  // @}
199 
200  /// Object for distribution of dofs.
201  std::shared_ptr<DOFHandlerMultiDim> dh_;
202 
203  /// Vector of solution data.
205 
207  std::shared_ptr<DOFHandlerMultiDim> dh_p0;
209  std::shared_ptr<Balance> balance_;
210  };
211 
212  enum DGVariant {
213  // Non-symmetric weighted interior penalty DG
215 
216  // Incomplete weighted interior penalty DG
218 
219  // Symmetric weighted interior penalty DG
220  symmetric = 1
221  };
222 
223  /**
224  * @brief Constructor.
225  * @param init_mesh computational mesh
226  * @param in_rec input record
227  */
228  TransportDG(Mesh &init_mesh, const Input::Record in_rec);
229  /**
230 
231  * @brief Declare input record type for the equation TransportDG.
232  */
233  static const Input::Type::Record & get_input_type();
234 
235  /**
236  * @brief Input type for the DG variant selection.
237  */
239 
240  /**
241  * @brief Initialize solution in the zero time.
242  */
243  void zero_time_step() override;
244 
246  { return false; }
247 
248  /**
249  * @brief Computes the solution in one time instant.
250  */
251  void update_solution() override;
252 
253  /**
254  * @brief Postprocesses the solution and writes to output file.
255  */
256  void output_data();
257 
258  /**
259  * @brief Destructor.
260  */
261  ~TransportDG() override;
262 
263  void initialize() override;
264 
266 
267  /// Return PETSc vector with solution for sbi-th component.
268  Vec get_component_vec(unsigned int sbi)
269  { return eq_data_->ls[sbi]->get_solution(); }
270 
271  /// Getter for P0 interpolation by FieldFE.
273  { return eq_data_->conc_fe;}
274 
275  /// Compute P0 interpolation of the solution (used in reaction term).
277 
278  void update_after_reactions(bool solution_changed);
279 
280  /// Access to balance object of Model
281  inline std::shared_ptr<Balance> balance() const {
282  return Model::balance_;
283  }
284 
285  inline typename Model::ModelEqFields &eq_fields() { return *eq_fields_; }
286 
287  inline typename Model::ModelEqData &eq_data() { return *eq_data_; }
288 
289 private:
290  /// Registrar of class to factory
291  static const int registrar;
292 
293  void preallocate();
294 
295  /**
296  * @brief Calculates the dispersivity (diffusivity) tensor from the velocity field.
297  *
298  * @param K The computed dispersivity tensor.
299  * @param velocity The velocity field (at quadrature points).
300  * @param Dm Molecular diffusivities.
301  * @param alphaL Longitudal dispersivities.
302  * @param alphaT Transversal dispersivities.
303  * @param porosity Porosities.
304  * @param cross_cut Cross-cuts of higher dimension.
305  */
306 // void calculate_dispersivity_tensor(arma::mat33 &K, const arma::vec3 &velocity,
307 // double Dm, double alphaL, double alphaT, double porosity,
308 // double cross_cut);
309 
310  /**
311  * @brief Sets the initial condition.
312  *
313  * This method just calls AssemblyDG::prepare_initial_condition() for each elements.
314  */
315  void set_initial_condition();
316 
317 
318 
320 
321 
322 
323  /// @name Physical parameters
324  // @{
325 
326  /// Fields for model parameters.
327  std::shared_ptr<EqFields> eq_fields_;
328 
329  /// Data for model parameters.
330  std::shared_ptr<EqData> eq_data_;
331 
332  // @}
333 
334 
335 
336  /// @name Solution of algebraic system
337  // @{
338 
339  /// Vector of right hand side.
341 
342  /// The stiffness matrix.
344 
345  /// The mass matrix.
347 
348  /// Mass from previous time instant (necessary when coefficients of mass matrix change in time).
350  // @}
351 
352 
353  /// @name Output to file
354  // @{
355 
356  /// Array for storing the output solution data.
357  //vector<double*> output_solution;
358 
359  /// Record with input specification.
361 
362  ofstream reg_stat_stream;
363 
364 
365  // @}
366 
367 
368  /// @name Auxiliary fields used during assembly
369  // @{
370 
371  /// Temporary values of increments due to retardation (e.g. sorption)
373 
374  // @}
375 
376 
377 
378 
379  /// @name Other
380  // @{
381 
382  /// Indicates whether matrices have been preallocated.
384 
386  // @}
387 
388  /// general assembly objects, hold assembly objects of appropriate dimension
394 };
395 
396 
397 
398 
399 
400 
401 #endif /* TRANSPORT_DG_HH_ */
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
Definition: asserts.hh:341
Provides the numbering of the finite element degrees of freedom on the computational mesh.
Definition: dofhandler.hh:151
Directing class of FieldValueCache.
Class template representing a field with values dependent on: point, element, and region.
Definition: field.hh:92
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
Generic class of assemblation.
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
Class for representation of a vector of fields of the same physical quantity.
Definition: multi_field.hh:87
The class for outputting data during time.
Definition: output_time.hh:51
Base class for quadrature rules on simplices in arbitrary dimensions.
Definition: quadrature.hh:48
Basic time management functionality for unsteady (and steady) solvers (class Equation).
FieldFEScalarVec conc_fe
std::shared_ptr< DOFHandlerMultiDim > dh_
Object for distribution of dofs.
std::vector< VectorMPI > output_vec
Vector of solution data.
LinSys ** ls_dt
Linear algebra system for the time derivative (actually it is used only for handling the matrix struc...
std::shared_ptr< DOFHandlerMultiDim > dh_p0
int dg_variant
DG variant ((non-)symmetric/incomplete.
TimeGovernor * time_
void set_time_governor(TimeGovernor *time)
std::shared_ptr< Balance > balance_
unsigned int dg_order
Polynomial order of finite elements.
LinSys ** ls
Linear algebra system for the transport equation.
std::vector< Vec > ret_vec
Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
unsigned int max_edg_sides
Maximal number of edge sides (evaluate from dim 1,2,3)
Field< 3, FieldValue< 3 >::Scalar > region_id
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
EquationOutput output_fields
Field< 3, FieldValue< 3 >::Scalar > subdomain
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
Transport with dispersion implemented using discontinuous Galerkin method.
std::vector< std::shared_ptr< FieldFE< 3, FieldValue< 3 >::Scalar > > > FieldFEScalarVec
void update_solution() override
Computes the solution in one time instant.
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
Definition: transport_dg.cc:78
std::vector< Vec > mass_vec
Mass from previous time instant (necessary when coefficients of mass matrix change in time).
bool allocation_done
Indicates whether matrices have been preallocated.
void output_region_statistics()
Vec get_component_vec(unsigned int sbi)
Return PETSc vector with solution for sbi-th component.
Input::Record input_rec
Array for storing the output solution data.
void zero_time_step() override
Initialize solution in the zero time.
std::shared_ptr< EqFields > eq_fields_
Fields for model parameters.
void initialize() override
GenericAssembly< MassAssemblyDim > * mass_assembly_
general assembly objects, hold assembly objects of appropriate dimension
static const Input::Type::Selection & get_dg_variant_selection_input_type()
Input type for the DG variant selection.
Definition: transport_dg.cc:52
GenericAssemblyBase * init_assembly_
vector< double > ret_sources
Temporary values of increments due to retardation (e.g. sorption)
FieldFEScalarVec & get_p0_interpolation()
Getter for P0 interpolation by FieldFE.
GenericAssembly< SourcesAssemblyDim > * sources_assembly_
GenericAssembly< BdrConditionAssemblyDim > * bdr_cond_assembly_
Model::ModelEqData & eq_data()
TransportDG(Mesh &init_mesh, const Input::Record in_rec)
Constructor.
~TransportDG() override
Destructor.
std::vector< Vec > rhs
Vector of right hand side.
bool init_projection
GenericAssembly< StiffnessAssemblyDim > * stiffness_assembly_
vector< double > ret_sources_prev
std::vector< Mat > stiffness_matrix
The stiffness matrix.
void compute_p0_interpolation()
Compute P0 interpolation of the solution (used in reaction term).
void output_data()
Postprocesses the solution and writes to output file.
static const int registrar
Registrar of class to factory.
ofstream reg_stat_stream
void calculate_cumulative_balance()
std::shared_ptr< Balance > balance() const
Access to balance object of Model.
std::vector< Mat > mass_matrix
The mass matrix.
void preallocate()
void set_initial_condition()
Calculates the dispersivity (diffusivity) tensor from the velocity field.
std::shared_ptr< EqData > eq_data_
Data for model parameters.
void update_after_reactions(bool solution_changed)
Model::ModelEqFields & eq_fields()
bool evaluate_time_constraint(double &)
Discontinuous Galerkin method for equation of transport with dispersion.
Discontinuous Galerkin method for equation of transport with dispersion.
Wrappers for linear systems based on MPIAIJ and MATIS format.
unsigned int IntDim
A dimension index type.
Definition: mixed.hh:19
Abstract linear system class.
Definition: balance.hh:40
Enum type UpdateFlags indicates which quantities are to be recomputed on each finite element cell.