2#ifndef DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
3#define DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
7#include<dune/geometry/referenceelements.hh>
10#include<dune/pdelab/common/quadraturerules.hh>
11#include <dune/pdelab/localoperator/defaultimp.hh>
13#include"../common/function.hh"
22 template<
typename GV,
typename RF>
54 using IntersectionType =
typename GV::Intersection;
57 template<
typename GV,
typename RF>
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);
177 bc_l (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
179 return asImp().bc_l(is,x,time);
184 bc_g (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
186 return asImp().bc_g(is,x,time);
190 typename Traits::RangeFieldType
191 g_l (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
193 return asImp().g_l(is,x,time);
197 typename Traits::RangeFieldType
198 g_g (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
200 return asImp().g_g(is,x,time);
204 typename Traits::RangeFieldType
205 j_l (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
207 return asImp().j_l(is,x,time);
211 typename Traits::RangeFieldType
212 j_g (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
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>
244 class TwoPhaseTwoPointFluxOperator
245 :
public NumericalJacobianSkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
246 public NumericalJacobianApplySkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
248 public NumericalJacobianBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
249 public NumericalJacobianApplyBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
251 public FullSkeletonPattern,
252 public FullVolumePattern,
253 public LocalOperatorDefaultFlags,
255 public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
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 };
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_)
279 template<
typename EG,
typename LFSV,
typename R>
280 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
283 const auto& cell = eg.entity();
286 auto geo = eg.geometry();
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>
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
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();
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();
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>
406 void alpha_boundary (
const IG& ig,
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();
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;
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>
475 void lambda_boundary (
const IG& ig,
const LFSV& lfsv, R& r_s)
const
478 auto geo = ig.geometry();
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();
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 };
621 using IntersectionIterator =
typename GV::IntersectionIterator;
628 typename TP::Traits::RangeFieldType time;
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();
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();
692 if (intersection.neighbor())
695 auto geo_outside = intersection.outside().geometry();
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()];
783 if (intersection.indexInInside()%2==0)
784 normalhat[intersection.indexInInside()/2] = -1.0;
786 normalhat[intersection.indexInInside()/2] = 1.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 };
856 using IntersectionIterator =
typename GV::IntersectionIterator;
863 typename TP::Traits::RangeFieldType time;
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();
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();
927 if (intersection.neighbor())
930 auto geo_outside = intersection.outside().geometry();
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()];
1020 if (intersection.indexInInside()%2==0)
1021 normalhat[intersection.indexInInside()/2] = -1.0;
1023 normalhat[intersection.indexInInside()/2] = 1.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));
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Dune::Intersection< GridImp, IntersectionImp > Intersection
Type of Intersection this IntersectionIterator points to.
Definition: intersectioniterator.hh:111
sparsity pattern generator
Definition: pattern.hh:14
leaf of a function tree
Definition: function.hh:302
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Default flags for all local operators.
Definition: flags.hh:19
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
Definition: twophaseccfv.hh:544
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:564
base class for parameter class
Definition: twophaseccfv.hh:70
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
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
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::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
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
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
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:131
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:137
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
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:76
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_g(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType p_g) const
gas phase viscosity
Definition: twophaseccfv.hh:123
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 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
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
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
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
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
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
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
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
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:848
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:613
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune namespace.
Definition: alignedallocator.hh:13
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
Convenience header that includes all available Raviart-Thomas local finite elements for cubes.
Static tag representing a codimension.
Definition: dimension.hh:24
R RangeType
range type
Definition: function.hh:62
traits class holding the function signature, same as in local function
Definition: function.hh:183
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:119
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:116
traits class for two phase parameter class
Definition: twophaseccfv.hh:24
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:44
GV GridViewType
the grid view
Definition: twophaseccfv.hh:26
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:35
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:50
@ dimDomain
dimension of the domain
Definition: twophaseccfv.hh:31
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:53