dune-pdelab  2.5-dev
twophaseccfv.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_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 
13 
14 #include"../common/function.hh"
15 #include"pattern.hh"
16 #include"flags.hh"
17 #include"idefault.hh"
18 
19 namespace Dune {
20  namespace PDELab {
21 
23  template<typename GV, typename RF>
25  {
27  using GridViewType = GV;
28 
30  enum {
32  dimDomain = GV::dimension
33  };
34 
36  using DomainFieldType = typename GV::Grid::ctype;
37 
39  using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>;
40 
42  using IntersectionDomainType = Dune::FieldVector<DomainFieldType,dimDomain-1>;
43 
45  using RangeFieldType = RF;
46 
48  using RangeType = Dune::FieldVector<RF,GV::dimensionworld>;
49 
52 
54  using ElementType = typename GV::Traits::template Codim<0>::Entity;
55  using IntersectionType = typename GV::Intersection;
56  };
57 
58  template<typename GV, typename RF>
60  {
63 
65  using PermTensorType = Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>;
66  };
67 
69  template<class T, class Imp>
71  {
72  public:
73  using Traits = T;
74 
76  typename Traits::RangeFieldType
77  phi (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
78  {
79  return asImp().phi(e,x);
80  }
81 
83  typename Traits::RangeFieldType
84  pc (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
85  typename Traits::RangeFieldType s_l) const
86  {
87  return asImp().pc(e,x,s_l);
88  }
89 
91  typename Traits::RangeFieldType
92  s_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
93  typename Traits::RangeFieldType pc) const
94  {
95  return asImp().s_l(e,x,pc);
96  }
97 
99  typename Traits::RangeFieldType
100  kr_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
101  typename Traits::RangeFieldType s_l) const
102  {
103  return asImp().kr_l(e,x,s_l);
104  }
105 
107  typename Traits::RangeFieldType
108  kr_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
109  typename Traits::RangeFieldType s_g) const
110  {
111  return asImp().kr_g(e,x,s_g);
112  }
113 
115  typename Traits::RangeFieldType
116  mu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
117  typename Traits::RangeFieldType p_l) const
118  {
119  return asImp().mu_l(e,x,p_l);
120  }
121 
123  typename Traits::RangeFieldType
124  mu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
125  typename Traits::RangeFieldType p_g) const
126  {
127  return asImp().mu_l(e,x,p_g);
128  }
129 
131  typename Traits::PermTensorType
132  k_abs (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
133  {
134  return asImp().k_abs(e,x);
135  }
136 
138  const typename Traits::RangeType& gravity () const
139  {
140  return asImp().gravity();
141  }
142 
144  template<typename E>
145  typename Traits::RangeFieldType
146  nu_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
147  typename Traits::RangeFieldType p_l) const
148  {
149  return asImp().nu_l(e,x,p_l);
150  }
151 
153  typename Traits::RangeFieldType
154  nu_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
155  typename Traits::RangeFieldType p_g) const
156  {
157  return asImp().nu_g(e,x,p_g);
158  }
159 
161  typename Traits::RangeFieldType
162  rho_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
163  typename Traits::RangeFieldType p_l) const
164  {
165  return asImp().rho_l(e,x,p_l);
166  }
167 
169  typename Traits::RangeFieldType
170  rho_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
171  typename Traits::RangeFieldType p_g) const
172  {
173  return asImp().rho_g(e,x,p_g);
174  }
175 
177  int
178  bc_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
179  {
180  return asImp().bc_l(is,x,time);
181  }
182 
184  int
185  bc_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
186  {
187  return asImp().bc_g(is,x,time);
188  }
189 
191  typename Traits::RangeFieldType
192  g_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
193  {
194  return asImp().g_l(is,x,time);
195  }
196 
198  typename Traits::RangeFieldType
199  g_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
200  {
201  return asImp().g_g(is,x,time);
202  }
203 
205  typename Traits::RangeFieldType
206  j_l (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
207  {
208  return asImp().j_l(is,x,time);
209  }
210 
212  typename Traits::RangeFieldType
213  j_g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, typename Traits::RangeFieldType time) const
214  {
215  return asImp().j_g(is,x,time);
216  }
217 
219  typename Traits::RangeFieldType
220  q_l (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
221  typename Traits::RangeFieldType time) const
222  {
223  return asImp().q_l(e,x,time);
224  }
225 
227  typename Traits::RangeFieldType
228  q_g (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
229  typename Traits::RangeFieldType time) const
230  {
231  return asImp().q_g(e,x,time);
232  }
233 
234  private:
235  Imp& asImp () {return static_cast<Imp &> (*this);}
236  const Imp& asImp () const {return static_cast<const Imp &>(*this);}
237  };
238 
239 
240  // a local operator for solving the two-phase flow in pressure-pressure formulation
241  // with two-point flux approximation
242  // TP : parameter class, see above
243  // V : Vector holding last time step
244  template<typename TP>
246  : public NumericalJacobianSkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
247  public NumericalJacobianApplySkeleton<TwoPhaseTwoPointFluxOperator<TP> >,
248 
249  public NumericalJacobianBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
250  public NumericalJacobianApplyBoundary<TwoPhaseTwoPointFluxOperator<TP> >,
251 
252  public FullSkeletonPattern,
253  public FullVolumePattern,
255 
256  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
257  {
258  enum { dim = TP::Traits::GridViewType::dimension };
259  enum { liquid = 0 };
260  enum { gas = 1 };
261 
262  using Real = typename TP::Traits::RangeFieldType;
263  public:
264  // pattern assembly flags
265  enum { doPatternVolume = true };
266  enum { doPatternSkeleton = true };
267 
268  // residual assembly flags
269  enum { doAlphaSkeleton = true };
270  enum { doAlphaBoundary = true };
271  enum { doLambdaVolume = true };
272  enum { doLambdaBoundary = true };
273 
275  TwoPhaseTwoPointFluxOperator (const TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
276  : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
277  {}
278 
279  // volume integral depending only on test functions
280  template<typename EG, typename LFSV, typename R>
281  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
282  {
283  // Reference to cell
284  const auto& cell = eg.entity();
285 
286  // get geometry
287  auto geo = eg.geometry();
288 
289  // cell geometry
290  auto ref_el = referenceElement(geo);
291  auto cell_center_local = ref_el.position(0,0);
292  auto cell_volume = geo.volume();
293 
294  // contribution from source term
295  r.accumulate(lfsv, liquid, -scale_l * tp.q_l(cell,cell_center_local,time) * cell_volume);
296  r.accumulate(lfsv, gas, -scale_g * tp.q_g(cell,cell_center_local,time) * cell_volume);
297  }
298 
299  // skeleton integral depending on test and ansatz functions
300  // each face is only visited ONCE!
301  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
302  void alpha_skeleton (const IG& ig,
303  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
304  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
305  R& r_s, R& r_n) const
306  {
307  // Define types
308  using namespace Indices;
309  using PLSpace = TypeTree::Child<LFSV,_0>;
310  using RF = typename PLSpace::Traits::FiniteElementType::
311  Traits::LocalBasisType::Traits::RangeFieldType;
312 
313  // References to inside and outside cells
314  const auto& cell_inside = ig.inside();
315  const auto& cell_outside = ig.outside();
316 
317  // get geometries
318  auto geo = ig.geometry();
319  auto geo_inside = cell_inside.geometry();
320  auto geo_outside = cell_outside.geometry();
321 
322  // cell geometries
323  auto ref_el_inside = referenceElement(geo_inside);
324  auto ref_el_outside = referenceElement(geo_outside);
325  auto inside_cell_center_local = ref_el_inside.position(0,0);
326  auto outside_cell_center_local = ref_el_outside.position(0,0);
327  auto inside_cell_center_global = geo_inside.center();
328  auto outside_cell_center_global = geo_outside.center();
329 
330  // distance of cell centers
331  auto d = outside_cell_center_global;
332  d -= inside_cell_center_global;
333  auto distance = d.two_norm();
334 
335  // face geometry
336  auto ref_el = referenceElement(geo);
337  auto face_local = ref_el.position(0,0);
338  auto face_volume = geo.volume();
339 
340  // absolute permeability
341  auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
342  auto k_abs_outside = tp.k_abs(cell_outside,outside_cell_center_local);
343 
344  // gravity times normal
345  auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
346 
347  // liquid phase calculation
348  auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
349  auto rho_l_outside = tp.rho_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
350  auto w_l = (x_s(lfsu_s,liquid)-x_n(lfsu_n,liquid))/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction
351  RF pc_upwind, s_l_upwind, s_g_upwind;
352  auto nu_l = aavg(tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)),
353  tp.nu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid)));
354  if (w_l>=0) // upwind capillary pressure on face
355  {
356  pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
357  s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
358  }
359  else
360  {
361  pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
362  s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
363  }
364  s_g_upwind = 1-s_l_upwind;
365  auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l_upwind)/
366  tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
367  auto lambda_l_outside = tp.kr_l(cell_outside,outside_cell_center_local,s_l_upwind)/
368  tp.mu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid));
369  auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
370 
371  r_s.accumulate(lfsv_s,liquid, scale_l * nu_l * sigma_l * w_l * face_volume);
372  r_n.accumulate(lfsv_n,liquid, -scale_l * nu_l * sigma_l * w_l * face_volume);
373 
374  // gas phase calculation
375  auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
376  auto rho_g_outside = tp.rho_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
377  auto w_g = (x_s(lfsu_s,gas)-x_n(lfsu_n,gas))/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction
378  auto nu_g = aavg(tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)),
379  tp.nu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas)));
380  if (w_l*w_g<=0) // new evaluation necessary only if signs differ
381  {
382  if (w_g>=0) // upwind capillary pressure on face
383  {
384  pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid);
385  s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind);
386  }
387  else
388  {
389  pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid);
390  s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind);
391  }
392  s_g_upwind = 1-s_l_upwind;
393  }
394  auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g_upwind)/
395  tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
396  auto lambda_g_outside = tp.kr_g(cell_outside,outside_cell_center_local,s_g_upwind)/
397  tp.mu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas));
398  auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
399 
400  r_s.accumulate(lfsv_s, gas, scale_g * nu_g * sigma_g * w_g * face_volume);
401  r_n.accumulate(lfsv_n, gas, -scale_g * nu_g * sigma_g * w_g * face_volume);
402  }
403 
404  // skeleton integral depending on test and ansatz functions
405  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
406  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
407  void alpha_boundary (const IG& ig,
408  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
409  R& r_s) const
410  {
411  // References to inside cell
412  const auto& cell_inside = ig.inside();
413 
414  // get geometries
415  auto geo = ig.geometry();
416  auto geo_inside = cell_inside.geometry();
417 
418  // face geometry
419  auto ref_el = referenceElement(geo);
420  auto face_local = ref_el.position(0,0);
421  auto face_volume = geo.volume();
422 
423  // evaluate boundary condition type
424  auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
425  auto 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  auto ref_el_inside = referenceElement(geo_inside);
430  auto inside_cell_center_local = ref_el_inside.position(0,0);
431  auto inside_cell_center_global = geo_inside.center();
432 
433  // distance of cell center to boundary
434  auto d = geo.global(face_local);
435  d -= inside_cell_center_global;
436  auto distance = d.two_norm();
437 
438  // absolute permeability
439  auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local);
440 
441  // gravity times normal
442  auto gn = tp.gravity()*ig.unitOuterNormal(face_local);
443 
444  // liquid phase Dirichlet boundary
445  if (bc_l==1)
446  {
447  auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
448  auto g_l = tp.g_l(ig.intersection(),face_local,time);
449  auto w_l = (x_s(lfsu_s,liquid)-g_l)/distance + rho_l_inside*gn;
450  auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
451  auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l)/
452  tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid));
453  auto sigma_l = lambda_l_inside*k_abs_inside;
454  auto nu_l = tp.nu_l(cell_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  auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
462  auto g_g = tp.g_g(ig.intersection(),face_local,time);
463  auto w_g = (x_s(lfsu_s,gas)-g_g)/distance + rho_g_inside*gn;
464  auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid));
465  auto s_g = 1-s_l;
466  auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g)/
467  tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas));
468  auto sigma_g = lambda_g_inside*k_abs_inside;
469  auto nu_g = tp.nu_g(cell_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  // get geometries
479  auto geo = ig.geometry();
480 
481  // face geometry
482  auto ref_el = referenceElement(geo);
483  auto face_local = ref_el.position(0,0);
484  auto face_volume = geo.volume();
485 
486  // evaluate boundary condition type
487  auto bc_l = tp.bc_l(ig.intersection(),face_local,time);
488  auto bc_g = tp.bc_g(ig.intersection(),face_local,time);
489  if (bc_l!=0 && bc_g!=0) return; // no Neumann boundary conditions
490 
491  // liquid phase Neumann boundary
492  if (bc_l==0)
493  {
494  auto j_l = tp.j_l(ig.intersection(),face_local,time);
495  r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume);
496  }
497 
498  // gas phase Neumann boundary
499  if (bc_g==0)
500  {
501  auto j_g = tp.j_g(ig.intersection(),face_local,time);
502  r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume);
503  }
504  }
505 
507  void setTime (typename TP::Traits::RangeFieldType t)
508  {
509  time = t;
510  }
511 
512  private:
513  const TP& tp; // two phase parameter class
514  typename TP::Traits::RangeFieldType time;
515  Real scale_l, scale_g;
516 
517  template<typename T>
518  T aavg (T a, T b) const
519  {
520  return 0.5*(a+b);
521  }
522 
523  template<typename T>
524  T havg (T a, T b) const
525  {
526  T eps = 1E-30;
527  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
528  }
529  };
530 
531 
538  template<class TP>
540  : public NumericalJacobianVolume<TwoPhaseOnePointTemporalOperator<TP> >,
541  public NumericalJacobianApplyVolume<TwoPhaseOnePointTemporalOperator<TP> >,
542  public FullVolumePattern,
544  public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType>
545  {
546  enum { dim = TP::Traits::GridViewType::dimension };
547  enum { liquid = 0 };
548  enum { gas = 1 };
549 
550  using Real = typename TP::Traits::RangeFieldType;
551 
552  public:
553  // pattern assembly flags
554  enum { doPatternVolume = true };
555 
556  // residual assembly flags
557  enum { doAlphaVolume = true };
558 
559  TwoPhaseOnePointTemporalOperator (TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0)
560  : tp(tp_), scale_l(scale_l_), scale_g(scale_g_)
561  {
562  }
563 
565  void setTime (typename TP::Traits::RangeFieldType t)
566  {
567  time = t;
568  }
569 
570  // volume integral depending on test and ansatz functions
571  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
572  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
573  {
574  // Reference to cell
575  const auto& cell = eg.entity();
576 
577  // get geometry
578  auto geo = eg.geometry();
579 
580  // cell geometry
581  auto ref_el = referenceElement(geo);
582  auto cell_center_local = ref_el.position(0,0);
583  auto cell_volume = geo.volume();
584 
585  auto phi = tp.phi(cell,cell_center_local);
586  auto s_l = tp.s_l(cell,cell_center_local,x(lfsu,gas)-x(lfsu,liquid));
587 
588  r.accumulate(lfsv, liquid, scale_l * phi * s_l * tp.nu_l(cell,cell_center_local,x(lfsu,liquid)) * cell_volume);
589  r.accumulate(lfsv, gas, scale_g * phi * (1-s_l) * tp.nu_g(cell,cell_center_local,x(lfsu,gas)) * cell_volume);
590  }
591 
592  private:
593  TP& tp;
594  typename TP::Traits::RangeFieldType time;
595  Real scale_l, scale_g;
596  };
597 
598 
607  template<typename TP, typename PL, typename PG>
608  class V_l
609  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
610  typename PL::Traits::RangeFieldType,
611  PL::Traits::GridViewType::dimension,
612  Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
613  V_l<TP,PL,PG> >
614  {
615  // extract useful types
616  using GV = typename PL::Traits::GridViewType;
617  using DF = typename GV::Grid::ctype;
618  using RF = typename PL::Traits::RangeFieldType;
619  using RangeType = typename PL::Traits::RangeType;
620  enum { dim = PL::Traits::GridViewType::dimension };
621  using Element = typename GV::Traits::template Codim<0>::Entity;
622  using IntersectionIterator = typename GV::IntersectionIterator;
623  using Intersection = typename IntersectionIterator::Intersection;
624 
625  const TP& tp;
626  const PL& pl;
627  const PG& pg;
628  Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
629  typename TP::Traits::RangeFieldType time;
630 
631 
632  using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
633 
634  public:
636 
638 
639  V_l (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
640 
641  // set time where operator is to be evaluated (i.e. end of the time intervall)
642  void set_time (typename TP::Traits::RangeFieldType time_)
643  {
644  time = time_;
645  }
646 
647  inline void evaluate (const typename Traits::ElementType& e,
648  const typename Traits::DomainType& x,
649  typename Traits::RangeType& y) const
650  {
651  // get geometries
652  auto geo = e.geometry();
653 
654  // face geometry
655  auto ref_el = referenceElement(geo);
656  auto face_local = ref_el.position(0,0);
657  auto face_volume = geo.volume();
658 
659  // cell geometry
660  auto inside_cell_center_local = ref_el.position(0,0);
661  auto inside_cell_center_global = geo.global(inside_cell_center_local);
662 
663  // absolute permeability
664  auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
665 
666  // pressure evaluation
667  typename PL::Traits::RangeType pl_inside, pg_inside;
668  pl.evaluate(e,inside_cell_center_local,pl_inside);
669  pg.evaluate(e,inside_cell_center_local,pg_inside);
670 
671  // density
672  auto nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside);
673 
674  // for coefficient computation
675  RF vn[2*dim]; // normal velocities
676  RF coeff[2*dim]; // RT0 coefficient
677  auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear
678  auto determinant = B.determinant();
679 
680  // loop over cell neighbors
681  for (const auto& intersection : intersections(pl.getGridView(),e))
682  {
683  // set to zero for processor boundary
684  vn[intersection.indexInInside()] = 0.0;
685 
686  // get geometry
687  auto geo_intersection = intersection.geometry();
688 
689  // face geometry
690  const auto& face_local = referenceElement(geo_intersection).position(0,0);
691 
692  // interior face
693  if (intersection.neighbor())
694  {
695  // get geometry
696  auto geo_outside = intersection.outside().geometry();
697  auto ref_el_outside = referenceElement(geo_outside);
698  auto outside_cell_center_local = ref_el_outside.position(0,0);
699  auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
700 
701  // distance of cell centers
702  auto d = outside_cell_center_global;
703  d -= inside_cell_center_global;
704  auto distance = d.two_norm();
705 
706  // absolute permeability
707  auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
708 
709  // gravity times normal
710  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
711 
712  // pressure evaluation
713  typename PL::Traits::RangeType pl_outside, pg_outside;
714  pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
715  pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
716 
717  // density
718  auto nu_l_outside = tp.nu_l(*(intersection.outside()),outside_cell_center_local,pg_outside);
719 
720  // liquid phase calculation
721  auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
722  auto rho_l_outside = tp.rho_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
723  auto w_l = (pl_inside-pl_outside)/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction
724  RF pc_upwind, s_l_upwind, s_g_upwind;
725  if (w_l>=0) // upwind capillary pressure on face
726  {
727  pc_upwind = pg_inside-pl_inside;
728  s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
729  }
730  else
731  {
732  pc_upwind = pg_outside-pl_outside;
733  s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
734  }
735  s_g_upwind = 1-s_l_upwind;
736  auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l_upwind)/
737  tp.mu_l(e,inside_cell_center_local,pl_inside);
738  auto lambda_l_outside = tp.kr_l(*(intersection.outside()),outside_cell_center_local,s_l_upwind)/
739  tp.mu_l(*(intersection.outside()),outside_cell_center_local,pl_outside);
740  auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside);
741  auto nu_l = aavg(nu_l_inside,nu_l_outside);
742 
743  // set coefficient
744  vn[intersection.indexInInside()] = nu_l * sigma_l * w_l;
745  }
746 
747  // boundary face
748  if (intersection.boundary())
749  {
750  // distance of cell center to boundary
751  auto d = geo_intersection.global(face_local);
752  d -= inside_cell_center_global;
753  auto distance = d.two_norm();
754 
755  // evaluate boundary condition type
756  auto bc_l = tp.bc_l(intersection,face_local,time);
757 
758  // liquid phase Dirichlet boundary
759  if (bc_l==1)
760  {
761  auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside);
762  auto g_l = tp.g_l(intersection,face_local,time);
763  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
764  auto w_l = (pl_inside-g_l)/distance + rho_l_inside*gn;
765  auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
766  auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l)/
767  tp.mu_l(e,inside_cell_center_local,pl_inside);
768  auto sigma_l = lambda_l_inside*k_abs_inside;
769  vn[intersection.indexInInside()] = nu_l_inside * sigma_l * w_l;
770  }
771 
772  // liquid phase Neumann boundary
773  if (bc_l==0)
774  {
775  auto j_l = tp.j_l(intersection,face_local,time);
776  vn[intersection.indexInInside()] = j_l;
777  }
778  }
779 
780  // compute coefficient
781  auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element
782  vstar *= vn[intersection.indexInInside()];
783  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
784  if (intersection.indexInInside()%2==0)
785  normalhat[intersection.indexInInside()/2] = -1.0;
786  else
787  normalhat[intersection.indexInInside()/2] = 1.0;
788  Dune::FieldVector<DF,dim> vstarhat(0);
789  B.umtv(vstar,vstarhat); // Piola backward transformation
790  vstarhat *= determinant;
791  coeff[intersection.indexInInside()] = vstarhat*normalhat;
792  }
793 
794  // compute velocity on reference element
795  std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
796  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
797  typename Traits::RangeType yhat(0);
798  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
799  yhat.axpy(coeff[i],rt0vectors[i]);
800 
801  // apply Piola transformation
802  B.invert();
803  y = 0;
804  B.umtv(yhat,y);
805  y /= determinant;
806 
807  // std::cout << "vn= " ;
808  // for (int i=0; i<2*dim; i++) std::cout << vn[i] << " ";
809  // std::cout << std::endl;
810  // std::cout << "V_l: x=" << x << " y=" << y << std::endl;
811  }
812 
813  inline const typename Traits::GridViewType& getGridView ()
814  {
815  return pl.getGridView();
816  }
817 
818  private:
819 
820  template<typename T>
821  T aavg (T a, T b) const
822  {
823  return 0.5*(a+b);
824  }
825 
826  template<typename T>
827  T havg (T a, T b) const
828  {
829  T eps = 1E-30;
830  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
831  }
832  };
833 
842  template<typename TP, typename PL, typename PG>
843  class V_g
844  : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType,
845  typename PL::Traits::RangeFieldType,
846  PL::Traits::GridViewType::dimension,
847  Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >,
848  V_g<TP,PL,PG> >
849  {
850  // extract useful types
851  using GV = typename PL::Traits::GridViewType;
852  using DF = typename GV::Grid::ctype;
853  using RF = typename PL::Traits::RangeFieldType;
854  using RangeType = typename PL::Traits::RangeType;
855  enum { dim = PL::Traits::GridViewType::dimension };
856  using Element = typename GV::Traits::template Codim<0>::Entity;
857  using IntersectionIterator = typename GV::IntersectionIterator;
858  using Intersection = typename IntersectionIterator::Intersection;
859 
860  const TP& tp;
861  const PL& pl;
862  const PG& pg;
863  Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe;
864  typename TP::Traits::RangeFieldType time;
865 
866 
867  using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType;
868 
869  public:
871 
873 
874  V_g (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {}
875 
876  // set time where operator is to be evaluated (i.e. end of the time intervall)
877  void set_time (typename TP::Traits::RangeFieldType time_)
878  {
879  time = time_;
880  }
881 
882  inline void evaluate (const typename Traits::ElementType& e,
883  const typename Traits::DomainType& x,
884  typename Traits::RangeType& y) const
885  {
886  // get geometries
887  auto geo = e.geometry();
888 
889  // face geometry
890  auto ref_el = referenceElement(geo);
891  auto face_local = ref_el.position(0,0);
892  auto face_volume = geo.volume();
893 
894  // cell geometry
895  auto inside_cell_center_local = ref_el.position(0,0);
896  auto inside_cell_center_global = geo.global(inside_cell_center_local);
897 
898  // absolute permeability
899  auto k_abs_inside = tp.k_abs(e,inside_cell_center_local);
900 
901  // pressure evaluation
902  typename PL::Traits::RangeType pl_inside, pg_inside;
903  pl.evaluate(e,inside_cell_center_local,pl_inside);
904  pg.evaluate(e,inside_cell_center_local,pg_inside);
905 
906  // density evaluation
907  auto nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside);
908 
909  // for coefficient computation
910  RF vn[2*dim]; // normal velocities
911  RF coeff[2*dim]; // RT0 coefficient
912  auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear
913  auto determinant = B.determinant();
914 
915  // loop over cell neighbors
916  for (const auto& intersection : intersections(pl.getGridView(),e))
917  {
918  // set to zero for processor boundary
919  vn[intersection.indexInInside()] = 0.0;
920 
921  // get geometry
922  auto geo_intersection = intersection.geometry();
923 
924  // face geometry
925  const auto& face_local = referenceElement(geo_intersection).position(0,0);
926 
927  // interior face
928  if (intersection.neighbor())
929  {
930  // get geometry
931  auto geo_outside = intersection.outside().geometry();
932  auto ref_el_outside = referenceElement(geo_outside);
933  auto outside_cell_center_local = ref_el_outside.position(0,0);
934  auto outside_cell_center_global = geo_outside.global(outside_cell_center_local);
935 
936  // distance of cell centers
937  auto d = outside_cell_center_global;
938  d -= inside_cell_center_global;
939  auto distance = d.two_norm();
940 
941  // absolute permeability
942  auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local);
943 
944  // gravity times normal
945  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
946 
947  // pressure evaluation
948  typename PL::Traits::RangeType pl_outside, pg_outside;
949  pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside);
950  pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside);
951 
952  // needed for both phases
953  RF pc_upwind, s_l_upwind, s_g_upwind;
954 
955  // gas phase calculation
956  auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
957  auto rho_g_outside = tp.rho_g(e,outside_cell_center_local,pg_outside);
958  auto w_g = (pg_inside-pg_outside)/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction
959  if (w_g>=0) // upwind capillary pressure on face
960  {
961  pc_upwind = pg_inside-pl_inside;
962  s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind);
963  }
964  else
965  {
966  pc_upwind = pg_outside-pl_outside;
967  s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind);
968  }
969  s_g_upwind = 1-s_l_upwind;
970  auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g_upwind)/
971  tp.mu_g(e,inside_cell_center_local,pg_inside);
972  auto lambda_g_outside = tp.kr_g(*(intersection.outside()),outside_cell_center_local,s_g_upwind)/
973  tp.mu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
974  auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside);
975 
976  auto nu_g_outside = tp.nu_g(*(intersection.outside()),outside_cell_center_local,pg_outside);
977 
978  // set coefficient
979  vn[intersection.indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g;
980  }
981 
982  // boundary face
983  if (intersection.boundary())
984  {
985  // distance of cell center to boundary
986  auto d = geo_intersection.global(face_local);
987  d -= inside_cell_center_global;
988  auto distance = d.two_norm();
989 
990  // evaluate boundary condition type
991  auto bc_g = tp.bc_g(intersection,face_local,time);
992 
993  // gas phase Dirichlet boundary
994  if (bc_g==1)
995  {
996  auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside);
997  auto g_g = tp.g_g(intersection,face_local,time);
998  auto gn = tp.gravity()*intersection.unitOuterNormal(face_local);
999  auto w_g = (pg_inside-g_g)/distance + rho_g_inside*gn;
1000  auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside);
1001  auto s_g = 1-s_l;
1002  auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g)/
1003  tp.mu_g(e,inside_cell_center_local,pg_inside);
1004  auto sigma_g = lambda_g_inside*k_abs_inside;
1005 
1006  vn[intersection.indexInInside()] = nu_g_inside * sigma_g * w_g;
1007  }
1008 
1009  // gas phase Neumann boundary
1010  if (bc_g==0)
1011  {
1012  auto j_g = tp.j_g(intersection,face_local,time);
1013  vn[intersection.indexInInside()] = j_g; /* /nu_g_inside*/;
1014  }
1015  }
1016 
1017  // compute coefficient
1018  auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element
1019  vstar *= vn[intersection.indexInInside()];
1020  Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element
1021  if (intersection.indexInInside()%2==0)
1022  normalhat[intersection.indexInInside()/2] = -1.0;
1023  else
1024  normalhat[intersection.indexInInside()/2] = 1.0;
1025  Dune::FieldVector<DF,dim> vstarhat(0);
1026  B.umtv(vstar,vstarhat); // Piola backward transformation
1027  vstarhat *= determinant;
1028  coeff[intersection.indexInInside()] = vstarhat*normalhat;
1029  }
1030 
1031  // compute velocity on reference element
1032  std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size());
1033  rt0fe.localBasis().evaluateFunction(x,rt0vectors);
1034  typename Traits::RangeType yhat(0);
1035  for (unsigned int i=0; i<rt0fe.localBasis().size(); i++)
1036  yhat.axpy(coeff[i],rt0vectors[i]);
1037 
1038  // apply Piola transformation
1039  B.invert();
1040  y = 0;
1041  B.umtv(yhat,y);
1042  y /= determinant;
1043  }
1044 
1045  inline const typename Traits::GridViewType& getGridView ()
1046  {
1047  return pl.getGridView();
1048  }
1049 
1050  private:
1051 
1052  template<typename T>
1053  T aavg (T a, T b) const
1054  {
1055  return 0.5*(a+b);
1056  }
1057 
1058  template<typename T>
1059  T havg (T a, T b) const
1060  {
1061  T eps = 1E-30;
1062  return 2.0/(1.0/(a+eps) + 1.0/(b+eps));
1063  }
1064  };
1065 
1066 
1067  }
1068 }
1069 
1070 #endif // DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH
const IG & ig
Definition: constraints.hh:148
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r_s) const
Definition: twophaseccfv.hh:476
static const int dim
Definition: adaptivity.hh:83
T Traits
Definition: twophaseccfv.hh:73
leaf of a function tree
Definition: function.hh:575
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:220
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
sparsity pattern generator
Definition: pattern.hh:29
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:407
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:228
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
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:170
Provide velocity field for gas phase.
Definition: twophaseccfv.hh:843
dimension of the domain
Definition: twophaseccfv.hh:32
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:117
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: twophaseccfv.hh:36
Definition: twophaseccfv.hh:245
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:877
traits class for two phase parameter class
Definition: twophaseccfv.hh:24
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: twophaseccfv.hh:42
RangeFieldType PermTensorType
permeability tensor type
Definition: twophaseccfv.hh:51
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary() ...
Definition: numericaljacobianapply.hh:285
V_g(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:874
V_l(const TP &tp_, const PL &pl_, const PG &pg_)
Definition: twophaseccfv.hh:639
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:813
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume() ...
Definition: numericaljacobianapply.hh:32
TwoPhaseTwoPointFluxOperator(const TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
constructor: pass parameter object
Definition: twophaseccfv.hh:275
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:199
Provide velocity field for liquid phase.
Definition: twophaseccfv.hh:608
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:185
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:100
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: twophaseccfv.hh:48
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: twophaseccfv.hh:39
const Traits::RangeType & gravity() const
gravity vector
Definition: twophaseccfv.hh:138
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:162
const Traits::GridViewType & getGridView()
Definition: twophaseccfv.hh:1045
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:114
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:281
Implements linear and nonlinear versions of jacobian_apply_skeleton() based on alpha_skeleton() ...
Definition: numericaljacobianapply.hh:180
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:647
IntersectionType
Enum describing the type of an intersection.
Definition: intersectiontype.hh:14
void set_time(typename TP::Traits::RangeFieldType time_)
Definition: twophaseccfv.hh:642
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Definition: twophaseccfv.hh:882
RF RangeFieldType
Export type for range field.
Definition: twophaseccfv.hh:45
Traits::PermTensorType k_abs(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
absolute permeability (scalar!)
Definition: twophaseccfv.hh:132
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:178
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:92
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:206
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: twophaseccfv.hh:572
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: twophaseccfv.hh:54
Traits::RangeFieldType phi(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
porosity
Definition: twophaseccfv.hh:77
GV GridViewType
the grid view
Definition: twophaseccfv.hh:27
Dune::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:48
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:108
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
sparsity pattern generator
Definition: pattern.hh:13
const Entity & e
Definition: localfunctionspace.hh:111
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:565
TwoPhaseOnePointTemporalOperator(TP &tp_, Real scale_l_=1.0, Real scale_g_=1.0)
Definition: twophaseccfv.hh:559
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:84
typename GV::Intersection IntersectionType
Definition: twophaseccfv.hh:55
base class for parameter class
Definition: twophaseccfv.hh:70
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:192
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:213
void setTime(typename TP::Traits::RangeFieldType t)
set time for subsequent evaluation
Definition: twophaseccfv.hh:507
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:116
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:146
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
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:154
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
Definition: twophaseccfv.hh:539
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:124
traits class holding the function signature, same as in local function
Definition: function.hh:175
ReferenceElementWrapper< ReferenceElement< typename Geometry::ctype, Geometry::mydimension > > referenceElement(const Geometry &geo)
Returns the reference element for the given geometry.
Definition: referenceelements.hh:144
Default flags for all local operators.
Definition: flags.hh:18
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:302