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> 14 #include"../common/function.hh" 23 template<
typename GV,
typename RF>
39 using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>;
48 using RangeType = Dune::FieldVector<RF,GV::dimensionworld>;
54 using ElementType =
typename GV::Traits::template Codim<0>::Entity;
58 template<
typename GV,
typename RF>
65 using PermTensorType = Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>;
69 template<
class T,
class Imp>
76 typename Traits::RangeFieldType
77 phi (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const 79 return asImp().phi(e,x);
83 typename Traits::RangeFieldType
84 pc (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
85 typename Traits::RangeFieldType s_l)
const 87 return asImp().pc(e,x,s_l);
91 typename Traits::RangeFieldType
92 s_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
93 typename Traits::RangeFieldType pc)
const 95 return asImp().s_l(e,x,pc);
99 typename Traits::RangeFieldType
100 kr_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
101 typename Traits::RangeFieldType s_l)
const 103 return asImp().kr_l(e,x,s_l);
107 typename Traits::RangeFieldType
108 kr_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
109 typename Traits::RangeFieldType s_g)
const 111 return asImp().kr_g(e,x,s_g);
115 typename Traits::RangeFieldType
116 mu_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
117 typename Traits::RangeFieldType p_l)
const 119 return asImp().mu_l(e,x,p_l);
123 typename Traits::RangeFieldType
124 mu_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
125 typename Traits::RangeFieldType p_g)
const 127 return asImp().mu_l(e,x,p_g);
131 typename Traits::PermTensorType
132 k_abs (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const 134 return asImp().k_abs(e,x);
138 const typename Traits::RangeType&
gravity ()
const 140 return asImp().gravity();
145 typename Traits::RangeFieldType
146 nu_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
147 typename Traits::RangeFieldType p_l)
const 149 return asImp().nu_l(e,x,p_l);
153 typename Traits::RangeFieldType
154 nu_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
155 typename Traits::RangeFieldType p_g)
const 157 return asImp().nu_g(e,x,p_g);
161 typename Traits::RangeFieldType
162 rho_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
163 typename Traits::RangeFieldType p_l)
const 165 return asImp().rho_l(e,x,p_l);
169 typename Traits::RangeFieldType
170 rho_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
171 typename Traits::RangeFieldType p_g)
const 173 return asImp().rho_g(e,x,p_g);
180 return asImp().bc_l(is,x,time);
187 return asImp().bc_g(is,x,time);
191 typename Traits::RangeFieldType
194 return asImp().g_l(is,x,time);
198 typename Traits::RangeFieldType
201 return asImp().g_g(is,x,time);
205 typename Traits::RangeFieldType
208 return asImp().j_l(is,x,time);
212 typename Traits::RangeFieldType
215 return asImp().j_g(is,x,time);
219 typename Traits::RangeFieldType
220 q_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
221 typename Traits::RangeFieldType time)
const 223 return asImp().q_l(e,x,time);
227 typename Traits::RangeFieldType
228 q_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
229 typename Traits::RangeFieldType time)
const 231 return asImp().q_g(e,x,time);
235 Imp& asImp () {
return static_cast<Imp &
> (*this);}
236 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
244 template<
typename TP>
258 enum {
dim = TP::Traits::GridViewType::dimension };
262 using Real =
typename TP::Traits::RangeFieldType;
265 enum { doPatternVolume =
true };
266 enum { doPatternSkeleton =
true };
269 enum { doAlphaSkeleton =
true };
270 enum { doAlphaBoundary =
true };
271 enum { doLambdaVolume =
true };
272 enum { doLambdaBoundary =
true };
276 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
280 template<
typename EG,
typename LFSV,
typename R>
284 const auto& cell = eg.entity();
287 auto geo = eg.geometry();
291 auto cell_center_local = ref_el.position(0,0);
292 auto cell_volume = geo.volume();
295 r.accumulate(lfsv, liquid, -scale_l * tp.q_l(cell,cell_center_local,time) * cell_volume);
296 r.accumulate(lfsv, gas, -scale_g * tp.q_g(cell,cell_center_local,time) * cell_volume);
301 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
303 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
304 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
305 R& r_s, R& r_n)
const 308 using namespace Indices;
309 using PLSpace = TypeTree::Child<LFSV,_0>;
310 using RF =
typename PLSpace::Traits::FiniteElementType::
311 Traits::LocalBasisType::Traits::RangeFieldType;
314 const auto& cell_inside = ig.inside();
315 const auto& cell_outside = ig.outside();
318 auto geo = ig.geometry();
319 auto geo_inside = cell_inside.geometry();
320 auto geo_outside = cell_outside.geometry();
325 auto inside_cell_center_local = ref_el_inside.position(0,0);
326 auto outside_cell_center_local = ref_el_outside.position(0,0);
327 auto inside_cell_center_global = geo_inside.center();
328 auto outside_cell_center_global = geo_outside.center();
331 auto d = outside_cell_center_global;
332 d -= inside_cell_center_global;
333 auto distance = d.two_norm();
337 auto face_local = ref_el.position(0,0);
338 auto face_volume = geo.volume();
341 auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
342 auto k_abs_outside = tp.k_abs(cell_outside,outside_cell_center_local);
345 auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
348 auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
349 auto rho_l_outside = tp.rho_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
350 auto w_l = (x_s(lfsu_s,liquid)-x_n(lfsu_n,liquid))/distance + aavg(rho_l_inside,rho_l_outside)*gn;
351 RF pc_upwind, s_l_upwind, s_g_upwind;
352 auto nu_l = aavg(tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)),
353 tp.nu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid)));
356 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
357 s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
361 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
362 s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
364 s_g_upwind = 1-s_l_upwind;
365 auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l_upwind)/
366 tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
367 auto lambda_l_outside = tp.kr_l(cell_outside,outside_cell_center_local,s_l_upwind)/
368 tp.mu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
369 auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
371 r_s.accumulate(lfsv_s,liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
372 r_n.accumulate(lfsv_n,liquid, -scale_l * nu_l * sigma_l * w_l * face_volume);
375 auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
376 auto rho_g_outside = tp.rho_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
377 auto w_g = (x_s(lfsu_s,gas)-x_n(lfsu_n,gas))/distance + aavg(rho_g_inside,rho_g_outside)*gn;
378 auto nu_g = aavg(tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)),
379 tp.nu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas)));
384 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
385 s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
389 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
390 s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
392 s_g_upwind = 1-s_l_upwind;
394 auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g_upwind)/
395 tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
396 auto lambda_g_outside = tp.kr_g(cell_outside,outside_cell_center_local,s_g_upwind)/
397 tp.mu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
398 auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
400 r_s.accumulate(lfsv_s, gas, scale_g * nu_g * sigma_g * w_g * face_volume);
401 r_n.accumulate(lfsv_n, gas, -scale_g * nu_g * sigma_g * w_g * face_volume);
406 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
408 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
412 const auto& cell_inside = ig.inside();
415 auto geo = ig.geometry();
416 auto geo_inside = cell_inside.geometry();
420 auto face_local = ref_el.position(0,0);
421 auto face_volume = geo.volume();
424 auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
425 auto bc_g = tp.bc_g(ig.intersection(),face_local,time);
426 if (bc_l!=1 && bc_g!=1)
return;
430 auto inside_cell_center_local = ref_el_inside.position(0,0);
431 auto inside_cell_center_global = geo_inside.center();
434 auto d = geo.global(face_local);
435 d -= inside_cell_center_global;
436 auto distance = d.two_norm();
439 auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
442 auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
447 auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
448 auto g_l = tp.g_l(ig.intersection(),face_local,time);
449 auto w_l = (x_s(lfsu_s,liquid)-g_l)/distance + rho_l_inside*gn;
450 auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
451 auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l)/
452 tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
453 auto sigma_l = lambda_l_inside*k_abs_inside;
454 auto nu_l = tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
455 r_s.accumulate(lfsv_s, liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
461 auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
462 auto g_g = tp.g_g(ig.intersection(),face_local,time);
463 auto w_g = (x_s(lfsu_s,gas)-g_g)/distance + rho_g_inside*gn;
464 auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
466 auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g)/
467 tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
468 auto sigma_g = lambda_g_inside*k_abs_inside;
469 auto nu_g = tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
470 r_s.accumulate(lfsv_s, gas, scale_l * nu_g * sigma_g * w_g * face_volume);
475 template<
typename IG,
typename LFSV,
typename R>
479 auto geo = ig.geometry();
483 auto face_local = ref_el.position(0,0);
484 auto face_volume = geo.volume();
487 auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
488 auto bc_g = tp.bc_g(ig.intersection(),face_local,time);
489 if (bc_l!=0 && bc_g!=0)
return;
494 auto j_l = tp.j_l(ig.intersection(),face_local,time);
495 r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
501 auto j_g = tp.j_g(ig.intersection(),face_local,time);
502 r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
507 void setTime (
typename TP::Traits::RangeFieldType t)
514 typename TP::Traits::RangeFieldType time;
515 Real scale_l, scale_g;
518 T aavg (T a, T b)
const 524 T havg (T a, T b)
const 527 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
546 enum {
dim = TP::Traits::GridViewType::dimension };
550 using Real =
typename TP::Traits::RangeFieldType;
554 enum { doPatternVolume =
true };
557 enum { doAlphaVolume =
true };
560 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
565 void setTime (
typename TP::Traits::RangeFieldType t)
571 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
572 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 575 const auto& cell = eg.entity();
578 auto geo = eg.geometry();
582 auto cell_center_local = ref_el.position(0,0);
583 auto cell_volume = geo.volume();
585 auto phi = tp.phi(cell,cell_center_local);
586 auto s_l = tp.s_l(cell,cell_center_local,x(lfsu,gas)-x(lfsu,liquid));
588 r.accumulate(lfsv, liquid, scale_l * phi * s_l * tp.nu_l(cell,cell_center_local,x(lfsu,liquid)) * cell_volume);
589 r.accumulate(lfsv, gas, scale_g * phi * (1-s_l) * tp.nu_g(cell,cell_center_local,x(lfsu,gas)) * cell_volume);
594 typename TP::Traits::RangeFieldType time;
595 Real scale_l, scale_g;
607 template<
typename TP,
typename PL,
typename PG>
610 typename PL::Traits::RangeFieldType,
611 PL::Traits::GridViewType::dimension,
612 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
616 using GV =
typename PL::Traits::GridViewType;
617 using DF =
typename GV::Grid::ctype;
618 using RF =
typename PL::Traits::RangeFieldType;
619 using RangeType =
typename PL::Traits::RangeType;
620 enum {
dim = PL::Traits::GridViewType::dimension };
621 using Element =
typename GV::Traits::template Codim<0>::Entity;
622 using IntersectionIterator =
typename GV::IntersectionIterator;
623 using Intersection =
typename IntersectionIterator::Intersection;
628 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
629 typename TP::Traits::RangeFieldType time;
632 using RT0RangeType =
typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
639 V_l (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
642 void set_time (
typename TP::Traits::RangeFieldType time_)
652 auto geo = e.geometry();
656 auto face_local = ref_el.position(0,0);
657 auto face_volume = geo.volume();
660 auto inside_cell_center_local = ref_el.position(0,0);
661 auto inside_cell_center_global = geo.global(inside_cell_center_local);
664 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
667 typename PL::Traits::RangeType pl_inside, pg_inside;
668 pl.evaluate(e,inside_cell_center_local,pl_inside);
669 pg.evaluate(e,inside_cell_center_local,pg_inside);
672 auto nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
677 auto B = geo.jacobianInverseTransposed(x);
678 auto determinant = B.determinant();
681 for (
const auto& intersection : intersections(pl.getGridView(),
e))
684 vn[intersection.indexInInside()] = 0.0;
687 auto geo_intersection = intersection.geometry();
693 if (intersection.neighbor())
696 auto geo_outside = intersection.outside().geometry();
698 auto outside_cell_center_local = ref_el_outside.position(0,0);
699 auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
702 auto d = outside_cell_center_global;
703 d -= inside_cell_center_global;
704 auto distance = d.two_norm();
707 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
710 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
713 typename PL::Traits::RangeType pl_outside, pg_outside;
714 pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
715 pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
718 auto nu_l_outside = tp.nu_l(*(intersection.outside()),outside_cell_center_local,pg_outside);
721 auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
722 auto rho_l_outside = tp.rho_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
723 auto w_l = (pl_inside-pl_outside)/distance + aavg(rho_l_inside,rho_l_outside)*gn;
724 RF pc_upwind, s_l_upwind, s_g_upwind;
727 pc_upwind = pg_inside-pl_inside;
728 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
732 pc_upwind = pg_outside-pl_outside;
733 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
735 s_g_upwind = 1-s_l_upwind;
736 auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l_upwind)/
737 tp.mu_l(e,inside_cell_center_local,pl_inside);
738 auto lambda_l_outside = tp.kr_l(*(intersection.outside()),outside_cell_center_local,s_l_upwind)/
739 tp.mu_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
740 auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
741 auto nu_l = aavg(nu_l_inside,nu_l_outside);
744 vn[intersection.indexInInside()] = nu_l * sigma_l * w_l;
748 if (intersection.boundary())
751 auto d = geo_intersection.global(face_local);
752 d -= inside_cell_center_global;
753 auto distance = d.two_norm();
756 auto bc_l = tp.bc_l(intersection,face_local,time);
761 auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
762 auto g_l = tp.g_l(intersection,face_local,time);
763 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
764 auto w_l = (pl_inside-g_l)/distance + rho_l_inside*gn;
765 auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
766 auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l)/
767 tp.mu_l(e,inside_cell_center_local,pl_inside);
768 auto sigma_l = lambda_l_inside*k_abs_inside;
769 vn[intersection.indexInInside()] = nu_l_inside * sigma_l * w_l;
775 auto j_l = tp.j_l(intersection,face_local,time);
776 vn[intersection.indexInInside()] = j_l;
781 auto vstar = intersection.unitOuterNormal(face_local);
782 vstar *= vn[intersection.indexInInside()];
783 Dune::FieldVector<RF,dim> normalhat(0);
784 if (intersection.indexInInside()%2==0)
785 normalhat[intersection.indexInInside()/2] = -1.0;
787 normalhat[intersection.indexInInside()/2] = 1.0;
788 Dune::FieldVector<DF,dim> vstarhat(0);
789 B.umtv(vstar,vstarhat);
790 vstarhat *= determinant;
791 coeff[intersection.indexInInside()] = vstarhat*normalhat;
795 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
796 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
798 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
799 yhat.axpy(coeff[i],rt0vectors[i]);
815 return pl.getGridView();
821 T aavg (T a, T b)
const 827 T havg (T a, T b)
const 830 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
842 template<
typename TP,
typename PL,
typename PG>
845 typename PL::Traits::RangeFieldType,
846 PL::Traits::GridViewType::dimension,
847 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
851 using GV =
typename PL::Traits::GridViewType;
852 using DF =
typename GV::Grid::ctype;
853 using RF =
typename PL::Traits::RangeFieldType;
854 using RangeType =
typename PL::Traits::RangeType;
855 enum {
dim = PL::Traits::GridViewType::dimension };
856 using Element =
typename GV::Traits::template Codim<0>::Entity;
857 using IntersectionIterator =
typename GV::IntersectionIterator;
858 using Intersection =
typename IntersectionIterator::Intersection;
863 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
864 typename TP::Traits::RangeFieldType time;
867 using RT0RangeType =
typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
874 V_g (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
877 void set_time (
typename TP::Traits::RangeFieldType time_)
887 auto geo = e.geometry();
891 auto face_local = ref_el.position(0,0);
892 auto face_volume = geo.volume();
895 auto inside_cell_center_local = ref_el.position(0,0);
896 auto inside_cell_center_global = geo.global(inside_cell_center_local);
899 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
902 typename PL::Traits::RangeType pl_inside, pg_inside;
903 pl.evaluate(e,inside_cell_center_local,pl_inside);
904 pg.evaluate(e,inside_cell_center_local,pg_inside);
907 auto nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
912 auto B = geo.jacobianInverseTransposed(x);
913 auto determinant = B.determinant();
916 for (
const auto& intersection : intersections(pl.getGridView(),
e))
919 vn[intersection.indexInInside()] = 0.0;
922 auto geo_intersection = intersection.geometry();
928 if (intersection.neighbor())
931 auto geo_outside = intersection.outside().geometry();
933 auto outside_cell_center_local = ref_el_outside.position(0,0);
934 auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
937 auto d = outside_cell_center_global;
938 d -= inside_cell_center_global;
939 auto distance = d.two_norm();
942 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
945 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
948 typename PL::Traits::RangeType pl_outside, pg_outside;
949 pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
950 pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
953 RF pc_upwind, s_l_upwind, s_g_upwind;
956 auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
957 auto rho_g_outside = tp.rho_g(e,outside_cell_center_local,pg_outside);
958 auto w_g = (pg_inside-pg_outside)/distance + aavg(rho_g_inside,rho_g_outside)*gn;
961 pc_upwind = pg_inside-pl_inside;
962 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
966 pc_upwind = pg_outside-pl_outside;
967 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
969 s_g_upwind = 1-s_l_upwind;
970 auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g_upwind)/
971 tp.mu_g(e,inside_cell_center_local,pg_inside);
972 auto lambda_g_outside = tp.kr_g(*(intersection.outside()),outside_cell_center_local,s_g_upwind)/
973 tp.mu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
974 auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
976 auto nu_g_outside = tp.nu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
979 vn[intersection.indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
983 if (intersection.boundary())
986 auto d = geo_intersection.global(face_local);
987 d -= inside_cell_center_global;
988 auto distance = d.two_norm();
991 auto bc_g = tp.bc_g(intersection,face_local,time);
996 auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
997 auto g_g = tp.g_g(intersection,face_local,time);
998 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
999 auto w_g = (pg_inside-g_g)/distance + rho_g_inside*gn;
1000 auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
1002 auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g)/
1003 tp.mu_g(e,inside_cell_center_local,pg_inside);
1004 auto sigma_g = lambda_g_inside*k_abs_inside;
1006 vn[intersection.indexInInside()] = nu_g_inside * sigma_g * w_g;
1012 auto j_g = tp.j_g(intersection,face_local,time);
1013 vn[intersection.indexInInside()] = j_g; ;
1018 auto vstar = intersection.unitOuterNormal(face_local);
1019 vstar *= vn[intersection.indexInInside()];
1020 Dune::FieldVector<RF,dim> normalhat(0);
1021 if (intersection.indexInInside()%2==0)
1022 normalhat[intersection.indexInInside()/2] = -1.0;
1024 normalhat[intersection.indexInInside()/2] = 1.0;
1025 Dune::FieldVector<DF,dim> vstarhat(0);
1026 B.umtv(vstar,vstarhat);
1027 vstarhat *= determinant;
1028 coeff[intersection.indexInInside()] = vstarhat*normalhat;
1032 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1033 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1035 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1036 yhat.axpy(coeff[i],rt0vectors[i]);
1047 return pl.getGridView();
1052 template<
typename T>
1053 T aavg (T a, T b)
const 1058 template<
typename T>
1059 T havg (T a, T b)
const 1062 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
1070 #endif // DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH const IG & ig
Definition: constraints.hh:148
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:108
static const int dim
Definition: adaptivity.hh:83
const Entity & e
Definition: localfunctionspace.hh:111
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:843
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
leaf of a function tree
Definition: function.hh:298
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: twophaseccfv.hh:39
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:882
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:124
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:49
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:220
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
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:178
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:170
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:813
TwoPhaseOnePointTemporalOperator(TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
Definition: twophaseccfv.hh:559
IntersectionType
Enum describing the type of an intersection.
Definition: intersectiontype.hh:14
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:100
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:36
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:199
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:118
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:84
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:115
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:192
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:162
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:507
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:565
dimension of the domain
Definition: twophaseccfv.hh:32
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:154
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:877
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: twophaseccfv.hh:42
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:213
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r_s) const
Definition: twophaseccfv.hh:476
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:51
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:54
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:572
V_l(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:639
Definition: twophaseccfv.hh:245
T Traits
Definition: twophaseccfv.hh:73
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:138
V_g(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:874
R RangeType
range type
Definition: function.hh:61
sparsity pattern generator
Definition: pattern.hh:29
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:116
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
GV GridViewType
the grid view
Definition: twophaseccfv.hh:27
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:647
traits class for two phase parameter class
Definition: twophaseccfv.hh:24
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:642
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:77
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:45
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:132
sparsity pattern generator
Definition: pattern.hh:13
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:407
typename GV::Intersection IntersectionType
Definition: twophaseccfv.hh:55
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:146
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:281
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:228
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:608
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
base class for parameter class
Definition: twophaseccfv.hh:70
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: twophaseccfv.hh:48
Default flags for all local operators.
Definition: flags.hh:18
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:92
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:1045
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:206
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
Definition: twophaseccfv.hh:539
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
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:302
Definition: twophaseccfv.hh:59
TwoPhaseTwoPointFluxOperator(const TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
constructor: pass parameter object
Definition: twophaseccfv.hh:275
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:185
traits class holding the function signature, same as in local function
Definition: function.hh:176