From 9ef601496459f753b575da822db1c07dce14aed5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Mon, 27 Aug 2012 12:20:03 +0200 Subject: [PATCH 01/10] Addes support for initialisation of three phases using SWOF and SGOF --- opm/core/utility/initState_impl.hpp | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index c50db8a91..f66704348 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -512,11 +512,12 @@ namespace Opm State& state) { const int num_phases = props.numPhases(); - if (num_phases != 2) { - THROW("initStateFromDeck(): currently handling only two-phase scenarios."); - } + state.init(grid, num_phases); if (deck.hasField("EQUIL")) { + if (num_phases != 2) { + THROW("initStateFromDeck(): currently handling only two-phase scenarios."); + } // Set saturations depending on oil-water contact. const EQUIL& equil= deck.getEQUIL(); if (equil.equil.size() != 1) { @@ -535,11 +536,22 @@ namespace Opm const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& p_deck = deck.getFloatingPointValue("PRESSURE"); 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]; - s[2*c] = sw_deck[c_deck]; - s[2*c + 1] = 1.0 - s[2*c]; - p[c] = p_deck[c_deck]; + if(num_phases == 2){ + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c] = sw_deck[c_deck]; + s[2*c + 1] = 1.0 - s[2*c]; + p[c] = p_deck[c_deck]; + } + }else{ + const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c] = sw_deck[c_deck]; + s[2*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); + s[2*c + 2] = sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } } } else { THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE."); From cf9b9cdf20c0150059a944ae3f7a218b2d9ca568 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Mon, 27 Aug 2012 12:22:32 +0200 Subject: [PATCH 02/10] Started work on supporting 3 phases for wellreport. --- opm/core/utility/miscUtilities.cpp | 45 +++++++++++++++++++----------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 16c59b17e..cc7ff1276 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -574,9 +574,10 @@ namespace Opm { int nw = well_bhp.size(); ASSERT(nw == wells.number_of_wells); - if (props.numPhases() != 2) { - THROW("WellReport for now assumes two phase flow."); - } + int np = props.numPhases(); + //if (props.numPhases() != 2) { + // THROW("WellReport for now assumes two phase flow."); + //} const double* visc = props.viscosity(); std::vector data_now; data_now.reserve(1 + 3*nw); @@ -586,7 +587,8 @@ namespace Opm double well_rate_total = 0.0; double well_rate_water = 0.0; for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { - const double perf_rate = well_perfrates[perf]*(unit::day/unit::second); + const double perf_rate = unit::convert::to(well_perfrates[perf], + unit::cubic(unit::meter)/unit::day); well_rate_total += perf_rate; if (perf_rate > 0.0) { // Injection. @@ -594,11 +596,14 @@ namespace Opm } else { // Production. const int cell = wells.well_cells[perf]; - double mob[2]; + double mob[np]; props.relperm(1, &saturation[2*cell], &cell, mob, 0); - mob[0] /= visc[0]; - mob[1] /= visc[1]; - const double fracflow = mob[0]/(mob[0] + mob[1]); + double tmob=0; + for(int i=0; i < np; ++i){ + mob[i] /= visc[i]; + tmob += mob[i]; + } + const double fracflow = mob[0]/tmob; well_rate_water += perf_rate*fracflow; } } @@ -627,9 +632,10 @@ namespace Opm // TODO: refactor, since this is almost identical to the other push(). int nw = well_bhp.size(); ASSERT(nw == wells.number_of_wells); - if (props.numPhases() != 2) { - THROW("WellReport for now assumes two phase flow."); - } + int np = props.numPhases(); + //if (props.numPhases() != 2) { + // THROW("WellReport for now assumes two phase flow."); + //} std::vector data_now; data_now.reserve(1 + 3*nw); data_now.push_back(time/unit::day); @@ -640,20 +646,25 @@ namespace Opm for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { const double perf_rate = well_perfrates[perf]*(unit::day/unit::second); well_rate_total += perf_rate; - if (perf_rate > 0.0) { + if (perf_rate > 0.0) { // Injection. well_rate_water += perf_rate*wells.comp_frac[0]; } else { // Production. const int cell = wells.well_cells[perf]; - double mob[2]; + double mob[np]; props.relperm(1, &s[2*cell], &cell, mob, 0); - double visc[2]; + double visc[np]; props.viscosity(1, &p[cell], &z[2*cell], &cell, visc, 0); - mob[0] /= visc[0]; - mob[1] /= visc[1]; - const double fracflow = mob[0]/(mob[0] + mob[1]); + double tmob=0; + for(int i=0; i < np; ++i){ + mob[i] /= visc[i]; + tmob += mob[i]; + } + const double fracflow = mob[0]/(tmob); well_rate_water += perf_rate*fracflow; + //const double fracflow = mob[0]/(tmob); + //well_rate_water += perf_rate*fracflow; } } data_now.push_back(well_rate_total); From 94c04f343d5187701bc1461fa90825af69ca9849 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Tue, 28 Aug 2012 14:27:14 +0200 Subject: [PATCH 03/10] Added param initializer on fluids with param to addjust table length. Moved internal class SatFunc to SatFuncStone2. --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 16 ++++++++- opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 7 ++-- opm/core/fluid/SaturationPropsFromDeck.hpp | 33 +++++-------------- 3 files changed, 29 insertions(+), 27 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 004e6f88f..ae6a61306 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -18,7 +18,7 @@ */ #include - +#include namespace Opm { @@ -34,6 +34,20 @@ namespace Opm } } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const parameter::ParameterGroup& param) + { + rock_.init(deck, grid); + pvt_.init(deck); + satprops_.init(deck, grid, param); + + if (pvt_.numPhases() != satprops_.numPhases()) { + THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); + } + } + BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck() { } diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index a2a6c3de9..1d0a0d380 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -26,6 +26,7 @@ #include #include #include +#include struct UnstructuredGrid; @@ -43,8 +44,10 @@ namespace Opm /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const UnstructuredGrid& grid); - + const UnstructuredGrid& grid); + BlackoilPropertiesFromDeck(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const parameter::ParameterGroup& param); /// Destructor. virtual ~BlackoilPropertiesFromDeck(); diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index b83a12711..2c6a89a06 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -19,10 +19,11 @@ #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED - +#include #include #include #include +#include #include struct UnstructuredGrid; @@ -44,6 +45,10 @@ namespace Opm void init(const EclipseGridParser& deck, const UnstructuredGrid& grid); + void init(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const parameter::ParameterGroup& param); + /// \return P, the number of phases. int numPhases() const; @@ -87,31 +92,11 @@ namespace Opm double* smax) const; private: - PhaseUsage phase_usage_; - class SatFuncSet - { - public: - void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg); - void evalKr(const double* s, double* kr) const; - void evalKrDeriv(const double* s, double* kr, double* dkrds) const; - void evalPc(const double* s, double* pc) const; - void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - UniformTableLinear krw_; - UniformTableLinear krow_; - UniformTableLinear pcow_; - UniformTableLinear krg_; - UniformTableLinear krog_; - UniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) - }; - std::vector satfuncset_; + PhaseUsage phase_usage_; + std::vector satfuncset_; std::vector cell_to_func_; // = SATNUM - 1 - const SatFuncSet& funcForCell(const int cell) const; + const SatFuncStone2& funcForCell(const int cell) const; }; From 6852be422c458f6deb43ef6fc22eba7b561a02e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Tue, 28 Aug 2012 16:41:06 +0200 Subject: [PATCH 04/10] Added new SatFuncSimple fluid. Introduced a simple fluid which has no problem with strange black oil behavior. Intended for testing, but for now it is used in SaturationPropsFromDeck. --- opm/core/fluid/SaturationPropsFromDeck.hpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index 2c6a89a06..bd70d8029 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include struct UnstructuredGrid; @@ -92,11 +93,12 @@ namespace Opm double* smax) const; private: - PhaseUsage phase_usage_; - std::vector satfuncset_; + PhaseUsage phase_usage_; + typedef SatFuncSimple satfunc_t; + std::vector satfuncset_; std::vector cell_to_func_; // = SATNUM - 1 - const SatFuncStone2& funcForCell(const int cell) const; + const satfunc_t& funcForCell(const int cell) const; }; From 2dede29f20824b5c2f089b91950c7fe6d474c7b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Fri, 31 Aug 2012 17:01:07 +0200 Subject: [PATCH 05/10] Introduced posibility to change number of sample points for pvt. Did change the PVTW calculation so derivatives are exact. Extended the test functions for pvt and relperm --- opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 3 ++- opm/core/fluid/blackoil/BlackoilPvtProperties.hpp | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index ae6a61306..67943fc76 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -39,7 +39,8 @@ namespace Opm const parameter::ParameterGroup& param) { rock_.init(deck, grid); - pvt_.init(deck); + int samples = param.getDefault("dead_tab_size", 1025); + pvt_.init(deck, samples); satprops_.init(deck, grid, param); if (pvt_.numPhases() != satprops_.numPhases()) { diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index 627f266b7..99fd5b8f4 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp @@ -47,7 +47,7 @@ namespace Opm BlackoilPvtProperties(); /// Initialize from deck. - void init(const EclipseGridParser& deck); + void init(const EclipseGridParser& deck,const int samples = 16); /// Number of active phases. int numPhases() const; From 913054c47369279b1d074d71f08757fa48ee9536 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Sep 2012 13:54:50 +0200 Subject: [PATCH 06/10] Added more checks in 3-phase init code. --- opm/core/utility/initState_impl.hpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index f66704348..9a635e93f 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -512,11 +512,10 @@ namespace Opm State& state) { const int num_phases = props.numPhases(); - state.init(grid, num_phases); if (deck.hasField("EQUIL")) { if (num_phases != 2) { - THROW("initStateFromDeck(): currently handling only two-phase scenarios."); + THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios."); } // Set saturations depending on oil-water contact. const EQUIL& equil= deck.getEQUIL(); @@ -536,22 +535,27 @@ namespace Opm const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& p_deck = deck.getFloatingPointValue("PRESSURE"); const int num_cells = grid.number_of_cells; - if(num_phases == 2){ + if (num_phases == 2) { for (int c = 0; c < num_cells; ++c) { int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; s[2*c] = sw_deck[c_deck]; s[2*c + 1] = 1.0 - s[2*c]; p[c] = p_deck[c_deck]; } - }else{ + } else if (num_phases == 3) { + if (!deck.hasField("SGAS")) { + THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found)."); + } const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); - for (int c = 0; c < num_cells; ++c) { + for (int c = 0; c < num_cells; ++c) { int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; s[2*c] = sw_deck[c_deck]; s[2*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); s[2*c + 2] = sg_deck[c_deck]; p[c] = p_deck[c_deck]; } + } else { + THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases."); } } else { THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE."); From c2d41a6639c0f24663ccdf0c3a6ed38c92853d06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Sep 2012 14:04:52 +0200 Subject: [PATCH 07/10] Whitespace cleanup and adding a check for #phases <= 3. --- opm/core/utility/miscUtilities.cpp | 340 ++++++++++++++--------------- 1 file changed, 170 insertions(+), 170 deletions(-) diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 0b8d20103..e1658b83a 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -40,14 +40,14 @@ namespace Opm /// @param[out] porevol the pore volume by cell. void computePorevolume(const UnstructuredGrid& grid, const double* porosity, - std::vector& porevol) + std::vector& porevol) { - int num_cells = grid.number_of_cells; - porevol.resize(num_cells); - std::transform(porosity, porosity + num_cells, - grid.cell_volumes, - porevol.begin(), - std::multiplies()); + int num_cells = grid.number_of_cells; + porevol.resize(num_cells); + std::transform(porosity, porosity + num_cells, + grid.cell_volumes, + porevol.begin(), + std::multiplies()); } @@ -98,20 +98,20 @@ namespace Opm /// For each phase p, we compute /// sat_vol_p = sum_i s_p_i pv_i void computeSaturatedVol(const std::vector& pv, - const std::vector& s, - double* sat_vol) + const std::vector& s, + double* sat_vol) { - const int num_cells = pv.size(); - const int np = s.size()/pv.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and pv vectors do not match."); - } - std::fill(sat_vol, sat_vol + np, 0.0); - for (int c = 0; c < num_cells; ++c) { - for (int p = 0; p < np; ++p) { - sat_vol[p] += pv[c]*s[np*c + p]; - } - } + const int num_cells = pv.size(); + const int np = s.size()/pv.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and pv vectors do not match."); + } + std::fill(sat_vol, sat_vol + np, 0.0); + for (int c = 0; c < num_cells; ++c) { + for (int p = 0; p < np; ++p) { + sat_vol[p] += pv[c]*s[np*c + p]; + } + } } @@ -123,28 +123,28 @@ namespace Opm /// For each phase p, we compute /// aver_sat_p = (sum_i s_p_i pv_i) / (sum_i pv_i). void computeAverageSat(const std::vector& pv, - const std::vector& s, - double* aver_sat) + const std::vector& s, + double* aver_sat) { - const int num_cells = pv.size(); - const int np = s.size()/pv.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and pv vectors do not match."); - } - double tot_pv = 0.0; - // Note that we abuse the output array to accumulate the - // saturated pore volumes. - std::fill(aver_sat, aver_sat + np, 0.0); - for (int c = 0; c < num_cells; ++c) { - tot_pv += pv[c]; - for (int p = 0; p < np; ++p) { - aver_sat[p] += pv[c]*s[np*c + p]; - } - } - // Must divide by pore volumes to get saturations. - for (int p = 0; p < np; ++p) { - aver_sat[p] /= tot_pv; - } + const int num_cells = pv.size(); + const int np = s.size()/pv.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and pv vectors do not match."); + } + double tot_pv = 0.0; + // Note that we abuse the output array to accumulate the + // saturated pore volumes. + std::fill(aver_sat, aver_sat + np, 0.0); + for (int c = 0; c < num_cells; ++c) { + tot_pv += pv[c]; + for (int p = 0; p < np; ++p) { + aver_sat[p] += pv[c]*s[np*c + p]; + } + } + // Must divide by pore volumes to get saturations. + for (int p = 0; p < np; ++p) { + aver_sat[p] /= tot_pv; + } } @@ -161,38 +161,38 @@ namespace Opm /// where P = s.size()/src.size(). /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const IncompPropertiesInterface& props, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced) + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); - } - std::fill(injected, injected + np, 0.0); - std::fill(produced, produced + np, 0.0); - const double* visc = props.viscosity(); - std::vector mob(np); - for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; - const double* sat = &s[np*c]; - props.relperm(1, sat, &c, &mob[0], 0); - double totmob = 0.0; - for (int p = 0; p < np; ++p) { - mob[p] /= visc[p]; - totmob += mob[p]; - } - for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; - } - } - } + const int num_cells = src.size(); + const int np = s.size()/src.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and src vectors do not match."); + } + std::fill(injected, injected + np, 0.0); + std::fill(produced, produced + np, 0.0); + const double* visc = props.viscosity(); + std::vector mob(np); + for (int c = 0; c < num_cells; ++c) { + if (src[c] > 0.0) { + injected[0] += src[c]*dt; + } else if (src[c] < 0.0) { + const double flux = -src[c]*dt; + const double* sat = &s[np*c]; + props.relperm(1, sat, &c, &mob[0], 0); + double totmob = 0.0; + for (int p = 0; p < np; ++p) { + mob[p] /= visc[p]; + totmob += mob[p]; + } + for (int p = 0; p < np; ++p) { + produced[p] += (mob[p]/totmob)*flux; + } + } + } } @@ -203,9 +203,9 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::IncompPropertiesInterface& props, - const std::vector& cells, - const std::vector& s, - std::vector& totmob) + const std::vector& cells, + const std::vector& s, + std::vector& totmob) { std::vector pmobc; @@ -231,10 +231,10 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::IncompPropertiesInterface& props, - const std::vector& cells, - const std::vector& s, - std::vector& totmob, - std::vector& omega) + const std::vector& cells, + const std::vector& s, + std::vector& totmob, + std::vector& omega) { std::vector pmobc; @@ -331,33 +331,33 @@ namespace Opm /// (+) positive inflow of first phase (water) /// (-) negative total outflow of both phases void computeTransportSource(const UnstructuredGrid& grid, - const std::vector& src, - const std::vector& faceflux, - const double inflow_frac, + const std::vector& src, + const std::vector& faceflux, + const double inflow_frac, const Wells* wells, const std::vector& well_perfrates, - std::vector& transport_src) + std::vector& transport_src) { - int nc = grid.number_of_cells; - transport_src.resize(nc); + int nc = grid.number_of_cells; + transport_src.resize(nc); // Source term and boundary contributions. - for (int c = 0; c < nc; ++c) { - transport_src[c] = 0.0; - transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; - for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) { - int f = grid.cell_faces[hf]; - const int* f2c = &grid.face_cells[2*f]; - double bdy_influx = 0.0; - if (f2c[0] == c && f2c[1] == -1) { - bdy_influx = -faceflux[f]; - } else if (f2c[0] == -1 && f2c[1] == c) { - bdy_influx = faceflux[f]; - } - if (bdy_influx != 0.0) { - transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; - } - } - } + for (int c = 0; c < nc; ++c) { + transport_src[c] = 0.0; + transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; + for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* f2c = &grid.face_cells[2*f]; + double bdy_influx = 0.0; + if (f2c[0] == c && f2c[1] == -1) { + bdy_influx = -faceflux[f]; + } else if (f2c[0] == -1 && f2c[1] == c) { + bdy_influx = faceflux[f]; + } + if (bdy_influx != 0.0) { + transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; + } + } + } // Well contributions. if (wells) { @@ -392,52 +392,52 @@ namespace Opm /// @param[in] face_flux signed per-face fluxes /// @param[out] cell_velocity the estimated velocities. void estimateCellVelocity(const UnstructuredGrid& grid, - const std::vector& face_flux, - std::vector& cell_velocity) + const std::vector& face_flux, + std::vector& cell_velocity) { - const int dim = grid.dimensions; - cell_velocity.clear(); - cell_velocity.resize(grid.number_of_cells*dim, 0.0); - for (int face = 0; face < grid.number_of_faces; ++face) { - int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; - const double* fc = &grid.face_centroids[face*dim]; - double flux = face_flux[face]; - for (int i = 0; i < 2; ++i) { - if (c[i] >= 0) { - const double* cc = &grid.cell_centroids[c[i]*dim]; - for (int d = 0; d < dim; ++d) { - double v_contrib = fc[d] - cc[d]; - v_contrib *= flux/grid.cell_volumes[c[i]]; - cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; - } - } - } - } + const int dim = grid.dimensions; + cell_velocity.clear(); + cell_velocity.resize(grid.number_of_cells*dim, 0.0); + for (int face = 0; face < grid.number_of_faces; ++face) { + int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; + const double* fc = &grid.face_centroids[face*dim]; + double flux = face_flux[face]; + for (int i = 0; i < 2; ++i) { + if (c[i] >= 0) { + const double* cc = &grid.cell_centroids[c[i]*dim]; + for (int d = 0; d < dim; ++d) { + double v_contrib = fc[d] - cc[d]; + v_contrib *= flux/grid.cell_volumes[c[i]]; + cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; + } + } + } + } } /// Extract a vector of water saturations from a vector of /// interleaved water and oil saturations. void toWaterSat(const std::vector& sboth, - std::vector& sw) + std::vector& sw) { - int num = sboth.size()/2; - sw.resize(num); - for (int i = 0; i < num; ++i) { - sw[i] = sboth[2*i]; - } + int num = sboth.size()/2; + sw.resize(num); + for (int i = 0; i < num; ++i) { + sw[i] = sboth[2*i]; + } } /// Make a a vector of interleaved water and oil saturations from /// a vector of water saturations. void toBothSat(const std::vector& sw, - std::vector& sboth) + std::vector& sboth) { - int num = sw.size(); - sboth.resize(2*num); - for (int i = 0; i < num; ++i) { - sboth[2*i] = sw[i]; - sboth[2*i + 1] = 1.0 - sw[i]; - } + int num = sw.size(); + sboth.resize(2*num); + for (int i = 0; i < num; ++i) { + sboth[2*i] = sw[i]; + sboth[2*i + 1] = 1.0 - sw[i]; + } } @@ -450,30 +450,30 @@ namespace Opm if (np != 2) { THROW("wellsToSrc() requires a 2 phase case."); } - src.resize(num_cells); - for (int w = 0; w < wells.number_of_wells; ++w) { + src.resize(num_cells); + for (int w = 0; w < wells.number_of_wells; ++w) { const int cur = wells.ctrls[w]->current; - if (wells.ctrls[w]->num != 1) { - MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored."); - } - if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) { - THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE."); - } - if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) { - THROW("In wellsToSrc(): well has multiple perforations."); - } + if (wells.ctrls[w]->num != 1) { + MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored."); + } + if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) { + THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE."); + } + if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) { + THROW("In wellsToSrc(): well has multiple perforations."); + } for (int p = 0; p < np; ++p) { if (wells.ctrls[w]->distr[np*cur + p] != 1.0) { THROW("In wellsToSrc(): well not controlled on total rate."); } } - double flow = wells.ctrls[w]->target[cur]; + double flow = wells.ctrls[w]->target[cur]; if (wells.type[w] == INJECTOR) { flow *= wells.comp_frac[np*w + 0]; // Obtaining water rate for inflow source. } - const double cell = wells.well_cells[wells.well_connpos[w]]; - src[cell] = flow; - } + const double cell = wells.well_cells[wells.well_connpos[w]]; + src[cell] = flow; + } } @@ -594,9 +594,10 @@ namespace Opm int nw = well_bhp.size(); ASSERT(nw == wells.number_of_wells); int np = props.numPhases(); - //if (props.numPhases() != 2) { - // THROW("WellReport for now assumes two phase flow."); - //} + const int max_np = 3; + if (np > max_np) { + THROW("WellReport for now assumes #phases <= " << max_np); + } const double* visc = props.viscosity(); std::vector data_now; data_now.reserve(1 + 3*nw); @@ -606,8 +607,8 @@ namespace Opm double well_rate_total = 0.0; double well_rate_water = 0.0; for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { - const double perf_rate = unit::convert::to(well_perfrates[perf], - unit::cubic(unit::meter)/unit::day); + const double perf_rate = unit::convert::to(well_perfrates[perf], + unit::cubic(unit::meter)/unit::day); well_rate_total += perf_rate; if (perf_rate > 0.0) { // Injection. @@ -615,13 +616,13 @@ namespace Opm } else { // Production. const int cell = wells.well_cells[perf]; - double mob[np]; + double mob[max_np]; props.relperm(1, &saturation[2*cell], &cell, mob, 0); - double tmob=0; - for(int i=0; i < np; ++i){ + double tmob = 0; + for(int i = 0; i < np; ++i) { mob[i] /= visc[i]; tmob += mob[i]; - } + } const double fracflow = mob[0]/tmob; well_rate_water += perf_rate*fracflow; } @@ -652,9 +653,10 @@ namespace Opm int nw = well_bhp.size(); ASSERT(nw == wells.number_of_wells); int np = props.numPhases(); - //if (props.numPhases() != 2) { - // THROW("WellReport for now assumes two phase flow."); - //} + const int max_np = 3; + if (np > max_np) { + THROW("WellReport for now assumes #phases <= " << max_np); + } std::vector data_now; data_now.reserve(1 + 3*nw); data_now.push_back(time/unit::day); @@ -665,25 +667,23 @@ namespace Opm for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { const double perf_rate = well_perfrates[perf]*(unit::day/unit::second); well_rate_total += perf_rate; - if (perf_rate > 0.0) { + if (perf_rate > 0.0) { // Injection. well_rate_water += perf_rate*wells.comp_frac[0]; } else { // Production. const int cell = wells.well_cells[perf]; - double mob[np]; + double mob[max_np]; props.relperm(1, &s[2*cell], &cell, mob, 0); - double visc[np]; + double visc[max_np]; props.viscosity(1, &p[cell], &z[2*cell], &cell, visc, 0); - double tmob=0; - for(int i=0; i < np; ++i){ + double tmob = 0; + for(int i = 0; i < np; ++i) { mob[i] /= visc[i]; tmob += mob[i]; - } + } const double fracflow = mob[0]/(tmob); well_rate_water += perf_rate*fracflow; - //const double fracflow = mob[0]/(tmob); - //well_rate_water += perf_rate*fracflow; } } data_now.push_back(well_rate_total); From 489501b49b3669548a3e6657e51de156692f87a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Sep 2012 14:16:54 +0200 Subject: [PATCH 08/10] Documented new constructor. Fixed layout. --- opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index 1d0a0d380..e6b9c5b41 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -45,9 +45,18 @@ namespace Opm /// to logical cartesian indices consistent with the deck. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid); + + /// Initialize from deck, grid and parameters. + /// \param deck Deck input parser + /// \param grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + /// \param param Parameters. Accepted parameters include: + /// tab_size (200) number of uniform sample points for saturation tables. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, const parameter::ParameterGroup& param); + /// Destructor. virtual ~BlackoilPropertiesFromDeck(); From 4e1647bb6272fa98ba57936d6761b8b39578f90b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Sep 2012 15:07:03 +0200 Subject: [PATCH 09/10] Formatting fixes. --- opm/core/fluid/blackoil/BlackoilPvtProperties.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index 99fd5b8f4..191c6fa89 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp @@ -47,7 +47,7 @@ namespace Opm BlackoilPvtProperties(); /// Initialize from deck. - void init(const EclipseGridParser& deck,const int samples = 16); + void init(const EclipseGridParser& deck, const int samples = 16); /// Number of active phases. int numPhases() const; From e9c4c2499cc362facef490fae23178c2b72fe6b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Sep 2012 15:09:55 +0200 Subject: [PATCH 10/10] Documented parameters. --- opm/core/fluid/BlackoilPropertiesFromDeck.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index e6b9c5b41..b63c0af88 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -52,7 +52,8 @@ namespace Opm /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. /// \param param Parameters. Accepted parameters include: - /// tab_size (200) number of uniform sample points for saturation tables. + /// dead_tab_size (1025) number of uniform sample points for dead-oil pvt tables. + /// tab_size_kr (200) number of uniform sample points for saturation tables. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, const parameter::ParameterGroup& param);