DUNE PDELab (git)

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