Flow123d
transport_dg.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2007 Technical University of Liberec. All rights reserved.
4  *
5  * Please make a following refer to Flow123d on your project site if you use the program for any purpose,
6  * especially for academic research:
7  * Flow123d, Research Centre: Advanced Remedial Technologies, Technical University of Liberec, Czech Republic
8  *
9  * This program is free software; you can redistribute it and/or modify it under the terms
10  * of the GNU General Public License version 3 as published by the Free Software Foundation.
11  *
12  * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
13  * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14  * See the GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License along with this program; if not,
17  * write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 021110-1307, USA.
18  *
19  *
20  * $Id$
21  * $Revision$
22  * $LastChangedBy$
23  * $LastChangedDate$
24  *
25  * @file
26  * @brief Discontinuous Galerkin method for equation of transport with dispersion.
27  * @author Jan Stebel
28  */
29 
30 #include "petscmat.h"
31 #include <armadillo>
32 #include <boost/foreach.hpp>
33 #include "system/xio.h"
34 #include "system/sys_profiler.hh"
35 #include "transport/transport_dg.hh"
36 #include "quadrature/quadrature_lib.hh"
37 #include "fem/mapping_p1.hh"
38 #include "fem/fe_values.hh"
39 #include "fem/fe_p.hh"
40 #include "fem/fe_rt.hh"
41 #include "io/output.h"
42 #include "fields/field_fe.hh"
43 #include "mesh/boundaries.h"
44 #include "la/distribution.hh"
45 #include "input/accessors.hh"
46 #include "flow/darcy_flow_mh.hh"
47 #include "la/linsys_PETSC.hh"
48 #include "mesh/partitioning.hh"
49 #include "transport/advection_diffusion_model.hh"
50 #include "transport/concentration_model.hh"
51 #include "transport/heat_model.hh"
52 
53 using namespace Input::Type;
54 
55 template<class Model>
56 Selection TransportDG<Model>::dg_variant_selection_input_type
57  = Selection("DG_variant", "Type of penalty term.")
58  .add_value(non_symmetric, "non-symmetric", "non-symmetric weighted interior penalty DG method")
59  .add_value(incomplete, "incomplete", "incomplete weighted interior penalty DG method")
60  .add_value(symmetric, "symmetric", "symmetric weighted interior penalty DG method");
61 
62 template<class Model>
63 Selection TransportDG<Model>::EqData::bc_type_selection =
64  Selection("TransportDG_EqData_bc_Type")
65  .add_value(inflow, "inflow", "Dirichlet BC on inflow and homogeneous Neumann BC on outflow.")
66  .add_value(dirichlet, "dirichlet")
67  .add_value(neumann, "neumann")
68  .add_value(robin, "robin");
69 
70 template<class Model>
71 Selection TransportDG<Model>::EqData::output_selection =
72  Model::ModelEqData::get_output_selection_input_type("DG", "DG solver")
73  .copy_values(EqData().output_fields.make_output_field_selection("").close())
74  .close();
75 
76 template<class Model>
77 Record TransportDG<Model>::input_type
78  = Model::get_input_type("DG", "DG solver")
79  .declare_key("solver", LinSys_PETSC::input_type, Default::obligatory(),
80  "Linear solver for MH problem.")
81  .declare_key("input_fields", Array(TransportDG<Model>::EqData().make_field_descriptor_type(Model::ModelEqData::name() + "_DG")), IT::Default::obligatory(), "")
82  .declare_key("dg_variant", TransportDG<Model>::dg_variant_selection_input_type, Default("non-symmetric"),
83  "Variant of interior penalty discontinuous Galerkin method.")
84  .declare_key("dg_order", Integer(0,3), Default("1"),
85  "Polynomial order for finite element in DG method (order 0 is suitable if there is no diffusion/dispersion).")
86  .declare_key("output_fields", Array(EqData::output_selection),
87  Default(Model::ModelEqData::default_output_field()),
88  "List of fields to write to output file.");
89 
90 
91 
92 FEObjects::FEObjects(Mesh *mesh_, unsigned int fe_order)
93 {
94  unsigned int q_order;
95 
96  switch (fe_order)
97  {
98  case 0:
99  q_order = 0;
100  fe1_ = new FE_P_disc<0,1,3>;
101  fe2_ = new FE_P_disc<0,2,3>;
102  fe3_ = new FE_P_disc<0,3,3>;
103  break;
104 
105  case 1:
106  q_order = 2;
107  fe1_ = new FE_P_disc<1,1,3>;
108  fe2_ = new FE_P_disc<1,2,3>;
109  fe3_ = new FE_P_disc<1,3,3>;
110  break;
111 
112  case 2:
113  q_order = 4;
114  fe1_ = new FE_P_disc<2,1,3>;
115  fe2_ = new FE_P_disc<2,2,3>;
116  fe3_ = new FE_P_disc<2,3,3>;
117  break;
118 
119  case 3:
120  q_order = 6;
121  fe1_ = new FE_P_disc<3,1,3>;
122  fe2_ = new FE_P_disc<3,2,3>;
123  fe3_ = new FE_P_disc<3,3,3>;
124  break;
125 
126  default:
127  xprintf(PrgErr, "Unsupported polynomial order %d for finite elements in TransportDG ", fe_order);
128  break;
129  }
130 
131  fe_rt1_ = new FE_RT0<1,3>;
132  fe_rt2_ = new FE_RT0<2,3>;
133  fe_rt3_ = new FE_RT0<3,3>;
134 
135  q0_ = new QGauss<0>(q_order);
136  q1_ = new QGauss<1>(q_order);
137  q2_ = new QGauss<2>(q_order);
138  q3_ = new QGauss<3>(q_order);
139 
140  map0_ = new MappingP1<0,3>;
141  map1_ = new MappingP1<1,3>;
142  map2_ = new MappingP1<2,3>;
143  map3_ = new MappingP1<3,3>;
144 
145  dh_ = new DOFHandlerMultiDim(*mesh_);
146 
147  dh_->distribute_dofs(*fe1_, *fe2_, *fe3_);
148 }
149 
150 
151 FEObjects::~FEObjects()
152 {
153  delete fe1_;
154  delete fe2_;
155  delete fe3_;
156  delete fe_rt1_;
157  delete fe_rt2_;
158  delete fe_rt3_;
159  delete q0_;
160  delete q1_;
161  delete q2_;
162  delete q3_;
163  delete map0_;
164  delete map1_;
165  delete map2_;
166  delete map3_;
167  delete dh_;
168 }
169 
170 template<> FiniteElement<0,3> *FEObjects::fe<0>() { return 0; }
171 template<> FiniteElement<1,3> *FEObjects::fe<1>() { return fe1_; }
172 template<> FiniteElement<2,3> *FEObjects::fe<2>() { return fe2_; }
173 template<> FiniteElement<3,3> *FEObjects::fe<3>() { return fe3_; }
174 
175 template<> FiniteElement<0,3> *FEObjects::fe_rt<0>() { return 0; }
176 template<> FiniteElement<1,3> *FEObjects::fe_rt<1>() { return fe_rt1_; }
177 template<> FiniteElement<2,3> *FEObjects::fe_rt<2>() { return fe_rt2_; }
178 template<> FiniteElement<3,3> *FEObjects::fe_rt<3>() { return fe_rt3_; }
179 
180 template<> Quadrature<0> *FEObjects::q<0>() { return q0_; }
181 template<> Quadrature<1> *FEObjects::q<1>() { return q1_; }
182 template<> Quadrature<2> *FEObjects::q<2>() { return q2_; }
183 template<> Quadrature<3> *FEObjects::q<3>() { return q3_; }
184 
185 template<> Mapping<0,3> *FEObjects::mapping<0>() { return map0_; }
186 template<> Mapping<1,3> *FEObjects::mapping<1>() { return map1_; }
187 template<> Mapping<2,3> *FEObjects::mapping<2>() { return map2_; }
188 template<> Mapping<3,3> *FEObjects::mapping<3>() { return map3_; }
189 
190 DOFHandlerMultiDim *FEObjects::dh() { return dh_; }
191 
192 
193 template<class Model>
194 TransportDG<Model>::EqData::EqData() : Model::ModelEqData()
195 {
196  ADD_FIELD(fracture_sigma, "Coefficient of diffusive transfer through fractures (for each substance).", "1.0");
197  ADD_FIELD(dg_penalty, "Penalty parameter influencing the discontinuity of the solution (for each substance). "
198  "Its default value 1 is sufficient in most cases. Higher value diminishes the inter-element jumps.", "1.0");
199 
200  ADD_FIELD(bc_type,"Boundary condition type, possible values: inflow, dirichlet, neumann, robin.", "\"inflow\"" );
201  bc_type.input_selection(&bc_type_selection);
202 
203 // std::vector<FieldEnum> list; list.push_back(neumann);
204  ADD_FIELD(bc_flux,"Flux in Neumann boundary condition.", "0.0");
205 // bc_flux.disable_where(bc_type, { dirichlet, inflow });
206  ADD_FIELD(bc_robin_sigma,"Conductivity coefficient in Robin boundary condition.", "0.0");
207 // bc_robin_sigma.disable_where(bc_type, {dirichlet, inflow, neumann});
208 
209  // add all input fields to the output list
210  this->output_fields += *this;
211 }
212 
213 
214 
215 template<class Model>
216 TransportDG<Model>::TransportDG(Mesh & init_mesh, const Input::Record &in_rec)
217  : TransportBase(init_mesh, in_rec),
218  mass_matrix(0),
219  allocation_done(false)
220 {
221  // Check that Model is derived from AdvectionDiffusionModel.
222  static_assert(std::is_base_of<AdvectionDiffusionModel, Model>::value, "");
223 
224  this->eq_data_ = &data_;
225 
226  time_ = new TimeGovernor(in_rec.val<Input::Record>("time"));
227  time_->fix_dt_until_mark();
228 
229  // Read names of transported substances.
230  Model::set_component_names(subst_names_, in_rec);
231  n_subst_ = subst_names_.size();
232 
233  Input::Iterator<Input::Record> it = in_rec.find<Input::Record>("mass_balance");
234  if (it)
235  mass_balance_ = new MassBalance(this, *it);
236 
237  // Set up physical parameters.
238  data_.set_mesh(init_mesh);
239  data_.bc_type.n_comp(n_subst_);
240  data_.bc_flux.n_comp(n_subst_);
241  data_.bc_robin_sigma.n_comp(n_subst_);
242  data_.fracture_sigma.n_comp(n_subst_);
243  data_.dg_penalty.n_comp(n_subst_);
244  Model::init_data(n_subst_);
245  data_.set_input_list( in_rec.val<Input::Array>("input_fields") );
246  data_.set_limit_side(LimitSide::left);
247 
248 
249  // DG variant
250  dg_variant = in_rec.val<DGVariant>("dg_variant");
251 
252  // DG stabilization parameters on boundary edges
253  gamma.resize(n_subst_);
254  for (int sbi=0; sbi<n_subst_; sbi++)
255  gamma[sbi].resize(mesh_->boundary_.size());
256 
257 
258  // create finite element structures and distribute DOFs
259  dg_order = in_rec.val<unsigned int>("dg_order");
260  feo = new FEObjects(mesh_, dg_order);
261  DBGMSG("TDG: solution size %d\n", feo->dh()->n_global_dofs());
262 
263  // Resize coefficient arrays
264  int qsize = max(feo->q<0>()->size(), max(feo->q<1>()->size(), max(feo->q<2>()->size(), feo->q<3>()->size())));
265  int max_edg_sides = max(mesh_->max_edge_sides(1), max(mesh_->max_edge_sides(2), mesh_->max_edge_sides(3)));
266  mm_coef.resize(qsize);
267  ad_coef.resize(n_subst_);
268  dif_coef.resize(n_subst_);
269  for (int sbi=0; sbi<n_subst_; sbi++)
270  {
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(n_subst_);
279  dif_coef_edg[sd].resize(n_subst_);
280  for (int sbi=0; sbi<n_subst_; sbi++)
281  {
282  ad_coef_edg[sd][sbi].resize(qsize);
283  dif_coef_edg[sd][sbi].resize(qsize);
284  }
285  }
286 
287  // register output fields
288  output_rec = in_rec.val<Input::Record>("output_stream");
289  output_vec.resize(n_subst_);
290  output_solution.resize(n_subst_);
291  for (int sbi=0; sbi<n_subst_; sbi++)
292  {
293  // for each substance we allocate output array and vector
294  output_solution[sbi] = new double[feo->dh()->n_global_dofs()];
295  VecCreateSeqWithArray(PETSC_COMM_SELF, 1, feo->dh()->n_global_dofs(), output_solution[sbi], &output_vec[sbi]);
296  }
297  data_.output_field.init(subst_names_);
298  data_.output_field.set_mesh(*mesh_);
299  data_.output_fields.output_type(OutputTime::CORNER_DATA);
300 
301  for (int sbi=0; sbi<n_subst_; sbi++)
302  {
303  // create shared pointer to a FieldFE, pass FE data and push this FieldFE to output_field on all regions
304  std::shared_ptr<FieldFE<3, FieldValue<3>::Scalar> > output_field_ptr(new FieldFE<3, FieldValue<3>::Scalar>);
305  output_field_ptr->set_fe_data(feo->dh(), feo->mapping<1>(), feo->mapping<2>(), feo->mapping<3>(), &output_vec[sbi]);
306  data_.output_field[sbi].set_field(mesh_->region_db().get_region_set("ALL"), output_field_ptr, 0);
307  }
308  data_.output_fields.set_limit_side(LimitSide::left);
309  output_stream = OutputTime::create_output_stream(output_rec);
310  output_stream->add_admissible_field_names(in_rec.val<Input::Array>("output_fields"), data_.output_selection);
311 
312  // set time marks for writing the output
313  output_stream->mark_output_times(*time_);
314  //output_mark_type = this->mark_type() | time_->marks().type_fixed_time() | time_->marks().type_output();
315  //time_->marks().add_time_marks(0.0, output_rec.val<double>("time_step"), time_->end_time(), output_mark_type);
316 
317  // allocate matrix and vector structures
318 
319  ls = new LinSys*[n_subst_];
320  ls_dt = new LinSys_PETSC(feo->dh()->distr());
321  ( (LinSys_PETSC *)ls_dt )->set_from_input( in_rec.val<Input::Record>("solver") );
322  for (int sbi = 0; sbi < n_subst_; sbi++) {
323  ls[sbi] = new LinSys_PETSC(feo->dh()->distr());
324  ( (LinSys_PETSC *)ls[sbi] )->set_from_input( in_rec.val<Input::Record>("solver") );
325  ls[sbi]->set_solution(NULL);
326  }
327  stiffness_matrix = new Mat[n_subst_];
328  rhs = new Vec[n_subst_];
329 
330 }
331 
332 template<class Model>
333 TransportDG<Model>::~TransportDG()
334 {
335  delete time_;
336  delete ls_dt;
337 
338  if (feo->dh()->el_ds()->myp() == 0)
339  {
340  for (int i=0; i<n_subst_; i++)
341  {
342  VecDestroy(&output_vec[i]);
343  delete[] output_solution[i];
344  }
345  }
346 
347  for (int i=0; i<n_subst_; i++)
348  {
349  delete ls[i];
350  MatDestroy(&stiffness_matrix[i]);
351  VecDestroy(&rhs[i]);
352  }
353  delete[] ls;
354  delete[] stiffness_matrix;
355  delete[] rhs;
356  delete feo;
357  if (mass_balance_ != NULL) delete mass_balance_;
358 
359  gamma.clear();
360  subst_names_.clear();
361  delete output_stream;
362 }
363 
364 
365 template<class Model>
366 void TransportDG<Model>::output_vector_gather()
367 {
368  IS is;
369  VecScatter output_scatter;
370  int idx[] = { 0 };
371  for (int sbi=0; sbi<n_subst_; sbi++)
372  {
373  // gather solution to output_vec[sbi]
374  ISCreateBlock(PETSC_COMM_SELF, ls[sbi]->size(), 1, idx, PETSC_COPY_VALUES, &is);
375  VecScatterCreate(ls[sbi]->get_solution(), is, output_vec[sbi], PETSC_NULL, &output_scatter);
376  VecScatterBegin(output_scatter, ls[sbi]->get_solution(), output_vec[sbi], INSERT_VALUES, SCATTER_FORWARD);
377  VecScatterEnd(output_scatter, ls[sbi]->get_solution(), output_vec[sbi], INSERT_VALUES, SCATTER_FORWARD);
378  VecScatterDestroy(&(output_scatter));
379  ISDestroy(&(is));
380  }
381 }
382 
383 
384 
385 template<class Model>
386 void TransportDG<Model>::zero_time_step()
387 {
388  data_.set_time(*time_);
389 
390  // set initial conditions
391  set_initial_condition();
392  for (int sbi = 0; sbi < n_subst_; sbi++)
393  ( (LinSys_PETSC *)ls[sbi] )->set_initial_guess_nonzero();
394 
395  output_data();
396 }
397 
398 
399 template<class Model>
400 void TransportDG<Model>::update_solution()
401 {
402  START_TIMER("DG-ONE STEP");
403 
404  time_->next_time();
405  time_->view("TDG");
406 
407  START_TIMER("data reinit");
408  data_.set_time(*time_);
409  END_TIMER("data reinit");
410 
411  // check first time assembly - needs preallocation
412  if (!allocation_done)
413  {
414  // preallocate mass matrix
415  ls_dt->start_allocation();
416  assemble_mass_matrix();
417  mass_matrix = NULL;
418 
419  // preallocate system matrix
420  for (int i=0; i<n_subst_; i++)
421  {
422  ls[i]->start_allocation();
423  stiffness_matrix[i] = NULL;
424  rhs[i] = NULL;
425  }
426  assemble_stiffness_matrix();
427  set_sources();
428  set_boundary_conditions();
429 
430  allocation_done = true;
431  }
432 
433  // assemble mass matrix
434  if (mass_matrix == NULL ||
435  mass_matrix_changed())
436  {
437  ls_dt->start_add_assembly();
438  ls_dt->mat_zero_entries();
439  assemble_mass_matrix();
440  ls_dt->finish_assembly();
441  mass_matrix = ls_dt->get_matrix();
442  }
443 
444  // assemble stiffness matrix
445  if (stiffness_matrix[0] == NULL ||
446  stiffness_matrix_changed())
447  {
448  // new fluxes can change the location of Neumann boundary,
449  // thus stiffness matrix must be reassembled
450  for (int i=0; i<n_subst_; i++)
451  {
452  ls[i]->start_add_assembly();
453  ls[i]->mat_zero_entries();
454  }
455  assemble_stiffness_matrix();
456  for (int i=0; i<n_subst_; i++)
457  {
458  ls[i]->finish_assembly();
459 
460  if (stiffness_matrix[i] == NULL)
461  MatConvert(ls[i]->get_matrix(), MATSAME, MAT_INITIAL_MATRIX, &stiffness_matrix[i]);
462  else
463  MatCopy(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  rhs_changed())
470  {
471  for (int i=0; i<n_subst_; i++)
472  {
473  ls[i]->start_add_assembly();
474  ls[i]->rhs_zero_entries();
475  }
476  set_sources();
477  set_boundary_conditions();
478  for (int i=0; i<n_subst_; i++)
479  {
480  ls[i]->finish_assembly();
481 
482  VecDuplicate(ls[i]->get_rhs(), &rhs[i]);
483  VecCopy(ls[i]->get_rhs(), rhs[i]);
484  }
485  }
486 
487  Model::flux_changed = false;
488 
489 
490  /* Apply backward Euler time integration.
491  *
492  * Denoting A the stiffness matrix and M the mass matrix, the algebraic system at the k-th time level reads
493  *
494  * (1/dt M + A)u^k = f + 1/dt M.u^{k-1}
495  *
496  * Hence we modify at each time level the right hand side:
497  *
498  * f^k = f + 1/dt M u^{k-1},
499  *
500  * where f stands for the term stemming from the force and boundary conditions.
501  * Accordingly, we set
502  *
503  * A^k = A + 1/dt M.
504  *
505  */
506  Mat m;
507  START_TIMER("solve");
508  for (int i=0; i<n_subst_; i++)
509  {
510  MatConvert(stiffness_matrix[i], MATSAME, MAT_INITIAL_MATRIX, &m);
511  MatAXPY(m, 1./time_->dt(), mass_matrix, SUBSET_NONZERO_PATTERN);
512  ls[i]->set_matrix(m, DIFFERENT_NONZERO_PATTERN);
513  Vec y,w;
514  VecDuplicate(rhs[i], &y);
515  VecDuplicate(rhs[i], &w);
516  MatMult(mass_matrix, ls[i]->get_solution(), y);
517  VecWAXPY(w, 1./time_->dt(), y, rhs[i]);
518  ls[i]->set_rhs(w);
519 
520  VecDestroy(&y);
521  VecDestroy(&w);
522  MatDestroy(&m);
523 
524  ls[i]->solve();
525  }
526  END_TIMER("solve");
527 
528  if (mass_balance() != NULL)
529  mass_balance()->calculate(time_->t());
530 
531  END_TIMER("DG-ONE STEP");
532 }
533 
534 
535 
536 template<class Model>
537 void TransportDG<Model>::get_solution_vector(double *& vector, unsigned int & size)
538 {}
539 
540 
541 template<class Model>
542 void TransportDG<Model>::get_parallel_solution_vector(Vec & vector)
543 {}
544 
545 
546 template<class Model>
547 void TransportDG<Model>::set_velocity_field(const MH_DofHandler &dh)
548 {
549  // So far the velocity_vector contains zeros, so we ignore it.
550  // Instead we use the value Side.flux.
551 
552  mh_dh = &dh;
553  Model::flux_changed = true;
554 
555 }
556 
557 
558 template<class Model>
559 void TransportDG<Model>::output_data()
560 {
561  double *solution;
562  unsigned int dof_indices[max(feo->fe<1>()->n_dofs(), max(feo->fe<2>()->n_dofs(), feo->fe<3>()->n_dofs()))];
563 
564  if (!time_->is_current( time_->marks().type_output() )) return;
565 
566  START_TIMER("DG-OUTPUT");
567 
568  // gather the solution from all processors
569  output_vector_gather();
570  data_.output_fields.set_time(*time_);
571  data_.output_fields.output(output_stream);
572  output_stream->write_time_frame();
573 
574  if (mass_balance() != NULL)
575  mass_balance()->output(time_->t());
576 
577  END_TIMER("DG-OUTPUT");
578 }
579 
580 
581 template<class Model>
582 void TransportDG<Model>::assemble_mass_matrix()
583 {
584  START_TIMER("assemble_mass");
585  assemble_mass_matrix<1>();
586  assemble_mass_matrix<2>();
587  assemble_mass_matrix<3>();
588  END_TIMER("assemble_mass");
589 }
590 
591 
592 template<class Model> template<unsigned int dim>
593 void TransportDG<Model>::assemble_mass_matrix()
594 {
595  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(), update_values | update_JxW_values | update_quadrature_points);
596  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
597  unsigned int dof_indices[ndofs];
598  PetscScalar local_mass_matrix[ndofs*ndofs];
599 
600  // assemble integral over elements
601  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
602  {
603  typename DOFHandlerBase::CellIterator cell = mesh_->element(feo->dh()->el_index(i_cell));
604  if (cell->dim() != dim) continue;
605 
606  fe_values.reinit(cell);
607  feo->dh()->get_dof_indices(cell, dof_indices);
608  ElementAccessor<3> ele_acc = cell->element_accessor();
609 
610  Model::compute_mass_matrix_coefficient(fe_values.point_list(), ele_acc, mm_coef);
611 
612  // assemble the local stiffness and mass matrix
613  for (unsigned int i=0; i<ndofs; i++)
614  {
615  for (unsigned int j=0; j<ndofs; j++)
616  {
617  local_mass_matrix[i*ndofs+j] = 0;
618  for (unsigned int k=0; k<qsize; k++)
619  local_mass_matrix[i*ndofs+j] += mm_coef[k]*fe_values.shape_value(j,k)*fe_values.shape_value(i,k)*fe_values.JxW(k);
620  }
621  }
622 
623  ls_dt->mat_set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, local_mass_matrix);
624  }
625 }
626 
627 
628 
629 
630 template<class Model>
631 void TransportDG<Model>::assemble_stiffness_matrix()
632 {
633  START_TIMER("assemble_stiffness");
634  START_TIMER("assemble_volume_integrals");
635  assemble_volume_integrals<1>();
636  assemble_volume_integrals<2>();
637  assemble_volume_integrals<3>();
638  END_TIMER("assemble_volume_integrals");
639 
640  START_TIMER("assemble_fluxes_boundary");
641  assemble_fluxes_boundary<1>();
642  assemble_fluxes_boundary<2>();
643  assemble_fluxes_boundary<3>();
644  END_TIMER("assemble_fluxes_boundary");
645 
646  START_TIMER("assemble_fluxes_elem_elem");
647  assemble_fluxes_element_element<1>();
648  assemble_fluxes_element_element<2>();
649  assemble_fluxes_element_element<3>();
650  END_TIMER("assemble_fluxes_elem_elem");
651 
652  START_TIMER("assemble_fluxes_elem_side");
653  assemble_fluxes_element_side<1>();
654  assemble_fluxes_element_side<2>();
655  assemble_fluxes_element_side<3>();
656  END_TIMER("assemble_fluxes_elem_side");
657  END_TIMER("assemble_stiffness");
658 }
659 
660 
661 
662 template<class Model>
663 template<unsigned int dim>
664 void TransportDG<Model>::assemble_volume_integrals()
665 {
666  FEValues<dim,3> fv_rt(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe_rt<dim>(),
667  update_values | update_gradients);
668  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
669  update_values | update_gradients | update_JxW_values | update_quadrature_points);
670  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
671  unsigned int dof_indices[ndofs];
672  vector<arma::vec3> velocity(qsize);
673  vector<arma::vec> sources_sigma(qsize);
674  PetscScalar local_matrix[ndofs*ndofs];
675 
676  // assemble integral over elements
677  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
678  {
679  typename DOFHandlerBase::CellIterator cell = mesh_->element(feo->dh()->el_index(i_cell));
680  if (cell->dim() != dim) continue;
681 
682  fe_values.reinit(cell);
683  fv_rt.reinit(cell);
684  ElementAccessor<3> ele_acc = cell->element_accessor();
685  feo->dh()->get_dof_indices(cell, dof_indices);
686 
687  calculate_velocity(cell, velocity, fv_rt);
688  Model::compute_advection_diffusion_coefficients(fe_values.point_list(), velocity, ele_acc, ad_coef, dif_coef);
689  Model::compute_sources_sigma(fe_values.point_list(), ele_acc, sources_sigma);
690 
691  // assemble the local stiffness matrix
692  for (int sbi=0; sbi<n_subst_; sbi++)
693  {
694  for (unsigned int i=0; i<ndofs; i++)
695  for (unsigned int j=0; j<ndofs; j++)
696  local_matrix[i*ndofs+j] = 0;
697 
698  for (unsigned int k=0; k<qsize; k++)
699  {
700  for (unsigned int i=0; i<ndofs; i++)
701  {
702  arma::vec3 Kt_grad_i = dif_coef[sbi][k].t()*fe_values.shape_grad(i,k);
703  double ad_dot_grad_i = arma::dot(ad_coef[sbi][k], fe_values.shape_grad(i,k));
704 
705  for (unsigned int j=0; j<ndofs; j++)
706  local_matrix[i*ndofs+j] += (arma::dot(Kt_grad_i, fe_values.shape_grad(j,k))
707  -fe_values.shape_value(j,k)*ad_dot_grad_i
708  +sources_sigma[k][sbi]*fe_values.shape_value(j,k)*fe_values.shape_value(i,k))*fe_values.JxW(k);
709  }
710  }
711  ls[sbi]->mat_set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, local_matrix);
712  }
713  }
714 }
715 
716 
717 template<class Model>
718 void TransportDG<Model>::set_sources()
719 {
720  START_TIMER("assemble_sources");
721  set_sources<1>();
722  set_sources<2>();
723  set_sources<3>();
724  END_TIMER("assemble_sources");
725 }
726 
727 template<class Model>
728 template<unsigned int dim>
729 void TransportDG<Model>::set_sources()
730 {
731  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
732  update_values | update_JxW_values | update_quadrature_points);
733  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
734  vector<arma::vec> sources_conc(qsize), sources_density(qsize), sources_sigma(qsize);
735  unsigned int dof_indices[ndofs];
736  PetscScalar local_rhs[ndofs];
737  double source;
738 
739  // assemble integral over elements
740  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
741  {
742  typename DOFHandlerBase::CellIterator cell = mesh_->element(feo->dh()->el_index(i_cell));
743  if (cell->dim() != dim) continue;
744 
745  fe_values.reinit(cell);
746  feo->dh()->get_dof_indices(cell, dof_indices);
747 
748  Model::compute_source_coefficients(fe_values.point_list(), cell->element_accessor(), sources_conc, sources_density, sources_sigma);
749 
750  // assemble the local stiffness matrix
751  for (int sbi=0; sbi<n_subst_; sbi++)
752  {
753  for (unsigned int i=0; i<ndofs; i++)
754  local_rhs[i] = 0;
755 
756  // compute sources
757  for (unsigned int k=0; k<qsize; k++)
758  {
759  source = (sources_density[k][sbi] + sources_conc[k][sbi]*sources_sigma[k][sbi])*fe_values.JxW(k);
760 
761  for (unsigned int i=0; i<ndofs; i++)
762  local_rhs[i] += source*fe_values.shape_value(i,k);
763  }
764  ls[sbi]->rhs_set_values(ndofs, (int *)dof_indices, local_rhs);
765  }
766  }
767 }
768 
769 
770 
771 template<class Model>
772 template<unsigned int dim>
773 void TransportDG<Model>::assemble_fluxes_element_element()
774 {
775  vector<FESideValues<dim,3>*> fe_values;
776  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
777  update_values);
778  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size(),
779  n_max_sides = ad_coef_edg.size();
780  vector<unsigned int*> side_dof_indices;
781  PetscScalar local_matrix[ndofs*ndofs];
782  vector<vector<arma::vec3> > side_velocity(n_max_sides);
783  vector<arma::vec> dg_penalty(n_max_sides);
784  double gamma_l, omega[2], transport_flux;
785 
786  for (int sid=0; sid<n_max_sides; sid++)
787  {
788  side_dof_indices.push_back(new unsigned int[ndofs]);
789  fe_values.push_back(new FESideValues<dim,3>(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
790  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points));
791  }
792 
793  // assemble integral over sides
794  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
795  {
796  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
797  if (edg->n_sides < 2 || edg->side(0)->element()->dim() != dim) continue;
798 
799  for (int sid=0; sid<edg->n_sides; sid++)
800  {
801  typename DOFHandlerBase::CellIterator cell = edg->side(sid)->element();
802  ElementAccessor<3> ele_acc = cell->element_accessor();
803  feo->dh()->get_dof_indices(cell, side_dof_indices[sid]);
804  fe_values[sid]->reinit(cell, edg->side(sid)->el_idx());
805  fsv_rt.reinit(cell, edg->side(sid)->el_idx());
806  calculate_velocity(cell, side_velocity[sid], fsv_rt);
807  Model::compute_advection_diffusion_coefficients(fe_values[sid]->point_list(), side_velocity[sid], ele_acc, ad_coef_edg[sid], dif_coef_edg[sid]);
808  dg_penalty[sid] = data_.dg_penalty.value(cell->centre(), ele_acc);
809  }
810 
811  // fluxes and penalty
812  for (int sbi=0; sbi<n_subst_; sbi++)
813  {
814  vector<double> fluxes(edg->n_sides);
815  for (unsigned int sid=0; sid<edg->n_sides; sid++)
816  {
817  fluxes[sid] = 0;
818  for (unsigned int k=0; k<qsize; k++)
819  fluxes[sid] += arma::dot(ad_coef_edg[sid][sbi][k], fe_values[sid]->normal_vector(k))*fe_values[sid]->JxW(k);
820  fluxes[sid] /= edg->side(sid)->measure();
821  }
822 
823  for (int s1=0; s1<edg->n_sides; s1++)
824  {
825  for (int s2=s1+1; s2<edg->n_sides; s2++)
826  {
827  ASSERT(edg->side(s1)->valid(), "Invalid side of edge.");
828  if (!feo->dh()->el_is_local(edg->side(s1)->element().index())
829  && !feo->dh()->el_is_local(edg->side(s2)->element().index())) continue;
830 
831  arma::vec3 nv = fe_values[s1]->normal_vector(0);
832 
833  // set up the parameters for DG method
834  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);
835 
836  int sd[2];
837  sd[0] = s1;
838  sd[1] = s2;
839 
840 #define AVERAGE(i,k,side_id) (fe_values[sd[side_id]]->shape_value(i,k)*0.5)
841 #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])
842 #define JUMP(i,k,side_id) ((side_id==0?1:-1)*fe_values[sd[side_id]]->shape_value(i,k))
843 
844  // For selected pair of elements:
845  for (int n=0; n<2; n++)
846  {
847  if (!feo->dh()->el_is_local(edg->side(sd[n])->element().index())) continue;
848 
849  for (int m=0; m<2; m++)
850  {
851  for (unsigned int i=0; i<fe_values[sd[n]]->n_dofs(); i++)
852  for (unsigned int j=0; j<fe_values[sd[m]]->n_dofs(); j++)
853  local_matrix[i*fe_values[sd[m]]->n_dofs()+j] = 0;
854 
855  for (unsigned int k=0; k<qsize; k++)
856  {
857  double flux_times_JxW = transport_flux*fe_values[0]->JxW(k);
858  double gamma_times_JxW = gamma_l*fe_values[0]->JxW(k);
859 
860  for (unsigned int i=0; i<fe_values[sd[n]]->n_dofs(); i++)
861  {
862  double flux_JxW_jump_i = flux_times_JxW*JUMP(i,k,n);
863  double gamma_JxW_jump_i = gamma_times_JxW*JUMP(i,k,n);
864  double JxW_jump_i = fe_values[0]->JxW(k)*JUMP(i,k,n);
865  double JxW_var_wavg_i = fe_values[0]->JxW(k)*WAVERAGE(i,k,n)*dg_variant;
866 
867  for (unsigned int j=0; j<fe_values[sd[m]]->n_dofs(); j++)
868  {
869  int index = i*fe_values[sd[m]]->n_dofs()+j;
870 
871  // flux due to transport (applied on interior edges) (average times jump)
872  local_matrix[index] += flux_JxW_jump_i*AVERAGE(j,k,m);
873 
874  // penalty enforcing continuity across edges (applied on interior and Dirichlet edges) (jump times jump)
875  local_matrix[index] += gamma_JxW_jump_i*JUMP(j,k,m);
876 
877  // terms due to diffusion
878  local_matrix[index] -= WAVERAGE(j,k,m)*JxW_jump_i;
879  local_matrix[index] -= JUMP(j,k,m)*JxW_var_wavg_i;
880  }
881  }
882  }
883  ls[sbi]->mat_set_values(fe_values[sd[n]]->n_dofs(), (int *)side_dof_indices[sd[n]], fe_values[sd[m]]->n_dofs(), (int *)side_dof_indices[sd[m]], local_matrix);
884  }
885  }
886 #undef AVERAGE
887 #undef WAVERAGE
888 #undef JUMP
889  }
890  }
891  }
892  }
893 
894  for (unsigned int i=0; i<n_max_sides; i++)
895  {
896  delete fe_values[i];
897  delete[] side_dof_indices[i];
898  }
899 }
900 
901 
902 template<class Model>
903 template<unsigned int dim>
904 void TransportDG<Model>::assemble_fluxes_boundary()
905 {
906  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
907  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points);
908  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
909  update_values);
910  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
911  unsigned int side_dof_indices[ndofs];
912  PetscScalar local_matrix[ndofs*ndofs];
913  vector<arma::vec3> side_velocity;
914  vector<arma::vec> robin_sigma(qsize);
915  arma::vec dg_penalty;
916  double gamma_l;
917 
918  // assemble boundary integral
919  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
920  {
921  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
922  if (edg->n_sides > 1) continue;
923  // check spatial dimension
924  if (edg->side(0)->dim() != dim-1) continue;
925  // skip edges lying not on the boundary
926  if (edg->side(0)->cond() == NULL) continue;
927 
928  SideIter side = edg->side(0);
929  typename DOFHandlerBase::CellIterator cell = side->element();
930  ElementAccessor<3> ele_acc = cell->element_accessor();
931  feo->dh()->get_dof_indices(cell, side_dof_indices);
932  fe_values_side.reinit(cell, side->el_idx());
933  fsv_rt.reinit(cell, side->el_idx());
934 
935  calculate_velocity(cell, side_velocity, fsv_rt);
936  Model::compute_advection_diffusion_coefficients(fe_values_side.point_list(), side_velocity, ele_acc, ad_coef, dif_coef);
937  dg_penalty = data_.dg_penalty.value(cell->centre(), ele_acc);
938  arma::uvec bc_type = data_.bc_type.value(side->cond()->element()->centre(), side->cond()->element_accessor());
939  data_.bc_robin_sigma.value_list(fe_values_side.point_list(), side->cond()->element_accessor(), robin_sigma);
940 
941  for (int sbi=0; sbi<n_subst_; sbi++)
942  {
943  for (unsigned int i=0; i<ndofs; i++)
944  for (unsigned int j=0; j<ndofs; j++)
945  local_matrix[i*ndofs+j] = 0;
946 
947  // On Neumann boundaries we have only term from integrating by parts the advective term,
948  // on Dirichlet boundaries we additionally apply the penalty which enforces the prescribed value.
949  double side_flux = 0;
950  for (unsigned int k=0; k<qsize; k++)
951  side_flux += arma::dot(ad_coef[sbi][k], fe_values_side.normal_vector(k))*fe_values_side.JxW(k);
952  double transport_flux = side_flux/side->measure();
953 
954  if ((bc_type[sbi] == EqData::dirichlet) || (bc_type[sbi] == EqData::inflow ))
955  {
956  // set up the parameters for DG method
957  set_DG_parameters_boundary(side, qsize, dif_coef[sbi], transport_flux, fe_values_side.normal_vector(0), dg_penalty[sbi], gamma_l);
958  gamma[sbi][side->cond_idx()] = gamma_l;
959  if (bc_type[sbi] == EqData::dirichlet || side_flux < -mh_dh->precision())
960  transport_flux += gamma_l;
961  }
962 
963  // fluxes and penalty
964  for (unsigned int k=0; k<qsize; k++)
965  {
966  double flux_times_JxW;
967  if (bc_type[sbi] == EqData::robin)
968  flux_times_JxW = (transport_flux + robin_sigma[k][sbi])*fe_values_side.JxW(k);
969  else
970  flux_times_JxW = transport_flux*fe_values_side.JxW(k);
971 
972  for (unsigned int i=0; i<ndofs; i++)
973  for (unsigned int j=0; j<ndofs; j++)
974  {
975  // flux due to advection and penalty
976  local_matrix[i*ndofs+j] += flux_times_JxW*fe_values_side.shape_value(i,k)*fe_values_side.shape_value(j,k);
977 
978  // flux due to diffusion (only on dirichlet and inflow boundary)
979  if (bc_type[sbi] == EqData::dirichlet || (bc_type[sbi] == EqData::inflow && side_flux < -mh_dh->precision()))
980  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)
981  + 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
982  )*fe_values_side.JxW(k);
983  }
984  }
985  ls[sbi]->mat_set_values(ndofs, (int *)side_dof_indices, ndofs, (int *)side_dof_indices, local_matrix);
986  }
987  }
988 }
989 
990 
991 template<class Model>
992 template<unsigned int dim>
993 void TransportDG<Model>::assemble_fluxes_element_side()
994 {
995 
996  if (dim == 1) return;
997  FEValues<dim-1,3> fe_values_vb(*feo->mapping<dim-1>(), *feo->q<dim-1>(), *feo->fe<dim-1>(),
998  update_values | update_gradients | update_JxW_values | update_quadrature_points);
999  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1000  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points);
1001  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1002  update_values);
1003  FEValues<dim-1,3> fv_rt(*feo->mapping<dim-1>(), *feo->q<dim-1>(), *feo->fe_rt<dim-1>(),
1004  update_values);
1005 
1006  vector<FEValuesSpaceBase<3>*> fv_sb(2);
1007  const unsigned int ndofs = feo->fe<dim>()->n_dofs(); // number of local dofs
1008  const unsigned int qsize = feo->q<dim-1>()->size(); // number of quadrature points
1009  unsigned int side_dof_indices[2*ndofs], n_dofs[2];
1010  vector<arma::vec3> velocity_higher, velocity_lower;
1011  vector<arma::vec> frac_sigma(qsize);
1012  vector<double> csection_lower(qsize), csection_higher(qsize), mm_coef_lower(qsize), mm_coef_higher(qsize);
1013  PetscScalar local_matrix[4*ndofs*ndofs];
1014  double comm_flux[2][2];
1015 
1016  // index 0 = element with lower dimension,
1017  // index 1 = side of element with higher dimension
1018  fv_sb[0] = &fe_values_vb;
1019  fv_sb[1] = &fe_values_side;
1020 
1021  // assemble integral over sides
1022  for (int inb=0; inb<feo->dh()->n_loc_nb(); inb++)
1023  {
1024  Neighbour *nb = &mesh_->vb_neighbours_[feo->dh()->nb_index(inb)];
1025  // skip neighbours of different dimension
1026  if (nb->element()->dim() != dim-1) continue;
1027 
1028  typename DOFHandlerBase::CellIterator cell_sub = mesh_->element.full_iter(nb->element());
1029  feo->dh()->get_dof_indices(cell_sub, side_dof_indices);
1030  fe_values_vb.reinit(cell_sub);
1031  n_dofs[0] = fv_sb[0]->n_dofs();
1032 
1033  typename DOFHandlerBase::CellIterator cell = nb->side()->element();
1034  feo->dh()->get_dof_indices(cell, side_dof_indices+n_dofs[0]);
1035  fe_values_side.reinit(cell, nb->side()->el_idx());
1036  n_dofs[1] = fv_sb[1]->n_dofs();
1037 
1038  // Element id's for testing if they belong to local partition.
1039  int element_id[2];
1040  element_id[0] = cell_sub.index();
1041  element_id[1] = cell.index();
1042 
1043  fsv_rt.reinit(cell, nb->side()->el_idx());
1044  fv_rt.reinit(cell_sub);
1045  calculate_velocity(cell, velocity_higher, fsv_rt);
1046  calculate_velocity(cell_sub, velocity_lower, fv_rt);
1047  Model::compute_advection_diffusion_coefficients(fe_values_vb.point_list(), velocity_lower, cell_sub->element_accessor(), ad_coef_edg[0], dif_coef_edg[0]);
1048  Model::compute_advection_diffusion_coefficients(fe_values_vb.point_list(), velocity_higher, cell->element_accessor(), ad_coef_edg[1], dif_coef_edg[1]);
1049  Model::compute_mass_matrix_coefficient(fe_values_vb.point_list(), cell_sub->element_accessor(), mm_coef_lower);
1050  Model::compute_mass_matrix_coefficient(fe_values_vb.point_list(), cell->element_accessor(), mm_coef_higher);
1051  data_.cross_section.value_list(fe_values_vb.point_list(), cell_sub->element_accessor(), csection_lower);
1052  data_.cross_section.value_list(fe_values_vb.point_list(), cell->element_accessor(), csection_higher);
1053  data_.fracture_sigma.value_list(fe_values_vb.point_list(), cell_sub->element_accessor(), frac_sigma);
1054 
1055  for (int sbi=0; sbi<n_subst_; sbi++)
1056  {
1057  for (unsigned int i=0; i<n_dofs[0]+n_dofs[1]; i++)
1058  for (unsigned int j=0; j<n_dofs[0]+n_dofs[1]; j++)
1059  local_matrix[i*(n_dofs[0]+n_dofs[1])+j] = 0;
1060 
1061  // set transmission conditions
1062  for (unsigned int k=0; k<qsize; k++)
1063  {
1064  /* The communication flux has two parts:
1065  * - "diffusive" term containing sigma
1066  * - "advective" term representing usual upwind
1067  *
1068  * The calculation differs from the reference manual, since ad_coef and dif_coef have different meaning
1069  * than b and A in the manual.
1070  * In calculation of sigma there appears one more csection_lower in the denominator.
1071  */
1072  double sigma = frac_sigma[k][sbi]*arma::dot(dif_coef_edg[0][sbi][k]*fe_values_side.normal_vector(k),fe_values_side.normal_vector(k))*
1073  2*csection_higher[k]*csection_higher[k]/(csection_lower[k]*csection_lower[k]);
1074 
1075  // Since mm_coef_* contains cross section, we have to divide by it.
1076  double transport_flux = arma::dot(ad_coef_edg[1][sbi][k], fe_values_side.normal_vector(k));
1077  double por_lower_over_higher = mm_coef_lower[k]*csection_higher[k]/(mm_coef_higher[k]*csection_lower[k]);
1078 
1079  comm_flux[0][0] = (sigma-min(0.,transport_flux*por_lower_over_higher))*fv_sb[0]->JxW(k);
1080  comm_flux[0][1] = -(sigma-min(0.,transport_flux*por_lower_over_higher))*fv_sb[0]->JxW(k);
1081  comm_flux[1][0] = -(sigma+max(0.,transport_flux))*fv_sb[0]->JxW(k);
1082  comm_flux[1][1] = (sigma+max(0.,transport_flux))*fv_sb[0]->JxW(k);
1083 
1084  for (int n=0; n<2; n++)
1085  {
1086  if (!feo->dh()->el_is_local(element_id[n])) continue;
1087 
1088  for (unsigned int i=0; i<n_dofs[n]; i++)
1089  for (int m=0; m<2; m++)
1090  for (unsigned int j=0; j<n_dofs[m]; j++)
1091  local_matrix[(i+n*n_dofs[0])*(n_dofs[0]+n_dofs[1]) + m*n_dofs[0] + j] +=
1092  comm_flux[m][n]*fv_sb[m]->shape_value(j,k)*fv_sb[n]->shape_value(i,k);
1093  }
1094  }
1095  ls[sbi]->mat_set_values(n_dofs[0]+n_dofs[1], (int *)side_dof_indices, n_dofs[0]+n_dofs[1], (int *)side_dof_indices, local_matrix);
1096  }
1097  }
1098 
1099 }
1100 
1101 
1102 
1103 
1104 
1105 
1106 template<class Model>
1107 void TransportDG<Model>::set_boundary_conditions()
1108 {
1109  START_TIMER("assemble_bc");
1110  set_boundary_conditions<1>();
1111  set_boundary_conditions<2>();
1112  set_boundary_conditions<3>();
1113  END_TIMER("assemble_bc");
1114 }
1115 
1116 
1117 template<class Model>
1118 template<unsigned int dim>
1119 void TransportDG<Model>::set_boundary_conditions()
1120 {
1121  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1122  update_values | update_gradients | update_normal_vectors | update_side_JxW_values | update_quadrature_points);
1123  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1124  update_values);
1125  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
1126  unsigned int side_dof_indices[ndofs];
1127  double local_rhs[ndofs];
1128  vector<arma::vec> bc_values(qsize), bc_fluxes(qsize), bc_sigma(qsize);
1129  vector<arma::vec3> velocity;
1130 
1131  for (int i=0; i<qsize; i++)
1132  {
1133  bc_values[i].resize(n_subst_);
1134  bc_fluxes[i].resize(n_subst_);
1135  bc_sigma[i].resize(n_subst_);
1136  }
1137 
1138  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
1139  {
1140  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
1141  if (edg->n_sides > 1) continue;
1142  if (edg->side(0)->dim() != dim-1) continue;
1143  // skip edges lying not on the boundary
1144  if (edg->side(0)->cond() == NULL) continue;
1145 
1146  SideIter side = edg->side(0);
1147  typename DOFHandlerBase::CellIterator cell = mesh().element.full_iter(side->element());
1148  ElementAccessor<3> ele_acc = side->cond()->element_accessor();
1149 
1150  arma::uvec bc_type = data_.bc_type.value(side->cond()->element()->centre(), ele_acc);
1151 
1152  fe_values_side.reinit(cell, side->el_idx());
1153  fsv_rt.reinit(cell, side->el_idx());
1154  calculate_velocity(cell, velocity, fsv_rt);
1155 
1156  Model::compute_advection_diffusion_coefficients(fe_values_side.point_list(), velocity, side->element()->element_accessor(), ad_coef, dif_coef);
1157  Model::compute_dirichlet_bc(fe_values_side.point_list(), ele_acc, bc_values);
1158  data_.bc_flux.value_list(fe_values_side.point_list(), ele_acc, bc_fluxes);
1159  data_.bc_robin_sigma.value_list(fe_values_side.point_list(), ele_acc, bc_sigma);
1160 
1161  feo->dh()->get_dof_indices(cell, side_dof_indices);
1162 
1163  for (int sbi=0; sbi<n_subst_; sbi++)
1164  {
1165  for (unsigned int i=0; i<ndofs; i++) local_rhs[i] = 0;
1166 
1167  for (unsigned int k=0; k<qsize; k++)
1168  {
1169  double bc_term = 0;
1170  arma::vec3 bc_grad;
1171  bc_grad.zeros();
1172  if ((bc_type[sbi] == EqData::inflow && mh_dh->side_flux( *edg->side(0) ) < -mh_dh->precision())
1173  || (bc_type[sbi] == EqData::dirichlet))
1174  {
1175  bc_term = gamma[sbi][side->cond_idx()]*bc_values[k][sbi]*fe_values_side.JxW(k);
1176  bc_grad = -bc_values[k][sbi]*fe_values_side.JxW(k)*dg_variant*(arma::trans(dif_coef[sbi][k])*fe_values_side.normal_vector(k));
1177  }
1178  else if (bc_type[sbi] == EqData::neumann)
1179  {
1180  bc_term = bc_fluxes[k][sbi]*fe_values_side.JxW(k);
1181  }
1182  else if (bc_type[sbi] == EqData::robin)
1183  {
1184  bc_term = bc_sigma[k][sbi]*bc_values[k][sbi]*fe_values_side.JxW(k);
1185  }
1186 
1187  for (unsigned int i=0; i<ndofs; i++)
1188  local_rhs[i] += bc_term*fe_values_side.shape_value(i,k)
1189  + arma::dot(bc_grad,fe_values_side.shape_grad(i,k));
1190  }
1191  ls[sbi]->rhs_set_values(ndofs, (int *)side_dof_indices, local_rhs);
1192  }
1193  }
1194 }
1195 
1196 
1197 
1198 // TODO: The detection of side number from SideIter
1199 // in TransportDG::calculate_velocity()
1200 // should be done with help of class RefElement. This however requires
1201 // that the MH_DofHandler uses the node/side ordering defined in
1202 // the respective RefElement.
1203 template<class Model>
1204 template<unsigned int dim>
1205 void TransportDG<Model>::calculate_velocity(const typename DOFHandlerBase::CellIterator &cell, vector<arma::vec3> &velocity, FEValuesBase<dim,3> &fv)
1206 {
1207  std::map<const Node*, int> node_nums;
1208  for (unsigned int i=0; i<cell->n_nodes(); i++)
1209  node_nums[cell->node[i]] = i;
1210 
1211  velocity.resize(fv.n_points());
1212 
1213  for (unsigned int k=0; k<fv.n_points(); k++)
1214  {
1215  velocity[k].zeros();
1216  for (unsigned int sid=0; sid<cell->n_sides(); sid++)
1217  {
1218  if (cell->side(sid)->dim() != dim-1) continue;
1219  int num = dim*(dim+1)/2;
1220  for (unsigned int i=0; i<cell->side(sid)->n_nodes(); i++)
1221  num -= node_nums[cell->side(sid)->node(i)];
1222  velocity[k] += fv.shape_vector(num,k) * mh_dh->side_flux( *(cell->side(sid)) );
1223  }
1224  }
1225 }
1226 
1227 
1228 
1229 
1230 
1231 template<class Model>
1232 void TransportDG<Model>::set_DG_parameters_edge(const Edge &edg,
1233  const int s1,
1234  const int s2,
1235  const int K_size,
1236  const vector<arma::mat33> &K1,
1237  const vector<arma::mat33> &K2,
1238  const vector<double> &fluxes,
1239  const arma::vec3 &normal_vector,
1240  const double alpha1,
1241  const double alpha2,
1242  double &gamma,
1243  double *omega,
1244  double &transport_flux)
1245 {
1246  double delta[2];
1247  double h = 0;
1248  double local_alpha = 0;
1249 
1250  ASSERT(edg.side(s1)->valid(), "Invalid side of an edge.");
1251  SideIter s = edg.side(s1);
1252 
1253  // calculate the side diameter
1254  if (s->dim() == 0)
1255  {
1256  h = 1;
1257  }
1258  else
1259  {
1260  for (unsigned int i=0; i<s->n_nodes(); i++)
1261  for (unsigned int j=i+1; j<s->n_nodes(); j++)
1262  h = max(h, s->node(i)->distance(*s->node(j)));
1263  }
1264 
1265  // calculate the total in- and out-flux through the edge
1266  double pflux = 0, nflux = 0;
1267  for (int i=0; i<edg.n_sides; i++)
1268  {
1269  if (fluxes[i] > 0)
1270  pflux += fluxes[i];
1271  else
1272  nflux += fluxes[i];
1273  }
1274 
1275  // calculate the flux from s1 to s2
1276  if (fluxes[s2] > 0 && fluxes[s1] < 0 && s1 < s2)
1277  transport_flux = fluxes[s1]*fabs(fluxes[s2]/pflux);
1278  else if (fluxes[s2] < 0 && fluxes[s1] > 0 && s1 < s2)
1279  transport_flux = fluxes[s1]*fabs(fluxes[s2]/nflux);
1280  else if (s1==s2)
1281  transport_flux = fluxes[s1];
1282  else
1283  transport_flux = 0;
1284 
1285  gamma = 0.5*fabs(transport_flux);
1286 
1287 
1288  // determine local DG penalty
1289  local_alpha = max(alpha1, alpha2);
1290 
1291  if (s1 == s2)
1292  {
1293  omega[0] = 1;
1294 
1295  // delta is set to the average value of Kn.n on the side
1296  delta[0] = 0;
1297  for (unsigned int k=0; k<K_size; k++)
1298  delta[0] += dot(K1[k]*normal_vector,normal_vector);
1299  delta[0] /= K_size;
1300 
1301  gamma += local_alpha/h*(delta[0]);
1302  }
1303  else
1304  {
1305  delta[0] = 0;
1306  delta[1] = 0;
1307  for (unsigned int k=0; k<K_size; k++)
1308  {
1309  delta[0] += dot(K1[k]*normal_vector,normal_vector);
1310  delta[1] += dot(K2[k]*normal_vector,normal_vector);
1311  }
1312  delta[0] /= K_size;
1313  delta[1] /= K_size;
1314 
1315  double delta_sum = delta[0] + delta[1];
1316 
1317  if (delta_sum > numeric_limits<double>::epsilon())
1318  {
1319  omega[0] = delta[1]/delta_sum;
1320  omega[1] = delta[0]/delta_sum;
1321  gamma += local_alpha/h*(delta[0]*delta[1]/delta_sum);
1322  }
1323  else
1324  for (int i=0; i<2; i++) omega[i] = 0;
1325  }
1326 }
1327 
1328 
1329 
1330 
1331 
1332 
1333 template<class Model>
1334 void TransportDG<Model>::set_DG_parameters_boundary(const SideIter side,
1335  const int K_size,
1336  const vector<arma::mat33> &K,
1337  const double flux,
1338  const arma::vec3 &normal_vector,
1339  const double alpha,
1340  double &gamma)
1341 {
1342  double delta = 0, h = 0;
1343 
1344  // calculate the side diameter
1345  if (side->dim() == 0)
1346  {
1347  h = 1;
1348  }
1349  else
1350  {
1351  for (unsigned int i=0; i<side->n_nodes(); i++)
1352  for (unsigned int j=i+1; j<side->n_nodes(); j++)
1353  h = max(h, side->node(i)->distance( *side->node(j) ));
1354  }
1355 
1356  // delta is set to the average value of Kn.n on the side
1357  for (unsigned int k=0; k<K_size; k++)
1358  delta += dot(K[k]*normal_vector,normal_vector);
1359  delta /= K_size;
1360 
1361  gamma = 0.5*fabs(flux) + alpha/h*delta;
1362 }
1363 
1364 
1365 
1366 
1367 
1368 template<class Model>
1369 void TransportDG<Model>::set_initial_condition()
1370 {
1371  START_TIMER("set_init_cond");
1372  for (int sbi=0; sbi<n_subst_; sbi++)
1373  ls[sbi]->start_allocation();
1374  prepare_initial_condition<1>();
1375  prepare_initial_condition<2>();
1376  prepare_initial_condition<3>();
1377 
1378  for (int sbi=0; sbi<n_subst_; sbi++)
1379  ls[sbi]->start_add_assembly();
1380  prepare_initial_condition<1>();
1381  prepare_initial_condition<2>();
1382  prepare_initial_condition<3>();
1383 
1384  for (int sbi=0; sbi<n_subst_; sbi++)
1385  {
1386  ls[sbi]->finish_assembly();
1387  ls[sbi]->solve();
1388  }
1389  END_TIMER("set_init_cond");
1390 }
1391 
1392 template<class Model>
1393 template<unsigned int dim>
1394 void TransportDG<Model>::prepare_initial_condition()
1395 {
1396  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
1397  update_values | update_JxW_values | update_quadrature_points);
1398  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
1399  unsigned int dof_indices[ndofs];
1400  double matrix[ndofs*ndofs], rhs[ndofs];
1401  std::vector<arma::vec> init_values(qsize);
1402 
1403  for (int i=0; i<qsize; i++)
1404  init_values[i].resize(n_subst_);
1405 
1406  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
1407  {
1408  typename DOFHandlerBase::CellIterator elem = mesh_->element(feo->dh()->el_index(i_cell));
1409  if (elem->dim() != dim) continue;
1410 
1411  ElementAccessor<3> ele_acc = elem->element_accessor();
1412  feo->dh()->get_dof_indices(elem, dof_indices);
1413  fe_values.reinit(elem);
1414 
1415  Model::compute_init_cond(fe_values.point_list(), ele_acc, init_values);
1416 
1417  for (int sbi=0; sbi<n_subst_; sbi++)
1418  {
1419  for (unsigned int i=0; i<ndofs; i++)
1420  {
1421  rhs[i] = 0;
1422  for (unsigned int j=0; j<ndofs; j++)
1423  matrix[i*ndofs+j] = 0;
1424  }
1425 
1426  for (unsigned int k=0; k<qsize; k++)
1427  {
1428  double rhs_term = init_values[k](sbi)*fe_values.JxW(k);
1429 
1430  for (unsigned int i=0; i<ndofs; i++)
1431  {
1432  for (unsigned int j=0; j<ndofs; j++)
1433  matrix[i*ndofs+j] += fe_values.shape_value(i,k)*fe_values.shape_value(j,k)*fe_values.JxW(k);
1434 
1435  rhs[i] += fe_values.shape_value(i,k)*rhs_term;
1436  }
1437  }
1438  ls[sbi]->set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, matrix, rhs);
1439  }
1440  }
1441 }
1442 
1443 
1444 template<class Model>
1445 void TransportDG<Model>::calc_fluxes(vector<vector<double> > &bcd_balance, vector<vector<double> > &bcd_plus_balance, vector<vector<double> > &bcd_minus_balance)
1446 {
1447  calc_fluxes<1>(bcd_balance, bcd_plus_balance, bcd_minus_balance);
1448  calc_fluxes<2>(bcd_balance, bcd_plus_balance, bcd_minus_balance);
1449  calc_fluxes<3>(bcd_balance, bcd_plus_balance, bcd_minus_balance);
1450 }
1451 
1452 template<class Model>
1453 template<unsigned int dim>
1454 void TransportDG<Model>::calc_fluxes(vector<vector<double> > &bcd_balance, vector<vector<double> > &bcd_plus_balance, vector<vector<double> > &bcd_minus_balance)
1455 {
1456  FESideValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1457  update_values | update_gradients | update_normal_vectors | update_side_JxW_values | update_quadrature_points);
1458  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1459  update_values);
1460  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
1461  unsigned int dof_indices[ndofs];
1462  vector<arma::vec3> side_velocity(qsize);
1463  vector<arma::vec> bc_values(qsize);
1464  arma::vec3 conc_grad;
1465 
1466  for (int i=0; i<qsize; i++)
1467  bc_values[i].resize(n_subst_);
1468 
1469  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
1470  {
1471  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
1472  if (edg->n_sides > 1) continue;
1473  if (edg->side(0)->dim() != dim-1) continue;
1474  // skip edges lying not on the boundary
1475  if (edg->side(0)->cond() == NULL) continue;
1476 
1477  SideIter side = edg->side(0);
1478  DOFHandlerMultiDim::CellIterator cell = side->element();
1479  ElementAccessor<3> ele_acc = cell->element_accessor();
1480  Region r = side->cond()->element_accessor().region();
1481  if (! r.is_valid()) xprintf(Msg, "Invalid region, ele % d, edg: % d\n", side->cond()->bc_ele_idx_, side->cond()->edge_idx_);
1482  unsigned int bc_region_idx = r.boundary_idx();
1483 
1484  fe_values.reinit(cell, side->el_idx());
1485  fsv_rt.reinit(cell, side->el_idx());
1486  feo->dh()->get_dof_indices(cell, dof_indices);
1487  calculate_velocity(cell, side_velocity, fsv_rt);
1488  Model::compute_advection_diffusion_coefficients(fe_values.point_list(), side_velocity, ele_acc, ad_coef, dif_coef);
1489  Model::compute_dirichlet_bc(fe_values.point_list(), side->cond()->element_accessor(), bc_values);
1490  arma::uvec bc_type = data_.bc_type.value(side->cond()->element()->centre(), side->cond()->element_accessor());
1491 
1492  for (int sbi=0; sbi<n_subst_; sbi++)
1493  {
1494  double mass_flux = 0;
1495  double water_flux = 0;
1496  for (unsigned int k=0; k<qsize; k++)
1497  water_flux += arma::dot(ad_coef[sbi][k], fe_values.normal_vector(k))*fe_values.JxW(k);
1498  water_flux /= side->measure();
1499 
1500  for (unsigned int k=0; k<qsize; k++)
1501  {
1502  double conc = 0;
1503  conc_grad.zeros();
1504  for (unsigned int i=0; i<ndofs; i++)
1505  {
1506  conc += fe_values.shape_value(i,k)*ls[sbi]->get_solution_array()[dof_indices[i]-feo->dh()->loffset()];
1507  conc_grad += fe_values.shape_grad(i,k)*ls[sbi]->get_solution_array()[dof_indices[i]-feo->dh()->loffset()];
1508  }
1509 
1510  // flux due to advection
1511  mass_flux += water_flux*conc*fe_values.JxW(k);
1512 
1513  if (bc_type[sbi] == EqData::dirichlet || (bc_type[sbi] == EqData::inflow && water_flux*side->measure() < -mh_dh->precision()))
1514  {
1515  // flux due to diffusion
1516  mass_flux -= arma::dot(dif_coef[sbi][k]*conc_grad,fe_values.normal_vector(k))*fe_values.JxW(k);
1517 
1518  // the penalty term has to be added otherwise the mass balance will not hold
1519  mass_flux -= gamma[sbi][side->cond_idx()]*(bc_values[k][sbi] - conc)*fe_values.JxW(k);
1520  }
1521 
1522  }
1523  bcd_balance[sbi][bc_region_idx] += mass_flux;
1524  if (mass_flux >= 0) bcd_plus_balance[sbi][bc_region_idx] += mass_flux;
1525  else bcd_minus_balance[sbi][bc_region_idx] += mass_flux;
1526  }
1527  }
1528 
1529 }
1530 
1531 template<class Model>
1532 void TransportDG<Model>::calc_elem_sources(vector<vector<double> > &mass, vector<vector<double> > &src_balance)
1533 {
1534  calc_elem_sources<1>(mass, src_balance);
1535  calc_elem_sources<2>(mass, src_balance);
1536  calc_elem_sources<3>(mass, src_balance);
1537 }
1538 
1539 template<class Model>
1540 template<unsigned int dim>
1541 void TransportDG<Model>::calc_elem_sources(vector<vector<double> > &mass, vector<vector<double> > &src_balance)
1542 {
1543  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
1544  update_values | update_JxW_values | update_quadrature_points);
1545  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
1546  unsigned int dof_indices[ndofs];
1547  vector<arma::vec> sources_conc(qsize), sources_density(qsize), sources_sigma(qsize);
1548  double mass_sum, sources_sum, conc, conc_diff;
1549 
1550  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
1551  {
1552  typename DOFHandlerBase::CellIterator elem = mesh_->element(feo->dh()->el_index(i_cell));
1553  if (elem->dim() != dim) continue;
1554 
1555  Region r = elem->element_accessor().region();
1556  if (! r.is_valid()) xprintf(Msg, "Invalid region, ele % d\n", elem.index());
1557  unsigned int region_idx = r.bulk_idx();
1558 
1559  ElementAccessor<3> ele_acc = elem->element_accessor();
1560  fe_values.reinit(elem);
1561  feo->dh()->get_dof_indices(elem, dof_indices);
1562 
1563  Model::compute_mass_matrix_coefficient(fe_values.point_list(), ele_acc, mm_coef);
1564  Model::compute_source_coefficients(fe_values.point_list(), ele_acc, sources_conc, sources_density, sources_sigma);
1565 
1566  for (int sbi=0; sbi<n_subst_; sbi++)
1567  {
1568  mass_sum = 0;
1569  sources_sum = 0;
1570 
1571  for (unsigned int k=0; k<qsize; k++)
1572  {
1573  conc = 0;
1574  for (unsigned int i=0; i<ndofs; i++)
1575  conc += fe_values.shape_value(i,k)*ls[sbi]->get_solution_array()[dof_indices[i]-feo->dh()->loffset()];
1576 
1577  mass_sum += mm_coef[k]*conc*fe_values.JxW(k);
1578 
1579  conc_diff = sources_conc[k][sbi] - conc;
1580  sources_sum += (sources_density[k][sbi] + conc_diff*sources_sigma[k][sbi])*fe_values.JxW(k);
1581  }
1582 
1583  mass[sbi][region_idx] += mass_sum;
1584  src_balance[sbi][region_idx] += sources_sum;
1585  }
1586  }
1587 }
1588 
1589 
1590 
1591 
1592 
1593 template class TransportDG<ConcentrationTransportModel>;
1594 template class TransportDG<HeatTransferModel>;
1595 
1596 
1597 
1598 
Generated on Thu Apr 17 2014 01:28:45 for Flow123d by   doxygen 1.8.4