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 template<class Model>
536 void TransportDG<Model>::set_velocity_field(const MH_DofHandler &dh)
537 {
538  // So far the velocity_vector contains zeros, so we ignore it.
539  // Instead we use the value Side.flux.
540 
541  mh_dh = &dh;
542  Model::flux_changed = true;
543 
544 }
545 
546 
547 template<class Model>
548 void TransportDG<Model>::output_data()
549 {
550  double *solution;
551  unsigned int dof_indices[max(feo->fe<1>()->n_dofs(), max(feo->fe<2>()->n_dofs(), feo->fe<3>()->n_dofs()))];
552 
553  if (!time_->is_current( time_->marks().type_output() )) return;
554 
555  START_TIMER("DG-OUTPUT");
556 
557  // gather the solution from all processors
558  output_vector_gather();
559  data_.output_fields.set_time(*time_);
560  data_.output_fields.output(output_stream);
561  output_stream->write_time_frame();
562 
563  if (mass_balance() != NULL)
564  mass_balance()->output(time_->t());
565 
566  END_TIMER("DG-OUTPUT");
567 }
568 
569 
570 template<class Model>
571 void TransportDG<Model>::assemble_mass_matrix()
572 {
573  START_TIMER("assemble_mass");
574  assemble_mass_matrix<1>();
575  assemble_mass_matrix<2>();
576  assemble_mass_matrix<3>();
577  END_TIMER("assemble_mass");
578 }
579 
580 
581 template<class Model> template<unsigned int dim>
582 void TransportDG<Model>::assemble_mass_matrix()
583 {
584  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(), update_values | update_JxW_values | update_quadrature_points);
585  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
586  unsigned int dof_indices[ndofs];
587  PetscScalar local_mass_matrix[ndofs*ndofs];
588 
589  // assemble integral over elements
590  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
591  {
592  typename DOFHandlerBase::CellIterator cell = mesh_->element(feo->dh()->el_index(i_cell));
593  if (cell->dim() != dim) continue;
594 
595  fe_values.reinit(cell);
596  feo->dh()->get_dof_indices(cell, dof_indices);
597  ElementAccessor<3> ele_acc = cell->element_accessor();
598 
599  Model::compute_mass_matrix_coefficient(fe_values.point_list(), ele_acc, mm_coef);
600 
601  // assemble the local stiffness and mass matrix
602  for (unsigned int i=0; i<ndofs; i++)
603  {
604  for (unsigned int j=0; j<ndofs; j++)
605  {
606  local_mass_matrix[i*ndofs+j] = 0;
607  for (unsigned int k=0; k<qsize; k++)
608  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);
609  }
610  }
611 
612  ls_dt->mat_set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, local_mass_matrix);
613  }
614 }
615 
616 
617 
618 
619 template<class Model>
620 void TransportDG<Model>::assemble_stiffness_matrix()
621 {
622  START_TIMER("assemble_stiffness");
623  START_TIMER("assemble_volume_integrals");
624  assemble_volume_integrals<1>();
625  assemble_volume_integrals<2>();
626  assemble_volume_integrals<3>();
627  END_TIMER("assemble_volume_integrals");
628 
629  START_TIMER("assemble_fluxes_boundary");
630  assemble_fluxes_boundary<1>();
631  assemble_fluxes_boundary<2>();
632  assemble_fluxes_boundary<3>();
633  END_TIMER("assemble_fluxes_boundary");
634 
635  START_TIMER("assemble_fluxes_elem_elem");
636  assemble_fluxes_element_element<1>();
637  assemble_fluxes_element_element<2>();
638  assemble_fluxes_element_element<3>();
639  END_TIMER("assemble_fluxes_elem_elem");
640 
641  START_TIMER("assemble_fluxes_elem_side");
642  assemble_fluxes_element_side<1>();
643  assemble_fluxes_element_side<2>();
644  assemble_fluxes_element_side<3>();
645  END_TIMER("assemble_fluxes_elem_side");
646  END_TIMER("assemble_stiffness");
647 }
648 
649 
650 
651 template<class Model>
652 template<unsigned int dim>
653 void TransportDG<Model>::assemble_volume_integrals()
654 {
655  FEValues<dim,3> fv_rt(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe_rt<dim>(),
656  update_values | update_gradients);
657  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
658  update_values | update_gradients | update_JxW_values | update_quadrature_points);
659  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
660  unsigned int dof_indices[ndofs];
661  vector<arma::vec3> velocity(qsize);
662  vector<arma::vec> sources_sigma(qsize);
663  PetscScalar local_matrix[ndofs*ndofs];
664 
665  // assemble integral over elements
666  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
667  {
668  typename DOFHandlerBase::CellIterator cell = mesh_->element(feo->dh()->el_index(i_cell));
669  if (cell->dim() != dim) continue;
670 
671  fe_values.reinit(cell);
672  fv_rt.reinit(cell);
673  ElementAccessor<3> ele_acc = cell->element_accessor();
674  feo->dh()->get_dof_indices(cell, dof_indices);
675 
676  calculate_velocity(cell, velocity, fv_rt);
677  Model::compute_advection_diffusion_coefficients(fe_values.point_list(), velocity, ele_acc, ad_coef, dif_coef);
678  Model::compute_sources_sigma(fe_values.point_list(), ele_acc, sources_sigma);
679 
680  // assemble the local stiffness matrix
681  for (int sbi=0; sbi<n_subst_; sbi++)
682  {
683  for (unsigned int i=0; i<ndofs; i++)
684  for (unsigned int j=0; j<ndofs; j++)
685  local_matrix[i*ndofs+j] = 0;
686 
687  for (unsigned int k=0; k<qsize; k++)
688  {
689  for (unsigned int i=0; i<ndofs; i++)
690  {
691  arma::vec3 Kt_grad_i = dif_coef[sbi][k].t()*fe_values.shape_grad(i,k);
692  double ad_dot_grad_i = arma::dot(ad_coef[sbi][k], fe_values.shape_grad(i,k));
693 
694  for (unsigned int j=0; j<ndofs; j++)
695  local_matrix[i*ndofs+j] += (arma::dot(Kt_grad_i, fe_values.shape_grad(j,k))
696  -fe_values.shape_value(j,k)*ad_dot_grad_i
697  +sources_sigma[k][sbi]*fe_values.shape_value(j,k)*fe_values.shape_value(i,k))*fe_values.JxW(k);
698  }
699  }
700  ls[sbi]->mat_set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, local_matrix);
701  }
702  }
703 }
704 
705 
706 template<class Model>
707 void TransportDG<Model>::set_sources()
708 {
709  START_TIMER("assemble_sources");
710  set_sources<1>();
711  set_sources<2>();
712  set_sources<3>();
713  END_TIMER("assemble_sources");
714 }
715 
716 template<class Model>
717 template<unsigned int dim>
718 void TransportDG<Model>::set_sources()
719 {
720  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
721  update_values | update_JxW_values | update_quadrature_points);
722  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
723  vector<arma::vec> sources_conc(qsize), sources_density(qsize), sources_sigma(qsize);
724  unsigned int dof_indices[ndofs];
725  PetscScalar local_rhs[ndofs];
726  double source;
727 
728  // assemble integral over elements
729  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
730  {
731  typename DOFHandlerBase::CellIterator cell = mesh_->element(feo->dh()->el_index(i_cell));
732  if (cell->dim() != dim) continue;
733 
734  fe_values.reinit(cell);
735  feo->dh()->get_dof_indices(cell, dof_indices);
736 
737  Model::compute_source_coefficients(fe_values.point_list(), cell->element_accessor(), sources_conc, sources_density, sources_sigma);
738 
739  // assemble the local stiffness matrix
740  for (int sbi=0; sbi<n_subst_; sbi++)
741  {
742  for (unsigned int i=0; i<ndofs; i++)
743  local_rhs[i] = 0;
744 
745  // compute sources
746  for (unsigned int k=0; k<qsize; k++)
747  {
748  source = (sources_density[k][sbi] + sources_conc[k][sbi]*sources_sigma[k][sbi])*fe_values.JxW(k);
749 
750  for (unsigned int i=0; i<ndofs; i++)
751  local_rhs[i] += source*fe_values.shape_value(i,k);
752  }
753  ls[sbi]->rhs_set_values(ndofs, (int *)dof_indices, local_rhs);
754  }
755  }
756 }
757 
758 
759 
760 template<class Model>
761 template<unsigned int dim>
762 void TransportDG<Model>::assemble_fluxes_element_element()
763 {
764  vector<FESideValues<dim,3>*> fe_values;
765  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
766  update_values);
767  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size(),
768  n_max_sides = ad_coef_edg.size();
769  vector<unsigned int*> side_dof_indices;
770  PetscScalar local_matrix[ndofs*ndofs];
771  vector<vector<arma::vec3> > side_velocity(n_max_sides);
772  vector<arma::vec> dg_penalty(n_max_sides);
773  double gamma_l, omega[2], transport_flux;
774 
775  for (int sid=0; sid<n_max_sides; sid++)
776  {
777  side_dof_indices.push_back(new unsigned int[ndofs]);
778  fe_values.push_back(new FESideValues<dim,3>(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
779  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points));
780  }
781 
782  // assemble integral over sides
783  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
784  {
785  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
786  if (edg->n_sides < 2 || edg->side(0)->element()->dim() != dim) continue;
787 
788  for (int sid=0; sid<edg->n_sides; sid++)
789  {
790  typename DOFHandlerBase::CellIterator cell = edg->side(sid)->element();
791  ElementAccessor<3> ele_acc = cell->element_accessor();
792  feo->dh()->get_dof_indices(cell, side_dof_indices[sid]);
793  fe_values[sid]->reinit(cell, edg->side(sid)->el_idx());
794  fsv_rt.reinit(cell, edg->side(sid)->el_idx());
795  calculate_velocity(cell, side_velocity[sid], fsv_rt);
796  Model::compute_advection_diffusion_coefficients(fe_values[sid]->point_list(), side_velocity[sid], ele_acc, ad_coef_edg[sid], dif_coef_edg[sid]);
797  dg_penalty[sid] = data_.dg_penalty.value(cell->centre(), ele_acc);
798  }
799 
800  // fluxes and penalty
801  for (int sbi=0; sbi<n_subst_; sbi++)
802  {
803  vector<double> fluxes(edg->n_sides);
804  for (unsigned int sid=0; sid<edg->n_sides; sid++)
805  {
806  fluxes[sid] = 0;
807  for (unsigned int k=0; k<qsize; k++)
808  fluxes[sid] += arma::dot(ad_coef_edg[sid][sbi][k], fe_values[sid]->normal_vector(k))*fe_values[sid]->JxW(k);
809  fluxes[sid] /= edg->side(sid)->measure();
810  }
811 
812  for (int s1=0; s1<edg->n_sides; s1++)
813  {
814  for (int s2=s1+1; s2<edg->n_sides; s2++)
815  {
816  ASSERT(edg->side(s1)->valid(), "Invalid side of edge.");
817  if (!feo->dh()->el_is_local(edg->side(s1)->element().index())
818  && !feo->dh()->el_is_local(edg->side(s2)->element().index())) continue;
819 
820  arma::vec3 nv = fe_values[s1]->normal_vector(0);
821 
822  // set up the parameters for DG method
823  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);
824 
825  int sd[2];
826  sd[0] = s1;
827  sd[1] = s2;
828 
829 #define AVERAGE(i,k,side_id) (fe_values[sd[side_id]]->shape_value(i,k)*0.5)
830 #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])
831 #define JUMP(i,k,side_id) ((side_id==0?1:-1)*fe_values[sd[side_id]]->shape_value(i,k))
832 
833  // For selected pair of elements:
834  for (int n=0; n<2; n++)
835  {
836  if (!feo->dh()->el_is_local(edg->side(sd[n])->element().index())) continue;
837 
838  for (int m=0; m<2; m++)
839  {
840  for (unsigned int i=0; i<fe_values[sd[n]]->n_dofs(); i++)
841  for (unsigned int j=0; j<fe_values[sd[m]]->n_dofs(); j++)
842  local_matrix[i*fe_values[sd[m]]->n_dofs()+j] = 0;
843 
844  for (unsigned int k=0; k<qsize; k++)
845  {
846  double flux_times_JxW = transport_flux*fe_values[0]->JxW(k);
847  double gamma_times_JxW = gamma_l*fe_values[0]->JxW(k);
848 
849  for (unsigned int i=0; i<fe_values[sd[n]]->n_dofs(); i++)
850  {
851  double flux_JxW_jump_i = flux_times_JxW*JUMP(i,k,n);
852  double gamma_JxW_jump_i = gamma_times_JxW*JUMP(i,k,n);
853  double JxW_jump_i = fe_values[0]->JxW(k)*JUMP(i,k,n);
854  double JxW_var_wavg_i = fe_values[0]->JxW(k)*WAVERAGE(i,k,n)*dg_variant;
855 
856  for (unsigned int j=0; j<fe_values[sd[m]]->n_dofs(); j++)
857  {
858  int index = i*fe_values[sd[m]]->n_dofs()+j;
859 
860  // flux due to transport (applied on interior edges) (average times jump)
861  local_matrix[index] += flux_JxW_jump_i*AVERAGE(j,k,m);
862 
863  // penalty enforcing continuity across edges (applied on interior and Dirichlet edges) (jump times jump)
864  local_matrix[index] += gamma_JxW_jump_i*JUMP(j,k,m);
865 
866  // terms due to diffusion
867  local_matrix[index] -= WAVERAGE(j,k,m)*JxW_jump_i;
868  local_matrix[index] -= JUMP(j,k,m)*JxW_var_wavg_i;
869  }
870  }
871  }
872  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);
873  }
874  }
875 #undef AVERAGE
876 #undef WAVERAGE
877 #undef JUMP
878  }
879  }
880  }
881  }
882 
883  for (unsigned int i=0; i<n_max_sides; i++)
884  {
885  delete fe_values[i];
886  delete[] side_dof_indices[i];
887  }
888 }
889 
890 
891 template<class Model>
892 template<unsigned int dim>
893 void TransportDG<Model>::assemble_fluxes_boundary()
894 {
895  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
896  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points);
897  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
898  update_values);
899  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
900  unsigned int side_dof_indices[ndofs];
901  PetscScalar local_matrix[ndofs*ndofs];
902  vector<arma::vec3> side_velocity;
903  vector<arma::vec> robin_sigma(qsize);
904  arma::vec dg_penalty;
905  double gamma_l;
906 
907  // assemble boundary integral
908  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
909  {
910  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
911  if (edg->n_sides > 1) continue;
912  // check spatial dimension
913  if (edg->side(0)->dim() != dim-1) continue;
914  // skip edges lying not on the boundary
915  if (edg->side(0)->cond() == NULL) continue;
916 
917  SideIter side = edg->side(0);
918  typename DOFHandlerBase::CellIterator cell = side->element();
919  ElementAccessor<3> ele_acc = cell->element_accessor();
920  feo->dh()->get_dof_indices(cell, side_dof_indices);
921  fe_values_side.reinit(cell, side->el_idx());
922  fsv_rt.reinit(cell, side->el_idx());
923 
924  calculate_velocity(cell, side_velocity, fsv_rt);
925  Model::compute_advection_diffusion_coefficients(fe_values_side.point_list(), side_velocity, ele_acc, ad_coef, dif_coef);
926  dg_penalty = data_.dg_penalty.value(cell->centre(), ele_acc);
927  arma::uvec bc_type = data_.bc_type.value(side->cond()->element()->centre(), side->cond()->element_accessor());
928  data_.bc_robin_sigma.value_list(fe_values_side.point_list(), side->cond()->element_accessor(), robin_sigma);
929 
930  for (int sbi=0; sbi<n_subst_; sbi++)
931  {
932  for (unsigned int i=0; i<ndofs; i++)
933  for (unsigned int j=0; j<ndofs; j++)
934  local_matrix[i*ndofs+j] = 0;
935 
936  // On Neumann boundaries we have only term from integrating by parts the advective term,
937  // on Dirichlet boundaries we additionally apply the penalty which enforces the prescribed value.
938  double side_flux = 0;
939  for (unsigned int k=0; k<qsize; k++)
940  side_flux += arma::dot(ad_coef[sbi][k], fe_values_side.normal_vector(k))*fe_values_side.JxW(k);
941  double transport_flux = side_flux/side->measure();
942 
943  if ((bc_type[sbi] == EqData::dirichlet) || (bc_type[sbi] == EqData::inflow ))
944  {
945  // set up the parameters for DG method
946  set_DG_parameters_boundary(side, qsize, dif_coef[sbi], transport_flux, fe_values_side.normal_vector(0), dg_penalty[sbi], gamma_l);
947  gamma[sbi][side->cond_idx()] = gamma_l;
948  if (bc_type[sbi] == EqData::dirichlet || side_flux < -mh_dh->precision())
949  transport_flux += gamma_l;
950  }
951 
952  // fluxes and penalty
953  for (unsigned int k=0; k<qsize; k++)
954  {
955  double flux_times_JxW;
956  if (bc_type[sbi] == EqData::robin)
957  flux_times_JxW = (transport_flux + robin_sigma[k][sbi])*fe_values_side.JxW(k);
958  else
959  flux_times_JxW = transport_flux*fe_values_side.JxW(k);
960 
961  for (unsigned int i=0; i<ndofs; i++)
962  for (unsigned int j=0; j<ndofs; j++)
963  {
964  // flux due to advection and penalty
965  local_matrix[i*ndofs+j] += flux_times_JxW*fe_values_side.shape_value(i,k)*fe_values_side.shape_value(j,k);
966 
967  // flux due to diffusion (only on dirichlet and inflow boundary)
968  if (bc_type[sbi] == EqData::dirichlet || (bc_type[sbi] == EqData::inflow && side_flux < -mh_dh->precision()))
969  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)
970  + 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
971  )*fe_values_side.JxW(k);
972  }
973  }
974  ls[sbi]->mat_set_values(ndofs, (int *)side_dof_indices, ndofs, (int *)side_dof_indices, local_matrix);
975  }
976  }
977 }
978 
979 
980 template<class Model>
981 template<unsigned int dim>
982 void TransportDG<Model>::assemble_fluxes_element_side()
983 {
984 
985  if (dim == 1) return;
986  FEValues<dim-1,3> fe_values_vb(*feo->mapping<dim-1>(), *feo->q<dim-1>(), *feo->fe<dim-1>(),
987  update_values | update_gradients | update_JxW_values | update_quadrature_points);
988  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
989  update_values | update_gradients | update_side_JxW_values | update_normal_vectors | update_quadrature_points);
990  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
991  update_values);
992  FEValues<dim-1,3> fv_rt(*feo->mapping<dim-1>(), *feo->q<dim-1>(), *feo->fe_rt<dim-1>(),
993  update_values);
994 
995  vector<FEValuesSpaceBase<3>*> fv_sb(2);
996  const unsigned int ndofs = feo->fe<dim>()->n_dofs(); // number of local dofs
997  const unsigned int qsize = feo->q<dim-1>()->size(); // number of quadrature points
998  unsigned int side_dof_indices[2*ndofs], n_dofs[2];
999  vector<arma::vec3> velocity_higher, velocity_lower;
1000  vector<arma::vec> frac_sigma(qsize);
1001  vector<double> csection_lower(qsize), csection_higher(qsize), mm_coef_lower(qsize), mm_coef_higher(qsize);
1002  PetscScalar local_matrix[4*ndofs*ndofs];
1003  double comm_flux[2][2];
1004 
1005  // index 0 = element with lower dimension,
1006  // index 1 = side of element with higher dimension
1007  fv_sb[0] = &fe_values_vb;
1008  fv_sb[1] = &fe_values_side;
1009 
1010  // assemble integral over sides
1011  for (int inb=0; inb<feo->dh()->n_loc_nb(); inb++)
1012  {
1013  Neighbour *nb = &mesh_->vb_neighbours_[feo->dh()->nb_index(inb)];
1014  // skip neighbours of different dimension
1015  if (nb->element()->dim() != dim-1) continue;
1016 
1017  typename DOFHandlerBase::CellIterator cell_sub = mesh_->element.full_iter(nb->element());
1018  feo->dh()->get_dof_indices(cell_sub, side_dof_indices);
1019  fe_values_vb.reinit(cell_sub);
1020  n_dofs[0] = fv_sb[0]->n_dofs();
1021 
1022  typename DOFHandlerBase::CellIterator cell = nb->side()->element();
1023  feo->dh()->get_dof_indices(cell, side_dof_indices+n_dofs[0]);
1024  fe_values_side.reinit(cell, nb->side()->el_idx());
1025  n_dofs[1] = fv_sb[1]->n_dofs();
1026 
1027  // Element id's for testing if they belong to local partition.
1028  int element_id[2];
1029  element_id[0] = cell_sub.index();
1030  element_id[1] = cell.index();
1031 
1032  fsv_rt.reinit(cell, nb->side()->el_idx());
1033  fv_rt.reinit(cell_sub);
1034  calculate_velocity(cell, velocity_higher, fsv_rt);
1035  calculate_velocity(cell_sub, velocity_lower, fv_rt);
1036  Model::compute_advection_diffusion_coefficients(fe_values_vb.point_list(), velocity_lower, cell_sub->element_accessor(), ad_coef_edg[0], dif_coef_edg[0]);
1037  Model::compute_advection_diffusion_coefficients(fe_values_vb.point_list(), velocity_higher, cell->element_accessor(), ad_coef_edg[1], dif_coef_edg[1]);
1038  Model::compute_mass_matrix_coefficient(fe_values_vb.point_list(), cell_sub->element_accessor(), mm_coef_lower);
1039  Model::compute_mass_matrix_coefficient(fe_values_vb.point_list(), cell->element_accessor(), mm_coef_higher);
1040  data_.cross_section.value_list(fe_values_vb.point_list(), cell_sub->element_accessor(), csection_lower);
1041  data_.cross_section.value_list(fe_values_vb.point_list(), cell->element_accessor(), csection_higher);
1042  data_.fracture_sigma.value_list(fe_values_vb.point_list(), cell_sub->element_accessor(), frac_sigma);
1043 
1044  for (int sbi=0; sbi<n_subst_; sbi++)
1045  {
1046  for (unsigned int i=0; i<n_dofs[0]+n_dofs[1]; i++)
1047  for (unsigned int j=0; j<n_dofs[0]+n_dofs[1]; j++)
1048  local_matrix[i*(n_dofs[0]+n_dofs[1])+j] = 0;
1049 
1050  // set transmission conditions
1051  for (unsigned int k=0; k<qsize; k++)
1052  {
1053  /* The communication flux has two parts:
1054  * - "diffusive" term containing sigma
1055  * - "advective" term representing usual upwind
1056  *
1057  * The calculation differs from the reference manual, since ad_coef and dif_coef have different meaning
1058  * than b and A in the manual.
1059  * In calculation of sigma there appears one more csection_lower in the denominator.
1060  */
1061  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))*
1062  2*csection_higher[k]*csection_higher[k]/(csection_lower[k]*csection_lower[k]);
1063 
1064  // Since mm_coef_* contains cross section, we have to divide by it.
1065  double transport_flux = arma::dot(ad_coef_edg[1][sbi][k], fe_values_side.normal_vector(k));
1066  double por_lower_over_higher = mm_coef_lower[k]*csection_higher[k]/(mm_coef_higher[k]*csection_lower[k]);
1067 
1068  comm_flux[0][0] = (sigma-min(0.,transport_flux*por_lower_over_higher))*fv_sb[0]->JxW(k);
1069  comm_flux[0][1] = -(sigma-min(0.,transport_flux*por_lower_over_higher))*fv_sb[0]->JxW(k);
1070  comm_flux[1][0] = -(sigma+max(0.,transport_flux))*fv_sb[0]->JxW(k);
1071  comm_flux[1][1] = (sigma+max(0.,transport_flux))*fv_sb[0]->JxW(k);
1072 
1073  for (int n=0; n<2; n++)
1074  {
1075  if (!feo->dh()->el_is_local(element_id[n])) continue;
1076 
1077  for (unsigned int i=0; i<n_dofs[n]; i++)
1078  for (int m=0; m<2; m++)
1079  for (unsigned int j=0; j<n_dofs[m]; j++)
1080  local_matrix[(i+n*n_dofs[0])*(n_dofs[0]+n_dofs[1]) + m*n_dofs[0] + j] +=
1081  comm_flux[m][n]*fv_sb[m]->shape_value(j,k)*fv_sb[n]->shape_value(i,k);
1082  }
1083  }
1084  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);
1085  }
1086  }
1087 
1088 }
1089 
1090 
1091 
1092 
1093 
1094 
1095 template<class Model>
1096 void TransportDG<Model>::set_boundary_conditions()
1097 {
1098  START_TIMER("assemble_bc");
1099  set_boundary_conditions<1>();
1100  set_boundary_conditions<2>();
1101  set_boundary_conditions<3>();
1102  END_TIMER("assemble_bc");
1103 }
1104 
1105 
1106 template<class Model>
1107 template<unsigned int dim>
1108 void TransportDG<Model>::set_boundary_conditions()
1109 {
1110  FESideValues<dim,3> fe_values_side(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1111  update_values | update_gradients | update_normal_vectors | update_side_JxW_values | update_quadrature_points);
1112  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1113  update_values);
1114  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
1115  unsigned int side_dof_indices[ndofs];
1116  double local_rhs[ndofs];
1117  vector<arma::vec> bc_values(qsize), bc_fluxes(qsize), bc_sigma(qsize);
1118  vector<arma::vec3> velocity;
1119 
1120  for (int i=0; i<qsize; i++)
1121  {
1122  bc_values[i].resize(n_subst_);
1123  bc_fluxes[i].resize(n_subst_);
1124  bc_sigma[i].resize(n_subst_);
1125  }
1126 
1127  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
1128  {
1129  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
1130  if (edg->n_sides > 1) continue;
1131  if (edg->side(0)->dim() != dim-1) continue;
1132  // skip edges lying not on the boundary
1133  if (edg->side(0)->cond() == NULL) continue;
1134 
1135  SideIter side = edg->side(0);
1136  typename DOFHandlerBase::CellIterator cell = mesh().element.full_iter(side->element());
1137  ElementAccessor<3> ele_acc = side->cond()->element_accessor();
1138 
1139  arma::uvec bc_type = data_.bc_type.value(side->cond()->element()->centre(), ele_acc);
1140 
1141  fe_values_side.reinit(cell, side->el_idx());
1142  fsv_rt.reinit(cell, side->el_idx());
1143  calculate_velocity(cell, velocity, fsv_rt);
1144 
1145  Model::compute_advection_diffusion_coefficients(fe_values_side.point_list(), velocity, side->element()->element_accessor(), ad_coef, dif_coef);
1146  Model::compute_dirichlet_bc(fe_values_side.point_list(), ele_acc, bc_values);
1147  data_.bc_flux.value_list(fe_values_side.point_list(), ele_acc, bc_fluxes);
1148  data_.bc_robin_sigma.value_list(fe_values_side.point_list(), ele_acc, bc_sigma);
1149 
1150  feo->dh()->get_dof_indices(cell, side_dof_indices);
1151 
1152  for (int sbi=0; sbi<n_subst_; sbi++)
1153  {
1154  for (unsigned int i=0; i<ndofs; i++) local_rhs[i] = 0;
1155 
1156  for (unsigned int k=0; k<qsize; k++)
1157  {
1158  double bc_term = 0;
1159  arma::vec3 bc_grad;
1160  bc_grad.zeros();
1161  if ((bc_type[sbi] == EqData::inflow && mh_dh->side_flux( *edg->side(0) ) < -mh_dh->precision())
1162  || (bc_type[sbi] == EqData::dirichlet))
1163  {
1164  bc_term = gamma[sbi][side->cond_idx()]*bc_values[k][sbi]*fe_values_side.JxW(k);
1165  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));
1166  }
1167  else if (bc_type[sbi] == EqData::neumann)
1168  {
1169  bc_term = bc_fluxes[k][sbi]*fe_values_side.JxW(k);
1170  }
1171  else if (bc_type[sbi] == EqData::robin)
1172  {
1173  bc_term = bc_sigma[k][sbi]*bc_values[k][sbi]*fe_values_side.JxW(k);
1174  }
1175 
1176  for (unsigned int i=0; i<ndofs; i++)
1177  local_rhs[i] += bc_term*fe_values_side.shape_value(i,k)
1178  + arma::dot(bc_grad,fe_values_side.shape_grad(i,k));
1179  }
1180  ls[sbi]->rhs_set_values(ndofs, (int *)side_dof_indices, local_rhs);
1181  }
1182  }
1183 }
1184 
1185 
1186 
1187 // TODO: The detection of side number from SideIter
1188 // in TransportDG::calculate_velocity()
1189 // should be done with help of class RefElement. This however requires
1190 // that the MH_DofHandler uses the node/side ordering defined in
1191 // the respective RefElement.
1192 template<class Model>
1193 template<unsigned int dim>
1194 void TransportDG<Model>::calculate_velocity(const typename DOFHandlerBase::CellIterator &cell, vector<arma::vec3> &velocity, FEValuesBase<dim,3> &fv)
1195 {
1196  std::map<const Node*, int> node_nums;
1197  for (unsigned int i=0; i<cell->n_nodes(); i++)
1198  node_nums[cell->node[i]] = i;
1199 
1200  velocity.resize(fv.n_points());
1201 
1202  for (unsigned int k=0; k<fv.n_points(); k++)
1203  {
1204  velocity[k].zeros();
1205  for (unsigned int sid=0; sid<cell->n_sides(); sid++)
1206  {
1207  if (cell->side(sid)->dim() != dim-1) continue;
1208  int num = dim*(dim+1)/2;
1209  for (unsigned int i=0; i<cell->side(sid)->n_nodes(); i++)
1210  num -= node_nums[cell->side(sid)->node(i)];
1211  velocity[k] += fv.shape_vector(num,k) * mh_dh->side_flux( *(cell->side(sid)) );
1212  }
1213  }
1214 }
1215 
1216 
1217 
1218 
1219 
1220 template<class Model>
1221 void TransportDG<Model>::set_DG_parameters_edge(const Edge &edg,
1222  const int s1,
1223  const int s2,
1224  const int K_size,
1225  const vector<arma::mat33> &K1,
1226  const vector<arma::mat33> &K2,
1227  const vector<double> &fluxes,
1228  const arma::vec3 &normal_vector,
1229  const double alpha1,
1230  const double alpha2,
1231  double &gamma,
1232  double *omega,
1233  double &transport_flux)
1234 {
1235  double delta[2];
1236  double h = 0;
1237  double local_alpha = 0;
1238 
1239  ASSERT(edg.side(s1)->valid(), "Invalid side of an edge.");
1240  SideIter s = edg.side(s1);
1241 
1242  // calculate the side diameter
1243  if (s->dim() == 0)
1244  {
1245  h = 1;
1246  }
1247  else
1248  {
1249  for (unsigned int i=0; i<s->n_nodes(); i++)
1250  for (unsigned int j=i+1; j<s->n_nodes(); j++)
1251  h = max(h, s->node(i)->distance(*s->node(j)));
1252  }
1253 
1254  // calculate the total in- and out-flux through the edge
1255  double pflux = 0, nflux = 0;
1256  for (int i=0; i<edg.n_sides; i++)
1257  {
1258  if (fluxes[i] > 0)
1259  pflux += fluxes[i];
1260  else
1261  nflux += fluxes[i];
1262  }
1263 
1264  // calculate the flux from s1 to s2
1265  if (fluxes[s2] > 0 && fluxes[s1] < 0 && s1 < s2)
1266  transport_flux = fluxes[s1]*fabs(fluxes[s2]/pflux);
1267  else if (fluxes[s2] < 0 && fluxes[s1] > 0 && s1 < s2)
1268  transport_flux = fluxes[s1]*fabs(fluxes[s2]/nflux);
1269  else if (s1==s2)
1270  transport_flux = fluxes[s1];
1271  else
1272  transport_flux = 0;
1273 
1274  gamma = 0.5*fabs(transport_flux);
1275 
1276 
1277  // determine local DG penalty
1278  local_alpha = max(alpha1, alpha2);
1279 
1280  if (s1 == s2)
1281  {
1282  omega[0] = 1;
1283 
1284  // delta is set to the average value of Kn.n on the side
1285  delta[0] = 0;
1286  for (unsigned int k=0; k<K_size; k++)
1287  delta[0] += dot(K1[k]*normal_vector,normal_vector);
1288  delta[0] /= K_size;
1289 
1290  gamma += local_alpha/h*(delta[0]);
1291  }
1292  else
1293  {
1294  delta[0] = 0;
1295  delta[1] = 0;
1296  for (unsigned int k=0; k<K_size; k++)
1297  {
1298  delta[0] += dot(K1[k]*normal_vector,normal_vector);
1299  delta[1] += dot(K2[k]*normal_vector,normal_vector);
1300  }
1301  delta[0] /= K_size;
1302  delta[1] /= K_size;
1303 
1304  double delta_sum = delta[0] + delta[1];
1305 
1306  if (delta_sum > numeric_limits<double>::epsilon())
1307  {
1308  omega[0] = delta[1]/delta_sum;
1309  omega[1] = delta[0]/delta_sum;
1310  gamma += local_alpha/h*(delta[0]*delta[1]/delta_sum);
1311  }
1312  else
1313  for (int i=0; i<2; i++) omega[i] = 0;
1314  }
1315 }
1316 
1317 
1318 
1319 
1320 
1321 
1322 template<class Model>
1323 void TransportDG<Model>::set_DG_parameters_boundary(const SideIter side,
1324  const int K_size,
1325  const vector<arma::mat33> &K,
1326  const double flux,
1327  const arma::vec3 &normal_vector,
1328  const double alpha,
1329  double &gamma)
1330 {
1331  double delta = 0, h = 0;
1332 
1333  // calculate the side diameter
1334  if (side->dim() == 0)
1335  {
1336  h = 1;
1337  }
1338  else
1339  {
1340  for (unsigned int i=0; i<side->n_nodes(); i++)
1341  for (unsigned int j=i+1; j<side->n_nodes(); j++)
1342  h = max(h, side->node(i)->distance( *side->node(j) ));
1343  }
1344 
1345  // delta is set to the average value of Kn.n on the side
1346  for (unsigned int k=0; k<K_size; k++)
1347  delta += dot(K[k]*normal_vector,normal_vector);
1348  delta /= K_size;
1349 
1350  gamma = 0.5*fabs(flux) + alpha/h*delta;
1351 }
1352 
1353 
1354 
1355 
1356 
1357 template<class Model>
1358 void TransportDG<Model>::set_initial_condition()
1359 {
1360  START_TIMER("set_init_cond");
1361  for (int sbi=0; sbi<n_subst_; sbi++)
1362  ls[sbi]->start_allocation();
1363  prepare_initial_condition<1>();
1364  prepare_initial_condition<2>();
1365  prepare_initial_condition<3>();
1366 
1367  for (int sbi=0; sbi<n_subst_; sbi++)
1368  ls[sbi]->start_add_assembly();
1369  prepare_initial_condition<1>();
1370  prepare_initial_condition<2>();
1371  prepare_initial_condition<3>();
1372 
1373  for (int sbi=0; sbi<n_subst_; sbi++)
1374  {
1375  ls[sbi]->finish_assembly();
1376  ls[sbi]->solve();
1377  }
1378  END_TIMER("set_init_cond");
1379 }
1380 
1381 template<class Model>
1382 template<unsigned int dim>
1383 void TransportDG<Model>::prepare_initial_condition()
1384 {
1385  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
1386  update_values | update_JxW_values | update_quadrature_points);
1387  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
1388  unsigned int dof_indices[ndofs];
1389  double matrix[ndofs*ndofs], rhs[ndofs];
1390  std::vector<arma::vec> init_values(qsize);
1391 
1392  for (int i=0; i<qsize; i++)
1393  init_values[i].resize(n_subst_);
1394 
1395  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
1396  {
1397  typename DOFHandlerBase::CellIterator elem = mesh_->element(feo->dh()->el_index(i_cell));
1398  if (elem->dim() != dim) continue;
1399 
1400  ElementAccessor<3> ele_acc = elem->element_accessor();
1401  feo->dh()->get_dof_indices(elem, dof_indices);
1402  fe_values.reinit(elem);
1403 
1404  Model::compute_init_cond(fe_values.point_list(), ele_acc, init_values);
1405 
1406  for (int sbi=0; sbi<n_subst_; sbi++)
1407  {
1408  for (unsigned int i=0; i<ndofs; i++)
1409  {
1410  rhs[i] = 0;
1411  for (unsigned int j=0; j<ndofs; j++)
1412  matrix[i*ndofs+j] = 0;
1413  }
1414 
1415  for (unsigned int k=0; k<qsize; k++)
1416  {
1417  double rhs_term = init_values[k](sbi)*fe_values.JxW(k);
1418 
1419  for (unsigned int i=0; i<ndofs; i++)
1420  {
1421  for (unsigned int j=0; j<ndofs; j++)
1422  matrix[i*ndofs+j] += fe_values.shape_value(i,k)*fe_values.shape_value(j,k)*fe_values.JxW(k);
1423 
1424  rhs[i] += fe_values.shape_value(i,k)*rhs_term;
1425  }
1426  }
1427  ls[sbi]->set_values(ndofs, (int *)dof_indices, ndofs, (int *)dof_indices, matrix, rhs);
1428  }
1429  }
1430 }
1431 
1432 
1433 template<class Model>
1434 void TransportDG<Model>::calc_fluxes(vector<vector<double> > &bcd_balance, vector<vector<double> > &bcd_plus_balance, vector<vector<double> > &bcd_minus_balance)
1435 {
1436  calc_fluxes<1>(bcd_balance, bcd_plus_balance, bcd_minus_balance);
1437  calc_fluxes<2>(bcd_balance, bcd_plus_balance, bcd_minus_balance);
1438  calc_fluxes<3>(bcd_balance, bcd_plus_balance, bcd_minus_balance);
1439 }
1440 
1441 template<class Model>
1442 template<unsigned int dim>
1443 void TransportDG<Model>::calc_fluxes(vector<vector<double> > &bcd_balance, vector<vector<double> > &bcd_plus_balance, vector<vector<double> > &bcd_minus_balance)
1444 {
1445  FESideValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe<dim>(),
1446  update_values | update_gradients | update_normal_vectors | update_side_JxW_values | update_quadrature_points);
1447  FESideValues<dim,3> fsv_rt(*feo->mapping<dim>(), *feo->q<dim-1>(), *feo->fe_rt<dim>(),
1448  update_values);
1449  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim-1>()->size();
1450  unsigned int dof_indices[ndofs];
1451  vector<arma::vec3> side_velocity(qsize);
1452  vector<arma::vec> bc_values(qsize);
1453  arma::vec3 conc_grad;
1454 
1455  for (int i=0; i<qsize; i++)
1456  bc_values[i].resize(n_subst_);
1457 
1458  for (int iedg=0; iedg<feo->dh()->n_loc_edges(); iedg++)
1459  {
1460  Edge *edg = &mesh_->edges[feo->dh()->edge_index(iedg)];
1461  if (edg->n_sides > 1) continue;
1462  if (edg->side(0)->dim() != dim-1) continue;
1463  // skip edges lying not on the boundary
1464  if (edg->side(0)->cond() == NULL) continue;
1465 
1466  SideIter side = edg->side(0);
1467  DOFHandlerMultiDim::CellIterator cell = side->element();
1468  ElementAccessor<3> ele_acc = cell->element_accessor();
1469  Region r = side->cond()->element_accessor().region();
1470  if (! r.is_valid()) xprintf(Msg, "Invalid region, ele % d, edg: % d\n", side->cond()->bc_ele_idx_, side->cond()->edge_idx_);
1471  unsigned int bc_region_idx = r.boundary_idx();
1472 
1473  fe_values.reinit(cell, side->el_idx());
1474  fsv_rt.reinit(cell, side->el_idx());
1475  feo->dh()->get_dof_indices(cell, dof_indices);
1476  calculate_velocity(cell, side_velocity, fsv_rt);
1477  Model::compute_advection_diffusion_coefficients(fe_values.point_list(), side_velocity, ele_acc, ad_coef, dif_coef);
1478  Model::compute_dirichlet_bc(fe_values.point_list(), side->cond()->element_accessor(), bc_values);
1479  arma::uvec bc_type = data_.bc_type.value(side->cond()->element()->centre(), side->cond()->element_accessor());
1480 
1481  for (int sbi=0; sbi<n_subst_; sbi++)
1482  {
1483  double mass_flux = 0;
1484  double water_flux = 0;
1485  for (unsigned int k=0; k<qsize; k++)
1486  water_flux += arma::dot(ad_coef[sbi][k], fe_values.normal_vector(k))*fe_values.JxW(k);
1487  water_flux /= side->measure();
1488 
1489  for (unsigned int k=0; k<qsize; k++)
1490  {
1491  double conc = 0;
1492  conc_grad.zeros();
1493  for (unsigned int i=0; i<ndofs; i++)
1494  {
1495  conc += fe_values.shape_value(i,k)*ls[sbi]->get_solution_array()[dof_indices[i]-feo->dh()->loffset()];
1496  conc_grad += fe_values.shape_grad(i,k)*ls[sbi]->get_solution_array()[dof_indices[i]-feo->dh()->loffset()];
1497  }
1498 
1499  // flux due to advection
1500  mass_flux += water_flux*conc*fe_values.JxW(k);
1501 
1502  if (bc_type[sbi] == EqData::dirichlet || (bc_type[sbi] == EqData::inflow && water_flux*side->measure() < -mh_dh->precision()))
1503  {
1504  // flux due to diffusion
1505  mass_flux -= arma::dot(dif_coef[sbi][k]*conc_grad,fe_values.normal_vector(k))*fe_values.JxW(k);
1506 
1507  // the penalty term has to be added otherwise the mass balance will not hold
1508  mass_flux -= gamma[sbi][side->cond_idx()]*(bc_values[k][sbi] - conc)*fe_values.JxW(k);
1509  }
1510 
1511  }
1512  bcd_balance[sbi][bc_region_idx] += mass_flux;
1513  if (mass_flux >= 0) bcd_plus_balance[sbi][bc_region_idx] += mass_flux;
1514  else bcd_minus_balance[sbi][bc_region_idx] += mass_flux;
1515  }
1516  }
1517 
1518 }
1519 
1520 template<class Model>
1521 void TransportDG<Model>::calc_elem_sources(vector<vector<double> > &mass, vector<vector<double> > &src_balance)
1522 {
1523  calc_elem_sources<1>(mass, src_balance);
1524  calc_elem_sources<2>(mass, src_balance);
1525  calc_elem_sources<3>(mass, src_balance);
1526 }
1527 
1528 template<class Model>
1529 template<unsigned int dim>
1530 void TransportDG<Model>::calc_elem_sources(vector<vector<double> > &mass, vector<vector<double> > &src_balance)
1531 {
1532  FEValues<dim,3> fe_values(*feo->mapping<dim>(), *feo->q<dim>(), *feo->fe<dim>(),
1533  update_values | update_JxW_values | update_quadrature_points);
1534  const unsigned int ndofs = feo->fe<dim>()->n_dofs(), qsize = feo->q<dim>()->size();
1535  unsigned int dof_indices[ndofs];
1536  vector<arma::vec> sources_conc(qsize), sources_density(qsize), sources_sigma(qsize);
1537  double mass_sum, sources_sum, conc, conc_diff;
1538 
1539  for (int i_cell=0; i_cell<feo->dh()->el_ds()->lsize(); i_cell++)
1540  {
1541  typename DOFHandlerBase::CellIterator elem = mesh_->element(feo->dh()->el_index(i_cell));
1542  if (elem->dim() != dim) continue;
1543 
1544  Region r = elem->element_accessor().region();
1545  if (! r.is_valid()) xprintf(Msg, "Invalid region, ele % d\n", elem.index());
1546  unsigned int region_idx = r.bulk_idx();
1547 
1548  ElementAccessor<3> ele_acc = elem->element_accessor();
1549  fe_values.reinit(elem);
1550  feo->dh()->get_dof_indices(elem, dof_indices);
1551 
1552  Model::compute_mass_matrix_coefficient(fe_values.point_list(), ele_acc, mm_coef);
1553  Model::compute_source_coefficients(fe_values.point_list(), ele_acc, sources_conc, sources_density, sources_sigma);
1554 
1555  for (int sbi=0; sbi<n_subst_; sbi++)
1556  {
1557  mass_sum = 0;
1558  sources_sum = 0;
1559 
1560  for (unsigned int k=0; k<qsize; k++)
1561  {
1562  conc = 0;
1563  for (unsigned int i=0; i<ndofs; i++)
1564  conc += fe_values.shape_value(i,k)*ls[sbi]->get_solution_array()[dof_indices[i]-feo->dh()->loffset()];
1565 
1566  mass_sum += mm_coef[k]*conc*fe_values.JxW(k);
1567 
1568  conc_diff = sources_conc[k][sbi] - conc;
1569  sources_sum += (sources_density[k][sbi] + conc_diff*sources_sigma[k][sbi])*fe_values.JxW(k);
1570  }
1571 
1572  mass[sbi][region_idx] += mass_sum;
1573  src_balance[sbi][region_idx] += sources_sum;
1574  }
1575  }
1576 }
1577 
1578 
1579 
1580 
1581 
1582 template class TransportDG<ConcentrationTransportModel>;
1583 template class TransportDG<HeatTransferModel>;
1584 
1585 
1586 
1587 
Generated on Thu May 29 2014 23:14:49 for Flow123d by   doxygen 1.8.4