From a358da7afaf388fea0c480906ca8a0fbb93fad7c Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Tue, 10 Dec 2013 20:36:30 +0800 Subject: [PATCH] fix the polymer source bug, warnning: water initial saturation should not bt zero when running this simulator. --- examples/.sim_poly2p_fincomp_ad.cpp.swp | Bin 20480 -> 0 bytes examples/sim_2p_fincomp_ad.cpp | 28 +++--- examples/sim_poly2p_fincomp_ad.cpp | 21 ++-- .../FullyImplicitTwophasePolymerSolver.cpp | 91 +++++++++--------- .../FullyImplicitTwophasePolymerSolver.hpp | 12 ++- 5 files changed, 79 insertions(+), 73 deletions(-) delete mode 100644 examples/.sim_poly2p_fincomp_ad.cpp.swp diff --git a/examples/.sim_poly2p_fincomp_ad.cpp.swp b/examples/.sim_poly2p_fincomp_ad.cpp.swp deleted file mode 100644 index 00e8add3fc22b334f8bb9566a5eba808f776ff56..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 20480 zcmeHPYm6jS6)wORFC`KJ{tz#{WT~0yt?r(6VdD6u0YpQ5{XnDNxsR@@?&_IE zK#8eNzNzke&pqedb02ljxo3B|ee1$ecAfQJ1D{I_q9T5&%l?CHnBR>+eU-*m!$z|&Uhhc*??_U)b@xpCOH z9nW(6{T&6=1hoZf3shTRkoETLzu4G2+i4TQuEtgD(s!P!rmahB3)B{fE<^=X><|LUGu(w)(kdP+sQZxu{mj^+*r*k_y!1i>kxnW(KN`?iQ*X^!z zH;a>bFCHuh{2&7Zk@A5&5Oc_AN-?OW&*_?sx;~jW4*7~BtCcAhlenKA@`h&U({Ci1 znA$F)K(Qnfj(l8u-3d~A-A~;nWAv~?jEX32GF9t%X_MiPD4om;M$UZKcAYev7o%6u zktdnY({A^Wm&Qpy^W*3SzUr*|@gO;95^F#5_~w)x=92DK$}v4cep4vT!#)bNTrm#p_nnGkZsZ-f>(zv>>OS3fG7?bww zzL6R+HK?XDS|g?iY9Zi#Q?dMPN2$SowFuoVK##&9}=cgHHX z%%^>2)4W=?WpgyBc_|uo*owc-Qx?U{U3H=r&isgBp@iI7_5&<6gQkK~w6y3dC6*{4 zv~6|vji_8$#+5Acbxx9LRAO54OBl$pf&*+Hn`SPLG9S;y*iX4!gr#P*peN9!6&v)h zt=Ojft4_+X8+NSrv=B5c#HdbHwO4yoxpiB-f~0bT6q$)EazdW=9hb`)(@%KM-`s|P z-i)%ECMCUTq-rqYY`I*#%22$L_?`vAWuDk^Prl(Kr5i_S#!lf+GH|m&!o4G)-7Q+6 z;WV|xE*d#2Jh5DxewwPsM0ixz2y#tkM%b<+?99=f*q!m+qh0#W^5A=Sf;$F7YN2P< z7Ii_dU0t|YEEGMaww*#xsoCgyG#k+trroUGb3{X8RYAB}Bs}Z}em~{ZyF~{%p$XnF z?iw17RhU$P?9KprqaiQO16Y#=wp3pfJuha9q{60AqH8EdcN=1sm_kmR^_9S7UK7qdDS50;@pCFs=GGD zP7A@yaU8Hx&9DPl2NQG_w&fX5x?R`~M#0x3Nj&IFtD8=;g0lwxt~Xf^hYxMEighq& zoHnv#%b;-!9c)D=eUl8ZKPA)Cg3pw(3)!{7z~gKldAKbX2Z5F!x40a|*4Us~tNs4L zsz!rM>~gJPAJSI24?DIzOapS%ag#XRk@ghU@Xk~r_PkRXE8oeft>c_MlAM?0{sYcbr@@n9p_+yep^CowYaw64q*Sxibc|8Dwg!H4RR+e7WpyWIDToF`tcZ91MjY$zhbR6 z4D$aUfo~kcZzlhLM#Uf=hJXKQU=_Fo_!E5l$AB%M3tR&H8h-t^fiu8Ez`ei#_z-Y4 z@G|`SzXIO?9t0i$Qs5-e1}+1hgOC41;4a{Uz#iaf`1qd%2Ebw9UBC;_?+?K5fbRfz z0p#OZSnJsUwkFl+QxGB)Y3QH;Ri@IdNi*Izy{w78z+8AF@ZSsf90i?~ zyQXeq=$c2mfVVoqfTyNAND_+9h2WBZpAo$Cdy@jcoQTaFV4jM2_2NWE0(=Z&dyb)J zA_68NBBX6P$7W`H9mUXb)uC|F$`C!IfFF4kh}1~mqG%C4uW6ECR#LoA&Iioe1I*Rc z?hmkZJb0X_aPj*AW3i#8UH_L_;)oWr;{B zCK*N1b!1UC*U`$XvXw_SAf2`~h2S)-_%T!{nbXe-1bT{0n~leZVcii}2^q0QUpr;~xgD2Ix-ztH5dC zW56NcYT(cC<(~z<349SqflGm}BSygi8#oCt;2hwmxWoSva3637-~;CaqCaoYLqM%o zZGqYXwFPPmydn#1{{qyt6JCU+RH@fHQ>u!;K?+tUO0iacvLOnY)zP*}P zWL9;^0ew#QI%iGog2DZs=yZNeB@>K%RCW#J3%u%7D)Jj&qhu=5vV3HEZ4;`9|Dw#x zgU`QBQWX*E%%}0hs#5J(a#g9eoWny^RZ~QZ8cTJe26-mbNTSt9WjWbur1I?a$VF9Z zLenKrx%+P)?UTBCHTHAt|fSx|CxzzOCog(a?XZi9OdsX{B5&>l&oY_cBl4J2mcbFN~bU@ea90f zgy2-XM0HME!ULQ5+Py)9SJ*JQmR6q#K*$W*6T`45?f0ojW9W8|^I#EeZ$eXK92FR0 s(5hR=7=nxns8v?%71|+I(w+_x!}O;;CvjE8MANFM2 density(num_phases, 1000.0); std::vector viscosity(num_phases, 1.0*centi*Poise); - double porosity = 0.5; + viscosity[0] = 0.5 * centi * Poise; + viscosity[1] = 5 * centi * Poise; + double porosity = 0.35; double permeability = 10.0*milli*darcy; SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear; IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity, porosity, permeability, grid.dimensions, num_cells); std::vector omega; std::vector src(num_cells, 0.0); - src[0] = 1.; - src[num_cells-1] = -1.; + src[0] = 10. / day; + src[num_cells-1] = -10. / day; FlowBCManager bcs; LinearSolverUmfpack linsolver; FullyImplicitTwoPhaseSolver solver(grid, props, linsolver); std::vector porevol; Opm::computePorevolume(grid, props.porosity(), porevol); - const double dt = param.getDefault("dt", 0.1) * day; - const int num_time_steps = param.getDefault("nsteps", 20); + const double dt = param.getDefault("dt", 10) * day; + const int num_time_steps = param.getDefault("nsteps", 10); std::vector allcells(num_cells); for (int cell = 0; cell < num_cells; ++cell) { allcells[cell] = cell; @@ -66,10 +68,10 @@ try //initial sat for (int c = 0; c < num_cells; ++c) { - state.saturation()[2*c] = 0.2; - state.saturation()[2*c+1] = 0.8; + state.saturation()[2*c] = 0; + state.saturation()[2*c+1] = 1; } - std::vector p(num_cells, 200*Opm::unit::barsa); + std::vector p(num_cells, 100*Opm::unit::barsa); state.pressure() = p; // state.setFirstSat(allcells, props, TwophaseState::MinSat); std::ostringstream vtkfilename; diff --git a/examples/sim_poly2p_fincomp_ad.cpp b/examples/sim_poly2p_fincomp_ad.cpp index 80e5f3e6a..f1262e1f9 100644 --- a/examples/sim_poly2p_fincomp_ad.cpp +++ b/examples/sim_poly2p_fincomp_ad.cpp @@ -36,12 +36,12 @@ try } std::string deck_filename = param.get("deck_filename"); EclipseGridParser deck = EclipseGridParser(deck_filename); - int nx = param.getDefault("nx", 20); - int ny = param.getDefault("ny", 20); + int nx = param.getDefault("nx", 30); + int ny = param.getDefault("ny", 30); int nz = 1; - double dx = 2.0; - double dy = 2.0; - double dz = 0.5; + double dx = 10.0; + double dy = 1.0; + double dz = 1.0; GridManager grid_manager(nx, ny, nz, dx, dy, dz); const UnstructuredGrid& grid = *grid_manager.c_grid(); int num_cells = grid.number_of_cells; @@ -51,7 +51,7 @@ try std::vector density(num_phases, 1000.0); std::vector viscosity(num_phases, 1.0*centi*Poise); viscosity[0] = 0.5 * centi * Poise; - viscosity[0] = 5 * centi * Poise; + viscosity[1] = 5 * centi * Poise; double porosity = 0.35; double permeability = 10.0*milli*darcy; SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear; @@ -97,10 +97,10 @@ try std::vector src(num_cells, 0.0); std::vector src_polymer(num_cells); src[0] = 10. / day; - src[num_cells-1] = -1. / day; + src[num_cells-1] = -10. / day; - PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 30.0)*Opm::unit::day, - param.getDefault("poly_end_days", 80.0)*Opm::unit::day, + PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day, + param.getDefault("poly_end_days", 800.0)*Opm::unit::day, param.getDefault("poly_amount", polymer_props.cMax())); FlowBCManager bcs; LinearSolverUmfpack linsolver; @@ -120,7 +120,7 @@ try state.saturation()[2*c] = 0.2; state.saturation()[2*c+1] = 0.8; } - std::vector p(num_cells, 200*Opm::unit::barsa); + std::vector p(num_cells, 100*Opm::unit::barsa); state.pressure() = p; std::vector c(num_cells, 0.0); @@ -137,6 +137,7 @@ try Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); + dm["concentration"] = &state.concentration(); Opm::writeVtkData(grid, dm, vtkfile); } } diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index 8e0b90f48..f6ba27bcb 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -97,13 +97,8 @@ typedef Eigen::Array& src, const std::vector& polymer_inflow) - // const bool if_polymer_actived) { // Create the primary variables. const SolutionState state = variableState(x); @@ -256,44 +250,33 @@ typedef Eigen::Array& polymer_inflow) const - { - double amount = 0; - for (int i = 0; i < int(polymer_inflow.size()); ++i) { - amount += polymer_inflow[i]; - } - return amount; - } - - ADB FullyImplicitTwophasePolymerSolver::accumSource(const int phase, const std::vector& kr, - const std::vector& src - ) const + const std::vector& src) const { //extract the source to out and in source. std::vector outsrc; @@ -329,8 +312,7 @@ typedef Eigen::Array& kr, + polymerSource(const std::vector& kr, const std::vector& src, const std::vector& polymer_inflow_c, const SolutionState& state) const @@ -359,13 +341,15 @@ typedef Eigen::Array& kr) const + FullyImplicitTwophasePolymerSolver::computeFracFlow(int phase, + const std::vector& kr) const { const double* mus = fluid_.viscosity(); ADB mob_phase = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]); @@ -498,8 +482,25 @@ typedef Eigen::Array& kr, + const SolutionState& state) const + { + const double* mus = fluid_.viscosity(); + ADB water_mob = kr[0] / V::Constant(kr[0].size(), 1, mus[0]); + ADB poly_mob = state.concentration * mc * water_mob; + + 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 diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp index 65ad86e3b..6205f9bde 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -30,7 +30,6 @@ namespace Opm { const std::vector& src, const std::vector& polymer_inflow ); -// const bool if_polymer_actived); private: typedef AutoDiffBlock ADB; typedef ADB::V V; @@ -44,7 +43,6 @@ namespace Opm { ADB pressure; std::vector saturation; ADB concentration; -// ADB cmax; }; const UnstructuredGrid& grid_; const IncompPropsAdInterface& fluid_; @@ -65,7 +63,6 @@ namespace Opm { const PolymerState& x, const std::vector& src, const std::vector& polymer_inflow); -// const bool if_polymer_actived); V solveJacobianSystem() const; void updateState(const V& dx, PolymerState& x) const; @@ -85,6 +82,11 @@ namespace Opm { const V& trans, const std::vector& kr, const SolutionState& state) const; + ADB + computePolymerMassFlux(const V& trans, + const ADB& mc, + const std::vector& kr, + const SolutionState& state) const; double residualNorm() const; ADB @@ -103,8 +105,8 @@ namespace Opm { fluidDensity(const int phase) const; ADB transMult(const ADB& p) const; - double - polymerInjectedAmount(const std::vector& polymer_inflow) const; }; + + } // namespace Opm #endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED