From 4a81a0acf11f442de8f4fb8ef4a63d8d972dae4a Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 28 Jul 2015 17:28:51 +0200 Subject: [PATCH] SaturationPropsFromDeck: add the code to calculate all saturation properties using opm-material also remove the now obsolete code in that class. --- opm/core/props/BlackoilPropertiesFromDeck.cpp | 30 +- opm/core/props/BlackoilPropertiesFromDeck.hpp | 21 +- opm/core/props/IncompPropertiesFromDeck.cpp | 15 +- .../props/satfunc/SaturationPropsFromDeck.hpp | 85 +- .../satfunc/SaturationPropsFromDeck_impl.hpp | 770 +++--------------- 5 files changed, 177 insertions(+), 744 deletions(-) diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index d314849f7..91f529d71 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -28,7 +28,20 @@ namespace Opm const UnstructuredGrid& grid, bool init_rock) { - init(deck, eclState, grid.number_of_cells, grid.global_cell, grid.cartdims, + std::vector compressedToCartesianIdx(grid.number_of_cells); + for (unsigned cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) { + if (grid.global_cell) { + compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx]; + } + else { + compressedToCartesianIdx[cellIdx] = cellIdx; + } + } + + auto materialLawManager = std::make_shared(); + materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx); + + init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, grid.dimensions, init_rock); } @@ -38,7 +51,20 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock) { - init(deck, eclState, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, + std::vector compressedToCartesianIdx(grid.number_of_cells); + for (unsigned cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) { + if (grid.global_cell) { + compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx]; + } + else { + compressedToCartesianIdx[cellIdx] = cellIdx; + } + } + + auto materialLawManager = std::make_shared(); + materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx); + + init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, grid.dimensions, param, init_rock); } diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 8a4b1c339..487b01a4f 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -27,6 +27,8 @@ #include #include +#include + #include #include @@ -41,6 +43,8 @@ namespace Opm class BlackoilPropertiesFromDeck : public BlackoilPropertiesInterface { public: + typedef typename SaturationPropsFromDeck::MaterialLawManager MaterialLawManager; + /// Initialize from deck and grid. /// \param[in] deck Deck input parser /// \param[in] grid Grid to which property object applies, needed for the @@ -67,7 +71,6 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock=true); - template BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState, @@ -89,6 +92,18 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock=true); + template + BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck, + Opm::EclipseStateConstPtr eclState, + std::shared_ptr materialLawManager, + int number_of_cells, + const int* global_cell, + const int* cart_dims, + const CentroidIterator& begin_cell_centroids, + int dimension, + const parameter::ParameterGroup& param, + bool init_rock=true); + /// Destructor. virtual ~BlackoilPropertiesFromDeck(); @@ -240,6 +255,7 @@ namespace Opm template void init(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState, + std::shared_ptr materialLawManager, int number_of_cells, const int* global_cell, const int* cart_dims, @@ -249,6 +265,7 @@ namespace Opm template void init(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState, + std::shared_ptr materialLawManager, int number_of_cells, const int* global_cell, const int* cart_dims, @@ -256,9 +273,11 @@ namespace Opm int dimension, const parameter::ParameterGroup& param, bool init_rock); + RockFromDeck rock_; std::vector cellPvtRegionIdx_; BlackoilPvtProperties pvt_; + std::shared_ptr materialLawManager_; std::shared_ptr satprops_; mutable std::vector B_; mutable std::vector dB_; diff --git a/opm/core/props/IncompPropertiesFromDeck.cpp b/opm/core/props/IncompPropertiesFromDeck.cpp index b501582e1..a56298e37 100644 --- a/opm/core/props/IncompPropertiesFromDeck.cpp +++ b/opm/core/props/IncompPropertiesFromDeck.cpp @@ -32,7 +32,20 @@ namespace Opm { rock_.init(eclState, grid.number_of_cells, grid.global_cell, grid.cartdims); pvt_.init(deck); - satprops_.init(deck, eclState, grid); + auto materialLawManager = std::make_shared(); + + std::vector compressedToCartesianIdx(grid.number_of_cells); + for (unsigned cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) { + if (grid.global_cell) { + compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx]; + } + else { + compressedToCartesianIdx[cellIdx] = cellIdx; + } + } + materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx); + + satprops_.init(deck, eclState, materialLawManager, grid); if (pvt_.numPhases() != satprops_.numPhases()) { OPM_THROW(std::runtime_error, "IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp index 594ae5e70..fcaf76817 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp @@ -28,6 +28,10 @@ #include #include +#include +#include +#include + #include struct UnstructuredGrid; @@ -38,9 +42,23 @@ namespace Opm class SaturationPropsFromDeck : public SaturationPropsInterface { public: + typedef Opm::ThreePhaseMaterialTraits MaterialTraits; + typedef Opm::EclMaterialLawManager MaterialLawManager; + typedef MaterialLawManager::MaterialLaw MaterialLaw; + typedef MaterialLawManager::MaterialLawParams MaterialLawParams; + /// Default constructor. inline SaturationPropsFromDeck(); + /// Initialize from a MaterialLawManager object and a compressed to cartesian cell index map. + /// \param[in] materialLawManager An initialized MaterialLawManager object + inline void init(const PhaseUsage &phaseUsage, + std::shared_ptr materialLawManager); + + /// Initialize from deck and grid. /// \param[in] deck Deck input parser /// \param[in] grid Grid to which property object applies, needed for the @@ -48,6 +66,7 @@ namespace Opm /// to logical cartesian indices consistent with the deck. inline void init(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclipseState, + std::shared_ptr materialLawManager, const UnstructuredGrid& grid); /// Initialize from deck and grid. @@ -64,6 +83,7 @@ namespace Opm template inline void init(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclipseState, + std::shared_ptr materialLawManager, int number_of_cells, const int* global_cell, const T& begin_cell_centroids, @@ -127,69 +147,8 @@ namespace Opm double & swat); private: - // internal helper method for satRange() - template - void satRange_(const SaturationFunction& satFunc, - const int cellIdx, - const int* cells, - double* smin, - double* smax) const; - - PhaseUsage phase_usage_; - typedef Opm::SatFuncMultiplexer SatFunc; - std::vector satfunc_; - std::vector cell_to_func_; // = SATNUM - 1 - std::vector cell_to_func_imb_; - - bool do_eps_; // ENDSCALE is active - bool do_3pt_; // SCALECRS: YES~true NO~false - bool do_hyst_; // Keywords ISWL etc detected - std::vector eps_transf_; - std::vector eps_transf_hyst_; - std::vector sat_hyst_; - - template - void initEPS(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, - int number_of_cells, - const int* global_cell, - const T& begin_cell_centroids, - int dimensions, - const std::vector& eps_kw, - std::vector& eps_transf); - template - void initEPSKey(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, - int number_of_cells, - const int* global_cell, - const T& begin_cell_centroids, - int dimensions, - const std::string& keyword, - std::vector& scaleparam); - void initEPSParam(const int cell, - EPSTransforms::Transform& data, - const bool oil, - const double sl_tab, - const double scr_tab, - const double su_tab, - const double sxcr_tab, - const double s0_tab, - const double krsr_tab, - const double krmax_tab, - const double pcmax_tab, - const std::vector& sl, - const std::vector& scr, - const std::vector& su, - const std::vector& sxcr, - const std::vector& s0, - const std::vector& krsr, - const std::vector& krmax, - const std::vector& pcmax); - - bool columnIsMasked_(Opm::DeckConstPtr deck, - const std::string& keywordName, - int columnIdx) - { return deck->getKeyword(keywordName)->getRecord(columnIdx)->getItem(0)->getSIDouble(0) != -1.0; } + std::shared_ptr materialLawManager_; + PhaseUsage phaseUsage_; }; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index baa43198a..09cbae8e3 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -50,212 +50,42 @@ namespace Opm /// Initialize from deck. inline void SaturationPropsFromDeck::init(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, + Opm::EclipseStateConstPtr eclipseState, + std::shared_ptr materialLawManager, const UnstructuredGrid& grid) { - this->init(deck, eclipseState, grid.number_of_cells, + this->init(deck, eclipseState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cell_centroids, grid.dimensions); } + /// Initialize from deck. + inline + void SaturationPropsFromDeck::init(const PhaseUsage &phaseUsage, + std::shared_ptr materialLawManager) + { + phaseUsage_ = phaseUsage; + materialLawManager_ = materialLawManager; + } + /// Initialize from deck. template void SaturationPropsFromDeck::init(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, + Opm::EclipseStateConstPtr eclipseState, + std::shared_ptr materialLawManager, int number_of_cells, const int* global_cell, const T& begin_cell_centroids, int dimensions) { - phase_usage_ = phaseUsageFromDeck(deck); - - // Extract input data. - // Oil phase should be active. - if (!phase_usage_.phase_used[BlackoilPhases::Liquid]) { - OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active."); - } - - // Check SATOPTS status - bool hysteresis_switch = false; - if (deck->hasKeyword("SATOPTS")) { - const std::vector& satopts = deck->getKeyword("SATOPTS")->getStringData(); - for (size_t i = 0; i < satopts.size(); ++i) { - if (satopts[i] == std::string("HYSTER")) { - hysteresis_switch = true; - } else { - OPM_THROW(std::runtime_error, "Keyword SATOPTS: Switch " << satopts[i] << " not supported. "); - } - } - } - - // Obtain SATNUM, if it exists, and create cell_to_func_. - // Otherwise, let the cell_to_func_ mapping be just empty. - int satfuncs_expected = 1; - cell_to_func_.resize(number_of_cells, /*value=*/0); - if (deck->hasKeyword("SATNUM")) { - const std::vector& satnum = deck->getKeyword("SATNUM")->getIntData(); - satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); - const int num_cells = number_of_cells; - const int* gc = global_cell; - for (int cell = 0; cell < num_cells; ++cell) { - const int deck_pos = (gc == NULL) ? cell : gc[cell]; - cell_to_func_[cell] = satnum[deck_pos] - 1; - } - } - - // Find number of tables, check for consistency. - enum { Uninitialized = -1 }; - int num_tables = Uninitialized; - if (phase_usage_.phase_used[BlackoilPhases::Aqua]) { - num_tables = deck->getKeyword("SWOF")->size(); - if (num_tables < satfuncs_expected) { - OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected); - } - } - if (phase_usage_.phase_used[BlackoilPhases::Vapour]) { - int num_sgof_tables = deck->getKeyword("SGOF")->size(); - if (num_sgof_tables < satfuncs_expected) { - OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected); - } - if (num_tables == Uninitialized) { - num_tables = num_sgof_tables; - } else if (num_tables != num_sgof_tables) { - OPM_THROW(std::runtime_error, "Inconsistent number of tables in SWOF and SGOF."); - } - } - - // Initialize saturation function objects. - satfunc_.resize(num_tables); - SatFuncMultiplexer::SatFuncType satFuncType = SatFuncMultiplexer::Gwseg; - if (deck->hasKeyword("STONE2")) { - satFuncType = SatFuncMultiplexer::Stone2; - } - - for (int table = 0; table < num_tables; ++table) { - satfunc_[table].initFromDeck(eclipseState, table, satFuncType); - } - - // Check EHYSTR status - do_hyst_ = false; - if (hysteresis_switch && deck->hasKeyword("EHYSTR")) { - const int& relative_perm_hyst = deck->getKeyword("EHYSTR")->getRecord(0)->getItem(1)->getInt(0); - const std::string& limiting_hyst_flag = deck->getKeyword("EHYSTR")->getRecord(0)->getItem(4)->getString(0); - if (relative_perm_hyst != int(0)) { - OPM_THROW(std::runtime_error, "Keyword EHYSTR, item 2: Flag '" << relative_perm_hyst << "' found, only '0' is supported. "); - } - if (limiting_hyst_flag != std::string("KR")) { - OPM_THROW(std::runtime_error, "Keyword EHYSTR, item 5: Flag '" << limiting_hyst_flag << "' found, only 'KR' is supported. "); - } - if ( ! deck->hasKeyword("ENDSCALE")) { - // TODO When use of IMBNUM is implemented, this constraint will be lifted. - OPM_THROW(std::runtime_error, "Currently hysteris effects is only available through endpoint scaling."); - } - do_hyst_ = true; - } else if (hysteresis_switch) { - OPM_THROW(std::runtime_error, "Switch HYSTER of keyword SATOPTS is active, but keyword EHYSTR not found."); - } else if (deck->hasKeyword("EHYSTR")) { - OPM_THROW(std::runtime_error, "Found keyword EHYSTR, but switch HYSTER of keyword SATOPTS is not set."); - } - - // Saturation table scaling - do_eps_ = false; - do_3pt_ = false; - if (deck->hasKeyword("ENDSCALE")) { - Opm::EndscaleWrapper endscale(deck->getKeyword("ENDSCALE")); - if (endscale.directionSwitch() != std::string("NODIR")) { - OPM_THROW(std::runtime_error, - "SaturationPropsFromDeck::init() -- ENDSCALE: " - "Currently only 'NODIR' accepted."); - } - if (!endscale.isReversible()) { - OPM_THROW(std::runtime_error, - "SaturationPropsFromDeck::init() -- ENDSCALE: " - "Currently only 'REVERS' accepted."); - } - if (deck->hasKeyword("SCALECRS")) { - Opm::ScalecrsWrapper scalecrs(deck->getKeyword("SCALECRS")); - if (scalecrs.isEnabled()) { - do_3pt_ = true; - } - } - do_eps_ = true; - - // Make a consistency check of ENDNUM: #regions = NTENDP (ENDSCALE::3, TABDIMS::8)... - if (deck->hasKeyword("ENDNUM")) { - const std::vector& endnum = deck->getKeyword("ENDNUM")->getIntData(); - int endnum_regions = *std::max_element(endnum.begin(), endnum.end()); - if (endnum_regions > endscale.numEndscaleTables()) { - OPM_THROW(std::runtime_error, - "ENDNUM: Found " << endnum_regions << - " regions. Maximum allowed is " << endscale.numEndscaleTables() << - " (confer item 3 of keyword ENDSCALE)."); // TODO See also item 8 of TABDIMS ... - } - } - // TODO: ENPTVD/ENKRVD: Too few tables gives a cryptical message from parser, - // superfluous tables are ignored by the parser without any warning ... - - const std::vector eps_kw{"SWL", "SWU", "SWCR", "SGL", "SGU", "SGCR", "SOWCR", - "SOGCR", "KRW", "KRG", "KRO", "KRWR", "KRGR", "KRORW", "KRORG", "PCW", "PCG"}; - eps_transf_.resize(number_of_cells); - initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids, - dimensions, eps_kw, eps_transf_); - - if (do_hyst_) { - if (deck->hasKeyword("KRW") - || deck->hasKeyword("KRG") - || deck->hasKeyword("KRO") - || deck->hasKeyword("KRWR") - || deck->hasKeyword("KRGR") - || deck->hasKeyword("KRORW") - || deck->hasKeyword("KRORG") - || deck->hasKeyword("ENKRVD") - || deck->hasKeyword("IKRG") - || deck->hasKeyword("IKRO") - || deck->hasKeyword("IKRWR") - || deck->hasKeyword("IKRGR") - || deck->hasKeyword("IKRORW") - || deck->hasKeyword("IKRORG") ) { - OPM_THROW(std::runtime_error,"Currently hysteresis and relperm value scaling cannot be combined."); - } - - if (deck->hasKeyword("IMBNUM")) { - const std::vector& imbnum = deck->getKeyword("IMBNUM")->getIntData(); - int imbnum_regions = *std::max_element(imbnum.begin(), imbnum.end()); - if (imbnum_regions > num_tables) { - OPM_THROW(std::runtime_error, - "IMBNUM: Found " << imbnum_regions << - " regions. Maximum allowed is " << num_tables << - " (number of tables provided by SWOF/SGOF)."); - } - const int num_cells = number_of_cells; - cell_to_func_imb_.resize(num_cells); - const int* gc = global_cell; - for (int cell = 0; cell < num_cells; ++cell) { - const int deck_pos = (gc == NULL) ? cell : gc[cell]; - cell_to_func_imb_[cell] = imbnum[deck_pos] - 1; - } - // TODO: Make actual use of IMBNUM. For now we just consider the imbibition curve - // to be a scaled version of the drainage curve (confer Norne model). - } - - const std::vector eps_i_kw{"ISWL", "ISWU", "ISWCR", "ISGL", "ISGU", "ISGCR", "ISOWCR", - "ISOGCR", "IKRW", "IKRG", "IKRO", "IKRWR", "IKRGR", "IKRORW", "IKRORG", "IPCW", "IPCG"}; - eps_transf_hyst_.resize(number_of_cells); - sat_hyst_.resize(number_of_cells); - initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids, - dimensions, eps_i_kw, eps_transf_hyst_); - } - } + init(Opm::phaseUsageFromDeck(deck), materialLawManager); } - - - /// \return P, the number of phases. inline int SaturationPropsFromDeck::numPhases() const { - return phase_usage_.num_phases; + return phaseUsage_.num_phases; } @@ -280,32 +110,39 @@ namespace Opm { assert(cells != 0); - const int np = phase_usage_.num_phases; - ExplicitArraysFluidState fluidState(phase_usage_); - fluidState.setSaturationArray(s); - + const int np = BlackoilPhases::MaxNumPhases; if (dkrds) { -// #pragma omp parallel for + ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_); + fluidState.setSaturationArray(s); + + typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation; + Evaluation relativePerms[BlackoilPhases::MaxNumPhases]; for (int i = 0; i < n; ++i) { fluidState.setIndex(i); - if (do_hyst_) { - satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); - } else if (do_eps_) { - satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]])); - } else { - satfunc_[cell_to_func_[cells[i]]].evalKrDeriv(fluidState, kr + np*i, dkrds + np*np*i); + const auto& params = materialLawManager_->materialLawParams(cells[i]); + MaterialLaw::relativePermeabilities(relativePerms, params, fluidState); + + // copy the values calculated using opm-material to the target arrays + for (int krPhaseIdx = 0; krPhaseIdx < np; ++krPhaseIdx) { + kr[np*i + krPhaseIdx] = relativePerms[krPhaseIdx].value; + + for (int satPhaseIdx = 0; satPhaseIdx < np; ++satPhaseIdx) + dkrds[np*np*i + satPhaseIdx*np + krPhaseIdx] = relativePerms[krPhaseIdx].derivatives[satPhaseIdx]; } } } else { -// #pragma omp parallel for + ExplicitArraysFluidState fluidState(phaseUsage_); + fluidState.setSaturationArray(s); + + double relativePerms[BlackoilPhases::MaxNumPhases]; for (int i = 0; i < n; ++i) { fluidState.setIndex(i); - if (do_hyst_) { - satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); - } else if (do_eps_) { - satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i, &(eps_transf_[cells[i]])); - } else { - satfunc_[cell_to_func_[cells[i]]].evalKr(fluidState, kr + np*i); + const auto& params = materialLawManager_->materialLawParams(cells[i]); + MaterialLaw::relativePermeabilities(relativePerms, params, fluidState); + + // copy the values calculated using opm-material to the target arrays + for (int krPhaseIdx = 0; krPhaseIdx < np; ++krPhaseIdx) { + kr[np*i + krPhaseIdx] = relativePerms[krPhaseIdx]; } } } @@ -333,28 +170,41 @@ namespace Opm { assert(cells != 0); - const int np = phase_usage_.num_phases; - ExplicitArraysFluidState fluidState(phase_usage_); - fluidState.setSaturationArray(s); - + const int np = BlackoilPhases::MaxNumPhases; if (dpcds) { -// #pragma omp parallel for + ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_); + typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation; + fluidState.setSaturationArray(s); + + Evaluation capillaryPressures[BlackoilPhases::MaxNumPhases]; for (int i = 0; i < n; ++i) { fluidState.setIndex(i); - if (do_eps_) { - satfunc_[cell_to_func_[cells[i]]].evalPcDeriv(fluidState, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]])); - } else { - satfunc_[cell_to_func_[cells[i]]].evalPcDeriv(fluidState, pc + np*i, dpcds + np*np*i); + const auto& params = materialLawManager_->materialLawParams(cells[i]); + MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState); + + // copy the values calculated using opm-material to the target arrays + for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) { + double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0; + pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].value; + + for (int satPhaseIdx = 0; satPhaseIdx < np; ++satPhaseIdx) + dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].derivatives[satPhaseIdx]; } } } else { -// #pragma omp parallel for + ExplicitArraysFluidState fluidState(phaseUsage_); + fluidState.setSaturationArray(s); + + double capillaryPressures[BlackoilPhases::MaxNumPhases]; for (int i = 0; i < n; ++i) { fluidState.setIndex(i); - if (do_eps_) { - satfunc_[cell_to_func_[cells[i]]].evalPc(fluidState, pc + np*i, &(eps_transf_[cells[i]])); - } else { - satfunc_[cell_to_func_[cells[i]]].evalPc(fluidState, pc + np*i); + const auto& params = materialLawManager_->materialLawParams(cells[i]); + MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState); + + // copy the values calculated using opm-material to the target arrays + for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) { + double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0; + pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx]; } } } @@ -372,67 +222,26 @@ namespace Opm double* smin, double* smax) const { - for (int cellIdx = 0; cellIdx < n; ++cellIdx) { - const SatFuncMultiplexer& multiplexerSatFunc = satfunc_[cell_to_func_[cellIdx]]; - switch (multiplexerSatFunc.satFuncType()) { - case SatFuncMultiplexer::Gwseg: - satRange_(multiplexerSatFunc.getGwseg(), cellIdx, cells, smin, smax); - break; - case SatFuncMultiplexer::Stone2: - satRange_(multiplexerSatFunc.getStone2(), cellIdx, cells, smin, smax); - break; - case SatFuncMultiplexer::Simple: - satRange_(multiplexerSatFunc.getSimple(), cellIdx, cells, smin, smax); - break; - } + const int np = BlackoilPhases::MaxNumPhases; + for (int i = 0; i < n; ++i) { + const auto& scaledDrainageInfo = + materialLawManager_->oilWaterScaledEpsInfoDrainage(cells[i]); + + smin[np*i + BlackoilPhases::Aqua] = scaledDrainageInfo.Swl; + smax[np*i + BlackoilPhases::Aqua] = scaledDrainageInfo.Swu; + smin[np*i + BlackoilPhases::Vapour] = scaledDrainageInfo.Sgl; + smax[np*i + BlackoilPhases::Vapour] = scaledDrainageInfo.Sgu; + + smin[np*i + BlackoilPhases::Liquid] = 1.0; + smax[np*i + BlackoilPhases::Liquid] = 1.0; + smin[np*i + BlackoilPhases::Liquid] -= smax[np*i + BlackoilPhases::Aqua]; + smax[np*i + BlackoilPhases::Liquid] -= smin[np*i + BlackoilPhases::Aqua]; + smin[np*i + BlackoilPhases::Liquid] -= smax[np*i + BlackoilPhases::Vapour]; + smax[np*i + BlackoilPhases::Liquid] -= smin[np*i + BlackoilPhases::Vapour]; + smin[np*i + BlackoilPhases::Liquid] = std::max(0.0, smin[np*i + BlackoilPhases::Liquid]); } } - template - void SaturationPropsFromDeck::satRange_(const SaturationFunction& satFunc, - const int cellIdx, - const int* cells, - double* smin, - double* smax) const - { - assert(cells != 0); - const int np = phase_usage_.num_phases; - - if (do_eps_) { - const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; - const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; - const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; - - smin[np*cellIdx + opos] = 1.0; - smax[np*cellIdx + opos] = 1.0; - if (phase_usage_.phase_used[BlackoilPhases::Aqua]) { - smin[np*cellIdx + wpos] = eps_transf_[cells[cellIdx]].wat.doNotScale ? satFunc.smin_[wpos] - : eps_transf_[cells[cellIdx]].wat.smin; - smax[np*cellIdx + wpos] = eps_transf_[cells[cellIdx]].wat.doNotScale ? satFunc.smax_[wpos] - : eps_transf_[cells[cellIdx]].wat.smax; - smin[np*cellIdx + opos] -= smax[np*cellIdx + wpos]; - smax[np*cellIdx + opos] -= smin[np*cellIdx + wpos]; - } - if (phase_usage_.phase_used[BlackoilPhases::Vapour]) { - smin[np*cellIdx + gpos] = eps_transf_[cells[cellIdx]].gas.doNotScale ? satFunc.smin_[gpos] - : eps_transf_[cells[cellIdx]].gas.smin; - smax[np*cellIdx + gpos] = eps_transf_[cells[cellIdx]].gas.doNotScale ? satFunc.smax_[gpos] - : eps_transf_[cells[cellIdx]].gas.smax; - smin[np*cellIdx + opos] -= smax[np*cellIdx + gpos]; - smax[np*cellIdx + opos] -= smin[np*cellIdx + gpos]; - } - if (phase_usage_.phase_used[BlackoilPhases::Vapour] && phase_usage_.phase_used[BlackoilPhases::Aqua]) { - smin[np*cellIdx + opos] = std::max(0.0,smin[np*cellIdx + opos]); - } - } else { - for (int p = 0; p < np; ++p) { - smin[np*cellIdx + p] = satFunc.smin_[p]; - smax[np*cellIdx + p] = satFunc.smax_[p]; - } - } - } - - /// Update saturation state for the hysteresis tracking /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. @@ -443,13 +252,15 @@ namespace Opm { assert(cells != 0); - const int np = phase_usage_.num_phases; - if (do_hyst_) { -// #pragma omp parallel for + if (materialLawManager_->enableHysteresis()) { + ExplicitArraysFluidState fluidState(phaseUsage_); + fluidState.setSaturationArray(s); + for (int i = 0; i < n; ++i) { - satfunc_[cell_to_func_[cells[i]]].getSatFuncBase().updateSatHyst(s + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); + fluidState.setIndex(i); + materialLawManager_->updateHysteresis(fluidState, cells[i]); } - } + } } @@ -462,403 +273,8 @@ namespace Opm const double pcow, double& swat) { - if (phase_usage_.phase_used[BlackoilPhases::Aqua]) { - const double pc_low_threshold = 1.0e-8; - // TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74 - if (swat <= eps_transf_[cell].wat.smin) { - swat = eps_transf_[cell].wat.smin; - } else if (pcow < pc_low_threshold) { - swat = eps_transf_[cell].wat.smax; - } else { - const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; - const int max_np = BlackoilPhases::MaxNumPhases; - double s[max_np] = { 0.0 }; - s[wpos] = swat; - ExplicitArraysFluidState fluidState(phase_usage_); - fluidState.setSaturationArray(s); - fluidState.setIndex(0); - double pc[max_np] = { 0.0 }; - satfunc_[cell_to_func_[cell]].evalPc(fluidState, pc, &(eps_transf_[cell])); - if (pc[wpos] > pc_low_threshold) { - eps_transf_[cell].wat.pcFactor *= pcow/pc[wpos]; - } - } - } else { - OPM_THROW(std::runtime_error, "swatInitScaling: no water phase! "); - } + swat = materialLawManager_->applySwatinit(cell, pcow, swat); } - - - // Initialize saturation scaling parameters - template - void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, - int number_of_cells, - const int* global_cell, - const T& begin_cell_centroid, - int dimensions, - const std::vector& eps_kw, - std::vector& eps_transf) - { - std::vector > eps_vec(eps_kw.size()); - const std::vector dummy; - - for (size_t i = 0; i < eps_kw.size(); ++i) { - initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions, - eps_kw[i], eps_vec[i]); - } - - const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; - const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; - const bool oilWater = phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && !phase_usage_.phase_used[BlackoilPhases::Vapour]; - const bool oilGas = !phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && phase_usage_.phase_used[BlackoilPhases::Vapour]; - const bool threephase = phase_usage_.phase_used[BlackoilPhases::Aqua] && phase_usage_.phase_used[BlackoilPhases::Liquid] && phase_usage_.phase_used[BlackoilPhases::Vapour]; - - for (int cell = 0; cell < number_of_cells; ++cell) { - auto& satFuncBase = satfunc_[cell_to_func_[cell]].getSatFuncBase(); - if (threephase || oilWater) { - // ### krw - initEPSParam(cell, eps_transf[cell].wat, false, - satFuncBase.smin_[wpos], - satFuncBase.swcr_, - satFuncBase.smax_[wpos], - satFuncBase.sowcr_, - oilWater ? -1.0 : satFuncBase.smin_[gpos], - satFuncBase.krwr_, - satFuncBase.krwmax_, - satFuncBase.pcwmax_, - eps_vec[0], eps_vec[2], eps_vec[1], eps_vec[6], eps_vec[3], eps_vec[11], eps_vec[8], eps_vec[15]); - // ### krow - initEPSParam(cell, eps_transf[cell].watoil, true, - 0.0, - satFuncBase.sowcr_, - satFuncBase.smin_[wpos], - satFuncBase.swcr_, - oilWater ? -1.0 : satFuncBase.smin_[gpos], - satFuncBase.krorw_, - satFuncBase.kromax_, - 0.0, - eps_vec[0], eps_vec[6], eps_vec[0], eps_vec[2], eps_vec[3], eps_vec[13], eps_vec[10], dummy); - } - if (threephase || oilGas) { - // ### krg - initEPSParam(cell, eps_transf[cell].gas, false, - satFuncBase.smin_[gpos], - satFuncBase.sgcr_, - satFuncBase.smax_[gpos], - satFuncBase.sogcr_, - oilGas ? -1.0 : satFuncBase.smin_[wpos], - satFuncBase.krgr_, - satFuncBase.krgmax_, - satFuncBase.pcgmax_, - eps_vec[3], eps_vec[5], eps_vec[4], eps_vec[7], eps_vec[0], eps_vec[12], eps_vec[9], eps_vec[16]); - // ### krog - initEPSParam(cell, eps_transf[cell].gasoil, true, - 0.0, - satFuncBase.sogcr_, - satFuncBase.smin_[gpos], - satFuncBase.sgcr_, - oilGas ? -1.0 : satFuncBase.smin_[wpos], - satFuncBase.krorg_, - satFuncBase.kromax_, - 0.0, - eps_vec[3], eps_vec[7], eps_vec[3], eps_vec[5], eps_vec[0], eps_vec[14], eps_vec[10], dummy); - } - } - } - - // Initialize saturation scaling parameter - template - void SaturationPropsFromDeck::initEPSKey(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, - int number_of_cells, - const int* global_cell, - const T& begin_cell_centroid, - int dimensions, - const std::string& keyword, - std::vector& scaleparam) - { - const bool useAqua = phase_usage_.phase_used[BlackoilPhases::Aqua]; - const bool useLiquid = phase_usage_.phase_used[BlackoilPhases::Liquid]; - const bool useVapour = phase_usage_.phase_used[BlackoilPhases::Vapour]; - bool useKeyword = deck->hasKeyword(keyword); - bool useStateKeyword = eclipseState->hasDoubleGridProperty(keyword); - const std::map kw2tab = { - {"SWL", 1}, {"SWCR", 2}, {"SWU", 3}, {"SGL", 4}, - {"SGCR", 5}, {"SGU", 6}, {"SOWCR", 7}, {"SOGCR", 8}, - {"ISWL", 1}, {"ISWCR", 2}, {"ISWU", 3}, {"ISGL", 4}, - {"ISGCR", 5}, {"ISGU", 6}, {"ISOWCR", 7}, {"ISOGCR", 8}}; - bool hasENPTVD = deck->hasKeyword("ENPTVD"); - bool hasENKRVD = deck->hasKeyword("ENKRVD"); - int itab = 0; - std::vector > param_col; - std::vector > depth_col; - std::vector col_names; - std::shared_ptr tables = eclipseState->getTableManager(); - - // Active keyword assigned default values for each cell (in case of possible box-wise assignment) - if ((keyword[0] == 'S' && (useStateKeyword || hasENPTVD)) || (keyword[1] == 'S' && useStateKeyword) ) { - if (useAqua && (useStateKeyword || columnIsMasked_(deck, "ENPTVD", kw2tab.find(keyword)->second-1))) { - itab = kw2tab.find(keyword)->second; - scaleparam.resize(number_of_cells); - } - if (!useKeyword && itab > 0) { - const auto& enptvdTables = tables->getEnptvdTables(); - int num_tables = enptvdTables.size(); - param_col.resize(num_tables); - depth_col.resize(num_tables); - col_names.resize(9); - for (int table_num=0; table_num 0) { - const auto& enkrvdTables = tables->getEnkrvdTables(); - int num_tables = enkrvdTables.size(); - param_col.resize(num_tables); - depth_col.resize(num_tables); - col_names.resize(8); - for (int table_num=0; table_num val; - if (keyword[0] == 'S' || keyword[1] == 'S') { // Saturation from EclipseState - val = eclipseState->getDoubleGridProperty(keyword)->getData(); - } else { - val = deck->getKeyword(keyword)->getSIDoubleData(); //KR and PC directly from deck. - } - for (int c = 0; c < int(scaleparam.size()); ++c) { - const int deck_pos = (gc == NULL) ? c : gc[c]; - scaleparam[c] = val[deck_pos]; - } - } - - if (itab > 0) { - const int dim = dimensions; - std::vector endnum; - if ( deck->hasKeyword("ENDNUM")) { - const std::vector& e = - deck->getKeyword("ENDNUM")->getIntData(); - endnum.resize(number_of_cells); - const int* gc = global_cell; - for (int cell = 0; cell < number_of_cells; ++cell) { - const int deck_pos = (gc == NULL) ? cell : gc[cell]; - endnum[cell] = e[deck_pos] - 1; // Deck value zero prevents scaling via ENPTVD/ENKRVD - } - } - else { - // Default deck value is one - endnum.assign(number_of_cells, 0); - } - if (keyword[0] == 'S' || keyword[1] == 'S') { // From EclipseState - for (int cell = 0; cell < number_of_cells; ++cell) { - if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) { - double zc = UgGridHelpers - ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim), - dim-1); - if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval - scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc); - } - } else if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0) { - // As of 1/9-2014: Reflects remaining work on opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp ... - OPM_THROW(std::runtime_error, " -- Inconsistent EclipseState: '" << keyword << "' (ENPTVD)"); - } - } - } else { //KR and PC from deck. - for (int cell = 0; cell < number_of_cells; ++cell) { - if (endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) { - double zc = UgGridHelpers - ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim), - dim-1); - if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval - scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc); - } - } - } - } - } - -// std::cout << keyword << ":" << std::endl; -// for (int c = 0; c < int(scaleparam.size()); ++c) { -// std::cout << c << " " << scaleparam[c] << std::endl; -// } - - } - - // Saturation scaling - inline - void SaturationPropsFromDeck::initEPSParam(const int cell, - EPSTransforms::Transform& data, - const bool oil, // flag indicating krow/krog calculations - const double sl_tab, // minimum saturation (for krow/krog calculations this is normally zero) - const double scr_tab, // critical saturation - const double su_tab, // maximum saturation (for krow/krog calculations this is minimum water/gas saturation) - const double sxcr_tab, // second critical saturation (not used for 2pt scaling) - const double s0_tab, // threephase complementary minimum saturation (-1.0 indicates 2-phase) - const double krsr_tab, // relperm at displacing critical saturation - const double krmax_tab, // relperm at maximum saturation - const double pcmax_tab, // cap-pres at maximum saturation (zero => no scaling) - const std::vector& sl, // For krow/krog calculations this is not used - const std::vector& scr, - const std::vector& su, // For krow/krog calculations this is SWL/SGL - const std::vector& sxcr, - const std::vector& s0, - const std::vector& krsr, - const std::vector& krmax, - const std::vector& pcmax) // For krow/krog calculations this is not used - { - if (scr.empty() && su.empty() && (sxcr.empty() || !do_3pt_) && s0.empty()) { - data.doNotScale = true; - data.smin = sl_tab; - if (oil) { - data.smax = (s0_tab < 0.0) ? 1.0 - su_tab : 1.0 - su_tab - s0_tab; - } else { - data.smax = su_tab; - } - data.scr = scr_tab; - } else { - data.doNotScale = false; - data.do_3pt = do_3pt_; - double s_r; - if (s0_tab < 0.0) { // 2phase - s_r = 1.0-sxcr_tab; - if (do_3pt_) data.sr = sxcr.empty() ? s_r : 1.0-sxcr[cell]; - } else { // 3phase - s_r = 1.0-sxcr_tab-s0_tab; - if (do_3pt_)data.sr = 1.0 - (sxcr.empty() ? sxcr_tab : sxcr[cell]) - - (s0.empty() ? s0_tab : s0[cell]); - } - data.scr = scr.empty() ? scr_tab : scr[cell]; - double s_max = su_tab; - if (oil) { - data.smin = sl_tab; - if (s0_tab < 0.0) { // 2phase - s_max = 1.0 - su_tab; - data.smax = 1.0 - (su.empty() ? su_tab : su[cell]); - } else { // 3phase - s_max = 1.0 - su_tab - s0_tab; - data.smax = 1.0 - (su.empty() ? su_tab : su[cell]) - - (s0.empty() ? s0_tab : s0[cell]); - } - } else { - data.smin = sl.empty() ? sl_tab : sl[cell]; - data.smax = su.empty() ? su_tab : su[cell]; - } - if (do_3pt_) { - data.slope1 = (s_r-scr_tab)/(data.sr-data.scr); - data.slope2 = (s_max-s_r)/(data.smax-data.sr); - } else { - data.slope2 = data.slope1 = (s_max-scr_tab)/(data.smax-data.scr); - // Inv transform of tabulated critical displacing saturation to prepare for possible value scaling (krwr etc) - data.sr = data.scr + (s_r-scr_tab)*(data.smax-data.scr)/(s_max-scr_tab); - } - } - - data.doKrMax = !krmax.empty(); - data.doKrCrit = !krsr.empty(); - data.doSatInterp = false; - data.krsr = krsr.empty() ? krsr_tab : krsr[cell]; - data.krmax = krmax.empty() ? krmax_tab : krmax[cell]; - data.krSlopeCrit = data.krsr/krsr_tab; - data.krSlopeMax = data.krmax/krmax_tab; - if (data.doKrCrit) { - if (data.sr > data.smax-1.0e-6) { - //Ignore krsr and do two-point (one might consider combining krsr and krmax linearly between scr and smax ... ) - data.doKrCrit = false; - } else if (std::fabs(krmax_tab- krsr_tab) > 1.0e-6) { // interpolate kr - data.krSlopeMax = (data.krmax-data.krsr)/(krmax_tab-krsr_tab); - } else { // interpolate sat - data.doSatInterp = true; - data.krSlopeMax = (data.krmax-data.krsr)/(data.smax-data.sr); - } - } - - if (std::fabs(pcmax_tab) < 1.0e-8 || pcmax.empty() || pcmax_tab*pcmax[cell] < 0.0) { - data.pcFactor = 1.0; - } else { - data.pcFactor = pcmax[cell]/pcmax_tab; - } - - } - } // namespace Opm #endif // OPM_SATURATIONPROPSFROMDECK_IMPL_HEADER_INCLUDED