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