dune-pdelab  2.4-dev
twophaseccfv.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_TWOPHASEOP_HH
3 #define DUNE_PDELAB_TWOPHASEOP_HH
4 
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>
9 
11 
12 #include"../common/function.hh"
13 #include"pattern.hh"
14 #include"flags.hh"
15 #include"idefault.hh"
16 
17 namespace Dune {
18  namespace PDELab {
19 
21  template<typename GV, typename RF>
23  {
25  typedef GV GridViewType;
26 
28  enum {
30  dimDomain = GV::dimension
31  };
32 
34  typedef typename GV::Grid::ctype DomainFieldType;
35 
37  typedef Dune::FieldVector<DomainFieldType,dimDomain> DomainType;
38 
40  typedef Dune::FieldVector<DomainFieldType,dimDomain-1> IntersectionDomainType;
41 
43  typedef RF RangeFieldType;
44 
46  typedef Dune::FieldVector<RF,GV::dimensionworld> RangeType;
47 
49  typedef RangeFieldType PermTensorType;
50 
52  typedef typename GV::Traits::template Codim<0>::Entity ElementType;
53  typedef typename GV::Intersection IntersectionType;
54  };
55 
56  template<typename GV, typename RF>
58  {
61 
63  typedef Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain> PermTensorType;
64  };
65 
67  template<class T, class Imp>
69  {
70  public:
71  typedef T Traits;
72 
74  typename Traits::RangeFieldType
75  phi (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
76  {
77  return asImp().phi(e,x);
78  }
79 
81  typename Traits::RangeFieldType
82  pc (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
83  typename Traits::RangeFieldType s_l) const
84  {
85  return asImp().pc(e,x,s_l);
86  }
87 
89  typename Traits::RangeFieldType
90  s_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
91  typename Traits::RangeFieldType pc) const
92  {
93  return asImp().s_l(e,x,pc);
94  }
95 
97  typename Traits::RangeFieldType
98  kr_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
99  typename Traits::RangeFieldType s_l) const
100  {
101  return asImp().kr_l(e,x,s_l);
102  }
103 
105  typename Traits::RangeFieldType
106  kr_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
107  typename Traits::RangeFieldType s_g) const
108  {
109  return asImp().kr_g(e,x,s_g);
110  }
111 
113  typename Traits::RangeFieldType
114  mu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
115  typename Traits::RangeFieldType p_l) const
116  {
117  return asImp().mu_l(e,x,p_l);
118  }
119 
121  typename Traits::RangeFieldType
122  mu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
123  typename Traits::RangeFieldType p_g) const
124  {
125  return asImp().mu_l(e,x,p_g);
126  }
127 
129  typename Traits::PermTensorType
130  k_abs (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
131  {
132  return asImp().k_abs(e,x);
133  }
134 
136  const typename Traits::RangeType& gravity () const
137  {
138  return asImp().gravity();
139  }
140 
142  template<typename E>
143  typename Traits::RangeFieldType
144  nu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
145  typename Traits::RangeFieldType p_l) const
146  {
147  return asImp().nu_l(e,x,p_l);
148  }
149 
151  typename Traits::RangeFieldType
152  nu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
153  typename Traits::RangeFieldType p_g) const
154  {
155  return asImp().nu_g(e,x,p_g);
156  }
157 
159  typename Traits::RangeFieldType
160  rho_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
161  typename Traits::RangeFieldType p_l) const
162  {
163  return asImp().rho_l(e,x,p_l);
164  }
165 
167  typename Traits::RangeFieldType
168  rho_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
169  typename Traits::RangeFieldType p_g) const
170  {
171  return asImp().rho_g(e,x,p_g);
172  }
173 
175  int
176  bc_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
177  {
178  return asImp().bc_l(is,x,time);
179  }
180 
182  int
183  bc_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
184  {
185  return asImp().bc_g(is,x,time);
186  }
187 
189  typename Traits::RangeFieldType
190  g_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
191  {
192  return asImp().g_l(is,x,time);
193  }
194 
196  typename Traits::RangeFieldType
197  g_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
198  {
199  return asImp().g_g(is,x,time);
200  }
201 
203  typename Traits::RangeFieldType
204  j_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
205  {
206  return asImp().j_l(is,x,time);
207  }
208 
210  typename Traits::RangeFieldType
211  j_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
212  {
213  return asImp().j_g(is,x,time);
214  }
215 
217  typename Traits::RangeFieldType
218  q_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
219  typename Traits::RangeFieldType time) const
220  {
221  return asImp().q_l(e,x,time);
222  }
223 
225  typename Traits::RangeFieldType
226  q_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
227  typename Traits::RangeFieldType time) const
228  {
229  return asImp().q_g(e,x,time);
230  }
231 
232  private:
233  Imp& asImp () {return static_cast<Imp &> (*this);}
234  const Imp& asImp () const {return static_cast<const Imp &>(*this);}
235  };
236 
237 
238  // a local operator for solving the two-phase flow in pressure-pressure formulation
239  // with two-point flux approximation
240  // TP : parameter class, see above
241  // V : Vector holding last time step
242  template<typename TP>
244  : public NumericalJacobianSkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
245  public NumericalJacobianApplySkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
246 
247  public NumericalJacobianBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
248  public NumericalJacobianApplyBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
249 
250  public FullSkeletonPattern,
251  public FullVolumePattern,
253 
254  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
255  {
256  enum { dim = TP::Traits::GridViewType::dimension };
257  enum { liquid = 0 };
258  enum { gas = 1 };
259 
260  typedef typename TP::Traits::RangeFieldType Real;
261  public:
262  // pattern assembly flags
263  enum { doPatternVolume = true };
264  enum { doPatternSkeleton = true };
265 
266  // residual assembly flags
267  enum { doAlphaSkeleton = true };
268  enum { doAlphaBoundary = true };
269  enum { doLambdaVolume = true };
270  enum { doLambdaBoundary = true };
271 
273  TwoPhaseTwoPointFluxOperator (const TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
274  : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
275  {}
276 
277  // volume integral depending only on test functions
278  template<typename EG, typename LFSV, typename R>
279  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
280  {
281  // select the two components
282  typedef typename LFSV::template Child<liquid>::Type PLSpace;
283 
284  // domain and range field type
285  typedef typename PLSpace::Traits::FiniteElementType::
286  Traits::LocalBasisType::Traits::DomainFieldType DF;
287  typedef typename PLSpace::Traits::FiniteElementType::
288  Traits::LocalBasisType::Traits::RangeFieldType RF;
289 
290  // cell geometry
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();
294 
295  // contribution from source term
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);
298  }
299 
300  // skeleton integral depending on test and ansatz functions
301  // each face is only visited ONCE!
302  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
303  void alpha_skeleton (const IG& ig,
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
307  {
308  // select the two components
309  typedef typename LFSV::template Child<0>::Type PLSpace;
310 
311  // domain and range field type
312  typedef typename PLSpace::Traits::FiniteElementType::
313  Traits::LocalBasisType::Traits::DomainFieldType DF;
314  typedef typename PLSpace::Traits::FiniteElementType::
315  Traits::LocalBasisType::Traits::RangeFieldType RF;
316 
317  // cell geometries
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();
326 
327  // distance of cell centers
328  Dune::FieldVector<DF,dim> d(outside_cell_center_global);
329  d -= inside_cell_center_global;
330  RF distance = d.two_norm();
331 
332  // face geometry
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();
336 
337  // absolute permeability
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);
340 
341  // gravity times normal
342  RF gn = tp.gravity()*ig.unitOuterNormal(face_local);
343 
344  // liquid phase calculation
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; // determines direction
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)));
351  if (w_l>=0) // upwind capillary pressure on face
352  {
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);
355  }
356  else
357  {
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);
360  }
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);
367 
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);
370 
371  // gas phase calculation
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; // determines direction
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)));
377  if (w_l*w_g<=0) // new evaluation necessary only if signs differ
378  {
379  if (w_g>=0) // upwind capillary pressure on face
380  {
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);
383  }
384  else
385  {
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);
388  }
389  s_g_upwind = 1-s_l_upwind;
390  }
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);
396 
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);
399  }
400 
401  // skeleton integral depending on test and ansatz functions
402  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
403  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
404  void alpha_boundary (const IG& ig,
405  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
406  R& r_s) const
407  {
408  // select the two components
409  typedef typename LFSV::template Child<0>::Type PLSpace;
410 
411  // domain and range field type
412  typedef typename PLSpace::Traits::FiniteElementType::
413  Traits::LocalBasisType::Traits::DomainFieldType DF;
414  typedef typename PLSpace::Traits::FiniteElementType::
415  Traits::LocalBasisType::Traits::RangeFieldType RF;
416 
417  // face geometry
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();
421 
422  // evaluate boundary condition type
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; // no Dirichlet boundary conditions
426 
427  // cell geometry
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);
432 
433  // distance of cell center to boundary
434  Dune::FieldVector<DF,dim> d = ig.geometry().global(face_local);
435  d -= inside_cell_center_global;
436  RF distance = d.two_norm();
437 
438  // absolute permeability
439  RF k_abs_inside = tp.k_abs(*(ig.inside()),inside_cell_center_local);
440 
441  // gravity times normal
442  RF gn = tp.gravity()*ig.unitOuterNormal(face_local);
443 
444  // liquid phase Dirichlet boundary
445  if (bc_l==1)
446  {
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);
456  }
457 
458  // gas phase Dirichlet boundary
459  if (bc_g==1)
460  {
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));
465  RF s_g = 1-s_l;
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);
471  }
472  }
473 
474  // boundary integral independent of ansatz functions
475  template<typename IG, typename LFSV, typename R>
476  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r_s) const
477  {
478  // select the two components
479  typedef typename LFSV::template Child<0>::Type PLSpace;
480 
481  // domain and range field type
482  typedef typename PLSpace::Traits::FiniteElementType::
483  Traits::LocalBasisType::Traits::DomainFieldType DF;
484  typedef typename PLSpace::Traits::FiniteElementType::
485  Traits::LocalBasisType::Traits::RangeFieldType RF;
486 
487  // face geometry
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();
492 
493  // evaluate boundary condition type
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; // no Neumann boundary conditions
497 
498  // liquid phase Neumann boundary
499  if (bc_l==0)
500  {
501  RF j_l = tp.j_l(ig.intersection(),face_local,time);
502  r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
503  }
504 
505  // gas phase Neumann boundary
506  if (bc_g==0)
507  {
508  RF j_g = tp.j_g(ig.intersection(),face_local,time);
509  r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
510  }
511  }
512 
514  void setTime (typename TP::Traits::RangeFieldType t)
515  {
516  time = t;
517  }
518 
519  private:
520  const TP& tp; // two phase parameter class
521  typename TP::Traits::RangeFieldType time;
522  Real scale_l, scale_g;
523 
524  template<typename T>
525  T aavg (T a, T b) const
526  {
527  return 0.5*(a+b);
528  }
529 
530  template<typename T>
531  T havg (T a, T b) const
532  {
533  T eps = 1E-30;
534  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
535  }
536  };
537 
538 
545  template<class TP>
547  : public NumericalJacobianVolume<TwoPhaseOnePointTemporalOperator<TP> >,
548  public NumericalJacobianApplyVolume<TwoPhaseOnePointTemporalOperator<TP> >,
549  public FullVolumePattern,
551  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
552  {
553  enum { dim = TP::Traits::GridViewType::dimension };
554  enum { liquid = 0 };
555  enum { gas = 1 };
556 
557  typedef typename TP::Traits::RangeFieldType Real;
558 
559  public:
560  // pattern assembly flags
561  enum { doPatternVolume = true };
562 
563  // residual assembly flags
564  enum { doAlphaVolume = true };
565 
566  TwoPhaseOnePointTemporalOperator (TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
567  : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
568  {
569  }
570 
572  void setTime (typename TP::Traits::RangeFieldType t)
573  {
574  time = t;
575  }
576 
577  // volume integral depending on test and ansatz functions
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
580  {
581  // select the two components
582  typedef typename LFSV::template Child<0>::Type PLSpace;
583 
584  // domain and range field type
585  typedef typename PLSpace::Traits::FiniteElementType::
586  Traits::LocalBasisType::Traits::DomainFieldType DF;
587  typedef typename PLSpace::Traits::FiniteElementType::
588  Traits::LocalBasisType::Traits::RangeFieldType RF;
589 
590  // cell geometry
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();
594 
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));
597 
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);
600  }
601 
602  private:
603  TP& tp;
604  typename TP::Traits::RangeFieldType time;
605  Real scale_l, scale_g;
606  };
607 
608 
617  template<typename TP, typename PL, typename PG>
618  class V_l
619  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
620  typename PL::Traits::RangeFieldType,
621  PL::Traits::GridViewType::dimension,
622  Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
623  V_l<TP,PL,PG> >
624  {
625  // extract useful types
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;
634 
635  const TP& tp;
636  const PL& pl;
637  const PG& pg;
638  Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
639  typename TP::Traits::RangeFieldType time;
640 
641 
642  typedef typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType RT0RangeType;
643 
644  public:
646 
648 
649  V_l (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
650 
651  // set time where operator is to be evaluated (i.e. end of the time intervall)
652  void set_time (typename TP::Traits::RangeFieldType time_)
653  {
654  time = time_;
655  }
656 
657  inline void evaluate (const typename Traits::ElementType& e,
658  const typename Traits::DomainType& x,
659  typename Traits::RangeType& y) const
660  {
661  // cell geometry
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);
667 
668  // absolute permeability
669  RF k_abs_inside = tp.k_abs(e,inside_cell_center_local);
670 
671  // pressure evaluation
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);
675 
676  // density
677  RF nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
678 
679  // for coefficient computation
680  RF vn[2*dim]; // normal velocities
681  RF coeff[2*dim]; // RT0 coefficient
682  typename Traits::ElementType::Geometry::JacobianInverseTransposed
683  B = e.geometry().jacobianInverseTransposed(x); // the transformation. Assume it is linear
684  RF determinant = B.determinant();
685 
686  // loop over cell neighbors
687  IntersectionIterator endit = pl.getGridView().iend(e);
688  for (IntersectionIterator iit = pl.getGridView().ibegin(e); iit!=endit; ++iit)
689  {
690  // set to zero for processor boundary
691  vn[iit->indexInInside()] = 0.0;
692 
693  // face geometry
694  const Dune::FieldVector<DF,dim-1>&
695  face_local = Dune::ReferenceElements<DF,dim-1>::general(iit->geometry().type()).position(0,0);
696 
697 
698  // interior face
699  if (iit->neighbor())
700  {
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);
706 
707  // distance of cell centers
708  Dune::FieldVector<DF,dim> d(outside_cell_center_global);
709  d -= inside_cell_center_global;
710  RF distance = d.two_norm();
711 
712  // absolute permeability
713  RF k_abs_outside = tp.k_abs(*(iit->outside()),outside_cell_center_local);
714 
715  // gravity times normal
716  RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
717 
718  // pressure evaluation
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);
722 
723  // density
724  RF nu_l_outside = tp.nu_l(*(iit->outside()),outside_cell_center_local,pg_outside);
725 
726  // liquid phase calculation
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; // determines direction
730  RF pc_upwind, s_l_upwind, s_g_upwind;
731  if (w_l>=0) // upwind capillary pressure on face
732  {
733  pc_upwind = pg_inside-pl_inside;
734  s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
735  }
736  else
737  {
738  pc_upwind = pg_outside-pl_outside;
739  s_l_upwind = tp.s_l(*(iit->outside()),outside_cell_center_local,pc_upwind);
740  }
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);
748 
749  // set coefficient
750  vn[iit->indexInInside()] = nu_l * sigma_l * w_l;
751  }
752 
753  // boundary face
754  if (iit->boundary())
755  {
756  // distance of cell center to boundary
757  Dune::FieldVector<DF,dim> d = iit->geometry().global(face_local);
758  d -= inside_cell_center_global;
759  RF distance = d.two_norm();
760 
761  // evaluate boundary condition type
762  int bc_l = tp.bc_l(*iit,face_local,time);
763 
764  // liquid phase Dirichlet boundary
765  if (bc_l==1)
766  {
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;
776  }
777 
778  // liquid phase Neumann boundary
779  if (bc_l==0)
780  {
781  RF j_l = tp.j_l(*iit,face_local,time);
782  vn[iit->indexInInside()] = j_l;
783  }
784  }
785 
786  // compute coefficient
787  Dune::FieldVector<DF,dim> vstar=iit->unitOuterNormal(face_local); // normal on tranformef element
788  vstar *= vn[iit->indexInInside()];
789  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
790  if (iit->indexInInside()%2==0)
791  normalhat[iit->indexInInside()/2] = -1.0;
792  else
793  normalhat[iit->indexInInside()/2] = 1.0;
794  Dune::FieldVector<DF,dim> vstarhat(0);
795  B.umtv(vstar,vstarhat); // Piola backward transformation
796  vstarhat *= determinant;
797  coeff[iit->indexInInside()] = vstarhat*normalhat;
798  }
799 
800  // compute velocity on reference element
801  std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
802  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
803  typename Traits::RangeType yhat(0);
804  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
805  yhat.axpy(coeff[i],rt0vectors[i]);
806 
807  // apply Piola transformation
808  B.invert();
809  y = 0;
810  B.umtv(yhat,y);
811  y /= determinant;
812 
813  // std::cout << "vn= " ;
814  // for (int i=0; i<2*dim; i++) std::cout << vn[i] << " ";
815  // std::cout << std::endl;
816  // std::cout << "V_l: x=" << x << " y=" << y << std::endl;
817  }
818 
819  inline const typename Traits::GridViewType& getGridView ()
820  {
821  return pl.getGridView();
822  }
823 
824  private:
825 
826  template<typename T>
827  T aavg (T a, T b) const
828  {
829  return 0.5*(a+b);
830  }
831 
832  template<typename T>
833  T havg (T a, T b) const
834  {
835  T eps = 1E-30;
836  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
837  }
838  };
839 
848  template<typename TP, typename PL, typename PG>
849  class V_g
850  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
851  typename PL::Traits::RangeFieldType,
852  PL::Traits::GridViewType::dimension,
853  Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
854  V_g<TP,PL,PG> >
855  {
856  // extract useful types
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;
865 
866  const TP& tp;
867  const PL& pl;
868  const PG& pg;
869  Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
870  typename TP::Traits::RangeFieldType time;
871 
872 
873  typedef typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType RT0RangeType;
874 
875  public:
877 
879 
880  V_g (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
881 
882  // set time where operator is to be evaluated (i.e. end of the time intervall)
883  void set_time (typename TP::Traits::RangeFieldType time_)
884  {
885  time = time_;
886  }
887 
888  inline void evaluate (const typename Traits::ElementType& e,
889  const typename Traits::DomainType& x,
890  typename Traits::RangeType& y) const
891  {
892  // cell geometry
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);
898 
899  // absolute permeability
900  RF k_abs_inside = tp.k_abs(e,inside_cell_center_local);
901 
902  // pressure evaluation
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);
906 
907  // density evaluation
908  RF nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
909 
910  // for coefficient computation
911  RF vn[2*dim]; // normal velocities
912  RF coeff[2*dim]; // RT0 coefficient
913  typename Traits::ElementType::Geometry::JacobianInverseTransposed
914  B = e.geometry().jacobianInverseTransposed(x); // the transformation. Assume it is linear
915  RF determinant = B.determinant();
916 
917  // loop over cell neighbors
918  IntersectionIterator endit = pl.getGridView().iend(e);
919  for (IntersectionIterator iit = pl.getGridView().ibegin(e); iit!=endit; ++iit)
920  {
921  // set to zero for processor boundary
922  vn[iit->indexInInside()] = 0.0;
923 
924  // face geometry
925  const Dune::FieldVector<DF,dim-1>&
926  face_local = Dune::ReferenceElements<DF,dim-1>::general(iit->geometry().type()).position(0,0);
927 
928  // interior face
929  if (iit->neighbor())
930  {
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);
936 
937  // distance of cell centers
938  Dune::FieldVector<DF,dim> d(outside_cell_center_global);
939  d -= inside_cell_center_global;
940  RF distance = d.two_norm();
941 
942  // absolute permeability
943  RF k_abs_outside = tp.k_abs(*(iit->outside()),outside_cell_center_local);
944 
945  // gravity times normal
946  RF gn = tp.gravity()*iit->unitOuterNormal(face_local);
947 
948  // pressure evaluation
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);
952 
953  // needed for both phases
954  RF pc_upwind, s_l_upwind, s_g_upwind;
955 
956  // gas phase calculation
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; // determines direction
960  if (w_g>=0) // upwind capillary pressure on face
961  {
962  pc_upwind = pg_inside-pl_inside;
963  s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
964  }
965  else
966  {
967  pc_upwind = pg_outside-pl_outside;
968  s_l_upwind = tp.s_l(*(iit->outside()),outside_cell_center_local,pc_upwind);
969  }
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);
976 
977  RF nu_g_outside = tp.nu_g(*(iit->outside()),outside_cell_center_local,pg_outside);
978 
979  // set coefficient
980  vn[iit->indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
981  }
982 
983  // boundary face
984  if (iit->boundary())
985  {
986  // distance of cell center to boundary
987  Dune::FieldVector<DF,dim> d = iit->geometry().global(face_local);
988  d -= inside_cell_center_global;
989  RF distance = d.two_norm();
990 
991  // evaluate boundary condition type
992  int bc_g = tp.bc_g(*iit,face_local,time);
993 
994  // gas phase Dirichlet boundary
995  if (bc_g==1)
996  {
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);
1002  RF s_g = 1-s_l;
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;
1006 
1007  vn[iit->indexInInside()] = nu_g_inside * sigma_g * w_g;
1008  }
1009 
1010  // gas phase Neumann boundary
1011  if (bc_g==0)
1012  {
1013  RF j_g = tp.j_g(*iit,face_local,time);
1014  vn[iit->indexInInside()] = j_g; /* /nu_g_inside*/;
1015  }
1016  }
1017 
1018  // compute coefficient
1019  Dune::FieldVector<DF,dim> vstar=iit->unitOuterNormal(face_local); // normal on tranformef element
1020  vstar *= vn[iit->indexInInside()];
1021  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
1022  if (iit->indexInInside()%2==0)
1023  normalhat[iit->indexInInside()/2] = -1.0;
1024  else
1025  normalhat[iit->indexInInside()/2] = 1.0;
1026  Dune::FieldVector<DF,dim> vstarhat(0);
1027  B.umtv(vstar,vstarhat); // Piola backward transformation
1028  vstarhat *= determinant;
1029  coeff[iit->indexInInside()] = vstarhat*normalhat;
1030  }
1031 
1032  // compute velocity on reference element
1033  std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1034  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1035  typename Traits::RangeType yhat(0);
1036  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1037  yhat.axpy(coeff[i],rt0vectors[i]);
1038 
1039  // apply Piola transformation
1040  B.invert();
1041  y = 0;
1042  B.umtv(yhat,y);
1043  y /= determinant;
1044  }
1045 
1046  inline const typename Traits::GridViewType& getGridView ()
1047  {
1048  return pl.getGridView();
1049  }
1050 
1051  private:
1052 
1053  template<typename T>
1054  T aavg (T a, T b) const
1055  {
1056  return 0.5*(a+b);
1057  }
1058 
1059  template<typename T>
1060  T havg (T a, T b) const
1061  {
1062  T eps = 1E-30;
1063  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
1064  }
1065  };
1066 
1067 
1068  }
1069 }
1070 
1071 #endif
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
traits class for two phase parameter class
Definition: twophaseccfv.hh:22
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
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
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
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
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
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
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
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