2 #ifndef DUNE_PDELAB_TWOPHASEOP_HH
3 #define DUNE_PDELAB_TWOPHASEOP_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>
12 #include"../common/function.hh"
21 template<
typename GV,
typename RF>
37 typedef Dune::FieldVector<DomainFieldType,dimDomain>
DomainType;
46 typedef Dune::FieldVector<RF,GV::dimensionworld>
RangeType;
52 typedef typename GV::Traits::template Codim<0>::Entity
ElementType;
56 template<
typename GV,
typename RF>
63 typedef Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>
PermTensorType;
67 template<
class T,
class Imp>
74 typename Traits::RangeFieldType
75 phi (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
77 return asImp().phi(e,x);
81 typename Traits::RangeFieldType
82 pc (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
83 typename Traits::RangeFieldType
s_l)
const
85 return asImp().pc(e,x,s_l);
89 typename Traits::RangeFieldType
90 s_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
91 typename Traits::RangeFieldType
pc)
const
93 return asImp().s_l(e,x,pc);
97 typename Traits::RangeFieldType
98 kr_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
99 typename Traits::RangeFieldType
s_l)
const
101 return asImp().kr_l(e,x,s_l);
105 typename Traits::RangeFieldType
106 kr_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
107 typename Traits::RangeFieldType s_g)
const
109 return asImp().kr_g(e,x,s_g);
113 typename Traits::RangeFieldType
114 mu_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
115 typename Traits::RangeFieldType p_l)
const
117 return asImp().mu_l(e,x,p_l);
121 typename Traits::RangeFieldType
122 mu_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
123 typename Traits::RangeFieldType p_g)
const
125 return asImp().mu_l(e,x,p_g);
129 typename Traits::PermTensorType
130 k_abs (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x)
const
132 return asImp().k_abs(e,x);
136 const typename Traits::RangeType&
gravity ()
const
138 return asImp().gravity();
143 typename Traits::RangeFieldType
144 nu_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
145 typename Traits::RangeFieldType p_l)
const
147 return asImp().nu_l(e,x,p_l);
151 typename Traits::RangeFieldType
152 nu_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
153 typename Traits::RangeFieldType p_g)
const
155 return asImp().nu_g(e,x,p_g);
159 typename Traits::RangeFieldType
160 rho_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
161 typename Traits::RangeFieldType p_l)
const
163 return asImp().rho_l(e,x,p_l);
167 typename Traits::RangeFieldType
168 rho_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
169 typename Traits::RangeFieldType p_g)
const
171 return asImp().rho_g(e,x,p_g);
176 bc_l (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
178 return asImp().bc_l(is,x,time);
183 bc_g (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
185 return asImp().bc_g(is,x,time);
189 typename Traits::RangeFieldType
190 g_l (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
192 return asImp().g_l(is,x,time);
196 typename Traits::RangeFieldType
197 g_g (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
199 return asImp().g_g(is,x,time);
203 typename Traits::RangeFieldType
204 j_l (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
206 return asImp().j_l(is,x,time);
210 typename Traits::RangeFieldType
211 j_g (
const typename Traits::IntersectionType& is,
const typename Traits::IntersectionDomainType& x,
typename Traits::RangeFieldType time)
const
213 return asImp().j_g(is,x,time);
217 typename Traits::RangeFieldType
218 q_l (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
219 typename Traits::RangeFieldType time)
const
221 return asImp().q_l(e,x,time);
225 typename Traits::RangeFieldType
226 q_g (
const typename Traits::ElementType&
e,
const typename Traits::DomainType& x,
227 typename Traits::RangeFieldType time)
const
229 return asImp().q_g(e,x,time);
233 Imp& asImp () {
return static_cast<Imp &
> (*this);}
234 const Imp& asImp ()
const {
return static_cast<const Imp &
>(*this);}
242 template<
typename TP>
256 enum { dim = TP::Traits::GridViewType::dimension };
260 typedef typename TP::Traits::RangeFieldType Real;
274 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
278 template<
typename EG,
typename LFSV,
typename R>
282 typedef typename LFSV::template Child<liquid>::Type PLSpace;
285 typedef typename PLSpace::Traits::FiniteElementType::
286 Traits::LocalBasisType::Traits::DomainFieldType DF;
287 typedef typename PLSpace::Traits::FiniteElementType::
288 Traits::LocalBasisType::Traits::RangeFieldType RF;
291 const Dune::FieldVector<DF,dim>&
292 cell_center_local = Dune::ReferenceElements<DF,dim>::general(eg.geometry().type()).position(0,0);
293 RF cell_volume = eg.geometry().volume();
296 r.accumulate(lfsv, liquid, -scale_l * tp.q_l(eg.entity(),cell_center_local,time) * cell_volume);
297 r.accumulate(lfsv, gas, -scale_g * tp.q_g(eg.entity(),cell_center_local,time) * cell_volume);
302 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
304 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
305 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
306 R& r_s, R& r_n)
const
309 typedef typename LFSV::template Child<0>::Type PLSpace;
312 typedef typename PLSpace::Traits::FiniteElementType::
313 Traits::LocalBasisType::Traits::DomainFieldType DF;
314 typedef typename PLSpace::Traits::FiniteElementType::
315 Traits::LocalBasisType::Traits::RangeFieldType RF;
318 const Dune::FieldVector<DF,dim>&
319 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
320 const Dune::FieldVector<DF,dim>&
321 outside_cell_center_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
322 Dune::FieldVector<DF,IG::dimension>
323 inside_cell_center_global = ig.inside()->geometry().center();
324 Dune::FieldVector<DF,IG::dimension>
325 outside_cell_center_global = ig.outside()->geometry().center();
328 Dune::FieldVector<DF,dim> d(outside_cell_center_global);
329 d -= inside_cell_center_global;
330 RF distance = d.two_norm();
333 const Dune::FieldVector<DF,IG::dimension-1>&
334 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
335 RF face_volume = ig.geometry().volume();
338 RF k_abs_inside = tp.k_abs(*(ig.inside()),inside_cell_center_local);
339 RF k_abs_outside = tp.k_abs(*(ig.outside()),outside_cell_center_local);
342 RF gn = tp.gravity()*ig.unitOuterNormal(face_local);
345 RF rho_l_inside = tp.rho_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
346 RF rho_l_outside = tp.rho_l(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,liquid));
347 RF w_l = (x_s(lfsu_s,liquid)-x_n(lfsu_n,liquid))/distance + aavg(rho_l_inside,rho_l_outside)*gn;
348 RF pc_upwind, s_l_upwind, s_g_upwind;
349 RF nu_l = aavg(tp.nu_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid)),
350 tp.nu_l(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,liquid)));
353 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
354 s_l_upwind = tp.s_l(*(ig.inside()),inside_cell_center_local,pc_upwind);
358 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
359 s_l_upwind = tp.s_l(*(ig.outside()),outside_cell_center_local,pc_upwind);
361 s_g_upwind = 1-s_l_upwind;
362 RF lambda_l_inside = tp.kr_l(*(ig.inside()),inside_cell_center_local,s_l_upwind)/
363 tp.mu_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
364 RF lambda_l_outside = tp.kr_l(*(ig.outside()),outside_cell_center_local,s_l_upwind)/
365 tp.mu_l(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,liquid));
366 RF sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
368 r_s.accumulate(lfsv_s,liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
369 r_n.accumulate(lfsv_n,liquid, -scale_l * nu_l * sigma_l * w_l * face_volume);
372 RF rho_g_inside = tp.rho_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
373 RF rho_g_outside = tp.rho_g(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,gas));
374 RF w_g = (x_s(lfsu_s,gas)-x_n(lfsu_n,gas))/distance + aavg(rho_g_inside,rho_g_outside)*gn;
375 RF nu_g = aavg(tp.nu_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas)),
376 tp.nu_g(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,gas)));
381 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
382 s_l_upwind = tp.s_l(*(ig.inside()),inside_cell_center_local,pc_upwind);
386 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
387 s_l_upwind = tp.s_l(*(ig.outside()),outside_cell_center_local,pc_upwind);
389 s_g_upwind = 1-s_l_upwind;
391 RF lambda_g_inside = tp.kr_g(*(ig.inside()),inside_cell_center_local,s_g_upwind)/
392 tp.mu_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
393 RF lambda_g_outside = tp.kr_g(*(ig.outside()),outside_cell_center_local,s_g_upwind)/
394 tp.mu_g(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,gas));
395 RF sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
397 r_s.accumulate(lfsv_s, gas, scale_g * nu_g * sigma_g * w_g * face_volume);
398 r_n.accumulate(lfsv_n, gas, -scale_g * nu_g * sigma_g * w_g * face_volume);
403 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
405 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
409 typedef typename LFSV::template Child<0>::Type PLSpace;
412 typedef typename PLSpace::Traits::FiniteElementType::
413 Traits::LocalBasisType::Traits::DomainFieldType DF;
414 typedef typename PLSpace::Traits::FiniteElementType::
415 Traits::LocalBasisType::Traits::RangeFieldType RF;
418 const Dune::FieldVector<DF,dim-1>&
419 face_local = Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).position(0,0);
420 RF face_volume = ig.geometry().volume();
423 int bc_l = tp.bc_l(ig.intersection(),face_local,time);
424 int bc_g = tp.bc_g(ig.intersection(),face_local,time);
425 if (bc_l!=1 && bc_g!=1)
return;
428 const Dune::FieldVector<DF,dim>&
429 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
430 Dune::FieldVector<DF,dim>
431 inside_cell_center_global = ig.inside()->geometry().global(inside_cell_center_local);
434 Dune::FieldVector<DF,dim> d = ig.geometry().global(face_local);
435 d -= inside_cell_center_global;
436 RF distance = d.two_norm();
439 RF k_abs_inside = tp.k_abs(*(ig.inside()),inside_cell_center_local);
442 RF gn = tp.gravity()*ig.unitOuterNormal(face_local);
447 RF rho_l_inside = tp.rho_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
448 RF g_l = tp.g_l(ig.intersection(),face_local,time);
449 RF w_l = (x_s(lfsu_s,liquid)-g_l)/distance + rho_l_inside*gn;
450 RF s_l = tp.s_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
451 RF lambda_l_inside = tp.kr_l(*(ig.inside()),inside_cell_center_local,s_l)/
452 tp.mu_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
453 RF sigma_l = lambda_l_inside*k_abs_inside;
454 RF nu_l = tp.nu_l(*(ig.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 RF rho_g_inside = tp.rho_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
462 RF g_g = tp.g_g(ig.intersection(),face_local,time);
463 RF w_g = (x_s(lfsu_s,gas)-g_g)/distance + rho_g_inside*gn;
464 RF s_l = tp.s_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
466 RF lambda_g_inside = tp.kr_g(*(ig.inside()),inside_cell_center_local,s_g)/
467 tp.mu_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
468 RF sigma_g = lambda_g_inside*k_abs_inside;
469 RF nu_g = tp.nu_g(*(ig.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 typedef typename LFSV::template Child<0>::Type PLSpace;
482 typedef typename PLSpace::Traits::FiniteElementType::
483 Traits::LocalBasisType::Traits::DomainFieldType DF;
484 typedef typename PLSpace::Traits::FiniteElementType::
485 Traits::LocalBasisType::Traits::RangeFieldType RF;
488 const Dune::FieldVector<DF,dim-1>&
489 face_local = Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).position(0,0);
490 RF face_volume = ig.geometry().integrationElement(face_local)*
491 Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).volume();
494 int bc_l = tp.bc_l(ig.intersection(),face_local,time);
495 int bc_g = tp.bc_g(ig.intersection(),face_local,time);
496 if (bc_l!=0 && bc_g!=0)
return;
501 RF j_l = tp.j_l(ig.intersection(),face_local,time);
502 r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
508 RF j_g = tp.j_g(ig.intersection(),face_local,time);
509 r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
514 void setTime (
typename TP::Traits::RangeFieldType t)
521 typename TP::Traits::RangeFieldType time;
522 Real scale_l, scale_g;
525 T aavg (T a, T b)
const
531 T havg (T a, T b)
const
534 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
553 enum { dim = TP::Traits::GridViewType::dimension };
557 typedef typename TP::Traits::RangeFieldType Real;
567 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
572 void setTime (
typename TP::Traits::RangeFieldType t)
578 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
579 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
582 typedef typename LFSV::template Child<0>::Type PLSpace;
585 typedef typename PLSpace::Traits::FiniteElementType::
586 Traits::LocalBasisType::Traits::DomainFieldType DF;
587 typedef typename PLSpace::Traits::FiniteElementType::
588 Traits::LocalBasisType::Traits::RangeFieldType RF;
591 const Dune::FieldVector<DF,dim>&
592 cell_center_local = Dune::ReferenceElements<DF,dim>::general(eg.geometry().type()).position(0,0);
593 RF cell_volume = eg.geometry().volume();
595 RF phi = tp.phi(eg.entity(),cell_center_local);
596 RF s_l = tp.s_l(eg.entity(),cell_center_local,x(lfsu,gas)-x(lfsu,liquid));
598 r.accumulate(lfsv, liquid, scale_l * phi * s_l * tp.nu_l(eg.entity(),cell_center_local,x(lfsu,liquid)) * cell_volume);
599 r.accumulate(lfsv, gas, scale_g * phi * (1-s_l) * tp.nu_g(eg.entity(),cell_center_local,x(lfsu,gas)) * cell_volume);
604 typename TP::Traits::RangeFieldType time;
605 Real scale_l, scale_g;
617 template<
typename TP,
typename PL,
typename PG>
620 typename PL::Traits::RangeFieldType,
621 PL::Traits::GridViewType::dimension,
622 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
626 typedef typename PL::Traits::GridViewType GV;
627 typedef typename GV::Grid::ctype DF;
628 typedef typename PL::Traits::RangeFieldType RF;
629 typedef typename PL::Traits::RangeType RangeType;
630 enum { dim = PL::Traits::GridViewType::dimension };
631 typedef typename GV::Traits::template Codim<0>::Entity Element;
632 typedef typename GV::IntersectionIterator IntersectionIterator;
633 typedef typename IntersectionIterator::Intersection Intersection;
638 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
639 typename TP::Traits::RangeFieldType time;
642 typedef typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType RT0RangeType;
649 V_l (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
652 void set_time (
typename TP::Traits::RangeFieldType time_)
662 const Dune::FieldVector<DF,dim>&
663 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::
664 general(e.type()).position(0,0);
665 Dune::FieldVector<DF,dim>
666 inside_cell_center_global = e.geometry().global(inside_cell_center_local);
669 RF k_abs_inside = tp.k_abs(e,inside_cell_center_local);
672 typename PL::Traits::RangeType pl_inside, pg_inside;
673 pl.evaluate(e,inside_cell_center_local,pl_inside);
674 pg.evaluate(e,inside_cell_center_local,pg_inside);
677 RF nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
682 typename Traits::ElementType::Geometry::JacobianInverseTransposed
683 B = e.geometry().jacobianInverseTransposed(x);
684 RF determinant = B.determinant();
687 IntersectionIterator endit = pl.getGridView().iend(e);
688 for (IntersectionIterator iit = pl.getGridView().ibegin(e); iit!=endit; ++iit)
691 vn[iit->indexInInside()] = 0.0;
694 const Dune::FieldVector<DF,dim-1>&
695 face_local = Dune::ReferenceElements<DF,dim-1>::general(iit->geometry().type()).position(0,0);
701 const Dune::FieldVector<DF,dim>&
702 outside_cell_center_local = Dune::ReferenceElements<DF,dim>::
703 general(iit->outside()->type()).position(0,0);
704 Dune::FieldVector<DF,dim>
705 outside_cell_center_global = iit->outside()->geometry().global(outside_cell_center_local);
708 Dune::FieldVector<DF,dim> d(outside_cell_center_global);
709 d -= inside_cell_center_global;
710 RF distance = d.two_norm();
713 RF k_abs_outside = tp.k_abs(*(iit->outside()),outside_cell_center_local);
716 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
719 typename PL::Traits::RangeType pl_outside, pg_outside;
720 pl.evaluate(*(iit->outside()),outside_cell_center_local,pl_outside);
721 pg.evaluate(*(iit->outside()),outside_cell_center_local,pg_outside);
724 RF nu_l_outside = tp.nu_l(*(iit->outside()),outside_cell_center_local,pg_outside);
727 RF rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
728 RF rho_l_outside = tp.rho_l(*(iit->outside()),outside_cell_center_local,pl_outside);
729 RF w_l = (pl_inside-pl_outside)/distance + aavg(rho_l_inside,rho_l_outside)*gn;
730 RF pc_upwind, s_l_upwind, s_g_upwind;
733 pc_upwind = pg_inside-pl_inside;
734 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
738 pc_upwind = pg_outside-pl_outside;
739 s_l_upwind = tp.s_l(*(iit->outside()),outside_cell_center_local,pc_upwind);
741 s_g_upwind = 1-s_l_upwind;
742 RF lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l_upwind)/
743 tp.mu_l(e,inside_cell_center_local,pl_inside);
744 RF lambda_l_outside = tp.kr_l(*(iit->outside()),outside_cell_center_local,s_l_upwind)/
745 tp.mu_l(*(iit->outside()),outside_cell_center_local,pl_outside);
746 RF sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
747 RF nu_l = aavg(nu_l_inside,nu_l_outside);
750 vn[iit->indexInInside()] = nu_l * sigma_l * w_l;
757 Dune::FieldVector<DF,dim> d = iit->geometry().global(face_local);
758 d -= inside_cell_center_global;
759 RF distance = d.two_norm();
762 int bc_l = tp.bc_l(*iit,face_local,time);
767 RF rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
768 RF g_l = tp.g_l(*iit,face_local,time);
769 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
770 RF w_l = (pl_inside-g_l)/distance + rho_l_inside*gn;
771 RF s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
772 RF lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l)/
773 tp.mu_l(e,inside_cell_center_local,pl_inside);
774 RF sigma_l = lambda_l_inside*k_abs_inside;
775 vn[iit->indexInInside()] = nu_l_inside * sigma_l * w_l;
781 RF j_l = tp.j_l(*iit,face_local,time);
782 vn[iit->indexInInside()] = j_l;
787 Dune::FieldVector<DF,dim> vstar=iit->unitOuterNormal(face_local);
788 vstar *= vn[iit->indexInInside()];
789 Dune::FieldVector<RF,dim> normalhat(0);
790 if (iit->indexInInside()%2==0)
791 normalhat[iit->indexInInside()/2] = -1.0;
793 normalhat[iit->indexInInside()/2] = 1.0;
794 Dune::FieldVector<DF,dim> vstarhat(0);
795 B.umtv(vstar,vstarhat);
796 vstarhat *= determinant;
797 coeff[iit->indexInInside()] = vstarhat*normalhat;
801 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
802 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
804 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
805 yhat.axpy(coeff[i],rt0vectors[i]);
821 return pl.getGridView();
827 T aavg (T a, T b)
const
833 T havg (T a, T b)
const
836 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
848 template<
typename TP,
typename PL,
typename PG>
851 typename PL::Traits::RangeFieldType,
852 PL::Traits::GridViewType::dimension,
853 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
857 typedef typename PL::Traits::GridViewType GV;
858 typedef typename GV::Grid::ctype DF;
859 typedef typename PL::Traits::RangeFieldType RF;
860 typedef typename PL::Traits::RangeType RangeType;
861 enum { dim = PL::Traits::GridViewType::dimension };
862 typedef typename GV::Traits::template Codim<0>::Entity Element;
863 typedef typename GV::IntersectionIterator IntersectionIterator;
864 typedef typename IntersectionIterator::Intersection Intersection;
869 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
870 typename TP::Traits::RangeFieldType time;
873 typedef typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType RT0RangeType;
880 V_g (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
883 void set_time (
typename TP::Traits::RangeFieldType time_)
893 const Dune::FieldVector<DF,dim>&
894 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::
895 general(e.type()).position(0,0);
896 Dune::FieldVector<DF,dim>
897 inside_cell_center_global = e.geometry().global(inside_cell_center_local);
900 RF k_abs_inside = tp.k_abs(e,inside_cell_center_local);
903 typename PL::Traits::RangeType pl_inside, pg_inside;
904 pl.evaluate(e,inside_cell_center_local,pl_inside);
905 pg.evaluate(e,inside_cell_center_local,pg_inside);
908 RF nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
913 typename Traits::ElementType::Geometry::JacobianInverseTransposed
914 B = e.geometry().jacobianInverseTransposed(x);
915 RF determinant = B.determinant();
918 IntersectionIterator endit = pl.getGridView().iend(e);
919 for (IntersectionIterator iit = pl.getGridView().ibegin(e); iit!=endit; ++iit)
922 vn[iit->indexInInside()] = 0.0;
925 const Dune::FieldVector<DF,dim-1>&
926 face_local = Dune::ReferenceElements<DF,dim-1>::general(iit->geometry().type()).position(0,0);
931 const Dune::FieldVector<DF,dim>&
932 outside_cell_center_local = Dune::ReferenceElements<DF,dim>::
933 general(iit->outside()->type()).position(0,0);
934 Dune::FieldVector<DF,dim>
935 outside_cell_center_global = iit->outside()->geometry().global(outside_cell_center_local);
938 Dune::FieldVector<DF,dim> d(outside_cell_center_global);
939 d -= inside_cell_center_global;
940 RF distance = d.two_norm();
943 RF k_abs_outside = tp.k_abs(*(iit->outside()),outside_cell_center_local);
946 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
949 typename PL::Traits::RangeType pl_outside, pg_outside;
950 pl.evaluate(*(iit->outside()),outside_cell_center_local,pl_outside);
951 pg.evaluate(*(iit->outside()),outside_cell_center_local,pg_outside);
954 RF pc_upwind, s_l_upwind, s_g_upwind;
957 RF rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
958 RF rho_g_outside = tp.rho_g(e,outside_cell_center_local,pg_outside);
959 RF w_g = (pg_inside-pg_outside)/distance + aavg(rho_g_inside,rho_g_outside)*gn;
962 pc_upwind = pg_inside-pl_inside;
963 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
967 pc_upwind = pg_outside-pl_outside;
968 s_l_upwind = tp.s_l(*(iit->outside()),outside_cell_center_local,pc_upwind);
970 s_g_upwind = 1-s_l_upwind;
971 RF lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g_upwind)/
972 tp.mu_g(e,inside_cell_center_local,pg_inside);
973 RF lambda_g_outside = tp.kr_g(*(iit->outside()),outside_cell_center_local,s_g_upwind)/
974 tp.mu_g(*(iit->outside()),outside_cell_center_local,pg_outside);
975 RF sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
977 RF nu_g_outside = tp.nu_g(*(iit->outside()),outside_cell_center_local,pg_outside);
980 vn[iit->indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
987 Dune::FieldVector<DF,dim> d = iit->geometry().global(face_local);
988 d -= inside_cell_center_global;
989 RF distance = d.two_norm();
992 int bc_g = tp.bc_g(*iit,face_local,time);
997 RF rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
998 RF g_g = tp.g_g(*iit,face_local,time);
999 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
1000 RF w_g = (pg_inside-g_g)/distance + rho_g_inside*gn;
1001 RF s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
1003 RF lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g)/
1004 tp.mu_g(e,inside_cell_center_local,pg_inside);
1005 RF sigma_g = lambda_g_inside*k_abs_inside;
1007 vn[iit->indexInInside()] = nu_g_inside * sigma_g * w_g;
1013 RF j_g = tp.j_g(*iit,face_local,time);
1014 vn[iit->indexInInside()] = j_g; ;
1019 Dune::FieldVector<DF,dim> vstar=iit->unitOuterNormal(face_local);
1020 vstar *= vn[iit->indexInInside()];
1021 Dune::FieldVector<RF,dim> normalhat(0);
1022 if (iit->indexInInside()%2==0)
1023 normalhat[iit->indexInInside()/2] = -1.0;
1025 normalhat[iit->indexInInside()/2] = 1.0;
1026 Dune::FieldVector<DF,dim> vstarhat(0);
1027 B.umtv(vstar,vstarhat);
1028 vstarhat *= determinant;
1029 coeff[iit->indexInInside()] = vstarhat*normalhat;
1033 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1034 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1036 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1037 yhat.axpy(coeff[i],rt0vectors[i]);
1048 return pl.getGridView();
1053 template<
typename T>
1054 T aavg (T a, T b)
const
1059 template<
typename T>
1060 T havg (T a, T b)
const
1063 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:572
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:657
TwoPhaseParameterTraits< GV, RF > Base
Definition: twophaseccfv.hh:59
Definition: twophaseccfv.hh:264
traits class for two phase parameter class
Definition: twophaseccfv.hh:22
Definition: twophaseccfv.hh:564
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:218
Definition: twophaseccfv.hh:561
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:114
Definition: twophaseccfv.hh:546
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: twophaseccfv.hh:37
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:883
Dune::FieldMatrix< RangeFieldType, Base::dimDomain, Base::dimDomain > PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:63
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:226
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:176
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:1046
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:618
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:303
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:211
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:75
Dune::PDELab::GridFunctionTraits< GV, RF, dim, Dune::FieldVector< RF, dim > > Traits
Definition: twophaseccfv.hh:645
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
TwoPhaseOnePointTemporalOperator(TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
Definition: twophaseccfv.hh:566
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:168
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:82
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r_s) const
Definition: twophaseccfv.hh:476
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:114
TwoPhaseTwoPointFluxOperator(const TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
constructor: pass parameter object
Definition: twophaseccfv.hh:273
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:514
const IG & ig
Definition: constraints.hh:147
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:279
leaf of a function tree
Definition: function.hh:576
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:43
V_g(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:880
V_l(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:649
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: twophaseccfv.hh:57
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:204
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:106
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:183
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:48
Definition: twophaseccfv.hh:269
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:98
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:52
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:122
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:190
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:819
traits class holding the function signature, same as in local function
Definition: function.hh:175
GV GridViewType
the grid view
Definition: twophaseccfv.hh:25
T Traits
Definition: twophaseccfv.hh:71
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:152
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:90
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:49
Definition: adaptivity.hh:27
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:197
Definition: twophaseccfv.hh:267
Dune::PDELab::GridFunctionBase< Traits, V_l< TP, PL, PG > > BaseT
Definition: twophaseccfv.hh:647
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
Dune::PDELab::GridFunctionTraits< GV, RF, dim, Dune::FieldVector< RF, dim > > Traits
Definition: twophaseccfv.hh:876
GV::Intersection IntersectionType
Definition: twophaseccfv.hh:53
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:144
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: twophaseccfv.hh:40
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:117
Definition: twophaseccfv.hh:263
Base::RangeFieldType RangeFieldType
Definition: twophaseccfv.hh:60
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:404
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:652
sparsity pattern generator
Definition: pattern.hh:29
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:160
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:136
GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:34
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:130
R RangeType
range type
Definition: function.hh:60
Definition: twophaseccfv.hh:268
Definition: twophaseccfv.hh:243
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:579
base class for parameter class
Definition: twophaseccfv.hh:68
dimension of the domain
Definition: twophaseccfv.hh:30
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
Definition: twophaseccfv.hh:270
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:888
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: twophaseccfv.hh:46
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:849
Dune::PDELab::GridFunctionBase< Traits, V_g< TP, PL, PG > > BaseT
Definition: twophaseccfv.hh:878
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538