From f41b82da3d7db0bc50a12787fa32baa900ce604b Mon Sep 17 00:00:00 2001 From: osae Date: Thu, 29 Nov 2012 16:36:13 +0100 Subject: [PATCH 1/5] Scaling of relperm functions - oil/water systems. opm/core/eclipse/EclipseGridParser.cpp opm/core/eclipse/EclipseGridParser.hpp - New keywords: ENDSCALE SCALECRS SWCR SWL SWU SOWCR KRW KRWR KRO KRORW opm/core/eclipse/SpecialEclipseFields.hpp - Parsers for ENDSCALE and SCALECRS. opm/core/fluid/BlackoilPropertiesFromDeck.cpp - Consistency check: ENDSCALE implemented for SatFuncSimple only. opm/core/fluid/SatFuncGwseg.hpp opm/core/fluid/SatFuncSimple.hpp opm/core/fluid/SatFuncStone2.hpp - Accomodate "default" values for scalable parameters. - For SatFuncGwseg and SatFuncStone2 the associated functionality not yet supported and the variables are dummies to satisfy the compiler. opm/core/fluid/SatFuncSimple.cpp - Initialisation for scalable parameters. - Evaluation of relperms: Use (1-so) for evaluation of oil-relperms. (For scaled arguments sw and so do not necessarily add to one.) - TODO: SatFuncGwseg.cpp and SatFuncStone2.cpp for oil-water systems. opm/core/fluid/SaturationPropsFromDeck.hpp - Struct to accomodate cell-wise scaling factors. - Two flags indicating scaling and method. - Methods for parameter initialisation and scaled relperm computation. opm/core/fluid/SaturationPropsFromDeck_impl.hpp - Initialize scaling options and relevant cell-wise scaling factors. - Relperm evaluation modified for possible end point scaling. --- opm/core/eclipse/EclipseGridParser.cpp | 11 +- opm/core/eclipse/EclipseGridParser.hpp | 2 + opm/core/eclipse/SpecialEclipseFields.hpp | 103 +++++++ opm/core/fluid/BlackoilPropertiesFromDeck.cpp | 3 + opm/core/fluid/SatFuncGwseg.hpp | 12 + opm/core/fluid/SatFuncSimple.cpp | 59 +++- opm/core/fluid/SatFuncSimple.hpp | 12 + opm/core/fluid/SatFuncStone2.hpp | 12 + opm/core/fluid/SaturationPropsFromDeck.hpp | 19 ++ .../fluid/SaturationPropsFromDeck_impl.hpp | 283 +++++++++++++++++- 10 files changed, 505 insertions(+), 11 deletions(-) diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index 9c0f437d..fef38c44 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -97,7 +97,10 @@ namespace EclipseKeywords string("MULTPV"), string("PRESSURE"), string("SGAS"), string("SWAT"), string("SOIL"), string("RS"), string("DXV"), string("DYV"), string("DZV"), - string("DEPTHZ"), string("TOPS"), string("MAPAXES") + string("DEPTHZ"), string("TOPS"), string("MAPAXES"), + string("SWCR"), string("SWL"), string("SWU"), + string("SOWCR"), string("KRW"), string("KRWR"), + string("KRO"), string("KRORW") }; const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); @@ -114,7 +117,7 @@ namespace EclipseKeywords string("PLYVISC"), string("PLYROCK"), string("PLYADS"), string("PLYMAX"), string("TLMIXPAR"), string("WPOLYMER"), string("GRUPTREE"), string("GCONINJE"), string("GCONPROD"), - string("WGRUPCON"), + string("WGRUPCON"), string("ENDSCALE"), string("SCALECRS"), // The following fields only have a dummy implementation // that allows us to ignore them. string("SWFN"), @@ -554,7 +557,9 @@ void EclipseGridParser::convertToSI() key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" || key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" || key == "SGAS" || key == "SWAT" || key == "SOIL" || - key == "RS") { + key == "RS" || key == "SWCR" || key == "SWL" || + key == "SWU" || key == "SOWCR" || key == "KRW" || + key == "KRWR" || key == "KRORW" || key == "KRO") { unit = 1.0; do_convert = false; // Dimensionless keywords... } else if (key == "PRESSURE") { diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index e38f8b28..ab871c68 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -192,6 +192,8 @@ public: SPECIAL_FIELD(GCONINJE) SPECIAL_FIELD(GCONPROD) SPECIAL_FIELD(WGRUPCON) + SPECIAL_FIELD(ENDSCALE) + SPECIAL_FIELD(SCALECRS) // The following fields only have a dummy implementation // that allows us to ignore them. diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index 507073d3..70bff0a0 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -2206,6 +2206,109 @@ struct WPOLYMER : public SpecialBase } }; +struct ENDSCALE : public SpecialBase { + std::string dir_switch_; + std::string revers_switch_; + int n_tab_; + int n_node_; + + virtual std::string name() const { + return std::string("ENDSCALE"); + } + + virtual void read(std::istream & is) { + dir_switch_ = std::string("NODIR"); + revers_switch_ = std::string("REVERS"); + n_tab_ = 1; + n_node_ = 20; + while (!is.eof()) { + std::string tmp = readString(is); + while (tmp.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + tmp = readString(is); + } + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } + dir_switch_ = tmp; + + tmp = readString(is); + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } + revers_switch_ = tmp; + + is >> ignoreWhitespace; + if(is.peek() == int('/')) { + is >> ignoreLine; + return; + } + is >> n_tab_; + + is >> ignoreWhitespace; + if(is.peek() == int('/')) { + is >> ignoreLine; + return; + } + is >> n_node_; + + is >> ignoreLine; + break; + } + } + + virtual void write(std::ostream & os) { + os << name() << '\n'; + os << dir_switch_ << " " << revers_switch_ << " " + << n_tab_ << " " << n_node_ << '\n'; + } + + virtual void convertToSI(const EclipseUnits&) { + } +}; + + +struct SCALECRS : public SpecialBase { + std::string scalecrs_; + + virtual std::string name() const { + return std::string("SCALECRS"); + } + + virtual void read(std::istream & is) { + scalecrs_ = std::string("NO"); + while (!is.eof()) { + std::string tmp = readString(is); + while (tmp.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + tmp = readString(is); + } + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } else if (tmp[0] == 'Y') { + scalecrs_ = std::string("YES"); + } + + is >> ignoreLine; + break; + } + } + + virtual void write(std::ostream & os) { + os << name() << '\n'; + os << scalecrs_ << '\n'; + } + + virtual void convertToSI(const EclipseUnits&) { + } +}; + + // The following fields only have a dummy implementation // that allows us to ignore them. struct SWFN : public MultRec {}; diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 271f9395..d49dc5b9 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -50,6 +50,9 @@ namespace Opm // Unfortunate lack of pointer smartness here... const int sat_samples = param.getDefault("sat_tab_size", 200); std::string threephase_model = param.getDefault("threephase_model", "simple"); + if (deck.hasField("ENDSCALE") && threephase_model != "simple") { + THROW("Sorry, end point scaling currently available for the 'simple' model only."); + } if (sat_samples > 1) { if (threephase_model == "stone2") { SaturationPropsFromDeck* ptr diff --git a/opm/core/fluid/SatFuncGwseg.hpp b/opm/core/fluid/SatFuncGwseg.hpp index 1932e7f0..ae5b47a6 100644 --- a/opm/core/fluid/SatFuncGwseg.hpp +++ b/opm/core/fluid/SatFuncGwseg.hpp @@ -40,6 +40,12 @@ namespace Opm void evalPcDeriv(const double* s, double* pc, double* dpcds) const; double smin_[PhaseUsage::MaxNumPhases]; double smax_[PhaseUsage::MaxNumPhases]; + double krwmax_; // Max water relperm + double kromax_; // Max oil relperm + double swcr_; // Critical water saturation. + double krwr_; // Water relperm at critical oil-in-water saturation. + double sowcr_; // Critical oil-in-water saturation. + double krorw_; // Oil relperm at critical water saturation. private: PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. UniformTableLinear krw_; @@ -65,6 +71,12 @@ namespace Opm void evalPcDeriv(const double* s, double* pc, double* dpcds) const; double smin_[PhaseUsage::MaxNumPhases]; double smax_[PhaseUsage::MaxNumPhases]; + double krwmax_; // Max water relperm + double kromax_; // Max oil relperm + double swcr_; // Critical water saturation. + double krwr_; // Water relperm at critical oil-in-water saturation. + double sowcr_; // Critical oil-in-water saturation. + double krorw_; // Oil relperm at critical water saturation. private: PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. NonuniformTableLinear krw_; diff --git a/opm/core/fluid/SatFuncSimple.cpp b/opm/core/fluid/SatFuncSimple.cpp index 7e226a4a..37e74fff 100644 --- a/opm/core/fluid/SatFuncSimple.cpp +++ b/opm/core/fluid/SatFuncSimple.cpp @@ -54,6 +54,27 @@ namespace Opm smin_[phase_usage.phase_pos[Aqua]] = sw[0]; swmax = sw.back(); smax_[phase_usage.phase_pos[Aqua]] = sw.back(); + + krwmax_ = krw.back(); + kromax_ = krow.front(); + swcr_ = swmax; + sowcr_ = 1.0 - swco; + krwr_ = krw.back(); + krorw_ = krow.front(); + for (unsigned int i=1; i 0.0) { + swcr_ = sw[i-1]; + krorw_ = krow[i-1]; + break; + } + } + for (unsigned int i=sw.size()-2; i>=0; --i) { + if (krow[i]> 0.0) { + sowcr_ = 1.0 - sw[i+1]; + krwr_ = krw[i+1]; + break; + } + } } if (phase_usage.phase_used[Vapour]) { const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; @@ -102,7 +123,8 @@ namespace Opm int opos = phase_usage.phase_pos[Liquid]; double sw = s[wpos]; double krw = krw_(sw); - double krow = krow_(sw); + double so = s[opos]; + double krow = krow_(1.0-so); kr[wpos] = krw; kr[opos] = krow; } else { @@ -160,8 +182,9 @@ namespace Opm double sw = s[wpos]; double krw = krw_(sw); double dkrww = krw_.derivative(sw); - double krow = krow_(sw); - double dkrow = krow_.derivative(sw); + double so = s[opos]; + double krow = krow_(1.0-so); + double dkrow = krow_.derivative(1.0-so); kr[wpos] = krw; kr[opos] = krow; dkrds[wpos + wpos*np] = dkrww; @@ -253,6 +276,28 @@ namespace Opm smin_[phase_usage.phase_pos[Aqua]] = sw[0]; swmax = sw.back(); smax_[phase_usage.phase_pos[Aqua]] = sw.back(); + + krwmax_ = krw.back(); + kromax_ = krow.front(); + swcr_ = swmax; + sowcr_ = 1.0 - swco; + krwr_ = krw.back(); + krorw_ = krow.front(); + for (unsigned int i=1; i 0.0) { + swcr_ = sw[i-1]; + krorw_ = krow[i-1]; + break; + } + } + for (unsigned int i=sw.size()-2; i>=0; --i) { + if (krow[i]> 0.0) { + sowcr_ = 1.0 - sw[i+1]; + krwr_ = krw[i+1]; + break; + } + } + } if (phase_usage.phase_used[Vapour]) { const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; @@ -301,7 +346,8 @@ namespace Opm int opos = phase_usage.phase_pos[Liquid]; double sw = s[wpos]; double krw = krw_(sw); - double krow = krow_(sw); + double so = s[opos]; + double krow = krow_(1.0-so); kr[wpos] = krw; kr[opos] = krow; } else { @@ -359,8 +405,9 @@ namespace Opm double sw = s[wpos]; double krw = krw_(sw); double dkrww = krw_.derivative(sw); - double krow = krow_(sw); - double dkrow = krow_.derivative(sw); + double so = s[opos]; + double krow = krow_(1.0-so); + double dkrow = krow_.derivative(1.0-so); kr[wpos] = krw; kr[opos] = krow; dkrds[wpos + wpos*np] = dkrww; diff --git a/opm/core/fluid/SatFuncSimple.hpp b/opm/core/fluid/SatFuncSimple.hpp index 338a359e..1b6ba03b 100644 --- a/opm/core/fluid/SatFuncSimple.hpp +++ b/opm/core/fluid/SatFuncSimple.hpp @@ -40,6 +40,12 @@ namespace Opm void evalPcDeriv(const double* s, double* pc, double* dpcds) const; double smin_[PhaseUsage::MaxNumPhases]; double smax_[PhaseUsage::MaxNumPhases]; + double krwmax_; // Max water relperm + double kromax_; // Max oil relperm + double swcr_; // Critical water saturation. + double krwr_; // Water relperm at critical oil-in-water saturation. + double sowcr_; // Critical oil-in-water saturation. + double krorw_; // Oil relperm at critical water saturation. private: PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. UniformTableLinear krw_; @@ -65,6 +71,12 @@ namespace Opm void evalPcDeriv(const double* s, double* pc, double* dpcds) const; double smin_[PhaseUsage::MaxNumPhases]; double smax_[PhaseUsage::MaxNumPhases]; + double krwmax_; // Max water relperm + double kromax_; // Max oil relperm + double swcr_; // Critical water saturation. + double krwr_; // Water relperm at critical oil-in-water saturation. + double sowcr_; // Critical oil-in-water saturation. + double krorw_; // Oil relperm at critical water saturation. private: PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. NonuniformTableLinear krw_; diff --git a/opm/core/fluid/SatFuncStone2.hpp b/opm/core/fluid/SatFuncStone2.hpp index 532bf7ce..fc906126 100644 --- a/opm/core/fluid/SatFuncStone2.hpp +++ b/opm/core/fluid/SatFuncStone2.hpp @@ -40,6 +40,12 @@ namespace Opm void evalPcDeriv(const double* s, double* pc, double* dpcds) const; double smin_[PhaseUsage::MaxNumPhases]; double smax_[PhaseUsage::MaxNumPhases]; + double krwmax_; // Max water relperm + double kromax_; // Max oil relperm + double swcr_; // Critical water saturation. + double krwr_; // Water relperm at critical oil-in-water saturation. + double sowcr_; // Critical oil-in-water saturation. + double krorw_; // Oil relperm at critical water saturation. private: PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. UniformTableLinear krw_; @@ -65,6 +71,12 @@ namespace Opm void evalPcDeriv(const double* s, double* pc, double* dpcds) const; double smin_[PhaseUsage::MaxNumPhases]; double smax_[PhaseUsage::MaxNumPhases]; + double krwmax_; // Max water relperm + double kromax_; // Max oil relperm + double swcr_; // Critical water saturation. + double krwr_; // Water relperm at critical oil-in-water saturation. + double sowcr_; // Critical oil-in-water saturation. + double krorw_; // Oil relperm at critical water saturation. private: PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. NonuniformTableLinear krw_; diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index ea2acb34..1492c5a5 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -107,9 +107,28 @@ namespace Opm std::vector satfuncset_; std::vector cell_to_func_; // = SATNUM - 1 + struct { // End point scaling parameters + std::vector swl_; + std::vector swcr_; + std::vector swu_; + std::vector sowcr_; + std::vector krw_; + std::vector krwr_; + std::vector kro_; + std::vector krorw_; + } eps_; + bool do_eps_; // ENDSCALE is active + bool do_3pt_; // SCALECRS: YES~true NO~false + typedef SatFuncSet Funcs; const Funcs& funcForCell(const int cell) const; + + void initEPS(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam); + void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const; }; diff --git a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp index f7a80b45..6b87bc29 100644 --- a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp @@ -96,6 +96,61 @@ namespace Opm for (int table = 0; table < num_tables; ++table) { satfuncset_[table].init(deck, table, phase_usage_, samples); } + + // Saturation table scaling + do_eps_ = false; + do_3pt_ = false; + if (deck.hasField("ENDSCALE")) { + if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) { + THROW("Currently endpoint-scaling limited to oil-water systems without gas."); + } + if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) { + THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted."); + } + if (deck.getENDSCALE().revers_switch_ != std::string("REVERS")) { + THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'REVERS' accepted."); + } + if (deck.hasField("SCALECRS")) { + if (deck.getSCALECRS().scalecrs_ == std::string("YES")) { + do_3pt_ = true; + } + } + do_eps_ = true; + initEPS(deck, grid, std::string("SWCR"), eps_.swcr_); + initEPS(deck, grid, std::string("SWL"), eps_.swl_); + initEPS(deck, grid, std::string("SWU"), eps_.swu_); + initEPS(deck, grid, std::string("SOWCR"), eps_.sowcr_); + initEPS(deck, grid, std::string("KRW"), eps_.krw_); + initEPS(deck, grid, std::string("KRWR"), eps_.krwr_); + initEPS(deck, grid, std::string("KRO"), eps_.kro_); + initEPS(deck, grid, std::string("KRORW"), eps_.krorw_); + +/* + double ss[PhaseUsage::MaxNumPhases], kr[PhaseUsage::MaxNumPhases]; + int oldP = std::cout.precision(); + std::cout.precision(4); + for (unsigned int i=0; i<=100; ++i) { + ss[phase_usage_.phase_pos[Aqua]] = i*0.01; + ss[phase_usage_.phase_pos[Liquid]] = 1.0 - ss[phase_usage_.phase_pos[Aqua]]; + endScaling(ss, 15, kr); + std::cout << std::showpoint + << std::setw(10) << ss[phase_usage_.phase_pos[Aqua]] + << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]]; + endScaling(ss, 45, kr); + std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]]; + endScaling(ss, 75, kr); + std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]]; + endScaling(ss, 105, kr); + std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]] + << std::noshowpoint << std::endl; + } + std::cout.precision(oldP); +*/ + } } @@ -134,12 +189,20 @@ namespace Opm if (dkrds) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); + if (do_eps_) { + relpermEPS(s + np*i, cells[i], kr + np*i, dkrds + np*np*i); + } else { + funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); + } } } else { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); + if (do_eps_) { + relpermEPS(s + np*i, cells[i], kr + np*i); + } else { + funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); + } } } } @@ -214,7 +277,223 @@ namespace Opm return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; } + // Initialize saturation scaling parameter + template + void SaturationPropsFromDeck::initEPS(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam) + { + + if (deck.hasField(keyword)) { + // Active keyword assigned default values for each cell (in case of possible box-wise assignment) + scaleparam.resize(grid.number_of_cells); + int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + if (keyword == std::string("SWCR")) { + for (int i=0; i& val = deck.getFloatingPointValue(keyword); + for (int c = 0; c < int(scaleparam.size()); ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + scaleparam[c] = val[deck_pos]; + } + } + } + + // Saturation scaling + template + void SaturationPropsFromDeck::relpermEPS(const double *s, const int cell, double *kr, double *dkrds) const + { + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; + double ss[PhaseUsage::MaxNumPhases]; + + if (do_3pt_) { // Three-point scaling + // Transforms for water saturation + if (eps_.swcr_.empty() && eps_.swu_.empty()) { + ss[wpos] = s[wpos]; + } else { + double s_r = 1.0-funcForCell(cell).sowcr_; + double sr = eps_.sowcr_.empty() ? s_r : 1.0-eps_.sowcr_[cell]; + if (s[wpos] <= sr) { + double sw_cr = funcForCell(cell).swcr_; + double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell]; + ss[wpos] = (s[wpos] <= swcr) ? sw_cr : sw_cr+(s[wpos]-swcr)*(s_r-sw_cr)/(sr-swcr); + } else { + double sw_max = funcForCell(cell).smax_[wpos]; + double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell]; + ss[wpos] = (s[wpos] >= swmax) ? sw_max : s_r+(s[wpos]-sr)*(sw_max-s_r)/(swmax-sr); + } + } + // Transforms for oil saturation + if (eps_.sowcr_.empty() && eps_.swl_.empty()) { + ss[opos] = s[opos]; + } else { + double s_r = 1.0-funcForCell(cell).swcr_; + double sr = eps_.swcr_.empty() ? s_r : 1.0-eps_.swcr_[cell]; + if (s[opos] <= sr) { + double sow_cr = funcForCell(cell).sowcr_; + double sowcr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell]; + ss[opos] = (s[opos] <= sowcr) ? sow_cr : sow_cr+(s[opos]-sowcr)*(s_r-sow_cr)/(sr-sowcr); + } else { + double sow_max = funcForCell(cell).smax_[opos]; + double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]); + ss[opos] = (s[opos] >= sowmax) ? sow_max : s_r+(s[opos]-sr)*(sow_max-s_r)/(sowmax-sr); + } + } + } else { // Two-point scaling + // Transforms for water saturation + if (eps_.swcr_.empty() && eps_.swu_.empty()) { + ss[wpos] = s[wpos]; + } else { + double sw_cr = funcForCell(cell).swcr_; + double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell]; + if (s[wpos] <= swcr) { + ss[wpos] = sw_cr; + } else { + double sw_max = funcForCell(cell).smax_[wpos]; + double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell]; + ss[wpos] = (s[wpos] >= swmax) ? sw_max : sw_cr + (s[wpos]-swcr)*(sw_max-sw_cr)/(swmax-swcr); + } + } + // Transforms for oil saturation + if (eps_.sowcr_.empty() && eps_.swl_.empty()) { + ss[opos] = s[opos]; + } else { + double sow_cr = funcForCell(cell).sowcr_; + double socr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell]; + if (s[opos] <= socr) { + ss[opos] = sow_cr; + } else { + double sow_max = funcForCell(cell).smax_[opos]; + double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]); + ss[opos] = (s[opos] >= sowmax) ? sow_max : sow_cr + (s[opos]-socr) *(sow_max-sow_cr)/(sowmax-socr); + } + } + } + + // Evaluation of relperms + if (dkrds) { + THROW("Relperm derivatives not yet available in combination with end point scaling ..."); + funcForCell(cell).evalKrDeriv(ss, kr, dkrds); + } else { + // Assume: sw_cr -> krw=0 sw_max -> krw= + // sow_cr -> kro=0 sow_max -> kro= + funcForCell(cell).evalKr(ss, kr); + } + + // Scaling of relperms values + // - Water + if (eps_.krw_.empty() && eps_.krwr_.empty()) { // No value scaling + } else if (eps_.krwr_.empty()) { // Two-point + kr[wpos] *= (eps_.krw_[cell]/funcForCell(cell).krwmax_); + } else { + double swcr = eps_.swcr_.empty() ? funcForCell(cell).swcr_ : eps_.swcr_[cell]; + double swmax = eps_.swu_.empty() ? funcForCell(cell).smax_[wpos] : eps_.swu_[cell]; + double sr; + if (do_3pt_) { + sr = eps_.sowcr_.empty() ? 1.0-funcForCell(cell).sowcr_ : 1.0-eps_.sowcr_[cell]; + } else { + double sw_cr = funcForCell(cell).swcr_; + double sw_max = funcForCell(cell).smax_[wpos]; + double s_r = 1.0-funcForCell(cell).sowcr_; + sr = swcr + (s_r-sw_cr)*(swmax-swcr)/(sw_max-sw_cr); + } + if (s[wpos] <= swcr) { + kr[wpos] = 0.0; + } else if (sr > swmax-1.0e-6) { + if (do_3pt_) { //Ignore krw and do two-point? + kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_; + } else if (!eps_.kro_.empty()){ //Ignore krwr and do two-point + kr[wpos] *= eps_.krw_[cell]/funcForCell(cell).krwmax_; + } + } else if (s[wpos] <= sr) { + kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_; + } else if (s[wpos] <= swmax) { + double krw_max = funcForCell(cell).krwmax_; + double krw = eps_.krw_.empty() ? krw_max : eps_.krw_[cell]; + double krw_r = funcForCell(cell).krwr_; + double krwr = eps_.krwr_.empty() ? krw_r : eps_.krwr_[cell]; + if (std::fabs(krw_max- krw_r) > 1.0e-6) { + kr[wpos] = krwr + (kr[wpos]-krw_r)*(krw-krwr)/(krw_max-krw_r); + } else { + kr[wpos] = krwr + (krw-krwr)*(s[wpos]-sr)/(swmax-sr); + } + } else { + kr[wpos] = eps_.krw_.empty() ? funcForCell(cell).krwmax_ : eps_.krw_[cell]; + } + } + + // - Oil + if (eps_.kro_.empty() && eps_.krorw_.empty()) { // No value scaling + } else if (eps_.krorw_.empty()) { // Two-point scaling + kr[opos] *= (eps_.kro_[cell]/funcForCell(cell).kromax_); + } else { + double sowcr = eps_.sowcr_.empty() ? funcForCell(cell).sowcr_ : eps_.sowcr_[cell]; + double sowmax = eps_.swl_.empty() ? funcForCell(cell).smax_[opos] : 1.0-eps_.swl_[cell]; + double sr; + if (do_3pt_) { + sr = eps_.swcr_.empty() ? 1.0-funcForCell(cell).swcr_ : 1.0-eps_.swcr_[cell]; + } else { + double sow_cr = funcForCell(cell).sowcr_; + double sow_max = funcForCell(cell).smax_[opos]; + double s_r = 1.0-funcForCell(cell).swcr_; + sr = sowcr + (s_r-sow_cr)*(sowmax-sowcr)/(sow_max-sow_cr); + } + if (s[opos] <= sowcr) { + kr[opos] = 0.0; + } else if (sr > sowmax-1.0e-6) { + if (do_3pt_) { //Ignore kro and do two-point? + kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_; + } else if (!eps_.kro_.empty()){ //Ignore krowr and do two-point + kr[opos] *= eps_.kro_[cell]/funcForCell(cell).kromax_; + } + } else if (s[opos] <= sr) { + kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_; + } else if (s[opos] <= sowmax) { + double kro_max = funcForCell(cell).kromax_; + double kro = eps_.kro_.empty() ? kro_max : eps_.kro_[cell]; + double kro_rw = funcForCell(cell).krorw_; + double krorw = eps_.krorw_[cell]; + if (std::fabs(kro_max- kro_rw) > 1.0e-6) { + kr[opos] = krorw + (kr[opos]- kro_rw)*(kro-krorw)/(kro_max- kro_rw); + } else { + kr[opos] = krorw + (kro-krorw)*(s[opos]-sr)/(sowmax-sr); + } + } else { + kr[opos] = eps_.kro_.empty() ? funcForCell(cell).kromax_ : eps_.kro_[cell]; + } + } + } } // namespace Opm From d5edb42b31b51b130485ac1d3e94409055d2e03c Mon Sep 17 00:00:00 2001 From: osae Date: Fri, 30 Nov 2012 12:08:40 +0100 Subject: [PATCH 2/5] Removed tabs. --- opm/core/eclipse/SpecialEclipseFields.hpp | 1963 ++++++++++----------- 1 file changed, 979 insertions(+), 984 deletions(-) diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index 70bff0a0..655b92d4 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -72,9 +72,9 @@ struct SPECGRID : public SpecialBase SPECGRID() { - dimensions.resize(3,1); - numres = 1; - qrdial = 'F'; + dimensions.resize(3,1); + numres = 1; + qrdial = 'F'; } virtual ~SPECGRID() @@ -85,33 +85,33 @@ struct SPECGRID : public SpecialBase virtual void read(std::istream& is) { - const int ndim = 3; - std::vector data(ndim+1,1); - int nread = readDefaultedVectorData(is , data, ndim+1); - int nd = std::min(nread, ndim); - copy(data.begin(), data.begin()+nd, &dimensions[0]); - numres = data[ndim]; - std::string candidate; - is >> candidate; - if (candidate == "/") { - return; - } else { - qrdial = candidate[0]; - } - - if (ignoreSlashLine(is)) { - return; - } else { - THROW("End of file reading" << name()); - } + const int ndim = 3; + std::vector data(ndim+1,1); + int nread = readDefaultedVectorData(is , data, ndim+1); + int nd = std::min(nread, ndim); + copy(data.begin(), data.begin()+nd, &dimensions[0]); + numres = data[ndim]; + std::string candidate; + is >> candidate; + if (candidate == "/") { + return; + } else { + qrdial = candidate[0]; + } + + if (ignoreSlashLine(is)) { + return; + } else { + THROW("End of file reading" << name()); + } } - + virtual void write(std::ostream& os) const { - os << name() << std::endl; - os << dimensions[0] << " " << dimensions[1] << " " - << dimensions[2] << " " << numres << " " << qrdial << std::endl; - os << std::endl; + os << name() << std::endl; + os << dimensions[0] << " " << dimensions[1] << " " + << dimensions[2] << " " << numres << " " << qrdial << std::endl; + os << std::endl; } virtual void convertToSI(const EclipseUnits&) @@ -146,40 +146,40 @@ struct FAULTS : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string fltname; - is >> fltname; - if (fltname[0] == '/') { - is >> ignoreLine; - break; - } - while (fltname.find("--") == 0) { - // This line is a comment - is >> ignoreLine >> fltname; - } - FaultSegment fault_segment; - fault_segment.ijk_coord.resize(6); - fault_segment.fault_name = fltname; - int nread = readDefaultedVectorData(is, fault_segment.ijk_coord, 6); - if (nread != 6) { - THROW("Error reading fault_segment " << fltname); - } - is >> fault_segment.face; - faults.push_back(fault_segment); - ignoreSlashLine(is); - } + while(is) { + std::string fltname; + is >> fltname; + if (fltname[0] == '/') { + is >> ignoreLine; + break; + } + while (fltname.find("--") == 0) { + // This line is a comment + is >> ignoreLine >> fltname; + } + FaultSegment fault_segment; + fault_segment.ijk_coord.resize(6); + fault_segment.fault_name = fltname; + int nread = readDefaultedVectorData(is, fault_segment.ijk_coord, 6); + if (nread != 6) { + THROW("Error reading fault_segment " << fltname); + } + is >> fault_segment.face; + faults.push_back(fault_segment); + ignoreSlashLine(is); + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int)faults.size(); ++i) { - os << faults[i].fault_name << " "; - copy(faults[i].ijk_coord.begin(), faults[i].ijk_coord.end(), - std::ostream_iterator(os, " ")); - os << faults[i].face << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int)faults.size(); ++i) { + os << faults[i].fault_name << " "; + copy(faults[i].ijk_coord.begin(), faults[i].ijk_coord.end(), + std::ostream_iterator(os, " ")); + os << faults[i].face << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits&) @@ -196,7 +196,7 @@ struct MultfltLine double transmis_multiplier; // Transmissibility multiplier double diffusivity_multiplier; // Diffusivity multiplier; MultfltLine() : - fault_name(""), transmis_multiplier(1.0), diffusivity_multiplier(1.0) + fault_name(""), transmis_multiplier(1.0), diffusivity_multiplier(1.0) {} }; @@ -216,38 +216,38 @@ struct MULTFLT : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string fltname; - is >> fltname; - if (fltname[0] == '/') { - is >> ignoreLine; - break; - } - while (fltname == "--") { - // This line is a comment - is >> ignoreLine >> fltname; - } - MultfltLine multflt_line; - multflt_line.fault_name = fltname; - std::vector data(2,1.0); - if (readDefaultedVectorData(is, data, 2) == 2) { - ignoreSlashLine(is); - } - multflt_line.transmis_multiplier = data[0]; - multflt_line.diffusivity_multiplier = data[1]; - multflts.push_back(multflt_line); - } + while(is) { + std::string fltname; + is >> fltname; + if (fltname[0] == '/') { + is >> ignoreLine; + break; + } + while (fltname == "--") { + // This line is a comment + is >> ignoreLine >> fltname; + } + MultfltLine multflt_line; + multflt_line.fault_name = fltname; + std::vector data(2,1.0); + if (readDefaultedVectorData(is, data, 2) == 2) { + ignoreSlashLine(is); + } + multflt_line.transmis_multiplier = data[0]; + multflt_line.diffusivity_multiplier = data[1]; + multflts.push_back(multflt_line); + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int)multflts.size(); ++i) { - os << multflts[i].fault_name << " " - << multflts[i].transmis_multiplier << " " - << multflts[i].diffusivity_multiplier << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int)multflts.size(); ++i) { + os << multflts[i].fault_name << " " + << multflts[i].transmis_multiplier << " " + << multflts[i].diffusivity_multiplier << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits&) @@ -329,20 +329,20 @@ struct DATES : public SpecialBase { return std::string("DATES"); } virtual void read(std::istream& is) { - while(is) { - dates.push_back(readDate(is)); - is >> ignoreWhitespace; - if (is.peek() == int('/')) { - is >> ignoreLine; - break; - } - } + while(is) { + dates.push_back(readDate(is)); + is >> ignoreWhitespace; + if (is.peek() == int('/')) { + is >> ignoreLine; + break; + } + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - copy(dates.begin(), dates.end(), - std::ostream_iterator(os, "\n")); + os << name() << '\n'; + copy(dates.begin(), dates.end(), + std::ostream_iterator(os, "\n")); } virtual void convertToSI(const EclipseUnits&) {} @@ -357,40 +357,40 @@ struct DENSITY : public SpecialBase virtual void read(std::istream& is) { - while (!is.eof()) { - std::vector density(3,-1e100); - if (readDefaultedVectorData(is, density, 3) == 3) { - ignoreSlashLine(is); - } - densities_.push_back(density); + while (!is.eof()) { + std::vector density(3,-1e100); + if (readDefaultedVectorData(is, density, 3) == 3) { + ignoreSlashLine(is); + } + densities_.push_back(density); - int action = next_action(is); // 0:continue 1:return 2:throw - if (action == 1) { - return; // Alphabetic char. Read next keyword. - } else if (action == 2) { - THROW("Error reading DENSITY. Next character is " - << (char)is.peek()); - } - } + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading DENSITY. Next character is " + << (char)is.peek()); + } + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)densities_.size(); ++i) { - os << densities_[i][0] << " " << densities_[i][1] << " " - << densities_[i][2] << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int i=0; i<(int)densities_.size(); ++i) { + os << densities_[i][0] << " " << densities_[i][1] << " " + << densities_[i][2] << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)densities_.size(); ++i) { - densities_[i][0] *= units.density; - densities_[i][1] *= units.density; - densities_[i][2] *= units.density; - } + for (int i=0; i<(int)densities_.size(); ++i) { + densities_[i][0] *= units.density; + densities_[i][1] *= units.density; + densities_[i][2] *= units.density; + } } }; @@ -402,32 +402,32 @@ struct PVDG : public SpecialBase virtual void read(std::istream& is) { - readPvdTable(is, pvdg_, name(), 3); + readPvdTable(is, pvdg_, name(), 3); } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int rn=0; rn<(int)pvdg_.size(); ++rn) { - for (int i=0; i<(int)pvdg_[rn][0].size(); ++i) { - os << pvdg_[rn][0][i] << " " << pvdg_[rn][1][i] << " " - << pvdg_[rn][2][i] << '\n'; - } - os << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int rn=0; rn<(int)pvdg_.size(); ++rn) { + for (int i=0; i<(int)pvdg_[rn][0].size(); ++i) { + os << pvdg_[rn][0][i] << " " << pvdg_[rn][1][i] << " " + << pvdg_[rn][2][i] << '\n'; + } + os << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - double volfac = units.gasvol_r/units.gasvol_s; - for (int rn=0; rn<(int)pvdg_.size(); ++rn) { - for (int i=0; i<(int)pvdg_[rn][0].size(); ++i) { - pvdg_[rn][0][i] *= units.pressure; - pvdg_[rn][1][i] *= volfac; - pvdg_[rn][2][i] *= units.viscosity; - } - } + double volfac = units.gasvol_r/units.gasvol_s; + for (int rn=0; rn<(int)pvdg_.size(); ++rn) { + for (int i=0; i<(int)pvdg_[rn][0].size(); ++i) { + pvdg_[rn][0][i] *= units.pressure; + pvdg_[rn][1][i] *= volfac; + pvdg_[rn][2][i] *= units.viscosity; + } + } } }; @@ -439,32 +439,32 @@ struct PVDO : public SpecialBase virtual void read(std::istream& is) { - readPvdTable(is, pvdo_, name(), 3); + readPvdTable(is, pvdo_, name(), 3); } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int rn=0; rn<(int)pvdo_.size(); ++rn) { - for (int i=0; i<(int)pvdo_[rn][0].size(); ++i) { - os << pvdo_[rn][0][i] << " " << pvdo_[rn][1][i] << " " - << pvdo_[rn][2][i] << '\n'; - } - os << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int rn=0; rn<(int)pvdo_.size(); ++rn) { + for (int i=0; i<(int)pvdo_[rn][0].size(); ++i) { + os << pvdo_[rn][0][i] << " " << pvdo_[rn][1][i] << " " + << pvdo_[rn][2][i] << '\n'; + } + os << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - double volfac = units.liqvol_r/units.liqvol_s; - for (int rn=0; rn<(int)pvdo_.size(); ++rn) { - for (int i=0; i<(int)pvdo_[rn][0].size(); ++i) { - pvdo_[rn][0][i] *= units.pressure; - pvdo_[rn][1][i] *= volfac; - pvdo_[rn][2][i] *= units.viscosity; - } - } + double volfac = units.liqvol_r/units.liqvol_s; + for (int rn=0; rn<(int)pvdo_.size(); ++rn) { + for (int i=0; i<(int)pvdo_[rn][0].size(); ++i) { + pvdo_[rn][0][i] *= units.pressure; + pvdo_[rn][1][i] *= volfac; + pvdo_[rn][2][i] *= units.viscosity; + } + } } }; @@ -476,46 +476,46 @@ struct PVTG : public SpecialBase virtual void read(std::istream& is) { - readPvtTable(is, pvtg_, name()); + readPvtTable(is, pvtg_, name()); } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int rn=0; rn<(int)pvtg_.size(); ++rn) { - for (int i=0; i<(int)pvtg_[rn].size(); ++i) { - int nl = (pvtg_[rn][i].size()-1) / 3; - os << pvtg_[rn][i][0] << " " << pvtg_[rn][i][1] << " " - << pvtg_[rn][i][2] << " " << pvtg_[rn][i][3] << '\n'; - for (int j=1, n=3; j pvtw; - readVectorData(is, pvtw); - if (pvtw.size() == 4) { - pvtw.push_back(0.0); // Not used by frontsim - } - pvtw_.push_back(pvtw); + while (!is.eof()) { + std::vector pvtw; + readVectorData(is, pvtw); + if (pvtw.size() == 4) { + pvtw.push_back(0.0); // Not used by frontsim + } + pvtw_.push_back(pvtw); - int action = next_action(is); // 0:continue 1:return 2:throw - if (action == 1) { - return; // Alphabetic char. Read next keyword. - } else if (action == 2) { - THROW("Error reading PVTW. Next character is " - << (char)is.peek()); - } - } + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading PVTW. Next character is " + << (char)is.peek()); + } + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)pvtw_.size(); ++i) { - os << pvtw_[i][0] << " " << pvtw_[i][1] << " " << pvtw_[i][2] - << " " << pvtw_[i][3] << " " << pvtw_[i][4] << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int i=0; i<(int)pvtw_.size(); ++i) { + os << pvtw_[i][0] << " " << pvtw_[i][1] << " " << pvtw_[i][2] + << " " << pvtw_[i][3] << " " << pvtw_[i][4] << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - double volfac = units.liqvol_r/units.liqvol_s; - for (int i=0; i<(int)pvtw_.size(); ++i) { - pvtw_[i][0] *= units.pressure; - pvtw_[i][1] *= volfac; - pvtw_[i][2] *= units.compressibility; - pvtw_[i][3] *= units.viscosity; - pvtw_[i][4] *= units.compressibility; - } + double volfac = units.liqvol_r/units.liqvol_s; + for (int i=0; i<(int)pvtw_.size(); ++i) { + pvtw_[i][0] *= units.pressure; + pvtw_[i][1] *= volfac; + pvtw_[i][2] *= units.compressibility; + pvtw_[i][3] *= units.viscosity; + pvtw_[i][4] *= units.compressibility; + } } }; @@ -630,37 +630,37 @@ struct ROCK : public SpecialBase virtual void read(std::istream& is) { - while (!is.eof()) { - std::vector rock; - readVectorData(is, rock); - rock_compressibilities_.push_back(rock); + while (!is.eof()) { + std::vector rock; + readVectorData(is, rock); + rock_compressibilities_.push_back(rock); - int action = next_action(is); // 0:continue 1:return 2:throw - if (action == 1) { - return; // Alphabetic char. Read next keyword. - } else if (action == 2) { - THROW("Error reading ROCK. Next character is " - << (char)is.peek()); - } - } + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading ROCK. Next character is " + << (char)is.peek()); + } + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)rock_compressibilities_.size(); ++i) { - os << rock_compressibilities_[i][0] << " " - << rock_compressibilities_[i][1] << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int i=0; i<(int)rock_compressibilities_.size(); ++i) { + os << rock_compressibilities_[i][0] << " " + << rock_compressibilities_[i][1] << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)rock_compressibilities_.size(); ++i) { - rock_compressibilities_[i][0] *= units.pressure; - rock_compressibilities_[i][1] *= units.compressibility; - } + for (int i=0; i<(int)rock_compressibilities_.size(); ++i) { + rock_compressibilities_[i][0] *= units.pressure; + rock_compressibilities_[i][1] *= units.compressibility; + } } }; @@ -673,29 +673,29 @@ struct ROCKTAB : public SpecialBase virtual void read(std::istream& is) { - readPvdTable(is, rocktab_, name(), 3); + readPvdTable(is, rocktab_, name(), 3); } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int rn=0; rn<(int)rocktab_.size(); ++rn) { - for (int i=0; i<(int)rocktab_[rn][0].size(); ++i) { - os << rocktab_[rn][0][i] << " " << rocktab_[rn][1][i] << " " - << rocktab_[rn][2][i] << '\n'; - } - os << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int rn=0; rn<(int)rocktab_.size(); ++rn) { + for (int i=0; i<(int)rocktab_[rn][0].size(); ++i) { + os << rocktab_[rn][0][i] << " " << rocktab_[rn][1][i] << " " + << rocktab_[rn][2][i] << '\n'; + } + os << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int rn=0; rn<(int)rocktab_.size(); ++rn) { - for (int i=0; i<(int)rocktab_[rn][0].size(); ++i) { - rocktab_[rn][0][i] *= units.pressure; - } - } + for (int rn=0; rn<(int)rocktab_.size(); ++rn) { + for (int i=0; i<(int)rocktab_[rn][0].size(); ++i) { + rocktab_[rn][0][i] *= units.pressure; + } + } } }; @@ -710,24 +710,24 @@ struct SGOF : public SpecialBase virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int rn=0; rn<(int)sgof_.size(); ++rn) { - for (int i=0; i<(int)sgof_[rn][0].size(); ++i) { - os << sgof_[rn][0][i] << " " << sgof_[rn][1][i] << " " - << sgof_[rn][2][i] << " " << sgof_[rn][3][i] << '\n'; - } - os << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int rn=0; rn<(int)sgof_.size(); ++rn) { + for (int i=0; i<(int)sgof_[rn][0].size(); ++i) { + os << sgof_[rn][0][i] << " " << sgof_[rn][1][i] << " " + << sgof_[rn][2][i] << " " << sgof_[rn][3][i] << '\n'; + } + os << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int rn=0; rn<(int)sgof_.size(); ++rn) { - for (int i=0; i<(int)sgof_[rn][0].size(); ++i) { - sgof_[rn][3][i] *= units.pressure; - } - } + for (int rn=0; rn<(int)sgof_.size(); ++rn) { + for (int i=0; i<(int)sgof_[rn][0].size(); ++i) { + sgof_[rn][3][i] *= units.pressure; + } + } } }; @@ -741,24 +741,24 @@ struct SWOF : public SpecialBase virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int rn=0; rn<(int)swof_.size(); ++rn) { - for (int i=0; i<(int)swof_[rn][0].size(); ++i) { - os << swof_[rn][0][i] << " " << swof_[rn][1][i] << " " - << swof_[rn][2][i] << " " << swof_[rn][3][i] << '\n'; - } - os << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int rn=0; rn<(int)swof_.size(); ++rn) { + for (int i=0; i<(int)swof_[rn][0].size(); ++i) { + os << swof_[rn][0][i] << " " << swof_[rn][1][i] << " " + << swof_[rn][2][i] << " " << swof_[rn][3][i] << '\n'; + } + os << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int rn=0; rn<(int)swof_.size(); ++rn) { - for (int i=0; i<(int)swof_[rn][0].size(); ++i) { - swof_[rn][3][i] *= units.pressure; - } - } + for (int rn=0; rn<(int)swof_.size(); ++rn) { + for (int i=0; i<(int)swof_[rn][0].size(); ++i) { + swof_[rn][3][i] *= units.pressure; + } + } } }; @@ -804,78 +804,78 @@ struct WELSPECS : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wellname = readString(is); - if (wellname[0] == '/') { - is >> ignoreLine; - break; - } - while (wellname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wellname = readString(is); - } - WelspecsLine welspecs_line; - welspecs_line.name_ = wellname; - welspecs_line.group_ = readString(is); - std::vector int_data(2,1); - readDefaultedVectorData(is, int_data, 2); - welspecs_line.I_ = int_data[0]; - welspecs_line.J_ = int_data[1]; - std::vector double_data(1,-1.0); - readDefaultedVectorData(is, double_data, 1); - welspecs_line.datum_depth_BHP_ = double_data[0]; - welspecs_line.pref_phase_ = readString(is); + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WelspecsLine welspecs_line; + welspecs_line.name_ = wellname; + welspecs_line.group_ = readString(is); + std::vector int_data(2,1); + readDefaultedVectorData(is, int_data, 2); + welspecs_line.I_ = int_data[0]; + welspecs_line.J_ = int_data[1]; + std::vector double_data(1,-1.0); + readDefaultedVectorData(is, double_data, 1); + welspecs_line.datum_depth_BHP_ = double_data[0]; + welspecs_line.pref_phase_ = readString(is); - // HACK! Ignore items 7-13. - ignoreSlashLine(is); - welspecs.push_back(welspecs_line); + // HACK! Ignore items 7-13. + ignoreSlashLine(is); + welspecs.push_back(welspecs_line); - // double_data[0] = 0.0; - // readDefaultedVectorData(is, double_data, 1); - // welspecs_line.drain_rad_ = double_data[0]; - // welspecs_line.spec_inflow_ = readString(is); - // welspecs_line.shut_in_ = readString(is); - // welspecs_line.crossflow_ = readString(is); - // int_data[0] = 0; - // readDefaultedVectorData(is, int_data, 1); - // welspecs_line.pressure_table_number_ = int_data[0]; - // welspecs_line.density_calc_type_ = readString(is); - // int_data[0] = 0; - // readDefaultedVectorData(is, int_data, 1); - // welspecs_line.fluids_in_place_reg_numb_ = int_data[0]; - // welspecs.push_back(welspecs_line); - } + // double_data[0] = 0.0; + // readDefaultedVectorData(is, double_data, 1); + // welspecs_line.drain_rad_ = double_data[0]; + // welspecs_line.spec_inflow_ = readString(is); + // welspecs_line.shut_in_ = readString(is); + // welspecs_line.crossflow_ = readString(is); + // int_data[0] = 0; + // readDefaultedVectorData(is, int_data, 1); + // welspecs_line.pressure_table_number_ = int_data[0]; + // welspecs_line.density_calc_type_ = readString(is); + // int_data[0] = 0; + // readDefaultedVectorData(is, int_data, 1); + // welspecs_line.fluids_in_place_reg_numb_ = int_data[0]; + // welspecs.push_back(welspecs_line); + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int)welspecs.size(); ++i) { - os << welspecs[i].name_ << " " - << welspecs[i].group_ << " " - << welspecs[i].I_ << " " - << welspecs[i].J_ << " " - << welspecs[i].datum_depth_BHP_ << " " - << welspecs[i].pref_phase_ << " " - << welspecs[i].drain_rad_ << " " - << welspecs[i].spec_inflow_ << " " - << welspecs[i].shut_in_ << " " - << welspecs[i].crossflow_ << " " - << welspecs[i].pressure_table_number_ << " " - << welspecs[i].density_calc_type_ << " " - << welspecs[i].fluids_in_place_reg_numb_ << " " + os << name() << std::endl; + for (int i=0; i<(int)welspecs.size(); ++i) { + os << welspecs[i].name_ << " " + << welspecs[i].group_ << " " + << welspecs[i].I_ << " " + << welspecs[i].J_ << " " + << welspecs[i].datum_depth_BHP_ << " " + << welspecs[i].pref_phase_ << " " + << welspecs[i].drain_rad_ << " " + << welspecs[i].spec_inflow_ << " " + << welspecs[i].shut_in_ << " " + << welspecs[i].crossflow_ << " " + << welspecs[i].pressure_table_number_ << " " + << welspecs[i].density_calc_type_ << " " + << welspecs[i].fluids_in_place_reg_numb_ << " " << "/\n"; - } - os << '/' << std::endl; + } + os << '/' << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)welspecs.size(); ++i) { - welspecs[i].datum_depth_BHP_ *= units.length; - welspecs[i].drain_rad_ *= units.length; - } + for (int i=0; i<(int)welspecs.size(); ++i) { + welspecs[i].datum_depth_BHP_ *= units.length; + welspecs[i].drain_rad_ *= units.length; + } } }; @@ -896,17 +896,17 @@ struct CompdatLine // Default values CompdatLine() : - open_shut_flag_("OPEN"), + open_shut_flag_("OPEN"), sat_table_number_(0), connect_trans_fac_(0.0), diameter_(0.0), Kh_(-1.0), skin_factor_(0.0), D_factor_(-1e100), - penetration_direct_("Z"), + penetration_direct_("Z"), r0_(-1.0) { - grid_ind_.resize(4); + grid_ind_.resize(4); } }; @@ -926,69 +926,69 @@ struct COMPDAT : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wellname = readString(is); - if (wellname[0] == '/') { - is >> ignoreLine; - break; - } - while (wellname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wellname = readString(is); - } - CompdatLine compdat_line; - compdat_line.well_ = wellname; - readDefaultedVectorData(is, compdat_line.grid_ind_, 4); - compdat_line.open_shut_flag_ = readString(is); - std::vector int_data(1,-1); - readDefaultedVectorData(is, int_data, 1); - compdat_line.sat_table_number_ = int_data[0]; + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + CompdatLine compdat_line; + compdat_line.well_ = wellname; + readDefaultedVectorData(is, compdat_line.grid_ind_, 4); + compdat_line.open_shut_flag_ = readString(is); + std::vector int_data(1,-1); + readDefaultedVectorData(is, int_data, 1); + compdat_line.sat_table_number_ = int_data[0]; std::vector double_data(2, 0.0); int num_to_read = 2; int num_read = readDefaultedVectorData(is, double_data, num_to_read); compdat_line.connect_trans_fac_ = double_data[0]; compdat_line.diameter_ = double_data[1]; - // HACK! Ignore items 10-14. + // HACK! Ignore items 10-14. if (num_read == num_to_read) { ignoreSlashLine(is); } - compdat.push_back(compdat_line); - } + compdat.push_back(compdat_line); + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int)compdat.size(); ++i) { - os << compdat[i].well_ << " " - << compdat[i].grid_ind_[0] << " " - << compdat[i].grid_ind_[1] << " " - << compdat[i].grid_ind_[2] << " " - << compdat[i].grid_ind_[3] << " " - << compdat[i].open_shut_flag_ << " " - << compdat[i].sat_table_number_ << " " - << compdat[i].connect_trans_fac_ << " " - << compdat[i].diameter_ << " " - << compdat[i].Kh_ << " " - << compdat[i].skin_factor_ << " " - << compdat[i].D_factor_ << " " - << compdat[i].penetration_direct_ << " " - << compdat[i].r0_ << " " + os << name() << std::endl; + for (int i=0; i<(int)compdat.size(); ++i) { + os << compdat[i].well_ << " " + << compdat[i].grid_ind_[0] << " " + << compdat[i].grid_ind_[1] << " " + << compdat[i].grid_ind_[2] << " " + << compdat[i].grid_ind_[3] << " " + << compdat[i].open_shut_flag_ << " " + << compdat[i].sat_table_number_ << " " + << compdat[i].connect_trans_fac_ << " " + << compdat[i].diameter_ << " " + << compdat[i].Kh_ << " " + << compdat[i].skin_factor_ << " " + << compdat[i].D_factor_ << " " + << compdat[i].penetration_direct_ << " " + << compdat[i].r0_ << " " << "/\n"; - } - os << '/' << std::endl; + } + os << '/' << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)compdat.size(); ++i) { - compdat[i].connect_trans_fac_ *= units.transmissibility; - compdat[i].diameter_ *= units.length; - compdat[i].Kh_ *= units.permeability*units.length; - compdat[i].r0_ *= units.length; - } + for (int i=0; i<(int)compdat.size(); ++i) { + compdat[i].connect_trans_fac_ *= units.transmissibility; + compdat[i].diameter_ *= units.length; + compdat[i].Kh_ *= units.permeability*units.length; + compdat[i].r0_ *= units.length; + } } }; @@ -1007,7 +1007,7 @@ struct GconinjeLine // Default values GconinjeLine() : - surface_flow_max_rate_(-1.0E20), + surface_flow_max_rate_(-1.0E20), resv_flow_max_rate_(-1E20), reinjection_fraction_target_(1), voidage_replacement_fraction_(1) @@ -1031,61 +1031,61 @@ struct GCONINJE : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string groupname = readString(is); - if (groupname[0] == '/') { - is >> ignoreLine; - break; - } - while (groupname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - groupname = readString(is); - } - GconinjeLine gconinje_line; - gconinje_line.group_ = groupname; - gconinje_line.injector_type_ = readString(is); - gconinje_line.control_mode_ = readString(is); - std::vector double_data(10, -1.0E20); + while(is) { + std::string groupname = readString(is); + if (groupname[0] == '/') { + is >> ignoreLine; + break; + } + while (groupname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + groupname = readString(is); + } + GconinjeLine gconinje_line; + gconinje_line.group_ = groupname; + gconinje_line.injector_type_ = readString(is); + gconinje_line.control_mode_ = readString(is); + std::vector double_data(10, -1.0E20); const int num_to_read = 4; - int num_read = readDefaultedVectorData(is, double_data, num_to_read); - gconinje_line.surface_flow_max_rate_ = double_data[0]; + int num_read = readDefaultedVectorData(is, double_data, num_to_read); + gconinje_line.surface_flow_max_rate_ = double_data[0]; gconinje_line.resv_flow_max_rate_ = double_data[1]; gconinje_line.reinjection_fraction_target_ = double_data[2]; gconinje_line.voidage_replacement_fraction_ = double_data[3]; - // HACK! Ignore any further items + // HACK! Ignore any further items if (num_read == num_to_read) { ignoreSlashLine(is); } - gconinje.push_back(gconinje_line); + gconinje.push_back(gconinje_line); - } + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int) gconinje.size(); ++i) { - os << gconinje[i].group_ << " " - << gconinje[i].injector_type_ << " " - << gconinje[i].control_mode_ << " " - << gconinje[i].surface_flow_max_rate_ << " " - << gconinje[i].reinjection_fraction_target_ - << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int) gconinje.size(); ++i) { + os << gconinje[i].group_ << " " + << gconinje[i].injector_type_ << " " + << gconinje[i].control_mode_ << " " + << gconinje[i].surface_flow_max_rate_ << " " + << gconinje[i].reinjection_fraction_target_ + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int) gconinje.size(); ++i) { + for (int i=0; i<(int) gconinje.size(); ++i) { if (gconinje[i].injector_type_ == "GAS") { gconinje[i].surface_flow_max_rate_ *= units.gasvol_s/units.time; } else { gconinje[i].surface_flow_max_rate_ *= units.liqvol_s/units.time; } - } + } } }; @@ -1111,9 +1111,9 @@ struct WconinjeLine // Default values WconinjeLine() : - open_shut_flag_("OPEN"), surface_flow_max_rate_(-1.0E20), - reservoir_flow_max_rate_(-1.0E20), BHP_limit_(6891), THP_limit_(-1.0E20), - VFP_table_number_(0), concentration_(0.0) + open_shut_flag_("OPEN"), surface_flow_max_rate_(-1.0E20), + reservoir_flow_max_rate_(-1.0E20), BHP_limit_(6891), THP_limit_(-1.0E20), + VFP_table_number_(0), concentration_(0.0) { } }; @@ -1134,74 +1134,74 @@ struct WCONINJE : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wellname = readString(is); - if (wellname[0] == '/') { - is >> ignoreLine; - break; - } - while (wellname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wellname = readString(is); - } - WconinjeLine wconinje_line; - wconinje_line.well_ = wellname; - wconinje_line.injector_type_ = readString(is); - wconinje_line.open_shut_flag_ = readString(is); - wconinje_line.control_mode_ = readString(is); - std::vector double_data(6, -1.0E20); - double_data[2] = wconinje_line.BHP_limit_; - double_data[4] = wconinje_line.VFP_table_number_; - double_data[5] = wconinje_line.concentration_; + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WconinjeLine wconinje_line; + wconinje_line.well_ = wellname; + wconinje_line.injector_type_ = readString(is); + wconinje_line.open_shut_flag_ = readString(is); + wconinje_line.control_mode_ = readString(is); + std::vector double_data(6, -1.0E20); + double_data[2] = wconinje_line.BHP_limit_; + double_data[4] = wconinje_line.VFP_table_number_; + double_data[5] = wconinje_line.concentration_; const int num_to_read = 6; - int num_read = readDefaultedVectorData(is, double_data, num_to_read); - wconinje_line.surface_flow_max_rate_ = double_data[0]; - wconinje_line.reservoir_flow_max_rate_ = double_data[1]; - wconinje_line.BHP_limit_ = double_data[2]; - wconinje_line.THP_limit_ = double_data[3]; - wconinje_line.VFP_table_number_ = (int)double_data[4]; - wconinje_line.concentration_ = double_data[5]; - // HACK! Ignore any further items + int num_read = readDefaultedVectorData(is, double_data, num_to_read); + wconinje_line.surface_flow_max_rate_ = double_data[0]; + wconinje_line.reservoir_flow_max_rate_ = double_data[1]; + wconinje_line.BHP_limit_ = double_data[2]; + wconinje_line.THP_limit_ = double_data[3]; + wconinje_line.VFP_table_number_ = (int)double_data[4]; + wconinje_line.concentration_ = double_data[5]; + // HACK! Ignore any further items if (num_read == num_to_read) { ignoreSlashLine(is); } - wconinje.push_back(wconinje_line); - } + wconinje.push_back(wconinje_line); + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int) wconinje.size(); ++i) { - os << wconinje[i].well_ << " " - << wconinje[i].injector_type_ << " " - << wconinje[i].open_shut_flag_ << " " - << wconinje[i].control_mode_ << " " - << wconinje[i].surface_flow_max_rate_ << " " - << wconinje[i].reservoir_flow_max_rate_ << " " - << wconinje[i].BHP_limit_ << " " - << wconinje[i].THP_limit_ << " " - << wconinje[i].VFP_table_number_ << " " - << wconinje[i].concentration_ - << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int) wconinje.size(); ++i) { + os << wconinje[i].well_ << " " + << wconinje[i].injector_type_ << " " + << wconinje[i].open_shut_flag_ << " " + << wconinje[i].control_mode_ << " " + << wconinje[i].surface_flow_max_rate_ << " " + << wconinje[i].reservoir_flow_max_rate_ << " " + << wconinje[i].BHP_limit_ << " " + << wconinje[i].THP_limit_ << " " + << wconinje[i].VFP_table_number_ << " " + << wconinje[i].concentration_ + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int) wconinje.size(); ++i) { + for (int i=0; i<(int) wconinje.size(); ++i) { if (wconinje[i].injector_type_ == "GAS") { wconinje[i].surface_flow_max_rate_ *= units.gasvol_s/units.time; } else { wconinje[i].surface_flow_max_rate_ *= units.liqvol_s/units.time; } - wconinje[i].reservoir_flow_max_rate_ *= units.liqvol_r/units.time; - wconinje[i].BHP_limit_ *= units.pressure; - wconinje[i].THP_limit_ *= units.pressure; - wconinje[i].concentration_ *= units.gasvol_s/units.liqvol_s; // ??? @bsp 10 - } + wconinje[i].reservoir_flow_max_rate_ *= units.liqvol_r/units.time; + wconinje[i].BHP_limit_ *= units.pressure; + wconinje[i].THP_limit_ *= units.pressure; + wconinje[i].concentration_ *= units.gasvol_s/units.liqvol_s; // ??? @bsp 10 + } } }; @@ -1220,8 +1220,8 @@ struct GconprodLine std::string procedure_; // Procedure on exceeding a maximum rate limit // Default values GconprodLine() : - oil_max_rate_(-1.0E20), water_max_rate_(-1.0E20), - gas_max_rate_(-1.0E20), liquid_max_rate_(-1.0E20), + oil_max_rate_(-1.0E20), water_max_rate_(-1.0E20), + gas_max_rate_(-1.0E20), liquid_max_rate_(-1.0E20), resv_max_rate_(-1.0E20) { } @@ -1243,25 +1243,25 @@ struct GCONPROD : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string groupname = readString(is); - if (groupname[0] == '/') { - is >> ignoreLine; - break; - } - while (groupname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - groupname = readString(is); - } - GconprodLine gconprod_line; - gconprod_line.group_ = groupname; - gconprod_line.control_mode_ = readString(is); + while(is) { + std::string groupname = readString(is); + if (groupname[0] == '/') { + is >> ignoreLine; + break; + } + while (groupname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + groupname = readString(is); + } + GconprodLine gconprod_line; + gconprod_line.group_ = groupname; + gconprod_line.control_mode_ = readString(is); if (gconprod_line.control_mode_[gconprod_line.control_mode_.size() - 1] == '*') { gconprod_line.control_mode_ = "NONE"; } - std::vector double_data(4, -1.0E20); + std::vector double_data(4, -1.0E20); const int num_to_read = 4; int num_read = readDefaultedVectorData(is, double_data, num_to_read); gconprod_line.oil_max_rate_ = double_data[0]; @@ -1319,39 +1319,39 @@ struct GCONPROD : public SpecialBase } - } + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int) gconprod.size(); ++i) { - os << gconprod[i].group_ << " " - << gconprod[i].control_mode_ << " " - << gconprod[i].oil_max_rate_ << " " - << gconprod[i].water_max_rate_ << " " - << gconprod[i].gas_max_rate_ << " " - << gconprod[i].liquid_max_rate_ << " " - << gconprod[i].procedure_ << " " + os << name() << std::endl; + for (int i=0; i<(int) gconprod.size(); ++i) { + os << gconprod[i].group_ << " " + << gconprod[i].control_mode_ << " " + << gconprod[i].oil_max_rate_ << " " + << gconprod[i].water_max_rate_ << " " + << gconprod[i].gas_max_rate_ << " " + << gconprod[i].liquid_max_rate_ << " " + << gconprod[i].procedure_ << " " << gconprod[i].resv_max_rate_ - << std::endl; - } - os << std::endl; + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - double lrat = units.liqvol_s / units.time; - double grat = units.gasvol_s / units.time; - double resv = units.liqvol_r / units.time; + double lrat = units.liqvol_s / units.time; + double grat = units.gasvol_s / units.time; + double resv = units.liqvol_r / units.time; - for (int i=0; i<(int) gconprod.size(); ++i) { - gconprod[i].oil_max_rate_ *= lrat; - gconprod[i].water_max_rate_ *= lrat; - gconprod[i].gas_max_rate_ *= grat; - gconprod[i].liquid_max_rate_ *= lrat; + for (int i=0; i<(int) gconprod.size(); ++i) { + gconprod[i].oil_max_rate_ *= lrat; + gconprod[i].water_max_rate_ *= lrat; + gconprod[i].gas_max_rate_ *= grat; + gconprod[i].liquid_max_rate_ *= lrat; gconprod[i].resv_max_rate_ *= resv; - } + } } }; @@ -1383,63 +1383,63 @@ struct WGRUPCON : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wname = readString(is); - if (wname[0] == '/') { - is >> ignoreLine; - break; - } - while (wname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wname = readString(is); - } - WgrupconLine wgrupcon_line; - wgrupcon_line.well_ = wname; + while(is) { + std::string wname = readString(is); + if (wname[0] == '/') { + is >> ignoreLine; + break; + } + while (wname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wname = readString(is); + } + WgrupconLine wgrupcon_line; + wgrupcon_line.well_ = wname; std::string available = readString(is); if (available[available.size()-1] == '*') { available = "YES"; } wgrupcon_line.available_for_group_control_ = (available == "YES"); - std::vector double_data(1, 1.0E20); + std::vector double_data(1, 1.0E20); const int num_to_read = 1; - int num_read = readDefaultedVectorData(is, double_data, num_to_read); + int num_read = readDefaultedVectorData(is, double_data, num_to_read); wgrupcon_line.guide_rate_ = double_data[0]; std::string phase = readString(is); - if (phase[0] == '/') { - is >> ignoreLine; - break; - } - while (phase.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - phase = readString(is); - } + if (phase[0] == '/') { + is >> ignoreLine; + break; + } + while (phase.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + phase = readString(is); + } wgrupcon_line.phase_ = phase; - - wgrupcon.push_back(wgrupcon_line); - // HACK! Ignore any further items + + wgrupcon.push_back(wgrupcon_line); + // HACK! Ignore any further items if (num_read == num_to_read) { ignoreSlashLine(is); } - } + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int) wgrupcon.size(); ++i) { - os << wgrupcon[i].well_ << " " - << wgrupcon[i].available_for_group_control_ << " " - << wgrupcon[i].guide_rate_ << " " - << wgrupcon[i].phase_ << " " - << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int) wgrupcon.size(); ++i) { + os << wgrupcon[i].well_ << " " + << wgrupcon[i].available_for_group_control_ << " " + << wgrupcon[i].guide_rate_ << " " + << wgrupcon[i].phase_ << " " + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& units) @@ -1470,10 +1470,10 @@ struct WconprodLine // depending on units used, but it is recommended practice // to not leave the BHP limit defaulted. WconprodLine() : - open_shut_flag_("OPEN"), oil_max_rate_(-1.0E20), water_max_rate_(-1.0E20), - gas_max_rate_(-1.0E20), liquid_max_rate_(-1.0E20), - reservoir_flow_max_rate_(-1.0E20), BHP_limit_(1.0), THP_limit_(0.0), - VFP_table_number_(0), artif_lift_quantity_(0.0) + open_shut_flag_("OPEN"), oil_max_rate_(-1.0E20), water_max_rate_(-1.0E20), + gas_max_rate_(-1.0E20), liquid_max_rate_(-1.0E20), + reservoir_flow_max_rate_(-1.0E20), BHP_limit_(1.0), THP_limit_(0.0), + VFP_table_number_(0), artif_lift_quantity_(0.0) { } }; @@ -1494,81 +1494,81 @@ struct WCONPROD : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wellname = readString(is); - if (wellname[0] == '/') { - is >> ignoreLine; - break; - } - while (wellname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wellname = readString(is); - } - WconprodLine wconprod_line; - wconprod_line.well_ = wellname; - wconprod_line.open_shut_flag_ = readString(is); - wconprod_line.control_mode_ = readString(is); - std::vector double_data(14, -1.0E20); - double_data[5] = wconprod_line.BHP_limit_; - double_data[6] = wconprod_line.THP_limit_; - double_data[7] = wconprod_line.VFP_table_number_; - double_data[8] = wconprod_line.artif_lift_quantity_; + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WconprodLine wconprod_line; + wconprod_line.well_ = wellname; + wconprod_line.open_shut_flag_ = readString(is); + wconprod_line.control_mode_ = readString(is); + std::vector double_data(14, -1.0E20); + double_data[5] = wconprod_line.BHP_limit_; + double_data[6] = wconprod_line.THP_limit_; + double_data[7] = wconprod_line.VFP_table_number_; + double_data[8] = wconprod_line.artif_lift_quantity_; const int num_to_read = 14; - int num_read = readDefaultedVectorData(is, double_data, num_to_read); - wconprod_line.oil_max_rate_ = double_data[0]; - wconprod_line.water_max_rate_ = double_data[1]; - wconprod_line.gas_max_rate_ = double_data[2]; - wconprod_line.liquid_max_rate_ = double_data[3]; - wconprod_line.reservoir_flow_max_rate_ = double_data[4]; - wconprod_line.BHP_limit_ = double_data[5]; - wconprod_line.THP_limit_ = double_data[6]; - wconprod_line.VFP_table_number_ = (int)double_data[7]; - wconprod_line.artif_lift_quantity_ = double_data[8]; - wconprod.push_back(wconprod_line); + int num_read = readDefaultedVectorData(is, double_data, num_to_read); + wconprod_line.oil_max_rate_ = double_data[0]; + wconprod_line.water_max_rate_ = double_data[1]; + wconprod_line.gas_max_rate_ = double_data[2]; + wconprod_line.liquid_max_rate_ = double_data[3]; + wconprod_line.reservoir_flow_max_rate_ = double_data[4]; + wconprod_line.BHP_limit_ = double_data[5]; + wconprod_line.THP_limit_ = double_data[6]; + wconprod_line.VFP_table_number_ = (int)double_data[7]; + wconprod_line.artif_lift_quantity_ = double_data[8]; + wconprod.push_back(wconprod_line); if(num_to_read == num_read) { ignoreSlashLine(is); } - } + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int) wconprod.size(); ++i) { - os << wconprod[i].well_ << " " - << wconprod[i].open_shut_flag_ << " " - << wconprod[i].control_mode_ << " " - << wconprod[i].oil_max_rate_ << " " - << wconprod[i].water_max_rate_ << " " - << wconprod[i].gas_max_rate_ << " " - << wconprod[i].liquid_max_rate_ << " " - << wconprod[i].reservoir_flow_max_rate_ << " " - << wconprod[i].BHP_limit_ << " " - << wconprod[i].THP_limit_ << " " - << wconprod[i].VFP_table_number_ << " " - << wconprod[i].artif_lift_quantity_ - << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int) wconprod.size(); ++i) { + os << wconprod[i].well_ << " " + << wconprod[i].open_shut_flag_ << " " + << wconprod[i].control_mode_ << " " + << wconprod[i].oil_max_rate_ << " " + << wconprod[i].water_max_rate_ << " " + << wconprod[i].gas_max_rate_ << " " + << wconprod[i].liquid_max_rate_ << " " + << wconprod[i].reservoir_flow_max_rate_ << " " + << wconprod[i].BHP_limit_ << " " + << wconprod[i].THP_limit_ << " " + << wconprod[i].VFP_table_number_ << " " + << wconprod[i].artif_lift_quantity_ + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - double lrat = units.liqvol_s / units.time; - double grat = units.gasvol_s / units.time; - double resv = units.liqvol_r / units.time; - for (int i=0; i<(int) wconprod.size(); ++i) { - wconprod[i].oil_max_rate_ *= lrat; - wconprod[i].water_max_rate_ *= lrat; - wconprod[i].gas_max_rate_ *= grat; - wconprod[i].liquid_max_rate_ *= lrat; - wconprod[i].reservoir_flow_max_rate_ *= resv; - wconprod[i].BHP_limit_ *= units.pressure; - wconprod[i].THP_limit_ *= units.pressure; - } + double lrat = units.liqvol_s / units.time; + double grat = units.gasvol_s / units.time; + double resv = units.liqvol_r / units.time; + for (int i=0; i<(int) wconprod.size(); ++i) { + wconprod[i].oil_max_rate_ *= lrat; + wconprod[i].water_max_rate_ *= lrat; + wconprod[i].gas_max_rate_ *= grat; + wconprod[i].liquid_max_rate_ *= lrat; + wconprod[i].reservoir_flow_max_rate_ *= resv; + wconprod[i].BHP_limit_ *= units.pressure; + wconprod[i].THP_limit_ *= units.pressure; + } } }; @@ -1603,60 +1603,60 @@ struct WELTARG : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wellname = readString(is); - if (wellname[0] == '/') { - is >> ignoreLine; - break; - } - while (wellname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wellname = readString(is); - } - WeltargLine weltarg_line; - weltarg_line.well_ = wellname; - weltarg_line.control_change_ = readString(is); - is >> weltarg_line.new_value_; - ignoreSlashLine(is); - weltarg.push_back(weltarg_line); - } + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WeltargLine weltarg_line; + weltarg_line.well_ = wellname; + weltarg_line.control_change_ = readString(is); + is >> weltarg_line.new_value_; + ignoreSlashLine(is); + weltarg.push_back(weltarg_line); + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int)weltarg.size(); ++i) { - os << weltarg[i].well_ << " " - << weltarg[i].control_change_ << " " - << weltarg[i].new_value_ << " " - << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int)weltarg.size(); ++i) { + os << weltarg[i].well_ << " " + << weltarg[i].control_change_ << " " + << weltarg[i].new_value_ << " " + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - double lrat = units.liqvol_s / units.time; - double grat = units.gasvol_s / units.time; - double resv = units.liqvol_r / units.time; - for (int i=0; i<(int) weltarg.size(); ++i) { - if (weltarg[i].control_change_[0] == 'O' || - weltarg[i].control_change_[0] == 'W' || - weltarg[i].control_change_[0] == 'L') { - weltarg[i].new_value_ *= lrat; - } else if (weltarg[i].control_change_[0] == 'G') { - weltarg[i].new_value_ *= grat; - } else if (weltarg[i].control_change_[0] == 'R') { - weltarg[i].new_value_ *= resv; - } else if (weltarg[i].control_change_[0] == 'B' || - weltarg[i].control_change_[0] == 'T') { - weltarg[i].new_value_ *= units.pressure; - } else { - THROW("WELTARG. Unknown control or constraint " - << weltarg[i].control_change_[0]); - } - } + double lrat = units.liqvol_s / units.time; + double grat = units.gasvol_s / units.time; + double resv = units.liqvol_r / units.time; + for (int i=0; i<(int) weltarg.size(); ++i) { + if (weltarg[i].control_change_[0] == 'O' || + weltarg[i].control_change_[0] == 'W' || + weltarg[i].control_change_[0] == 'L') { + weltarg[i].new_value_ *= lrat; + } else if (weltarg[i].control_change_[0] == 'G') { + weltarg[i].new_value_ *= grat; + } else if (weltarg[i].control_change_[0] == 'R') { + weltarg[i].new_value_ *= resv; + } else if (weltarg[i].control_change_[0] == 'B' || + weltarg[i].control_change_[0] == 'T') { + weltarg[i].new_value_ *= units.pressure; + } else { + THROW("WELTARG. Unknown control or constraint " + << weltarg[i].control_change_[0]); + } + } } }; @@ -1684,34 +1684,34 @@ struct WELOPEN : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wellname = readString(is); - if (wellname[0] == '/') { - is >> ignoreLine; - break; - } - while (wellname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wellname = readString(is); - } - WelopenLine welopen_line; - welopen_line.well_ = wellname; - welopen_line.openshutflag_ = readString(is); - ignoreSlashLine(is); - welopen.push_back(welopen_line); - } + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WelopenLine welopen_line; + welopen_line.well_ = wellname; + welopen_line.openshutflag_ = readString(is); + ignoreSlashLine(is); + welopen.push_back(welopen_line); + } } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int)welopen.size(); ++i) { - os << welopen[i].well_ << " " - << welopen[i].openshutflag_ - << std::endl; - } - os << std::endl; + os << name() << std::endl; + for (int i=0; i<(int)welopen.size(); ++i) { + os << welopen[i].well_ << " " + << welopen[i].openshutflag_ + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& /*units*/) @@ -1768,56 +1768,56 @@ struct EQUIL : public SpecialBase virtual void read(std::istream& is) { - // Note. This function assumes that NTEQUL = 1, and reads only one line. - int num_to_read = 9; - std::vector data(num_to_read,0); - int num_read = readDefaultedVectorData(is, data, num_to_read); - if (num_read == num_to_read) { - ignoreSlashLine(is); - } + // Note. This function assumes that NTEQUL = 1, and reads only one line. + int num_to_read = 9; + std::vector data(num_to_read,0); + int num_read = readDefaultedVectorData(is, data, num_to_read); + if (num_read == num_to_read) { + ignoreSlashLine(is); + } - EquilLine equil_line; - equil_line.datum_depth_ = data[0]; - equil_line.datum_depth_pressure_ = data[1]; - equil_line.water_oil_contact_depth_ = data[2]; - equil_line.oil_water_cap_pressure_ = data[3]; - equil_line.gas_oil_contact_depth_ = data[4]; - equil_line.gas_oil_cap_pressure_ = data[5]; - equil_line.live_oil_table_index_ = int(data[6]); - equil_line.wet_gas_table_index_ = int(data[7]); - equil_line.N_ = int(data[8]); - equil.push_back(equil_line); + EquilLine equil_line; + equil_line.datum_depth_ = data[0]; + equil_line.datum_depth_pressure_ = data[1]; + equil_line.water_oil_contact_depth_ = data[2]; + equil_line.oil_water_cap_pressure_ = data[3]; + equil_line.gas_oil_contact_depth_ = data[4]; + equil_line.gas_oil_cap_pressure_ = data[5]; + equil_line.live_oil_table_index_ = int(data[6]); + equil_line.wet_gas_table_index_ = int(data[7]); + equil_line.N_ = int(data[8]); + equil.push_back(equil_line); } virtual void write(std::ostream& os) const { - os << name() << std::endl; - for (int i=0; i<(int)equil.size(); ++i) { - os << equil[i].datum_depth_ << " " - << equil[i].datum_depth_pressure_ << " " - << equil[i].water_oil_contact_depth_ << " " - << equil[i].oil_water_cap_pressure_ << " " - << equil[i].gas_oil_contact_depth_ << " " - << equil[i].gas_oil_cap_pressure_ << " " - << equil[i].live_oil_table_index_ << " " - << equil[i].wet_gas_table_index_ << " " - << equil[i].N_ << " " + os << name() << std::endl; + for (int i=0; i<(int)equil.size(); ++i) { + os << equil[i].datum_depth_ << " " + << equil[i].datum_depth_pressure_ << " " + << equil[i].water_oil_contact_depth_ << " " + << equil[i].oil_water_cap_pressure_ << " " + << equil[i].gas_oil_contact_depth_ << " " + << equil[i].gas_oil_cap_pressure_ << " " + << equil[i].live_oil_table_index_ << " " + << equil[i].wet_gas_table_index_ << " " + << equil[i].N_ << " " - << std::endl; - } - os << std::endl; + << std::endl; + } + os << std::endl; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)equil.size(); ++i) { - equil[i].datum_depth_ *= units.length; - equil[i].datum_depth_pressure_ *= units.pressure; - equil[i].water_oil_contact_depth_ *= units.length; - equil[i].oil_water_cap_pressure_ *= units.pressure; - equil[i].gas_oil_contact_depth_ *= units.length; - equil[i].gas_oil_cap_pressure_ *= units.pressure; - } + for (int i=0; i<(int)equil.size(); ++i) { + equil[i].datum_depth_ *= units.length; + equil[i].datum_depth_pressure_ *= units.pressure; + equil[i].water_oil_contact_depth_ *= units.length; + equil[i].oil_water_cap_pressure_ *= units.pressure; + equil[i].gas_oil_contact_depth_ *= units.length; + equil[i].gas_oil_cap_pressure_ *= units.pressure; + } } }; @@ -1829,44 +1829,44 @@ struct PVCDO : public SpecialBase virtual void read(std::istream& is) { - while (!is.eof()) { - std::vector pvcdo; - readVectorData(is, pvcdo); - if (pvcdo.size() == 4) { - pvcdo.push_back(0.0); - } - pvcdo_.push_back(pvcdo); + while (!is.eof()) { + std::vector pvcdo; + readVectorData(is, pvcdo); + if (pvcdo.size() == 4) { + pvcdo.push_back(0.0); + } + pvcdo_.push_back(pvcdo); - int action = next_action(is); // 0:continue 1:return 2:throw - if (action == 1) { - return; // Alphabetic char. Read next keyword. - } else if (action == 2) { - THROW("Error reading PVCDO. Next character is " - << (char)is.peek()); - } - } + int action = next_action(is); // 0:continue 1:return 2:throw + if (action == 1) { + return; // Alphabetic char. Read next keyword. + } else if (action == 2) { + THROW("Error reading PVCDO. Next character is " + << (char)is.peek()); + } + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)pvcdo_.size(); ++i) { - os << pvcdo_[i][0] << " " << pvcdo_[i][1] << " " << pvcdo_[i][2] - << " " << pvcdo_[i][3] << " " << pvcdo_[i][4] << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int i=0; i<(int)pvcdo_.size(); ++i) { + os << pvcdo_[i][0] << " " << pvcdo_[i][1] << " " << pvcdo_[i][2] + << " " << pvcdo_[i][3] << " " << pvcdo_[i][4] << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - double volfac = units.liqvol_r/units.liqvol_s; - for (int i=0; i<(int)pvcdo_.size(); ++i) { - pvcdo_[i][0] *= units.pressure; - pvcdo_[i][1] *= volfac; - pvcdo_[i][2] *= units.compressibility; - pvcdo_[i][3] *= units.viscosity; - pvcdo_[i][4] *= units.compressibility; - } + double volfac = units.liqvol_r/units.liqvol_s; + for (int i=0; i<(int)pvcdo_.size(); ++i) { + pvcdo_[i][0] *= units.pressure; + pvcdo_[i][1] *= volfac; + pvcdo_[i][2] *= units.compressibility; + pvcdo_[i][3] *= units.viscosity; + pvcdo_[i][4] *= units.compressibility; + } } }; @@ -1878,27 +1878,27 @@ struct TSTEP : public SpecialBase virtual void read(std::istream& is) { - std::vector tstep; - readVectorData(is, tstep); - if (!tstep.empty()) { - tstep_.insert(tstep_.end(), tstep.begin(), tstep.end()); - } + std::vector tstep; + readVectorData(is, tstep); + if (!tstep.empty()) { + tstep_.insert(tstep_.end(), tstep.begin(), tstep.end()); + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - copy(tstep_.begin(), tstep_.end(), - std::ostream_iterator(os, " ")); - os << '\n'; + os << name() << '\n'; + copy(tstep_.begin(), tstep_.end(), + std::ostream_iterator(os, " ")); + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - int num_steps = tstep_.size(); - for (int i = 0; i < num_steps; ++i) { - tstep_[i] *= units.time; - } + int num_steps = tstep_.size(); + for (int i = 0; i < num_steps; ++i) { + tstep_[i] *= units.time; + } } }; @@ -1909,43 +1909,43 @@ struct MultRec : public SpecialBase #ifdef VERBOSE std::cout << "(dummy implementation)" << std::endl; #endif - const std::ctype& ct = std::use_facet< std::ctype >(std::locale::classic()); - is >> ignoreSlashLine; + const std::ctype& ct = std::use_facet< std::ctype >(std::locale::classic()); + is >> ignoreSlashLine; while (!is.eof()) { is >> ignoreWhitespace; - std::streampos pos = is.tellg(); + std::streampos pos = is.tellg(); char c; is.get(c); - if (is.eof()) { - return; - } - if (ct.is(std::ctype_base::alpha, c)) { - std::string name; // Unquoted name or new keyword? - std::getline(is, name); - if (name.rfind('/') != std::string::npos) { - continue; // Unquoted name - } else { - is.seekg(pos); - break; // Read next keyword - } - } else if (ct.is(std::ctype_base::digit, c) || c== '.') { - is >> ignoreSlashLine; // Decimal digit. Ignore data. - continue; - } else if (c== '\'') { - is >> ignoreSlashLine; // Quote. Ignore data. - continue; - } else if(c == '-' && is.peek() == int('-')) { - is >> ignoreLine; // This line is a comment - continue; - } else if (c == '/' ) { - is >> ignoreLine; // This line is a null record. - continue; // (No data before slash) - } else { - is.putback(c); - std::string temp; - is >> temp; - std::cout << "READ ERROR! Next word is " << temp << std::endl; - } + if (is.eof()) { + return; + } + if (ct.is(std::ctype_base::alpha, c)) { + std::string name; // Unquoted name or new keyword? + std::getline(is, name); + if (name.rfind('/') != std::string::npos) { + continue; // Unquoted name + } else { + is.seekg(pos); + break; // Read next keyword + } + } else if (ct.is(std::ctype_base::digit, c) || c== '.') { + is >> ignoreSlashLine; // Decimal digit. Ignore data. + continue; + } else if (c== '\'') { + is >> ignoreSlashLine; // Quote. Ignore data. + continue; + } else if(c == '-' && is.peek() == int('-')) { + is >> ignoreLine; // This line is a comment + continue; + } else if (c == '/' ) { + is >> ignoreLine; // This line is a null record. + continue; // (No data before slash) + } else { + is.putback(c); + std::string temp; + is >> temp; + std::cout << "READ ERROR! Next word is " << temp << std::endl; + } } } @@ -1964,29 +1964,29 @@ struct PLYVISC : public SpecialBase virtual void read(std::istream& is) { - // Note. This function assumes that NTPVT = 1, and reads only one table. - std::vector plyvisc; - readVectorData(is, plyvisc); - for (int i=0; i<(int)plyvisc.size(); i+=2) { - concentration_.push_back(plyvisc[i]); - factor_.push_back(plyvisc[i+1]); - } + // Note. This function assumes that NTPVT = 1, and reads only one table. + std::vector plyvisc; + readVectorData(is, plyvisc); + for (int i=0; i<(int)plyvisc.size(); i+=2) { + concentration_.push_back(plyvisc[i]); + factor_.push_back(plyvisc[i+1]); + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)concentration_.size(); ++i) { - os << concentration_[i] << " " << factor_[i] << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int i=0; i<(int)concentration_.size(); ++i) { + os << concentration_[i] << " " << factor_[i] << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)concentration_.size(); ++i) { - concentration_[i] *= units.polymer_density; - } + for (int i=0; i<(int)concentration_.size(); ++i) { + concentration_[i] *= units.polymer_density; + } } }; @@ -1999,29 +1999,29 @@ struct PLYROCK : public SpecialBase virtual void read(std::istream& is) { - // Note. This function assumes that NTSFUN = 1, and reads only one line. - plyrock_.resize(5,-1e00); - plyrock_[3] = 1; // Default value - int num_to_read = 5; - int num_read = readDefaultedVectorData(is, plyrock_, num_to_read); - if (num_read == num_to_read) { - ignoreSlashLine(is); - } + // Note. This function assumes that NTSFUN = 1, and reads only one line. + plyrock_.resize(5,-1e00); + plyrock_[3] = 1; // Default value + int num_to_read = 5; + int num_read = readDefaultedVectorData(is, plyrock_, num_to_read); + if (num_read == num_to_read) { + ignoreSlashLine(is); + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)plyrock_.size(); ++i) { - os << plyrock_[i] << " "; - } - os << "\n\n"; + os << name() << '\n'; + for (int i=0; i<(int)plyrock_.size(); ++i) { + os << plyrock_[i] << " "; + } + os << "\n\n"; } virtual void convertToSI(const EclipseUnits& units) { - - plyrock_[2] *= units.polymer_density; + + plyrock_[2] *= units.polymer_density; } }; @@ -2035,29 +2035,29 @@ struct PLYADS : public SpecialBase virtual void read(std::istream& is) { - // Note. This function assumes that NTSFUN = 1, and reads only one table. - std::vector plyads; - readVectorData(is, plyads); - for (int i=0; i<(int)plyads.size(); i+=2) { - local_concentration_.push_back(plyads[i]); - adsorbed_concentration_.push_back(plyads[i+1]); - } + // Note. This function assumes that NTSFUN = 1, and reads only one table. + std::vector plyads; + readVectorData(is, plyads); + for (int i=0; i<(int)plyads.size(); i+=2) { + local_concentration_.push_back(plyads[i]); + adsorbed_concentration_.push_back(plyads[i+1]); + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)local_concentration_.size(); ++i) { - os << local_concentration_[i] << " " << adsorbed_concentration_[i] << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int i=0; i<(int)local_concentration_.size(); ++i) { + os << local_concentration_[i] << " " << adsorbed_concentration_[i] << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)local_concentration_.size(); ++i) { - local_concentration_[i] *= units.polymer_density; - } + for (int i=0; i<(int)local_concentration_.size(); ++i) { + local_concentration_[i] *= units.polymer_density; + } } }; @@ -2070,24 +2070,24 @@ struct PLYMAX : public SpecialBase virtual void read(std::istream& is) { - // Note. This function assumes that NTMISC = 1, and reads only one line. - plymax_.resize(2); - readVectorData(is, plymax_); + // Note. This function assumes that NTMISC = 1, and reads only one line. + plymax_.resize(2); + readVectorData(is, plymax_); } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)plymax_.size(); ++i) { - os << plymax_[i] << " "; - } - os << "\n\n"; + os << name() << '\n'; + for (int i=0; i<(int)plymax_.size(); ++i) { + os << plymax_[i] << " "; + } + os << "\n\n"; } virtual void convertToSI(const EclipseUnits& units) { - plymax_[0] *= units.polymer_density; - plymax_[1] *= units.polymer_density; + plymax_[0] *= units.polymer_density; + plymax_[1] *= units.polymer_density; } }; @@ -2100,25 +2100,25 @@ struct TLMIXPAR : public SpecialBase virtual void read(std::istream& is) { - // Note. This function assumes that NTMISC = 1, and reads only one record. - tlmixpar_.resize(2, -1e100); - int num_to_read = 2; - int num_read = readDefaultedVectorData(is, tlmixpar_, num_to_read); - if (tlmixpar_[1] < 0) { - tlmixpar_[1] = tlmixpar_[0]; - } - if (num_read == num_to_read) { - ignoreSlashLine(is); - } + // Note. This function assumes that NTMISC = 1, and reads only one record. + tlmixpar_.resize(2, -1e100); + int num_to_read = 2; + int num_read = readDefaultedVectorData(is, tlmixpar_, num_to_read); + if (tlmixpar_[1] < 0) { + tlmixpar_[1] = tlmixpar_[0]; + } + if (num_read == num_to_read) { + ignoreSlashLine(is); + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)tlmixpar_.size(); ++i) { - os << tlmixpar_[i] << " "; - } - os << "\n\n"; + os << name() << '\n'; + for (int i=0; i<(int)tlmixpar_.size(); ++i) { + os << tlmixpar_[i] << " "; + } + os << "\n\n"; } virtual void convertToSI(const EclipseUnits& units) @@ -2138,8 +2138,8 @@ struct WpolymerLine WpolymerLine() { - well_ = polymer_group_ = salt_group_ = ""; - polymer_concentration_ = salt_concentration_ = 0.0; + well_ = polymer_group_ = salt_group_ = ""; + polymer_concentration_ = salt_concentration_ = 0.0; } }; @@ -2151,58 +2151,58 @@ struct WPOLYMER : public SpecialBase virtual void read(std::istream& is) { - while(is) { - std::string wellname = readString(is); - if (wellname[0] == '/') { - is >> ignoreLine; - break; - } - while (wellname.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - wellname = readString(is); - } - WpolymerLine wpolymer_line; - wpolymer_line.well_ = wellname; - is >> wpolymer_line.polymer_concentration_; - is >> wpolymer_line.salt_concentration_; - std::string group = readString(is); - if (group[0] == '/') { - is >> ignoreLine; - } else { - wpolymer_line.polymer_group_ = group; - group = readString(is); - if (group[0] == '/') { - is >> ignoreLine; - } else { - wpolymer_line.salt_group_ = group; - is >> ignoreLine; - } - } - wpolymer_.push_back(wpolymer_line); - } + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WpolymerLine wpolymer_line; + wpolymer_line.well_ = wellname; + is >> wpolymer_line.polymer_concentration_; + is >> wpolymer_line.salt_concentration_; + std::string group = readString(is); + if (group[0] == '/') { + is >> ignoreLine; + } else { + wpolymer_line.polymer_group_ = group; + group = readString(is); + if (group[0] == '/') { + is >> ignoreLine; + } else { + wpolymer_line.salt_group_ = group; + is >> ignoreLine; + } + } + wpolymer_.push_back(wpolymer_line); + } } virtual void write(std::ostream& os) const { - os << name() << '\n'; - for (int i=0; i<(int)wpolymer_.size(); ++i) { - os << wpolymer_[i].well_ - << " " << wpolymer_[i].polymer_concentration_ - << " " << wpolymer_[i].salt_concentration_ - << " " << wpolymer_[i].polymer_group_ - << " " << wpolymer_[i].salt_group_; - os << '\n'; - } - os << '\n'; + os << name() << '\n'; + for (int i=0; i<(int)wpolymer_.size(); ++i) { + os << wpolymer_[i].well_ + << " " << wpolymer_[i].polymer_concentration_ + << " " << wpolymer_[i].salt_concentration_ + << " " << wpolymer_[i].polymer_group_ + << " " << wpolymer_[i].salt_group_; + os << '\n'; + } + os << '\n'; } virtual void convertToSI(const EclipseUnits& units) { - for (int i=0; i<(int)wpolymer_.size(); ++i) { - wpolymer_[i].polymer_concentration_ *= units.polymer_density; - wpolymer_[i].salt_concentration_ *= units.polymer_density; - } + for (int i=0; i<(int)wpolymer_.size(); ++i) { + wpolymer_[i].polymer_concentration_ *= units.polymer_density; + wpolymer_[i].salt_concentration_ *= units.polymer_density; + } } }; @@ -2222,39 +2222,35 @@ struct ENDSCALE : public SpecialBase { n_tab_ = 1; n_node_ = 20; while (!is.eof()) { - std::string tmp = readString(is); - while (tmp.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - tmp = readString(is); - } - if (tmp[0] == '/') { - is >> ignoreLine; - break; - } + std::string tmp = readString(is); + while (tmp.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + tmp = readString(is); + } + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } dir_switch_ = tmp; - tmp = readString(is); - if (tmp[0] == '/') { - is >> ignoreLine; - break; - } + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } revers_switch_ = tmp; - is >> ignoreWhitespace; if(is.peek() == int('/')) { is >> ignoreLine; return; } is >> n_tab_; - is >> ignoreWhitespace; if(is.peek() == int('/')) { is >> ignoreLine; return; } is >> n_node_; - is >> ignoreLine; break; } @@ -2281,19 +2277,18 @@ struct SCALECRS : public SpecialBase { virtual void read(std::istream & is) { scalecrs_ = std::string("NO"); while (!is.eof()) { - std::string tmp = readString(is); - while (tmp.find("--") == 0) { - // This line is a comment - is >> ignoreLine; - tmp = readString(is); - } - if (tmp[0] == '/') { - is >> ignoreLine; - break; - } else if (tmp[0] == 'Y') { + std::string tmp = readString(is); + while (tmp.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + tmp = readString(is); + } + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } else if (tmp[0] == 'Y') { scalecrs_ = std::string("YES"); } - is >> ignoreLine; break; } From 7510ded183a3b583573ee1a2ddd7a1584bfa5a9c Mon Sep 17 00:00:00 2001 From: osae Date: Fri, 7 Dec 2012 14:34:06 +0100 Subject: [PATCH 3/5] Fixed loop index and added parameter check. --- opm/core/fluid/SatFuncSimple.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/opm/core/fluid/SatFuncSimple.cpp b/opm/core/fluid/SatFuncSimple.cpp index 37e74fff..9c42ad30 100644 --- a/opm/core/fluid/SatFuncSimple.cpp +++ b/opm/core/fluid/SatFuncSimple.cpp @@ -46,6 +46,9 @@ namespace Opm const std::vector& krw = swof_table[table_num][1]; const std::vector& krow = swof_table[table_num][2]; const std::vector& pcow = swof_table[table_num][3]; + if (krw.front() != 0.0 || krow.back() != 0.0) { + THROW("Error SWOF data - non-zero krw(swco) and/or krow(1-sor)"); + } buildUniformMonotoneTable(sw, krw, samples, krw_); buildUniformMonotoneTable(sw, krow, samples, krow_); buildUniformMonotoneTable(sw, pcow, samples, pcow_); @@ -68,10 +71,10 @@ namespace Opm break; } } - for (unsigned int i=sw.size()-2; i>=0; --i) { - if (krow[i]> 0.0) { - sowcr_ = 1.0 - sw[i+1]; - krwr_ = krw[i+1]; + for (unsigned int i=sw.size()-1; i>=1; --i) { + if (krow[i-1]> 0.0) { + sowcr_ = 1.0 - sw[i]; + krwr_ = krw[i]; break; } } @@ -268,6 +271,9 @@ namespace Opm const std::vector& krw = swof_table[table_num][1]; const std::vector& krow = swof_table[table_num][2]; const std::vector& pcow = swof_table[table_num][3]; + if (krw.front() != 0.0 || krow.back() != 0.0) { + THROW("Error SWOF data - non-zero krw(swco) and/or krow(1-sor)"); + } krw_ = NonuniformTableLinear(sw, krw); krow_ = NonuniformTableLinear(sw, krow); pcow_ = NonuniformTableLinear(sw, pcow); @@ -290,10 +296,10 @@ namespace Opm break; } } - for (unsigned int i=sw.size()-2; i>=0; --i) { - if (krow[i]> 0.0) { - sowcr_ = 1.0 - sw[i+1]; - krwr_ = krw[i+1]; + for (unsigned int i=sw.size()-1; i>=1; --i) { + if (krow[i-1]> 0.0) { + sowcr_ = 1.0 - sw[i]; + krwr_ = krw[i]; break; } } From cf5eab879415f7c1363a278d5d8869d099639b90 Mon Sep 17 00:00:00 2001 From: osae Date: Mon, 17 Dec 2012 13:58:20 +0100 Subject: [PATCH 4/5] Loop index. --- opm/core/fluid/SatFuncSimple.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/core/fluid/SatFuncSimple.cpp b/opm/core/fluid/SatFuncSimple.cpp index 9c42ad30..c4d98156 100644 --- a/opm/core/fluid/SatFuncSimple.cpp +++ b/opm/core/fluid/SatFuncSimple.cpp @@ -64,14 +64,14 @@ namespace Opm sowcr_ = 1.0 - swco; krwr_ = krw.back(); krorw_ = krow.front(); - for (unsigned int i=1; i::size_type i=1; i 0.0) { swcr_ = sw[i-1]; krorw_ = krow[i-1]; break; } } - for (unsigned int i=sw.size()-1; i>=1; --i) { + for (std::vector::size_type i=sw.size()-1; i>=1; --i) { if (krow[i-1]> 0.0) { sowcr_ = 1.0 - sw[i]; krwr_ = krw[i]; @@ -289,14 +289,14 @@ namespace Opm sowcr_ = 1.0 - swco; krwr_ = krw.back(); krorw_ = krow.front(); - for (unsigned int i=1; i::size_type i=1; i 0.0) { swcr_ = sw[i-1]; krorw_ = krow[i-1]; break; } } - for (unsigned int i=sw.size()-1; i>=1; --i) { + for (std::vector::size_type i=sw.size()-1; i>=1; --i) { if (krow[i-1]> 0.0) { sowcr_ = 1.0 - sw[i]; krwr_ = krw[i]; From d1aada8d81a7605d092edce8097823690029c60e Mon Sep 17 00:00:00 2001 From: osae Date: Mon, 17 Dec 2012 14:02:30 +0100 Subject: [PATCH 5/5] Support for keywords ENPTVD and ENKRVD. --- opm/core/eclipse/EclipseGridParser.cpp | 1 + opm/core/eclipse/EclipseGridParser.hpp | 2 + opm/core/eclipse/EclipseGridParserHelpers.hpp | 16 +- opm/core/eclipse/SpecialEclipseFields.hpp | 176 +++++++++++++++++- .../fluid/SaturationPropsFromDeck_impl.hpp | 160 ++++++++++------ 5 files changed, 292 insertions(+), 63 deletions(-) diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index fef38c44..7c7962b7 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -118,6 +118,7 @@ namespace EclipseKeywords string("PLYMAX"), string("TLMIXPAR"), string("WPOLYMER"), string("GRUPTREE"), string("GCONINJE"), string("GCONPROD"), string("WGRUPCON"), string("ENDSCALE"), string("SCALECRS"), + string("ENPTVD"), string("ENKRVD"), // The following fields only have a dummy implementation // that allows us to ignore them. string("SWFN"), diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index ab871c68..f09280ee 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -194,6 +194,8 @@ public: SPECIAL_FIELD(WGRUPCON) SPECIAL_FIELD(ENDSCALE) SPECIAL_FIELD(SCALECRS) + SPECIAL_FIELD(ENPTVD) + SPECIAL_FIELD(ENKRVD) // The following fields only have a dummy implementation // that allows us to ignore them. diff --git a/opm/core/eclipse/EclipseGridParserHelpers.hpp b/opm/core/eclipse/EclipseGridParserHelpers.hpp index 7bc0e003..20893eac 100644 --- a/opm/core/eclipse/EclipseGridParserHelpers.hpp +++ b/opm/core/eclipse/EclipseGridParserHelpers.hpp @@ -378,7 +378,10 @@ namespace } // Replace default values -1 by linear interpolation - inline void insertDefaultValues(std::vector >& table, int ncol) + inline void insertDefaultValues(std::vector >& table, + int ncol, + double defaultOut = 0.0, + bool defaultInterpolation = true) { const int sz = table[0].size(); for (int k=1; k > > table_; + std::vector mask_; + + virtual std::string name() const { + return std::string("ENPTVD"); + } + + virtual void read(std::istream & is) { + table_.resize(5); + std::vector > sub_table(5); + while (!is.eof()) { + if(is.peek() == int('/')) { + if (sub_table[0].empty() && !(table_[0].empty())) { + is >> ignoreLine; + break; + } else { + THROW("Error reading ENPTVD data - none or incomplete table."); + } + } + std::vector data(9,-1.0); + int nread = readDefaultedVectorData(is, data, 9); + if (nread != 9) { + THROW("Error reading ENPTVD data - depth and 8 saturations pr line."); + } + if (data[0] == -1.0) { + THROW("Error reading ENPTVD data - depth can not be defaulted."); + } + if ((data[4] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0) || (data[8] != -1.0)) { + THROW("Error reading ENPTVD data - non-default values in column 5-7,9 not supported."); + } + sub_table[0].push_back(data[0]); //depth + sub_table[1].push_back(data[1]); //swl + sub_table[2].push_back(data[2]); //swcr + sub_table[3].push_back(data[3]); //swu + sub_table[4].push_back(data[7]); //sowcr + is >> ignoreWhitespace; + if(is.peek() == int('/')) { + is >> ignoreLine; + if (sub_table[0].size() >= 2) { + insertDefaultValues(sub_table, 5, -1.0, false); + std::vector >::iterator it_sub = sub_table.begin(); + for(std::vector > >::size_type i=0; i > >::size_type i=1; i >::size_type j=0; j >::size_type j=0; j::size_type k=0; k > >::size_type i=0; i >::size_type j=0; j::size_type k=0; k > > table_; + std::vector mask_; + + virtual std::string name() const { + return std::string("ENKRVD"); + } + + virtual void read(std::istream & is) { + table_.resize(5); + std::vector > sub_table(5); + while (!is.eof()) { + if(is.peek() == int('/')) { + if (sub_table[0].empty() && !(table_[0].empty())) { + is >> ignoreLine; + break; + } else { + THROW("Error reading ENKRVD data - none or incomplete table."); + } + } + std::vector data(8,-1.0); + int nread = readDefaultedVectorData(is, data, 8); + if (nread != 8) { + THROW("Error reading ENKRVD data - depth and 7 relperms pr line."); + } + if (data[0] == -1.0) { + THROW("Error reading ENKRVD data - depth can not be defaulted."); + } + if ((data[2] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0)) { + THROW("Error reading ENKRVD data - non-default values in column 3,6-7 not supported."); + } + sub_table[0].push_back(data[0]); //depth + sub_table[1].push_back(data[1]); //krw + sub_table[2].push_back(data[3]); //kro + sub_table[3].push_back(data[4]); //krw(sowcr) + sub_table[4].push_back(data[7]); //kro(swcr) + is >> ignoreWhitespace; + if(is.peek() == int('/')) { + is >> ignoreLine; + if (sub_table[0].size() >= 2) { + insertDefaultValues(sub_table, 5, -1.0, false); + std::vector >::iterator it_sub = sub_table.begin(); + for(std::vector > >::size_type i=0; i > >::size_type i=1; i >::size_type j=0; j >::size_type j=0; j::size_type k=0; k > >::size_type i=0; i >::size_type j=0; j::size_type k=0; k void SaturationPropsFromDeck::initEPS(const EclipseGridParser& deck, - const UnstructuredGrid& grid, - const std::string& keyword, - std::vector& scaleparam) + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam) { - - if (deck.hasField(keyword)) { + bool useKeyword = deck.hasField(keyword); + bool hasENPTVD = deck.hasField("ENPTVD"); + bool hasENKRVD = deck.hasField("ENKRVD"); + int itab = 0; + std::vector > > table_dummy; + std::vector > >& table = table_dummy; - // Active keyword assigned default values for each cell (in case of possible box-wise assignment) - scaleparam.resize(grid.number_of_cells); - int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua]; - if (keyword == std::string("SWCR")) { - for (int i=0; i 0) { + table = deck.getENPTVD().table_; + } + } else if (keyword[0] == 'K' && (useKeyword || hasENKRVD)) { + if (keyword == std::string("KRW")) { + if (useKeyword || deck.getENKRVD().mask_[0]) { + itab = 1; + scaleparam.resize(grid.number_of_cells); + for (int i=0; i 0) { + table = deck.getENKRVD().table_; + } + } + if (scaleparam.empty()) { + return; + } else if (useKeyword) { // Keyword values from deck std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl; const int* gc = grid.global_cell; @@ -326,8 +350,28 @@ namespace Opm const int deck_pos = (gc == NULL) ? c : gc[c]; scaleparam[c] = val[deck_pos]; } + } else { + std::cout << "--- Scaling parameter '" << keyword << "' assigned via "; + if (keyword[0] == 'S') + deck.getENPTVD().write(std::cout); + else + deck.getENKRVD().write(std::cout); + const double* cc = grid.cell_centroids; + const int dim = grid.dimensions; + for (int cell = 0; cell < grid.number_of_cells; ++cell) { + int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell]; + if (table[itab][jtab][0] != -1.0) { + std::vector& depth = table[0][jtab]; + std::vector& val = table[itab][jtab]; + double zc = cc[dim*cell+dim-1]; + if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval + scaleparam[cell] = linearInterpolation(depth, val, zc); + } + } + } } } + // Saturation scaling template