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