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 10 #include<dune/pdelab/common/quadraturerules.hh> 11 #include <dune/pdelab/localoperator/defaultimp.hh> 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 21 //! traits class for two phase parameter class 22 template<typename GV, typename RF> 23 struct TwoPhaseParameterTraits 24 { 25 //! \brief the grid view 26 using GridViewType = GV; 27 28 //! \brief Enum for domain dimension 29 enum { 30 //! \brief dimension of the domain 31 dimDomain = GV::dimension 32 }; 33 34 //! \brief Export type for domain field 35 using DomainFieldType = typename GV::Grid::ctype; 36 37 //! \brief domain type 38 using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>; 39 40 //! \brief domain type 41 using IntersectionDomainType = Dune::FieldVector<DomainFieldType,dimDomain-1>; 42 43 //! \brief Export type for range field 44 using RangeFieldType = RF; 45 46 //! \brief range type 47 using RangeType = Dune::FieldVector<RF,GV::dimensionworld>; 48 49 //! \brief permeability tensor type 50 using PermTensorType = RangeFieldType; 51 52 //! grid types 53 using ElementType = typename GV::Traits::template Codim<0>::Entity; 54 using IntersectionType = typename GV::Intersection; 55 }; 56 57 template<typename GV, typename RF> 58 struct TwoPhaseFullTensorParameterTraits : TwoPhaseParameterTraits<GV, RF> 59 { 60 using Base = TwoPhaseParameterTraits<GV, RF>; 61 using RangeFieldType = typename Base::RangeFieldType; 62 63 //! \brief permeability tensor type 64 using PermTensorType = Dune::FieldMatrix<RangeFieldType,Base::dimDomain,Base::dimDomain>; 65 }; 66 67 //! base class for parameter class 68 template<class T, class Imp> 69 class TwoPhaseParameterInterface 70 { 71 public: 72 using Traits = T; 73 74 //! porosity 75 typename Traits::RangeFieldType phi(const typename Traits::ElementType & e,const typename Traits::DomainType & x) const76 phi (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const 77 { 78 return asImp().phi(e,x); 79 } 80 81 //! capillary pressure function 82 typename Traits::RangeFieldType pc(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType s_l) const83 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 89 //! inverse capillary pressure function 90 typename Traits::RangeFieldType s_l(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType pc) const91 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 97 //! liquid phase relative permeability 98 typename Traits::RangeFieldType kr_l(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType s_l) const99 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 105 //! gas phase relative permeability 106 typename Traits::RangeFieldType kr_g(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType s_g) const107 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 113 //! liquid phase viscosity 114 typename Traits::RangeFieldType mu_l(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType p_l) const115 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 121 //! gas phase viscosity 122 typename Traits::RangeFieldType mu_g(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType p_g) const123 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 129 //! absolute permeability (scalar!) 130 typename Traits::PermTensorType k_abs(const typename Traits::ElementType & e,const typename Traits::DomainType & x) const131 k_abs (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const 132 { 133 return asImp().k_abs(e,x); 134 } 135 136 //! gravity vector gravity() const137 const typename Traits::RangeType& gravity () const 138 { 139 return asImp().gravity(); 140 } 141 142 //! liquid phase molar density 143 template<typename E> 144 typename Traits::RangeFieldType nu_l(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType p_l) const145 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 151 //! gas phase molar density 152 typename Traits::RangeFieldType nu_g(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType p_g) const153 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 159 //! liquid phase mass density 160 typename Traits::RangeFieldType rho_l(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType p_l) const161 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 167 //! gas phase mass density 168 typename Traits::RangeFieldType rho_g(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType p_g) const169 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 175 //! liquid phase boundary condition type 176 int bc_l(const typename Traits::IntersectionType & is,const typename Traits::IntersectionDomainType & x,typename Traits::RangeFieldType time) const177 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 182 //! gas phase boundary condition type 183 int bc_g(const typename Traits::IntersectionType & is,const typename Traits::IntersectionDomainType & x,typename Traits::RangeFieldType time) const184 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 189 //! liquid phase Dirichlet boundary condition 190 typename Traits::RangeFieldType g_l(const typename Traits::IntersectionType & is,const typename Traits::IntersectionDomainType & x,typename Traits::RangeFieldType time) const191 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 196 //! gas phase Dirichlet boundary condition 197 typename Traits::RangeFieldType g_g(const typename Traits::IntersectionType & is,const typename Traits::IntersectionDomainType & x,typename Traits::RangeFieldType time) const198 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 203 //! liquid phase Neumann boundary condition 204 typename Traits::RangeFieldType j_l(const typename Traits::IntersectionType & is,const typename Traits::IntersectionDomainType & x,typename Traits::RangeFieldType time) const205 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 210 //! gas phase Neumann boundary condition 211 typename Traits::RangeFieldType j_g(const typename Traits::IntersectionType & is,const typename Traits::IntersectionDomainType & x,typename Traits::RangeFieldType time) const212 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 217 //! liquid phase source term 218 typename Traits::RangeFieldType q_l(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType time) const219 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 225 //! gas phase source term 226 typename Traits::RangeFieldType q_g(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeFieldType time) const227 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: asImp()234 Imp& asImp () {return static_cast<Imp &> (*this);} asImp() const235 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> 244 class TwoPhaseTwoPointFluxOperator 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, 253 public LocalOperatorDefaultFlags, 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 using Real = typename TP::Traits::RangeFieldType; 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 273 //! constructor: pass parameter object TwoPhaseTwoPointFluxOperator(const TP & tp_,Real scale_l_=1.0,Real scale_g_=1.0)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> lambda_volume(const EG & eg,const LFSV & lfsv,R & r) const280 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const 281 { 282 // Reference to cell 283 const auto& cell = eg.entity(); 284 285 // get geometry 286 auto geo = eg.geometry(); 287 288 // cell geometry 289 auto ref_el = referenceElement(geo); 290 auto cell_center_local = ref_el.position(0,0); 291 auto cell_volume = geo.volume(); 292 293 // contribution from source term 294 r.accumulate(lfsv, liquid, -scale_l * tp.q_l(cell,cell_center_local,time) * cell_volume); 295 r.accumulate(lfsv, gas, -scale_g * tp.q_g(cell,cell_center_local,time) * cell_volume); 296 } 297 298 // skeleton integral depending on test and ansatz functions 299 // each face is only visited ONCE! 300 template<typename IG, typename LFSU, typename X, typename LFSV, typename R> 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) const301 void alpha_skeleton (const IG& ig, 302 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s, 303 const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n, 304 R& r_s, R& r_n) const 305 { 306 // Define types 307 using namespace Indices; 308 using PLSpace = TypeTree::Child<LFSV,_0>; 309 using RF = typename PLSpace::Traits::FiniteElementType:: 310 Traits::LocalBasisType::Traits::RangeFieldType; 311 312 // References to inside and outside cells 313 const auto& cell_inside = ig.inside(); 314 const auto& cell_outside = ig.outside(); 315 316 // get geometries 317 auto geo = ig.geometry(); 318 auto geo_inside = cell_inside.geometry(); 319 auto geo_outside = cell_outside.geometry(); 320 321 // cell geometries 322 auto ref_el_inside = referenceElement(geo_inside); 323 auto ref_el_outside = referenceElement(geo_outside); 324 auto inside_cell_center_local = ref_el_inside.position(0,0); 325 auto outside_cell_center_local = ref_el_outside.position(0,0); 326 auto inside_cell_center_global = geo_inside.center(); 327 auto outside_cell_center_global = geo_outside.center(); 328 329 // distance of cell centers 330 auto d = outside_cell_center_global; 331 d -= inside_cell_center_global; 332 auto distance = d.two_norm(); 333 334 // face geometry 335 auto ref_el = referenceElement(geo); 336 auto face_local = ref_el.position(0,0); 337 auto face_volume = geo.volume(); 338 339 // absolute permeability 340 auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local); 341 auto k_abs_outside = tp.k_abs(cell_outside,outside_cell_center_local); 342 343 // gravity times normal 344 auto gn = tp.gravity()*ig.unitOuterNormal(face_local); 345 346 // liquid phase calculation 347 auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)); 348 auto rho_l_outside = tp.rho_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid)); 349 auto w_l = (x_s(lfsu_s,liquid)-x_n(lfsu_n,liquid))/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction 350 RF pc_upwind, s_l_upwind, s_g_upwind; 351 auto nu_l = aavg(tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)), 352 tp.nu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid))); 353 if (w_l>=0) // upwind capillary pressure on face 354 { 355 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid); 356 s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind); 357 } 358 else 359 { 360 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid); 361 s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind); 362 } 363 s_g_upwind = 1-s_l_upwind; 364 auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l_upwind)/ 365 tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)); 366 auto lambda_l_outside = tp.kr_l(cell_outside,outside_cell_center_local,s_l_upwind)/ 367 tp.mu_l(cell_outside,outside_cell_center_local,x_n(lfsu_n,liquid)); 368 auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside); 369 370 r_s.accumulate(lfsv_s,liquid, scale_l * nu_l * sigma_l * w_l * face_volume); 371 r_n.accumulate(lfsv_n,liquid, -scale_l * nu_l * sigma_l * w_l * face_volume); 372 373 // gas phase calculation 374 auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)); 375 auto rho_g_outside = tp.rho_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas)); 376 auto w_g = (x_s(lfsu_s,gas)-x_n(lfsu_n,gas))/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction 377 auto nu_g = aavg(tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)), 378 tp.nu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas))); 379 if (w_l*w_g<=0) // new evaluation necessary only if signs differ 380 { 381 if (w_g>=0) // upwind capillary pressure on face 382 { 383 pc_upwind = x_s(lfsu_s,gas)-x_s(lfsu_s,liquid); 384 s_l_upwind = tp.s_l(cell_inside,inside_cell_center_local,pc_upwind); 385 } 386 else 387 { 388 pc_upwind = x_n(lfsu_n,gas)-x_n(lfsu_n,liquid); 389 s_l_upwind = tp.s_l(cell_outside,outside_cell_center_local,pc_upwind); 390 } 391 s_g_upwind = 1-s_l_upwind; 392 } 393 auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g_upwind)/ 394 tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)); 395 auto lambda_g_outside = tp.kr_g(cell_outside,outside_cell_center_local,s_g_upwind)/ 396 tp.mu_g(cell_outside,outside_cell_center_local,x_n(lfsu_n,gas)); 397 auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside); 398 399 r_s.accumulate(lfsv_s, gas, scale_g * nu_g * sigma_g * w_g * face_volume); 400 r_n.accumulate(lfsv_n, gas, -scale_g * nu_g * sigma_g * w_g * face_volume); 401 } 402 403 // skeleton integral depending on test and ansatz functions 404 // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations 405 template<typename IG, typename LFSU, typename X, typename LFSV, typename R> alpha_boundary(const IG & ig,const LFSU & lfsu_s,const X & x_s,const LFSV & lfsv_s,R & r_s) const406 void alpha_boundary (const IG& ig, 407 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s, 408 R& r_s) const 409 { 410 // References to inside cell 411 const auto& cell_inside = ig.inside(); 412 413 // get geometries 414 auto geo = ig.geometry(); 415 auto geo_inside = cell_inside.geometry(); 416 417 // face geometry 418 auto ref_el = referenceElement(geo); 419 auto face_local = ref_el.position(0,0); 420 auto face_volume = geo.volume(); 421 422 // evaluate boundary condition type 423 auto bc_l = tp.bc_l(ig.intersection(),face_local,time); 424 auto bc_g = tp.bc_g(ig.intersection(),face_local,time); 425 if (bc_l!=1 && bc_g!=1) return; // no Dirichlet boundary conditions 426 427 // cell geometry 428 auto ref_el_inside = referenceElement(geo_inside); 429 auto inside_cell_center_local = ref_el_inside.position(0,0); 430 auto inside_cell_center_global = geo_inside.center(); 431 432 // distance of cell center to boundary 433 auto d = geo.global(face_local); 434 d -= inside_cell_center_global; 435 auto distance = d.two_norm(); 436 437 // absolute permeability 438 auto k_abs_inside = tp.k_abs(cell_inside,inside_cell_center_local); 439 440 // gravity times normal 441 auto gn = tp.gravity()*ig.unitOuterNormal(face_local); 442 443 // liquid phase Dirichlet boundary 444 if (bc_l==1) 445 { 446 auto rho_l_inside = tp.rho_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)); 447 auto g_l = tp.g_l(ig.intersection(),face_local,time); 448 auto w_l = (x_s(lfsu_s,liquid)-g_l)/distance + rho_l_inside*gn; 449 auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid)); 450 auto lambda_l_inside = tp.kr_l(cell_inside,inside_cell_center_local,s_l)/ 451 tp.mu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)); 452 auto sigma_l = lambda_l_inside*k_abs_inside; 453 auto nu_l = tp.nu_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,liquid)); 454 r_s.accumulate(lfsv_s, liquid, scale_l * nu_l * sigma_l * w_l * face_volume); 455 } 456 457 // gas phase Dirichlet boundary 458 if (bc_g==1) 459 { 460 auto rho_g_inside = tp.rho_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)); 461 auto g_g = tp.g_g(ig.intersection(),face_local,time); 462 auto w_g = (x_s(lfsu_s,gas)-g_g)/distance + rho_g_inside*gn; 463 auto s_l = tp.s_l(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)-x_s(lfsu_s,liquid)); 464 auto s_g = 1-s_l; 465 auto lambda_g_inside = tp.kr_g(cell_inside,inside_cell_center_local,s_g)/ 466 tp.mu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)); 467 auto sigma_g = lambda_g_inside*k_abs_inside; 468 auto nu_g = tp.nu_g(cell_inside,inside_cell_center_local,x_s(lfsu_s,gas)); 469 r_s.accumulate(lfsv_s, gas, scale_l * nu_g * sigma_g * w_g * face_volume); 470 } 471 } 472 473 // boundary integral independent of ansatz functions 474 template<typename IG, typename LFSV, typename R> lambda_boundary(const IG & ig,const LFSV & lfsv,R & r_s) const475 void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r_s) const 476 { 477 // get geometries 478 auto geo = ig.geometry(); 479 480 // face geometry 481 auto ref_el = referenceElement(geo); 482 auto face_local = ref_el.position(0,0); 483 auto face_volume = geo.volume(); 484 485 // evaluate boundary condition type 486 auto bc_l = tp.bc_l(ig.intersection(),face_local,time); 487 auto bc_g = tp.bc_g(ig.intersection(),face_local,time); 488 if (bc_l!=0 && bc_g!=0) return; // no Neumann boundary conditions 489 490 // liquid phase Neumann boundary 491 if (bc_l==0) 492 { 493 auto j_l = tp.j_l(ig.intersection(),face_local,time); 494 r_s.accumulate(lfsv, liquid, scale_l * j_l * face_volume); 495 } 496 497 // gas phase Neumann boundary 498 if (bc_g==0) 499 { 500 auto j_g = tp.j_g(ig.intersection(),face_local,time); 501 r_s.accumulate(lfsv, gas, scale_g * j_g * face_volume); 502 } 503 } 504 505 //! set time for subsequent evaluation setTime(typename TP::Traits::RangeFieldType t)506 void setTime (typename TP::Traits::RangeFieldType t) 507 { 508 time = t; 509 } 510 511 private: 512 const TP& tp; // two phase parameter class 513 typename TP::Traits::RangeFieldType time; 514 Real scale_l, scale_g; 515 516 template<typename T> aavg(T a,T b) const517 T aavg (T a, T b) const 518 { 519 return 0.5*(a+b); 520 } 521 522 template<typename T> havg(T a,T b) const523 T havg (T a, T b) const 524 { 525 T eps = 1E-30; 526 return 2.0/(1.0/(a+eps) + 1.0/(b+eps)); 527 } 528 }; 529 530 531 /** a local operator for the storage operator 532 * 533 * \f{align*}{ 534 \int_\Omega c(x) uv dx 535 * \f} 536 */ 537 template<class TP> 538 class TwoPhaseOnePointTemporalOperator 539 : public NumericalJacobianVolume<TwoPhaseOnePointTemporalOperator<TP> >, 540 public NumericalJacobianApplyVolume<TwoPhaseOnePointTemporalOperator<TP> >, 541 public FullVolumePattern, 542 public LocalOperatorDefaultFlags, 543 public InstationaryLocalOperatorDefaultMethods<typename TP::Traits::RangeFieldType> 544 { 545 enum { dim = TP::Traits::GridViewType::dimension }; 546 enum { liquid = 0 }; 547 enum { gas = 1 }; 548 549 using Real = typename TP::Traits::RangeFieldType; 550 551 public: 552 // pattern assembly flags 553 enum { doPatternVolume = true }; 554 555 // residual assembly flags 556 enum { doAlphaVolume = true }; 557 TwoPhaseOnePointTemporalOperator(TP & tp_,Real scale_l_=1.0,Real scale_g_=1.0)558 TwoPhaseOnePointTemporalOperator (TP& tp_, Real scale_l_=1.0, Real scale_g_=1.0) 559 : tp(tp_), scale_l(scale_l_), scale_g(scale_g_) 560 { 561 } 562 563 //! set time for subsequent evaluation setTime(typename TP::Traits::RangeFieldType t)564 void setTime (typename TP::Traits::RangeFieldType t) 565 { 566 time = t; 567 } 568 569 // volume integral depending on test and ansatz functions 570 template<typename EG, typename LFSU, typename X, typename LFSV, typename R> alpha_volume(const EG & eg,const LFSU & lfsu,const X & x,const LFSV & lfsv,R & r) const571 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const 572 { 573 // Reference to cell 574 const auto& cell = eg.entity(); 575 576 // get geometry 577 auto geo = eg.geometry(); 578 579 // cell geometry 580 auto ref_el = referenceElement(geo); 581 auto cell_center_local = ref_el.position(0,0); 582 auto cell_volume = geo.volume(); 583 584 auto phi = tp.phi(cell,cell_center_local); 585 auto s_l = tp.s_l(cell,cell_center_local,x(lfsu,gas)-x(lfsu,liquid)); 586 587 r.accumulate(lfsv, liquid, scale_l * phi * s_l * tp.nu_l(cell,cell_center_local,x(lfsu,liquid)) * cell_volume); 588 r.accumulate(lfsv, gas, scale_g * phi * (1-s_l) * tp.nu_g(cell,cell_center_local,x(lfsu,gas)) * cell_volume); 589 } 590 591 private: 592 TP& tp; 593 typename TP::Traits::RangeFieldType time; 594 Real scale_l, scale_g; 595 }; 596 597 598 /** \brief Provide velocity field for liquid phase 599 600 Uses RT0 interpolation on a cell. 601 602 - T : provides TwoPhaseParameterInterface 603 - PL : P0 function for liquid phase pressure 604 - PG : P0 function for gas phase pressure 605 */ 606 template<typename TP, typename PL, typename PG> 607 class V_l 608 : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType, 609 typename PL::Traits::RangeFieldType, 610 PL::Traits::GridViewType::dimension, 611 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >, 612 V_l<TP,PL,PG> > 613 { 614 // extract useful types 615 using GV = typename PL::Traits::GridViewType; 616 using DF = typename GV::Grid::ctype; 617 using RF = typename PL::Traits::RangeFieldType; 618 using RangeType = typename PL::Traits::RangeType; 619 enum { dim = PL::Traits::GridViewType::dimension }; 620 using Element = typename GV::Traits::template Codim<0>::Entity; 621 using IntersectionIterator = typename GV::IntersectionIterator; 622 using Intersection = typename IntersectionIterator::Intersection; 623 624 const TP& tp; 625 const PL& pl; 626 const PG& pg; 627 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe; 628 typename TP::Traits::RangeFieldType time; 629 630 631 using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType; 632 633 public: 634 using Traits = Dune::PDELab::GridFunctionTraits<GV,RF,dim,Dune::FieldVector<RF,dim> >; 635 636 using BaseT = Dune::PDELab::GridFunctionBase<Traits,V_l<TP,PL,PG> >; 637 V_l(const TP & tp_,const PL & pl_,const PG & pg_)638 V_l (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {} 639 640 // set time where operator is to be evaluated (i.e. end of the time intervall) set_time(typename TP::Traits::RangeFieldType time_)641 void set_time (typename TP::Traits::RangeFieldType time_) 642 { 643 time = time_; 644 } 645 evaluate(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeType & y) const646 inline void evaluate (const typename Traits::ElementType& e, 647 const typename Traits::DomainType& x, 648 typename Traits::RangeType& y) const 649 { 650 // get geometries 651 auto geo = e.geometry(); 652 653 // face geometry 654 auto ref_el = referenceElement(geo); 655 auto face_local = ref_el.position(0,0); 656 auto face_volume = geo.volume(); 657 658 // cell geometry 659 auto inside_cell_center_local = ref_el.position(0,0); 660 auto inside_cell_center_global = geo.global(inside_cell_center_local); 661 662 // absolute permeability 663 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local); 664 665 // pressure evaluation 666 typename PL::Traits::RangeType pl_inside, pg_inside; 667 pl.evaluate(e,inside_cell_center_local,pl_inside); 668 pg.evaluate(e,inside_cell_center_local,pg_inside); 669 670 // density 671 auto nu_l_inside = tp.nu_l(e,inside_cell_center_local,pl_inside); 672 673 // for coefficient computation 674 RF vn[2*dim]; // normal velocities 675 RF coeff[2*dim]; // RT0 coefficient 676 auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear 677 auto determinant = B.determinant(); 678 679 // loop over cell neighbors 680 for (const auto& intersection : intersections(pl.getGridView(),e)) 681 { 682 // set to zero for processor boundary 683 vn[intersection.indexInInside()] = 0.0; 684 685 // get geometry 686 auto geo_intersection = intersection.geometry(); 687 688 // face geometry 689 const auto& face_local = referenceElement(geo_intersection).position(0,0); 690 691 // interior face 692 if (intersection.neighbor()) 693 { 694 // get geometry 695 auto geo_outside = intersection.outside().geometry(); 696 auto ref_el_outside = referenceElement(geo_outside); 697 auto outside_cell_center_local = ref_el_outside.position(0,0); 698 auto outside_cell_center_global = geo_outside.global(outside_cell_center_local); 699 700 // distance of cell centers 701 auto d = outside_cell_center_global; 702 d -= inside_cell_center_global; 703 auto distance = d.two_norm(); 704 705 // absolute permeability 706 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local); 707 708 // gravity times normal 709 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local); 710 711 // pressure evaluation 712 typename PL::Traits::RangeType pl_outside, pg_outside; 713 pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside); 714 pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside); 715 716 // density 717 auto nu_l_outside = tp.nu_l(*(intersection.outside()),outside_cell_center_local,pg_outside); 718 719 // liquid phase calculation 720 auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside); 721 auto rho_l_outside = tp.rho_l(*(intersection.outside()),outside_cell_center_local,pl_outside); 722 auto w_l = (pl_inside-pl_outside)/distance + aavg(rho_l_inside,rho_l_outside)*gn; // determines direction 723 RF pc_upwind, s_l_upwind, s_g_upwind; 724 if (w_l>=0) // upwind capillary pressure on face 725 { 726 pc_upwind = pg_inside-pl_inside; 727 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind); 728 } 729 else 730 { 731 pc_upwind = pg_outside-pl_outside; 732 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind); 733 } 734 s_g_upwind = 1-s_l_upwind; 735 auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l_upwind)/ 736 tp.mu_l(e,inside_cell_center_local,pl_inside); 737 auto lambda_l_outside = tp.kr_l(*(intersection.outside()),outside_cell_center_local,s_l_upwind)/ 738 tp.mu_l(*(intersection.outside()),outside_cell_center_local,pl_outside); 739 auto sigma_l = havg(lambda_l_inside*k_abs_inside,lambda_l_outside*k_abs_outside); 740 auto nu_l = aavg(nu_l_inside,nu_l_outside); 741 742 // set coefficient 743 vn[intersection.indexInInside()] = nu_l * sigma_l * w_l; 744 } 745 746 // boundary face 747 if (intersection.boundary()) 748 { 749 // distance of cell center to boundary 750 auto d = geo_intersection.global(face_local); 751 d -= inside_cell_center_global; 752 auto distance = d.two_norm(); 753 754 // evaluate boundary condition type 755 auto bc_l = tp.bc_l(intersection,face_local,time); 756 757 // liquid phase Dirichlet boundary 758 if (bc_l==1) 759 { 760 auto rho_l_inside = tp.rho_l(e,inside_cell_center_local,pl_inside); 761 auto g_l = tp.g_l(intersection,face_local,time); 762 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local); 763 auto w_l = (pl_inside-g_l)/distance + rho_l_inside*gn; 764 auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside); 765 auto lambda_l_inside = tp.kr_l(e,inside_cell_center_local,s_l)/ 766 tp.mu_l(e,inside_cell_center_local,pl_inside); 767 auto sigma_l = lambda_l_inside*k_abs_inside; 768 vn[intersection.indexInInside()] = nu_l_inside * sigma_l * w_l; 769 } 770 771 // liquid phase Neumann boundary 772 if (bc_l==0) 773 { 774 auto j_l = tp.j_l(intersection,face_local,time); 775 vn[intersection.indexInInside()] = j_l; 776 } 777 } 778 779 // compute coefficient 780 auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element 781 vstar *= vn[intersection.indexInInside()]; 782 Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element 783 if (intersection.indexInInside()%2==0) 784 normalhat[intersection.indexInInside()/2] = -1.0; 785 else 786 normalhat[intersection.indexInInside()/2] = 1.0; 787 Dune::FieldVector<DF,dim> vstarhat(0); 788 B.umtv(vstar,vstarhat); // Piola backward transformation 789 vstarhat *= determinant; 790 coeff[intersection.indexInInside()] = vstarhat*normalhat; 791 } 792 793 // compute velocity on reference element 794 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size()); 795 rt0fe.localBasis().evaluateFunction(x,rt0vectors); 796 typename Traits::RangeType yhat(0); 797 for (unsigned int i=0; i<rt0fe.localBasis().size(); i++) 798 yhat.axpy(coeff[i],rt0vectors[i]); 799 800 // apply Piola transformation 801 B.invert(); 802 y = 0; 803 B.umtv(yhat,y); 804 y /= determinant; 805 806 // std::cout << "vn= " ; 807 // for (int i=0; i<2*dim; i++) std::cout << vn[i] << " "; 808 // std::cout << std::endl; 809 // std::cout << "V_l: x=" << x << " y=" << y << std::endl; 810 } 811 getGridView()812 inline const typename Traits::GridViewType& getGridView () 813 { 814 return pl.getGridView(); 815 } 816 817 private: 818 819 template<typename T> aavg(T a,T b) const820 T aavg (T a, T b) const 821 { 822 return 0.5*(a+b); 823 } 824 825 template<typename T> havg(T a,T b) const826 T havg (T a, T b) const 827 { 828 T eps = 1E-30; 829 return 2.0/(1.0/(a+eps) + 1.0/(b+eps)); 830 } 831 }; 832 833 /** \brief Provide velocity field for gas phase 834 835 Uses RT0 interpolation on a cell. 836 837 - T : provides TwoPhaseParameterInterface 838 - PL : P0 function for liquid phase pressure 839 - PG : P0 function for gas phase pressure 840 */ 841 template<typename TP, typename PL, typename PG> 842 class V_g 843 : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename PL::Traits::GridViewType, 844 typename PL::Traits::RangeFieldType, 845 PL::Traits::GridViewType::dimension, 846 Dune::FieldVector<typename PL::Traits::RangeFieldType,PL::Traits::GridViewType::dimension> >, 847 V_g<TP,PL,PG> > 848 { 849 // extract useful types 850 using GV = typename PL::Traits::GridViewType; 851 using DF = typename GV::Grid::ctype; 852 using RF = typename PL::Traits::RangeFieldType; 853 using RangeType = typename PL::Traits::RangeType; 854 enum { dim = PL::Traits::GridViewType::dimension }; 855 using Element = typename GV::Traits::template Codim<0>::Entity; 856 using IntersectionIterator = typename GV::IntersectionIterator; 857 using Intersection = typename IntersectionIterator::Intersection; 858 859 const TP& tp; 860 const PL& pl; 861 const PG& pg; 862 Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0> rt0fe; 863 typename TP::Traits::RangeFieldType time; 864 865 866 using RT0RangeType = typename Dune::RaviartThomasCubeLocalFiniteElement<DF,RF,dim,0>::Traits::LocalBasisType::Traits::RangeType; 867 868 public: 869 using Traits = Dune::PDELab::GridFunctionTraits<GV,RF,dim,Dune::FieldVector<RF,dim> >; 870 871 using BaseT = Dune::PDELab::GridFunctionBase<Traits,V_g<TP,PL,PG> >; 872 V_g(const TP & tp_,const PL & pl_,const PG & pg_)873 V_g (const TP& tp_, const PL& pl_, const PG& pg_) : tp(tp_), pl(pl_), pg(pg_), time(0) {} 874 875 // set time where operator is to be evaluated (i.e. end of the time intervall) set_time(typename TP::Traits::RangeFieldType time_)876 void set_time (typename TP::Traits::RangeFieldType time_) 877 { 878 time = time_; 879 } 880 evaluate(const typename Traits::ElementType & e,const typename Traits::DomainType & x,typename Traits::RangeType & y) const881 inline void evaluate (const typename Traits::ElementType& e, 882 const typename Traits::DomainType& x, 883 typename Traits::RangeType& y) const 884 { 885 // get geometries 886 auto geo = e.geometry(); 887 888 // face geometry 889 auto ref_el = referenceElement(geo); 890 auto face_local = ref_el.position(0,0); 891 auto face_volume = geo.volume(); 892 893 // cell geometry 894 auto inside_cell_center_local = ref_el.position(0,0); 895 auto inside_cell_center_global = geo.global(inside_cell_center_local); 896 897 // absolute permeability 898 auto k_abs_inside = tp.k_abs(e,inside_cell_center_local); 899 900 // pressure evaluation 901 typename PL::Traits::RangeType pl_inside, pg_inside; 902 pl.evaluate(e,inside_cell_center_local,pl_inside); 903 pg.evaluate(e,inside_cell_center_local,pg_inside); 904 905 // density evaluation 906 auto nu_g_inside = tp.nu_g(e,inside_cell_center_local,pg_inside); 907 908 // for coefficient computation 909 RF vn[2*dim]; // normal velocities 910 RF coeff[2*dim]; // RT0 coefficient 911 auto B = geo.jacobianInverseTransposed(x); // the transformation. Assume it is linear 912 auto determinant = B.determinant(); 913 914 // loop over cell neighbors 915 for (const auto& intersection : intersections(pl.getGridView(),e)) 916 { 917 // set to zero for processor boundary 918 vn[intersection.indexInInside()] = 0.0; 919 920 // get geometry 921 auto geo_intersection = intersection.geometry(); 922 923 // face geometry 924 const auto& face_local = referenceElement(geo_intersection).position(0,0); 925 926 // interior face 927 if (intersection.neighbor()) 928 { 929 // get geometry 930 auto geo_outside = intersection.outside().geometry(); 931 auto ref_el_outside = referenceElement(geo_outside); 932 auto outside_cell_center_local = ref_el_outside.position(0,0); 933 auto outside_cell_center_global = geo_outside.global(outside_cell_center_local); 934 935 // distance of cell centers 936 auto d = outside_cell_center_global; 937 d -= inside_cell_center_global; 938 auto distance = d.two_norm(); 939 940 // absolute permeability 941 auto k_abs_outside = tp.k_abs(*(intersection.outside()),outside_cell_center_local); 942 943 // gravity times normal 944 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local); 945 946 // pressure evaluation 947 typename PL::Traits::RangeType pl_outside, pg_outside; 948 pl.evaluate(*(intersection.outside()),outside_cell_center_local,pl_outside); 949 pg.evaluate(*(intersection.outside()),outside_cell_center_local,pg_outside); 950 951 // needed for both phases 952 RF pc_upwind, s_l_upwind, s_g_upwind; 953 954 // gas phase calculation 955 auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside); 956 auto rho_g_outside = tp.rho_g(e,outside_cell_center_local,pg_outside); 957 auto w_g = (pg_inside-pg_outside)/distance + aavg(rho_g_inside,rho_g_outside)*gn; // determines direction 958 if (w_g>=0) // upwind capillary pressure on face 959 { 960 pc_upwind = pg_inside-pl_inside; 961 s_l_upwind = tp.s_l(e,inside_cell_center_local,pc_upwind); 962 } 963 else 964 { 965 pc_upwind = pg_outside-pl_outside; 966 s_l_upwind = tp.s_l(*(intersection.outside()),outside_cell_center_local,pc_upwind); 967 } 968 s_g_upwind = 1-s_l_upwind; 969 auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g_upwind)/ 970 tp.mu_g(e,inside_cell_center_local,pg_inside); 971 auto lambda_g_outside = tp.kr_g(*(intersection.outside()),outside_cell_center_local,s_g_upwind)/ 972 tp.mu_g(*(intersection.outside()),outside_cell_center_local,pg_outside); 973 auto sigma_g = havg(lambda_g_inside*k_abs_inside,lambda_g_outside*k_abs_outside); 974 975 auto nu_g_outside = tp.nu_g(*(intersection.outside()),outside_cell_center_local,pg_outside); 976 977 // set coefficient 978 vn[intersection.indexInInside()] = aavg(nu_g_inside,nu_g_outside) * sigma_g * w_g; 979 } 980 981 // boundary face 982 if (intersection.boundary()) 983 { 984 // distance of cell center to boundary 985 auto d = geo_intersection.global(face_local); 986 d -= inside_cell_center_global; 987 auto distance = d.two_norm(); 988 989 // evaluate boundary condition type 990 auto bc_g = tp.bc_g(intersection,face_local,time); 991 992 // gas phase Dirichlet boundary 993 if (bc_g==1) 994 { 995 auto rho_g_inside = tp.rho_g(e,inside_cell_center_local,pg_inside); 996 auto g_g = tp.g_g(intersection,face_local,time); 997 auto gn = tp.gravity()*intersection.unitOuterNormal(face_local); 998 auto w_g = (pg_inside-g_g)/distance + rho_g_inside*gn; 999 auto s_l = tp.s_l(e,inside_cell_center_local,pg_inside-pl_inside); 1000 auto s_g = 1-s_l; 1001 auto lambda_g_inside = tp.kr_g(e,inside_cell_center_local,s_g)/ 1002 tp.mu_g(e,inside_cell_center_local,pg_inside); 1003 auto sigma_g = lambda_g_inside*k_abs_inside; 1004 1005 vn[intersection.indexInInside()] = nu_g_inside * sigma_g * w_g; 1006 } 1007 1008 // gas phase Neumann boundary 1009 if (bc_g==0) 1010 { 1011 auto j_g = tp.j_g(intersection,face_local,time); 1012 vn[intersection.indexInInside()] = j_g; /* /nu_g_inside*/; 1013 } 1014 } 1015 1016 // compute coefficient 1017 auto vstar = intersection.unitOuterNormal(face_local); // normal on tranformef element 1018 vstar *= vn[intersection.indexInInside()]; 1019 Dune::FieldVector<RF,dim> normalhat(0); // normal on reference element 1020 if (intersection.indexInInside()%2==0) 1021 normalhat[intersection.indexInInside()/2] = -1.0; 1022 else 1023 normalhat[intersection.indexInInside()/2] = 1.0; 1024 Dune::FieldVector<DF,dim> vstarhat(0); 1025 B.umtv(vstar,vstarhat); // Piola backward transformation 1026 vstarhat *= determinant; 1027 coeff[intersection.indexInInside()] = vstarhat*normalhat; 1028 } 1029 1030 // compute velocity on reference element 1031 std::vector<RT0RangeType> rt0vectors(rt0fe.localBasis().size()); 1032 rt0fe.localBasis().evaluateFunction(x,rt0vectors); 1033 typename Traits::RangeType yhat(0); 1034 for (unsigned int i=0; i<rt0fe.localBasis().size(); i++) 1035 yhat.axpy(coeff[i],rt0vectors[i]); 1036 1037 // apply Piola transformation 1038 B.invert(); 1039 y = 0; 1040 B.umtv(yhat,y); 1041 y /= determinant; 1042 } 1043 getGridView()1044 inline const typename Traits::GridViewType& getGridView () 1045 { 1046 return pl.getGridView(); 1047 } 1048 1049 private: 1050 1051 template<typename T> aavg(T a,T b) const1052 T aavg (T a, T b) const 1053 { 1054 return 0.5*(a+b); 1055 } 1056 1057 template<typename T> havg(T a,T b) const1058 T havg (T a, T b) const 1059 { 1060 T eps = 1E-30; 1061 return 2.0/(1.0/(a+eps) + 1.0/(b+eps)); 1062 } 1063 }; 1064 1065 1066 } 1067 } 1068 1069 #endif // DUNE_PDELAB_LOCALOPERATOR_TWOPHASECCFV_HH 1070