From 0055ae1def5bcaabf0900f976c5adbc0214d2e2c Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Thu, 23 Jan 2014 08:48:42 +0800 Subject: [PATCH] output watercut for incomp polymer solver. --- opm/polymer/fullyimplicit/.utilities.cpp.swp | Bin 16384 -> 0 bytes .../FullyImplicitTwophasePolymerSolver.cpp | 35 ++++-------- .../FullyImplicitTwophasePolymerSolver.hpp | 7 ++- .../SimulatorFullyImplicitTwophasePolymer.cpp | 52 ++++++++++++++++-- opm/polymer/fullyimplicit/utilities.cpp | 2 +- 5 files changed, 62 insertions(+), 34 deletions(-) delete mode 100644 opm/polymer/fullyimplicit/.utilities.cpp.swp diff --git a/opm/polymer/fullyimplicit/.utilities.cpp.swp b/opm/polymer/fullyimplicit/.utilities.cpp.swp deleted file mode 100644 index b7460437f9637d580d00a4119b128a6abf20961e..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeHN%a0sK8E^6cCW#3o1VM5suM;z~kD1-Il)&smti6UfVzbu9UL!4=p{J*Mx9omV zRlVblv2w}*Bu+sbA`S>4F5CzQE?kfzB_dq-2T;W0fKbE%35g5*zUt@9&aC4T2vy5J z&vre&@2lVU)mPnB-Cx|jvc+C#FFLqB?Kpqy{bA?n2Y!9#p|3ekC<2jE`8}>a7dIDz z+s`juSpM=0r?l-PGQP<@-Qlt9q>9T9iiRh?vHJYh{#XVq1C{~H zfMvikU>UFsSOzQumI2FvWxz7OraGr(USbDZ~qUjW|+z6D$W z{`IKi{0cY%mVv)N;y7;ucYt-^_YKEMfX9J%Kj}C>0B!(J0%w4~e!_8n3;Yy#14w{1 zpaFdNamV={@Eu?W_#*J9haKm2U={f1L$C=r02Y7`u(9!7;48of*d+NS@Fvg&{s_+A z1KtL113qAQ1eL6I-5t53ijyT&s-+vLB$m3Tq_>xF>GwPy1p91mjhQ>?YQEx7)A4Jv zKx#F;hH6W(Y<7)#ahUXc>3VyM`zx%|VK<;pGRg6&#!)OJJBouebvjw))UTm>dw!{~nZ%rN_DZ0yepg6nZOw&r`)qvL6VtaAA zv{50et5wg=b`o~GX(V*FJ8+K~>SkE@>V*z&+B0-(#>JySa3A7lm^g@H$p>w=QxZx{ zvZ1Tkj5BJ_v^yvtW(LtfVOMTANE{+m)WSOb^t?EVl348x_fM5;#^J@OHX>(n!R0c- z7dc(#$yg~rh=)u>M)oN-3ZFG|_4+vBP|MD)u`4gXbotV??VZ*sE<=!29utbWOvl;{ zpwpAN1?nh>2Dkz`t>VVz=PnmvO!?FYPO&zc!`Pic})LU!Gp8xmc+yXDdWS$?R}+sHKfV^jKV>fMn5-6b{F=nE*S&CttB#&GfCajV44! z$=F0cnWna{ZuPFcxU+pDOCvIkvF~~U41(h2pj~- zGD{9!qzn)LwP(gOQe;F41?uv+*_9K~p!8mphJ7x3vEMUN)hH{qKEzAznW~Q|sJfI> zrmmY(E(Hu^X({^n)(z2AW!jPTp^YJm*LXAX=1ovnqDPx~0e5AyHg<}x$b*Ua8ZS|8 zSa^vdvS@PW)QmxSX@XPFkFH&U~)g%*qWfvd$4$C7Qb6*8zvSFz9fu|CYA zPO)Z{zOt3GmHF}>3g%i2xf>rjl7c)DLb5L;GE37WJhU)@eNsa;gl{IUbi+N^wm<3} zcI}W?B}|!#!@xwS+EjPUT!{uT>vKwETqcrhK41f$V9mqgz~YE0tY%%~DH=_#OGNvu z=1C&a#4`jaBy5_K(SZ3udK=$wPVKCypSh^BVgpwr-)otA8ZsvLGqGHBV_ywYL;~rO zv3JY(7LOzJ(WpahD*n-J0I`^6=Aj1mVdM&*mqg4mdNT^0vC>UMYH;g)+dmg{U2PPiC=%_`5F)^YPVazo)|c1~_r~nxkBB4N5gg2t;QZ1>fy0C z|FHYGa~;p!6Qz%|x=c?2TP{7MPO`+bRU%9S@>Zu@XC}2RzBA8>DV7McE}hnX^t<$; z+xW^dNrJJ?WUNlTTb82AHEC+kKtu?WPCsx}FK`j3x7s5yR8l&m^f{&OK#BoX(fj|; z@Gky2yu;J`e{p{Q4&L{F1NnWzTUPlalRbOXj|f)u9EDrlT0dy3^I@oI2#T8%eJ?7;Rf2+ f&$f7aFa;a@`JFk--R0QA6Z#kjSt4$FI2HI0YMeuL diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index c3e99d1b2..f24ac775a 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -150,7 +150,8 @@ namespace { step(const double dt, PolymerState& x, WellState& xw, - const std::vector& polymer_inflow) + const std::vector& polymer_inflow, + std::vector& src) { V pvol(grid_.number_of_cells); @@ -169,7 +170,7 @@ namespace { const double atol = 1.0e-12; const double rtol = 5.0e-8; const int maxit = 40; - assemble(pvdt, old_state, x, xw, polymer_inflow); + assemble(pvdt, x, xw, polymer_inflow, src); const double r0 = residualNorm(); int it = 0; @@ -181,7 +182,7 @@ namespace { const V dx = solveJacobianSystem(); updateState(dx, x, xw); - assemble(pvdt, old_state, x, xw, polymer_inflow); + assemble(pvdt, x, xw, polymer_inflow, src); const double r = residualNorm(); @@ -393,10 +394,10 @@ namespace { void FullyImplicitTwophasePolymerSolver:: assemble(const V& pvdt, - const SolutionState& old_state, const PolymerState& x, const WellState& xw, - const std::vector& polymer_inflow) + const std::vector& polymer_inflow, + std::vector& src) { // Create the primary variables. const SolutionState state = variableState(x, xw); @@ -404,27 +405,10 @@ namespace { // -------- Mass balance equations for water and oil -------- const V trans = subset(transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); - const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); - // const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); const ADB mc = computeMc(state); computeMassFlux(trans, mc, kr[1], krw_eff, state); - //const std::vector source = accumSource(kr[1], krw_eff, state.concentration, src, polymer_inflow); - // const std::vector source = polymerSource(); - // const double rho_r = polymer_props_ad_.rockDensity(); -// const V phi = V::Constant(pvdt.size(), 1, *fluid_.porosity()); - -// const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); -// residual_.mass_balance[0] = pvdt * (state.saturation[0] - old_state.saturation[0]) -// + ops_.div * mflux[0]; -// residual_.mass_balance[1] = pvdt * (state.saturation[1] - old_state.saturation[1]) -// + ops_.div * mflux[1]; - // Mass balance equation for polymer -// residual_.mass_balance[2] = pvdt * (state.saturation[0] * state.concentration -// - old_state.saturation[0] * old_state.concentration) * (1. - dead_pore_vol) -// + pvdt * rho_r * (1. - phi) / phi * ads - // + ops_.div * mflux[2]; residual_.mass_balance[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0]) + ops_.div*rq_[0].mflux; residual_.mass_balance[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0]) @@ -516,14 +500,17 @@ namespace { well_contribs[phase] = superset(perf_flux, well_cells, nc); // DUMP(well_contribs[phase]); residual_.mass_balance[phase] += well_contribs[phase]; + for (int i = 0; i < nc; ++i) { + src[i] += well_contribs[phase].value()[i]; + } } // well rates contribs to polymer mass balance eqn. // for injection wells. const V polyin = Eigen::Map(& polymer_inflow[0], nc); const V poly_in_perf = subset(polyin, well_cells); - const V poly_c_cell = subset(state.concentration, well_cells).value(); - const V poly_c = producer.select(poly_c_cell, poly_in_perf); + const V poly_mc_cell = subset(mc, well_cells).value(); + const V poly_c = producer.select(poly_mc_cell, poly_in_perf); residual_.mass_balance[2] += superset(well_perf_rates[0] * poly_c, well_cells, nc); // Set the well flux equation diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp index d96b19750..525b80a7d 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -29,7 +29,8 @@ namespace Opm { void step(const double dt, PolymerState& state, WellState& well_state, - const std::vector& polymer_inflow); + const std::vector& polymer_inflow, + std::vector& src); private: typedef AutoDiffBlock ADB; typedef ADB::V V; @@ -88,10 +89,10 @@ namespace Opm { const WellState& xw); void assemble(const V& pvdt, - const SolutionState& old_state, const PolymerState& x, const WellState& xw, - const std::vector& polymer_inflow); + const std::vector& polymer_inflow, + std::vector& src); V solveJacobianSystem() const; void updateState(const V& dx, PolymerState& x, diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp index 34d19a8ea..291207b08 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -146,6 +147,7 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); dm["concentration"] = &state.concentration(); std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); @@ -162,6 +164,7 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); dm["concentration"] = &state.concentration(); std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); @@ -191,7 +194,6 @@ namespace Opm -/* static void outputWaterCut(const Opm::Watercut& watercut, const std::string& output_dir) { @@ -203,6 +205,7 @@ namespace Opm } watercut.write(os); } +/* static void outputWellReport(const Opm::WellReport& wellreport, const std::string& output_dir) { @@ -298,8 +301,9 @@ namespace Opm std::vector porevol; Opm::computePorevolume(grid_, props_.porosity(), porevol); - // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); std::vector polymer_inflow_c(grid_.number_of_cells); + std::vector transport_src(grid_.number_of_cells); // Main simulation loop. Opm::time::StopWatch solver_timer; @@ -307,6 +311,10 @@ namespace Opm Opm::time::StopWatch step_timer; Opm::time::StopWatch total_timer; total_timer.start(); + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); #if 0 // These must be changed for three-phase. double init_surfvol[2] = { 0.0 }; @@ -333,7 +341,7 @@ namespace Opm outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); } outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); + // outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); } @@ -348,7 +356,7 @@ namespace Opm polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c); solver_timer.start(); std::vector initial_pressure = state.pressure(); - solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c); + solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src); // Stop timer and report. solver_timer.stop(); @@ -378,6 +386,38 @@ namespace Opm } } } while (!well_control_passed); + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double polyinj = 0; + double polyprod = 0; + + Opm::computeInjectedProduced(props_, polymer_props_, + state, + transport_src, polymer_inflow_c, timer.currentStepLength(), + injected, produced, + polyinj, polyprod); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + std::cout.precision(5); + const int width = 18; + std::cout << "\nMass balance report.\n"; + std::cout << " Injected reservoir volumes: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] << std::endl; + std::cout << " Produced reservoir volumes: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] << std::endl; + std::cout << " Total inj reservoir volumes: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] << std::endl; + std::cout << " Total prod reservoir volumes: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] << std::endl; // Update pore volumes if rock is compressible. @@ -440,9 +480,9 @@ namespace Opm 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 0 + outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); if (wells_) { outputWellReport(wellreport, output_dir_); } diff --git a/opm/polymer/fullyimplicit/utilities.cpp b/opm/polymer/fullyimplicit/utilities.cpp index bb03ff1cd..2008795a0 100644 --- a/opm/polymer/fullyimplicit/utilities.cpp +++ b/opm/polymer/fullyimplicit/utilities.cpp @@ -144,7 +144,7 @@ namespace Opm const V inv_muw_eff = polymer_props.effectiveInvWaterVisc(c, mus); std::vector mob(np); mob[0] = krw_eff * inv_muw_eff; - mob[1] = kr[1] / mu[1]; + mob[1] = kr[1] / mus[1]; const V watmob_c = src_selector.select(mob[0], one); const V oilmob_c = src_selector.select(mob[1], zero);