From 8592ca825a138ff8a23d6bc073795c5135bc321c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 21 May 2013 23:54:30 +0200 Subject: [PATCH 01/12] Added wellRates() member. --- opm/core/simulator/WellState.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index f1edc5d7a..97c71e401 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -54,6 +54,7 @@ namespace Opm } perfrates_.resize(wells->well_connpos[nw], 0.0); perfpress_.resize(wells->well_connpos[nw], -1e100); + wellrates_.resize(wells->well_connpos[nw] * wells->number_of_phases, 0.0); } } @@ -61,6 +62,10 @@ namespace Opm std::vector& bhp() { return bhp_; } const std::vector& bhp() const { return bhp_; } + /// One rate per well and phase. + std::vector& wellRates() { return wellrates_; } + const std::vector& wellRates() const { return wellrates_; } + /// One rate per well connection. std::vector& perfRates() { return perfrates_; } const std::vector& perfRates() const { return perfrates_; } @@ -71,6 +76,7 @@ namespace Opm private: std::vector bhp_; + std::vector wellrates_; std::vector perfrates_; std::vector perfpress_; }; From 1cf0bb465b07d25563c355422618cc06ba59c945 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 May 2013 09:20:46 +0200 Subject: [PATCH 02/12] Fix bug in upwinding code. Accidental usage of std::vector's operator< discovered. --- opm/core/pressure/CompressibleTpfa.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 3dcb20440..436399029 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -429,7 +429,7 @@ namespace Opm for (int phase = 0; phase < np; ++phase) { int upwindc = -1; if (c[0] >=0 && c[1] >= 0) { - upwindc = (pot[0] < pot[1]) ? c[1] : c[0]; + upwindc = (pot[0][phase] < pot[1][phase]) ? c[1] : c[0]; } else { upwindc = (c[0] >= 0) ? c[0] : c[1]; } From 125c57346116e904e457a8cc68bdbe4f2299ab93 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 May 2013 09:21:41 +0200 Subject: [PATCH 03/12] Added assert to guard against wrong usage. --- opm/core/simulator/BlackoilState.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp index 153c2eed5..caf570e5c 100644 --- a/opm/core/simulator/BlackoilState.hpp +++ b/opm/core/simulator/BlackoilState.hpp @@ -22,6 +22,7 @@ #include #include +#include #include namespace Opm @@ -62,6 +63,7 @@ namespace Opm const int n = cells.size(); std::vector smin(num_phases_*n); std::vector smax(num_phases_*n); + ASSERT(n > 0); props.satRange(n, &cells[0], &smin[0], &smax[0]); const double* svals = (es == MinSat) ? &smin[0] : &smax[0]; for (int ci = 0; ci < n; ++ci) { From 6e9a46b34de9ff179abfe3ea359173fb645267b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 May 2013 12:53:06 +0200 Subject: [PATCH 04/12] Added gasoilratio() to BlackoilState. --- opm/core/simulator/BlackoilState.hpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp index caf570e5c..974deffd8 100644 --- a/opm/core/simulator/BlackoilState.hpp +++ b/opm/core/simulator/BlackoilState.hpp @@ -48,6 +48,7 @@ namespace Opm // but use available phase information instead. sat_[num_phases*cell + 1] = 1.0; } + gor_.resize(g.number_of_cells, 0.0); } enum ExtremalSat { MinSat, MaxSat }; @@ -60,10 +61,13 @@ namespace Opm const Opm::BlackoilPropertiesInterface& props, ExtremalSat es) { + if (cells.empty()) { + return; + } const int n = cells.size(); + ASSERT(n > 0); std::vector smin(num_phases_*n); std::vector smax(num_phases_*n); - ASSERT(n > 0); props.satRange(n, &cells[0], &smin[0], &smax[0]); const double* svals = (es == MinSat) ? &smin[0] : &smax[0]; for (int ci = 0; ci < n; ++ci) { @@ -83,12 +87,14 @@ namespace Opm std::vector& faceflux () { return flux_ ; } std::vector& surfacevol () { return surfvol_; } std::vector& saturation () { return sat_ ; } + std::vector& gasoilratio () { return gor_ ; } const std::vector& pressure () const { return press_ ; } const std::vector& facepressure() const { return fpress_; } const std::vector& faceflux () const { return flux_ ; } const std::vector& surfacevol () const { return surfvol_; } const std::vector& saturation () const { return sat_ ; } + const std::vector& gasoilratio () const { return gor_ ; } private: int num_phases_; @@ -97,6 +103,7 @@ namespace Opm std::vector flux_ ; std::vector surfvol_; std::vector sat_ ; + std::vector gor_ ; }; } // namespace Opm From a0a9c482eac144c20a14ac9ae944e283ef220a45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 May 2013 15:44:07 +0200 Subject: [PATCH 05/12] Created new initialization routine, using RS from deck. --- opm/core/simulator/initState_impl.hpp | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index b6492bb98..ed9d4fd97 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -599,6 +599,29 @@ namespace Opm } + /// Initialize a blackoil state from input deck. + template + void initBlackoilStateFromDeck(const UnstructuredGrid& grid, + const Props& props, + const EclipseGridParser& deck, + const double gravity, + State& state) + { + initStateFromDeck(grid, props, deck, gravity, state); + initBlackoilSurfvol(grid, props, state); + if (deck.hasField("RS")) { + const std::vector& rs_deck = deck.getFloatingPointValue("RS"); + const int num_cells = grid.number_of_cells; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + state.gasoilratio()[c] = rs_deck[c_deck]; + } + } else { + THROW("Temporarily, we require the RS field."); + } + } + + } // namespace Opm From c280f2b0dd9620d203a5aad7d2e181a35ee22184 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 30 May 2013 11:03:08 +0200 Subject: [PATCH 06/12] Initialize rate-controlled well bhp with safety factor. Safety factor is 1.01 (INJECTOR) or 0.99 (PRODUCER), similar to mrst's ad-fi/utils/initWellSolLocal.m > initialize(). --- opm/core/simulator/WellState.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 97c71e401..3915ed0e0 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -46,7 +46,8 @@ namespace Opm if ((ctrl->current < 0) || // SHUT (ctrl->type[ctrl->current] != BHP)) { const int cell = wells->well_cells[wells->well_connpos[w]]; - bhp_[w] = state.pressure()[cell]; + const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99; + bhp_[w] = safety_factor*state.pressure()[cell]; } else { bhp_[w] = ctrl->target[ctrl->current]; From ad4c9a47e04eeeb72060052bb89e3ec6780c80c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 21:58:30 +0200 Subject: [PATCH 07/12] Refined well state initialization. For SURFACE_RATE controlled wells, initialize wellRates() to match. --- opm/core/simulator/WellState.hpp | 35 ++++++++++++++++++++++++-------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 3915ed0e0..b8b7a31d0 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -31,31 +31,48 @@ namespace Opm { public: /// Allocate and initialize if wells is non-null. + /// Also tries to give useful initial values to the bhp() and + /// wellRates() fields, depending on controls. The + /// perfRates() field is filled with zero, and perfPress() + /// with -1e100. template void init(const Wells* wells, const State& state) { if (wells) { const int nw = wells->number_of_wells; + const int np = wells->number_of_phases; bhp_.resize(nw); - // Initialize bhp to be target pressure - // if bhp-controlled well, otherwise set - // to pressure in first perforation cell. + wellrates_.resize(nw * np, 0.0); for (int w = 0; w < nw; ++w) { const WellControls* ctrl = wells->ctrls[w]; - + // Initialize bhp to be target pressure if + // bhp-controlled well, otherwise set to a little + // above or below (depending on if the well is an + // injector or producer) pressure in first perforation + // cell. if ((ctrl->current < 0) || // SHUT (ctrl->type[ctrl->current] != BHP)) { - const int cell = wells->well_cells[wells->well_connpos[w]]; + const int first_cell = wells->well_cells[wells->well_connpos[w]]; const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99; - bhp_[w] = safety_factor*state.pressure()[cell]; - } - else { + bhp_[w] = safety_factor*state.pressure()[first_cell]; + } else { bhp_[w] = ctrl->target[ctrl->current]; } + // Initialize well rates to match controls if type is SURFACE_RATE + if ((ctrl->current >= 0) || // open well + (ctrl->type[ctrl->current] != SURFACE_RATE)) { + const double rate_target = ctrl->target[ctrl->current]; + for (int p = 0; p < np; ++p) { + const double phase_distr = ctrl->distr[np * ctrl->current + p]; + wellrates_[np*w + p] = rate_target * phase_distr; + } + } } + // The perforation rates and perforation pressures are + // not expected to be consistent with bhp_ and wellrates_ + // after init(). perfrates_.resize(wells->well_connpos[nw], 0.0); perfpress_.resize(wells->well_connpos[nw], -1e100); - wellrates_.resize(wells->well_connpos[nw] * wells->number_of_phases, 0.0); } } From 175cecacdaa3442d663c720341badf6b73495d03 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 22:19:43 +0200 Subject: [PATCH 08/12] Bugfix in well rate init. Do not always try to initialize, also initialize proper phase rates. --- opm/core/simulator/WellState.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index b8b7a31d0..eb4e583df 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -59,11 +59,12 @@ namespace Opm bhp_[w] = ctrl->target[ctrl->current]; } // Initialize well rates to match controls if type is SURFACE_RATE - if ((ctrl->current >= 0) || // open well + if ((ctrl->current >= 0) && // open well (ctrl->type[ctrl->current] != SURFACE_RATE)) { const double rate_target = ctrl->target[ctrl->current]; for (int p = 0; p < np; ++p) { - const double phase_distr = ctrl->distr[np * ctrl->current + p]; + const double phase_distr = ctrl->distr[np * ctrl->current + p] + * wells->comp_frac[np * w + p]; wellrates_[np*w + p] = rate_target * phase_distr; } } From 5d457ff70816792ad384262c6d8358bd869ece8f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 23:30:43 +0200 Subject: [PATCH 09/12] Ensures well rate initialization actually happens. Do not use the well's comp_frac member, only rely on the control's distr member for initialization. This forced a change to WellsManager's initialization of the distr member. --- opm/core/simulator/WellState.hpp | 5 ++--- opm/core/wells/WellsManager.cpp | 24 ++++++++++++++++++++++-- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index eb4e583df..9849bd680 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -60,11 +60,10 @@ namespace Opm } // Initialize well rates to match controls if type is SURFACE_RATE if ((ctrl->current >= 0) && // open well - (ctrl->type[ctrl->current] != SURFACE_RATE)) { + (ctrl->type[ctrl->current] == SURFACE_RATE)) { const double rate_target = ctrl->target[ctrl->current]; for (int p = 0; p < np; ++p) { - const double phase_distr = ctrl->distr[np * ctrl->current + p] - * wells->comp_frac[np * w + p]; + const double phase_distr = ctrl->distr[np * ctrl->current + p]; wellrates_[np*w + p] = rate_target * phase_distr; } } diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 328abd495..a32545595 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -505,13 +505,33 @@ namespace Opm int control_pos[5] = { -1, -1, -1, -1, -1 }; if (ok && wci_line.surface_flow_max_rate_ >= 0.0) { control_pos[InjectionControl::RATE] = w_->ctrls[wix]->num; - const double distr[3] = { 1.0, 1.0, 1.0 }; + double distr[3] = { 0.0, 0.0, 0.0 }; + if (wci_line.injector_type_ == "WATER") { + distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (wci_line.injector_type_ == "OIL") { + distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (wci_line.injector_type_ == "GAS") { + distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } else { + THROW("Injector type " << wci_line.injector_type_ << " not supported." + "WellsManager only supports WATER, OIL and GAS injector types."); + } ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_, distr, wix, w_); } if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) { control_pos[InjectionControl::RESV] = w_->ctrls[wix]->num; - const double distr[3] = { 1.0, 1.0, 1.0 }; + double distr[3] = { 0.0, 0.0, 0.0 }; + if (wci_line.injector_type_ == "WATER") { + distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (wci_line.injector_type_ == "OIL") { + distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (wci_line.injector_type_ == "GAS") { + distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } else { + THROW("Injector type " << wci_line.injector_type_ << " not supported." + "WellsManager only supports WATER, OIL and GAS injector types."); + } ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_, distr, wix, w_); } From 55e58deea626f5468a75912eb1c3aad1a1ed9e00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Jun 2013 14:33:18 +0200 Subject: [PATCH 10/12] Add transMult(), poroMultDeriv() and transMultDeriv(). --- opm/core/props/rock/RockCompressibility.cpp | 36 +++++++++++++++++++-- opm/core/props/rock/RockCompressibility.hpp | 10 ++++++ 2 files changed, 44 insertions(+), 2 deletions(-) diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index b290f9be7..fd072798c 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -42,12 +42,16 @@ namespace Opm { if (deck.hasField("ROCKTAB")) { const table_t& rt = deck.getROCKTAB().rocktab_; - int n = rt[0][0].size(); + if (rt.size() != 1) { + THROW("Can only handle a single region in ROCKTAB."); + } + const int n = rt[0][0].size(); p_.resize(n); poromult_.resize(n); for (int i = 0; i < n; ++i) { p_[i] = rt[0][0][i]; poromult_[i] = rt[0][1][i]; + transmult_[i] = rt[0][2][i]; } } else if (deck.hasField("ROCK")) { const ROCK& r = deck.getROCK(); @@ -70,11 +74,39 @@ namespace Opm const double cpnorm = rock_comp_*(pressure - pref_); return (1.0 + cpnorm + 0.5*cpnorm*cpnorm); } else { - // return Opm::linearInterpolation(p_, poromult_, pressure); return Opm::linearInterpolation(p_, poromult_, pressure); } } + double RockCompressibility::poroMultDeriv(double pressure) const + { + if (p_.empty()) { + // Approximating poro multiplier with a quadratic curve, + // we must use its derivative. + return rock_comp_ + 2 * rock_comp_ * rock_comp_ * (pressure - pref_); + } else { + return Opm::linearInterpolationDerivative(p_, poromult_, pressure); + } + } + + double RockCompressibility::transMult(double pressure) const + { + if (p_.empty()) { + return 1.0; + } else { + return Opm::linearInterpolation(p_, transmult_, pressure); + } + } + + double RockCompressibility::transMultDeriv(double pressure) const + { + if (p_.empty()) { + return 0.0; + } else { + return Opm::linearInterpolationDerivative(p_, transmult_, pressure); + } + } + double RockCompressibility::rockComp(double pressure) const { if (p_.empty()) { diff --git a/opm/core/props/rock/RockCompressibility.hpp b/opm/core/props/rock/RockCompressibility.hpp index a5c9a983b..03567caed 100644 --- a/opm/core/props/rock/RockCompressibility.hpp +++ b/opm/core/props/rock/RockCompressibility.hpp @@ -47,12 +47,22 @@ namespace Opm /// Porosity multiplier. double poroMult(double pressure) const; + /// Derivative of porosity multiplier with respect to pressure. + double poroMultDeriv(double pressure) const; + + /// Transmissibility multiplier. + double transMult(double pressure) const; + + /// Derivative of transmissibility multiplier with respect to pressure. + double transMultDeriv(double pressure) const; + /// Rock compressibility = (d poro / d p)*(1 / poro). double rockComp(double pressure) const; private: std::vector p_; std::vector poromult_; + std::vector transmult_; double pref_; double rock_comp_; }; From d7d80e6225283f1affc1afefe046f85cc7a75b73 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 5 Jun 2013 09:40:38 +0200 Subject: [PATCH 11/12] A test that compares the blackoil fluid based on the [p,r] interface with the blackoil fluid based on the [p,z] interface --- tests/testFluid.DATA | 46 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 tests/testFluid.DATA diff --git a/tests/testFluid.DATA b/tests/testFluid.DATA new file mode 100644 index 000000000..6bc3dabc0 --- /dev/null +++ b/tests/testFluid.DATA @@ -0,0 +1,46 @@ +OIL +WATER +GAS + +PVTO +-- Rs Pbub Bo Vo + .0 14.7 1.0000 1.20 / + .165 400. 1.0120 1.17 / + .335 800. 1.0255 1.14 / + .500 1200. 1.0380 1.11 / + .665 1600. 1.0510 1.08 / + .828 2000. 1.0630 1.06 / + .985 2400. 1.0750 1.03 / + 1.130 2800. 1.0870 1.00 / + 1.270 3200. 1.0985 .98 / + 1.390 3600. 1.1100 .95 / + 1.500 4000. 1.1200 .94 + 5000. 1.1189 .94 / + / + +PVDG +-- Pg Bg Vg + 14.7 178.08 .0125 + 400. 5.4777 .0130 + 800. 2.7392 .0135 + 1200. 1.8198 .0140 + 1600. 1.3648 .0145 + 2000. 1.0957 .0150 + 2400. 0.9099 .0155 + 2800. 0.7799 .0160 + 3200. 0.6871 .0165 + 3600. 0.6035 .0170 + 4000. 0.5432 .0175 / + +PVTW +--Depth Bw Comp Vw Cv + 3600. 1.0034 1.0E-6 0.96 0.0 / + + +DENSITY +-- Oil Water Gas + 44.98 63.01 0.0702 / + + + + From a1b0a0af8351e5416b65d1a1207fc16e2109272c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 5 Jun 2013 12:24:23 +0200 Subject: [PATCH 12/12] Mark deck as FIELD units. The numbers in the deck are more indicative of FIELD unit conventions than METRIC unit conventions, so allow the input parser to interpret the data in that manner. --- tests/testFluid.DATA | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/testFluid.DATA b/tests/testFluid.DATA index 6bc3dabc0..78c598cca 100644 --- a/tests/testFluid.DATA +++ b/tests/testFluid.DATA @@ -2,6 +2,8 @@ OIL WATER GAS +FIELD + PVTO -- Rs Pbub Bo Vo .0 14.7 1.0000 1.20 / @@ -40,7 +42,3 @@ PVTW DENSITY -- Oil Water Gas 44.98 63.01 0.0702 / - - - -