From 6d5b90df548214d4993fe132861d2403e63b4e8c Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Tue, 17 Dec 2013 20:28:32 +0800 Subject: [PATCH] modified the simulator more general. --- .../.FullyImplicitTwoPhaseSolver.cpp.swp | Bin 36864 -> 0 bytes .../SimulatorFullyImplicitTwophase.cpp | 279 +++++++++++++++--- .../SimulatorFullyImplicitTwophase.hpp | 9 +- 3 files changed, 251 insertions(+), 37 deletions(-) delete mode 100644 opm/polymer/fullyimplicit/.FullyImplicitTwoPhaseSolver.cpp.swp diff --git a/opm/polymer/fullyimplicit/.FullyImplicitTwoPhaseSolver.cpp.swp b/opm/polymer/fullyimplicit/.FullyImplicitTwoPhaseSolver.cpp.swp deleted file mode 100644 index 776e79ab91a6a8d8336cc313dcba0b0fa42293bb..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 36864 zcmeI533Mb!dB>Z>fB_q0ZWAZ%<-PStD~+TzV5~i0PrQWniaoTzdKpGDEopXVrpG-q zlGe*|B#;0m=1Rae;dmGVCMKNb#s*#>#Ch<5B;4k(frK+aILsCD`|9Yb?iopI8(v-> z)9?K?Gu>5Ref3qno;_hfg;I#}?z{=vdrp2u53R{a7PP)NlKJwp*_st2J8nTBX(*+bO?x znP2-C@^Shh-w*1o!!#4k2%yXOW&vMVr1efsNoH;O^k3 z_YQ*3gLi<-!C`PQm;!$a9tj=@&H(>;uON6cxE8zwyZ|hLHQ=G(EbtKUANQnuZ~|Ni zwu1YEn@fdb=op zOWbdA)3KZ8vHUj|`DXJpyN#JBDaErTMMGyho^q%2^pZbhip*sX4-chTN5a!fd@; zs}_-uiV|uUhDQ9-hI~Secm^7)WC$!%BjIA)YL~QZ9nT}7x;vNDs{YN2NDk^V-CDi6 zrC#4AzYCH|cr)EdB^vUx%_Zf963NO7wazjBLA};&tB9vkA^qwSGbGs*~bP)k$W=Hp zgmz2jWBD?y+L4cf^vhB;YIWwrQOcOr-#xW1+2!&(M(E|*6wVdZ&q~znByp|UK2)tW zAT2k;b)ozi4JQset~EVN5v{G0TGcD9S(41{Lq~GCX5%C*NQNV+qSmRkazC_qcmt0% z_R;&yRqz5A;ktFyKdY@WRZC@%7NMk5;=a%iEcYodizFQ15H3<~b5ZCo3jH))xy_8( zDbNDNcJ0y#Ex1WFgZbr2dRm$ZV$+k8PmYqf&}=bWM#2e_q|`L4A_9@C$ez!r*{(4I zI^J9=MI)_AMp(?=eRD{C)QIx5m-?+!A5VF)k}#hcDSFRDV>or>J&>wl-xl`U+k1I!=tovoyehk{R$}&B`P)!r%mWo@FqxXRUD>EO>)v6uCjYt4J!PwTC-B`R-nm z{%m6e%E+tDDE5{uhJ^qk9P7JMST>VPvPDFclSluc^8QnPKid8CZ!!Ji!l;f$=qucn zGpY{NuRPU$N=mFF#yYLp@>nM+H`^uodeZeF+yZU}lI~Mra1A6dkib9!0|^WyFp$7N0s{#QBruS`Kmr2^3?wj+!2h8H za2#aJQi>I^r#Ty-u`FUv-JtgWYf1A(38ZKLKWkWjFMIq;z)`RXJP14x+=`$7%iv4k zdhk-P51bDk2_6K#gJ1tepaG`99&j!Qz;*cXUk+XdO5iDA6W9pO1*e0X@!@|R+z2iM z&jLHa+2H5+@oxg}25$y008a%aupUf+hk^%yGr@hq8Q?zPbZ`o|1z*3!09+1c!DGOg z;NIXCeE-*gH-ZjW2Y!fAyWFW3VQrk$4mur^d@xi#Ja z@%;7m5R8q5+r()uF09ZRpztZo&PGYptVE&q(uZPxR(ADbU<*-Mm2q7?>CN4#mGQbe zp9L9qijA*+OdHQtnL;!a2S<#F9q&a-rX7}>RZV%e9xp=37OZQBIW}ySB&TvPZ3^4n znRe8{{w}smTNOOG6rkNmo{I%mggRz7ER(IcS>>G<87?$iVz$H}=ao&y9ZQ(zkk&#u zsg(&rXe;6PRZt9B4%$utn(@^3tyxp?3F1#K;>E{ZW~|z~%dL&M^qaLGx>OXOMrjr% zM({L&)=&{caZ+0<^ty-;)7o%mUi?W~GJN;k3c85tP+1$rxyn>5J_Zz$sIbq4!EFawc91b z^dTh@aWJitU^HyU-J}vZaR^blWZ3FwP$^d`-NrtZvY<$hgh!K+u$@@9l+F`|wb`)H zj63qA7&Tg*V}+r-^H!3sQ-K*Mx`p8eUi_ur;#&*V&alxD?J}NpC01Ln(J-!8OM0C2 zJk9zpR}oNPIhRk}(WK^zD=KchlDm=ZPWZHySEJC#K-8Ew&{F@Ok`pBG0qUQhnD&;s&G*IA94Ai2@Cir+>p_3cRHjZ$R{ zswhF#mn6)i?Ga*B*v=SA?bG3~QfJkr6Fj0ul5$pQWdTtixboYn0kd3GBf*<3q-bIo zn$c1Rxk~^84PD@C+eR@Udqqmfsod#ItlO%6C^Wg(tqZ-l1Nz0*xYLOnqw{fHT1wBn zZqj$WInv$`oe-teW^pS|?bVLeO~Ovvi@7Iik>KSYSNU_P)ZE#YSy8IE4ChjqR%#M8 z$hXN!)!JnwLj=d3CgZWmOg+Kam}!+eUhi}rPxlYM&mq{&Ox$F!#g88@&WOu=N`+T8 zWv;gsMSA;EW#-Kzx1|NtJZKvh?cx)kN}D6So9UyvsR>lcHoptE5g- zmNb=}U0gDRyQ^7~#cZaFTA3`*c-z-RmT%Ra+9g)OEaU10skv=Moy08R z<^q*1x^8T$EZt4Dg{7RSX(oL-n>seU`Z~&9r^>fazt)lEW)It}B!yi&wo_(35o8&4 zDOYCGrt4vnZmXh{MHpJeiA4NiHoaMmZMu74g{7N84F{>lQf6N5mui*8;PEt%En$^U z%|OqTq9>jewYnWuYAQ1GLw?d3QTzX`=>4nA{(q0zP=1NM{*6H5{`Z0r@WYefHnq09Sym;2z+U?7=1O{>eb%>wnJP`+eXtZ~)vFypcWklfeDK+u3Iy0q26V z!5@K#g731|{uH!)g^rzd@xi8LRO3WT1UINg#}wmPNg5_&g!l88pV*3YN_BGlq` zq|~f^uJ^fd8(qk*g`8ohoXkZXZ+onH5+;1Yx*VoqxoPsa!W@HHcIB5`d!}}m_FuSn z&p{npmvoyQLLEzyL^+LwYqph}Vcf)~Fwc&jrNh<>r$Pz#l64$$s8}G@cnU@co(Q8l zRqGv>9(#m%!u4WhbhGC7kbJFlF&8vD-m*o=$=*aIq|i(T4?6uq5s#QTL_lKv*mq#>&Mn(^l=g1f zzeAT@q~|C~{K>Y8?Em3Qc&r37mI?}ZN#upl1-a$ zvSA*}Zw0mZq>H6FAfvt$+>ZgJ&^$`0X|<(JY@qk4`$e&hR}a~jH1Y-;uLIlauli{V zEY+Gv)IpY9C(E%ZtifdxSGyI?# z1TEGfF*~G&p`l_~;voq;9@fI@#2xl$l$vcOq;;Sz-9$YPc`i|VHDl#s*`jxy+UTY0 zF2w}Cz%;XjF-yHHC$K5I-kMtW=)GyS|5S$QZSPOuo(gHZI0*Im=Gly(UxNtIKU(Oyosl9HV2fvS}(YRSpX>=sA| zu`dp@3d6Rec3x=O+-An(M>owrHAFZdW_I7YQ*>fAtV?_Toah-X*f$qZL?j zr0pI0IbhdVKG7SVHF?ySmZfqYHAh6JmB*0qM;|rA45fsk+FOok2-`y(>!im@2)!E}mXI1Hvsh3aE2|JY zoVB#5Nn6H2sVIvDHQjFtr}pmIeqh^%4X(crlZHIOO=*gCKeyEuw&7K*xy|l7zlqoo5O3zFxR*UzYPkgOQ z;hAC`>(X|XeM`=C{*#y#^zer|S(LqVyPS5x5Udx!Ib)O8WzzQRqrTVo)!v%Pbd#;8 zEsE{g(56XJ7fG*Y#?{{%&)8I^YrVX)!Q)$WYqm|6Rjo1Zk8-x%vt33S zUaM%g_qWz@eUwXjr>m4}RN>ds`u(?&Uz#2(l)22zQ@P!;$x2G>qO3|D+bt*6usm1h%x5QEtk#;W)Q@%MWj%%)DXf=y zm>%?Ui!RY-b3Rqwh4swI89x{?=QYJ@DDnRfVX1$U>{;>uAJH4e;2rGqF9y@#Vc@sy z^TqFfE$DzJf%V|tK>Yo3=Kq!8x!|eb0!oy;Ir5O-VR<1o(m3u zhl4u-@$vsN_z3tg_z<`OTn|17t^*A)2Q~of=l?DCf?t4J!Oy|Zz)!(X!1uvh!JEJ{ zz)^4pI2F7Vo5EATH24MfgiisnFFX#P@#_c81S^Yr)fjoZr6yJP7~#9F5{fz>Ge40d!~wPe9d z9(G&I=TX`#-Ll6F8ZjKQo5rBMDXe0)RHq>6kAY&CKWPG#P*s+QI*HRZmT7rNz0AB* zWy`p)Rc>lC`$$U0Uk7GLw-=YEtnITpl;C6ti=8kfzDbVcrPWEf?a!;v9%8g%1uZ1h z4&zxwL?im~3HW%Ua%Dc06h$XBTTI-HkV$6cA#R6^g0g8Z_~wjfaB^h)R5_cgqNdlM;&NNr#HUA*Vkqz7qd)3lFJ-y8$vEq(p1a#U z^w|o-HTrZNo-+AQDk4=rc(DfY2(#r{ol)FlEpr6cWK5G=>o{lPS{+9aG$tn%$;Xu| z7!1nIeaG6Js8JZ24A%sXz$&)9G^I>u=Oj@JQ-S8MXU6ntn(un^ilNFbbj^lfM0Abb zXn0+?4*|I)VKl{<^TNrlW;I%x>Tu?FMBlA7n|l<)6`U4j%r9daTq=8RbCG!vz$e9Vt?DM z;)m>MS-aGA2X=*d2DPJFjgr~|!4td#M*M8)xK6pzg~&QB%>)5_UiEUT9nn$R1^P{g zm&G6%LQH5m?>!$nI4HVAl;o)*SVzvLhMFa?jCoS*7Rx-EkT=#S!YA~Hl##PE*W=#+Sfd2eY_&`#`v)lgXK=* z{EW^Mu_!_NN*MNBMS2zBY>APJd6U2Uk_xQMVk6~ut;tWvdsa89G-}#Zp8CE)UEa_! zz5{7ev!7&*@3*SR+tE}*aW$G5ViP{GX+_Gh&N9J>OU1L)nxnFjmhi}Q)uoz7OCmKa zrLxr$U6vk#>q9V*D#@jX;Ix(}ZQipVRt{b^FiAyu-dC%)rDxh}If6XCtrXJ$v;Hie zEV*n?$^H_7sH%!Ph~6EiW5uN!uK9r`lEUk`JPbGY zcH4e*`j^(jRBv(8simBa*>d#mSYldegyKvz$BQ!X(U1`JD7T_!Rfn%p0nbI}o5hVm z*ch2hHD+shlL16r&+e06H`mHZbj4hnu`wcI#Fe7Po9R}p!%3%nXy_MxQ;)ynfynSS zU3w7oVRx+xHuJRj9VCnLB%p}!O11xIi>~&Qvj4a5*Z&ZE{u{vc;ML#+_!Dq1@D29- zuLBkE1TYMK4Ur#%e*ss57lO;dA`m;kDd3~*`)>g62k!@8gx2CoJ3&cDZjUEpl+ zBWwXb1RndAHvq96yb@diUJ5P;&jZf_^R*z121TnjD+ zr-SdX=a>Ee(?9|q2TlWDWAA?)OoId9e&9?X@ACU1_#F5w_zd_Ykn;dHfaiff2MzFO za8K|#V-t9Xu?sZ$jKDeIKHyYv3b+&a4R(QFfS-c@0@s1JfTw}IU;u<;t3d!LReS0^5rM`$gt_ZC2%5)^W4Vny~w7a`u+WgJ9jg zERR5{u1xhHG*~SrLFL9#QqHCtx+d7!^NTL#J!Cz*OsL*?D zr~l1{OL>l*u1F;+@9Ze4hj_&G0f|C)ZX}#n_pF~nvH3*cDr#s%MFJ$#gdMq>vab)!ZH#6Zfm#Qu-xd(q8J-j}`yb zkJE;e)l+IwvsqF5EX*db)wh>v`Rcx~VoNIFW@3>&-8n6d(+|C?$i%g%?0A%BmiIal zYvT42%197hl5RC+dnCs%ip`klY2UIVg{sSSANiRlT^ODjItgl84w@i;XJpol#Hc)z z3rm(HVh(aFt8J5s<2SQxAgMB!@`9Uiplqy%ovYnWRZZpo#3eIKeo_uBc|Gz~y-jTC zNwlfu{vst+uRS$2A)bOS^1Zh~yRF4qvuZ;wyseFdLd-|ys<#0;pwjj@J~avJmX!#1 z#vD0qi!oMI*;0P0hExigKy^%t>b1U#>Q!%gQ!(>RiCQ6@vKFOLl9+A9bbo05VICFG zo8^T$`IgCYsdLMD!`5+PvTf4A`4tYGB+%fl>+6}saGg<%YTsae_tqTTQ<7Eg!Y>p& zkAuiXuIw1Nw7{#&?66m^BlqgS3jCeN7H++>EU2g*_0~U87C2Vq6=?l+hIuPmpbF*B zba<&zNzUOW1uFkFn<{-tbI?5?el)Q@ud$3lHN9^PqiPClbkiBoiPP(5cZ{CT-CDz* zN8FXD8_>1^Oi`|xE*5I-iiGVE{j(uR33@%XlazQJn5@cGXWh;B+&Mba*m}%bv7XeY z<_jyKJG`L`t*V=cebm-!R5ii8;*q7TTywaKLwYra-Wo=;dzeu}W}5+@@cd$>rmk?Rqv0Dsw*Z(wElfiqR8lYvTBiCf1C zYg@|dpgzzX5(MvT`Q{9O*>Lu#E_qnSFLJ=N&m4td*(&?jMT4|#gkC|$)n zMPr5*_YfI9(+Zbw^;(0Sdv&PTZE_}S$5JaP=(l2KmQR*I)#Z@hSet7WxFs_s^T8S5+wALa0XKm! zgD(NG3w!~59$W)vz*#`v5AZGa_7Wo??*(Xq7|1*ScY{ZP)4(U#+y53Bf#(G5;e7=md_}3U8dx5 zs(P72g}SGrZW#c~O`ZvYdp&37?10wk``dTzVJ@Y9?GLkL^K^rDber>EW`gK@{ELWr zH|kV<%pxj3*bgb)oUMdA>owu;3q>?N+0_C;^zU*cv;n_ID5?6_D-5+d)?ax^aQ_hH zUhSr%cY7q7(Qdzw!)8v&6`ho`>Z?a0%l8ZEMLl)+ou53jF)tp+neMNf7Ox-OEs&;` zOAd2-=NgyqnEuj-Rpe%sQ(TB6nbQ!u31ykdqyJwA)fDRl$L^CQ=&}dZL}<5PVfYdMM|#&L5B3Nmxfz#aGMfU)v&)J^TMTZnn*)49`foCfvAD+%WuI!>71ZIM>m8 zcuFriIDDm2-cjn|t~o4bgg2gZadVAZ`iS!g+t2r76*-O0!K%MyaUJ?)z z4#x`!l=~=`$B!+N;5L=gn-PT*&4t90zyUpEL1HyBK|QXW7U$_)nWj>86zLH(vxE!J z_uxwyH@4uO=Py~Zv|PYvwN%pRX= zBQj&%loZ(p9haxeHtB6`PEFgCuVfANiAOj&W(%jsN%4;;jz^=9S|*WnQb|O}U^F_h zKA$S1s9b*|gDNW3=C)I%LeiW}htg3=-CCfSv9OA_HccYZ@-EdFVv*vXJ#qBR9s$bi z487#?6=P-^@FgYKAAwD$XeyHqtEGnTp&)gaf40Ku8udmWWr;NRv>hR(p^j@%1wY7A z6{cgXlsRzNm}3afGB!7FfN4?4i4%EYo6tEq|A$lpT|dQZ{&{PMu?OPX8uW+R!aDfanXAgwp^Se;fX?;(7%?M$c!=zY3-A+ zqp95Xr6fYh>3aTG78d-@&LnQ^*Z%^H|Ey&G7n#}F+J7DvYVU?8bzZF1liB{4vNW^3 zp#l(m)$EwavLwIcvTV5E3pgyGfNS&zqHjG$bOt13{zY0g3+csF)SE>u?+KDXtE&IE h#{ZbX48-eGuyfaZ4V&%*)LA7sps!*m&;3i9{{;vV@NNJA diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp index 377a42265..be908b0a4 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp @@ -50,18 +50,23 @@ namespace Opm { + + class SimulatorFullyImplicitTwophase::Impl { public: Impl(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const IncompPropsAdInterface& props, + WellsManager& wells_manager, LinearSolverInterface& linsolver, - std::vector& src); + std::vector& src, + const double* gravity); SimulatorReport run(SimulatorTimer& timer, TwophaseState& state, - std::vector& src); + std::vector& src, + WellState& well_state); private: @@ -71,9 +76,15 @@ namespace Opm std::string output_dir_; int output_interval_; // Parameters for well control + bool check_well_controls_; + int max_well_control_iterations_; // Observed objects. const UnstructuredGrid& grid_; const IncompPropsAdInterface& props_; + WellsManager& wells_manager_; + const Wells* wells_; + const std::vector& src_; + const double* gravity_; // Solvers FullyImplicitTwoPhaseSolver solver_; // Misc. data @@ -86,10 +97,12 @@ namespace Opm SimulatorFullyImplicitTwophase::SimulatorFullyImplicitTwophase(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const IncompPropsAdInterface& props, + WellsManager& wells_manager, LinearSolverInterface& linsolver, - std::vector& src) + std::vector& src, + const double* gravity) { - pimpl_.reset(new Impl(param, grid, props, linsolver, src)); + pimpl_.reset(new Impl(param, grid, props, wells_manager, linsolver, src, gravity)); } @@ -98,9 +111,10 @@ namespace Opm SimulatorReport SimulatorFullyImplicitTwophase::run(SimulatorTimer& timer, TwophaseState& state, - std::vector& src) + std::vector& src, + WellState& well_state) { - return pimpl_->run(timer, state, src); + return pimpl_->run(timer, state, src, well_state); } @@ -133,6 +147,75 @@ namespace Opm dm["velocity"] = &cell_velocity; Opm::writeVtkData(grid, dm, vtkfile); } + + + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + static void outputWellStateMatlab(const Opm::WellState& well_state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["bhp"] = &well_state.bhp(); + dm["wellrates"] = &well_state.wellRates(); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error,"Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + + /* static void outputWaterCut(const Opm::Watercut& watercut, const std::string& output_dir) @@ -145,20 +228,36 @@ namespace Opm } watercut.write(os); } - - */ - + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + wellreport.write(os); + } + */ SimulatorFullyImplicitTwophase::Impl::Impl(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const IncompPropsAdInterface& props, + WellsManager& wells_manager, LinearSolverInterface& linsolver, - std::vector& src) + std::vector& src, + const double* gravity) : grid_(grid), props_(props), - solver_(grid_, props_, linsolver) + wells_manager_(wells_manager), + wells_ (wells_manager.c_wells()), + src_ (src), + gravity_(gravity), + solver_(grid_, props_, *wells_manager.c_wells(), linsolver, gravity_) + { // For output. output_ = param.getDefault("output", true); @@ -175,6 +274,9 @@ namespace Opm } output_interval_ = param.getDefault("output_interval", 1); } + // Well control related init. + check_well_controls_ = param.getDefault("check_well_controls", false); + max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); // Misc init. const int num_cells = grid.number_of_cells; @@ -184,16 +286,16 @@ namespace Opm } } - - - SimulatorReport SimulatorFullyImplicitTwophase::Impl::run(SimulatorTimer& timer, TwophaseState& state, - std::vector& src) + std::vector& src, + WellState& well_state) { + // Initialisation. std::vector porevol; - computePorevolume(grid_, props_.porosity(), porevol); + Opm::computePorevolume(grid_, props_.porosity(), porevol); + // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); std::vector initial_porevol = porevol; @@ -203,14 +305,27 @@ namespace Opm Opm::time::StopWatch step_timer; Opm::time::StopWatch total_timer; total_timer.start(); +#if 0 // These must be changed for three-phase. + double init_surfvol[2] = { 0.0 }; + double inplace_surfvol[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol); + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); +#endif std::vector fractional_flows; + std::vector well_resflows_phase; + if (wells_) { + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + } std::fstream tstep_os; if (output_) { std::string filename = output_dir_ + "/step_timing.param"; tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); } - for (; !timer.done(); ++timer) { + while (!timer.done()) { // Report timestep and (optionally) write state to disk. step_timer.start(); timer.report(std::cout); @@ -218,34 +333,125 @@ namespace Opm if (output_vtk_) { outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); + } - + SimulatorReport sreport; - // Run solver. - solver_timer.start(); - std::vector initial_pressure = state.pressure(); - solver_.step(timer.currentStepLength(), state, src); + bool well_control_passed = !check_well_controls_; + int well_control_iteration = 0; + do { + // Run solver. + solver_timer.start(); + std::vector initial_pressure = state.pressure(); + solver_.step(timer.currentStepLength(), state, src, well_state); - // Stop timer and report. - solver_timer.stop(); - const double st = solver_timer.secsSinceStart(); - std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; - stime += st; - sreport.pressure_time = st; + // Stop timer and report. + solver_timer.stop(); + const double st = solver_timer.secsSinceStart(); + std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; + stime += st; + sreport.pressure_time = st; + + // Optionally, check if well controls are satisfied. + if (check_well_controls_) { + Opm::computePhaseFlowRatesPerWell(*wells_, + well_state.perfRates(), + fractional_flows, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { + OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); + } + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; + } + } + } while (!well_control_passed); + + // Update pore volumes if rock is compressible. + initial_porevol = porevol; + + // The reports below are geared towards two phases only. +#if 0 + // Report mass balances. + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + Opm::computeInjectedProduced(props_, state, transport_src, stepsize, + injected, produced); + Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + std::cout.precision(5); + const int width = 18; + std::cout << "\nMass balance report.\n"; + std::cout << " Injected surface volumes: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] << std::endl; + std::cout << " Produced surface volumes: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] << std::endl; + std::cout << " Total inj surface volumes: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] << std::endl; + std::cout << " Total prod surface volumes: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] << std::endl; + const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0], + init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] }; + std::cout << " Initial - inplace + inj - prod: " + << std::setw(width) << balance[0] + << std::setw(width) << balance[1] + << std::endl; + std::cout << " Relative mass error: " + << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0]) + << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1]) + << std::endl; + std::cout.precision(8); + + // Make well reports. + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, + state.pressure(), state.surfacevol(), state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } +#endif sreport.total_time = step_timer.secsSinceStart(); if (output_) { sreport.reportParam(tstep_os); - } - } - if (output_) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); +#if 0 + outputWaterCut(watercut, output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } +#endif + tstep_os.close(); } - // outputWaterCut(watercut, output_dir_); - tstep_os.close(); + + // advance to next timestep before reporting at this location + ++timer; + + // write an output file for later inspection } total_timer.stop(); @@ -258,4 +464,7 @@ namespace Opm } + + + } // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp index 0dac9641c..c70a65e77 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp @@ -32,6 +32,8 @@ namespace Opm class LinearSolverInterface; class SimulatorTimer; class TwophaseState; + class WellsManager; + class WellState; struct SimulatorReport; /// Class collecting all necessary components for a two-phase simulation. @@ -60,8 +62,10 @@ namespace Opm SimulatorFullyImplicitTwophase(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const IncompPropsAdInterface& props, + WellsManager& wells_manager, LinearSolverInterface& linsolver, - std::vector& src); + std::vector& src, + const double* gravity); /// Run the simulation. /// This will run succesive timesteps until timer.done() is true. It will @@ -72,7 +76,8 @@ namespace Opm /// \return simulation report, with timing data SimulatorReport run(SimulatorTimer& timer, TwophaseState& state, - std::vector& src); + std::vector& src, + WellState& well_state); private: class Impl;