2 #ifndef DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH 3 #define DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH 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> 13 #include"../common/function.hh" 22 template<
typename GV,
typename RF>
38 using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>;
47 using RangeType = Dune::FieldVector<RF,GV::dimensionworld>;
53 using ElementType =
typename GV::Traits::template Codim<0>::Entity;
57 template<
typename GV,
typename RF>
64 using PermTensorType = Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>;
68 template<
class T,
class Imp>
75 typename Traits::RangeFieldType
76 phi (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const 78 return asImp().phi(e,x);
82 typename Traits::RangeFieldType
83 pc (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
84 typename Traits::RangeFieldType s_l)
const 86 return asImp().pc(e,x,s_l);
90 typename Traits::RangeFieldType
91 s_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
92 typename Traits::RangeFieldType pc)
const 94 return asImp().s_l(e,x,pc);
98 typename Traits::RangeFieldType
99 kr_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
100 typename Traits::RangeFieldType s_l)
const 102 return asImp().kr_l(e,x,s_l);
106 typename Traits::RangeFieldType
107 kr_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
108 typename Traits::RangeFieldType s_g)
const 110 return asImp().kr_g(e,x,s_g);
114 typename Traits::RangeFieldType
115 mu_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
116 typename Traits::RangeFieldType p_l)
const 118 return asImp().mu_l(e,x,p_l);
122 typename Traits::RangeFieldType
123 mu_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
124 typename Traits::RangeFieldType p_g)
const 126 return asImp().mu_l(e,x,p_g);
130 typename Traits::PermTensorType
131 k_abs (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const 133 return asImp().k_abs(e,x);
137 const typename Traits::RangeType&
gravity ()
const 139 return asImp().gravity();
144 typename Traits::RangeFieldType
145 nu_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
146 typename Traits::RangeFieldType p_l)
const 148 return asImp().nu_l(e,x,p_l);
152 typename Traits::RangeFieldType
153 nu_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
154 typename Traits::RangeFieldType p_g)
const 156 return asImp().nu_g(e,x,p_g);
160 typename Traits::RangeFieldType
161 rho_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
162 typename Traits::RangeFieldType p_l)
const 164 return asImp().rho_l(e,x,p_l);
168 typename Traits::RangeFieldType
169 rho_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
170 typename Traits::RangeFieldType p_g)
const 172 return asImp().rho_g(e,x,p_g);
179 return asImp().bc_l(is,x,time);
186 return asImp().bc_g(is,x,time);
190 typename Traits::RangeFieldType
193 return asImp().g_l(is,x,time);
197 typename Traits::RangeFieldType
200 return asImp().g_g(is,x,time);
204 typename Traits::RangeFieldType
207 return asImp().j_l(is,x,time);
211 typename Traits::RangeFieldType
214 return asImp().j_g(is,x,time);
218 typename Traits::RangeFieldType
219 q_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
220 typename Traits::RangeFieldType time)
const 222 return asImp().q_l(e,x,time);
226 typename Traits::RangeFieldType
227 q_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
228 typename Traits::RangeFieldType time)
const 230 return asImp().q_g(e,x,time);
234 Imp& asImp () {
return static_cast<Imp &
> (*this);}
235 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
243 template<
typename TP>
257 enum {
dim = TP::Traits::GridViewType::dimension };
261 using Real =
typename TP::Traits::RangeFieldType;
264 enum { doPatternVolume =
true };
265 enum { doPatternSkeleton =
true };
268 enum { doAlphaSkeleton =
true };
269 enum { doAlphaBoundary =
true };
270 enum { doLambdaVolume =
true };
271 enum { doLambdaBoundary =
true };
275 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
279 template<
typename EG,
typename LFSV,
typename R>
283 const auto& cell = eg.entity();
286 auto geo = eg.geometry();
289 auto ref_el = referenceElement(geo);
290 auto cell_center_local = ref_el.position(0,0);
291 auto cell_volume = geo.volume();
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);
300 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
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 307 using namespace Indices;
308 using PLSpace = TypeTree::Child<LFSV,_0>;
309 using RF =
typename PLSpace::Traits::FiniteElementType::
310 Traits::LocalBasisType::Traits::RangeFieldType;
313 const auto& cell_inside = ig.inside();
314 const auto& cell_outside = ig.outside();
317 auto geo = ig.geometry();
318 auto geo_inside = cell_inside.geometry();
319 auto geo_outside = cell_outside.geometry();
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();
330 auto d = outside_cell_center_global;
331 d -= inside_cell_center_global;
332 auto distance = d.two_norm();
335 auto ref_el = referenceElement(geo);
336 auto face_local = ref_el.position(0,0);
337 auto face_volume = geo.volume();
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);
344 auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
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;
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)));
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);
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);
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);
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);
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;
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)));
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);
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);
391 s_g_upwind = 1-s_l_upwind;
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);
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);
405 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
407 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
411 const auto& cell_inside = ig.inside();
414 auto geo = ig.geometry();
415 auto geo_inside = cell_inside.geometry();
418 auto ref_el = referenceElement(geo);
419 auto face_local = ref_el.position(0,0);
420 auto face_volume = geo.volume();
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;
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();
433 auto d = geo.global(face_local);
434 d -= inside_cell_center_global;
435 auto distance = d.two_norm();
438 auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
441 auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
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);
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));
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);
474 template<
typename IG,
typename LFSV,
typename R>
478 auto geo = ig.geometry();
481 auto ref_el = referenceElement(geo);
482 auto face_local = ref_el.position(0,0);
483 auto face_volume = geo.volume();
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;
493 auto j_l = tp.j_l(ig.intersection(),face_local,time);
494 r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
500 auto j_g = tp.j_g(ig.intersection(),face_local,time);
501 r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
506 void setTime (
typename TP::Traits::RangeFieldType t)
513 typename TP::Traits::RangeFieldType time;
514 Real scale_l, scale_g;
517 T aavg (T a, T b)
const 523 T havg (T a, T b)
const 526 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
545 enum {
dim = TP::Traits::GridViewType::dimension };
549 using Real =
typename TP::Traits::RangeFieldType;
553 enum { doPatternVolume =
true };
556 enum { doAlphaVolume =
true };
559 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
564 void setTime (
typename TP::Traits::RangeFieldType t)
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 574 const auto& cell = eg.entity();
577 auto geo = eg.geometry();
580 auto ref_el = referenceElement(geo);
581 auto cell_center_local = ref_el.position(0,0);
582 auto cell_volume = geo.volume();
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));
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);
593 typename TP::Traits::RangeFieldType time;
594 Real scale_l, scale_g;
606 template<
typename TP,
typename PL,
typename PG>
609 typename PL::Traits::RangeFieldType,
610 PL::Traits::GridViewType::dimension,
611 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
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;
627 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
628 typename TP::Traits::RangeFieldType time;
631 using RT0RangeType =
typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
638 V_l (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
641 void set_time (
typename TP::Traits::RangeFieldType time_)
651 auto geo = e.geometry();
654 auto ref_el = referenceElement(geo);
655 auto face_local = ref_el.position(0,0);
656 auto face_volume = geo.volume();
659 auto inside_cell_center_local = ref_el.position(0,0);
660 auto inside_cell_center_global = geo.global(inside_cell_center_local);
663 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
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);
671 auto nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
676 auto B = geo.jacobianInverseTransposed(x);
677 auto determinant = B.determinant();
680 for (
const auto& intersection : intersections(pl.getGridView(),
e))
683 vn[intersection.indexInInside()] = 0.0;
686 auto geo_intersection = intersection.geometry();
689 const auto& face_local = referenceElement(geo_intersection).position(0,0);
692 if (intersection.neighbor())
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);
701 auto d = outside_cell_center_global;
702 d -= inside_cell_center_global;
703 auto distance = d.two_norm();
706 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
709 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
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);
717 auto nu_l_outside = tp.nu_l(*(intersection.outside()),outside_cell_center_local,pg_outside);
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;
723 RF pc_upwind, s_l_upwind, s_g_upwind;
726 pc_upwind = pg_inside-pl_inside;
727 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
731 pc_upwind = pg_outside-pl_outside;
732 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
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);
743 vn[intersection.indexInInside()] = nu_l * sigma_l * w_l;
747 if (intersection.boundary())
750 auto d = geo_intersection.global(face_local);
751 d -= inside_cell_center_global;
752 auto distance = d.two_norm();
755 auto bc_l = tp.bc_l(intersection,face_local,time);
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;
774 auto j_l = tp.j_l(intersection,face_local,time);
775 vn[intersection.indexInInside()] = j_l;
780 auto vstar = intersection.unitOuterNormal(face_local);
781 vstar *= vn[intersection.indexInInside()];
782 Dune::FieldVector<RF,dim> normalhat(0);
783 if (intersection.indexInInside()%2==0)
784 normalhat[intersection.indexInInside()/2] = -1.0;
786 normalhat[intersection.indexInInside()/2] = 1.0;
787 Dune::FieldVector<DF,dim> vstarhat(0);
788 B.umtv(vstar,vstarhat);
789 vstarhat *= determinant;
790 coeff[intersection.indexInInside()] = vstarhat*normalhat;
794 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
795 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
797 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
798 yhat.axpy(coeff[i],rt0vectors[i]);
814 return pl.getGridView();
820 T aavg (T a, T b)
const 826 T havg (T a, T b)
const 829 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
841 template<
typename TP,
typename PL,
typename PG>
844 typename PL::Traits::RangeFieldType,
845 PL::Traits::GridViewType::dimension,
846 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
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;
862 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
863 typename TP::Traits::RangeFieldType time;
866 using RT0RangeType =
typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
873 V_g (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
876 void set_time (
typename TP::Traits::RangeFieldType time_)
886 auto geo = e.geometry();
889 auto ref_el = referenceElement(geo);
890 auto face_local = ref_el.position(0,0);
891 auto face_volume = geo.volume();
894 auto inside_cell_center_local = ref_el.position(0,0);
895 auto inside_cell_center_global = geo.global(inside_cell_center_local);
898 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
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);
906 auto nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
911 auto B = geo.jacobianInverseTransposed(x);
912 auto determinant = B.determinant();
915 for (
const auto& intersection : intersections(pl.getGridView(),
e))
918 vn[intersection.indexInInside()] = 0.0;
921 auto geo_intersection = intersection.geometry();
924 const auto& face_local = referenceElement(geo_intersection).position(0,0);
927 if (intersection.neighbor())
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);
936 auto d = outside_cell_center_global;
937 d -= inside_cell_center_global;
938 auto distance = d.two_norm();
941 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
944 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
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);
952 RF pc_upwind, s_l_upwind, s_g_upwind;
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;
960 pc_upwind = pg_inside-pl_inside;
961 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
965 pc_upwind = pg_outside-pl_outside;
966 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
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);
975 auto nu_g_outside = tp.nu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
978 vn[intersection.indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
982 if (intersection.boundary())
985 auto d = geo_intersection.global(face_local);
986 d -= inside_cell_center_global;
987 auto distance = d.two_norm();
990 auto bc_g = tp.bc_g(intersection,face_local,time);
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);
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;
1005 vn[intersection.indexInInside()] = nu_g_inside * sigma_g * w_g;
1011 auto j_g = tp.j_g(intersection,face_local,time);
1012 vn[intersection.indexInInside()] = j_g; ;
1017 auto vstar = intersection.unitOuterNormal(face_local);
1018 vstar *= vn[intersection.indexInInside()];
1019 Dune::FieldVector<RF,dim> normalhat(0);
1020 if (intersection.indexInInside()%2==0)
1021 normalhat[intersection.indexInInside()/2] = -1.0;
1023 normalhat[intersection.indexInInside()/2] = 1.0;
1024 Dune::FieldVector<DF,dim> vstarhat(0);
1025 B.umtv(vstar,vstarhat);
1026 vstarhat *= determinant;
1027 coeff[intersection.indexInInside()] = vstarhat*normalhat;
1031 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1032 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1034 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1035 yhat.axpy(coeff[i],rt0vectors[i]);
1046 return pl.getGridView();
1051 template<
typename T>
1052 T aavg (T a, T b)
const 1057 template<
typename T>
1058 T havg (T a, T b)
const 1061 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
1069 #endif // DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH dimension of the domain
Definition: twophaseccfv.hh:31
R RangeType
range type
Definition: function.hh:61
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:44
TwoPhaseTwoPointFluxOperator(const TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
constructor: pass parameter object
Definition: twophaseccfv.hh:274
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:876
static const int dim
Definition: adaptivity.hh:84
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: twophaseccfv.hh:38
const Entity & e
Definition: localfunctionspace.hh:120
leaf of a function tree
Definition: function.hh:298
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:29
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
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:571
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
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
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:881
traits class for two phase parameter class
Definition: twophaseccfv.hh:23
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
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
sparsity pattern generator
Definition: pattern.hh:13
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
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:607
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
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:1044
base class for parameter class
Definition: twophaseccfv.hh:69
IntersectionType
Enum describing the type of an intersection.
Definition: intersectiontype.hh:14
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
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:137
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:812
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:49
V_l(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:638
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:118
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:564
Definition: twophaseccfv.hh:538
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:506
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
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:50
T Traits
Definition: twophaseccfv.hh:72
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
traits class holding the function signature, same as in local function
Definition: function.hh:176
TwoPhaseOnePointTemporalOperator(TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
Definition: twophaseccfv.hh:558
V_g(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:873
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: twophaseccfv.hh:406
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
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:35
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
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
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
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:646
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:53
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:641
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r_s) const
Definition: twophaseccfv.hh:475
const IG & ig
Definition: constraints.hh:148
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:131
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:280
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:76
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: twophaseccfv.hh:47
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
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:115
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: twophaseccfv.hh:301
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
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
typename GV::Intersection IntersectionType
Definition: twophaseccfv.hh:54
Definition: twophaseccfv.hh:58
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: twophaseccfv.hh:41
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
GV GridViewType
the grid view
Definition: twophaseccfv.hh:26
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
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:842
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
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
Definition: twophaseccfv.hh:244