Flow123d  master-f44eb46
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"
30 #include "la/linsys_PETSC.hh"
31 #include "coupling/balance.hh"
35 #include "transport/heat_model.hh"
36 #include "transport/assembly_dg.hh"
37 
38 #include "fields/multi_field.hh"
39 #include "fields/generic_field.hh"
40 #include "input/factory.hh"
42 #include "mesh/accessors.hh"
43 
44 FLOW123D_FORCE_LINK_IN_CHILD(concentrationTransportModel)
46 
47 
48 
49 using namespace Input::Type;
50 
51 template<class Model>
53  return Selection("DG_variant", "Type of penalty term.")
54  .add_value(non_symmetric, "non-symmetric", "non-symmetric weighted interior penalty DG method")
55  .add_value(incomplete, "incomplete", "incomplete weighted interior penalty DG method")
56  .add_value(symmetric, "symmetric", "symmetric weighted interior penalty DG method")
57  .close();
58 }
59 
60 /*
61 * Should be removed
62 template<class Model>
63 const Selection & TransportDG<Model>::EqData::get_output_selection() {
64  return Model::ModelEqData::get_output_selection_input_type(
65  "DG",
66  "Implicit in time Discontinuous Galerkin solver")
67  .copy_values(EqData().make_output_field_selection("").close())
68  ConvectionTransport::EqData().output_fields
69  .make_output_field_selection(
70  "ConvectionTransport_output_fields",
71  "Selection of output fields for Convection Solute Transport model.")
72  .close()),
73  .close();
74 }
75 */
76 
77 template<class Model>
79  std::string equation_name = std::string(Model::ModelEqData::name()) + "_DG";
80  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  this->add_coords_field();
145  this->set_default_fieldset();
146 }
147 
148 
149 
150 template<typename Model>
152  : Model(init_mesh, in_rec),
153  input_rec(in_rec),
154  allocation_done(false),
155  mass_assembly_(nullptr)
156 {
157  // Can not use name() + "constructor" here, since START_TIMER only accepts const char *
158  // due to constexpr optimization.
159  START_TIMER(Model::ModelEqData::name());
160  // Check that Model is derived from AdvectionDiffusionModel.
162 
163  eq_data_ = make_shared<EqData>();
164  eq_fields_ = make_shared<EqFields>();
165  this->eq_fieldset_ = eq_fields_;
166  Model::init_balance(in_rec);
167 
168 
169  // Set up physical parameters.
170  eq_fields_->set_mesh(init_mesh);
171  eq_fields_->region_id = GenericField<3>::region_id(*Model::mesh_);
172  eq_fields_->subdomain = GenericField<3>::subdomain(*Model::mesh_);
173 
174 
175  // DG data parameters
176  eq_data_->dg_variant = in_rec.val<DGVariant>("dg_variant");
177  eq_data_->dg_order = in_rec.val<unsigned int>("dg_order");
178 
179  Model::init_from_input(in_rec);
180 
181  MixedPtr<FE_P_disc> fe(eq_data_->dg_order);
182  shared_ptr<DiscreteSpace> ds = make_shared<EqualOrderDiscreteSpace>(Model::mesh_, fe);
183  eq_data_->dh_ = make_shared<DOFHandlerMultiDim>(*Model::mesh_);
184  eq_data_->dh_->distribute_dofs(ds);
185  //DebugOut().fmt("TDG: solution size {}\n", eq_data_->dh_->n_global_dofs());
186 
187 }
188 
189 
190 template<class Model>
192 {
193  eq_fields_->set_components(eq_data_->substances_.names());
194  eq_fields_->set_input_list( input_rec.val<Input::Array>("input_fields"), *(Model::time_) );
195  eq_data_->set_time_governor(Model::time_);
196  eq_data_->balance_ = this->balance();
197  eq_fields_->initialize();
198 
199  // Resize coefficient arrays
200  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)));
201  ret_sources.resize(eq_data_->n_substances());
202  ret_sources_prev.resize(eq_data_->n_substances());
203 
204  Input::Array user_fields_arr;
205  if (input_rec.opt_val("user_fields", user_fields_arr)) {
206  this->init_user_fields(user_fields_arr, eq_fields_->output_fields);
207  }
208 
209  eq_data_->output_vec.resize(eq_data_->n_substances());
210  eq_fields_->output_field.set_components(eq_data_->substances_.names());
211  eq_fields_->output_field.set_mesh(*Model::mesh_);
212  eq_fields_->output_fields.set_mesh(*Model::mesh_);
213  eq_fields_->output_type(OutputTime::CORNER_DATA);
214 
215  eq_fields_->output_field.setup_components();
216  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
217  {
218  // create shared pointer to a FieldFE, pass FE data and push this FieldFE to output_field on all regions
219  auto output_field_ptr = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_);
220  eq_fields_->output_field[sbi].set(output_field_ptr, 0);
221  eq_data_->output_vec[sbi] = output_field_ptr->vec();
222  }
223 
224  // set time marks for writing the output
225  eq_fields_->output_fields.initialize(Model::output_stream_, Model::mesh_, input_rec.val<Input::Record>("output"), this->time());
226 
227  // equation default PETSc solver options
228  std::string petsc_default_opts;
229  if (eq_data_->dh_->distr()->np() == 1)
230  petsc_default_opts = "-ksp_type bcgs -pc_type ilu -pc_factor_levels 2 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
231  else
232  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";
233 
234  // allocate matrix and vector structures
235  eq_data_->ls = new LinSys*[eq_data_->n_substances()];
236  eq_data_->ls_dt = new LinSys*[eq_data_->n_substances()];
237  eq_data_->conc_fe.resize(eq_data_->n_substances());
238 
239  MixedPtr<FE_P_disc> fe(0);
240  shared_ptr<DiscreteSpace> ds = make_shared<EqualOrderDiscreteSpace>(Model::mesh_, fe);
241  eq_data_->dh_p0 = make_shared<DOFHandlerMultiDim>(*Model::mesh_);
242  eq_data_->dh_p0->distribute_dofs(ds);
243 
244  stiffness_matrix.resize(eq_data_->n_substances(), nullptr);
245  mass_matrix.resize(eq_data_->n_substances(), nullptr);
246  rhs.resize(eq_data_->n_substances(), nullptr);
247  mass_vec.resize(eq_data_->n_substances(), nullptr);
248  eq_data_->ret_vec.resize(eq_data_->n_substances(), nullptr);
249 
250  for (unsigned int sbi = 0; sbi < eq_data_->n_substances(); sbi++) {
251  eq_data_->ls[sbi] = new LinSys_PETSC(eq_data_->dh_->distr().get(), petsc_default_opts);
252  ( (LinSys_PETSC *)eq_data_->ls[sbi] )->set_from_input( input_rec.val<Input::Record>("solver") );
253  eq_data_->ls[sbi]->set_solution(eq_data_->output_vec[sbi].petsc_vec());
254 
255  eq_data_->ls_dt[sbi] = new LinSys_PETSC(eq_data_->dh_->distr().get(), petsc_default_opts);
256  ( (LinSys_PETSC *)eq_data_->ls_dt[sbi] )->set_from_input( input_rec.val<Input::Record>("solver") );
257 
258  eq_data_->conc_fe[sbi] = create_field_fe< 3, FieldValue<3>::Scalar >(eq_data_->dh_p0);
259 
260  VecDuplicate(eq_data_->ls[sbi]->get_solution(), &eq_data_->ret_vec[sbi]);
261  }
262 
263 
264  init_projection = input_rec.val<bool>("init_projection");
265 
266  // create assemblation object, finite element structures and distribute DOFs
267  mass_assembly_ = new GenericAssembly< MassAssemblyDim >(eq_fields_.get(), eq_data_.get());
268  stiffness_assembly_ = new GenericAssembly< StiffnessAssemblyDim >(eq_fields_.get(), eq_data_.get());
269  sources_assembly_ = new GenericAssembly< SourcesAssemblyDim >(eq_fields_.get(), eq_data_.get());
270  bdr_cond_assembly_ = new GenericAssembly< BdrConditionAssemblyDim >(eq_fields_.get(), eq_data_.get());
271 
272  if(init_projection)
273  init_assembly_ = new GenericAssembly< InitProjectionAssemblyDim >(eq_fields_.get(), eq_data_.get());
274  else
275  init_assembly_ = new GenericAssembly< InitConditionAssemblyDim >(eq_fields_.get(), eq_data_.get());
276 
277  // initialization of balance object
278  Model::balance_->allocate(eq_data_->dh_->distr()->lsize(), mass_assembly_->eval_points()->max_size());
279 
280  eq_fields_->init_condition.setup_components();
281  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
282  {
283  eq_fields_->init_condition[sbi].add_factory( std::make_shared<FieldFE<3, FieldValue<3>::Scalar>::NativeFactory>(sbi, eq_data_->dh_));
284  }
285 }
286 
287 
288 template<class Model>
290 {
291  delete Model::time_;
292 
293  if (rhs.size() > 0) {
294  // initialize called
295 
296  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
297  {
298  if (eq_data_->ls != nullptr) {
299  delete eq_data_->ls[i];
300  delete eq_data_->ls_dt[i];
301  }
302 
303  if (stiffness_matrix.size() > 0) {
304  if (stiffness_matrix[i])
305  chkerr(MatDestroy(&stiffness_matrix[i]));
306  if (mass_matrix[i])
307  chkerr(MatDestroy(&mass_matrix[i]));
308  if (rhs[i])
309  chkerr(VecDestroy(&rhs[i]));
310  if (mass_vec[i])
311  chkerr(VecDestroy(&mass_vec[i]));
312  if (eq_data_->ret_vec[i])
313  chkerr(VecDestroy(&eq_data_->ret_vec[i]));
314  }
315  }
316  if (eq_data_->ls != nullptr) {
317  delete[] eq_data_->ls;
318  delete[] eq_data_->ls_dt;
319  eq_data_->ls = nullptr;
320  }
321  //delete[] stiffness_matrix;
322  //delete[] mass_matrix;
323  //delete[] rhs;
324  //delete[] mass_vec;
325  //delete[] ret_vec;
326 
327  if (mass_assembly_ != nullptr) {
328  delete mass_assembly_;
329  delete stiffness_assembly_;
330  delete sources_assembly_;
331  delete bdr_cond_assembly_;
332  delete init_assembly_;
333  }
334  }
335 
336 
337 }
338 
339 
340 template<class Model>
342 {
343  START_TIMER(Model::ModelEqData::name());
344  eq_fields_->mark_input_times( *(Model::time_) );
345  eq_fields_->set_time(Model::time_->step(), LimitSide::left);
346  std::stringstream ss; // print warning message with table of uninitialized fields
347  if ( FieldCommon::print_message_table(ss, "transport DG") ) {
348  WarningOut() << ss.str();
349  }
350 
351 
352  // set initial conditions
353  set_initial_condition();
354  for (unsigned int sbi = 0; sbi < eq_data_->n_substances(); sbi++)
355  ( (LinSys_PETSC *)eq_data_->ls[sbi] )->set_initial_guess_nonzero();
356 
357  // check first time assembly - needs preallocation
358  if (!allocation_done) preallocate();
359 
360  // after preallocation we assemble the matrices and vectors required for mass balance
361  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
362  {
363  Model::balance_->calculate_instant(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
364 
365  // add sources due to sorption
366  ret_sources_prev[sbi] = 0;
367  }
368 
369  output_data();
370 }
371 
372 
373 template<class Model>
375 {
376  // preallocate system matrix
377  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
378  {
379  // preallocate system matrix
380  eq_data_->ls[i]->start_allocation();
381  stiffness_matrix[i] = NULL;
382  rhs[i] = NULL;
383 
384  // preallocate mass matrix
385  eq_data_->ls_dt[i]->start_allocation();
386  mass_matrix[i] = NULL;
387  VecZeroEntries(eq_data_->ret_vec[i]);
388  }
389  stiffness_assembly_->assemble(eq_data_->dh_);
390  mass_assembly_->assemble(eq_data_->dh_);
391  sources_assembly_->assemble(eq_data_->dh_);
392  bdr_cond_assembly_->assemble(eq_data_->dh_);
393  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
394  {
395  VecAssemblyBegin(eq_data_->ret_vec[i]);
396  VecAssemblyEnd(eq_data_->ret_vec[i]);
397  }
398 
399  allocation_done = true;
400 }
401 
402 
403 
404 template<class Model>
406 {
407  START_TIMER("DG-ONE STEP");
408 
409  Model::time_->next_time();
410  Model::time_->view("TDG");
411 
412  START_TIMER("data reinit");
413  eq_fields_->set_time(Model::time_->step(), LimitSide::left);
414  END_TIMER("data reinit");
415 
416  // assemble mass matrix
417  if (mass_matrix[0] == NULL || eq_fields_->subset(FieldFlag::in_time_term).changed() )
418  {
419  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
420  {
421  eq_data_->ls_dt[i]->start_add_assembly();
422  eq_data_->ls_dt[i]->mat_zero_entries();
423  VecZeroEntries(eq_data_->ret_vec[i]);
424  }
425  mass_assembly_->assemble(eq_data_->dh_);
426  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
427  {
428  eq_data_->ls_dt[i]->finish_assembly();
429  VecAssemblyBegin(eq_data_->ret_vec[i]);
430  VecAssemblyEnd(eq_data_->ret_vec[i]);
431  // construct mass_vec for initial time
432  if (mass_matrix[i] == NULL)
433  {
434  VecDuplicate(eq_data_->ls[i]->get_solution(), &mass_vec[i]);
435  MatMult(*(eq_data_->ls_dt[i]->get_matrix()), eq_data_->ls[i]->get_solution(), mass_vec[i]);
436  MatConvert(*( eq_data_->ls_dt[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &mass_matrix[i]);
437  }
438  else
439  MatCopy(*( eq_data_->ls_dt[i]->get_matrix() ), mass_matrix[i], DIFFERENT_NONZERO_PATTERN);
440  }
441  }
442 
443  // assemble stiffness matrix
444  if (stiffness_matrix[0] == NULL
445  || eq_fields_->subset(FieldFlag::in_main_matrix).changed()
446  || eq_fields_->flow_flux.changed())
447  {
448  // new fluxes can change the location of Neumann boundary,
449  // thus stiffness matrix must be reassembled
450  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
451  {
452  eq_data_->ls[i]->start_add_assembly();
453  eq_data_->ls[i]->mat_zero_entries();
454  }
455  stiffness_assembly_->assemble(eq_data_->dh_);
456  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
457  {
458  eq_data_->ls[i]->finish_assembly();
459 
460  if (stiffness_matrix[i] == NULL)
461  MatConvert(*( eq_data_->ls[i]->get_matrix() ), MATSAME, MAT_INITIAL_MATRIX, &stiffness_matrix[i]);
462  else
463  MatCopy(*( eq_data_->ls[i]->get_matrix() ), stiffness_matrix[i], DIFFERENT_NONZERO_PATTERN);
464  }
465  }
466 
467  // assemble right hand side (due to sources and boundary conditions)
468  if (rhs[0] == NULL
469  || eq_fields_->subset(FieldFlag::in_rhs).changed()
470  || eq_fields_->flow_flux.changed())
471  {
472  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
473  {
474  eq_data_->ls[i]->start_add_assembly();
475  eq_data_->ls[i]->rhs_zero_entries();
476  }
477  sources_assembly_->assemble(eq_data_->dh_);
478  bdr_cond_assembly_->assemble(eq_data_->dh_);
479  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
480  {
481  eq_data_->ls[i]->finish_assembly();
482 
483  if (rhs[i] == nullptr) VecDuplicate(*( eq_data_->ls[i]->get_rhs() ), &rhs[i]);
484  VecCopy(*( eq_data_->ls[i]->get_rhs() ), rhs[i]);
485  }
486  }
487 
488  /* Apply backward Euler time integration.
489  *
490  * Denoting A the stiffness matrix and M the mass matrix, the algebraic system at the k-th time level reads
491  *
492  * (1/dt M + A)u^k = f + 1/dt M.u^{k-1}
493  *
494  * Hence we modify at each time level the right hand side:
495  *
496  * f^k = f + 1/dt M u^{k-1},
497  *
498  * where f stands for the term stemming from the force and boundary conditions.
499  * Accordingly, we set
500  *
501  * A^k = A + 1/dt M.
502  *
503  */
504  Mat m;
505  START_TIMER("solve");
506  for (unsigned int i=0; i<eq_data_->n_substances(); i++)
507  {
508  MatConvert(stiffness_matrix[i], MATSAME, MAT_INITIAL_MATRIX, &m);
509  MatAXPY(m, 1./Model::time_->dt(), mass_matrix[i], SUBSET_NONZERO_PATTERN);
510  eq_data_->ls[i]->set_matrix(m, DIFFERENT_NONZERO_PATTERN);
511  Vec w;
512  VecDuplicate(rhs[i], &w);
513  VecWAXPY(w, 1./Model::time_->dt(), mass_vec[i], rhs[i]);
514  eq_data_->ls[i]->set_rhs(w);
515 
516  VecDestroy(&w);
517  chkerr(MatDestroy(&m));
518 
519  eq_data_->ls[i]->solve();
520 
521  // update mass_vec due to possible changes in mass matrix
522  MatMult(*(eq_data_->ls_dt[i]->get_matrix()), eq_data_->ls[i]->get_solution(), mass_vec[i]);
523  }
524  END_TIMER("solve");
525 
526  // Possibly output matrices for debug reasons.
527  // for (unsigned int i=0; i<eq_data_->n_substances(); i++){
528  // string conc_name = eq_data_->substances().names()[i] + "_" + std::to_string(eq_data_->time_->step().index());
529  // eq_data_->ls[i]->view("stiff_" + conc_name);
530  // eq_data_->ls_dt[i]->view("mass_" + conc_name);
531  // }
532 
533  calculate_cumulative_balance();
534 
535  END_TIMER("DG-ONE STEP");
536 }
537 
538 
539 template<class Model>
541 {
542  // calculate element averages of solution
543  for (auto cell : eq_data_->dh_->own_range() )
544  {
545  LocDofVec loc_dof_indices = cell.get_loc_dof_indices();
546  unsigned int n_dofs=loc_dof_indices.n_rows;
547 
548  DHCellAccessor dh_p0_cell = eq_data_->dh_p0->cell_accessor_from_element(cell.elm_idx());
549  IntIdx dof_p0 = dh_p0_cell.get_loc_dof_indices()[0];
550 
551  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
552  {
553  eq_data_->conc_fe[sbi]->vec().set(dof_p0, 0);
554 
555  for (unsigned int j=0; j<n_dofs; ++j)
556  eq_data_->conc_fe[sbi]->vec().add( dof_p0, eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]] );
557 
558  eq_data_->conc_fe[sbi]->vec().set( dof_p0, max(eq_data_->conc_fe[sbi]->vec().get(dof_p0)/n_dofs, 0.) );
559  }
560  }
561 }
562 
563 
564 
565 
566 template<class Model>
568 {
569  //if (!Model::time_->is_current( Model::time_->marks().type_output() )) return;
570 
571 
572  START_TIMER("DG-OUTPUT");
573 
574  // gather the solution from all processors
575  eq_fields_->output_fields.set_time( this->time().step(), LimitSide::left);
576  //if (eq_fields_->output_fields.is_field_output_time(eq_fields_->output_field, this->time().step()) )
577  eq_fields_->output_fields.output(this->time().step());
578 
579  Model::output_data();
580 
581  START_TIMER("TOS-balance");
582  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
583  Model::balance_->calculate_instant(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
584  Model::balance_->output();
585  END_TIMER("TOS-balance");
586 
587  END_TIMER("DG-OUTPUT");
588 }
589 
590 
591 template<class Model>
593 {
594  if (Model::balance_->cumulative())
595  {
596  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
597  {
598  Model::balance_->calculate_cumulative(eq_data_->subst_idx_[sbi], eq_data_->ls[sbi]->get_solution());
599 
600  // update source increment due to retardation
601  VecDot(eq_data_->ret_vec[sbi], eq_data_->ls[sbi]->get_solution(), &ret_sources[sbi]);
602 
603  Model::balance_->add_cumulative_source(eq_data_->subst_idx_[sbi], (ret_sources[sbi]-ret_sources_prev[sbi])/Model::time_->dt());
604  ret_sources_prev[sbi] = ret_sources[sbi];
605  }
606  }
607 }
608 
609 
610 
611 
612 template<class Model>
614 {
615  if(init_projection)
616  {
617  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
618  eq_data_->ls[sbi]->start_allocation();
619 
620  init_assembly_->assemble(eq_data_->dh_);
621 
622  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
623  eq_data_->ls[sbi]->start_add_assembly();
624 
625  init_assembly_->assemble(eq_data_->dh_);
626 
627  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); sbi++)
628  {
629  eq_data_->ls[sbi]->finish_assembly();
630  eq_data_->ls[sbi]->solve();
631  }
632  }
633  else
634  init_assembly_->assemble(eq_data_->dh_);
635 }
636 
637 
638 template<class Model>
640 {
641  if (solution_changed)
642  {
643  for (auto cell : eq_data_->dh_->own_range() )
644  {
645  LocDofVec loc_dof_indices = cell.get_loc_dof_indices();
646  unsigned int n_dofs=loc_dof_indices.n_rows;
647 
648  DHCellAccessor dh_p0_cell = eq_data_->dh_p0->cell_accessor_from_element(cell.elm_idx());
649  IntIdx dof_p0 = dh_p0_cell.get_loc_dof_indices()[0];
650 
651  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
652  {
653  double old_average = 0;
654  for (unsigned int j=0; j<n_dofs; ++j)
655  old_average += eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]];
656  old_average /= n_dofs;
657 
658  for (unsigned int j=0; j<n_dofs; ++j)
659  eq_data_->ls[sbi]->get_solution_array()[loc_dof_indices[j]]
660  += eq_data_->conc_fe[sbi]->vec().get(dof_p0) - old_average;
661  }
662  }
663  }
664  // update mass_vec for the case that mass matrix changes in next time step
665  for (unsigned int sbi=0; sbi<eq_data_->n_substances(); ++sbi)
666  MatMult(*(eq_data_->ls_dt[sbi]->get_matrix()), eq_data_->ls[sbi]->get_solution(), mass_vec[sbi]);
667 }
668 
669 
670 
671 
672 
673 
675 template class TransportDG<HeatTransferModel>;
676 
677 
678 
679 
Discontinuous Galerkin method for equation of transport with dispersion.
Cell accessor allow iterate over DOF handler cells.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
static Input::Type::Record & user_fields_template(std::string equation_name)
Template Record with common key user_fields for derived equations.
Definition: equation.cc:46
static bool print_message_table(ostream &stream, std::string equation_name)
Definition: field_common.cc:95
FieldCommon & description(const string &description)
FieldCommon & flags_add(FieldFlag::Flags::Mask mask)
FieldCommon & flags(FieldFlag::Flags::Mask mask)
FieldCommon & name(const string &name)
FieldCommon & units(const UnitSI &units)
Set basic units of the field.
FieldCommon & input_default(const string &input_default)
static constexpr Mask in_time_term
A field is part of time term of the equation.
Definition: field_flag.hh:47
static constexpr Mask equation_external_output
Match an output field, that can be also copy of other field.
Definition: field_flag.hh:58
static constexpr Mask in_rhs
A field is part of the right hand side of the equation.
Definition: field_flag.hh:51
static constexpr Mask in_main_matrix
A field is part of main "stiffness matrix" of the equation.
Definition: field_flag.hh:49
Generic class of assemblation.
static auto region_id(Mesh &mesh) -> IndexField
static auto subdomain(Mesh &mesh) -> IndexField
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
const Ret val(const string &key) const
Class for declaration of inputs sequences.
Definition: type_base.hh:339
Class for declaration of the input of type Bool.
Definition: type_base.hh:452
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
Class for declaration of the integral input data.
Definition: type_base.hh:483
Record type proxy class.
Definition: type_record.hh:182
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & copy_keys(const Record &other)
Copy keys from other record.
Definition: type_record.cc:216
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
Template for classes storing finite set of named values.
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.
const Selection & close() const
Close the Selection, no more values can be added.
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
Definition: mesh.h:362
Field< 3, FieldValue< 3 >::Scalar > region_id
MultiField< 3, FieldValue< 3 >::Scalar > dg_penalty
Penalty enforcing inter-element continuity of solution (for each substance).
EquationOutput output_fields
Field< 3, FieldValue< 3 >::Scalar > subdomain
MultiField< 3, FieldValue< 3 >::Scalar > fracture_sigma
Transition parameter for diffusive transfer on fractures (for each substance).
Transport with dispersion implemented using discontinuous Galerkin method.
void update_solution() override
Computes the solution in one time instant.
static const Input::Type::Record & get_input_type()
Declare input record type for the equation TransportDG.
Definition: transport_dg.cc:78
bool allocation_done
Indicates whether matrices have been preallocated.
Input::Record input_rec
Array for storing the output solution data.
void zero_time_step() override
Initialize solution in the zero time.
std::shared_ptr< EqFields > eq_fields_
Fields for model parameters.
void initialize() override
GenericAssembly< MassAssemblyDim > * mass_assembly_
general assembly objects, hold assembly objects of appropriate dimension
static const Input::Type::Selection & get_dg_variant_selection_input_type()
Input type for the DG variant selection.
Definition: transport_dg.cc:52
TransportDG(Mesh &init_mesh, const Input::Record in_rec)
Constructor.
~TransportDG() override
Destructor.
void compute_p0_interpolation()
Compute P0 interpolation of the solution (used in reaction term).
void output_data()
Postprocesses the solution and writes to output file.
void calculate_cumulative_balance()
void preallocate()
void set_initial_condition()
Calculates the dispersivity (diffusivity) tensor from the velocity field.
std::shared_ptr< EqData > eq_data_
Data for model parameters.
void update_after_reactions(bool solution_changed)
static UnitSI & dimensionless()
Returns dimensionless unit.
Definition: unit_si.cc:55
Discontinuous Galerkin method for equation of transport with dispersion.
Definitions of basic Lagrangean finite elements with polynomial shape functions.
Definitions of Raviart-Thomas finite elements.
Class FEValues calculates finite element data on the actual cells such as shape function values,...
#define FLOW123D_FORCE_LINK_IN_CHILD(x)
Definition: global_defs.h:104
Discontinuous Galerkin method for equation of transport with dispersion.
int IntIdx
Definition: index_types.hh:25
arma::Col< IntIdx > LocDofVec
Definition: index_types.hh:28
static constexpr bool value
Definition: json.hpp:87
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
Definitions of particular quadrature rules on simplices.
#define END_TIMER(tag)
Ends a timer with specified tag.
#define START_TIMER(tag)
Starts a timer with specified tag.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142
Discontinuous Galerkin method for equation of transport with dispersion.