DUNE PDELab (git)

twophaseccfv.hh
1// -*- tab-width: 4; indent-tabs-mode: nil -*-
2#ifndef DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
3#define DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
4
5#include<dune/common/exceptions.hh>
6#include<dune/common/fvector.hh>
7#include<dune/geometry/referenceelements.hh>
8#include<dune/localfunctions/raviartthomas/raviartthomascube.hh>
9
10#include<dune/pdelab/common/quadraturerules.hh>
11#include <dune/pdelab/localoperator/defaultimp.hh>
12
13#include"../common/function.hh"
14#include"pattern.hh"
15#include"flags.hh"
16#include"idefault.hh"
17
18namespace Dune {
19 namespace PDELab {
20
22 template<typename GV, typename RF>
23 struct TwoPhaseParameterTraits
24 {
26 using GridViewType = GV;
27
29 enum {
31 dimDomain = GV::dimension
32 };
33
35 using DomainFieldType = typename GV::Grid::ctype;
36
38 using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>;
39
41 using IntersectionDomainType = Dune::FieldVector<DomainFieldType,dimDomain-1>;
42
44 using RangeFieldType = RF;
45
47 using RangeType = Dune::FieldVector<RF,GV::dimensionworld>;
48
50 using PermTensorType = RangeFieldType;
51
53 using ElementType = typename GV::Traits::template Codim<0>::Entity;
54 using IntersectionType = typename GV::Intersection;
55 };
56
57 template<typename GV, typename RF>
58 struct TwoPhaseFullTensorParameterTraits : TwoPhaseParameterTraits<GV, RF>
59 {
60 using Base = TwoPhaseParameterTraits<GV, RF>;
61 using RangeFieldType = typename Base::RangeFieldType;
62
64 using PermTensorType = Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>;
65 };
66
68 template<class T, class Imp>
69 class TwoPhaseParameterInterface
70 {
71 public:
72 using Traits = T;
73
75 typename Traits::RangeFieldType
76 phi (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
77 {
78 return asImp().phi(e,x);
79 }
80
82 typename Traits::RangeFieldType
83 pc (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
84 typename Traits::RangeFieldType s_l) const
85 {
86 return asImp().pc(e,x,s_l);
87 }
88
90 typename Traits::RangeFieldType
91 s_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
92 typename Traits::RangeFieldType pc) const
93 {
94 return asImp().s_l(e,x,pc);
95 }
96
98 typename Traits::RangeFieldType
99 kr_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
100 typename Traits::RangeFieldType s_l) const
101 {
102 return asImp().kr_l(e,x,s_l);
103 }
104
106 typename Traits::RangeFieldType
107 kr_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
108 typename Traits::RangeFieldType s_g) const
109 {
110 return asImp().kr_g(e,x,s_g);
111 }
112
114 typename Traits::RangeFieldType
115 mu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
116 typename Traits::RangeFieldType p_l) const
117 {
118 return asImp().mu_l(e,x,p_l);
119 }
120
122 typename Traits::RangeFieldType
123 mu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
124 typename Traits::RangeFieldType p_g) const
125 {
126 return asImp().mu_l(e,x,p_g);
127 }
128
130 typename Traits::PermTensorType
131 k_abs (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
132 {
133 return asImp().k_abs(e,x);
134 }
135
137 const typename Traits::RangeType& gravity () const
138 {
139 return asImp().gravity();
140 }
141
143 template<typename E>
144 typename Traits::RangeFieldType
145 nu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
146 typename Traits::RangeFieldType p_l) const
147 {
148 return asImp().nu_l(e,x,p_l);
149 }
150
152 typename Traits::RangeFieldType
153 nu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
154 typename Traits::RangeFieldType p_g) const
155 {
156 return asImp().nu_g(e,x,p_g);
157 }
158
160 typename Traits::RangeFieldType
161 rho_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
162 typename Traits::RangeFieldType p_l) const
163 {
164 return asImp().rho_l(e,x,p_l);
165 }
166
168 typename Traits::RangeFieldType
169 rho_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
170 typename Traits::RangeFieldType p_g) const
171 {
172 return asImp().rho_g(e,x,p_g);
173 }
174
176 int
177 bc_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
178 {
179 return asImp().bc_l(is,x,time);
180 }
181
183 int
184 bc_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
185 {
186 return asImp().bc_g(is,x,time);
187 }
188
190 typename Traits::RangeFieldType
191 g_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
192 {
193 return asImp().g_l(is,x,time);
194 }
195
197 typename Traits::RangeFieldType
198 g_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
199 {
200 return asImp().g_g(is,x,time);
201 }
202
204 typename Traits::RangeFieldType
205 j_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
206 {
207 return asImp().j_l(is,x,time);
208 }
209
211 typename Traits::RangeFieldType
212 j_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
213 {
214 return asImp().j_g(is,x,time);
215 }
216
218 typename Traits::RangeFieldType
219 q_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
220 typename Traits::RangeFieldType time) const
221 {
222 return asImp().q_l(e,x,time);
223 }
224
226 typename Traits::RangeFieldType
227 q_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
228 typename Traits::RangeFieldType time) const
229 {
230 return asImp().q_g(e,x,time);
231 }
232
233 private:
234 Imp& asImp () {return static_cast<Imp &> (*this);}
235 const Imp& asImp () const {return static_cast<const Imp &>(*this);}
236 };
237
238
239 // a local operator for solving the two-phase flow in pressure-pressure formulation
240 // with two-point flux approximation
241 // TP : parameter class, see above
242 // V : Vector holding last time step
243 template<typename TP>
244 class TwoPhaseTwoPointFluxOperator
245 : public NumericalJacobianSkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
246 public NumericalJacobianApplySkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
247
248 public NumericalJacobianBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
249 public NumericalJacobianApplyBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
250
251 public FullSkeletonPattern,
252 public FullVolumePattern,
253 public LocalOperatorDefaultFlags,
254
255 public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
256 {
257 enum { dim = TP::Traits::GridViewType::dimension };
258 enum { liquid = 0 };
259 enum { gas = 1 };
260
261 using Real = typename TP::Traits::RangeFieldType;
262 public:
263 // pattern assembly flags
264 enum { doPatternVolume = true };
265 enum { doPatternSkeleton = true };
266
267 // residual assembly flags
268 enum { doAlphaSkeleton = true };
269 enum { doAlphaBoundary = true };
270 enum { doLambdaVolume = true };
271 enum { doLambdaBoundary = true };
272
274 TwoPhaseTwoPointFluxOperator (const TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
275 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
276 {}
277
278 // volume integral depending only on test functions
279 template<typename EG, typename LFSV, typename R>
280 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
281 {
282 // Reference to cell
283 const auto& cell = eg.entity();
284
285 // get geometry
286 auto geo = eg.geometry();
287
288 // cell geometry
289 auto ref_el = referenceElement(geo);
290 auto cell_center_local = ref_el.position(0,0);
291 auto cell_volume = geo.volume();
292
293 // contribution from source term
294 r.accumulate(lfsv, liquid, -scale_l * tp.q_l(cell,cell_center_local,time) * cell_volume);
295 r.accumulate(lfsv, gas, -scale_g * tp.q_g(cell,cell_center_local,time) * cell_volume);
296 }
297
298 // skeleton integral depending on test and ansatz functions
299 // each face is only visited ONCE!
300 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
301 void alpha_skeleton (const IG& ig,
302 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
303 const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
304 R& r_s, R& r_n) const
305 {
306 // Define types
307 using namespace Indices;
308 using PLSpace = TypeTree::Child<LFSV,_0>;
309 using RF = typename PLSpace::Traits::FiniteElementType::
310 Traits::LocalBasisType::Traits::RangeFieldType;
311
312 // References to inside and outside cells
313 const auto& cell_inside = ig.inside();
314 const auto& cell_outside = ig.outside();
315
316 // get geometries
317 auto geo = ig.geometry();
318 auto geo_inside = cell_inside.geometry();
319 auto geo_outside = cell_outside.geometry();
320
321 // cell geometries
322 auto ref_el_inside = referenceElement(geo_inside);
323 auto ref_el_outside = referenceElement(geo_outside);
324 auto inside_cell_center_local = ref_el_inside.position(0,0);
325 auto outside_cell_center_local = ref_el_outside.position(0,0);
326 auto inside_cell_center_global = geo_inside.center();
327 auto outside_cell_center_global = geo_outside.center();
328
329 // distance of cell centers
330 auto d = outside_cell_center_global;
331 d -= inside_cell_center_global;
332 auto distance = d.two_norm();
333
334 // face geometry
335 auto ref_el = referenceElement(geo);
336 auto face_local = ref_el.position(0,0);
337 auto face_volume = geo.volume();
338
339 // absolute permeability
340 auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
341 auto k_abs_outside = tp.k_abs(cell_outside,outside_cell_center_local);
342
343 // gravity times normal
344 auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
345
346 // liquid phase calculation
347 auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
348 auto rho_l_outside = tp.rho_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
349 auto w_l = (x_s(lfsu_s,liquid)-x_n(lfsu_n,liquid))/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction
350 RF pc_upwind, s_l_upwind, s_g_upwind;
351 auto nu_l = aavg(tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)),
352 tp.nu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid)));
353 if (w_l>=0) // upwind capillary pressure on face
354 {
355 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
356 s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
357 }
358 else
359 {
360 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
361 s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
362 }
363 s_g_upwind = 1-s_l_upwind;
364 auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l_upwind)/
365 tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
366 auto lambda_l_outside = tp.kr_l(cell_outside,outside_cell_center_local,s_l_upwind)/
367 tp.mu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
368 auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
369
370 r_s.accumulate(lfsv_s,liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
371 r_n.accumulate(lfsv_n,liquid, -scale_l * nu_l * sigma_l * w_l * face_volume);
372
373 // gas phase calculation
374 auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
375 auto rho_g_outside = tp.rho_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
376 auto w_g = (x_s(lfsu_s,gas)-x_n(lfsu_n,gas))/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction
377 auto nu_g = aavg(tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)),
378 tp.nu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas)));
379 if (w_l*w_g<=0) // new evaluation necessary only if signs differ
380 {
381 if (w_g>=0) // upwind capillary pressure on face
382 {
383 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
384 s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
385 }
386 else
387 {
388 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
389 s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
390 }
391 s_g_upwind = 1-s_l_upwind;
392 }
393 auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g_upwind)/
394 tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
395 auto lambda_g_outside = tp.kr_g(cell_outside,outside_cell_center_local,s_g_upwind)/
396 tp.mu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
397 auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
398
399 r_s.accumulate(lfsv_s, gas, scale_g * nu_g * sigma_g * w_g * face_volume);
400 r_n.accumulate(lfsv_n, gas, -scale_g * nu_g * sigma_g * w_g * face_volume);
401 }
402
403 // skeleton integral depending on test and ansatz functions
404 // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
405 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
406 void alpha_boundary (const IG& ig,
407 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
408 R& r_s) const
409 {
410 // References to inside cell
411 const auto& cell_inside = ig.inside();
412
413 // get geometries
414 auto geo = ig.geometry();
415 auto geo_inside = cell_inside.geometry();
416
417 // face geometry
418 auto ref_el = referenceElement(geo);
419 auto face_local = ref_el.position(0,0);
420 auto face_volume = geo.volume();
421
422 // evaluate boundary condition type
423 auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
424 auto bc_g = tp.bc_g(ig.intersection(),face_local,time);
425 if (bc_l!=1 && bc_g!=1) return; // no Dirichlet boundary conditions
426
427 // cell geometry
428 auto ref_el_inside = referenceElement(geo_inside);
429 auto inside_cell_center_local = ref_el_inside.position(0,0);
430 auto inside_cell_center_global = geo_inside.center();
431
432 // distance of cell center to boundary
433 auto d = geo.global(face_local);
434 d -= inside_cell_center_global;
435 auto distance = d.two_norm();
436
437 // absolute permeability
438 auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
439
440 // gravity times normal
441 auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
442
443 // liquid phase Dirichlet boundary
444 if (bc_l==1)
445 {
446 auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
447 auto g_l = tp.g_l(ig.intersection(),face_local,time);
448 auto w_l = (x_s(lfsu_s,liquid)-g_l)/distance + rho_l_inside*gn;
449 auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
450 auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l)/
451 tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
452 auto sigma_l = lambda_l_inside*k_abs_inside;
453 auto nu_l = tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
454 r_s.accumulate(lfsv_s, liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
455 }
456
457 // gas phase Dirichlet boundary
458 if (bc_g==1)
459 {
460 auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
461 auto g_g = tp.g_g(ig.intersection(),face_local,time);
462 auto w_g = (x_s(lfsu_s,gas)-g_g)/distance + rho_g_inside*gn;
463 auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
464 auto s_g = 1-s_l;
465 auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g)/
466 tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
467 auto sigma_g = lambda_g_inside*k_abs_inside;
468 auto nu_g = tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
469 r_s.accumulate(lfsv_s, gas, scale_l * nu_g * sigma_g * w_g * face_volume);
470 }
471 }
472
473 // boundary integral independent of ansatz functions
474 template<typename IG, typename LFSV, typename R>
475 void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r_s) const
476 {
477 // get geometries
478 auto geo = ig.geometry();
479
480 // face geometry
481 auto ref_el = referenceElement(geo);
482 auto face_local = ref_el.position(0,0);
483 auto face_volume = geo.volume();
484
485 // evaluate boundary condition type
486 auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
487 auto bc_g = tp.bc_g(ig.intersection(),face_local,time);
488 if (bc_l!=0 && bc_g!=0) return; // no Neumann boundary conditions
489
490 // liquid phase Neumann boundary
491 if (bc_l==0)
492 {
493 auto j_l = tp.j_l(ig.intersection(),face_local,time);
494 r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
495 }
496
497 // gas phase Neumann boundary
498 if (bc_g==0)
499 {
500 auto j_g = tp.j_g(ig.intersection(),face_local,time);
501 r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
502 }
503 }
504
506 void setTime (typename TP::Traits::RangeFieldType t)
507 {
508 time = t;
509 }
510
511 private:
512 const TP& tp; // two phase parameter class
513 typename TP::Traits::RangeFieldType time;
514 Real scale_l, scale_g;
515
516 template<typename T>
517 T aavg (T a, T b) const
518 {
519 return 0.5*(a+b);
520 }
521
522 template<typename T>
523 T havg (T a, T b) const
524 {
525 T eps = 1E-30;
526 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
527 }
528 };
529
530
537 template<class TP>
538 class TwoPhaseOnePointTemporalOperator
539 : public NumericalJacobianVolume<TwoPhaseOnePointTemporalOperator<TP> >,
540 public NumericalJacobianApplyVolume<TwoPhaseOnePointTemporalOperator<TP> >,
541 public FullVolumePattern,
542 public LocalOperatorDefaultFlags,
543 public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
544 {
545 enum { dim = TP::Traits::GridViewType::dimension };
546 enum { liquid = 0 };
547 enum { gas = 1 };
548
549 using Real = typename TP::Traits::RangeFieldType;
550
551 public:
552 // pattern assembly flags
553 enum { doPatternVolume = true };
554
555 // residual assembly flags
556 enum { doAlphaVolume = true };
557
558 TwoPhaseOnePointTemporalOperator (TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
559 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
560 {
561 }
562
564 void setTime (typename TP::Traits::RangeFieldType t)
565 {
566 time = t;
567 }
568
569 // volume integral depending on test and ansatz functions
570 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
571 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
572 {
573 // Reference to cell
574 const auto& cell = eg.entity();
575
576 // get geometry
577 auto geo = eg.geometry();
578
579 // cell geometry
580 auto ref_el = referenceElement(geo);
581 auto cell_center_local = ref_el.position(0,0);
582 auto cell_volume = geo.volume();
583
584 auto phi = tp.phi(cell,cell_center_local);
585 auto s_l = tp.s_l(cell,cell_center_local,x(lfsu,gas)-x(lfsu,liquid));
586
587 r.accumulate(lfsv, liquid, scale_l * phi * s_l * tp.nu_l(cell,cell_center_local,x(lfsu,liquid)) * cell_volume);
588 r.accumulate(lfsv, gas, scale_g * phi * (1-s_l) * tp.nu_g(cell,cell_center_local,x(lfsu,gas)) * cell_volume);
589 }
590
591 private:
592 TP& tp;
593 typename TP::Traits::RangeFieldType time;
594 Real scale_l, scale_g;
595 };
596
597
606 template<typename TP, typename PL, typename PG>
607 class V_l
608 : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
609 typename PL::Traits::RangeFieldType,
610 PL::Traits::GridViewType::dimension,
611 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
612 V_l<TP,PL,PG> >
613 {
614 // extract useful types
615 using GV = typename PL::Traits::GridViewType;
616 using DF = typename GV::Grid::ctype;
617 using RF = typename PL::Traits::RangeFieldType;
618 using RangeType = typename PL::Traits::RangeType;
619 enum { dim = PL::Traits::GridViewType::dimension };
620 using Element = typename GV::Traits::template Codim<0>::Entity;
621 using IntersectionIterator = typename GV::IntersectionIterator;
622 using Intersection = typename IntersectionIterator::Intersection;
623
624 const TP& tp;
625 const PL& pl;
626 const PG& pg;
627 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
628 typename TP::Traits::RangeFieldType time;
629
630
631 using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
632
633 public:
634 using Traits = Dune::PDELab::GridFunctionTraits<GV,RF,dim,Dune::FieldVector<RF,dim> >;
635
636 using BaseT = Dune::PDELab::GridFunctionBase<Traits,V_l<TP,PL,PG> >;
637
638 V_l (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
639
640 // set time where operator is to be evaluated (i.e. end of the time intervall)
641 void set_time (typename TP::Traits::RangeFieldType time_)
642 {
643 time = time_;
644 }
645
646 inline void evaluate (const typename Traits::ElementType& e,
647 const typename Traits::DomainType& x,
648 typename Traits::RangeType& y) const
649 {
650 // get geometries
651 auto geo = e.geometry();
652
653 // face geometry
654 auto ref_el = referenceElement(geo);
655 auto face_local = ref_el.position(0,0);
656 auto face_volume = geo.volume();
657
658 // cell geometry
659 auto inside_cell_center_local = ref_el.position(0,0);
660 auto inside_cell_center_global = geo.global(inside_cell_center_local);
661
662 // absolute permeability
663 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
664
665 // pressure evaluation
666 typename PL::Traits::RangeType pl_inside, pg_inside;
667 pl.evaluate(e,inside_cell_center_local,pl_inside);
668 pg.evaluate(e,inside_cell_center_local,pg_inside);
669
670 // density
671 auto nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
672
673 // for coefficient computation
674 RF vn[2*dim]; // normal velocities
675 RF coeff[2*dim]; // RT0 coefficient
676 auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear
677 auto determinant = B.determinant();
678
679 // loop over cell neighbors
680 for (const auto& intersection : intersections(pl.getGridView(),e))
681 {
682 // set to zero for processor boundary
683 vn[intersection.indexInInside()] = 0.0;
684
685 // get geometry
686 auto geo_intersection = intersection.geometry();
687
688 // face geometry
689 const auto& face_local = referenceElement(geo_intersection).position(0,0);
690
691 // interior face
692 if (intersection.neighbor())
693 {
694 // get geometry
695 auto geo_outside = intersection.outside().geometry();
696 auto ref_el_outside = referenceElement(geo_outside);
697 auto outside_cell_center_local = ref_el_outside.position(0,0);
698 auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
699
700 // distance of cell centers
701 auto d = outside_cell_center_global;
702 d -= inside_cell_center_global;
703 auto distance = d.two_norm();
704
705 // absolute permeability
706 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
707
708 // gravity times normal
709 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
710
711 // pressure evaluation
712 typename PL::Traits::RangeType pl_outside, pg_outside;
713 pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
714 pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
715
716 // density
717 auto nu_l_outside = tp.nu_l(*(intersection.outside()),outside_cell_center_local,pg_outside);
718
719 // liquid phase calculation
720 auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
721 auto rho_l_outside = tp.rho_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
722 auto w_l = (pl_inside-pl_outside)/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction
723 RF pc_upwind, s_l_upwind, s_g_upwind;
724 if (w_l>=0) // upwind capillary pressure on face
725 {
726 pc_upwind = pg_inside-pl_inside;
727 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
728 }
729 else
730 {
731 pc_upwind = pg_outside-pl_outside;
732 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
733 }
734 s_g_upwind = 1-s_l_upwind;
735 auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l_upwind)/
736 tp.mu_l(e,inside_cell_center_local,pl_inside);
737 auto lambda_l_outside = tp.kr_l(*(intersection.outside()),outside_cell_center_local,s_l_upwind)/
738 tp.mu_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
739 auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
740 auto nu_l = aavg(nu_l_inside,nu_l_outside);
741
742 // set coefficient
743 vn[intersection.indexInInside()] = nu_l * sigma_l * w_l;
744 }
745
746 // boundary face
747 if (intersection.boundary())
748 {
749 // distance of cell center to boundary
750 auto d = geo_intersection.global(face_local);
751 d -= inside_cell_center_global;
752 auto distance = d.two_norm();
753
754 // evaluate boundary condition type
755 auto bc_l = tp.bc_l(intersection,face_local,time);
756
757 // liquid phase Dirichlet boundary
758 if (bc_l==1)
759 {
760 auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
761 auto g_l = tp.g_l(intersection,face_local,time);
762 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
763 auto w_l = (pl_inside-g_l)/distance + rho_l_inside*gn;
764 auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
765 auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l)/
766 tp.mu_l(e,inside_cell_center_local,pl_inside);
767 auto sigma_l = lambda_l_inside*k_abs_inside;
768 vn[intersection.indexInInside()] = nu_l_inside * sigma_l * w_l;
769 }
770
771 // liquid phase Neumann boundary
772 if (bc_l==0)
773 {
774 auto j_l = tp.j_l(intersection,face_local,time);
775 vn[intersection.indexInInside()] = j_l;
776 }
777 }
778
779 // compute coefficient
780 auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element
781 vstar *= vn[intersection.indexInInside()];
782 Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
783 if (intersection.indexInInside()%2==0)
784 normalhat[intersection.indexInInside()/2] = -1.0;
785 else
786 normalhat[intersection.indexInInside()/2] = 1.0;
787 Dune::FieldVector<DF,dim> vstarhat(0);
788 B.umtv(vstar,vstarhat); // Piola backward transformation
789 vstarhat *= determinant;
790 coeff[intersection.indexInInside()] = vstarhat*normalhat;
791 }
792
793 // compute velocity on reference element
794 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
795 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
796 typename Traits::RangeType yhat(0);
797 for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
798 yhat.axpy(coeff[i],rt0vectors[i]);
799
800 // apply Piola transformation
801 B.invert();
802 y = 0;
803 B.umtv(yhat,y);
804 y /= determinant;
805
806 // std::cout << "vn= " ;
807 // for (int i=0; i<2*dim; i++) std::cout << vn[i] << " ";
808 // std::cout << std::endl;
809 // std::cout << "V_l: x=" << x << " y=" << y << std::endl;
810 }
811
812 inline const typename Traits::GridViewType& getGridView ()
813 {
814 return pl.getGridView();
815 }
816
817 private:
818
819 template<typename T>
820 T aavg (T a, T b) const
821 {
822 return 0.5*(a+b);
823 }
824
825 template<typename T>
826 T havg (T a, T b) const
827 {
828 T eps = 1E-30;
829 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
830 }
831 };
832
841 template<typename TP, typename PL, typename PG>
842 class V_g
843 : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
844 typename PL::Traits::RangeFieldType,
845 PL::Traits::GridViewType::dimension,
846 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
847 V_g<TP,PL,PG> >
848 {
849 // extract useful types
850 using GV = typename PL::Traits::GridViewType;
851 using DF = typename GV::Grid::ctype;
852 using RF = typename PL::Traits::RangeFieldType;
853 using RangeType = typename PL::Traits::RangeType;
854 enum { dim = PL::Traits::GridViewType::dimension };
855 using Element = typename GV::Traits::template Codim<0>::Entity;
856 using IntersectionIterator = typename GV::IntersectionIterator;
857 using Intersection = typename IntersectionIterator::Intersection;
858
859 const TP& tp;
860 const PL& pl;
861 const PG& pg;
862 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
863 typename TP::Traits::RangeFieldType time;
864
865
866 using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
867
868 public:
869 using Traits = Dune::PDELab::GridFunctionTraits<GV,RF,dim,Dune::FieldVector<RF,dim> >;
870
871 using BaseT = Dune::PDELab::GridFunctionBase<Traits,V_g<TP,PL,PG> >;
872
873 V_g (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
874
875 // set time where operator is to be evaluated (i.e. end of the time intervall)
876 void set_time (typename TP::Traits::RangeFieldType time_)
877 {
878 time = time_;
879 }
880
881 inline void evaluate (const typename Traits::ElementType& e,
882 const typename Traits::DomainType& x,
883 typename Traits::RangeType& y) const
884 {
885 // get geometries
886 auto geo = e.geometry();
887
888 // face geometry
889 auto ref_el = referenceElement(geo);
890 auto face_local = ref_el.position(0,0);
891 auto face_volume = geo.volume();
892
893 // cell geometry
894 auto inside_cell_center_local = ref_el.position(0,0);
895 auto inside_cell_center_global = geo.global(inside_cell_center_local);
896
897 // absolute permeability
898 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
899
900 // pressure evaluation
901 typename PL::Traits::RangeType pl_inside, pg_inside;
902 pl.evaluate(e,inside_cell_center_local,pl_inside);
903 pg.evaluate(e,inside_cell_center_local,pg_inside);
904
905 // density evaluation
906 auto nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
907
908 // for coefficient computation
909 RF vn[2*dim]; // normal velocities
910 RF coeff[2*dim]; // RT0 coefficient
911 auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear
912 auto determinant = B.determinant();
913
914 // loop over cell neighbors
915 for (const auto& intersection : intersections(pl.getGridView(),e))
916 {
917 // set to zero for processor boundary
918 vn[intersection.indexInInside()] = 0.0;
919
920 // get geometry
921 auto geo_intersection = intersection.geometry();
922
923 // face geometry
924 const auto& face_local = referenceElement(geo_intersection).position(0,0);
925
926 // interior face
927 if (intersection.neighbor())
928 {
929 // get geometry
930 auto geo_outside = intersection.outside().geometry();
931 auto ref_el_outside = referenceElement(geo_outside);
932 auto outside_cell_center_local = ref_el_outside.position(0,0);
933 auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
934
935 // distance of cell centers
936 auto d = outside_cell_center_global;
937 d -= inside_cell_center_global;
938 auto distance = d.two_norm();
939
940 // absolute permeability
941 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
942
943 // gravity times normal
944 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
945
946 // pressure evaluation
947 typename PL::Traits::RangeType pl_outside, pg_outside;
948 pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
949 pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
950
951 // needed for both phases
952 RF pc_upwind, s_l_upwind, s_g_upwind;
953
954 // gas phase calculation
955 auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
956 auto rho_g_outside = tp.rho_g(e,outside_cell_center_local,pg_outside);
957 auto w_g = (pg_inside-pg_outside)/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction
958 if (w_g>=0) // upwind capillary pressure on face
959 {
960 pc_upwind = pg_inside-pl_inside;
961 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
962 }
963 else
964 {
965 pc_upwind = pg_outside-pl_outside;
966 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
967 }
968 s_g_upwind = 1-s_l_upwind;
969 auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g_upwind)/
970 tp.mu_g(e,inside_cell_center_local,pg_inside);
971 auto lambda_g_outside = tp.kr_g(*(intersection.outside()),outside_cell_center_local,s_g_upwind)/
972 tp.mu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
973 auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
974
975 auto nu_g_outside = tp.nu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
976
977 // set coefficient
978 vn[intersection.indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
979 }
980
981 // boundary face
982 if (intersection.boundary())
983 {
984 // distance of cell center to boundary
985 auto d = geo_intersection.global(face_local);
986 d -= inside_cell_center_global;
987 auto distance = d.two_norm();
988
989 // evaluate boundary condition type
990 auto bc_g = tp.bc_g(intersection,face_local,time);
991
992 // gas phase Dirichlet boundary
993 if (bc_g==1)
994 {
995 auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
996 auto g_g = tp.g_g(intersection,face_local,time);
997 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
998 auto w_g = (pg_inside-g_g)/distance + rho_g_inside*gn;
999 auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
1000 auto s_g = 1-s_l;
1001 auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g)/
1002 tp.mu_g(e,inside_cell_center_local,pg_inside);
1003 auto sigma_g = lambda_g_inside*k_abs_inside;
1004
1005 vn[intersection.indexInInside()] = nu_g_inside * sigma_g * w_g;
1006 }
1007
1008 // gas phase Neumann boundary
1009 if (bc_g==0)
1010 {
1011 auto j_g = tp.j_g(intersection,face_local,time);
1012 vn[intersection.indexInInside()] = j_g; /* /nu_g_inside*/;
1013 }
1014 }
1015
1016 // compute coefficient
1017 auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element
1018 vstar *= vn[intersection.indexInInside()];
1019 Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
1020 if (intersection.indexInInside()%2==0)
1021 normalhat[intersection.indexInInside()/2] = -1.0;
1022 else
1023 normalhat[intersection.indexInInside()/2] = 1.0;
1024 Dune::FieldVector<DF,dim> vstarhat(0);
1025 B.umtv(vstar,vstarhat); // Piola backward transformation
1026 vstarhat *= determinant;
1027 coeff[intersection.indexInInside()] = vstarhat*normalhat;
1028 }
1029
1030 // compute velocity on reference element
1031 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1032 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1033 typename Traits::RangeType yhat(0);
1034 for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1035 yhat.axpy(coeff[i],rt0vectors[i]);
1036
1037 // apply Piola transformation
1038 B.invert();
1039 y = 0;
1040 B.umtv(yhat,y);
1041 y /= determinant;
1042 }
1043
1044 inline const typename Traits::GridViewType& getGridView ()
1045 {
1046 return pl.getGridView();
1047 }
1048
1049 private:
1050
1051 template<typename T>
1052 T aavg (T a, T b) const
1053 {
1054 return 0.5*(a+b);
1055 }
1056
1057 template<typename T>
1058 T havg (T a, T b) const
1059 {
1060 T eps = 1E-30;
1061 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
1062 }
1063 };
1064
1065
1066 }
1067}
1068
1069#endif // DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:117
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Dune::IntersectionIterator::Intersection
Dune::Intersection< GridImp, IntersectionImp > Intersection
Type of Intersection this IntersectionIterator points to.
Definition: intersectioniterator.hh:111
Dune::PDELab::FullVolumePattern
sparsity pattern generator
Definition: pattern.hh:14
Dune::PDELab::GridFunctionBase
leaf of a function tree
Definition: function.hh:302
Dune::PDELab::InstationaryLocalOperatorDefaultMethods
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Dune::PDELab::LocalOperatorDefaultFlags
Default flags for all local operators.
Definition: flags.hh:19
Dune::PDELab::NumericalJacobianApplyVolume
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Dune::PDELab::NumericalJacobianVolume
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
Dune::PDELab::TwoPhaseOnePointTemporalOperator
Definition: twophaseccfv.hh:544
Dune::PDELab::TwoPhaseOnePointTemporalOperator::setTime
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:564
Dune::PDELab::TwoPhaseParameterInterface
base class for parameter class
Definition: twophaseccfv.hh:70
Dune::PDELab::TwoPhaseParameterInterface::bc_g
int bc_g(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
gas phase boundary condition type
Definition: twophaseccfv.hh:184
Dune::PDELab::TwoPhaseParameterInterface::nu_g
Traits::RangeFieldType nu_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_g) const
gas phase molar density
Definition: twophaseccfv.hh:153
Dune::PDELab::TwoPhaseParameterInterface::pc
Traits::RangeFieldType pc(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType s_l) const
capillary pressure function
Definition: twophaseccfv.hh:83
Dune::PDELab::TwoPhaseParameterInterface::j_l
Traits::RangeFieldType j_l(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
liquid phase Neumann boundary condition
Definition: twophaseccfv.hh:205
Dune::PDELab::TwoPhaseParameterInterface::g_l
Traits::RangeFieldType g_l(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
liquid phase Dirichlet boundary condition
Definition: twophaseccfv.hh:191
Dune::PDELab::TwoPhaseParameterInterface::kr_g
Traits::RangeFieldType kr_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType s_g) const
gas phase relative permeability
Definition: twophaseccfv.hh:107
Dune::PDELab::TwoPhaseParameterInterface::k_abs
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:131
Dune::PDELab::TwoPhaseParameterInterface::gravity
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:137
Dune::PDELab::TwoPhaseParameterInterface::nu_l
Traits::RangeFieldType nu_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_l) const
liquid phase molar density
Definition: twophaseccfv.hh:145
Dune::PDELab::TwoPhaseParameterInterface::phi
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:76
Dune::PDELab::TwoPhaseParameterInterface::bc_l
int bc_l(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
liquid phase boundary condition type
Definition: twophaseccfv.hh:177
Dune::PDELab::TwoPhaseParameterInterface::mu_g
Traits::RangeFieldType mu_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_g) const
gas phase viscosity
Definition: twophaseccfv.hh:123
Dune::PDELab::TwoPhaseParameterInterface::rho_g
Traits::RangeFieldType rho_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_g) const
gas phase mass density
Definition: twophaseccfv.hh:169
Dune::PDELab::TwoPhaseParameterInterface::kr_l
Traits::RangeFieldType kr_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType s_l) const
liquid phase relative permeability
Definition: twophaseccfv.hh:99
Dune::PDELab::TwoPhaseParameterInterface::q_g
Traits::RangeFieldType q_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType time) const
gas phase source term
Definition: twophaseccfv.hh:227
Dune::PDELab::TwoPhaseParameterInterface::mu_l
Traits::RangeFieldType mu_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_l) const
liquid phase viscosity
Definition: twophaseccfv.hh:115
Dune::PDELab::TwoPhaseParameterInterface::q_l
Traits::RangeFieldType q_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType time) const
liquid phase source term
Definition: twophaseccfv.hh:219
Dune::PDELab::TwoPhaseParameterInterface::rho_l
Traits::RangeFieldType rho_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_l) const
liquid phase mass density
Definition: twophaseccfv.hh:161
Dune::PDELab::TwoPhaseParameterInterface::s_l
Traits::RangeFieldType s_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType pc) const
inverse capillary pressure function
Definition: twophaseccfv.hh:91
Dune::PDELab::TwoPhaseParameterInterface::j_g
Traits::RangeFieldType j_g(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
gas phase Neumann boundary condition
Definition: twophaseccfv.hh:212
Dune::PDELab::TwoPhaseParameterInterface::g_g
Traits::RangeFieldType g_g(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &x, typename Traits::RangeFieldType time) const
gas phase Dirichlet boundary condition
Definition: twophaseccfv.hh:198
Dune::PDELab::V_g
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:848
Dune::PDELab::V_l
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:613
exceptions.hh
A few common exception classes.
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune::referenceElement
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune::Indices
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
raviartthomascube.hh
Convenience header that includes all available Raviart-Thomas local finite elements for cubes.
Dune::Codim
Static tag representing a codimension.
Definition: dimension.hh:24
Dune::PDELab::GridFunctionTraits
traits class holding the function signature, same as in local function
Definition: function.hh:183
Dune::PDELab::PowerCompositeGridFunctionTraits::ElementType
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:119
Dune::PDELab::PowerCompositeGridFunctionTraits::GridViewType
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:116
Dune::PDELab::TwoPhaseParameterTraits
traits class for two phase parameter class
Definition: twophaseccfv.hh:24
Dune::PDELab::TwoPhaseParameterTraits::RangeFieldType
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:44
Dune::PDELab::TwoPhaseParameterTraits::GridViewType
GV GridViewType
the grid view
Definition: twophaseccfv.hh:26
Dune::PDELab::TwoPhaseParameterTraits::DomainFieldType
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:35
Dune::PDELab::TwoPhaseParameterTraits::PermTensorType
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:50
Dune::PDELab::TwoPhaseParameterTraits::dimDomain
@ dimDomain
dimension of the domain
Definition: twophaseccfv.hh:31
Dune::PDELab::TwoPhaseParameterTraits::ElementType
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:53
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)