dune-pdelab  2.4.1
convectiondiffusionccfv.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH
4 
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/typetraits.hh>
8 #include<dune/geometry/referenceelements.hh>
9 
18 
19 namespace Dune {
20  namespace PDELab {
21 
37  template<typename TP>
39  :
40  // public NumericalJacobianSkeleton<ConvectionDiffusionCCFV<TP> >,
41  // public NumericalJacobianBoundary<ConvectionDiffusionCCFV<TP> >,
42  // public NumericalJacobianVolume<ConvectionDiffusionCCFV<TP> >,
43  public FullSkeletonPattern,
44  public FullVolumePattern,
46  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
47  {
48  using BCType = typename ConvectionDiffusionBoundaryConditions::Type;
49 
50  public:
51  // pattern assembly flags
52  enum { doPatternVolume = true };
53  enum { doPatternSkeleton = true };
54 
55  // residual assembly flags
56  enum { doAlphaVolume = true };
57  enum { doAlphaSkeleton = true };
58  enum { doAlphaBoundary = true };
59  enum { doLambdaVolume = true };
60  enum { doLambdaSkeleton = false };
61  enum { doLambdaBoundary = false };
62 
63  ConvectionDiffusionCCFV (TP& param_) : param(param_)
64  {}
65 
66  // volume integral depending on test and ansatz functions
67  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
68  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
69  {
70  // cell center
71  auto geo = eg.geometry();
72  auto ref_el = referenceElement(geo);
73  auto local_inside = ref_el.position(0,0);
74 
75  // evaluate reaction term
76  auto c = param.c(eg.entity(),local_inside);
77 
78  // and accumulate
79  r.accumulate(lfsu,0,(c*x(lfsu,0))*geo.volume());
80  }
81 
82  // jacobian of volume term
83  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
84  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
85  M& mat) const
86  {
87  // cell center
88  auto geo = eg.geometry();
89  auto ref_el = referenceElement(geo);
90  auto local_inside = ref_el.position(0,0);
91 
92  // evaluate reaction term
93  auto c = param.c(eg.entity(),local_inside);
94 
95  // and accumulate
96  mat.accumulate(lfsu,0,lfsu,0,c*geo.volume());
97  }
98 
99  // skeleton integral depending on test and ansatz functions
100  // each face is only visited ONCE!
101  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
102  void alpha_skeleton (const IG& ig,
103  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
104  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
105  R& r_s, R& r_n) const
106  {
107  // define types
108  using RF = typename LFSU::Traits::FiniteElementType::
109  Traits::LocalBasisType::Traits::RangeFieldType;
110 
111  // dimensions
112  const auto dim = IG::dimension;
113 
114  // get cell entities from both sides of the intersection
115  auto cell_inside = ig.inside();
116  auto cell_outside = ig.outside();
117 
118  // get geometries
119  auto geo = ig.geometry();
120  auto geo_inside = cell_inside.geometry();
121  auto geo_outside = cell_outside.geometry();
122 
123  // get geometry of intersection in local coordinates of neighbor cells
124  auto geo_in_inside = ig.geometryInInside();
125 
126  // center in face's reference element
127  auto ref_el = referenceElement(geo);
128  auto face_local = ref_el.position(0,0);
129 
130  // face volume for integration
131  auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
132 
133  // cell centers in references elements
134  auto ref_el_inside = referenceElement(geo_inside);
135  auto ref_el_outside = referenceElement(geo_outside);
136  auto local_inside = ref_el_inside.position(0,0);
137  auto local_outside = ref_el_outside.position(0,0);
138 
139  // evaluate diffusion coefficient from either side and take harmonic average
140  auto tensor_inside = param.A(cell_inside,local_inside);
141  auto tensor_outside = param.A(cell_outside,local_outside);
142  auto n_F = ig.centerUnitOuterNormal();
143  Dune::FieldVector<RF,dim> An_F;
144  tensor_inside.mv(n_F,An_F);
145  auto k_inside = n_F*An_F;
146  tensor_outside.mv(n_F,An_F);
147  auto k_outside = n_F*An_F;
148  auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
149 
150  // evaluate convective term
151  auto iplocal_s = geo_in_inside.global(face_local);
152  auto b = param.b(cell_inside,iplocal_s);
153  auto vn = b*n_F;
154  auto u_upwind=0;
155  if (vn>=0) u_upwind = x_s(lfsu_s,0); else u_upwind = x_n(lfsu_n,0);
156 
157  // cell centers in global coordinates
158  auto global_inside = geo_inside.global(local_inside);
159  auto global_outside = geo_outside.global(local_outside);
160 
161  // distance between the two cell centers
162  global_inside -= global_outside;
163  auto distance = global_inside.two_norm();
164 
165  // contribution to residual on inside element, other residual is computed by symmetric call
166  r_s.accumulate(lfsu_s,0,(u_upwind*vn)*face_volume+k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
167  r_n.accumulate(lfsu_n,0,-(u_upwind*vn)*face_volume-k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
168  }
169 
170  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
171  void jacobian_skeleton (const IG& ig,
172  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
173  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
174  M& mat_ss, M& mat_sn,
175  M& mat_ns, M& mat_nn) const
176  {
177  // define types
178  using RF = typename LFSU::Traits::FiniteElementType::
179  Traits::LocalBasisType::Traits::RangeFieldType;
180 
181  // dimensions
182  const auto dim = IG::dimension;
183 
184  // get cell entities from both sides of the intersection
185  auto cell_inside = ig.inside();
186  auto cell_outside = ig.outside();
187 
188  // get geometries
189  auto geo = ig.geometry();
190  auto geo_inside = cell_inside.geometry();
191  auto geo_outside = cell_outside.geometry();
192 
193  // get geometry of intersection in local coordinates of neighbor cells
194  auto geo_in_inside = ig.geometryInInside();
195 
196  // center in face's reference element
197  auto ref_el = referenceElement(geo);
198  auto face_local = ref_el.position(0,0);
199 
200  // face volume for integration
201  auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
202 
203  // cell centers in references elements
204  auto ref_el_inside = referenceElement(geo_inside);
205  auto ref_el_outside = referenceElement(geo_outside);
206  auto local_inside = ref_el_inside.position(0,0);
207  auto local_outside = ref_el_outside.position(0,0);
208 
209  // evaluate diffusion coefficient from either side and take harmonic average
210  auto tensor_inside = param.A(cell_inside,local_inside);
211  auto tensor_outside = param.A(cell_outside,local_outside);
212  auto n_F = ig.centerUnitOuterNormal();
213  Dune::FieldVector<RF,dim> An_F;
214  tensor_inside.mv(n_F,An_F);
215  auto k_inside = n_F*An_F;
216  tensor_outside.mv(n_F,An_F);
217  auto k_outside = n_F*An_F;
218  auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
219 
220  // evaluate convective term
221  auto iplocal_s = geo_in_inside.global(face_local);
222  auto b = param.b(cell_inside,iplocal_s);
223  auto vn = b*n_F;
224 
225  // cell centers in global coordinates
226  auto global_inside = geo_inside.global(local_inside);
227  auto global_outside = geo_outside.global(local_outside);
228 
229  // distance between the two cell centers
230  global_inside -= global_outside;
231  auto distance = global_inside.two_norm();
232 
233  // contribution to residual on inside element, other residual is computed by symmetric call
234  mat_ss.accumulate(lfsu_s,0,lfsu_s,0, k_avg*face_volume/distance );
235  mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -k_avg*face_volume/distance );
236  mat_sn.accumulate(lfsu_s,0,lfsu_n,0, -k_avg*face_volume/distance );
237  mat_nn.accumulate(lfsu_n,0,lfsu_n,0, k_avg*face_volume/distance );
238  if (vn>=0)
239  {
240  mat_ss.accumulate(lfsu_s,0,lfsu_s,0, vn*face_volume );
241  mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -vn*face_volume );
242  }
243  else
244  {
245  mat_sn.accumulate(lfsu_s,0,lfsu_n,0, vn*face_volume );
246  mat_nn.accumulate(lfsu_n,0,lfsu_n,0, -vn*face_volume );
247  }
248  }
249 
250 
251  // post skeleton: compute time step allowable for cell; to be done later
252  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
253  void alpha_volume_post_skeleton(const EG& eg, const LFSU& lfsu, const X& x,
254  const LFSV& lfsv, R& r) const
255  {
256  if (!first_stage) return; // time step calculation is only done in first stage
257 
258  // cell center
259  auto geo = eg.geometry();
260  auto ref_el = referenceElement(geo);
261  auto local_inside = ref_el.position(0,0);
262 
263  // compute optimal dt for this cell
264  auto cellcapacity = param.d(eg.entity(),local_inside)*geo.volume();
265  auto celldt = cellcapacity/(cellinflux+1E-30);
266  dtmin = std::min(dtmin,celldt);
267  }
268 
269 
270  // skeleton integral depending on test and ansatz functions
271  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
272  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
273  void alpha_boundary (const IG& ig,
274  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
275  R& r_s) const
276  {
277  // define types
278  using RF = typename LFSU::Traits::FiniteElementType::
279  Traits::LocalBasisType::Traits::RangeFieldType;
280 
281  // dimensions
282  const auto dim = IG::dimension;
283 
284  // get cell entities from both sides of the intersection
285  auto cell_inside = ig.inside();
286 
287  // get geometries
288  auto geo = ig.geometry();
289  auto geo_inside = cell_inside.geometry();
290 
291  // get geometry of intersection in local coordinates of neighbor cells
292  auto geo_in_inside = ig.geometryInInside();
293 
294  // center in face's reference element
295  auto ref_el = referenceElement(geo);
296  auto face_local = ref_el.position(0,0);
297 
298  // face volume for integration
299  auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
300 
301  // cell centers in references elements
302  auto ref_el_inside = referenceElement(geo_inside);
303  auto local_inside = ref_el_inside.position(0,0);
304 
305  // evaluate boundary condition type
306  auto bctype = param.bctype(ig.intersection(),face_local);
307 
309  {
310  // Dirichlet boundary
311  // distance between cell center and face center
312  auto global_inside = geo_inside.global(local_inside);
313  auto global_outside = geo.global(face_local);
314  global_inside -= global_outside;
315  auto distance = global_inside.two_norm();
316 
317  // evaluate diffusion coefficient
318  auto tensor_inside = param.A(cell_inside,local_inside);
319  auto n_F = ig.centerUnitOuterNormal();
320  Dune::FieldVector<RF,dim> An_F;
321  tensor_inside.mv(n_F,An_F);
322  auto k_inside = n_F*An_F;
323 
324  // evaluate boundary condition function
325  auto iplocal_s = geo_in_inside.global(face_local);
326  auto g = param.g(cell_inside,iplocal_s);
327 
328  // velocity needed for convection term
329  auto b = param.b(cell_inside,iplocal_s);
330  auto n = ig.centerUnitOuterNormal();
331 
332  // contribution to residual on inside element, assumes that Dirichlet boundary is inflow
333  r_s.accumulate(lfsu_s,0,(b*n)*g*face_volume + k_inside*(x_s(lfsu_s,0)-g)*face_volume/distance);
334 
335  return;
336  }
337 
339  {
340  // Neumann boundary
341  // evaluate flux boundary condition
342 
343  //evaluate boundary function
344  auto j = param.j(ig.intersection(),face_local);
345 
346  // contribution to residual on inside element
347  r_s.accumulate(lfsu_s,0,j*face_volume);
348 
349  return;
350  }
351 
353  {
354  // evaluate velocity field and outer unit normal
355  auto iplocal_s = geo_in_inside.global(face_local);
356  auto b = param.b(cell_inside,iplocal_s);
357  auto n = ig.centerUnitOuterNormal();
358 
359  // evaluate outflow boundary condition
360  auto o = param.o(ig.intersection(),face_local);
361 
362  // integrate o
363  r_s.accumulate(lfsu_s,0,((b*n)*x_s(lfsu_s,0) + o)*face_volume);
364 
365  return;
366  }
367  }
368 
369  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
370  void jacobian_boundary (const IG& ig,
371  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
372  M& mat_ss) const
373  {
374  // define types
375  using RF = typename LFSU::Traits::FiniteElementType::
376  Traits::LocalBasisType::Traits::RangeFieldType;
377 
378  // dimensions
379  const auto dim = IG::dimension;
380 
381  // get cell entities from both sides of the intersection
382  auto cell_inside = ig.inside();
383 
384  // get geometries
385  auto geo = ig.geometry();
386  auto geo_inside = cell_inside.geometry();
387 
388  // get geometry of intersection in local coordinates of neighbor cells
389  auto geo_in_inside = ig.geometryInInside();
390 
391  // center in face's reference element
392  auto ref_el = referenceElement(geo);
393  auto face_local = ref_el.position(0,0);
394 
395  // face volume for integration
396  auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
397 
398  // cell centers in references elements
399  auto ref_el_inside = referenceElement(geo_inside);
400  auto local_inside = ref_el_inside.position(0,0);
401 
402  // evaluate boundary condition type
403  auto bctype = param.bctype(ig.intersection(),face_local);
404 
406  {
407  // Dirichlet boundary
408  // distance between cell center and face center
409  auto global_inside = geo_inside.global(local_inside);
410  auto global_outside = geo.global(face_local);
411  global_inside -= global_outside;
412  auto distance = global_inside.two_norm();
413 
414  // evaluate diffusion coefficient
415  auto tensor_inside = param.A(cell_inside,local_inside);
416  auto n_F = ig.centerUnitOuterNormal();
417  Dune::FieldVector<RF,dim> An_F;
418  tensor_inside.mv(n_F,An_F);
419  auto k_inside = n_F*An_F;
420 
421  // contribution to residual on inside element
422  mat_ss.accumulate(lfsu_s,0,lfsv_s,0, k_inside*face_volume/distance );
423 
424  return;
425  }
426 
428  {
429  // evaluate velocity field and outer unit normal
430  auto iplocal_s = geo_in_inside.global(face_local);
431  auto b = param.b(cell_inside,iplocal_s);
432  auto n = ig.centerUnitOuterNormal();
433 
434  // integrate o
435  mat_ss.accumulate(lfsu_s,0,lfsv_s,0, (b*n)*face_volume );
436 
437  return;
438  }
439  }
440 
441  // volume integral depending only on test functions
442  template<typename EG, typename LFSV, typename R>
443  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
444  {
445  // cell center
446  auto geo = eg.geometry();
447  auto ref_el = referenceElement(geo);
448  auto local_inside = ref_el.position(0,0);
449 
450  // evaluate source and sink term
451  auto f = param.f(eg.entity(),local_inside);
452 
453  r.accumulate(lfsv,0,-f*geo.volume());
454  }
455 
457  void setTime (typename TP::Traits::RangeFieldType t)
458  {
459  param.setTime(t);
460  }
461 
463  void preStep (typename TP::Traits::RangeFieldType time, typename TP::Traits::RangeFieldType dt,
464  int stages)
465  {
466  }
467 
469  void preStage (typename TP::Traits::RangeFieldType time, int r)
470  {
471  if (r==1)
472  {
473  first_stage = true;
474  dtmin = 1E100;
475  }
476  else first_stage = false;
477  }
478 
480  void postStage ()
481  {
482  }
483 
485  typename TP::Traits::RangeFieldType suggestTimestep (typename TP::Traits::RangeFieldType dt) const
486  {
487  return dtmin;
488  }
489 
490  private:
491  TP& param;
492  bool first_stage;
493  mutable typename TP::Traits::RangeFieldType dtmin; // accumulate minimum dt here
494  mutable typename TP::Traits::RangeFieldType cellinflux;
495  };
496 
497 
498 
499 
506  template<class TP>
508  :
509  // public NumericalJacobianApplyVolume<ConvectionDiffusionCCFVTemporalOperator<TP> >,
510  public FullVolumePattern,
512  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
513  {
514  public:
515  // pattern assembly flags
516  enum { doPatternVolume = true };
517 
518  // residual assembly flags
519  enum { doAlphaVolume = true };
520 
522  : param(param_)
523  {
524  }
525 
526  // volume integral depending on test and ansatz functions
527  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
528  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
529  {
530  // cell center
531  auto geo = eg.geometry();
532  auto ref_el = referenceElement(geo);
533  auto local_inside = ref_el.position(0,0);
534 
535  // capacity term
536  auto capacity = param.d(eg.entity(),local_inside);
537 
538  // residual contribution
539  r.accumulate(lfsu,0,capacity*x(lfsu,0)*geo.volume());
540  }
541 
542  // jacobian of volume term
543  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
544  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
545  M& mat) const
546  {
547  // cell center
548  auto geo = eg.geometry();
549  auto ref_el = referenceElement(geo);
550  auto local_inside = ref_el.position(0,0);
551 
552  // capacity term
553  auto capacity = param.d(eg.entity(),local_inside);
554 
555  // residual contribution
556  mat.accumulate(lfsu,0,lfsu,0,capacity*geo.volume());
557  }
558 
559  private:
560  TP& param;
561  };
562 
563 
565  } // namespace PDELab
566 } // namespace Dune
567 
568 #endif
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:68
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:528
Definition: convectiondiffusionccfv.hh:52
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionccfv.hh:273
Definition: convectiondiffusionccfv.hh:59
sparsity pattern generator
Definition: pattern.hh:13
Definition: convectiondiffusionccfv.hh:58
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:253
const IG & ig
Definition: constraints.hh:148
void preStage(typename TP::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: convectiondiffusionccfv.hh:469
static const int dim
Definition: adaptivity.hh:83
Definition: convectiondiffusionccfv.hh:60
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Definition: adaptivity.hh:27
Definition: convectiondiffusionparameter.hh:65
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: convectiondiffusionccfv.hh:457
Type
Definition: convectiondiffusionparameter.hh:65
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:84
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: convectiondiffusionccfv.hh:102
sparsity pattern generator
Definition: pattern.hh:29
void preStep(typename TP::Traits::RangeFieldType time, typename TP::Traits::RangeFieldType dt, int stages)
to be called once before each time step
Definition: convectiondiffusionccfv.hh:463
Default flags for all local operators.
Definition: flags.hh:18
ConvectionDiffusionCCFVTemporalOperator(TP &param_)
Definition: convectiondiffusionccfv.hh:521
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:443
Definition: convectiondiffusionparameter.hh:65
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, M &mat_ss, M &mat_sn, M &mat_ns, M &mat_nn) const
Definition: convectiondiffusionccfv.hh:171
Definition: convectiondiffusionccfv.hh:61
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: convectiondiffusionccfv.hh:507
Definition: convectiondiffusionparameter.hh:65
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: convectiondiffusionccfv.hh:485
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:544
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_ss) const
Definition: convectiondiffusionccfv.hh:370
ConvectionDiffusionCCFV(TP &param_)
Definition: convectiondiffusionccfv.hh:63
Definition: convectiondiffusionccfv.hh:38
Definition: convectiondiffusionccfv.hh:56
Definition: convectiondiffusionccfv.hh:53
Definition: convectiondiffusionccfv.hh:57
void postStage()
to be called once at the end of each stage
Definition: convectiondiffusionccfv.hh:480