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/common/static_assert.hh>
8 #include<dune/geometry/referenceelements.hh>
9 #include<dune/localfunctions/raviartthomas/raviartthomascube.hh>
13 #include"../common/function.hh"
22 template<
typename GV,
typename RF>
38 typedef Dune::FieldVector<DomainFieldType,dimDomain>
DomainType;
47 typedef Dune::FieldVector<RF,GV::dimensionworld>
RangeType;
53 typedef typename GV::Traits::template Codim<0>::Entity
ElementType;
57 template<
typename GV,
typename RF>
64 typedef Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>
PermTensorType;
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>
257 enum { dim = TP::Traits::GridViewType::dimension };
261 typedef typename TP::Traits::RangeFieldType Real;
275 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
279 template<
typename EG,
typename LFSV,
typename R>
283 typedef typename LFSV::template Child<liquid>::Type PLSpace;
286 typedef typename PLSpace::Traits::FiniteElementType::
287 Traits::LocalBasisType::Traits::DomainFieldType DF;
288 typedef typename PLSpace::Traits::FiniteElementType::
289 Traits::LocalBasisType::Traits::RangeFieldType RF;
292 const Dune::FieldVector<DF,dim>&
293 cell_center_local = Dune::ReferenceElements<DF,dim>::general(eg.geometry().type()).position(0,0);
294 RF cell_volume = eg.geometry().volume();
297 r.accumulate(lfsv, liquid, -scale_l * tp.q_l(eg.entity(),cell_center_local,time) * cell_volume);
298 r.accumulate(lfsv, gas, -scale_g * tp.q_g(eg.entity(),cell_center_local,time) * cell_volume);
303 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
305 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
306 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
307 R& r_s, R& r_n)
const
310 typedef typename LFSV::template Child<0>::Type PLSpace;
313 typedef typename PLSpace::Traits::FiniteElementType::
314 Traits::LocalBasisType::Traits::DomainFieldType DF;
315 typedef typename PLSpace::Traits::FiniteElementType::
316 Traits::LocalBasisType::Traits::RangeFieldType RF;
319 const Dune::FieldVector<DF,dim>&
320 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
321 const Dune::FieldVector<DF,dim>&
322 outside_cell_center_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
323 Dune::FieldVector<DF,IG::dimension>
324 inside_cell_center_global = ig.inside()->geometry().center();
325 Dune::FieldVector<DF,IG::dimension>
326 outside_cell_center_global = ig.outside()->geometry().center();
329 Dune::FieldVector<DF,dim> d(outside_cell_center_global);
330 d -= inside_cell_center_global;
331 RF distance = d.two_norm();
334 const Dune::FieldVector<DF,IG::dimension-1>&
335 face_local = Dune::ReferenceElements<DF,IG::dimension-1>::general(ig.geometry().type()).position(0,0);
336 RF face_volume = ig.geometry().volume();
339 RF k_abs_inside = tp.k_abs(*(ig.inside()),inside_cell_center_local);
340 RF k_abs_outside = tp.k_abs(*(ig.outside()),outside_cell_center_local);
343 RF gn = tp.gravity()*ig.unitOuterNormal(face_local);
346 RF rho_l_inside = tp.rho_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
347 RF rho_l_outside = tp.rho_l(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,liquid));
348 RF w_l = (x_s(lfsu_s,liquid)-x_n(lfsu_n,liquid))/distance + aavg(rho_l_inside,rho_l_outside)*gn;
349 RF pc_upwind, s_l_upwind, s_g_upwind;
350 RF nu_l = aavg(tp.nu_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid)),
351 tp.nu_l(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,liquid)));
354 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
355 s_l_upwind = tp.s_l(*(ig.inside()),inside_cell_center_local,pc_upwind);
359 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
360 s_l_upwind = tp.s_l(*(ig.outside()),outside_cell_center_local,pc_upwind);
362 s_g_upwind = 1-s_l_upwind;
363 RF lambda_l_inside = tp.kr_l(*(ig.inside()),inside_cell_center_local,s_l_upwind)/
364 tp.mu_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
365 RF lambda_l_outside = tp.kr_l(*(ig.outside()),outside_cell_center_local,s_l_upwind)/
366 tp.mu_l(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,liquid));
367 RF sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
369 r_s.accumulate(lfsv_s,liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
370 r_n.accumulate(lfsv_n,liquid, -scale_l * nu_l * sigma_l * w_l * face_volume);
373 RF rho_g_inside = tp.rho_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
374 RF rho_g_outside = tp.rho_g(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,gas));
375 RF w_g = (x_s(lfsu_s,gas)-x_n(lfsu_n,gas))/distance + aavg(rho_g_inside,rho_g_outside)*gn;
376 RF nu_g = aavg(tp.nu_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas)),
377 tp.nu_g(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,gas)));
382 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
383 s_l_upwind = tp.s_l(*(ig.inside()),inside_cell_center_local,pc_upwind);
387 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
388 s_l_upwind = tp.s_l(*(ig.outside()),outside_cell_center_local,pc_upwind);
390 s_g_upwind = 1-s_l_upwind;
392 RF lambda_g_inside = tp.kr_g(*(ig.inside()),inside_cell_center_local,s_g_upwind)/
393 tp.mu_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
394 RF lambda_g_outside = tp.kr_g(*(ig.outside()),outside_cell_center_local,s_g_upwind)/
395 tp.mu_g(*(ig.outside()),outside_cell_center_local,x_n(lfsu_n,gas));
396 RF sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
398 r_s.accumulate(lfsv_s, gas, scale_g * nu_g * sigma_g * w_g * face_volume);
399 r_n.accumulate(lfsv_n, gas, -scale_g * nu_g * sigma_g * w_g * face_volume);
404 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
406 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
410 typedef typename LFSV::template Child<0>::Type PLSpace;
413 typedef typename PLSpace::Traits::FiniteElementType::
414 Traits::LocalBasisType::Traits::DomainFieldType DF;
415 typedef typename PLSpace::Traits::FiniteElementType::
416 Traits::LocalBasisType::Traits::RangeFieldType RF;
419 const Dune::FieldVector<DF,dim-1>&
420 face_local = Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).position(0,0);
421 RF face_volume = ig.geometry().volume();
424 int bc_l = tp.bc_l(ig.intersection(),face_local,time);
425 int bc_g = tp.bc_g(ig.intersection(),face_local,time);
426 if (bc_l!=1 && bc_g!=1)
return;
429 const Dune::FieldVector<DF,dim>&
430 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
431 Dune::FieldVector<DF,dim>
432 inside_cell_center_global = ig.inside()->geometry().global(inside_cell_center_local);
435 Dune::FieldVector<DF,dim> d = ig.geometry().global(face_local);
436 d -= inside_cell_center_global;
437 RF distance = d.two_norm();
440 RF k_abs_inside = tp.k_abs(*(ig.inside()),inside_cell_center_local);
443 RF gn = tp.gravity()*ig.unitOuterNormal(face_local);
448 RF rho_l_inside = tp.rho_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
449 RF g_l = tp.g_l(ig.intersection(),face_local,time);
450 RF w_l = (x_s(lfsu_s,liquid)-g_l)/distance + rho_l_inside*gn;
451 RF s_l = tp.s_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
452 RF lambda_l_inside = tp.kr_l(*(ig.inside()),inside_cell_center_local,s_l)/
453 tp.mu_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
454 RF sigma_l = lambda_l_inside*k_abs_inside;
455 RF nu_l = tp.nu_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,liquid));
456 r_s.accumulate(lfsv_s, liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
462 RF rho_g_inside = tp.rho_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
463 RF g_g = tp.g_g(ig.intersection(),face_local,time);
464 RF w_g = (x_s(lfsu_s,gas)-g_g)/distance + rho_g_inside*gn;
465 RF s_l = tp.s_l(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
467 RF lambda_g_inside = tp.kr_g(*(ig.inside()),inside_cell_center_local,s_g)/
468 tp.mu_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
469 RF sigma_g = lambda_g_inside*k_abs_inside;
470 RF nu_g = tp.nu_g(*(ig.inside()),inside_cell_center_local,x_s(lfsu_s,gas));
471 r_s.accumulate(lfsv_s, gas, scale_l * nu_g * sigma_g * w_g * face_volume);
476 template<
typename IG,
typename LFSV,
typename R>
480 typedef typename LFSV::template Child<0>::Type PLSpace;
483 typedef typename PLSpace::Traits::FiniteElementType::
484 Traits::LocalBasisType::Traits::DomainFieldType DF;
485 typedef typename PLSpace::Traits::FiniteElementType::
486 Traits::LocalBasisType::Traits::RangeFieldType RF;
489 const Dune::FieldVector<DF,dim-1>&
490 face_local = Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).position(0,0);
491 RF face_volume = ig.geometry().integrationElement(face_local)*
492 Dune::ReferenceElements<DF,dim-1>::general(ig.geometry().type()).volume();
495 int bc_l = tp.bc_l(ig.intersection(),face_local,time);
496 int bc_g = tp.bc_g(ig.intersection(),face_local,time);
497 if (bc_l!=0 && bc_g!=0)
return;
502 RF j_l = tp.j_l(ig.intersection(),face_local,time);
503 r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
509 RF j_g = tp.j_g(ig.intersection(),face_local,time);
510 r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
515 void setTime (
typename TP::Traits::RangeFieldType t)
522 typename TP::Traits::RangeFieldType time;
523 Real scale_l, scale_g;
526 T aavg (T a, T b)
const
532 T havg (T a, T b)
const
535 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
554 enum { dim = TP::Traits::GridViewType::dimension };
558 typedef typename TP::Traits::RangeFieldType Real;
568 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
573 void setTime (
typename TP::Traits::RangeFieldType t)
579 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
580 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
583 typedef typename LFSV::template Child<0>::Type PLSpace;
586 typedef typename PLSpace::Traits::FiniteElementType::
587 Traits::LocalBasisType::Traits::DomainFieldType DF;
588 typedef typename PLSpace::Traits::FiniteElementType::
589 Traits::LocalBasisType::Traits::RangeFieldType RF;
592 const Dune::FieldVector<DF,dim>&
593 cell_center_local = Dune::ReferenceElements<DF,dim>::general(eg.geometry().type()).position(0,0);
594 RF cell_volume = eg.geometry().volume();
596 RF phi = tp.phi(eg.entity(),cell_center_local);
597 RF s_l = tp.s_l(eg.entity(),cell_center_local,x(lfsu,gas)-x(lfsu,liquid));
599 r.accumulate(lfsv, liquid, scale_l * phi * s_l * tp.nu_l(eg.entity(),cell_center_local,x(lfsu,liquid)) * cell_volume);
600 r.accumulate(lfsv, gas, scale_g * phi * (1-s_l) * tp.nu_g(eg.entity(),cell_center_local,x(lfsu,gas)) * cell_volume);
605 typename TP::Traits::RangeFieldType time;
606 Real scale_l, scale_g;
618 template<
typename TP,
typename PL,
typename PG>
621 typename PL::Traits::RangeFieldType,
622 PL::Traits::GridViewType::dimension,
623 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
627 typedef typename PL::Traits::GridViewType GV;
628 typedef typename GV::Grid::ctype DF;
629 typedef typename PL::Traits::RangeFieldType RF;
630 typedef typename PL::Traits::RangeType RangeType;
631 enum { dim = PL::Traits::GridViewType::dimension };
632 typedef typename GV::Traits::template Codim<0>::Entity Element;
633 typedef typename GV::IntersectionIterator IntersectionIterator;
634 typedef typename IntersectionIterator::Intersection Intersection;
639 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
640 typename TP::Traits::RangeFieldType time;
643 typedef typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType RT0RangeType;
650 V_l (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
653 void set_time (
typename TP::Traits::RangeFieldType time_)
663 const Dune::FieldVector<DF,dim>&
664 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::
665 general(e.type()).position(0,0);
666 Dune::FieldVector<DF,dim>
667 inside_cell_center_global = e.geometry().global(inside_cell_center_local);
670 RF k_abs_inside = tp.k_abs(e,inside_cell_center_local);
673 typename PL::Traits::RangeType pl_inside, pg_inside;
674 pl.evaluate(e,inside_cell_center_local,pl_inside);
675 pg.evaluate(e,inside_cell_center_local,pg_inside);
678 RF nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
683 typename Traits::ElementType::Geometry::JacobianInverseTransposed
684 B = e.geometry().jacobianInverseTransposed(x);
685 RF determinant = B.determinant();
688 IntersectionIterator endit = pl.getGridView().iend(e);
689 for (IntersectionIterator iit = pl.getGridView().ibegin(e); iit!=endit; ++iit)
692 vn[iit->indexInInside()] = 0.0;
695 const Dune::FieldVector<DF,dim-1>&
696 face_local = Dune::ReferenceElements<DF,dim-1>::general(iit->geometry().type()).position(0,0);
702 const Dune::FieldVector<DF,dim>&
703 outside_cell_center_local = Dune::ReferenceElements<DF,dim>::
704 general(iit->outside()->type()).position(0,0);
705 Dune::FieldVector<DF,dim>
706 outside_cell_center_global = iit->outside()->geometry().global(outside_cell_center_local);
709 Dune::FieldVector<DF,dim> d(outside_cell_center_global);
710 d -= inside_cell_center_global;
711 RF distance = d.two_norm();
714 RF k_abs_outside = tp.k_abs(*(iit->outside()),outside_cell_center_local);
717 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
720 typename PL::Traits::RangeType pl_outside, pg_outside;
721 pl.evaluate(*(iit->outside()),outside_cell_center_local,pl_outside);
722 pg.evaluate(*(iit->outside()),outside_cell_center_local,pg_outside);
725 RF nu_l_outside = tp.nu_l(*(iit->outside()),outside_cell_center_local,pg_outside);
728 RF rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
729 RF rho_l_outside = tp.rho_l(*(iit->outside()),outside_cell_center_local,pl_outside);
730 RF w_l = (pl_inside-pl_outside)/distance + aavg(rho_l_inside,rho_l_outside)*gn;
731 RF pc_upwind, s_l_upwind, s_g_upwind;
734 pc_upwind = pg_inside-pl_inside;
735 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
739 pc_upwind = pg_outside-pl_outside;
740 s_l_upwind = tp.s_l(*(iit->outside()),outside_cell_center_local,pc_upwind);
742 s_g_upwind = 1-s_l_upwind;
743 RF lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l_upwind)/
744 tp.mu_l(e,inside_cell_center_local,pl_inside);
745 RF lambda_l_outside = tp.kr_l(*(iit->outside()),outside_cell_center_local,s_l_upwind)/
746 tp.mu_l(*(iit->outside()),outside_cell_center_local,pl_outside);
747 RF sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
748 RF nu_l = aavg(nu_l_inside,nu_l_outside);
751 vn[iit->indexInInside()] = nu_l * sigma_l * w_l;
758 Dune::FieldVector<DF,dim> d = iit->geometry().global(face_local);
759 d -= inside_cell_center_global;
760 RF distance = d.two_norm();
763 int bc_l = tp.bc_l(*iit,face_local,time);
768 RF rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
769 RF g_l = tp.g_l(*iit,face_local,time);
770 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
771 RF w_l = (pl_inside-g_l)/distance + rho_l_inside*gn;
772 RF s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
773 RF lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l)/
774 tp.mu_l(e,inside_cell_center_local,pl_inside);
775 RF sigma_l = lambda_l_inside*k_abs_inside;
776 vn[iit->indexInInside()] = nu_l_inside * sigma_l * w_l;
782 RF j_l = tp.j_l(*iit,face_local,time);
783 vn[iit->indexInInside()] = j_l;
788 Dune::FieldVector<DF,dim> vstar=iit->unitOuterNormal(face_local);
789 vstar *= vn[iit->indexInInside()];
790 Dune::FieldVector<RF,dim> normalhat(0);
791 if (iit->indexInInside()%2==0)
792 normalhat[iit->indexInInside()/2] = -1.0;
794 normalhat[iit->indexInInside()/2] = 1.0;
795 Dune::FieldVector<DF,dim> vstarhat(0);
796 B.umtv(vstar,vstarhat);
797 vstarhat *= determinant;
798 coeff[iit->indexInInside()] = vstarhat*normalhat;
802 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
803 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
805 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
806 yhat.axpy(coeff[i],rt0vectors[i]);
822 return pl.getGridView();
828 T aavg (T a, T b)
const
834 T havg (T a, T b)
const
837 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
849 template<
typename TP,
typename PL,
typename PG>
852 typename PL::Traits::RangeFieldType,
853 PL::Traits::GridViewType::dimension,
854 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
858 typedef typename PL::Traits::GridViewType GV;
859 typedef typename GV::Grid::ctype DF;
860 typedef typename PL::Traits::RangeFieldType RF;
861 typedef typename PL::Traits::RangeType RangeType;
862 enum { dim = PL::Traits::GridViewType::dimension };
863 typedef typename GV::Traits::template Codim<0>::Entity Element;
864 typedef typename GV::IntersectionIterator IntersectionIterator;
865 typedef typename IntersectionIterator::Intersection Intersection;
870 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
871 typename TP::Traits::RangeFieldType time;
874 typedef typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType RT0RangeType;
881 V_g (
const TP& tp_,
const PL& pl_,
const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
884 void set_time (
typename TP::Traits::RangeFieldType time_)
894 const Dune::FieldVector<DF,dim>&
895 inside_cell_center_local = Dune::ReferenceElements<DF,dim>::
896 general(e.type()).position(0,0);
897 Dune::FieldVector<DF,dim>
898 inside_cell_center_global = e.geometry().global(inside_cell_center_local);
901 RF k_abs_inside = tp.k_abs(e,inside_cell_center_local);
904 typename PL::Traits::RangeType pl_inside, pg_inside;
905 pl.evaluate(e,inside_cell_center_local,pl_inside);
906 pg.evaluate(e,inside_cell_center_local,pg_inside);
909 RF nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
914 typename Traits::ElementType::Geometry::JacobianInverseTransposed
915 B = e.geometry().jacobianInverseTransposed(x);
916 RF determinant = B.determinant();
919 IntersectionIterator endit = pl.getGridView().iend(e);
920 for (IntersectionIterator iit = pl.getGridView().ibegin(e); iit!=endit; ++iit)
923 vn[iit->indexInInside()] = 0.0;
926 const Dune::FieldVector<DF,dim-1>&
927 face_local = Dune::ReferenceElements<DF,dim-1>::general(iit->geometry().type()).position(0,0);
932 const Dune::FieldVector<DF,dim>&
933 outside_cell_center_local = Dune::ReferenceElements<DF,dim>::
934 general(iit->outside()->type()).position(0,0);
935 Dune::FieldVector<DF,dim>
936 outside_cell_center_global = iit->outside()->geometry().global(outside_cell_center_local);
939 Dune::FieldVector<DF,dim> d(outside_cell_center_global);
940 d -= inside_cell_center_global;
941 RF distance = d.two_norm();
944 RF k_abs_outside = tp.k_abs(*(iit->outside()),outside_cell_center_local);
947 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
950 typename PL::Traits::RangeType pl_outside, pg_outside;
951 pl.evaluate(*(iit->outside()),outside_cell_center_local,pl_outside);
952 pg.evaluate(*(iit->outside()),outside_cell_center_local,pg_outside);
955 RF pc_upwind, s_l_upwind, s_g_upwind;
958 RF rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
959 RF rho_g_outside = tp.rho_g(e,outside_cell_center_local,pg_outside);
960 RF w_g = (pg_inside-pg_outside)/distance + aavg(rho_g_inside,rho_g_outside)*gn;
963 pc_upwind = pg_inside-pl_inside;
964 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
968 pc_upwind = pg_outside-pl_outside;
969 s_l_upwind = tp.s_l(*(iit->outside()),outside_cell_center_local,pc_upwind);
971 s_g_upwind = 1-s_l_upwind;
972 RF lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g_upwind)/
973 tp.mu_g(e,inside_cell_center_local,pg_inside);
974 RF lambda_g_outside = tp.kr_g(*(iit->outside()),outside_cell_center_local,s_g_upwind)/
975 tp.mu_g(*(iit->outside()),outside_cell_center_local,pg_outside);
976 RF sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
978 RF nu_g_outside = tp.nu_g(*(iit->outside()),outside_cell_center_local,pg_outside);
981 vn[iit->indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
988 Dune::FieldVector<DF,dim> d = iit->geometry().global(face_local);
989 d -= inside_cell_center_global;
990 RF distance = d.two_norm();
993 int bc_g = tp.bc_g(*iit,face_local,time);
998 RF rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
999 RF g_g = tp.g_g(*iit,face_local,time);
1000 RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
1001 RF w_g = (pg_inside-g_g)/distance + rho_g_inside*gn;
1002 RF s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
1004 RF lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g)/
1005 tp.mu_g(e,inside_cell_center_local,pg_inside);
1006 RF sigma_g = lambda_g_inside*k_abs_inside;
1008 vn[iit->indexInInside()] = nu_g_inside * sigma_g * w_g;
1014 RF j_g = tp.j_g(*iit,face_local,time);
1015 vn[iit->indexInInside()] = j_g; ;
1020 Dune::FieldVector<DF,dim> vstar=iit->unitOuterNormal(face_local);
1021 vstar *= vn[iit->indexInInside()];
1022 Dune::FieldVector<RF,dim> normalhat(0);
1023 if (iit->indexInInside()%2==0)
1024 normalhat[iit->indexInInside()/2] = -1.0;
1026 normalhat[iit->indexInInside()/2] = 1.0;
1027 Dune::FieldVector<DF,dim> vstarhat(0);
1028 B.umtv(vstar,vstarhat);
1029 vstarhat *= determinant;
1030 coeff[iit->indexInInside()] = vstarhat*normalhat;
1034 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1035 rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1037 for (
unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1038 yhat.axpy(coeff[i],rt0vectors[i]);
1049 return pl.getGridView();
1054 template<
typename T>
1055 T aavg (T a, T b)
const
1060 template<
typename T>
1061 T havg (T a, T b)
const
1064 return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
traits class for two phase parameter class
Definition: twophaseccfv.hh:23
Dune::FieldMatrix< RangeFieldType, Base::dimDomain, Base::dimDomain > PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:64
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:580
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:137
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:76
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:53
sparsity pattern generator
Definition: pattern.hh:30
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
TwoPhaseTwoPointFluxOperator(const TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
constructor: pass parameter object
Definition: twophaseccfv.hh:274
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:658
Dune::PDELab::GridFunctionTraits< GV, RF, dim, Dune::FieldVector< RF, dim > > Traits
Definition: twophaseccfv.hh:646
Definition: twophaseccfv.hh:264
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
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: twophaseccfv.hh:47
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:118
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
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:405
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:653
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:1047
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 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
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:619
V_l(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:650
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
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
Definition: twophaseccfv.hh:271
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:884
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
dimension of the domain
Definition: twophaseccfv.hh:31
base class for parameter class
Definition: twophaseccfv.hh:69
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:49
Definition: twophaseccfv.hh:565
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:304
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
TwoPhaseParameterTraits< GV, RF > Base
Definition: twophaseccfv.hh:60
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r_s) const
Definition: twophaseccfv.hh:477
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:820
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:573
Dune::PDELab::GridFunctionBase< Traits, V_g< TP, PL, PG > > BaseT
Definition: twophaseccfv.hh:879
R RangeType
range type
Definition: function.hh:61
const IG & ig
Definition: common/constraints.hh:146
leaf of a function tree
Definition: function.hh:577
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
Definition: twophaseccfv.hh:265
GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:35
Definition: twophaseccfv.hh:547
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: twophaseccfv.hh:41
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 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
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:850
Dune::PDELab::GridFunctionBase< Traits, V_l< TP, PL, PG > > BaseT
Definition: twophaseccfv.hh:648
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:44
V_g(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:881
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 kr_l(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType s_l) const
liquid phase relative permeability
Definition: twophaseccfv.hh:99
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
GV::Intersection IntersectionType
Definition: twophaseccfv.hh:54
Definition: twophaseccfv.hh:269
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:131
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::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 class holding the function signature, same as in local function
Definition: function.hh:176
Base::RangeFieldType RangeFieldType
Definition: twophaseccfv.hh:61
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: twophaseccfv.hh:270
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:115
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: twophaseccfv.hh:38
const EG & eg
Definition: common/constraints.hh:277
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:889
Definition: twophaseccfv.hh:244
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
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
T Traits
Definition: twophaseccfv.hh:72
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
Definition: twophaseccfv.hh:562
Dune::PDELab::GridFunctionTraits< GV, RF, dim, Dune::FieldVector< RF, dim > > Traits
Definition: twophaseccfv.hh:877
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:515
const E & e
Definition: interpolate.hh:172
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
TwoPhaseOnePointTemporalOperator(TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
Definition: twophaseccfv.hh:567
Definition: twophaseccfv.hh:58
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
Definition: twophaseccfv.hh:268
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
GV GridViewType
the grid view
Definition: twophaseccfv.hh:26
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:50