Flow123d  JS_before_hm-1728-g5d7ebd453
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/index_types.hh"
20 #include "system/sys_profiler.hh"
22 
23 #include "io/output_time.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"
31 #include "la/linsys_PETSC.hh"
32 #include "coupling/balance.hh"
36 #include "transport/heat_model.hh"
37 #include "transport/assembly_dg.hh"
38 
39 #include "fields/multi_field.hh"
40 #include "fields/generic_field.hh"
41 #include "input/factory.hh"
43 #include "mesh/accessors.hh"
44 
45 FLOW123D_FORCE_LINK_IN_CHILD(concentrationTransportModel)
47 
48 
49 
50 using namespace Input::Type;
51 
52 template<class Model>
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>
80  std::string equation_name = std::string(Model::ModelEqData::name()) + "_DG";
81  return Model::get_input_type("DG", "Discontinuous Galerkin (DG) solver")
83  "Solver for the linear system.")
84  .declare_key("input_fields", Array(
86  .make_field_descriptor_type(equation_name)),
88  "Input fields of the equation.")
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("init_projection", Bool(), Default("true"),
94  "If true, use DG projection of the initial condition field."
95  "Otherwise, evaluate initial condition field directly (well suited for reading native data).")
96  .declare_key("output",
97  EqFields().output_fields.make_output_type(equation_name, ""),
98  IT::Default("{ \"fields\": [ " + Model::ModelEqData::default_output_field() + "] }"),
99  "Specification of output fields and output times.")
100  .close();
101 }
102 
103 template<class Model>
105  Input::register_class< TransportDG<Model>, Mesh &, const Input::Record>(std::string(Model::ModelEqData::name()) + "_DG") +
107 
108 
109 
110 template<class Model>
112 {
113  *this+=fracture_sigma
114  .name("fracture_sigma")
115  .description(
116  "Coefficient of diffusive transfer through fractures (for each substance).")
118  .input_default("1.0")
120 
121  *this+=dg_penalty
122  .name("dg_penalty")
123  .description(
124  "Penalty parameter influencing the discontinuity of the solution (for each substance). "
125  "Its default value 1 is sufficient in most cases. Higher value diminishes the inter-element jumps.")
127  .input_default("1.0")
129 
130  *this += region_id.name("region_id")
133  .description("Region ids.");
134 
135  *this += subdomain.name("subdomain")
138  .description("Subdomain ids of the domain decomposition.");
139 
140 
141  // add all input fields to the output list
142  output_fields += *this;
143 
144 }
145 
146 
147 
148 // return the ratio of longest and shortest edge
149 template<class Model>
151 {
152  double h_max = 0, h_min = numeric_limits<double>::infinity();
153  for (unsigned int i=0; i<e->n_nodes(); i++)
154  for (unsigned int j=i+1; j<e->n_nodes(); j++)
155  {
156  double dist = arma::norm(*e.node(i) - *e.node(j));
157  h_max = max(h_max, dist);
158  h_min = min(h_min, dist);
159  }
160  return h_max/h_min;
161 }
162 
163 
164 
165 template<class Model>
167  const int K_size,
168  const vector<arma::mat33> &K,
169  const double flux,
170  const arma::vec3 &normal_vector,
171  const double alpha,
172  double &gamma)
173 {
174  double delta = 0, h = 0;
175 
176  // calculate the side diameter
177  if (side.dim() == 0)
178  {
179  h = 1;
180  }
181  else
182  {
183  for (unsigned int i=0; i<side.n_nodes(); i++)
184  for (unsigned int j=i+1; j<side.n_nodes(); j++) {
185  double dist = arma::norm(*side.node(i) - *side.node(j));
186  h = max(h, dist);
187  }
188 
189  }
190 
191  // delta is set to the average value of Kn.n on the side
192  for (int k=0; k<K_size; k++)
193  delta += dot(K[k]*normal_vector,normal_vector);
194  delta /= K_size;
195 
196  gamma = 0.5*fabs(flux) + alpha/h*delta*elem_anisotropy(side.element());
197 }
198 
199 
200 
201 template<typename Model>
203  : Model(init_mesh, in_rec),
204  input_rec(in_rec),
205  allocation_done(false)
206 {
207  // Can not use name() + "constructor" here, since START_TIMER only accepts const char *
208  // due to constexpr optimization.
209  START_TIMER(Model::ModelEqData::name());
210  // Check that Model is derived from AdvectionDiffusionModel.
212 
213  eq_data_ = make_shared<EqData>();
214  eq_fields_ = make_shared<EqFields>();
215  eq_fields_->add_coords_field();
216  this->eq_fieldset_ = eq_fields_.get();
217  Model::init_balance(in_rec);
218 
219 
220  // Set up physical parameters.
221  eq_fields_->set_mesh(init_mesh);
222  eq_fields_->region_id = GenericField<3>::region_id(*Model::mesh_);
223  eq_fields_->subdomain = GenericField<3>::subdomain(*Model::mesh_);
224 
225 
226  // DG data parameters
227  eq_data_->dg_variant = in_rec.val<DGVariant>("dg_variant");
228  eq_data_->dg_order = in_rec.val<unsigned int>("dg_order");
229 
230  Model::init_from_input(in_rec);
231 
232  MixedPtr<FE_P_disc> fe(eq_data_->dg_order);
233  shared_ptr<DiscreteSpace> ds = make_shared<EqualOrderDiscreteSpace>(Model::mesh_, fe);
234  eq_data_->dh_ = make_shared<DOFHandlerMultiDim>(*Model::mesh_);
235  eq_data_->dh_->distribute_dofs(ds);
236  //DebugOut().fmt("TDG: solution size {}\n", eq_data_->dh_->n_global_dofs());
237 
238 }
239 
240 
241 template<class Model>
243 {
244  eq_fields_->set_components(eq_data_->substances_.names());
245  eq_fields_->set_input_list( input_rec.val<Input::Array>("input_fields"), *(Model::time_) );
246  eq_data_->set_time_governor(Model::time_);
247  eq_data_->balance_ = this->balance();
248  eq_fields_->initialize();
249 
250  // DG stabilization parameters on boundary edges
251  eq_data_->gamma.resize(eq_data_->n_substances());
252  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
253  eq_data_->gamma[sbi].resize(Model::mesh_->boundary_.size());
254 
255  // Resize coefficient arrays
256  eq_data_->max_edg_sides = max(Model::mesh_->max_edge_sides(1), max(Model::mesh_->max_edge_sides(2), Model::mesh_->max_edge_sides(3)));
257  ret_sources.resize(eq_data_->n_substances());
258  ret_sources_prev.resize(eq_data_->n_substances());
259 
260  eq_data_->output_vec.resize(eq_data_->n_substances());
261  eq_fields_->output_field.set_components(eq_data_->substances_.names());
262  eq_fields_->output_field.set_mesh(*Model::mesh_);
263  eq_fields_->output_type(OutputTime::CORNER_DATA);
264 
265  eq_fields_->output_field.setup_components();
266  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
267  {
268  // create shared pointer to a FieldFE, pass FE data and push this FieldFE to output_field on all regions
269  auto output_field_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_);
270  eq_fields_->output_field[sbi].set(output_field_ptr, 0);
271  eq_data_->output_vec[sbi] = output_field_ptr->vec();
272  }
273 
274  // set time marks for writing the output
275  eq_fields_->output_fields.initialize(Model::output_stream_, Model::mesh_, input_rec.val<Input::Record>("output"), this->time());
276 
277  // equation default PETSc solver options
278  std::string petsc_default_opts;
279  if (eq_data_->dh_->distr()->np() == 1)
280  petsc_default_opts = "-ksp_type bcgs -pc_type ilu -pc_factor_levels 2 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
281  else
282  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";
283 
284  // allocate matrix and vector structures
285  eq_data_->ls = new LinSys*[eq_data_->n_substances()];
286  eq_data_->ls_dt = new LinSys*[eq_data_->n_substances()];
287  eq_data_->conc_fe.resize(eq_data_->n_substances());
288 
289  MixedPtr<FE_P_disc> fe(0);
290  shared_ptr<DiscreteSpace> ds = make_shared<EqualOrderDiscreteSpace>(Model::mesh_, fe);
291  eq_data_->dh_p0 = make_shared<DOFHandlerMultiDim>(*Model::mesh_);
292  eq_data_->dh_p0->distribute_dofs(ds);
293 
294  stiffness_matrix.resize(eq_data_->n_substances(), nullptr);
295  mass_matrix.resize(eq_data_->n_substances(), nullptr);
296  rhs.resize(eq_data_->n_substances(), nullptr);
297  mass_vec.resize(eq_data_->n_substances(), nullptr);
298  eq_data_->ret_vec.resize(eq_data_->n_substances(), nullptr);
299 
300  for (unsigned int sbi = 0; sbi < eq_data_->n_substances(); sbi++) {
301  eq_data_->ls[sbi] = new LinSys_PETSC(eq_data_->dh_->distr().get(), petsc_default_opts);
302  ( (LinSys_PETSC *)eq_data_->ls[sbi] )->set_from_input( input_rec.val<Input::Record>("solver") );
303  eq_data_->ls[sbi]->set_solution(eq_data_->output_vec[sbi].petsc_vec());
304 
305  eq_data_->ls_dt[sbi] = new LinSys_PETSC(eq_data_->dh_->distr().get(), petsc_default_opts);
306  ( (LinSys_PETSC *)eq_data_->ls_dt[sbi] )->set_from_input( input_rec.val<Input::Record>("solver") );
307 
308  eq_data_->conc_fe[sbi] = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_p0);
309 
310  VecDuplicate(eq_data_->ls[sbi]->get_solution(), &eq_data_->ret_vec[sbi]);
311  }
312 
313 
314  init_projection = input_rec.val<bool>("init_projection");
315 
316  // create assemblation object, finite element structures and distribute DOFs
317  mass_assembly_ = new GenericAssembly< MassAssemblyDim >(eq_fields_.get(), eq_data_.get());
318  stiffness_assembly_ = new GenericAssembly< StiffnessAssemblyDim >(eq_fields_.get(), eq_data_.get());
319  sources_assembly_ = new GenericAssembly< SourcesAssemblyDim >(eq_fields_.get(), eq_data_.get());
320  bdr_cond_assembly_ = new GenericAssembly< BdrConditionAssemblyDim >(eq_fields_.get(), eq_data_.get());
321 
322  if(init_projection)
323  init_assembly_ = new GenericAssembly< InitProjectionAssemblyDim >(eq_fields_.get(), eq_data_.get());
324  else
325  init_assembly_ = new GenericAssembly< InitConditionAssemblyDim >(eq_fields_.get(), eq_data_.get());
326 
327  // initialization of balance object
328  Model::balance_->allocate(eq_data_->dh_->distr()->lsize(), mass_assembly_->eval_points()->max_size());
329 
330  int qsize = mass_assembly_->eval_points()->max_size();
331  eq_data_->dif_coef.resize(eq_data_->n_substances());
332  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
333  {
334  eq_data_->dif_coef[sbi].resize(qsize);
335  }
336 
337  eq_fields_->init_condition.setup_components();
338  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
339  {
340  eq_fields_->init_condition[sbi].add_factory( std::make_shared<FieldFE<3, FieldValue<3>::Scalar>::NativeFactory>(sbi, eq_data_->dh_));
341  }
342 }
343 
344 
345 template<class Model>
347 {
348  delete Model::time_;
349 
350  if (eq_data_->gamma.size() > 0) {
351  // initialize called
352 
353  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
354  {
355  delete eq_data_->ls[i];
356  delete eq_data_->ls_dt[i];
357 
358  if (stiffness_matrix[i])
359  chkerr(MatDestroy(&stiffness_matrix[i]));
360  if (mass_matrix[i])
361  chkerr(MatDestroy(&mass_matrix[i]));
362  if (rhs[i])
363  chkerr(VecDestroy(&rhs[i]));
364  if (mass_vec[i])
365  chkerr(VecDestroy(&mass_vec[i]));
366  if (eq_data_->ret_vec[i])
367  chkerr(VecDestroy(&eq_data_->ret_vec[i]));
368  }
369  delete[] eq_data_->ls;
370  delete[] eq_data_->ls_dt;
371  //delete[] stiffness_matrix;
372  //delete[] mass_matrix;
373  //delete[] rhs;
374  //delete[] mass_vec;
375  //delete[] ret_vec;
376 
377  delete mass_assembly_;
378  delete stiffness_assembly_;
379  delete sources_assembly_;
380  delete bdr_cond_assembly_;
381  delete init_assembly_;
382  }
383 
384 }
385 
386 
387 template<class Model>
389 {
390  START_TIMER(Model::ModelEqData::name());
391  eq_fields_->mark_input_times( *(Model::time_) );
392  eq_fields_->set_time(Model::time_->step(), LimitSide::left);
393  std::stringstream ss; // print warning message with table of uninitialized fields
394  if ( FieldCommon::print_message_table(ss, "transport DG") ) {
395  WarningOut() << ss.str();
396  }
397 
398 
399  // set initial conditions
400  set_initial_condition();
401  for (unsigned int sbi = 0; sbi < eq_data_->n_substances(); sbi++)
402  ( (LinSys_PETSC *)eq_data_->ls[sbi] )->set_initial_guess_nonzero();
403 
404  // check first time assembly - needs preallocation
405  if (!allocation_done) preallocate();
406 
407  // after preallocation we assemble the matrices and vectors required for mass balance
408  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
409  {
410  Model::balance_->calculate_instant(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
411 
412  // add sources due to sorption
413  ret_sources_prev[sbi] = 0;
414  }
415 
416  output_data();
417 }
418 
419 
420 template<class Model>
422 {
423  // preallocate system matrix
424  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
425  {
426  // preallocate system matrix
427  eq_data_->ls[i]->start_allocation();
428  stiffness_matrix[i] = NULL;
429  rhs[i] = NULL;
430 
431  // preallocate mass matrix
432  eq_data_->ls_dt[i]->start_allocation();
433  mass_matrix[i] = NULL;
434  VecZeroEntries(eq_data_->ret_vec[i]);
435  }
436  stiffness_assembly_->assemble(eq_data_->dh_);
437  mass_assembly_->assemble(eq_data_->dh_);
438  sources_assembly_->assemble(eq_data_->dh_);
439  bdr_cond_assembly_->assemble(eq_data_->dh_);
440  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
441  {
442  VecAssemblyBegin(eq_data_->ret_vec[i]);
443  VecAssemblyEnd(eq_data_->ret_vec[i]);
444  }
445 
446  allocation_done = true;
447 }
448 
449 
450 
451 template<class Model>
453 {
454  START_TIMER("DG-ONE STEP");
455 
456  Model::time_->next_time();
457  Model::time_->view("TDG");
458 
459  START_TIMER("data reinit");
460  eq_fields_->set_time(Model::time_->step(), LimitSide::left);
461  END_TIMER("data reinit");
462 
463  // assemble mass matrix
464  if (mass_matrix[0] == NULL || eq_fields_->subset(FieldFlag::in_time_term).changed() )
465  {
466  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
467  {
468  eq_data_->ls_dt[i]->start_add_assembly();
469  eq_data_->ls_dt[i]->mat_zero_entries();
470  VecZeroEntries(eq_data_->ret_vec[i]);
471  }
472  mass_assembly_->assemble(eq_data_->dh_);
473  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
474  {
475  eq_data_->ls_dt[i]->finish_assembly();
476  VecAssemblyBegin(eq_data_->ret_vec[i]);
477  VecAssemblyEnd(eq_data_->ret_vec[i]);
478  // construct mass_vec for initial time
479  if (mass_matrix[i] == NULL)
480  {
481  VecDuplicate(eq_data_->ls[i]->get_solution(), &mass_vec[i]);
482  MatMult(*(eq_data_->ls_dt[i]->get_matrix()), eq_data_->ls[i]->get_solution(), mass_vec[i]);
483  MatConvert(*( eq_data_->ls_dt[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &mass_matrix[i]);
484  }
485  else
486  MatCopy(*( eq_data_->ls_dt[i]->get_matrix() ), mass_matrix[i], DIFFERENT_NONZERO_PATTERN);
487  }
488  }
489 
490  // assemble stiffness matrix
491  if (stiffness_matrix[0] == NULL
492  || eq_fields_->subset(FieldFlag::in_main_matrix).changed()
493  || eq_fields_->flow_flux.changed())
494  {
495  // new fluxes can change the location of Neumann boundary,
496  // thus stiffness matrix must be reassembled
497  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
498  {
499  eq_data_->ls[i]->start_add_assembly();
500  eq_data_->ls[i]->mat_zero_entries();
501  }
502  stiffness_assembly_->assemble(eq_data_->dh_);
503  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
504  {
505  eq_data_->ls[i]->finish_assembly();
506 
507  if (stiffness_matrix[i] == NULL)
508  MatConvert(*( eq_data_->ls[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &stiffness_matrix[i]);
509  else
510  MatCopy(*( eq_data_->ls[i]->get_matrix() ), stiffness_matrix[i], DIFFERENT_NONZERO_PATTERN);
511  }
512  }
513 
514  // assemble right hand side (due to sources and boundary conditions)
515  if (rhs[0] == NULL
516  || eq_fields_->subset(FieldFlag::in_rhs).changed()
517  || eq_fields_->flow_flux.changed())
518  {
519  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
520  {
521  eq_data_->ls[i]->start_add_assembly();
522  eq_data_->ls[i]->rhs_zero_entries();
523  }
524  sources_assembly_->assemble(eq_data_->dh_);
525  bdr_cond_assembly_->assemble(eq_data_->dh_);
526  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
527  {
528  eq_data_->ls[i]->finish_assembly();
529 
530  if (rhs[i] == nullptr) VecDuplicate(*( eq_data_->ls[i]->get_rhs() ), &rhs[i]);
531  VecCopy(*( eq_data_->ls[i]->get_rhs() ), rhs[i]);
532  }
533  }
534 
535  /* Apply backward Euler time integration.
536  *
537  * Denoting A the stiffness matrix and M the mass matrix, the algebraic system at the k-th time level reads
538  *
539  * (1/dt M + A)u^k = f + 1/dt M.u^{k-1}
540  *
541  * Hence we modify at each time level the right hand side:
542  *
543  * f^k = f + 1/dt M u^{k-1},
544  *
545  * where f stands for the term stemming from the force and boundary conditions.
546  * Accordingly, we set
547  *
548  * A^k = A + 1/dt M.
549  *
550  */
551  Mat m;
552  START_TIMER("solve");
553  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
554  {
555  MatConvert(stiffness_matrix[i], MATSAME, MAT_INITIAL_MATRIX, &m);
556  MatAXPY(m, 1./Model::time_->dt(), mass_matrix[i], SUBSET_NONZERO_PATTERN);
557  eq_data_->ls[i]->set_matrix(m, DIFFERENT_NONZERO_PATTERN);
558  Vec w;
559  VecDuplicate(rhs[i], &w);
560  VecWAXPY(w, 1./Model::time_->dt(), mass_vec[i], rhs[i]);
561  eq_data_->ls[i]->set_rhs(w);
562 
563  VecDestroy(&w);
564  chkerr(MatDestroy(&m));
565 
566  eq_data_->ls[i]->solve();
567 
568  // update mass_vec due to possible changes in mass matrix
569  MatMult(*(eq_data_->ls_dt[i]->get_matrix()), eq_data_->ls[i]->get_solution(), mass_vec[i]);
570  }
571  END_TIMER("solve");
572 
573  // Possibly output matrices for debug reasons.
574  // for (unsigned int i=0; i<eq_data_->n_substances(); i++){
575  // string conc_name = eq_data_->substances().names()[i] + "_" + std::to_string(eq_data_->time_->step().index());
576  // eq_data_->ls[i]->view("stiff_" + conc_name);
577  // eq_data_->ls_dt[i]->view("mass_" + conc_name);
578  // }
579 
580  calculate_cumulative_balance();
581 
582  END_TIMER("DG-ONE STEP");
583 }
584 
585 
586 template<class Model>
588 {
589  // calculate element averages of solution
590  for (auto cell : eq_data_->dh_->own_range() )
591  {
592  LocDofVec loc_dof_indices = cell.get_loc_dof_indices();
593  unsigned int n_dofs=loc_dof_indices.n_rows;
594 
595  DHCellAccessor dh_p0_cell = eq_data_->dh_p0->cell_accessor_from_element(cell.elm_idx());
596  IntIdx dof_p0 = dh_p0_cell.get_loc_dof_indices()[0];
597 
598  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
599  {
600  eq_data_->conc_fe[sbi]->vec().set(dof_p0, 0);
601 
602  for (unsigned int j=0; j<n_dofs; ++j)
603  eq_data_->conc_fe[sbi]->vec().add( dof_p0, eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]] );
604 
605  eq_data_->conc_fe[sbi]->vec().set( dof_p0, max(eq_data_->conc_fe[sbi]->vec().get(dof_p0)/n_dofs, 0.) );
606  }
607  }
608 }
609 
610 
611 
612 
613 template<class Model>
615 {
616  //if (!Model::time_->is_current( Model::time_->marks().type_output() )) return;
617 
618 
619  START_TIMER("DG-OUTPUT");
620 
621  // gather the solution from all processors
622  eq_fields_->output_fields.set_time( this->time().step(), LimitSide::left);
623  //if (eq_fields_->output_fields.is_field_output_time(eq_fields_->output_field, this->time().step()) )
624  eq_fields_->output_fields.output(this->time().step());
625 
626  Model::output_data();
627 
628  START_TIMER("TOS-balance");
629  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
630  Model::balance_->calculate_instant(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
631  Model::balance_->output();
632  END_TIMER("TOS-balance");
633 
634  END_TIMER("DG-OUTPUT");
635 }
636 
637 
638 template<class Model>
640 {
641  if (Model::balance_->cumulative())
642  {
643  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
644  {
645  Model::balance_->calculate_cumulative(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
646 
647  // update source increment due to retardation
648  VecDot(eq_data_->ret_vec[sbi], eq_data_->ls[sbi]->get_solution(), &ret_sources[sbi]);
649 
650  Model::balance_->add_cumulative_source(eq_data_->subst_idx_[sbi], (ret_sources[sbi]-ret_sources_prev[sbi])/Model::time_->dt());
651  ret_sources_prev[sbi] = ret_sources[sbi];
652  }
653  }
654 }
655 
656 
657 
658 
659 template<class Model>
661 {
662  if(init_projection)
663  {
664  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
665  eq_data_->ls[sbi]->start_allocation();
666 
667  init_assembly_->assemble(eq_data_->dh_);
668 
669  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
670  eq_data_->ls[sbi]->start_add_assembly();
671 
672  init_assembly_->assemble(eq_data_->dh_);
673 
674  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
675  {
676  eq_data_->ls[sbi]->finish_assembly();
677  eq_data_->ls[sbi]->solve();
678  }
679  }
680  else
681  init_assembly_->assemble(eq_data_->dh_);
682 }
683 
684 
685 template<class Model>
687 {
688  el_4_loc = Model::mesh_->get_el_4_loc();
689  el_ds = Model::mesh_->get_el_ds();
690 }
691 
692 
693 template<class Model>
695 {
696  if (solution_changed)
697  {
698  for (auto cell : eq_data_->dh_->own_range() )
699  {
700  LocDofVec loc_dof_indices = cell.get_loc_dof_indices();
701  unsigned int n_dofs=loc_dof_indices.n_rows;
702 
703  DHCellAccessor dh_p0_cell = eq_data_->dh_p0->cell_accessor_from_element(cell.elm_idx());
704  IntIdx dof_p0 = dh_p0_cell.get_loc_dof_indices()[0];
705 
706  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
707  {
708  double old_average = 0;
709  for (unsigned int j=0; j<n_dofs; ++j)
710  old_average += eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]];
711  old_average /= n_dofs;
712 
713  for (unsigned int j=0; j<n_dofs; ++j)
714  eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]]
715  += eq_data_->conc_fe[sbi]->vec().get(dof_p0) - old_average;
716  }
717  }
718  }
719  // update mass_vec for the case that mass matrix changes in next time step
720  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
721  MatMult(*(eq_data_->ls_dt[sbi]->get_matrix()), eq_data_->ls[sbi]->get_solution(), mass_vec[sbi]);
722 }
723 
724 template<class Model>
726 {
727  return Model::mesh_->get_row_4_el();
728 }
729 
730 
731 
732 
733 
734 
736 template class TransportDG<HeatTransferModel>;
737 
738 
739 
740 
FieldCommon::units
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
Definition: field_common.hh:150
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
LocDofVec
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
Input::Type::Bool
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
fe_rt.hh
Definitions of Raviart-Thomas finite elements.
UnitSI::dimensionless
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
TransportDG::EqFields::dg_penalty
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
Definition: transport_dg.hh:153
factory.hh
LinSys
Definition: la_linsys_new.hh:169
transport_dg.hh
Discontinuous Galerkin method for equation of transport with dispersion.
TransportDG::input_rec
Input::Record input_rec
Array for storing the output solution data.
Definition: transport_dg.hh:396
Input::Type::Integer
Class for declaration of the integral input data.
Definition: type_base.hh:483
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:686
IntIdx
int IntIdx
Definition: index_types.hh:25
EquationBase::time
TimeGovernor & time()
Definition: equation.hh:149
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,...
value
static constexpr bool value
Definition: json.hpp:87
chkerr
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142
FLOW123D_FORCE_LINK_IN_CHILD
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:157
LinSys_PETSC
Definition: linsys_PETSC.hh:43
TransportDG::EqFields::EqFields
EqFields()
Definition: transport_dg.cc:111
Side::dim
unsigned int dim() const
Returns dimension of the side, that is dimension of the element minus one.
Definition: accessors_impl.hh:213
std::vector
Definition: doxy_dummy_defs.hh:7
TransportDG::update_solution
void update_solution() override
Computes the solution in one time instant.
Definition: transport_dg.cc:452
ElementAccessor< 3 >
arma::vec3
Definition: doxy_dummy_defs.hh:17
FieldCommon::flags
FieldCommon & flags(FieldFlag::Flags::Mask mask)
Definition: field_common.hh:189
advection_diffusion_model.hh
Discontinuous Galerkin method for equation of transport with dispersion.
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()
Calculates the dispersivity (diffusivity) tensor from the velocity field.
Definition: transport_dg.cc:660
TransportDG::EqData::set_DG_parameters_boundary
void set_DG_parameters_boundary(Side 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:166
index_types.hh
TransportDG::calculate_cumulative_balance
void calculate_cumulative_balance()
Definition: transport_dg.cc:639
TransportDG::compute_p0_interpolation
void compute_p0_interpolation()
Compute P0 interpolation of the solution (used in reaction term).
Definition: transport_dg.cc:587
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
TransportDG::EqFields::region_id
Field< 3, FieldValue< 3 >::Scalar > region_id
Definition: transport_dg.hh:154
TransportDG::TransportDG
TransportDG(Mesh &init_mesh, const Input::Record in_rec)
Constructor.
Definition: transport_dg.cc:202
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:195
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
Distribution
Definition: distribution.hh:50
EquationBase::balance
std::shared_ptr< Balance > balance() const
Definition: equation.hh:185
TransportDG::output_data
void output_data()
Postprocesses the solution and writes to output file.
Definition: transport_dg.cc:614
LimitSide::left
@ left
dh_cell_accessor.hh
accessors.hh
TransportDG::eq_fields_
std::shared_ptr< EqFields > eq_fields_
Fields for model parameters.
Definition: transport_dg.hh:363
Input::Record
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
sys_profiler.hh
generic_assembly.hh
quadrature_lib.hh
Definitions of particular quadrature rules on simplices.
TransportDG::get_row_4_el
LongIdx * get_row_4_el()
Definition: transport_dg.cc:725
TransportDG::initialize
void initialize() override
Definition: transport_dg.cc:242
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
TransportDG::preallocate
void preallocate()
Definition: transport_dg.cc:421
Side
Definition: accessors.hh:361
GenericField::subdomain
static auto subdomain(Mesh &mesh) -> IndexField
Definition: generic_field.impl.hh:58
TransportOperatorSplitting::output_data
void output_data() override
Write computed fields.
Definition: transport_operator_splitting.cc:248
Input::Type::Default::obligatory
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
FieldCommon::print_message_table
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:96
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:503
Side::element
ElementAccessor< 3 > element() const
Returns iterator to the element of the side.
Definition: accessors_impl.hh:233
TransportDG::allocation_done
bool allocation_done
Indicates whether matrices have been preallocated.
Definition: transport_dg.hh:419
TransportDG::EqFields
Definition: transport_dg.hh:147
Input::Type::Selection
Template for classes storing finite set of named values.
Definition: type_selection.hh:65
Side::n_nodes
unsigned int n_nodes() const
Returns number of nodes of the side.
Definition: accessors.hh:407
ElementAccessor::node
NodeAccessor< 3 > node(unsigned int ni) const
Definition: accessors.hh:200
Input::Type::Record::close
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
heat_model.hh
Discontinuous Galerkin method for equation of transport with dispersion.
generic_field.hh
Input::Type
Definition: balance.hh:41
Input::Type::Record
Record type proxy class.
Definition: type_record.hh:182
TransportDG
Transport with dispersion implemented using discontinuous Galerkin method.
Definition: transport_dg.hh:134
TransportDG::EqFields::fracture_sigma
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
Definition: transport_dg.hh:152
GenericField::region_id
static auto region_id(Mesh &mesh) -> IndexField
Definition: generic_field.impl.hh:39
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:137
LongIdx
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
FieldFE
Definition: field.hh:63
TransportDG::eq_data_
std::shared_ptr< EqData > eq_data_
Data for model parameters.
Definition: transport_dg.hh:366
TransportDG::DGVariant
DGVariant
Definition: transport_dg.hh:244
Mesh
Definition: mesh.h:77
multi_field.hh
Input::Type::Array
Class for declaration of inputs sequences.
Definition: type_base.hh:339
TransportDG::update_after_reactions
void update_after_reactions(bool solution_changed)
Definition: transport_dg.cc:694
Model
Definition: field_model.hh:304
DHCellAccessor
Cell accessor allow iterate over DOF handler cells.
Definition: dh_cell_accessor.hh:43
Input::Array
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
WarningOut
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:270
fe_value_handler.hh
TransportDG::EqFields::subdomain
Field< 3, FieldValue< 3 >::Scalar > subdomain
Definition: transport_dg.hh:155
OutputTime::CORNER_DATA
@ CORNER_DATA
Definition: output_time.hh:110
MixedPtr< FE_P_disc >
DHCellAccessor::get_loc_dof_indices
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
Definition: dh_cell_accessor.hh:88
assembly_dg.hh
concentration_model.hh
Discontinuous Galerkin method for equation of transport with dispersion.
EquationBase::eq_fieldset_
FieldSet * eq_fieldset_
Definition: equation.hh:227
TransportDG::zero_time_step
void zero_time_step() override
Initialize solution in the zero time.
Definition: transport_dg.cc:388
TransportDG::~TransportDG
~TransportDG() override
Destructor.
Definition: transport_dg.cc:346
balance.hh
FieldCommon::description
FieldCommon & description(const string &description)
Definition: field_common.hh:125
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
Returns node for given local index i on the side.
Definition: accessors_impl.hh:227
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
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
START_TIMER
#define START_TIMER(tag)
Starts a timer with specified tag.
Definition: sys_profiler.hh:115
TransportDG::EqData::elem_anisotropy
double elem_anisotropy(ElementAccessor< 3 > e) const
Compute and return anisotropy of given element.
Definition: transport_dg.cc:150
GenericAssembly
Generic class of assemblation.
Definition: generic_assembly.hh:70
TransportDG::EqFields::output_fields
EquationOutput output_fields
Definition: transport_dg.hh:157
LinSys_PETSC::get_input_type
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
END_TIMER
#define END_TIMER(tag)
Ends a timer with specified tag.
Definition: sys_profiler.hh:149
Element::n_nodes
unsigned int n_nodes() const
Definition: elements.h:126
FieldCommon::name
FieldCommon & name(const string &name)
Definition: field_common.hh:118
FieldValue
Definition: field_values.hh:645