diff --git a/opm/models/blackoil/blackoilbioeffectsparams.cpp b/opm/models/blackoil/blackoilbioeffectsparams.cpp index 8b23ad3db68..23bed6538a2 100644 --- a/opm/models/blackoil/blackoilbioeffectsparams.cpp +++ b/opm/models/blackoil/blackoilbioeffectsparams.cpp @@ -34,12 +34,10 @@ #include -namespace Opm { +namespace { -template -template -void BlackOilBioeffectsParams:: -initFromState(const EclipseState& eclState) +template +void verifyState(const Opm::EclipseState& eclState) { // some sanity checks: // if biofilm is enabled, the BIOFILM keyword must be present, @@ -70,6 +68,18 @@ initFromState(const EclipseState& eclState) "contains the MICP keyword"); } } +} + +} + +namespace Opm { + +template +template +void BlackOilBioeffectsParams:: +initFromState(const EclipseState& eclState) +{ + verifyState(eclState); if (!eclState.runspec().micp() && !eclState.runspec().biof()) return; // bioeffects are supposed to be disabled diff --git a/opm/models/blackoil/blackoilbrineparams.cpp b/opm/models/blackoil/blackoilbrineparams.cpp index 7bf963c70ca..c05aee1a52c 100644 --- a/opm/models/blackoil/blackoilbrineparams.cpp +++ b/opm/models/blackoil/blackoilbrineparams.cpp @@ -36,27 +36,37 @@ #include #include -namespace Opm { +namespace { -template template -void BlackOilBrineParams:: -initFromState(const EclipseState& eclState) +void verifyState(const Opm::EclipseState& eclState) { // some sanity checks: if brine are enabled, the BRINE keyword must be // present, if brine are disabled the keyword must not be present. if constexpr (enableBrine) { - if (!eclState.runspec().phases().active(Phase::BRINE)) { + if (!eclState.runspec().phases().active(Opm::Phase::BRINE)) { throw std::runtime_error("Non-trivial brine treatment requested at compile time, but " "the deck does not contain the BRINE keyword"); } } else { - if (eclState.runspec().phases().active(Phase::BRINE)) { + if (eclState.runspec().phases().active(Opm::Phase::BRINE)) { throw std::runtime_error("Brine treatment disabled at compile time, but the deck " "contains the BRINE keyword"); } } +} + +} + +namespace Opm { + +template +template +void BlackOilBrineParams:: +initFromState(const EclipseState& eclState) +{ + verifyState(eclState); if (!eclState.runspec().phases().active(Phase::BRINE)) { return; // brine treatment is supposed to be disabled diff --git a/opm/models/blackoil/blackoilextboparams.cpp b/opm/models/blackoil/blackoilextboparams.cpp index b79e532d119..cfcbc5742a0 100644 --- a/opm/models/blackoil/blackoilextboparams.cpp +++ b/opm/models/blackoil/blackoilextboparams.cpp @@ -37,38 +37,65 @@ #include #include -namespace Opm { +namespace { -template template -void BlackOilExtboParams:: -initFromState(const EclipseState& eclState) +void verifyState(const Opm::EclipseState& eclState) { // some sanity checks: if extended BO is enabled, the PVTSOL keyword must be // present, if extended BO is disabled the keyword must not be present. if constexpr (enableExtbo) { - if (!eclState.runspec().phases().active(Phase::ZFRACTION)) { + if (!eclState.runspec().phases().active(Opm::Phase::ZFRACTION)) { throw std::runtime_error("Extended black oil treatment requested at compile " "time, but the deck does not contain the PVTSOL keyword"); } } else { - if (!enableExtbo && eclState.runspec().phases().active(Phase::ZFRACTION)) { + if (eclState.runspec().phases().active(Opm::Phase::ZFRACTION)) { throw std::runtime_error("Extended black oil treatment disabled at compile time, but the deck " "contains the PVTSOL keyword"); } } +} + +template +std::vector +parseSdensity(const Opm::EclipseState& eclState, + const std::size_t numPvtRegions) +{ + const auto& sdensityTables = eclState.getTableManager().getSolventDensityTables(); + if (sdensityTables.size() == numPvtRegions) { + std::vector zReferenceDensity(numPvtRegions); + for (std::size_t regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) { + const Scalar rhoRefS = sdensityTables[regionIdx].getSolventDensityColumn().front(); + zReferenceDensity[regionIdx] = rhoRefS; + } + return zReferenceDensity; + } + else { + throw std::runtime_error("Extbo: kw SDENSITY is missing or not aligned with NTPVT\n"); + } +} +} + +namespace Opm { + +template +template +void BlackOilExtboParams:: +initFromState(const EclipseState& eclState) +{ + verifyState(eclState); if (!eclState.runspec().phases().active(Phase::ZFRACTION)) { return; // solvent treatment is supposed to be disabled } // pvt properties from kw PVTSOL: - const auto& tableManager = eclState.getTableManager(); const auto& pvtsolTables = tableManager.getPvtsolTables(); - std::size_t numPvtRegions = pvtsolTables.size(); + const std::size_t numPvtRegions = pvtsolTables.size(); BO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); BG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); @@ -120,7 +147,7 @@ initFromState(const EclipseState& eclState) PBUB_RV_[regionIdx].appendXPos(ZCO2); const auto& underSaturatedTable = pvtsolTable.getUnderSaturatedTable(outerIdx); - std::size_t numRows = underSaturatedTable.numRows(); + const std::size_t numRows = underSaturatedTable.numRows(); Scalar bo0 = 0.0; Scalar po0 = 0.0; @@ -138,7 +165,7 @@ initFromState(const EclipseState& eclState) if (bo0 > bo) { // This is undersaturated oil-phase for ZCO2 <= zLim ... // Here we assume tabulated bo to decay beyond boiling point if (extractCmpFromPvt) { - Scalar cmpFactor = (bo - bo0) / (po - po0); + const Scalar cmpFactor = (bo - bo0) / (po - po0); oilCmp[outerIdx] = cmpFactor; zLim_[regionIdx] = ZCO2; //std::cout << "### cmpFactorOil: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl; @@ -147,9 +174,9 @@ initFromState(const EclipseState& eclState) } else if (bo0 == bo) { // This is undersaturated gas-phase for ZCO2 > zLim ... // Here we assume tabulated bo to be constant extrapolated beyond dew point if (innerIdx+1 < numRows && ZCO2<1.0 && extractCmpFromPvt) { - Scalar rvNxt = underSaturatedTable.get("RV", innerIdx + 1) + innerIdx * 1.0e-10; - Scalar bgNxt = underSaturatedTable.get("B_G", innerIdx + 1); - Scalar cmpFactor = (bgNxt - bg) / (rvNxt - rv); + const Scalar rvNxt = underSaturatedTable.get("RV", innerIdx + 1) + innerIdx * 1.0e-10; + const Scalar bgNxt = underSaturatedTable.get("B_G", innerIdx + 1); + const Scalar cmpFactor = (bgNxt - bg) / (rvNxt - rv); gasCmp[outerIdx] = cmpFactor; //std::cout << "### cmpFactorGas: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl; } @@ -189,17 +216,8 @@ initFromState(const EclipseState& eclState) gasCmp_[regionIdx].setXYContainers(zArg, gasCmp, /*sortInput=*/false); } - // Reference density for pure z-component taken from kw SDENSITY - const auto& sdensityTables = eclState.getTableManager().getSolventDensityTables(); - if (sdensityTables.size() == numPvtRegions) { - zReferenceDensity_.resize(numPvtRegions); - for (unsigned regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) { - Scalar rhoRefS = sdensityTables[regionIdx].getSolventDensityColumn().front(); - zReferenceDensity_[regionIdx] = rhoRefS; - } - } - else - throw std::runtime_error("Extbo: kw SDENSITY is missing or not aligned with NTPVT\n"); + // Reference density for pure z-component taken from kw SDENSITY + zReferenceDensity_ = parseSdensity(eclState, numPvtRegions); } #define INSTANTIATE_TYPE(T) \ diff --git a/opm/models/blackoil/blackoilpolymerparams.cpp b/opm/models/blackoil/blackoilpolymerparams.cpp index c76cece8d6e..4008cfa4930 100644 --- a/opm/models/blackoil/blackoilpolymerparams.cpp +++ b/opm/models/blackoil/blackoilpolymerparams.cpp @@ -31,6 +31,7 @@ copyright holders. #include #include +#include #include #include #include @@ -51,38 +52,32 @@ convertVecToVec(const std::vector>& input) return output; } -} - -namespace Opm { - -template template -void BlackOilPolymerParams:: -initFromState(const EclipseState& eclState) +void verifyState(const Opm::EclipseState& eclState) { // some sanity checks: if polymers are enabled, the POLYMER keyword must be // present, if polymers are disabled the keyword must not be present. if constexpr (enablePolymer) { - if (!eclState.runspec().phases().active(Phase::POLYMER)) { + if (!eclState.runspec().phases().active(Opm::Phase::POLYMER)) { throw std::runtime_error("Non-trivial polymer treatment requested at compile time, but " "the deck does not contain the POLYMER keyword"); } } else { - if (eclState.runspec().phases().active(Phase::POLYMER)) { + if (eclState.runspec().phases().active(Opm::Phase::POLYMER)) { throw std::runtime_error("Polymer treatment disabled at compile time, but the deck " "contains the POLYMER keyword"); } } if constexpr (enablePolymerMolarWeight) { - if (!eclState.runspec().phases().active(Phase::POLYMW)) { + if (!eclState.runspec().phases().active(Opm::Phase::POLYMW)) { throw std::runtime_error("Polymer molecular weight tracking is enabled at compile time, but " "the deck does not contain the POLYMW keyword"); } } else { - if (eclState.runspec().phases().active(Phase::POLYMW)) { + if (eclState.runspec().phases().active(Opm::Phase::POLYMW)) { throw std::runtime_error("Polymer molecular weight tracking is disabled at compile time, but the deck " "contains the POLYMW keyword"); } @@ -92,7 +87,18 @@ initFromState(const EclipseState& eclState) throw std::runtime_error("Polymer molecular weight tracking is enabled while polymer treatment " "is disabled at compile time"); } +} + +} + +namespace Opm { +template +template +void BlackOilPolymerParams:: +initFromState(const EclipseState& eclState) +{ + verifyState(eclState); if (!eclState.runspec().phases().active(Phase::POLYMER)) { return; // polymer treatment is supposed to be disabled } @@ -102,7 +108,79 @@ initFromState(const EclipseState& eclState) unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); setNumSatRegions(numSatRegions); - // initialize the objects which deal with the PLYROCK keyword + // initialize the objects which deal with the PLYROCK keyword + processPlyrock(eclState); + + // initialize the objects which deal with the PLYADS keyword + processPlyads(eclState); + + const unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + plyviscViscosityMultiplierTable_.resize(numPvtRegions); + + // initialize the objects which deal with the PLYVISC keyword + processPlyvisc(eclState); + + // initialize the objects which deal with the PLYMAX keyword + processPlymax(eclState); + + processPlmixpar(eclState); + + hasPlyshlog_ = tableManager.hasTables("PLYSHLOG"); + hasShrate_ = tableManager.useShrate(); + + if ((hasPlyshlog_ || hasShrate_) && enablePolymerMolarWeight) { + OpmLog::warning("PLYSHLOG and SHRATE should not be used in POLYMW runs, they will have no effect.\n"); + } + + if (hasPlyshlog_ && !enablePolymerMolarWeight) { + processPlyshlog(eclState); + } + + if (hasShrate_ && !enablePolymerMolarWeight) { + if (!hasPlyshlog_) { + throw std::runtime_error("PLYSHLOG must be specified if SHRATE is used in POLYMER runs\n"); + } + processShrate(eclState); + } + + if constexpr (enablePolymerMolarWeight) { + processPlyvmh(eclState); + processPlymwinj(eclState); + processSkprwat(eclState); + processSkprpoly(eclState); + } +} + +template +void BlackOilPolymerParams:: +setNumSatRegions(unsigned numRegions) +{ + plyrockDeadPoreVolume_.resize(numRegions); + plyrockResidualResistanceFactor_.resize(numRegions); + plyrockRockDensityFactor_.resize(numRegions); + plyrockAdsorbtionIndex_.resize(numRegions); + plyrockMaxAdsorbtion_.resize(numRegions); + plyadsAdsorbedPolymer_.resize(numRegions); +} + +template +void BlackOilPolymerParams:: +setNumMixRegions(unsigned numRegions, bool enablePolymerMolarWeight) +{ + plymaxMaxConcentration_.resize(numRegions); + plymixparToddLongstaff_.resize(numRegions); + + if (enablePolymerMolarWeight) { + plyvmhCoefficients_.resize(numRegions); + } +} + +template +void BlackOilPolymerParams:: +processPlyrock(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); const auto& plyrockTables = tableManager.getPlyrockTables(); if (!plyrockTables.empty()) { assert(numSatRegions == plyrockTables.size()); @@ -119,8 +197,14 @@ initFromState(const EclipseState& eclState) else { throw std::runtime_error("PLYROCK must be specified in POLYMER runs\n"); } +} - // initialize the objects which deal with the PLYADS keyword +template +void BlackOilPolymerParams:: +processPlyads(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); const auto& plyadsTables = tableManager.getPlyadsTables(); if (!plyadsTables.empty()) { assert(numSatRegions == plyadsTables.size()); @@ -135,13 +219,16 @@ initFromState(const EclipseState& eclState) else { throw std::runtime_error("PLYADS must be specified in POLYMER runs\n"); } +} - - unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); - plyviscViscosityMultiplierTable_.resize(numPvtRegions); - - // initialize the objects which deal with the PLYVISC keyword +template +template +void BlackOilPolymerParams:: +processPlyvisc(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); const auto& plyviscTables = tableManager.getPlyviscTables(); + const unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); if (!plyviscTables.empty()) { // different viscosity model is used for POLYMW if constexpr (enablePolymerMolarWeight) { @@ -159,11 +246,19 @@ initFromState(const EclipseState& eclState) } } } - else if constexpr (!enablePolymerMolarWeight) { - throw std::runtime_error("PLYVISC must be specified in POLYMER runs\n"); + else { + if constexpr (!enablePolymerMolarWeight) { + throw std::runtime_error("PLYVISC must be specified in POLYMER runs\n"); + } } +} - // initialize the objects which deal with the PLYMAX keyword +template +template +void BlackOilPolymerParams:: +processPlymax(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); const auto& plymaxTables = tableManager.getPlymaxTables(); const unsigned numMixRegions = plymaxTables.size(); setNumMixRegions(numMixRegions, enablePolymerMolarWeight); @@ -176,193 +271,203 @@ initFromState(const EclipseState& eclState) else { throw std::runtime_error("PLYMAX must be specified in POLYMER runs\n"); } +} - if (!eclState.getTableManager().getPlmixparTable().empty()) { +template +template +void BlackOilPolymerParams:: +processPlmixpar(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const auto& plmixparTable = tableManager.getPlmixparTable(); + const auto& plymaxTables = tableManager.getPlymaxTables(); + const unsigned numMixRegions = plymaxTables.size(); + if (!plmixparTable.empty()) { if constexpr (enablePolymerMolarWeight) { OpmLog::warning("PLMIXPAR should not be used in POLYMW runs, it will have no effect.\n"); } else { - const auto& plmixparTable = eclState.getTableManager().getPlmixparTable(); // initialize the objects which deal with the PLMIXPAR keyword for (unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++mixRegionIdx) { plymixparToddLongstaff_[mixRegionIdx] = plmixparTable[mixRegionIdx].todd_langstaff; } } } - else if constexpr (!enablePolymerMolarWeight) { - throw std::runtime_error("PLMIXPAR must be specified in POLYMER runs\n"); - } - - hasPlyshlog_ = eclState.getTableManager().hasTables("PLYSHLOG"); - hasShrate_ = eclState.getTableManager().useShrate(); - - if ((hasPlyshlog_ || hasShrate_) && enablePolymerMolarWeight) { - OpmLog::warning("PLYSHLOG and SHRATE should not be used in POLYMW runs, they will have no effect.\n"); - } - - if (hasPlyshlog_ && !enablePolymerMolarWeight) { - const auto& plyshlogTables = tableManager.getPlyshlogTables(); - assert(numPvtRegions == plyshlogTables.size()); - plyshlogShearEffectRefMultiplier_.resize(numPvtRegions); - plyshlogShearEffectRefLogVelocity_.resize(numPvtRegions); - for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { - const auto& plyshlogTable = plyshlogTables.template getTable(pvtRegionIdx); - - Scalar plyshlogRefPolymerConcentration = plyshlogTable.getRefPolymerConcentration(); - auto waterVelocity = plyshlogTable.getWaterVelocityColumn().vectorCopy(); - auto shearMultiplier = plyshlogTable.getShearMultiplierColumn().vectorCopy(); - - // do the unit version here for the waterVelocity - UnitSystem unitSystem = eclState.getDeckUnitSystem(); - double siFactor = hasShrate_? unitSystem.parse("1/Time").getSIScaling() : unitSystem.parse("Length/Time").getSIScaling(); - for (std::size_t i = 0; i < waterVelocity.size(); ++i) { - waterVelocity[i] *= siFactor; - // for plyshlog the input must be stored as logarithms - // the interpolation is then done the log-space. - waterVelocity[i] = std::log(waterVelocity[i]); - } - - Scalar refViscMult = plyviscViscosityMultiplierTable_[pvtRegionIdx].eval(plyshlogRefPolymerConcentration, /*extrapolate=*/true); - // convert the table using referece conditions - for (std::size_t i = 0; i < waterVelocity.size(); ++i) { - shearMultiplier[i] *= refViscMult; - shearMultiplier[i] -= 1; - shearMultiplier[i] /= (refViscMult - 1); - } - plyshlogShearEffectRefMultiplier_[pvtRegionIdx].resize(waterVelocity.size()); - plyshlogShearEffectRefLogVelocity_[pvtRegionIdx].resize(waterVelocity.size()); - - for (std::size_t i = 0; i < waterVelocity.size(); ++i) { - plyshlogShearEffectRefMultiplier_[pvtRegionIdx][i] = shearMultiplier[i]; - plyshlogShearEffectRefLogVelocity_[pvtRegionIdx][i] = waterVelocity[i]; - } + else { + if constexpr (!enablePolymerMolarWeight) { + throw std::runtime_error("PLMIXPAR must be specified in POLYMER runs\n"); } } +} - if (hasShrate_ && !enablePolymerMolarWeight) { - if (!hasPlyshlog_) { - throw std::runtime_error("PLYSHLOG must be specified if SHRATE is used in POLYMER runs\n"); - } - const auto& shrateTable = eclState.getTableManager().getShrateTable(); - shrate_.resize(numPvtRegions); - for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { - if (shrateTable.empty()) { - shrate_[pvtRegionIdx] = 4.8; //default; - } - else if (shrateTable.size() == numPvtRegions) { - shrate_[pvtRegionIdx] = shrateTable[pvtRegionIdx].rate; - } - else { - throw std::runtime_error("SHRATE must either have 0 or number of NUMPVT entries\n"); - } +template +void BlackOilPolymerParams:: +processPlyshlog(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const auto& plyshlogTables = tableManager.getPlyshlogTables(); + const unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + assert(numPvtRegions == plyshlogTables.size()); + plyshlogShearEffectRefMultiplier_.resize(numPvtRegions); + plyshlogShearEffectRefLogVelocity_.resize(numPvtRegions); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { + const auto& plyshlogTable = plyshlogTables.template getTable(pvtRegionIdx); + + const Scalar plyshlogRefPolymerConcentration = plyshlogTable.getRefPolymerConcentration(); + auto waterVelocity = plyshlogTable.getWaterVelocityColumn().vectorCopy(); + auto shearMultiplier = plyshlogTable.getShearMultiplierColumn().vectorCopy(); + + // do the unit version here for the waterVelocity + UnitSystem unitSystem = eclState.getDeckUnitSystem(); + const double siFactor = hasShrate_? unitSystem.parse("1/Time").getSIScaling() + : unitSystem.parse("Length/Time").getSIScaling(); + + // for plyshlog the input must be stored as logarithms + // the interpolation is then done in log-space. + std::ranges::transform(waterVelocity, waterVelocity.begin(), + [siFactor](const double input) + { return std::log(input * siFactor); }); + + const Scalar refViscMult = plyviscViscosityMultiplierTable_[pvtRegionIdx]. + eval(plyshlogRefPolymerConcentration, /*extrapolate=*/true); + + // convert the table using reference conditions + std::ranges::transform(shearMultiplier, shearMultiplier.begin(), + [refViscMult](const double input) + { return (input * refViscMult - 1.0) / (refViscMult - 1.0); }); + + plyshlogShearEffectRefMultiplier_[pvtRegionIdx].resize(waterVelocity.size()); + plyshlogShearEffectRefLogVelocity_[pvtRegionIdx].resize(waterVelocity.size()); + + for (std::size_t i = 0; i < waterVelocity.size(); ++i) { + plyshlogShearEffectRefMultiplier_[pvtRegionIdx][i] = shearMultiplier[i]; + plyshlogShearEffectRefLogVelocity_[pvtRegionIdx][i] = waterVelocity[i]; } } +} - if constexpr (enablePolymerMolarWeight) { - const auto& plyvmhTable = eclState.getTableManager().getPlyvmhTable(); - if (!plyvmhTable.empty()) { - assert(plyvmhTable.size() == numMixRegions); - for (std::size_t regionIdx = 0; regionIdx < numMixRegions; ++regionIdx) { - plyvmhCoefficients_[regionIdx].k_mh = plyvmhTable[regionIdx].k_mh; - plyvmhCoefficients_[regionIdx].a_mh = plyvmhTable[regionIdx].a_mh; - plyvmhCoefficients_[regionIdx].gamma = plyvmhTable[regionIdx].gamma; - plyvmhCoefficients_[regionIdx].kappa = plyvmhTable[regionIdx].kappa; - } +template +void BlackOilPolymerParams:: +processShrate(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const auto& shrateTable = tableManager.getShrateTable(); + const unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + shrate_.resize(numPvtRegions); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { + if (shrateTable.empty()) { + shrate_[pvtRegionIdx] = 4.8; //default; } - else { - throw std::runtime_error("PLYVMH keyword must be specified in POLYMW rus \n"); + else if (shrateTable.size() == numPvtRegions) { + shrate_[pvtRegionIdx] = shrateTable[pvtRegionIdx].rate; } - - // handling PLYMWINJ keyword - const auto& plymwinjTables = tableManager.getPlymwinjTables(); - for (const auto& table : plymwinjTables) { - const int tableNumber = table.first; - const auto& plymwinjtable = table.second; - const std::vector& throughput = plymwinjtable.getThroughputs(); - const std::vector& watervelocity = plymwinjtable.getVelocities(); - const std::vector>& molecularweight = plymwinjtable.getMoleWeights(); - if constexpr (std::is_same_v) { - const std::vector tp(throughput.begin(), throughput.end()); - const std::vector wv(watervelocity.begin(), watervelocity.end()); - const auto mw = convertVecToVec(molecularweight); - TabulatedTwoDFunction tablefunc(tp, wv, mw, true, false); - plymwinjTables_[tableNumber] = std::move(tablefunc); - } else { - TabulatedTwoDFunction tablefunc(throughput, watervelocity, molecularweight, true, false); - plymwinjTables_[tableNumber] = std::move(tablefunc); - } + else { + throw std::runtime_error("SHRATE must either have 0 or number of NUMPVT entries\n"); } + } +} - // handling SKPRWAT keyword - const auto& skprwatTables = tableManager.getSkprwatTables(); - for (const auto& table : skprwatTables) { - const int tableNumber = table.first; - const auto& skprwattable = table.second; - const std::vector& throughput = skprwattable.getThroughputs(); - const std::vector& watervelocity = skprwattable.getVelocities(); - const std::vector>& skinpressure = skprwattable.getSkinPressures(); - if constexpr (std::is_same_v) { - const std::vector tp(throughput.begin(), throughput.end()); - const std::vector wv(watervelocity.begin(), watervelocity.end()); - const auto sp = convertVecToVec(skinpressure); - TabulatedTwoDFunction tablefunc(tp, wv, sp, true, false); - skprwatTables_[tableNumber] = std::move(tablefunc); - } else { - TabulatedTwoDFunction tablefunc(throughput, watervelocity, skinpressure, true, false); - skprwatTables_[tableNumber] = std::move(tablefunc); - } +template +void BlackOilPolymerParams:: +processPlyvmh(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const auto& plyvmhTable = tableManager.getPlyvmhTable(); + const auto& plymaxTables = tableManager.getPlymaxTables(); + const unsigned numMixRegions = plymaxTables.size(); + if (!plyvmhTable.empty()) { + assert(plyvmhTable.size() == numMixRegions); + for (std::size_t regionIdx = 0; regionIdx < numMixRegions; ++regionIdx) { + plyvmhCoefficients_[regionIdx].k_mh = plyvmhTable[regionIdx].k_mh; + plyvmhCoefficients_[regionIdx].a_mh = plyvmhTable[regionIdx].a_mh; + plyvmhCoefficients_[regionIdx].gamma = plyvmhTable[regionIdx].gamma; + plyvmhCoefficients_[regionIdx].kappa = plyvmhTable[regionIdx].kappa; } + } + else { + throw std::runtime_error("PLYVMH keyword must be specified in POLYMW runs\n"); + } +} - // handling SKPRPOLY keyword - const auto& skprpolyTables = tableManager.getSkprpolyTables(); - for (const auto& table : skprpolyTables) { - const int tableNumber = table.first; - const auto& skprpolytable = table.second; - const std::vector& throughput = skprpolytable.getThroughputs(); - const std::vector& watervelocity = skprpolytable.getVelocities(); - const std::vector>& skinpressure = skprpolytable.getSkinPressures(); - const double refPolymerConcentration = skprpolytable.referenceConcentration(); - if constexpr (std::is_same_v) { - const std::vector tp(throughput.begin(), throughput.end()); - const std::vector wv(watervelocity.begin(), watervelocity.end()); - const auto sp = convertVecToVec(skinpressure); - SkprpolyTable tablefunc { - refPolymerConcentration, - TabulatedTwoDFunction(tp, wv, sp, true, false) - }; - skprpolyTables_[tableNumber] = std::move(tablefunc); - } else { - SkprpolyTable tablefunc { - refPolymerConcentration, - TabulatedTwoDFunction(throughput, watervelocity, skinpressure, true, false) - }; - skprpolyTables_[tableNumber] = std::move(tablefunc); - } +template +void BlackOilPolymerParams:: +processPlymwinj(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const auto& plymwinjTables = tableManager.getPlymwinjTables(); + for (const auto& table : plymwinjTables) { + const int tableNumber = table.first; + const auto& plymwinjtable = table.second; + const std::vector& throughput = plymwinjtable.getThroughputs(); + const std::vector& watervelocity = plymwinjtable.getVelocities(); + const std::vector>& molecularweight = plymwinjtable.getMoleWeights(); + if constexpr (std::is_same_v) { + const std::vector tp(throughput.begin(), throughput.end()); + const std::vector wv(watervelocity.begin(), watervelocity.end()); + const auto mw = convertVecToVec(molecularweight); + TabulatedTwoDFunction tablefunc(tp, wv, mw, true, false); + plymwinjTables_[tableNumber] = std::move(tablefunc); + } else { + TabulatedTwoDFunction tablefunc(throughput, watervelocity, molecularweight, true, false); + plymwinjTables_[tableNumber] = std::move(tablefunc); } } } template void BlackOilPolymerParams:: -setNumSatRegions(unsigned numRegions) +processSkprwat(const EclipseState& eclState) { - plyrockDeadPoreVolume_.resize(numRegions); - plyrockResidualResistanceFactor_.resize(numRegions); - plyrockRockDensityFactor_.resize(numRegions); - plyrockAdsorbtionIndex_.resize(numRegions); - plyrockMaxAdsorbtion_.resize(numRegions); - plyadsAdsorbedPolymer_.resize(numRegions); + const auto& tableManager = eclState.getTableManager(); + const auto& skprwatTables = tableManager.getSkprwatTables(); + for (const auto& table : skprwatTables) { + const int tableNumber = table.first; + const auto& skprwattable = table.second; + const std::vector& throughput = skprwattable.getThroughputs(); + const std::vector& watervelocity = skprwattable.getVelocities(); + const std::vector>& skinpressure = skprwattable.getSkinPressures(); + if constexpr (std::is_same_v) { + const std::vector tp(throughput.begin(), throughput.end()); + const std::vector wv(watervelocity.begin(), watervelocity.end()); + const auto sp = convertVecToVec(skinpressure); + TabulatedTwoDFunction tablefunc(tp, wv, sp, true, false); + skprwatTables_[tableNumber] = std::move(tablefunc); + } else { + TabulatedTwoDFunction tablefunc(throughput, watervelocity, skinpressure, true, false); + skprwatTables_[tableNumber] = std::move(tablefunc); + } + } } template void BlackOilPolymerParams:: -setNumMixRegions(unsigned numRegions, bool enablePolymerMolarWeight) +processSkprpoly(const EclipseState& eclState) { - plymaxMaxConcentration_.resize(numRegions); - plymixparToddLongstaff_.resize(numRegions); - - if (enablePolymerMolarWeight) { - plyvmhCoefficients_.resize(numRegions); + const auto& tableManager = eclState.getTableManager(); + const auto& skprpolyTables = tableManager.getSkprpolyTables(); + for (const auto& table : skprpolyTables) { + const int tableNumber = table.first; + const auto& skprpolytable = table.second; + const std::vector& throughput = skprpolytable.getThroughputs(); + const std::vector& watervelocity = skprpolytable.getVelocities(); + const std::vector>& skinpressure = skprpolytable.getSkinPressures(); + const double refPolymerConcentration = skprpolytable.referenceConcentration(); + if constexpr (std::is_same_v) { + const std::vector tp(throughput.begin(), throughput.end()); + const std::vector wv(watervelocity.begin(), watervelocity.end()); + const auto sp = convertVecToVec(skinpressure); + SkprpolyTable tablefunc { + refPolymerConcentration, + TabulatedTwoDFunction(tp, wv, sp, true, false) + }; + skprpolyTables_[tableNumber] = std::move(tablefunc); + } else { + SkprpolyTable tablefunc { + refPolymerConcentration, + TabulatedTwoDFunction(throughput, watervelocity, skinpressure, true, false) + }; + skprpolyTables_[tableNumber] = std::move(tablefunc); + } } } @@ -382,6 +487,7 @@ setPlyrock(unsigned satRegionIdx, plyrockMaxAdsorbtion_[satRegionIdx] = plyrockMaxAdsorbtion; } + #define INSTANTIATE_TYPE(T) \ template struct BlackOilPolymerParams; \ template void BlackOilPolymerParams::initFromState(const EclipseState&); \ diff --git a/opm/models/blackoil/blackoilpolymerparams.hpp b/opm/models/blackoil/blackoilpolymerparams.hpp index 56ac139f6ca..b062c7e0273 100644 --- a/opm/models/blackoil/blackoilpolymerparams.hpp +++ b/opm/models/blackoil/blackoilpolymerparams.hpp @@ -49,32 +49,6 @@ struct BlackOilPolymerParams { template void initFromState(const EclipseState& eclState); - /*! - * \brief Specify the number of satuation regions. - * - * This must be called before setting the PLYROCK and PLYADS of any region. - */ - void setNumSatRegions(unsigned numRegions); - - /*! - * \brief Specify the number of mix regions. - * - * This must be called before setting the PLYMAC and PLMIXPAR of any region. - */ - void setNumMixRegions(unsigned numRegions, bool enablePolymerMolarWeight); - - /*! - * \brief Specify the polymer rock properties a single region. - * - * The index of specified here must be in range [0, numSatRegions) - */ - void setPlyrock(unsigned satRegionIdx, - const Scalar& plyrockDeadPoreVolume, - const Scalar& plyrockResidualResistanceFactor, - const Scalar& plyrockRockDensityFactor, - const Scalar& plyrockAdsorbtionIndex, - const Scalar& plyrockMaxAdsorbtion); - // a struct containing the constants to calculate polymer viscosity // based on Mark-Houwink equation and Huggins equation, the constants are provided // by the keyword PLYVMH @@ -110,6 +84,91 @@ struct BlackOilPolymerParams { std::map skprwatTables_{}; std::map skprpolyTables_{}; + +private: + /*! + * \brief Specify the number of saturation regions. + * + * This must be called before setting the PLYROCK and PLYADS of any region. + */ + void setNumSatRegions(unsigned numRegions); + + /*! + * \brief Specify the number of mix regions. + * + * This must be called before setting the PLYMAX and PLMIXPAR of any region. + */ + void setNumMixRegions(unsigned numRegions, bool enablePolymerMolarWeight); + + /*! + * \brief Specify the polymer rock properties a single region. + * + * The index of specified here must be in range [0, numSatRegions) + */ + void setPlyrock(unsigned satRegionIdx, + const Scalar& plyrockDeadPoreVolume, + const Scalar& plyrockResidualResistanceFactor, + const Scalar& plyrockRockDensityFactor, + const Scalar& plyrockAdsorbtionIndex, + const Scalar& plyrockMaxAdsorbtion); + + /*! + * \brief Process the Plyrock data. + */ + void processPlyrock(const EclipseState& eclState); + + /*! + * \brief Process the Plyads data. + */ + void processPlyads(const EclipseState& eclState); + + /*! + * \brief Process the Plyvisc data. + */ + template + void processPlyvisc(const EclipseState& eclState); + + /*! + * \brief Process the Plymax data. + */ + template + void processPlymax(const EclipseState& eclState); + + /*! + * \brief Process the Plmixpar data. + */ + template + void processPlmixpar(const EclipseState& eclState); + + /*! + * \brief Process the Plyshlog data. + */ + void processPlyshlog(const EclipseState& eclState); + + /*! + * \brief Process the Shrate data. + */ + void processShrate(const EclipseState& eclState); + + /*! + * \brief Process the Plyvmh data. + */ + void processPlyvmh(const EclipseState& eclState); + + /*! + * \brief Process the Plymwinj data. + */ + void processPlymwinj(const EclipseState& eclState); + + /*! + * \brief Process the Skprwat data. + */ + void processSkprwat(const EclipseState& eclState); + + /*! + * \brief Process the Skprpoly data. + */ + void processSkprpoly(const EclipseState& eclState); }; } // namespace Opm diff --git a/opm/models/blackoil/blackoilsolventparams.cpp b/opm/models/blackoil/blackoilsolventparams.cpp index 465a920df86..cc6d9e258b8 100644 --- a/opm/models/blackoil/blackoilsolventparams.cpp +++ b/opm/models/blackoil/blackoilsolventparams.cpp @@ -34,31 +34,83 @@ #include #include -namespace Opm { +namespace { -template template -void BlackOilSolventParams:: -initFromState(const EclipseState& eclState, const Schedule& schedule) +void verifyState(const Opm::EclipseState& eclState) { // some sanity checks: if solvents are enabled, the SOLVENT keyword must be // present, if solvents are disabled the keyword must not be present. if constexpr (enableSolvent) { - if (!eclState.runspec().phases().active(Phase::SOLVENT)) { + if (!eclState.runspec().phases().active(Opm::Phase::SOLVENT)) { throw std::runtime_error("Non-trivial solvent treatment requested at compile " "time, but the deck does not contain the SOLVENT keyword"); } } else { - if (eclState.runspec().phases().active(Phase::SOLVENT)) { + if (eclState.runspec().phases().active(Opm::Phase::SOLVENT)) { throw std::runtime_error("Solvent treatment disabled at compile time, but the deck " "contains the SOLVENT keyword"); } } +} + +} + +namespace Opm { +template +template +void BlackOilSolventParams:: +initFromState(const EclipseState& eclState, const Schedule& schedule) +{ + verifyState(eclState); if (!eclState.runspec().phases().active(Phase::SOLVENT)) { return; // solvent treatment is supposed to be disabled } + setupPvts(eclState, schedule); + processSsfn(eclState); + + const auto& tableManager = eclState.getTableManager(); + + // initialize the objects needed for miscible solvent and oil simulations + isMiscible_ = false; + if (!tableManager.getMiscTables().empty()) { + isMiscible_ = true; + processSof2(eclState); + processMisc(eclState); + processPmisc(eclState); + processMsfn(eclState); + processSorwmis(eclState); + processSgcwmis(eclState); + processTlmixpar(eclState); + processTlpmixpa(eclState); + } +} + +template +void BlackOilSolventParams:: +setNumSatRegions(unsigned numRegions) +{ + ssfnKrg_.resize(numRegions); + ssfnKrs_.resize(numRegions); +} + +template +void BlackOilSolventParams:: +setMsfn(unsigned satRegionIdx, + const TabulatedFunction& msfnKrsg, + const TabulatedFunction& msfnKro) +{ + msfnKrsg_[satRegionIdx] = msfnKrsg; + msfnKro_[satRegionIdx] = msfnKro; +} + +template +void BlackOilSolventParams:: +setupPvts(const EclipseState& eclState, + const Schedule& schedule) +{ co2sol_ = eclState.runspec().co2Sol(); h2sol_ = eclState.runspec().h2Sol(); @@ -70,20 +122,28 @@ initFromState(const EclipseState& eclState, const Schedule& schedule) if (co2sol_) { co2GasPvt_.initFromState(eclState, schedule); brineCo2Pvt_.initFromState(eclState, schedule); - } else { + } + else { h2GasPvt_.initFromState(eclState, schedule); brineH2Pvt_.initFromState(eclState, schedule); } if (eclState.getSimulationConfig().hasDISGASW()) { rsSolw_active_ = true; } - } else + } + else { solventPvt_.initFromState(eclState, schedule); + } +} +template +void BlackOilSolventParams:: +processSsfn(const EclipseState& eclState) +{ const auto& tableManager = eclState.getTableManager(); // initialize the objects which deal with the SSFN keyword const auto& ssfnTables = tableManager.getSsfnTables(); - unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + const unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); setNumSatRegions(numSatRegions); for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { const auto& ssfnTable = ssfnTables.template getTable(satRegionIdx); @@ -94,226 +154,254 @@ initFromState(const EclipseState& eclState, const Schedule& schedule) ssfnTable.getSolventRelPermMultiplierColumn(), /*sortInput=*/true); } +} - // initialize the objects needed for miscible solvent and oil simulations - isMiscible_ = false; - if (!eclState.getTableManager().getMiscTables().empty()) { - isMiscible_ = true; - - unsigned numMiscRegions = 1; - - // misicible hydrocabon relative permeability wrt water - const auto& sof2Tables = tableManager.getSof2Tables(); - if (!sof2Tables.empty()) { - // resize the attributes of the object - sof2Krn_.resize(numSatRegions); - for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { - const auto& sof2Table = sof2Tables.template getTable(satRegionIdx); - sof2Krn_[satRegionIdx].setXYContainers(sof2Table.getSoColumn(), - sof2Table.getKroColumn(), - /*sortInput=*/true); - } - } - else if (eclState.runspec().phases().active(Phase::OIL)) { - throw std::runtime_error("SOF2 must be specified in MISCIBLE (SOLVENT and OIL) runs\n"); +template +void BlackOilSolventParams:: +processSof2(const EclipseState& eclState) +{ + // miscible hydrocarbon relative permeability wrt water + const auto& tableManager = eclState.getTableManager(); + const auto& sof2Tables = tableManager.getSof2Tables(); + const unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + if (!sof2Tables.empty()) { + // resize the attributes of the object + sof2Krn_.resize(numSatRegions); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { + const auto& sof2Table = sof2Tables.template getTable(satRegionIdx); + sof2Krn_[satRegionIdx].setXYContainers(sof2Table.getSoColumn(), + sof2Table.getKroColumn(), + /*sortInput=*/true); } + } + else if (eclState.runspec().phases().active(Phase::OIL)) { + throw std::runtime_error("SOF2 must be specified in MISCIBLE (SOLVENT and OIL) runs\n"); + } +} - const auto& miscTables = tableManager.getMiscTables(); - if (!miscTables.empty()) { - assert(numMiscRegions == miscTables.size()); - - // resize the attributes of the object - misc_.resize(numMiscRegions); - for (unsigned miscRegionIdx = 0; miscRegionIdx < numMiscRegions; ++miscRegionIdx) { - const auto& miscTable = miscTables.template getTable(miscRegionIdx); +template +void BlackOilSolventParams:: +processMisc(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const auto& miscTables = tableManager.getMiscTables(); + const unsigned numMiscRegions = 1; + if (!miscTables.empty()) { + assert(numMiscRegions == miscTables.size()); - // solventFraction = Ss / (Ss + Sg); - const auto& solventFraction = miscTable.getSolventFractionColumn(); - const auto& misc = miscTable.getMiscibilityColumn(); - misc_[miscRegionIdx].setXYContainers(solventFraction, misc); - } - } - else { - throw std::runtime_error("MISC must be specified in MISCIBLE (SOLVENT) runs\n"); + // resize the attributes of the object + misc_.resize(numMiscRegions); + for (unsigned miscRegionIdx = 0; miscRegionIdx < numMiscRegions; ++miscRegionIdx) { + const auto& miscTable = miscTables.template getTable(miscRegionIdx); + + // solventFraction = Ss / (Ss + Sg); + const auto& solventFraction = miscTable.getSolventFractionColumn(); + const auto& misc = miscTable.getMiscibilityColumn(); + misc_[miscRegionIdx].setXYContainers(solventFraction, misc); } + } + else { + throw std::runtime_error("MISC must be specified in MISCIBLE (SOLVENT) runs\n"); + } +} - // resize the attributes of the object - pmisc_.resize(numMiscRegions); - const auto& pmiscTables = tableManager.getPmiscTables(); - if (!pmiscTables.empty()) { - assert(numMiscRegions == pmiscTables.size()); +template +void BlackOilSolventParams:: +processPmisc(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const unsigned numMiscRegions = 1; - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { - const auto& pmiscTable = pmiscTables.template getTable(regionIdx); + // resize the attributes of the object + pmisc_.resize(numMiscRegions); + const auto& pmiscTables = tableManager.getPmiscTables(); + if (!pmiscTables.empty()) { + assert(numMiscRegions == pmiscTables.size()); - // Copy data - const auto& po = pmiscTable.getOilPhasePressureColumn(); - const auto& pmisc = pmiscTable.getMiscibilityColumn(); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + const auto& pmiscTable = pmiscTables.template getTable(regionIdx); - pmisc_[regionIdx].setXYContainers(po, pmisc); - } + // Copy data + const auto& po = pmiscTable.getOilPhasePressureColumn(); + const auto& pmisc = pmiscTable.getMiscibilityColumn(); + + pmisc_[regionIdx].setXYContainers(po, pmisc); } - else { - std::vector x = {0.0,1.0e20}; - std::vector y = {1.0,1.0}; - TabulatedFunction constant = TabulatedFunction(2, x, y); - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { - pmisc_[regionIdx] = constant; - } + } + else { + const std::vector x = {0.0,1.0e20}; + const std::vector y = {1.0,1.0}; + const TabulatedFunction constant = TabulatedFunction(2, x, y); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + pmisc_[regionIdx] = constant; } + } +} - // miscible relative permeability multipleiers - msfnKrsg_.resize(numSatRegions); - msfnKro_.resize(numSatRegions); - const auto& msfnTables = tableManager.getMsfnTables(); - if (!msfnTables.empty()) { - assert(numSatRegions == msfnTables.size()); - - for (unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) { - const MsfnTable& msfnTable = msfnTables.template getTable(regionIdx); - - // Copy data - // Ssg = Ss + Sg; - const auto& Ssg = msfnTable.getGasPhaseFractionColumn(); - const auto& krsg = msfnTable.getGasSolventRelpermMultiplierColumn(); - const auto& kro = msfnTable.getOilRelpermMultiplierColumn(); - - msfnKrsg_[regionIdx].setXYContainers(Ssg, krsg); - msfnKro_[regionIdx].setXYContainers(Ssg, kro); - } +template +void BlackOilSolventParams:: +processMsfn(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + + // miscible relative permeability multipleiers + msfnKrsg_.resize(numSatRegions); + msfnKro_.resize(numSatRegions); + const auto& msfnTables = tableManager.getMsfnTables(); + if (!msfnTables.empty()) { + assert(numSatRegions == msfnTables.size()); + + for (unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) { + const MsfnTable& msfnTable = msfnTables.template getTable(regionIdx); + + // Copy data + // Ssg = Ss + Sg; + const auto& Ssg = msfnTable.getGasPhaseFractionColumn(); + const auto& krsg = msfnTable.getGasSolventRelpermMultiplierColumn(); + const auto& kro = msfnTable.getOilRelpermMultiplierColumn(); + + msfnKrsg_[regionIdx].setXYContainers(Ssg, krsg); + msfnKro_[regionIdx].setXYContainers(Ssg, kro); } - else { - std::vector x = {0.0,1.0}; - std::vector y = {1.0,0.0}; - TabulatedFunction unit = TabulatedFunction(2, x, x); - TabulatedFunction invUnit = TabulatedFunction(2, x, y); - - for (unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) { - setMsfn(regionIdx, unit, invUnit); - } + } + else { + const std::vector x = {0.0,1.0}; + const std::vector y = {1.0,0.0}; + const TabulatedFunction unit = TabulatedFunction(2, x, x); + const TabulatedFunction invUnit = TabulatedFunction(2, x, y); + + for (unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) { + setMsfn(regionIdx, unit, invUnit); } - // resize the attributes of the object - sorwmis_.resize(numMiscRegions); - const auto& sorwmisTables = tableManager.getSorwmisTables(); - if (!sorwmisTables.empty()) { - assert(numMiscRegions == sorwmisTables.size()); + } +} - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { - const auto& sorwmisTable = sorwmisTables.template getTable(regionIdx); +template +void BlackOilSolventParams:: +processSorwmis(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const unsigned numMiscRegions = 1; + sorwmis_.resize(numMiscRegions); + const auto& sorwmisTables = tableManager.getSorwmisTables(); + if (!sorwmisTables.empty()) { + assert(numMiscRegions == sorwmisTables.size()); + + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + const auto& sorwmisTable = sorwmisTables.template getTable(regionIdx); - // Copy data - const auto& sw = sorwmisTable.getWaterSaturationColumn(); - const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn(); + // Copy data + const auto& sw = sorwmisTable.getWaterSaturationColumn(); + const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn(); - sorwmis_[regionIdx].setXYContainers(sw, sorwmis); - } + sorwmis_[regionIdx].setXYContainers(sw, sorwmis); } - else { - // default - std::vector x = {0.0,1.0}; - std::vector y = {0.0,0.0}; - TabulatedFunction zero = TabulatedFunction(2, x, y); - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { - sorwmis_[regionIdx] = zero; - } + } + else { + // default + const std::vector x = {0.0,1.0}; + const std::vector y = {0.0,0.0}; + const TabulatedFunction zero = TabulatedFunction(2, x, y); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + sorwmis_[regionIdx] = zero; } + } +} - // resize the attributes of the object - sgcwmis_.resize(numMiscRegions); - const auto& sgcwmisTables = tableManager.getSgcwmisTables(); - if (!sgcwmisTables.empty()) { - assert(numMiscRegions == sgcwmisTables.size()); - - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { - const auto& sgcwmisTable = sgcwmisTables.template getTable(regionIdx); +template +void BlackOilSolventParams:: +processSgcwmis(const EclipseState& eclState) +{ + const auto& tableManager = eclState.getTableManager(); + const unsigned numMiscRegions = 1; - // Copy data - const auto& sw = sgcwmisTable.getWaterSaturationColumn(); - const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn(); + // resize the attributes of the object + sgcwmis_.resize(numMiscRegions); + const auto& sgcwmisTables = tableManager.getSgcwmisTables(); + if (!sgcwmisTables.empty()) { + assert(numMiscRegions == sgcwmisTables.size()); - sgcwmis_[regionIdx].setXYContainers(sw, sgcwmis); - } - } - else { - // default - std::vector x = {0.0,1.0}; - std::vector y = {0.0,0.0}; - TabulatedFunction zero = TabulatedFunction(2, x, y); - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) - sgcwmis_[regionIdx] = zero; - } + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + const auto& sgcwmisTable = sgcwmisTables.template getTable(regionIdx); - const auto& tlmixpar = eclState.getTableManager().getTLMixpar(); - if (!tlmixpar.empty()) { - // resize the attributes of the object - tlMixParamViscosity_.resize(numMiscRegions); - tlMixParamDensity_.resize(numMiscRegions); - - assert(numMiscRegions == tlmixpar.size()); - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { - const auto& tlp = tlmixpar[regionIdx]; - tlMixParamViscosity_[regionIdx] = tlp.viscosity_parameter; - tlMixParamDensity_[regionIdx] = tlp.density_parameter; - } - } - else - throw std::runtime_error("TLMIXPAR must be specified in MISCIBLE (SOLVENT) runs\n"); + // Copy data + const auto& sw = sgcwmisTable.getWaterSaturationColumn(); + const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn(); - // resize the attributes of the object - tlPMixTable_.resize(numMiscRegions); - if (!eclState.getTableManager().getTlpmixpaTables().empty()) { - const auto& tlpmixparTables = tableManager.getTlpmixpaTables(); - if (!tlpmixparTables.empty()) { - assert(numMiscRegions == tlpmixparTables.size()); - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { - const auto& tlpmixparTable = tlpmixparTables.template getTable(regionIdx); - - // Copy data - const auto& po = tlpmixparTable.getOilPhasePressureColumn(); - const auto& tlpmixpa = tlpmixparTable.getMiscibilityColumn(); - - tlPMixTable_[regionIdx].setXYContainers(po, tlpmixpa); - } - } - else { - // if empty keyword. Try to use the pmisc table as default. - if (!pmisc_.empty()) { - tlPMixTable_ = pmisc_; - } - else { - throw std::invalid_argument("If the pressure dependent TL values in " - "TLPMIXPA is defaulted (no entries), then " - "the PMISC tables must be specified."); - } - } + sgcwmis_[regionIdx].setXYContainers(sw, sgcwmis); } - else { - // default - std::vector x = {0.0,1.0e20}; - std::vector y = {1.0,1.0}; - TabulatedFunction ones = TabulatedFunction(2, x, y); - for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) - tlPMixTable_[regionIdx] = ones; + } + else { + // default + const std::vector x = {0.0,1.0}; + const std::vector y = {0.0,0.0}; + const TabulatedFunction zero = TabulatedFunction(2, x, y); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + sgcwmis_[regionIdx] = zero; } } } template void BlackOilSolventParams:: -setNumSatRegions(unsigned numRegions) +processTlmixpar(const EclipseState& eclState) { - ssfnKrg_.resize(numRegions); - ssfnKrs_.resize(numRegions); + const auto& tableManager = eclState.getTableManager(); + const unsigned numMiscRegions = 1; + const auto& tlmixpar = tableManager.getTLMixpar(); + if (!tlmixpar.empty()) { + // resize the attributes of the object + tlMixParamViscosity_.resize(numMiscRegions); + tlMixParamDensity_.resize(numMiscRegions); + + assert(numMiscRegions == tlmixpar.size()); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + const auto& tlp = tlmixpar[regionIdx]; + tlMixParamViscosity_[regionIdx] = tlp.viscosity_parameter; + tlMixParamDensity_[regionIdx] = tlp.density_parameter; + } + } + else { + throw std::runtime_error("TLMIXPAR must be specified in MISCIBLE (SOLVENT) runs\n"); + } } template void BlackOilSolventParams:: -setMsfn(unsigned satRegionIdx, - const TabulatedFunction& msfnKrsg, - const TabulatedFunction& msfnKro) +processTlpmixpa(const EclipseState& eclState) { - msfnKrsg_[satRegionIdx] = msfnKrsg; - msfnKro_[satRegionIdx] = msfnKro; + const auto& tableManager = eclState.getTableManager(); + const unsigned numMiscRegions = 1; + const auto& tlpmixparTables = tableManager.getTlpmixpaTables(); + + // resize the attributes of the object + tlPMixTable_.resize(numMiscRegions); + if (!tableManager.getTlpmixpaTables().empty()) { + assert(numMiscRegions == tlpmixparTables.size()); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + const auto& tlpmixparTable = tlpmixparTables.template getTable(regionIdx); + + // Copy data + const auto& po = tlpmixparTable.getOilPhasePressureColumn(); + const auto& tlpmixpa = tlpmixparTable.getMiscibilityColumn(); + + tlPMixTable_[regionIdx].setXYContainers(po, tlpmixpa); + } + } + else { + OpmLog::warning("The pressure dependent TL values in " + "TLPMIXPA is missing or defaulted," + "using dummy values"); + + // default + const std::vector x = {0.0,1.0e20}; + const std::vector y = {1.0,1.0}; + const TabulatedFunction ones = TabulatedFunction(2, x, y); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + tlPMixTable_[regionIdx] = ones; + } + } } #define INSTANTIATE_TYPE(T) \ diff --git a/opm/models/blackoil/blackoilsolventparams.hpp b/opm/models/blackoil/blackoilsolventparams.hpp index 3d2ba1c397a..3562f691513 100644 --- a/opm/models/blackoil/blackoilsolventparams.hpp +++ b/opm/models/blackoil/blackoilsolventparams.hpp @@ -80,8 +80,9 @@ struct BlackOilSolventParams bool co2sol_ = false; bool h2sol_ = false; +private: /*! - * \brief Specify the number of satuation regions. + * \brief Specify the number of saturation regions. * * This must be called before setting the SSFN of any region. */ @@ -95,6 +96,57 @@ struct BlackOilSolventParams void setMsfn(unsigned satRegionIdx, const TabulatedFunction& msfnKrsg, const TabulatedFunction& msfnKro); + + /*! + * \brief Setup active pvt members. + */ + void setupPvts(const EclipseState& eclState, + const Schedule& schedule); + + /*! + * \brief Process the SSFN data. + */ + void processSsfn(const EclipseState& eclState); + + /*! + * \brief Process the SOF2 data. + */ + void processSof2(const EclipseState& eclState); + + /*! + * \brief Process the MISC data. + */ + void processMisc(const EclipseState& eclState); + + /*! + * \brief Process the PMISC data. + */ + void processPmisc(const EclipseState& eclState); + + /*! + * \brief Process the MSFN data. + */ + void processMsfn(const EclipseState& eclState); + + /*! + * \brief Process the SORWMIS data. + */ + void processSorwmis(const EclipseState& eclState); + + /*! + * \brief Process the SGCWMIS data. + */ + void processSgcwmis(const EclipseState& eclState); + + /*! + * \brief Process the TLMIXPAR data. + */ + void processTlmixpar(const EclipseState& eclState); + + /*! + * \brief Process the TLPMIXPA data. + */ + void processTlpmixpa(const EclipseState& eclState); }; } // namespace Opm