From 0558184439efbc6aeee1cbe6f808651fde11167c Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Thu, 12 Dec 2013 20:29:40 +0800 Subject: [PATCH] modify the function for the polymer term. --- .../.IncompPropsAdFromDeck.cpp.swp | Bin 0 -> 16384 bytes .../FullyImplicitTwophasePolymerSolver.cpp | 93 +++++++++--------- opm/polymer/fullyimplicit/PolymerPropsAd.cpp | 10 ++ opm/polymer/fullyimplicit/PolymerPropsAd.hpp | 1 + 4 files changed, 55 insertions(+), 49 deletions(-) create mode 100644 opm/polymer/fullyimplicit/.IncompPropsAdFromDeck.cpp.swp diff --git a/opm/polymer/fullyimplicit/.IncompPropsAdFromDeck.cpp.swp b/opm/polymer/fullyimplicit/.IncompPropsAdFromDeck.cpp.swp new file mode 100644 index 0000000000000000000000000000000000000000..5a92cb6fba3072c52aea393cf7b815ab985ff255 GIT binary patch literal 16384 zcmeI3TWB3c7{|x^wY9WZ@lqV4hI7*NoFppR%>`4_+G0aPX;ZMo!)A9*+-7gi?n#KD zim2^VQ4qlgUqn!R@CMcw`ywd5ct-^BMT_92555$Hir{}{cF)c^Jts{SM05`PvbULU zzWL8(cjw#84ovTwq<2^DB>1f(Syt-4!}f5Kr+aPR89GJPIb5GhcDy z_?#u01V#a)K%u}Ia>vl%fH=Oj=N7tQ{r=*T*)s|l1&jhl0i%FXz$jo8FbWt2{__f? z@rC4Z)bhftrVH8jg0}5N+41^pyV!O;J2V@kfKk9GU=%P47zK<1MggOMQNSo*6fg=H z1^$H!*e)Tx@aqP#0D$;^{{H{_HH4f1Pl6dR2sVM6z*=xMII)(H0L+2=K?Pg_K3+q} zA+Qm=eKjGkf|o!5=70?z02{$N@XJ+%`~cnoZ-8UqC^!OkgB{>@a5MP+NLVg2ZgU`Wd;8XBCcm^;~0}p{2unAlTe!ZNKpTSSyNAMna z4p?9-xCh(@PG3gI2jD1p0?dG|U=uiXDar)L!HZxYxC2}ZPF@24!7JbxI0}w{r@>Po z1S8-|@X5u5d<0$wj{^q`g0k+4=A{9~$3}=bx2o2+5x_PrQP^Qbx0^q__%4u?NXa+f!X>YRFk3SJTTv$$*V&<}X z<_L21QI)Mo!kHg49wZgc62=!;DU7vRwAFcx=E(~KMI{Rz;7}4(aJ)^*QqnEjX)Hcs zsi8nwaw+VmNzL~CifgBKsoYN;=KD$WXr4EB-R8b#P!>UPNf-*_z83%0#r&FCI?(iSmiPE`50L3j!z1|{{e3|5eDyuK7@9>V&w6li>Hphvr_ZzB ziu4_NEgC!PviT&t-J4^fWsP%gFO9fSW7cQ=G&?!56xxC3bUof1El%2pB0ffCJF9le z4R;w{({y@O=`irbtUugsgvASKt&%i&9iwN0g^8DX#hkM&$DPchOJFq0PqCDy{ns z&lQgmOtw+Ci4qSWwTifxGkeza|5_bevukb?lWsmQtKuU5SLKYGS#veTbl#;LXd`b| zt9KJg=Qk33@KA7&IA&I$DbrHB-HyOAmOgJ$4QrvJF+)h!KEiQ{$h~bgadLwkiDJq_^XsL%HeQ{^fd!95@>Cx^yc@dQM z2>&2ck_wy1HOn}S>griFL47Ia61@mHwP@B!VT;l|Q|I4@@!dME?~mhLVKDecy(H8FP`3l=hEi=So2H2-xmqI&`V1$KC|wp zbKGCk|FN0!!`X*9t^7&V4oc?#HVHZ#ZaM=H3r)dZkaQYjo^3 zmeuieNtxU3_+FeaVQI?7r|Z0z%9@g8f;$yUV%!Yqo;biqvR>%;4VTdoFM#*^4d4vyoCYT| zo8Jk$=4TWz3K#{90!9I&fKk9GU=%P47zO?j1rF0{we`JDyml4~3HhQ<77S8rW35T@ z_tfTkc*}x278ornUu9@newR@%M0G@{OkaxcY`t(PFHdYm-h;CumNY{*kBVT=V9ia2 zYOZI`(LTh$u9|&ZXvm8Zba+A#{fHbM10BRTEwq`I-gx%y>eRUu= zA?rKY<_T68V&PIr1gb>*0W^ZJz1G5pdJ7i|SSk&iP$51UBit+kS*lLrco?lm_c8}T zGK=nrc*^aNMqJc2oTGL8GdOLFsi}2=+*z@Og>giy0<20)udC3nLi7}3zi_K4J$G%W zRmi4@Dy0@n8EOrt47KZbb{p&MJ1Q2d5nPpVLDjH>of}q2jP*KgDTHT_fNTW BYpVbN literal 0 HcmV?d00001 diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index a101c5512..a68110425 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -226,6 +226,19 @@ typedef Eigen::Array c.value()(i)) ? cmax(i) : c.value()(i); + } + + return ADB::constant(cmax); + } @@ -243,9 +256,12 @@ typedef Eigen::Array kr = computeRelPerm(state); - for (int phase = 0; phase < fluid_.numPhases(); ++phase) { - const ADB mflux = computeMassFlux(phase, trans, kr, state); - ADB source = accumSource(phase, kr, src); + + const ADB cmax = computeCmax(state.concentration); + const ADB krw_eff = polymer_props_.effectiveRelPerm(c, cmax, kr[0], state.saturation[0]); + + const std::vector mflux = computeMassFlux(trans, kr, state); + const std::vector source = accumSource(phase, kr, src); residual_[phase] = pvdt*(state.saturation[phase] - old_state.saturation[phase]) + ops_.div*mflux - source; @@ -255,17 +271,19 @@ typedef Eigen::Array& kr, - const std::vector& src) const + std::vector + FullyImplicitTwophasePolymerSolver::accumSource(const std::vector& kr, + const std::vector& src, + const std::vector& polymer_inflow_c, + const SolutionState& state) const { //extract the source to out and in source. std::vector outsrc; @@ -286,6 +304,7 @@ typedef Eigen::Array(& src[0], grid_.number_of_cells); const V outSrc = Eigen::Map(& outsrc[0], grid_.number_of_cells); const V inSrc = Eigen::Map(& insrc[0], grid_.number_of_cells); + const V polyin = Eigen::Map(& polymer_inflow_c[0], grid_.number_of_cells); // compute the out-fracflow. ADB f_out = computeFracFlow(phase, kr); // compute the in-fracflow. @@ -336,7 +355,7 @@ typedef Eigen::Array FullyImplicitTwophasePolymerSolver::computeFracFlow(int phase, const std::vector& kr) const { @@ -452,54 +471,30 @@ typedef Eigen::Array& kr , - const SolutionState& state ) const + std::vector + FullyImplicitTwophasePolymerSolver::computeMassFlux(const V& trans, + const ADB& mc, + const std::vector& kr , + const SolutionState& state ) const { -// const ADB tr_mult = transMult(state.pressure); const double* mus = fluid_.viscosity(); - ADB mob = ADB::null(); - if (phase == 0) { - ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus); -// for (int i = 0; i < grid_.number_of_cells; ++i) { -// std::cout << state.concentration.value()(i) << " " << 1./inv_wat_eff_vis.value()(i) << std::endl; -// std::cout << "water effective vis\n"<<1./inv_wat_eff_v1./inv_wat_eff_vis.value()is.value()< mflux(2, ADB::null()); + ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus); + ADB wat_mob = kr[0] * inv_wat_eff_vis; + ADB oil_mob = kr[1] / V::Constant(kr[1].size(), 1, mus[1]); + ADB poly_mob = mc * kr[0] * inv_wat_eff_vis; + + const ADB dp = ops_.ngrad * state.pressure; const ADB head = trans * dp; - UpwindSelector upwind(grid_, ops_, head.value()); - return upwind.select(mob) * head; + mflux.push_back(upwind.select(wat_mob)*head); + mflux.push_back(upwind.select(oil_mob)*head); + mflux.push_back(upwind.select(poly_mob)*head); + return mflux; } - ADB - FullyImplicitTwophasePolymerSolver::computePolymerMassFlux(const V& trans, - const ADB& mc, - const std::vector& kr, - const SolutionState& state) const - - { - const double* mus = fluid_.viscosity(); - ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus); - ADB poly_mob = mc * kr[0] * inv_wat_eff_vis; - - const ADB dp = ops_.ngrad * state.pressure; - const ADB head = trans * dp; - - UpwindSelector upwind(grid_, ops_, head.value()); - - return upwind.select(poly_mob) * head; - - } - double FullyImplicitTwophasePolymerSolver::residualNorm() const diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp index f9ce55952..d1b984ac4 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp @@ -163,6 +163,16 @@ namespace Opm { } */ + double + PolymerPropsAd::rockDensity() const + { + return polymer_props_.rockDensity(); + } + + + + + PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props) : polymer_props_ (polymer_props) { diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp index b964a241d..5c899ab40 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -63,6 +63,7 @@ namespace Opm { const ADB& muWat, const ADB& muPolyEff) const; */ + double rockDensity() const; typedef AutoDiffBlock ADB; typedef ADB::V V; PolymerPropsAd(const PolymerProperties& polymer_props);