Flow123d  release_1.8.2-2141-g34b6400
transport_dg.cc
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.cc
15  * @brief Discontinuous Galerkin method for equation of transport with dispersion.
16  * @author Jan Stebel
17  */
18 
19 #include "system/sys_profiler.hh"
20 #include "transport/transport_dg.hh"
21 
22 #include "io/output_time.hh"
23 #include "quadrature/quadrature_lib.hh"
24 #include "fem/mapping_p1.hh"
25 #include "fem/fe_values.hh"
26 #include "fem/fe_p.hh"
27 #include "fem/fe_rt.hh"
28 #include "fields/field_fe.hh"
29 #include "la/linsys_PETSC.hh"
30 #include "transport/advection_diffusion_model.hh"
31 #include "transport/concentration_model.hh"
32 #include "transport/heat_model.hh"
33 #include "coupling/balance.hh"
34 
35 #include "fields/multi_field.hh"
36 #include "fields/generic_field.hh"
37 #include "input/factory.hh"
38 #include "io/equation_output.hh"
39 
40 FLOW123D_FORCE_LINK_IN_CHILD(concentrationTransportModel);
41 FLOW123D_FORCE_LINK_IN_CHILD(heatModel);
42 
43 
44 
45 using namespace Input::Type;
46 
47 template<class Model>
48 const Selection & TransportDG<Model>::get_dg_variant_selection_input_type() {
49  return Selection("DG_variant", "Type of penalty term.")
50  .add_value(non_symmetric, "non-symmetric", "non-symmetric weighted interior penalty DG method")
51  .add_value(incomplete, "incomplete", "incomplete weighted interior penalty DG method")
52  .add_value(symmetric, "symmetric", "symmetric weighted interior penalty DG method")
53  .close();
54 }
55 
56 /*
57  * Should be removed
58 template<class Model>
59 const Selection & TransportDG<Model>::EqData::get_output_selection() {
60  return Model::ModelEqData::get_output_selection_input_type(
61  "DG",
62  "Implicit in time Discontinuous Galerkin solver")
63  .copy_values(EqData().make_output_field_selection("").close())
64  ConvectionTransport::EqData().output_fields
65  .make_output_field_selection(
66  "ConvectionTransport_output_fields",
67  "Selection of output fields for Convection Solute Transport model.")
68  .close()),
69  .close();
70 }
71 */
72 
73 template<class Model>
74 const Record & TransportDG<Model>::get_input_type() {
75  std::string equation_name = std::string(Model::ModelEqData::name()) + "_DG";
76  return Model::get_input_type("DG", "DG solver")
77  .declare_key("solver", LinSys_PETSC::get_input_type(), Default::obligatory(),
78  "Linear solver for MH problem.")
79  .declare_key("input_fields", Array(
80  TransportDG<Model>::EqData()
81  .make_field_descriptor_type(equation_name)),
82  IT::Default::obligatory(),
83  "Input fields of the equation.")
84  .declare_key("dg_variant", TransportDG<Model>::get_dg_variant_selection_input_type(), Default("\"non-symmetric\""),
85  "Variant of interior penalty discontinuous Galerkin method.")
86  .declare_key("dg_order", Integer(0,3), Default("1"),
87  "Polynomial order for finite element in DG method (order 0 is suitable if there is no diffusion/dispersion).")
88 /*
89  .declare_key("output_fields",
90  Array(
91  // Get selection name and description from the model
92  Model::ModelEqData::get_output_selection()
93  // EqData contains both TransportDG and model specific fields.
94  .copy_values(
95  EqData().make_output_field_selection("DG_output_fields","Auxiliary Selection")
96  .close())
97  .close()),
98  Default(Model::ModelEqData::default_output_field()),
99  "List of fields to write to output file.")
100 */
101  .declare_key("output",
102  EqData().output_fields.make_output_type(equation_name, ""),
103  IT::Default("{ \"fields\": [ " + Model::ModelEqData::default_output_field() + "] }"),
104  "Setting of the field output.")
105  .close();
106 }
107 
108 template<class Model>
109 const int TransportDG<Model>::registrar =
110  Input::register_class< TransportDG<Model>, Mesh &, const Input::Record>(std::string(Model::ModelEqData::name()) + "_DG") +
111  TransportDG<Model>::get_input_type().size();
112 
113 
114 
115 
116 FEObjects::FEObjects(Mesh *mesh_, unsigned int fe_order)
117 {
118  unsigned int q_order;
119 
120  switch (fe_order)
121  {
122  case 0:
123  q_order = 0;
124  fe1_ = new FE_P_disc<0,1,3>;
125  fe2_ = new FE_P_disc<0,2,3>;
126  fe3_ = new FE_P_disc<0,3,3>;
127  break;
128 
129  case 1:
130  q_order = 2;
131  fe1_ = new FE_P_disc<1,1,3>;
132  fe2_ = new FE_P_disc<1,2,3>;
133  fe3_ = new FE_P_disc<1,3,3>;
134  break;
135 
136  case 2:
137  q_order = 4;
138  fe1_ = new FE_P_disc<2,1,3>;
139  fe2_ = new FE_P_disc<2,2,3>;
140  fe3_ = new FE_P_disc<2,3,3>;
141  break;
142 
143  case 3:
144  q_order = 6;
145  fe1_ = new FE_P_disc<3,1,3>;
146  fe2_ = new FE_P_disc<3,2,3>;
147  fe3_ = new FE_P_disc<3,3,3>;
148  break;
149 
150  default:
151  q_order=0;
152  xprintf(PrgErr, "Unsupported polynomial order %d for finite elements in TransportDG ", fe_order);
153  break;
154  }
155 
156  fe_rt1_ = new FE_RT0<1,3>;
157  fe_rt2_ = new FE_RT0<2,3>;
158  fe_rt3_ = new FE_RT0<3,3>;
159 
160  q0_ = new QGauss<0>(q_order);
161  q1_ = new QGauss<1>(q_order);
162  q2_ = new QGauss<2>(q_order);
163  q3_ = new QGauss<3>(q_order);
164 
165  map0_ = new MappingP1<0,3>;
166  map1_ = new MappingP1<1,3>;
167  map2_ = new MappingP1<2,3>;
168  map3_ = new MappingP1<3,3>;
169 
170  dh_ = new DOFHandlerMultiDim(*mesh_);
171 
172  dh_->distribute_dofs(*fe1_, *fe2_, *fe3_);
173 }
174 
175 
176 FEObjects::~FEObjects()
177 {
178  delete fe1_;
179  delete fe2_;
180  delete fe3_;
181  delete fe_rt1_;
182  delete fe_rt2_;
183  delete fe_rt3_;
184  delete q0_;
185  delete q1_;
186  delete q2_;
187  delete q3_;
188  delete map0_;
189  delete map1_;
190  delete map2_;
191  delete map3_;
192  delete dh_;
193 }
194 
195 template<> FiniteElement<0,3> *FEObjects::fe<0>() { return 0; }
196 template<> FiniteElement<1,3> *FEObjects::fe<1>() { return fe1_; }
197 template<> FiniteElement<2,3> *FEObjects::fe<2>() { return fe2_; }
198 template<> FiniteElement<3,3> *FEObjects::fe<3>() { return fe3_; }
199 
200 template<> FiniteElement<0,3> *FEObjects::fe_rt<0>() { return 0; }
201 template<> FiniteElement<1,3> *FEObjects::fe_rt<1>() { return fe_rt1_; }
202 template<> FiniteElement<2,3> *FEObjects::fe_rt<2>() { return fe_rt2_; }
203 template<> FiniteElement<3,3> *FEObjects::fe_rt<3>() { return fe_rt3_; }
204 
205 template<> Quadrature<0> *FEObjects::q<0>() { return q0_; }
206 template<> Quadrature<1> *FEObjects::q<1>() { return q1_; }
207 template<> Quadrature<2> *FEObjects::q<2>() { return q2_; }
208 template<> Quadrature<3> *FEObjects::q<3>() { return q3_; }
209 
210 template<> Mapping<0,3> *FEObjects::mapping<0>() { return map0_; }
211 template<> Mapping<1,3> *FEObjects::mapping<1>() { return map1_; }
212 template<> Mapping<2,3> *FEObjects::mapping<2>() { return map2_; }
213 template<> Mapping<3,3> *FEObjects::mapping<3>() { return map3_; }
214 
215 DOFHandlerMultiDim *FEObjects::dh() { return dh_; }
216 
217 
218 template<class Model>
219 TransportDG<Model>::EqData::EqData() : Model::ModelEqData()
220 {
221  *this+=fracture_sigma
222  .name("fracture_sigma")
223  .description(
224  "Coefficient of diffusive transfer through fractures (for each substance).")
225  .units( UnitSI::dimensionless() )
226  .input_default("1.0")
227  .flags_add(FieldFlag::in_main_matrix);
228 
229  *this+=dg_penalty
230  .name("dg_penalty")
231  .description(
232  "Penalty parameter influencing the discontinuity of the solution (for each substance). "
233  "Its default value 1 is sufficient in most cases. Higher value diminishes the inter-element jumps.")
234  .units( UnitSI::dimensionless() )
235  .input_default("1.0")
236  .flags_add(FieldFlag::in_rhs & FieldFlag::in_main_matrix);
237 
238  *this += region_id.name("region_id")
239  .units( UnitSI::dimensionless())
240  .flags(FieldFlag::equation_external_output);
241 
242 
243  // add all input fields to the output list
244  output_fields += *this;
245 
246 }
247 
248 template<class Model>
249 TransportDG<Model>::TransportDG(Mesh & init_mesh, const Input::Record in_rec)
250  : Model(init_mesh, in_rec),
251  input_rec(in_rec),
252  allocation_done(false)
253 {
254  // Can not use name() + "constructor" here, since START_TIMER only accepts const char *
255  // due to constexpr optimization.
256  START_TIMER(Model::ModelEqData::name());
257  // Check that Model is derived from AdvectionDiffusionModel.
258  static_assert(std::is_base_of<AdvectionDiffusionModel, Model>::value, "");
259 
260  this->eq_data_ = &data_;
261 
262 
263  // Set up physical parameters.
264  data_.set_mesh(init_mesh);
265  data_.region_id = GenericField<3>::region_id(*Model::mesh_);
266 
267 
268  // DG variant and order
269  dg_variant = in_rec.val<DGVariant>("dg_variant");
270  dg_order = in_rec.val<unsigned int>("dg_order");
271 
272  Model::init_from_input(in_rec);
273 
274  // create finite element structures and distribute DOFs
275  feo = new FEObjects(Model::mesh_, dg_order);
276  //DebugOut().fmt("TDG: solution size {}\n", feo->dh()->n_global_dofs());
277 
278 }
279 
280 
281 template<class Model>
282 void TransportDG<Model>::initialize()
283 {
284  data_.set_components(Model::substances_.names());
285  data_.set_input_list( input_rec.val<Input::Array>("input_fields") );
286 
287  // DG stabilization parameters on boundary edges
288  gamma.resize(Model::n_substances());
289  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
290  gamma[sbi].resize(Model::mesh_->boundary_.size());
291 
292  // Resize coefficient arrays
293  int qsize = max(feo->q<0>()->size(), max(feo->q<1>()->size(), max(feo->q<2>()->size(), feo->q<3>()->size())));
294  int max_edg_sides = max(Model::mesh_->max_edge_sides(1), max(Model::mesh_->max_edge_sides(2), Model::mesh_->max_edge_sides(3)));
295  mm_coef.resize(qsize);
296  ret_coef.resize(Model::n_substances());
297  sorption_sources.resize(Model::n_substances());
298  ad_coef.resize(Model::n_substances());
299  dif_coef.resize(Model::n_substances());
300  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
301  {
302  ret_coef[sbi].resize(qsize);
303  ad_coef[sbi].resize(qsize);
304  dif_coef[sbi].resize(qsize);
305  }
306  ad_coef_edg.resize(max_edg_sides);
307  dif_coef_edg.resize(max_edg_sides);
308  for (int sd=0; sd<max_edg_sides; sd++)
309  {
310  ad_coef_edg[sd].resize(Model::n_substances());
311  dif_coef_edg[sd].resize(Model::n_substances());
312  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
313  {
314  ad_coef_edg[sd][sbi].resize(qsize);
315  dif_coef_edg[sd][sbi].resize(qsize);
316  }
317  }
318 
319  output_vec.resize(Model::n_substances());
320  //output_solution.resize(Model::n_substances());
321  int rank;
322  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
323  unsigned int output_vector_size= (rank==0)?feo->dh()->n_global_dofs():0;
324  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
325  {
326  // for each substance we allocate output array and vector
327  //output_solution[sbi] = new double[feo->dh()->n_global_dofs()];
328  VecCreateSeq(PETSC_COMM_SELF, output_vector_size, &output_vec[sbi]);
329  }
330  data_.output_field.set_components(Model::substances_.names());
331  data_.output_field.set_mesh(*Model::mesh_);
332  data_.output_type(OutputTime::CORNER_DATA);
333 
334  data_.output_field.setup_components();
335  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
336  {
337  // create shared pointer to a FieldFE, pass FE data and push this FieldFE to output_field on all regions
338  std::shared_ptr<FieldFE<3, FieldValue<3>::Scalar> > output_field_ptr(new FieldFE<3, FieldValue<3>::Scalar>);
339  output_field_ptr->set_fe_data(feo->dh(), feo->mapping<1>(), feo->mapping<2>(), feo->mapping<3>(), &output_vec[sbi]);
340  data_.output_field[sbi].set_field(Model::mesh_->region_db().get_region_set("ALL"), output_field_ptr, 0);
341  }
342 
343  // set time marks for writing the output
344  data_.output_fields.initialize(Model::output_stream_, input_rec.val<Input::Record>("output"), this->time());
345 
346  // allocate matrix and vector structures
347  ls = new LinSys*[Model::n_substances()];
348  ls_dt = new LinSys*[Model::n_substances()];
349  solution_elem_ = new double*[Model::n_substances()];
350  for (unsigned int sbi = 0; sbi < Model::n_substances(); sbi++) {
351  ls[sbi] = new LinSys_PETSC(feo->dh()->distr());
352  ( (LinSys_PETSC *)ls[sbi] )->set_from_input( input_rec.val<Input::Record>("solver") );
353  ls[sbi]->set_solution(NULL);
354 
355  ls_dt[sbi] = new LinSys_PETSC(feo->dh()->distr());
356  ( (LinSys_PETSC *)ls_dt[sbi] )->set_from_input( input_rec.val<Input::Record>("solver") );
357  solution_elem_[sbi] = new double[Model::mesh_->get_el_ds()->lsize()];
358  }
359  stiffness_matrix = new Mat[Model::n_substances()];
360  mass_matrix = new Mat[Model::n_substances()];
361  rhs = new Vec[Model::n_substances()];
362  mass_vec = new Vec[Model::n_substances()];
363 
364 
365  // initialization of balance object
366  if (Model::balance_ != nullptr)
367  {
368  Model::balance_->allocate(feo->dh()->distr()->lsize(),
369  max(feo->fe<1>()->n_dofs(), max(feo->fe<2>()->n_dofs(), feo->fe<3>()->n_dofs())));
370  }
371 
372 }
373 
374 
375 template<class Model>
376 TransportDG<Model>::~TransportDG()
377 {
378  delete Model::time_;
379 
380  if (gamma.size() > 0) {
381  // initialize called
382 
383  for (auto &vec : output_vec) VecDestroy(&vec);
384 
385  for (unsigned int i=0; i<Model::n_substances(); i++)
386  {
387  delete ls[i];
388  delete[] solution_elem_[i];
389  delete ls_dt[i];
390  MatDestroy(&stiffness_matrix[i]);
391  MatDestroy(&mass_matrix[i]);
392  VecDestroy(&rhs[i]);
393  VecDestroy(&mass_vec[i]);
394  }
395  delete[] ls;
396  delete[] solution_elem_;
397  delete[] ls_dt;
398  delete[] stiffness_matrix;
399  delete[] mass_matrix;
400  delete[] rhs;
401  delete[] mass_vec;
402  delete feo;
403 
404  }
405 
406 }
407 
408 
409 template<class Model>
410 void TransportDG<Model>::output_vector_gather()
411 {
412  VecScatter output_scatter;
413  VecScatterCreateToZero(ls[0]->get_solution(), &output_scatter, PETSC_NULL);
414  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
415  {
416  // gather solution to output_vec[sbi]
417  VecScatterBegin(output_scatter, ls[sbi]->get_solution(), output_vec[sbi], INSERT_VALUES, SCATTER_FORWARD);
418  VecScatterEnd(output_scatter, ls[sbi]->get_solution(), output_vec[sbi], INSERT_VALUES, SCATTER_FORWARD);
419  }
420  VecScatterDestroy(&(output_scatter));
421 
422 }
423 
424 
425 
426 template<class Model>
427 void TransportDG<Model>::zero_time_step()
428 {
429  START_TIMER(Model::ModelEqData::name());
430  data_.mark_input_times( *(Model::time_) );
431  data_.set_time(Model::time_->step(), LimitSide::left);
432 
433 
434  // set initial conditions
435  set_initial_condition();
436  for (unsigned int sbi = 0; sbi < Model::n_substances(); sbi++)
437  ( (LinSys_PETSC *)ls[sbi] )->set_initial_guess_nonzero();
438 
439  // check first time assembly - needs preallocation
440  if (!allocation_done) preallocate();
441 
442  // after preallocation we assemble the matrices and vectors required for mass balance
443  if (Model::balance_ != nullptr)
444  {
445  for (unsigned int sbi=0; sbi<Model::n_substances(); ++sbi)
446  {
447  Model::balance_->calculate_mass(Model::subst_idx[sbi], ls[sbi]->get_solution());
448  Model::balance_->calculate_source(Model::subst_idx[sbi], ls[sbi]->get_solution());
449  Model::balance_->calculate_flux(Model::subst_idx[sbi], ls[sbi]->get_solution());
450 
451  // add sources due to sorption
452  vector<double> masses(Model::mesh_->region_db().bulk_size());
453  double mass = 0;
454  Model::balance_->calculate_mass(Model::subst_idx[sbi], ls[sbi]->get_solution(), masses);
455  for (auto reg_mass : masses)
456  mass += reg_mass;
457 
458  sorption_sources[sbi] = mass;
459  }
460  }
461 
462  output_data();
463 }
464 
465 
466 template<class Model>
467 void TransportDG<Model>::preallocate()
468 {
469  // preallocate system matrix
470  for (unsigned int i=0; i<Model::n_substances(); i++)
471  {
472  // preallocate system matrix
473  ls[i]->start_allocation();
474  stiffness_matrix[i] = NULL;
475  rhs[i] = NULL;
476 
477  // preallocate mass matrix
478  ls_dt[i]->start_allocation();
479  mass_matrix[i] = NULL;
480  }
481  assemble_stiffness_matrix();
482  assemble_mass_matrix();
483  set_sources();
484  set_boundary_conditions();
485 
486  allocation_done = true;
487 }
488 
489 
490 
491 template<class Model>
492 void TransportDG<Model>::update_solution()
493 {
494  START_TIMER("DG-ONE STEP");
495 
496  Model::time_->next_time();
497  Model::time_->view("TDG");
498 
499  START_TIMER("data reinit");
500  data_.set_time(Model::time_->step(), LimitSide::left);
501  END_TIMER("data reinit");
502 
503  // assemble mass matrix
504  if (mass_matrix[0] == NULL || data_.subset(FieldFlag::in_time_term).changed() )
505  {
506  for (unsigned int i=0; i<Model::n_substances(); i++)
507  {
508  ls_dt[i]->start_add_assembly();
509  ls_dt[i]->mat_zero_entries();
510  }
511  assemble_mass_matrix();
512  for (unsigned int i=0; i<Model::n_substances(); i++)
513  {
514  ls_dt[i]->finish_assembly();
515  // construct mass_vec for initial time
516  if (mass_matrix[i] == NULL)
517  {
518  VecDuplicate(ls[i]->get_solution(), &mass_vec[i]);
519  MatMult(*(ls_dt[i]->get_matrix()), ls[i]->get_solution(), mass_vec[i]);
520  if (Model::balance_ != nullptr && Model::balance_->cumulative())
521  {
522  double total_mass = 0;
523  VecSum(mass_vec[i], &total_mass);
524  sorption_sources[i] -= total_mass;
525  }
526  MatConvert(*( ls_dt[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &mass_matrix[i]);
527  }
528  else
529  MatCopy(*( ls_dt[i]->get_matrix() ), mass_matrix[i], DIFFERENT_NONZERO_PATTERN);
530  }
531  }
532 
533  // assemble stiffness matrix
534  if (stiffness_matrix[0] == NULL
535  || data_.subset(FieldFlag::in_main_matrix).changed()
536  || Model::flux_changed)
537  {
538  // new fluxes can change the location of Neumann boundary,
539  // thus stiffness matrix must be reassembled
540  for (unsigned int i=0; i<Model::n_substances(); i++)
541  {
542  ls[i]->start_add_assembly();
543  ls[i]->mat_zero_entries();
544  }
545  assemble_stiffness_matrix();
546  for (unsigned int i=0; i<Model::n_substances(); i++)
547  {
548  ls[i]->finish_assembly();
549 
550  if (stiffness_matrix[i] == NULL)
551  MatConvert(*( ls[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &stiffness_matrix[i]);
552  else
553  MatCopy(*( ls[i]->get_matrix() ), stiffness_matrix[i], DIFFERENT_NONZERO_PATTERN);
554  }
555  }
556 
557  // assemble right hand side (due to sources and boundary conditions)
558  if (rhs[0] == NULL
559  || data_.subset(FieldFlag::in_rhs).changed()
560  || Model::flux_changed)
561  {
562  for (unsigned int i=0; i<Model::n_substances(); i++)
563  {
564  ls[i]->start_add_assembly();
565  ls[i]->rhs_zero_entries();
566  }
567  set_sources();
568  set_boundary_conditions();
569  for (unsigned int i=0; i<Model::n_substances(); i++)
570  {
571  ls[i]->finish_assembly();
572 
573  if (rhs[i] == nullptr) VecDuplicate(*( ls[i]->get_rhs() ), &rhs[i]);
574  VecCopy(*( ls[i]->get_rhs() ), rhs[i]);
575  }
576  }
577 
578  Model::flux_changed = false;
579 
580 
581  /* Apply backward Euler time integration.
582  *
583  * Denoting A the stiffness matrix and M the mass matrix, the algebraic system at the k-th time level reads
584  *
585  * (1/dt M + A)u^k = f + 1/dt M.u^{k-1}
586  *
587  * Hence we modify at each time level the right hand side:
588  *
589  * f^k = f + 1/dt M u^{k-1},
590  *
591  * where f stands for the term stemming from the force and boundary conditions.
592  * Accordingly, we set
593  *
594  * A^k = A + 1/dt M.
595  *
596  */
597  Mat m;
598  START_TIMER("solve");
599  for (unsigned int i=0; i<Model::n_substances(); i++)
600  {
601  MatConvert(stiffness_matrix[i], MATSAME, MAT_INITIAL_MATRIX, &m);
602  MatAXPY(m, 1./Model::time_->dt(), mass_matrix[i], SUBSET_NONZERO_PATTERN);
603  ls[i]->set_matrix(m, DIFFERENT_NONZERO_PATTERN);
604  Vec w;
605  VecDuplicate(rhs[i], &w);
606  VecWAXPY(w, 1./Model::time_->dt(), mass_vec[i], rhs[i]);
607  ls[i]->set_rhs(w);
608 
609  VecDestroy(&w);
610  MatDestroy(&m);
611 
612  ls[i]->solve();
613 
614  // update mass_vec due to possible changes in mass matrix
615  MatMult(*(ls_dt[i]->get_matrix()), ls[i]->get_solution(), mass_vec[i]);
616  }
617  END_TIMER("solve");
618 
619  calculate_cumulative_balance();
620 
621  END_TIMER("DG-ONE STEP");
622 }
623 
624 
625 template<class Model>
626 void TransportDG<Model>::calculate_concentration_matrix()
627 {
628  // calculate element averages of solution
629  for (unsigned int i_cell=0; i_cell<Model::mesh_->get_el_ds()->lsize(); i_cell++)
630  {
631  typename DOFHandlerBase::CellIterator elem = Model::mesh_->element(feo->dh()->el_index(i_cell));
632 
633  unsigned int n_dofs;
634  switch (elem->dim())
635  {
636  case 1:
637  n_dofs = feo->fe<1>()->n_dofs();
638  break;
639  case 2:
640  n_dofs = feo->fe<2>()->n_dofs();
641  break;
642  case 3:
643  n_dofs = feo->fe<3>()->n_dofs();
644  break;
645  }
646 
647  unsigned int dof_indices[n_dofs];
648  feo->dh()->get_dof_indices(elem, dof_indices);
649 
650  for (unsigned int sbi=0; sbi<Model::n_substances(); ++sbi)
651  {
652  solution_elem_[sbi][i_cell] = 0;
653 
654  for (unsigned int j=0; j<n_dofs; ++j)
655  solution_elem_[sbi][i_cell] += ls[sbi]->get_solution_array()[dof_indices[j]-feo->dh()->distr()->begin()];
656 
657  solution_elem_[sbi][i_cell] = max(solution_elem_[sbi][i_cell]/n_dofs, 0.);
658  }
659  }
660 }
661 
662 
663 
664 
665 template<class Model>
666 void TransportDG<Model>::output_data()
667 {
668  //if (!Model::time_->is_current( Model::time_->marks().type_output() )) return;
669 
670 
671  START_TIMER("DG-OUTPUT");
672 
673  // gather the solution from all processors
674  data_.output_fields.set_time( this->time().step(), LimitSide::left);
675  //if (data_.output_fields.is_field_output_time(data_.output_field, this->time().step()) )
676  output_vector_gather();
677  data_.output_fields.output(this->time().step());
678 
679  Model::output_data();
680 
681  END_TIMER("DG-OUTPUT");
682 }
683 
684 
685 template<class Model>
686 void TransportDG<Model>::calculate_cumulative_balance()
687 {
688  if (Model::balance_ != nullptr && Model::balance_->cumulative())
689  {
690  for (unsigned int sbi=0; sbi<Model::n_substances(); ++sbi)
691  {
692  Model::balance_->calculate_cumulative_sources(Model::subst_idx[sbi], ls[sbi]->get_solution(), Model::time_->dt());
693  Model::balance_->calculate_cumulative_fluxes(Model::subst_idx[sbi], ls[sbi]->get_solution(), Model::time_->dt());
694 
695  // add sources due to sorption
696  vector<double> masses(Model::mesh_->region_db().bulk_size());
697  double mass = 0;
698  Model::balance_->calculate_mass(Model::subst_idx[sbi], ls[sbi]->get_solution(), masses);
699  for (auto reg_mass : masses)
700  mass += reg_mass;
701  double total_mass = 0;
702  VecSum(mass_vec[sbi], &total_mass);
703 
704  Model::balance_->add_cumulative_source(Model::subst_idx[sbi], mass-total_mass-sorption_sources[sbi]);
705  sorption_sources[sbi] = mass - total_mass;
706  }
707  }
708 }
709 
710 
711 template<class Model>
712 void TransportDG<Model>::calculate_instant_balance()
713 {
714  if (Model::balance_ != nullptr)
715  {
716  for (unsigned int sbi=0; sbi<Model::n_substances(); ++sbi)
717  {
718  Model::balance_->calculate_mass(Model::subst_idx[sbi], ls[sbi]->get_solution());
719  Model::balance_->calculate_source(Model::subst_idx[sbi], ls[sbi]->get_solution());
720  Model::balance_->calculate_flux(Model::subst_idx[sbi], ls[sbi]->get_solution());
721  }
722  }
723 }
724 
725 
726 template<class Model>
727 void TransportDG<Model>::assemble_mass_matrix()
728 {
729  START_TIMER("assemble_mass");
730  if (Model::balance_ != nullptr)
731  Model::balance_->start_mass_assembly(Model::subst_idx);
732  assemble_mass_matrix<1>();
733  assemble_mass_matrix<2>();
734  assemble_mass_matrix<3>();
735  if (Model::balance_ != nullptr)
736  Model::balance_->finish_mass_assembly(Model::subst_idx);
737  END_TIMER("assemble_mass");
738 }
739 
740 
741 template<class Model> template<unsigned int dim>
742 void TransportDG<Model>::assemble_mass_matrix()
743 {
744  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(), update_values | update_JxW_values | update_quadrature_points);
745  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
746  vector<int> dof_indices(ndofs);
747  PetscScalar local_mass_matrix[ndofs*ndofs];
748  vector<PetscScalar> local_mass_balance_vector(ndofs);
749 
750  // assemble integral over elements
751  for (unsigned int i_cell=0; i_cell<Model::mesh_->get_el_ds()->lsize(); i_cell++)
752  {
753  typename DOFHandlerBase::CellIterator cell = Model::mesh_->element(feo->dh()->el_index(i_cell));
754  if (cell->dim() != dim) continue;
755 
756  fe_values.reinit(cell);
757  feo->dh()->get_dof_indices(cell, (unsigned int*)&(dof_indices[0]));
758  ElementAccessor<3> ele_acc = cell->element_accessor();
759 
760  Model::compute_mass_matrix_coefficient(fe_values.point_list(), ele_acc, mm_coef);
761  Model::compute_retardation_coefficient(fe_values.point_list(), ele_acc, ret_coef);
762 
763  for (unsigned int sbi=0; sbi<Model::n_substances(); ++sbi)
764  {
765  // assemble the local mass matrix
766  for (unsigned int i=0; i<ndofs; i++)
767  {
768  for (unsigned int j=0; j<ndofs; j++)
769  {
770  local_mass_matrix[i*ndofs+j] = 0;
771  for (unsigned int k=0; k<qsize; k++)
772  local_mass_matrix[i*ndofs+j] += (mm_coef[k]+ret_coef[sbi][k])*fe_values.shape_value(j,k)*fe_values.shape_value(i,k)*fe_values.JxW(k);
773  }
774  }
775 
776  if (Model::balance_ != nullptr)
777  {
778  for (unsigned int i=0; i<ndofs; i++)
779  {
780  local_mass_balance_vector[i] = 0;
781  for (unsigned int k=0; k<qsize; k++)
782  local_mass_balance_vector[i] += mm_coef[k]*fe_values.shape_value(i,k)*fe_values.JxW(k);
783  }
784  }
785 
786  Model::balance_->add_mass_matrix_values(Model::subst_idx[sbi], ele_acc.region().bulk_idx(), dof_indices, local_mass_balance_vector);
787  ls_dt[sbi]->mat_set_values(ndofs, &(dof_indices[0]), ndofs, &(dof_indices[0]), local_mass_matrix);
788  }
789  }
790 }
791 
792 
793 
794 
795 template<class Model>
796 void TransportDG<Model>::assemble_stiffness_matrix()
797 {
798  START_TIMER("assemble_stiffness");
799  START_TIMER("assemble_volume_integrals");
800  assemble_volume_integrals<1>();
801  assemble_volume_integrals<2>();
802  assemble_volume_integrals<3>();
803  END_TIMER("assemble_volume_integrals");
804 
805  START_TIMER("assemble_fluxes_boundary");
806  assemble_fluxes_boundary<1>();
807  assemble_fluxes_boundary<2>();
808  assemble_fluxes_boundary<3>();
809  END_TIMER("assemble_fluxes_boundary");
810 
811  START_TIMER("assemble_fluxes_elem_elem");
812  assemble_fluxes_element_element<1>();
813  assemble_fluxes_element_element<2>();
814  assemble_fluxes_element_element<3>();
815  END_TIMER("assemble_fluxes_elem_elem");
816 
817  START_TIMER("assemble_fluxes_elem_side");
818  assemble_fluxes_element_side<1>();
819  assemble_fluxes_element_side<2>();
820  assemble_fluxes_element_side<3>();
821  END_TIMER("assemble_fluxes_elem_side");
822  END_TIMER("assemble_stiffness");
823 }
824 
825 
826 
827 template<class Model>
828 template<unsigned int dim>
829 void TransportDG<Model>::assemble_volume_integrals()
830 {
831  FEValues<dim,3> fv_rt(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe_rt<dim>(),
832  update_values | update_gradients);
833  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
834  update_values | update_gradients | update_JxW_values | update_quadrature_points);
835  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
836  unsigned int dof_indices[ndofs];
837  vector<arma::vec3> velocity(qsize);
838  vector<arma::vec> sources_sigma(qsize, arma::vec(Model::n_substances()));
839  PetscScalar local_matrix[ndofs*ndofs];
840 
841  // assemble integral over elements
842  for (unsigned int i_cell=0; i_cell<Model::mesh_->get_el_ds()->lsize(); i_cell++)
843  {
844  typename DOFHandlerBase::CellIterator cell = Model::mesh_->element(feo->dh()->el_index(i_cell));
845  if (cell->dim() != dim) continue;
846 
847  fe_values.reinit(cell);
848  fv_rt.reinit(cell);
849  ElementAccessor<3> ele_acc = cell->element_accessor();
850  feo->dh()->get_dof_indices(cell, dof_indices);
851 
852  calculate_velocity(cell, velocity, fv_rt);
853  Model::compute_advection_diffusion_coefficients(fe_values.point_list(), velocity, ele_acc, ad_coef, dif_coef);
854  Model::compute_sources_sigma(fe_values.point_list(), ele_acc, sources_sigma);
855 
856  // assemble the local stiffness matrix
857  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
858  {
859  for (unsigned int i=0; i<ndofs; i++)
860  for (unsigned int j=0; j<ndofs; j++)
861  local_matrix[i*ndofs+j] = 0;
862 
863  for (unsigned int k=0; k<qsize; k++)
864  {
865  for (unsigned int i=0; i<ndofs; i++)
866  {
867  arma::vec3 Kt_grad_i = dif_coef[sbi][k].t()*fe_values.shape_grad(i,k);
868  double ad_dot_grad_i = arma::dot(ad_coef[sbi][k], fe_values.shape_grad(i,k));
869 
870  for (unsigned int j=0; j<ndofs; j++)
871  local_matrix[i*ndofs+j] += (arma::dot(Kt_grad_i, fe_values.shape_grad(j,k))
872  -fe_values.shape_value(j,k)*ad_dot_grad_i
873  +sources_sigma[k][sbi]*fe_values.shape_value(j,k)*fe_values.shape_value(i,k))*fe_values.JxW(k);
874  }
875  }
876  ls[sbi]->mat_set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, local_matrix);
877  }
878  }
879 }
880 
881 
882 template<class Model>
883 void TransportDG<Model>::set_sources()
884 {
885  START_TIMER("assemble_sources");
886  if (Model::balance_ != nullptr)
887  Model::balance_->start_source_assembly(Model::subst_idx);
888  set_sources<1>();
889  set_sources<2>();
890  set_sources<3>();
891  if (Model::balance_ != nullptr)
892  Model::balance_->finish_source_assembly(Model::subst_idx);
893  END_TIMER("assemble_sources");
894 }
895 
896 template<class Model>
897 template<unsigned int dim>
898 void TransportDG<Model>::set_sources()
899 {
900  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
901  update_values | update_JxW_values | update_quadrature_points);
902  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
903  vector<arma::vec> sources_conc(qsize, arma::vec(Model::n_substances())),
904  sources_density(qsize, arma::vec(Model::n_substances())),
905  sources_sigma(qsize, arma::vec(Model::n_substances()));
906  vector<int> dof_indices(ndofs);
907  PetscScalar local_rhs[ndofs];
908  vector<PetscScalar> local_source_balance_vector(ndofs), local_source_balance_rhs(ndofs);
909  double source;
910 
911  // assemble integral over elements
912  for (unsigned int i_cell=0; i_cell<Model::mesh_->get_el_ds()->lsize(); i_cell++)
913  {
914  typename DOFHandlerBase::CellIterator cell = Model::mesh_->element(feo->dh()->el_index(i_cell));
915  if (cell->dim() != dim) continue;
916 
917  fe_values.reinit(cell);
918  feo->dh()->get_dof_indices(cell, (unsigned int *)&(dof_indices[0]));
919 
920  Model::compute_source_coefficients(fe_values.point_list(), cell->element_accessor(), sources_conc, sources_density, sources_sigma);
921 
922  // assemble the local stiffness matrix
923  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
924  {
925  fill_n(local_rhs, ndofs, 0);
926  local_source_balance_vector.assign(ndofs, 0);
927  local_source_balance_rhs.assign(ndofs, 0);
928 
929  // compute sources
930  for (unsigned int k=0; k<qsize; k++)
931  {
932  source = (sources_density[k][sbi] + sources_conc[k][sbi]*sources_sigma[k][sbi])*fe_values.JxW(k);
933 
934  for (unsigned int i=0; i<ndofs; i++)
935  local_rhs[i] += source*fe_values.shape_value(i,k);
936  }
937  ls[sbi]->rhs_set_values(ndofs, &(dof_indices[0]), local_rhs);
938 
939  if (Model::balance_ != nullptr)
940  {
941  for (unsigned int i=0; i<ndofs; i++)
942  {
943  for (unsigned int k=0; k<qsize; k++)
944  local_source_balance_vector[i] -= sources_sigma[k][sbi]*fe_values.shape_value(i,k)*fe_values.JxW(k);
945 
946  local_source_balance_rhs[i] += local_rhs[i];
947  }
948  Model::balance_->add_source_matrix_values(Model::subst_idx[sbi], cell->region().bulk_idx(), dof_indices, local_source_balance_vector);
949  Model::balance_->add_source_vec_values(Model::subst_idx[sbi], cell->region().bulk_idx(), dof_indices, local_source_balance_rhs);
950  }
951  }
952  }
953 }
954 
955 
956 
957 template<class Model>
958 template<unsigned int dim>
959 void TransportDG<Model>::assemble_fluxes_element_element()
960 {
961  vector<FESideValues<dim,3>*> fe_values;
962  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
963  update_values);
964  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size(),
965  n_max_sides = ad_coef_edg.size();
966  vector<unsigned int*> side_dof_indices;
967  PetscScalar local_matrix[ndofs*ndofs];
968  vector<vector<arma::vec3> > side_velocity(n_max_sides);
969  vector<arma::vec> dg_penalty(n_max_sides);
970  double gamma_l, omega[2], transport_flux;
971 
972  for (unsigned int sid=0; sid<n_max_sides; sid++)
973  {
974  side_dof_indices.push_back(new unsigned int[ndofs]);
975  fe_values.push_back(new FESideValues<dim,3>(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
976  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points));
977  }
978 
979  // assemble integral over sides
980  for (unsigned int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
981  {
982  Edge *edg = &Model::mesh_->edges[feo->dh()->edge_index(iedg)];
983  if (edg->n_sides < 2 || edg->side(0)->element()->dim() != dim) continue;
984 
985  for (int sid=0; sid<edg->n_sides; sid++)
986  {
987  typename DOFHandlerBase::CellIterator cell = edg->side(sid)->element();
988  ElementAccessor<3> ele_acc = cell->element_accessor();
989  feo->dh()->get_dof_indices(cell, side_dof_indices[sid]);
990  fe_values[sid]->reinit(cell, edg->side(sid)->el_idx());
991  fsv_rt.reinit(cell, edg->side(sid)->el_idx());
992  calculate_velocity(cell, side_velocity[sid], fsv_rt);
993  Model::compute_advection_diffusion_coefficients(fe_values[sid]->point_list(), side_velocity[sid], ele_acc, ad_coef_edg[sid], dif_coef_edg[sid]);
994  dg_penalty[sid] = data_.dg_penalty.value(cell->centre(), ele_acc);
995  }
996 
997  // fluxes and penalty
998  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
999  {
1000  vector<double> fluxes(edg->n_sides);
1001  for (int sid=0; sid<edg->n_sides; sid++)
1002  {
1003  fluxes[sid] = 0;
1004  for (unsigned int k=0; k<qsize; k++)
1005  fluxes[sid] += arma::dot(ad_coef_edg[sid][sbi][k], fe_values[sid]->normal_vector(k))*fe_values[sid]->JxW(k);
1006  fluxes[sid] /= edg->side(sid)->measure();
1007  }
1008 
1009  for (int s1=0; s1<edg->n_sides; s1++)
1010  {
1011  for (int s2=s1+1; s2<edg->n_sides; s2++)
1012  {
1013  OLD_ASSERT(edg->side(s1)->valid(), "Invalid side of edge.");
1014  if (!feo->dh()->el_is_local(edg->side(s1)->element().index())
1015  && !feo->dh()->el_is_local(edg->side(s2)->element().index())) continue;
1016 
1017  arma::vec3 nv = fe_values[s1]->normal_vector(0);
1018 
1019  // set up the parameters for DG method
1020  set_DG_parameters_edge(*edg, s1, s2, qsize, dif_coef_edg[s1][sbi], dif_coef_edg[s2][sbi], fluxes, fe_values[0]->normal_vector(0), dg_penalty[s1][sbi], dg_penalty[s2][sbi], gamma_l, omega, transport_flux);
1021 
1022  int sd[2];
1023  sd[0] = s1;
1024  sd[1] = s2;
1025 
1026 #define AVERAGE(i,k,side_id) (fe_values[sd[side_id]]->shape_value(i,k)*0.5)
1027 #define WAVERAGE(i,k,side_id) (arma::dot(dif_coef_edg[sd[side_id]][sbi][k]*fe_values[sd[side_id]]->shape_grad(i,k),nv)*omega[side_id])
1028 #define JUMP(i,k,side_id) ((side_id==0?1:-1)*fe_values[sd[side_id]]->shape_value(i,k))
1029 
1030  // For selected pair of elements:
1031  for (int n=0; n<2; n++)
1032  {
1033  if (!feo->dh()->el_is_local(edg->side(sd[n])->element().index())) continue;
1034 
1035  for (int m=0; m<2; m++)
1036  {
1037  for (unsigned int i=0; i<fe_values[sd[n]]->n_dofs(); i++)
1038  for (unsigned int j=0; j<fe_values[sd[m]]->n_dofs(); j++)
1039  local_matrix[i*fe_values[sd[m]]->n_dofs()+j] = 0;
1040 
1041  for (unsigned int k=0; k<qsize; k++)
1042  {
1043  double flux_times_JxW = transport_flux*fe_values[0]->JxW(k);
1044  double gamma_times_JxW = gamma_l*fe_values[0]->JxW(k);
1045 
1046  for (unsigned int i=0; i<fe_values[sd[n]]->n_dofs(); i++)
1047  {
1048  double flux_JxW_jump_i = flux_times_JxW*JUMP(i,k,n);
1049  double gamma_JxW_jump_i = gamma_times_JxW*JUMP(i,k,n);
1050  double JxW_jump_i = fe_values[0]->JxW(k)*JUMP(i,k,n);
1051  double JxW_var_wavg_i = fe_values[0]->JxW(k)*WAVERAGE(i,k,n)*dg_variant;
1052 
1053  for (unsigned int j=0; j<fe_values[sd[m]]->n_dofs(); j++)
1054  {
1055  int index = i*fe_values[sd[m]]->n_dofs()+j;
1056 
1057  // flux due to transport (applied on interior edges) (average times jump)
1058  local_matrix[index] += flux_JxW_jump_i*AVERAGE(j,k,m);
1059 
1060  // penalty enforcing continuity across edges (applied on interior and Dirichlet edges) (jump times jump)
1061  local_matrix[index] += gamma_JxW_jump_i*JUMP(j,k,m);
1062 
1063  // terms due to diffusion
1064  local_matrix[index] -= WAVERAGE(j,k,m)*JxW_jump_i;
1065  local_matrix[index] -= JUMP(j,k,m)*JxW_var_wavg_i;
1066  }
1067  }
1068  }
1069  ls[sbi]->mat_set_values(fe_values[sd[n]]->n_dofs(), (int *)side_dof_indices[sd[n]], fe_values[sd[m]]->n_dofs(), (int *)side_dof_indices[sd[m]], local_matrix);
1070  }
1071  }
1072 #undef AVERAGE
1073 #undef WAVERAGE
1074 #undef JUMP
1075  }
1076  }
1077  }
1078  }
1079 
1080  for (unsigned int i=0; i<n_max_sides; i++)
1081  {
1082  delete fe_values[i];
1083  delete[] side_dof_indices[i];
1084  }
1085 }
1086 
1087 
1088 template<class Model>
1089 template<unsigned int dim>
1090 void TransportDG<Model>::assemble_fluxes_boundary()
1091 {
1092  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1093  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points);
1094  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1095  update_values);
1096  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
1097  unsigned int side_dof_indices[ndofs];
1098  PetscScalar local_matrix[ndofs*ndofs];
1099  vector<arma::vec3> side_velocity;
1100  vector<double> robin_sigma(qsize);
1101  vector<double> csection(qsize);
1102  arma::vec dg_penalty;
1103  double gamma_l;
1104 
1105  // assemble boundary integral
1106  for (unsigned int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
1107  {
1108  Edge *edg = &Model::mesh_->edges[feo->dh()->edge_index(iedg)];
1109  if (edg->n_sides > 1) continue;
1110  // check spatial dimension
1111  if (edg->side(0)->dim() != dim-1) continue;
1112  // skip edges lying not on the boundary
1113  if (edg->side(0)->cond() == NULL) continue;
1114 
1115  SideIter side = edg->side(0);
1116  typename DOFHandlerBase::CellIterator cell = side->element();
1117  ElementAccessor<3> ele_acc = cell->element_accessor();
1118  feo->dh()->get_dof_indices(cell, side_dof_indices);
1119  fe_values_side.reinit(cell, side->el_idx());
1120  fsv_rt.reinit(cell, side->el_idx());
1121 
1122  calculate_velocity(cell, side_velocity, fsv_rt);
1123  Model::compute_advection_diffusion_coefficients(fe_values_side.point_list(), side_velocity, ele_acc, ad_coef, dif_coef);
1124  dg_penalty = data_.dg_penalty.value(cell->centre(), ele_acc);
1125  arma::uvec bc_type;
1126  Model::get_bc_type(side->cond()->element_accessor(), bc_type);
1127  data_.cross_section.value_list(fe_values_side.point_list(), ele_acc, csection);
1128 
1129  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
1130  {
1131  for (unsigned int i=0; i<ndofs; i++)
1132  for (unsigned int j=0; j<ndofs; j++)
1133  local_matrix[i*ndofs+j] = 0;
1134 
1135  // On Neumann boundaries we have only term from integrating by parts the advective term,
1136  // on Dirichlet boundaries we additionally apply the penalty which enforces the prescribed value.
1137  double side_flux = 0;
1138  for (unsigned int k=0; k<qsize; k++)
1139  side_flux += arma::dot(ad_coef[sbi][k], fe_values_side.normal_vector(k))*fe_values_side.JxW(k);
1140  double transport_flux = side_flux/side->measure();
1141 
1142  if (bc_type[sbi] == AdvectionDiffusionModel::abc_dirichlet)
1143  {
1144  // set up the parameters for DG method
1145  set_DG_parameters_boundary(side, qsize, dif_coef[sbi], transport_flux, fe_values_side.normal_vector(0), dg_penalty[sbi], gamma_l);
1146  gamma[sbi][side->cond_idx()] = gamma_l;
1147  transport_flux += gamma_l;
1148  }
1149 
1150  // fluxes and penalty
1151  for (unsigned int k=0; k<qsize; k++)
1152  {
1153  double flux_times_JxW;
1154  if (bc_type[sbi] == AdvectionDiffusionModel::abc_total_flux)
1155  {
1156  Model::get_flux_bc_sigma(sbi, fe_values_side.point_list(), side->cond()->element_accessor(), robin_sigma);
1157  flux_times_JxW = csection[k]*robin_sigma[k]*fe_values_side.JxW(k);
1158  }
1159  else if (bc_type[sbi] == AdvectionDiffusionModel::abc_diffusive_flux)
1160  {
1161  Model::get_flux_bc_sigma(sbi, fe_values_side.point_list(), side->cond()->element_accessor(), robin_sigma);
1162  flux_times_JxW = (transport_flux + csection[k]*robin_sigma[k])*fe_values_side.JxW(k);
1163  }
1164  else if (bc_type[sbi] == AdvectionDiffusionModel::abc_inflow && side_flux < 0)
1165  flux_times_JxW = 0;
1166  else
1167  flux_times_JxW = transport_flux*fe_values_side.JxW(k);
1168 
1169  for (unsigned int i=0; i<ndofs; i++)
1170  {
1171  for (unsigned int j=0; j<ndofs; j++)
1172  {
1173  // flux due to advection and penalty
1174  local_matrix[i*ndofs+j] += flux_times_JxW*fe_values_side.shape_value(i,k)*fe_values_side.shape_value(j,k);
1175 
1176  // flux due to diffusion (only on dirichlet and inflow boundary)
1177  if (bc_type[sbi] == AdvectionDiffusionModel::abc_dirichlet)
1178  local_matrix[i*ndofs+j] -= (arma::dot(dif_coef[sbi][k]*fe_values_side.shape_grad(j,k),fe_values_side.normal_vector(k))*fe_values_side.shape_value(i,k)
1179  + arma::dot(dif_coef[sbi][k]*fe_values_side.shape_grad(i,k),fe_values_side.normal_vector(k))*fe_values_side.shape_value(j,k)*dg_variant
1180  )*fe_values_side.JxW(k);
1181  }
1182  }
1183  }
1184 
1185  ls[sbi]->mat_set_values(ndofs, (int *)side_dof_indices, ndofs, (int *)side_dof_indices, local_matrix);
1186  }
1187  }
1188 }
1189 
1190 
1191 template<class Model>
1192 template<unsigned int dim>
1193 void TransportDG<Model>::assemble_fluxes_element_side()
1194 {
1195 
1196  if (dim == 1) return;
1197  FEValues<dim-1,3> fe_values_vb(*feo->mapping<dim-1>(), *feo->q<dim-1>(), *feo->fe<dim-1>(),
1198  update_values | update_gradients | update_JxW_values | update_quadrature_points);
1199  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1200  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points);
1201  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1202  update_values);
1203  FEValues<dim-1,3> fv_rt(*feo->mapping<dim-1>(), *feo->q<dim-1>(), *feo->fe_rt<dim-1>(),
1204  update_values);
1205 
1206  vector<FEValuesSpaceBase<3>*> fv_sb(2);
1207  const unsigned int ndofs = feo->fe<dim>()->n_dofs(); // number of local dofs
1208  const unsigned int qsize = feo->q<dim-1>()->size(); // number of quadrature points
1209  unsigned int side_dof_indices[2*ndofs], n_dofs[2];
1210  vector<arma::vec3> velocity_higher, velocity_lower;
1211  vector<arma::vec> frac_sigma(qsize, arma::vec(Model::n_substances()));
1212  vector<double> csection_lower(qsize), csection_higher(qsize), por_lower(qsize), por_higher(qsize);
1213  PetscScalar local_matrix[4*ndofs*ndofs];
1214  double comm_flux[2][2];
1215 
1216  // index 0 = element with lower dimension,
1217  // index 1 = side of element with higher dimension
1218  fv_sb[0] = &fe_values_vb;
1219  fv_sb[1] = &fe_values_side;
1220 
1221  // assemble integral over sides
1222  for (unsigned int inb=0; inb<feo->dh()->n_loc_nb(); inb++)
1223  {
1224  Neighbour *nb = &Model::mesh_->vb_neighbours_[feo->dh()->nb_index(inb)];
1225  // skip neighbours of different dimension
1226  if (nb->element()->dim() != dim-1) continue;
1227 
1228  typename DOFHandlerBase::CellIterator cell_sub = Model::mesh_->element.full_iter(nb->element());
1229  feo->dh()->get_dof_indices(cell_sub, side_dof_indices);
1230  fe_values_vb.reinit(cell_sub);
1231  n_dofs[0] = fv_sb[0]->n_dofs();
1232 
1233  typename DOFHandlerBase::CellIterator cell = nb->side()->element();
1234  feo->dh()->get_dof_indices(cell, side_dof_indices+n_dofs[0]);
1235  fe_values_side.reinit(cell, nb->side()->el_idx());
1236  n_dofs[1] = fv_sb[1]->n_dofs();
1237 
1238  // Element id's for testing if they belong to local partition.
1239  int element_id[2];
1240  element_id[0] = cell_sub.index();
1241  element_id[1] = cell.index();
1242 
1243  fsv_rt.reinit(cell, nb->side()->el_idx());
1244  fv_rt.reinit(cell_sub);
1245  calculate_velocity(cell, velocity_higher, fsv_rt);
1246  calculate_velocity(cell_sub, velocity_lower, fv_rt);
1247  Model::compute_advection_diffusion_coefficients(fe_values_vb.point_list(), velocity_lower, cell_sub->element_accessor(), ad_coef_edg[0], dif_coef_edg[0]);
1248  Model::compute_advection_diffusion_coefficients(fe_values_vb.point_list(), velocity_higher, cell->element_accessor(), ad_coef_edg[1], dif_coef_edg[1]);
1249  data_.cross_section.value_list(fe_values_vb.point_list(), cell_sub->element_accessor(), csection_lower);
1250  data_.cross_section.value_list(fe_values_vb.point_list(), cell->element_accessor(), csection_higher);
1251  data_.fracture_sigma.value_list(fe_values_vb.point_list(), cell_sub->element_accessor(), frac_sigma);
1252  data_.porosity.value_list(fe_values_vb.point_list(), cell_sub->element_accessor(), por_lower);
1253  data_.porosity.value_list(fe_values_vb.point_list(), cell->element_accessor(), por_higher);
1254 
1255  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
1256  {
1257  for (unsigned int i=0; i<n_dofs[0]+n_dofs[1]; i++)
1258  for (unsigned int j=0; j<n_dofs[0]+n_dofs[1]; j++)
1259  local_matrix[i*(n_dofs[0]+n_dofs[1])+j] = 0;
1260 
1261  // set transmission conditions
1262  for (unsigned int k=0; k<qsize; k++)
1263  {
1264  /* The communication flux has two parts:
1265  * - "diffusive" term containing sigma
1266  * - "advective" term representing usual upwind
1267  *
1268  * The calculation differs from the reference manual, since ad_coef and dif_coef have different meaning
1269  * than b and A in the manual.
1270  * In calculation of sigma there appears one more csection_lower in the denominator.
1271  */
1272  double sigma = frac_sigma[k][sbi]*arma::dot(dif_coef_edg[0][sbi][k]*fe_values_side.normal_vector(k),fe_values_side.normal_vector(k))*
1273  2*csection_higher[k]*csection_higher[k]/(csection_lower[k]*csection_lower[k]);
1274 
1275  double transport_flux = arma::dot(ad_coef_edg[1][sbi][k], fe_values_side.normal_vector(k));
1276  double por_lower_over_higher = por_lower[k]/por_higher[k];
1277 
1278  comm_flux[0][0] = (sigma-min(0.,transport_flux*por_lower_over_higher))*fv_sb[0]->JxW(k);
1279  comm_flux[0][1] = -(sigma-min(0.,transport_flux*por_lower_over_higher))*fv_sb[0]->JxW(k);
1280  comm_flux[1][0] = -(sigma+max(0.,transport_flux))*fv_sb[0]->JxW(k);
1281  comm_flux[1][1] = (sigma+max(0.,transport_flux))*fv_sb[0]->JxW(k);
1282 
1283  for (int n=0; n<2; n++)
1284  {
1285  if (!feo->dh()->el_is_local(element_id[n])) continue;
1286 
1287  for (unsigned int i=0; i<n_dofs[n]; i++)
1288  for (int m=0; m<2; m++)
1289  for (unsigned int j=0; j<n_dofs[m]; j++)
1290  local_matrix[(i+n*n_dofs[0])*(n_dofs[0]+n_dofs[1]) + m*n_dofs[0] + j] +=
1291  comm_flux[m][n]*fv_sb[m]->shape_value(j,k)*fv_sb[n]->shape_value(i,k);
1292  }
1293  }
1294  ls[sbi]->mat_set_values(n_dofs[0]+n_dofs[1], (int *)side_dof_indices, n_dofs[0]+n_dofs[1], (int *)side_dof_indices, local_matrix);
1295  }
1296  }
1297 
1298 }
1299 
1300 
1301 
1302 
1303 
1304 
1305 template<class Model>
1306 void TransportDG<Model>::set_boundary_conditions()
1307 {
1308  START_TIMER("assemble_bc");
1309  if (Model::balance_ != nullptr)
1310  Model::balance_->start_flux_assembly(Model::subst_idx);
1311  set_boundary_conditions<1>();
1312  set_boundary_conditions<2>();
1313  set_boundary_conditions<3>();
1314  if (Model::balance_ != nullptr)
1315  Model::balance_->finish_flux_assembly(Model::subst_idx);
1316  END_TIMER("assemble_bc");
1317 }
1318 
1319 
1320 template<class Model>
1321 template<unsigned int dim>
1322 void TransportDG<Model>::set_boundary_conditions()
1323 {
1324  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1325  update_values | update_gradients | update_normal_vectors | update_side_JxW_values | update_quadrature_points);
1326  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1327  update_values);
1328  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
1329  vector<int> side_dof_indices(ndofs);
1330  unsigned int loc_b=0;
1331  double local_rhs[ndofs];
1332  vector<PetscScalar> local_flux_balance_vector(ndofs);
1333  PetscScalar local_flux_balance_rhs;
1334  vector<arma::vec> bc_values(qsize, arma::vec(Model::n_substances()));
1335  vector<double> bc_fluxes(qsize),
1336  bc_sigma(qsize),
1337  bc_ref_values(qsize);
1338  vector<double> csection(qsize);
1339  vector<arma::vec3> velocity;
1340 
1341  for (unsigned int loc_el = 0; loc_el < Model::mesh_->get_el_ds()->lsize(); loc_el++)
1342  {
1343  ElementFullIter elm = Model::mesh_->element(feo->dh()->el_index(loc_el));
1344  if (elm->boundary_idx_ == nullptr) continue;
1345 
1346  FOR_ELEMENT_SIDES(elm,si)
1347  {
1348  Edge *edg = elm->side(si)->edge();
1349  if (edg->n_sides > 1) continue;
1350  // skip edges lying not on the boundary
1351  if (edg->side(0)->cond() == NULL) continue;
1352 
1353  if (edg->side(0)->dim() != dim-1)
1354  {
1355  if (edg->side(0)->cond() != nullptr) ++loc_b;
1356  continue;
1357  }
1358 
1359  SideIter side = edg->side(0);
1360  typename DOFHandlerBase::CellIterator cell = Model::mesh_->element.full_iter(side->element());
1361  ElementAccessor<3> ele_acc = side->cond()->element_accessor();
1362 
1363  arma::uvec bc_type;
1364  Model::get_bc_type(ele_acc, bc_type);
1365 
1366  fe_values_side.reinit(cell, side->el_idx());
1367  fsv_rt.reinit(cell, side->el_idx());
1368  calculate_velocity(cell, velocity, fsv_rt);
1369 
1370  Model::compute_advection_diffusion_coefficients(fe_values_side.point_list(), velocity, side->element()->element_accessor(), ad_coef, dif_coef);
1371  data_.cross_section.value_list(fe_values_side.point_list(), side->element()->element_accessor(), csection);
1372  // The b.c. data are fetched for all possible b.c. types since we allow
1373  // different bc_type for each substance.
1374  data_.bc_dirichlet_value.value_list(fe_values_side.point_list(), ele_acc, bc_values);
1375 
1376  feo->dh()->get_dof_indices(cell, (unsigned int *)&(side_dof_indices[0]));
1377 
1378  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
1379  {
1380  fill_n(local_rhs, ndofs, 0);
1381  local_flux_balance_vector.assign(ndofs, 0);
1382  local_flux_balance_rhs = 0;
1383 
1384  double side_flux = 0;
1385  for (unsigned int k=0; k<qsize; k++)
1386  side_flux += arma::dot(ad_coef[sbi][k], fe_values_side.normal_vector(k))*fe_values_side.JxW(k);
1387  double transport_flux = side_flux/side->measure();
1388 
1389  if (bc_type[sbi] == AdvectionDiffusionModel::abc_inflow && side_flux < 0)
1390  {
1391  for (unsigned int k=0; k<qsize; k++)
1392  {
1393  double bc_term = -transport_flux*bc_values[k][sbi]*fe_values_side.JxW(k);
1394  for (unsigned int i=0; i<ndofs; i++)
1395  local_rhs[i] += bc_term*fe_values_side.shape_value(i,k);
1396  }
1397  if (Model::balance_ != nullptr)
1398  for (unsigned int i=0; i<ndofs; i++)
1399  local_flux_balance_rhs -= local_rhs[i];
1400  }
1401  else if (bc_type[sbi] == AdvectionDiffusionModel::abc_dirichlet)
1402  {
1403  for (unsigned int k=0; k<qsize; k++)
1404  {
1405  double bc_term = gamma[sbi][side->cond_idx()]*bc_values[k][sbi]*fe_values_side.JxW(k);
1406  arma::vec3 bc_grad = -bc_values[k][sbi]*fe_values_side.JxW(k)*dg_variant*(arma::trans(dif_coef[sbi][k])*fe_values_side.normal_vector(k));
1407  for (unsigned int i=0; i<ndofs; i++)
1408  local_rhs[i] += bc_term*fe_values_side.shape_value(i,k)
1409  + arma::dot(bc_grad,fe_values_side.shape_grad(i,k));
1410  }
1411  if (Model::balance_ != nullptr)
1412  {
1413  for (unsigned int k=0; k<qsize; k++)
1414  {
1415  for (unsigned int i=0; i<ndofs; i++)
1416  {
1417  local_flux_balance_vector[i] += (arma::dot(ad_coef[sbi][k], fe_values_side.normal_vector(k))*fe_values_side.shape_value(i,k)
1418  - arma::dot(dif_coef[sbi][k]*fe_values_side.shape_grad(i,k),fe_values_side.normal_vector(k))
1419  + gamma[sbi][side->cond_idx()]*fe_values_side.shape_value(i,k))*fe_values_side.JxW(k);
1420  }
1421  }
1422  if (Model::time_->tlevel() > 0)
1423  for (unsigned int i=0; i<ndofs; i++)
1424  local_flux_balance_rhs -= local_rhs[i];
1425  }
1426  }
1427  else if (bc_type[sbi] == AdvectionDiffusionModel::abc_total_flux)
1428  {
1429  Model::get_flux_bc_data(sbi, fe_values_side.point_list(), ele_acc, bc_fluxes, bc_sigma, bc_ref_values);
1430  for (unsigned int k=0; k<qsize; k++)
1431  {
1432  double bc_term = csection[k]*(bc_sigma[k]*bc_ref_values[k]+bc_fluxes[k])*fe_values_side.JxW(k);
1433  for (unsigned int i=0; i<ndofs; i++)
1434  local_rhs[i] += bc_term*fe_values_side.shape_value(i,k);
1435  }
1436 
1437  if (Model::balance_ != nullptr)
1438  {
1439  for (unsigned int i=0; i<ndofs; i++)
1440  {
1441  for (unsigned int k=0; k<qsize; k++)
1442  local_flux_balance_vector[i] += csection[k]*bc_sigma[k]*fe_values_side.JxW(k)*fe_values_side.shape_value(i,k);
1443  local_flux_balance_rhs -= local_rhs[i];
1444  }
1445  }
1446  }
1447  else if (bc_type[sbi] == AdvectionDiffusionModel::abc_diffusive_flux)
1448  {
1449  Model::get_flux_bc_data(sbi, fe_values_side.point_list(), ele_acc, bc_fluxes, bc_sigma, bc_ref_values);
1450  for (unsigned int k=0; k<qsize; k++)
1451  {
1452  double bc_term = csection[k]*(bc_sigma[k]*bc_ref_values[k]+bc_fluxes[k])*fe_values_side.JxW(k);
1453  for (unsigned int i=0; i<ndofs; i++)
1454  local_rhs[i] += bc_term*fe_values_side.shape_value(i,k);
1455  }
1456 
1457  if (Model::balance_ != nullptr)
1458  {
1459  for (unsigned int i=0; i<ndofs; i++)
1460  {
1461  for (unsigned int k=0; k<qsize; k++)
1462  local_flux_balance_vector[i] += csection[k]*(arma::dot(ad_coef[sbi][k], fe_values_side.normal_vector(k)) + bc_sigma[k])*fe_values_side.JxW(k)*fe_values_side.shape_value(i,k);
1463  local_flux_balance_rhs -= local_rhs[i];
1464  }
1465  }
1466  }
1467  else if (bc_type[sbi] == AdvectionDiffusionModel::abc_inflow && side_flux >= 0)
1468  {
1469  if (Model::balance_ != nullptr)
1470  {
1471  for (unsigned int k=0; k<qsize; k++)
1472  {
1473  for (unsigned int i=0; i<ndofs; i++)
1474  local_flux_balance_vector[i] += arma::dot(ad_coef[sbi][k], fe_values_side.normal_vector(k))*fe_values_side.JxW(k)*fe_values_side.shape_value(i,k);
1475  }
1476  }
1477  }
1478  ls[sbi]->rhs_set_values(ndofs, &(side_dof_indices[0]), local_rhs);
1479 
1480  if (Model::balance_ != nullptr)
1481  {
1482  Model::balance_->add_flux_matrix_values(Model::subst_idx[sbi], loc_b, side_dof_indices, local_flux_balance_vector);
1483  Model::balance_->add_flux_vec_value(Model::subst_idx[sbi], loc_b, local_flux_balance_rhs);
1484  }
1485  }
1486  ++loc_b;
1487  }
1488  }
1489 }
1490 
1491 
1492 
1493 template<class Model>
1494 template<unsigned int dim>
1495 void TransportDG<Model>::calculate_velocity(const typename DOFHandlerBase::CellIterator &cell, vector<arma::vec3> &velocity, FEValuesBase<dim,3> &fv)
1496 {
1497  OLD_ASSERT(cell->dim() == dim, "Element dimension mismatch!");
1498 
1499  velocity.resize(fv.n_points());
1500 
1501  for (unsigned int k=0; k<fv.n_points(); k++)
1502  {
1503  velocity[k].zeros();
1504  for (unsigned int sid=0; sid<cell->n_sides(); sid++)
1505  velocity[k] += fv.shape_vector(sid,k) * Model::mh_dh->side_flux( *(cell->side(sid)) );
1506  }
1507 }
1508 
1509 
1510 
1511 
1512 
1513 template<class Model>
1514 void TransportDG<Model>::set_DG_parameters_edge(const Edge &edg,
1515  const int s1,
1516  const int s2,
1517  const int K_size,
1518  const vector<arma::mat33> &K1,
1519  const vector<arma::mat33> &K2,
1520  const vector<double> &fluxes,
1521  const arma::vec3 &normal_vector,
1522  const double alpha1,
1523  const double alpha2,
1524  double &gamma,
1525  double *omega,
1526  double &transport_flux)
1527 {
1528  double delta[2];
1529  double h = 0;
1530  double local_alpha = 0;
1531 
1532  OLD_ASSERT(edg.side(s1)->valid(), "Invalid side of an edge.");
1533  SideIter s = edg.side(s1);
1534 
1535  // calculate the side diameter
1536  if (s->dim() == 0)
1537  {
1538  h = 1;
1539  }
1540  else
1541  {
1542  for (unsigned int i=0; i<s->n_nodes(); i++)
1543  for (unsigned int j=i+1; j<s->n_nodes(); j++)
1544  h = max(h, s->node(i)->distance(*s->node(j)));
1545  }
1546 
1547  // calculate the total in- and out-flux through the edge
1548  double pflux = 0, nflux = 0;
1549  for (int i=0; i<edg.n_sides; i++)
1550  {
1551  if (fluxes[i] > 0)
1552  pflux += fluxes[i];
1553  else
1554  nflux += fluxes[i];
1555  }
1556 
1557  // calculate the flux from s1 to s2
1558  if (fluxes[s2] > 0 && fluxes[s1] < 0 && s1 < s2)
1559  transport_flux = fluxes[s1]*fabs(fluxes[s2]/pflux);
1560  else if (fluxes[s2] < 0 && fluxes[s1] > 0 && s1 < s2)
1561  transport_flux = fluxes[s1]*fabs(fluxes[s2]/nflux);
1562  else if (s1==s2)
1563  transport_flux = fluxes[s1];
1564  else
1565  transport_flux = 0;
1566 
1567  gamma = 0.5*fabs(transport_flux);
1568 
1569 
1570  // determine local DG penalty
1571  local_alpha = max(alpha1, alpha2);
1572 
1573  if (s1 == s2)
1574  {
1575  omega[0] = 1;
1576 
1577  // delta is set to the average value of Kn.n on the side
1578  delta[0] = 0;
1579  for (int k=0; k<K_size; k++)
1580  delta[0] += dot(K1[k]*normal_vector,normal_vector);
1581  delta[0] /= K_size;
1582 
1583  gamma += local_alpha/h*(delta[0]);
1584  }
1585  else
1586  {
1587  delta[0] = 0;
1588  delta[1] = 0;
1589  for (int k=0; k<K_size; k++)
1590  {
1591  delta[0] += dot(K1[k]*normal_vector,normal_vector);
1592  delta[1] += dot(K2[k]*normal_vector,normal_vector);
1593  }
1594  delta[0] /= K_size;
1595  delta[1] /= K_size;
1596 
1597  double delta_sum = delta[0] + delta[1];
1598 
1599  if (delta_sum > numeric_limits<double>::epsilon())
1600  {
1601  omega[0] = delta[1]/delta_sum;
1602  omega[1] = delta[0]/delta_sum;
1603  gamma += local_alpha/h*(delta[0]*delta[1]/delta_sum);
1604  }
1605  else
1606  for (int i=0; i<2; i++) omega[i] = 0;
1607  }
1608 }
1609 
1610 
1611 
1612 
1613 
1614 
1615 template<class Model>
1616 void TransportDG<Model>::set_DG_parameters_boundary(const SideIter side,
1617  const int K_size,
1618  const vector<arma::mat33> &K,
1619  const double flux,
1620  const arma::vec3 &normal_vector,
1621  const double alpha,
1622  double &gamma)
1623 {
1624  double delta = 0, h = 0;
1625 
1626  // calculate the side diameter
1627  if (side->dim() == 0)
1628  {
1629  h = 1;
1630  }
1631  else
1632  {
1633  for (unsigned int i=0; i<side->n_nodes(); i++)
1634  for (unsigned int j=i+1; j<side->n_nodes(); j++)
1635  h = max(h, side->node(i)->distance( *side->node(j) ));
1636  }
1637 
1638  // delta is set to the average value of Kn.n on the side
1639  for (int k=0; k<K_size; k++)
1640  delta += dot(K[k]*normal_vector,normal_vector);
1641  delta /= K_size;
1642 
1643  gamma = 0.5*fabs(flux) + alpha/h*delta;
1644 }
1645 
1646 
1647 
1648 
1649 
1650 template<class Model>
1651 void TransportDG<Model>::set_initial_condition()
1652 {
1653  START_TIMER("set_init_cond");
1654  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
1655  ls[sbi]->start_allocation();
1656  prepare_initial_condition<1>();
1657  prepare_initial_condition<2>();
1658  prepare_initial_condition<3>();
1659 
1660  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
1661  ls[sbi]->start_add_assembly();
1662  prepare_initial_condition<1>();
1663  prepare_initial_condition<2>();
1664  prepare_initial_condition<3>();
1665 
1666  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
1667  {
1668  ls[sbi]->finish_assembly();
1669  ls[sbi]->solve();
1670  }
1671  END_TIMER("set_init_cond");
1672 }
1673 
1674 template<class Model>
1675 template<unsigned int dim>
1676 void TransportDG<Model>::prepare_initial_condition()
1677 {
1678  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
1679  update_values | update_JxW_values | update_quadrature_points);
1680  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
1681  unsigned int dof_indices[ndofs];
1682  double matrix[ndofs*ndofs], rhs[ndofs];
1683  std::vector<arma::vec> init_values(qsize);
1684 
1685  for (unsigned int k=0; k<qsize; k++)
1686  init_values[k].resize(Model::n_substances());
1687 
1688  for (unsigned int i_cell=0; i_cell<Model::mesh_->get_el_ds()->lsize(); i_cell++)
1689  {
1690  typename DOFHandlerBase::CellIterator elem = Model::mesh_->element(feo->dh()->el_index(i_cell));
1691  if (elem->dim() != dim) continue;
1692 
1693  ElementAccessor<3> ele_acc = elem->element_accessor();
1694  feo->dh()->get_dof_indices(elem, dof_indices);
1695  fe_values.reinit(elem);
1696 
1697  Model::compute_init_cond(fe_values.point_list(), ele_acc, init_values);
1698 
1699  for (unsigned int sbi=0; sbi<Model::n_substances(); sbi++)
1700  {
1701  for (unsigned int i=0; i<ndofs; i++)
1702  {
1703  rhs[i] = 0;
1704  for (unsigned int j=0; j<ndofs; j++)
1705  matrix[i*ndofs+j] = 0;
1706  }
1707 
1708  for (unsigned int k=0; k<qsize; k++)
1709  {
1710  double rhs_term = init_values[k](sbi)*fe_values.JxW(k);
1711 
1712  for (unsigned int i=0; i<ndofs; i++)
1713  {
1714  for (unsigned int j=0; j<ndofs; j++)
1715  matrix[i*ndofs+j] += fe_values.shape_value(i,k)*fe_values.shape_value(j,k)*fe_values.JxW(k);
1716 
1717  rhs[i] += fe_values.shape_value(i,k)*rhs_term;
1718  }
1719  }
1720  ls[sbi]->set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, matrix, rhs);
1721  }
1722  }
1723 }
1724 
1725 
1726 template<class Model>
1727 void TransportDG<Model>::get_par_info(int * &el_4_loc, Distribution * &el_ds)
1728 {
1729  el_4_loc = Model::mesh_->get_el_4_loc();
1730  el_ds = Model::mesh_->get_el_ds();
1731 }
1732 
1733 
1734 template<class Model>
1735 void TransportDG<Model>::update_after_reactions(bool solution_changed)
1736 {
1737  if (solution_changed)
1738  {
1739  for (unsigned int i_cell=0; i_cell<Model::mesh_->get_el_ds()->lsize(); i_cell++)
1740  {
1741  typename DOFHandlerBase::CellIterator elem = Model::mesh_->element(feo->dh()->el_index(i_cell));
1742 
1743  unsigned int n_dofs;
1744  switch (elem->dim())
1745  {
1746  case 1:
1747  n_dofs = feo->fe<1>()->n_dofs();
1748  break;
1749  case 2:
1750  n_dofs = feo->fe<2>()->n_dofs();
1751  break;
1752  case 3:
1753  n_dofs = feo->fe<3>()->n_dofs();
1754  break;
1755  }
1756 
1757  unsigned int dof_indices[n_dofs];
1758  feo->dh()->get_dof_indices(elem, dof_indices);
1759 
1760  for (unsigned int sbi=0; sbi<Model::n_substances(); ++sbi)
1761  {
1762  double old_average = 0;
1763  for (unsigned int j=0; j<n_dofs; ++j)
1764  old_average += ls[sbi]->get_solution_array()[dof_indices[j]-feo->dh()->distr()->begin()];
1765  old_average /= n_dofs;
1766 
1767  for (unsigned int j=0; j<n_dofs; ++j)
1768  ls[sbi]->get_solution_array()[dof_indices[j]-feo->dh()->distr()->begin()] += solution_elem_[sbi][i_cell] - old_average;
1769  }
1770  }
1771  }
1772  // update mass_vec for the case that mass matrix changes in next time step
1773  for (unsigned int sbi=0; sbi<Model::n_substances(); ++sbi)
1774  MatMult(*(ls_dt[sbi]->get_matrix()), ls[sbi]->get_solution(), mass_vec[sbi]);
1775 }
1776 
1777 template<class Model>
1778 int *TransportDG<Model>::get_row_4_el()
1779 {
1780  return Model::mesh_->get_row_4_el();
1781 }
1782 
1783 
1784 
1785 
1786 
1787 
1788 
1789 
1790 template class TransportDG<ConcentrationTransportModel>;
1791 template class TransportDG<HeatTransferModel>;
1792 
1793 
1794 
1795 
EquationBase::time
TimeGovernor & time()
Definition: equation.hh:146
TransportDG::input_rec
Input::Record input_rec
Record with input specification.
Definition: transport_dg.hh:466
mapping_p1.hh
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell...
update_values
Shape function values.
Definition: update_flags.hh:87
EquationBase::eq_data_
FieldSet * eq_data_
Definition: equation.hh:230
FEValuesBase::JxW
double JxW(const unsigned int point_no)
Return the product of Jacobian determinant and the quadrature weight at given quadrature point...
Definition: fe_values.hh:292
TransportDG::assemble_fluxes_element_element
void assemble_fluxes_element_element()
Assembles the fluxes between elements of the same dimension.
Definition: transport_dg.cc:959
TransportDG::set_sources
void set_sources()
Assembles the right hand side due to volume sources.
Definition: transport_dg.cc:883
TransportDG::set_boundary_conditions
void set_boundary_conditions()
Assembles the r.h.s. components corresponding to the Dirichlet boundary conditions.
Definition: transport_dg.cc:1306
update_side_JxW_values
Transformed quadrature weight for cell sides.
Definition: update_flags.hh:153
Input::Array
Accessor to input data conforming to declared Array.
Definition: accessors.hh:552
TransportDG::assemble_fluxes_element_side
void assemble_fluxes_element_side()
Assembles the fluxes between elements of different dimensions.
Definition: transport_dg.cc:1193
FieldFlag::in_main_matrix
static constexpr Mask in_main_matrix
A field is part of main &quot;stiffness matrix&quot; of the equation.
Definition: field_flag.hh:49
linsys_PETSC.hh
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
TransportDG::calculate_concentration_matrix
void calculate_concentration_matrix()
Definition: transport_dg.cc:626
TransportDG
Transport with dispersion implemented using discontinuous Galerkin method.
Definition: transport_dg.hh:130
TransportDG::EqData::region_id
Field< 3, FieldValue< 3 >::Integer > region_id
Definition: transport_dg.hh:141
Input::Type::Default
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:50
FieldCommon::flags_add
FieldCommon & flags_add(FieldFlag::Flags::Mask mask)
Definition: field_common.hh:153
Node::distance
double distance(const Node &n2) const
Definition: nodes.hh:86
LinSys_PETSC
Definition: linsys_PETSC.hh:31
TransportDG::assemble_mass_matrix
void assemble_mass_matrix()
Assembles the mass matrix.
Definition: transport_dg.cc:727
Side::cond
Boundary * cond() const
Definition: side_impl.hh:70
TransportDG::update_solution
void update_solution() override
Computes the solution in one time instant.
Definition: transport_dg.cc:492
Input::Type::Default::obligatory
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:99
TransportDG::dg_variant
int dg_variant
DG variant ((non-)symmetric/incomplete.
Definition: transport_dg.hh:420
TransportDG::prepare_initial_condition
void prepare_initial_condition()
Assembles the auxiliary linear system to calculate the initial solution as L^2-projection of the pres...
Definition: transport_dg.cc:1676
Mesh
Definition: mesh.h:95
generic_field.hh
Fields computed from the mesh data.
TransportDG::set_initial_condition
void set_initial_condition()
Sets the initial condition.
Definition: transport_dg.cc:1651
fe_values.hh
Class FEValues calculates finite element data on the actual cells such as shape function values...
FOR_ELEMENT_SIDES
#define FOR_ELEMENT_SIDES(i, j)
Definition: elements.h:189
flow::FullIterator::index
int index() const
Definition: sys_vector.hh:78
Edge::n_sides
int n_sides
Definition: edges.h:36
Edge
Definition: edges.h:26
Side::valid
bool valid() const
Definition: side_impl.hh:86
LinSys_PETSC::set_from_input
void set_from_input(const Input::Record in_rec)
Definition: linsys_PETSC.cc:428
AVERAGE
#define AVERAGE(i, k, side_id)
TransportDG::assemble_stiffness_matrix
void assemble_stiffness_matrix()
Assembles the stiffness matrix.
Definition: transport_dg.cc:796
concentration_model.hh
Discontinuous Galerkin method for equation of transport with dispersion.
Input::Type::Integer
Class for declaration of the integral input data.
Definition: type_base.hh:465
FieldCommon::units
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
Definition: field_common.hh:118
TransportDG::EqData::dg_penalty
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
Definition: transport_dg.hh:140
TransportDG::get_dg_variant_selection_input_type
static const Input::Type::Selection & get_dg_variant_selection_input_type()
Input type for the DG variant selection.
Definition: transport_dg.cc:48
FEValuesBase::n_points
unsigned int n_points()
Returns the number of quadrature points.
Definition: fe_values.hh:330
Input::Type::Array
Class for declaration of inputs sequences.
Definition: type_base.hh:321
TransportDG::EqData::output_fields
EquationOutput output_fields
Definition: transport_dg.hh:143
TransportDG::calculate_cumulative_balance
void calculate_cumulative_balance()
Definition: transport_dg.cc:686
TransportDG::assemble_fluxes_boundary
void assemble_fluxes_boundary()
Assembles the fluxes on the boundary.
Definition: transport_dg.cc:1090
QGauss
Symmetric Gauss-Legendre quadrature formulae on simplices.
Definition: darcy_flow_mh.hh:69
WAVERAGE
#define WAVERAGE(i, k, side_id)
sys_profiler.hh
TransportDG::TransportDG
TransportDG(Mesh &init_mesh, const Input::Record in_rec)
Constructor.
Definition: transport_dg.cc:249
Element::dim
unsigned int dim() const
Definition: element_impls.hh:26
FEValuesBase::normal_vector
arma::vec::fixed< spacedim > normal_vector(unsigned int point_no)
Returns the normal vector to a side at given quadrature point.
Definition: fe_values.hh:322
OLD_ASSERT
#define OLD_ASSERT(...)
Definition: global_defs.h:131
TransportDG::initialize
void initialize() override
Definition: transport_dg.cc:282
FE_P_disc
Discontinuous Lagrangean finite element on dim dimensional simplex.
Definition: fe_p.hh:202
update_quadrature_points
Transformed quadrature points.
Definition: update_flags.hh:102
equation_output.hh
TransportDG::output_data
void output_data()
Postprocesses the solution and writes to output file.
Definition: transport_dg.cc:666
TransportDG::EqData::fracture_sigma
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
Definition: transport_dg.hh:139
TransportDG::preallocate
void preallocate()
Definition: transport_dg.cc:467
Side::dim
unsigned int dim() const
Definition: side_impl.hh:37
TransportDG::feo
FEObjects * feo
Finite element objects.
Definition: transport_dg.hh:414
FieldFlag::in_time_term
static constexpr Mask in_time_term
A field is part of time term of the equation.
Definition: field_flag.hh:47
FieldCommon::input_default
FieldCommon & input_default(const string &input_default)
Definition: field_common.hh:105
FE_RT0
Raviart-Thomas element of order 0.
Definition: fe_rt.hh:33
fe_p.hh
Definitions of basic Lagrangean finite elements with polynomial shape functions.
FieldFlag::equation_external_output
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
Definition: field_flag.hh:58
Neighbour::side
SideIter side()
Definition: neighbours_impl.hh:26
Input::Type::Selection::add_value
Selection & add_value(const int value, const std::string &key, const std::string &description="")
Adds one new value with name given by key to the Selection.
Definition: type_selection.cc:49
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:277
Input::Record::val
const Ret val(const string &key) const
Definition: accessors_impl.hh:30
GenericField::region_id
static auto region_id(Mesh &mesh) -> IndexField
Definition: generic_field.impl.hh:28
xprintf
#define xprintf(...)
Definition: system.hh:87
START_TIMER
#define START_TIMER(tag)
Starts a timer with specified tag.
Definition: sys_profiler.hh:112
FEValuesBase::fe
FiniteElement< dim, spacedim > * fe
The used finite element.
Definition: fe_values.hh:383
DOFHandlerMultiDim
Provides the numbering of the finite element degrees of freedom on the computational mesh...
Definition: dofhandler.hh:248
EquationBase::mesh_
Mesh * mesh_
Definition: equation.hh:221
Input::Type::Record::declare_key
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:468
TransportDG::dg_order
unsigned int dg_order
Polynomial order of finite elements.
Definition: transport_dg.hh:423
TransportDG::output_vector_gather
void output_vector_gather()
Definition: transport_dg.cc:410
FieldFlag::in_rhs
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Definition: field_flag.hh:51
MappingP1
Affine mapping between reference and actual cell.
Definition: mapping_p1.hh:53
update_gradients
Shape function gradients.
Definition: update_flags.hh:95
FEValuesBase::shape_vector
arma::vec::fixed< spacedim > shape_vector(const unsigned int function_no, const unsigned int point_no)
Return the value of the function_no-th shape function at the point_no-th quadrature point...
Definition: fe_values.hh:254
output_time.hh
Neighbour::element
ElementIter element()
Definition: neighbours_impl.hh:42
Side::measure
double measure() const
Definition: sides.cc:29
FEValuesBase::shape_value
double shape_value(const unsigned int function_no, const unsigned int point_no)
Return the value of the function_no-th shape function at the point_no-th quadrature point...
Definition: fe_values.hh:226
FEObjects::dh
DOFHandlerMultiDim * dh()
Definition: transport_dg.cc:215
ElementAccessor::region
Region region() const
Definition: accessors.hh:98
update_normal_vectors
Normal vectors.
Definition: update_flags.hh:124
transport_dg.hh
Discontinuous Galerkin method for equation of transport with dispersion.
MPI_Comm_rank
#define MPI_Comm_rank
Definition: mpi.h:236
epsilon
const double epsilon
Definition: mathfce.h:23
TransportOperatorSplitting::output_data
void output_data() override
Write computed fields.
Definition: transport_operator_splitting.cc:232
FieldCommon::description
FieldCommon & description(const string &description)
Definition: field_common.hh:93
TransportDG::data_
EqData data_
Field data for model parameters.
Definition: transport_dg.hh:405
TransportDG::allocation_done
bool allocation_done
Indicates whether matrices have been preallocated.
Definition: transport_dg.hh:499
multi_field.hh
Side::cond_idx
unsigned int cond_idx() const
Definition: side_impl.hh:75
TransportDG::get_input_type
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
Definition: transport_dg.cc:74
advection_diffusion_model.hh
Discontinuous Galerkin method for equation of transport with dispersion.
Input::Type::Selection::close
const Selection & close() const
Close the Selection, no more values can be added.
Definition: type_selection.cc:62
Side::element
ElementFullIter element() const
Definition: side_impl.hh:52
TransportDG::update_after_reactions
void update_after_reactions(bool solution_changed)
Definition: transport_dg.cc:1735
PrgErr
Definition: system.hh:59
FESideValues::reinit
void reinit(ElementFullIter &cell, unsigned int sid)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:225
FEValues::reinit
void reinit(ElementFullIter &cell)
Update cell-dependent data (gradients, Jacobians etc.)
Definition: fe_values.cc:157
RegionIdx::bulk_idx
unsigned int bulk_idx() const
Returns index of the region in the bulk set.
Definition: region.hh:90
LinSys
Abstract linear system class.
Definition: la_linsys_new.hh:169
quadrature_lib.hh
Definitions of particular quadrature rules on simplices.
FieldCommon::name
FieldCommon & name(const string &name)
Definition: field_common.hh:86
END_TIMER
#define END_TIMER(tag)
Ends a timer with specified tag.
Definition: sys_profiler.hh:146
heat_model.hh
Discontinuous Galerkin method for equation of transport with dispersion.
FEValuesBase::point_list
vector< arma::vec::fixed< spacedim > > & point_list()
Return coordinates of all quadrature points in the actual cell system.
Definition: fe_values.hh:311
LinSys_PETSC::get_input_type
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:31
Side::el_idx
unsigned int el_idx() const
Definition: side_impl.hh:81
TransportDG::calculate_instant_balance
void calculate_instant_balance()
Definition: transport_dg.cc:712
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:171
LimitSide::left
TransportDG::assemble_volume_integrals
void assemble_volume_integrals()
Assembles the volume integrals into the stiffness matrix.
Definition: transport_dg.cc:829
FieldCommon::flags
FieldCommon & flags(FieldFlag::Flags::Mask mask)
Definition: field_common.hh:147
FiniteElement
Abstract class for the description of a general finite element on a reference simplex in dim dimensio...
Definition: dofhandler.hh:29
FEValuesBase
Base class for FEValues and FESideValues.
Definition: fe_values.hh:32
Side::node
const Node * node(unsigned int i) const
Definition: side_impl.hh:46
TransportDG::get_row_4_el
int * get_row_4_el()
Definition: transport_dg.cc:1778
TransportDG::calculate_velocity
void calculate_velocity(const ElementFullIter &cell, std::vector< arma::vec3 > &velocity, FEValuesBase< dim, 3 > &fv)
Calculates the velocity field on a given dim dimensional cell.
field_fe.hh
FEValuesBase::shape_grad
arma::vec::fixed< spacedim > shape_grad(const unsigned int function_no, const unsigned int point_no)
Return the gradient of the function_no-th shape function at the point_no-th quadrature point...
Definition: fe_values.hh:239
UnitSI::dimensionless
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:53
TransportDG::~TransportDG
~TransportDG()
Destructor.
Definition: transport_dg.cc:376
TransportDG::set_DG_parameters_edge
void set_DG_parameters_edge(const Edge &edg, const int s1, const int s2, const int K_size, const std::vector< arma::mat33 > &K1, const std::vector< arma::mat33 > &K2, const std::vector< double > &fluxes, const arma::vec3 &normal_vector, const double alpha1, const double alpha2, double &gamma, double *omega, double &transport_flux)
Sets up some parameters of the DG method for two sides of an edge.
Definition: transport_dg.cc:1514
Edge::side
SideIter side(const unsigned int i) const
Definition: edges.h:31
JUMP
#define JUMP(i, k, side_id)
Input::Type::Selection
Template for classes storing finite set of named values.
Definition: type_selection.hh:63
TransportDG::set_DG_parameters_boundary
void set_DG_parameters_boundary(const SideIter side, const int K_size, const std::vector< arma::mat33 > &K, const double flux, const arma::vec3 &normal_vector, const double alpha, double &gamma)
Sets up parameters of the DG method on a given boundary edge.
Definition: transport_dg.cc:1616
fe_rt.hh
Definitions of Raviart-Thomas finite elements.
TransportDG::get_par_info
void get_par_info(int *&el_4_loc, Distribution *&el_ds)
Definition: transport_dg.cc:1727
FieldValue_
Definition: field_values.hh:231
FLOW123D_FORCE_LINK_IN_CHILD
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:180
Boundary::element_accessor
ElementAccessor< 3 > element_accessor()
Definition: boundaries.cc:47
TransportDG::zero_time_step
void zero_time_step() override
Initialize solution in the zero time.
Definition: transport_dg.cc:427
FEObjects::FEObjects
FEObjects(Mesh *mesh_, unsigned int fe_order)
Definition: transport_dg.cc:116
Side::n_nodes
unsigned int n_nodes() const
Definition: side_impl.hh:33
update_JxW_values
Transformed quadrature weights.
Definition: update_flags.hh:114
Generated by   doxygen 1.8.5