From a15481dc0f78db51e76b1276dc4c4f5d88db593d Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Sep 2022 11:41:39 +0200 Subject: [PATCH 1/6] import code for initializing brine module from eclipsestate this allows us to keep this code in a separate compile unit, limiting the amount of times it has to be built. --- CMakeLists_files.cmake | 1 + ebos/eclblackoilmoduleinit.cc | 106 ++++++++++++++++++++++++++++++++++ ebos/eclblackoilmoduleinit.hh | 36 ++++++++++++ ebos/eclproblem.hh | 4 +- 4 files changed, 146 insertions(+), 1 deletion(-) create mode 100644 ebos/eclblackoilmoduleinit.cc create mode 100644 ebos/eclblackoilmoduleinit.hh diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 347a1e9a91d..284026a1d49 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -25,6 +25,7 @@ list (APPEND MAIN_SOURCE_FILES ebos/collecttoiorank.cc ebos/eclactionhandler.cc + ebos/eclblackoilmoduleinit.cc ebos/eclgenericcpgridvanguard.cc ebos/eclgenericoutputblackoilmodule.cc ebos/eclgenericproblem.cc diff --git a/ebos/eclblackoilmoduleinit.cc b/ebos/eclblackoilmoduleinit.cc new file mode 100644 index 00000000000..ce685a91b8e --- /dev/null +++ b/ebos/eclblackoilmoduleinit.cc @@ -0,0 +1,106 @@ +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace Opm { + +template +BlackOilBrineParams +setupBrineParams(bool enableBrine, + const EclipseState& eclState) +{ + BlackOilBrineParams params; + + // some sanity checks: if brine are enabled, the BRINE keyword must be + // present, if brine are disabled the keyword must not be present. + if (enableBrine && !eclState.runspec().phases().active(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 (!enableBrine && eclState.runspec().phases().active(Phase::BRINE)) { + throw std::runtime_error("Brine treatment disabled at compile time, but the deck " + "contains the BRINE keyword"); + } + + if (!eclState.runspec().phases().active(Phase::BRINE)) + return params; // brine treatment is supposed to be disabled + + const auto& tableManager = eclState.getTableManager(); + + unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + params.referencePressure_.resize(numPvtRegions); + + const auto& pvtwsaltTables = tableManager.getPvtwSaltTables(); + + // initialize the objects which deal with the BDENSITY keyword + const auto& bdensityTables = tableManager.getBrineDensityTables(); + if (!bdensityTables.empty()) { + params.bdensityTable_.resize(numPvtRegions); + assert(numPvtRegions == bdensityTables.size()); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) { + const auto& bdensityTable = bdensityTables[pvtRegionIdx]; + const auto& pvtwsaltTable = pvtwsaltTables[pvtRegionIdx]; + const auto& c = pvtwsaltTable.getSaltConcentrationColumn(); + params.bdensityTable_[pvtRegionIdx].setXYContainers(c, bdensityTable); + } + } + + if constexpr (enableSaltPrecipitation) { + const TableContainer& permfactTables = tableManager.getPermfactTables(); + params.permfactTable_.resize(numPvtRegions); + for (size_t i = 0; i < permfactTables.size(); ++i) { + const PermfactTable& permfactTable = permfactTables.getTable(i); + params.permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn()); + } + + const TableContainer& saltsolTables = tableManager.getSaltsolTables(); + if (!saltsolTables.empty()) { + params.saltsolTable_.resize(numPvtRegions); + assert(numPvtRegions == saltsolTables.size()); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) { + const SaltsolTable& saltsolTable = saltsolTables.getTable(pvtRegionIdx ); + params.saltsolTable_[pvtRegionIdx] = saltsolTable.getSaltsolColumn().front(); + } + } + } + + return params; +} + +template BlackOilBrineParams +setupBrineParams(bool, const EclipseState&); +template BlackOilBrineParams +setupBrineParams(bool, const EclipseState&); + +} diff --git a/ebos/eclblackoilmoduleinit.hh b/ebos/eclblackoilmoduleinit.hh new file mode 100644 index 00000000000..08c8b063c84 --- /dev/null +++ b/ebos/eclblackoilmoduleinit.hh @@ -0,0 +1,36 @@ +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ + +#ifndef ECL_BLACKOILMODULE_INIT_HH +#define ECL_BLACKOILMODULE_INIT_HH + +namespace Opm { + +class EclipseState; +template struct BlackOilBrineParams; + +//! \brief Setup parameters for brine module from an EclipseState. +template +BlackOilBrineParams setupBrineParams(bool enableBrine, + const EclipseState& eclState); +} + +#endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 1400f97d3ad..794182bd5d4 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -49,6 +49,7 @@ #endif #include "eclactionhandler.hh" +#include "eclblackoilmoduleinit.hh" #include "eclequilinitializer.hh" #include "eclwriter.hh" #include "ecloutputblackoilmodule.hh" @@ -783,10 +784,11 @@ public: this->model().addOutputModule(new VtkEclTracerModule(simulator)); // Tell the black-oil extensions to initialize their internal data structures const auto& vanguard = simulator.vanguard(); + BrineModule::setParams(setupBrineParams(enableBrine, + vanguard.eclState())); SolventModule::initFromState(vanguard.eclState(), vanguard.schedule()); PolymerModule::initFromState(vanguard.eclState()); FoamModule::initFromState(vanguard.eclState()); - BrineModule::initFromState(vanguard.eclState()); ExtboModule::initFromState(vanguard.eclState()); MICPModule::initFromState(vanguard.eclState()); From 0acddf5c674d973c1408d272861b1bdab5a35d8c Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Sep 2022 11:41:39 +0200 Subject: [PATCH 2/6] import code for initializing extbo module from eclipsestate this allows us to keep this code in a separate compile unit, limiting the amount of times it has to be built. --- ebos/eclblackoilmoduleinit.cc | 175 +++++++++++++++++++++++++++++++++- ebos/eclblackoilmoduleinit.hh | 6 ++ ebos/eclproblem.hh | 2 +- 3 files changed, 181 insertions(+), 2 deletions(-) diff --git a/ebos/eclblackoilmoduleinit.cc b/ebos/eclblackoilmoduleinit.cc index ce685a91b8e..58c466b35fd 100644 --- a/ebos/eclblackoilmoduleinit.cc +++ b/ebos/eclblackoilmoduleinit.cc @@ -23,14 +23,24 @@ #include #include +#include +#include #include #include +#include #include -#include +#include #include +#include +#include +#include +#include +#include #include +#include +#include #include namespace Opm { @@ -98,9 +108,172 @@ setupBrineParams(bool enableBrine, return params; } +template +BlackOilExtboParams setupExtboParams(bool enableExtbo, + const EclipseState& eclState) +{ + BlackOilExtboParams params; + // 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 (enableExtbo && !eclState.runspec().phases().active(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)) + throw std::runtime_error("Extended black oil treatment disabled at compile time, but the deck " + "contains the PVTSOL keyword"); + + if (!eclState.runspec().phases().active(Phase::ZFRACTION)) + return params; // solvent treatment is supposed to be disabled + + // pvt properties from kw PVTSOL: + + const auto& tableManager = eclState.getTableManager(); + const auto& pvtsolTables = tableManager.getPvtsolTables(); + + size_t numPvtRegions = pvtsolTables.size(); + + using Tabulated2DFunction = typename BlackOilExtboParams::Tabulated2DFunction; + + params.BO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.BG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.X_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.Y_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.VISCO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.VISCG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + + params.PBUB_RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + params.PBUB_RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + + params.zLim_.resize(numPvtRegions); + + const bool extractCmpFromPvt = true; //: Default values used in [*] + params.oilCmp_.resize(numPvtRegions); + params.gasCmp_.resize(numPvtRegions); + + for (unsigned regionIdx = 0; regionIdx < numPvtRegions; ++ regionIdx) { + const auto& pvtsolTable = pvtsolTables[regionIdx]; + + const auto& saturatedTable = pvtsolTable.getSaturatedTable(); + assert(saturatedTable.numRows() > 1); + + std::vector oilCmp(saturatedTable.numRows(), -4.0e-9); //Default values used in [*] + std::vector gasCmp(saturatedTable.numRows(), -0.08); //-------------"------------- + params.zLim_[regionIdx] = 0.7; //-------------"------------- + std::vector zArg(saturatedTable.numRows(), 0.0); + + for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) { + Scalar ZCO2 = saturatedTable.get("ZCO2", outerIdx); + + zArg[outerIdx] = ZCO2; + + params.BO_[regionIdx].appendXPos(ZCO2); + params.BG_[regionIdx].appendXPos(ZCO2); + + params.RS_[regionIdx].appendXPos(ZCO2); + params.RV_[regionIdx].appendXPos(ZCO2); + + params.X_[regionIdx].appendXPos(ZCO2); + params.Y_[regionIdx].appendXPos(ZCO2); + + params.VISCO_[regionIdx].appendXPos(ZCO2); + params.VISCG_[regionIdx].appendXPos(ZCO2); + + params.PBUB_RS_[regionIdx].appendXPos(ZCO2); + params.PBUB_RV_[regionIdx].appendXPos(ZCO2); + + const auto& underSaturatedTable = pvtsolTable.getUnderSaturatedTable(outerIdx); + size_t numRows = underSaturatedTable.numRows(); + + Scalar bo0 = 0.0; + Scalar po0 = 0.0; + for (unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) { + Scalar po = underSaturatedTable.get("P", innerIdx); + Scalar bo = underSaturatedTable.get("B_O", innerIdx); + Scalar bg = underSaturatedTable.get("B_G", innerIdx); + Scalar rs = underSaturatedTable.get("RS", innerIdx)+innerIdx*1.0e-10; + Scalar rv = underSaturatedTable.get("RV", innerIdx)+innerIdx*1.0e-10; + Scalar xv = underSaturatedTable.get("XVOL", innerIdx); + Scalar yv = underSaturatedTable.get("YVOL", innerIdx); + Scalar mo = underSaturatedTable.get("MU_O", innerIdx); + Scalar mg = underSaturatedTable.get("MU_G", innerIdx); + + 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); + oilCmp[outerIdx] = cmpFactor; + params.zLim_[regionIdx] = ZCO2; + //std::cout << "### cmpFactorOil: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl; + } + break; + } 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); + gasCmp[outerIdx] = cmpFactor; + //std::cout << "### cmpFactorGas: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl; + } + + params.BO_[regionIdx].appendSamplePoint(outerIdx,po,bo); + params.BG_[regionIdx].appendSamplePoint(outerIdx,po,bg); + params.RS_[regionIdx].appendSamplePoint(outerIdx,po,rs); + params.RV_[regionIdx].appendSamplePoint(outerIdx,po,rv); + params.X_[regionIdx].appendSamplePoint(outerIdx,po,xv); + params.Y_[regionIdx].appendSamplePoint(outerIdx,po,yv); + params.VISCO_[regionIdx].appendSamplePoint(outerIdx,po,mo); + params.VISCG_[regionIdx].appendSamplePoint(outerIdx,po,mg); + break; + } + + bo0 = bo; + po0 = po; + + params.BO_[regionIdx].appendSamplePoint(outerIdx,po,bo); + params.BG_[regionIdx].appendSamplePoint(outerIdx,po,bg); + + params.RS_[regionIdx].appendSamplePoint(outerIdx,po,rs); + params.RV_[regionIdx].appendSamplePoint(outerIdx,po,rv); + + params.X_[regionIdx].appendSamplePoint(outerIdx,po,xv); + params.Y_[regionIdx].appendSamplePoint(outerIdx,po,yv); + + params.VISCO_[regionIdx].appendSamplePoint(outerIdx,po,mo); + params.VISCG_[regionIdx].appendSamplePoint(outerIdx,po,mg); + + // rs,rv -> pressure + params.PBUB_RS_[regionIdx].appendSamplePoint(outerIdx, rs, po); + params.PBUB_RV_[regionIdx].appendSamplePoint(outerIdx, rv, po); + } + } + params.oilCmp_[regionIdx].setXYContainers(zArg, oilCmp, /*sortInput=*/false); + params.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) { + params.zReferenceDensity_.resize(numPvtRegions); + for (unsigned regionIdx = 0; regionIdx < numPvtRegions; ++ regionIdx) { + Scalar rhoRefS = sdensityTables[regionIdx].getSolventDensityColumn().front(); + params.zReferenceDensity_[regionIdx] = rhoRefS; + } + } + else + throw std::runtime_error("Extbo: kw SDENSITY is missing or not aligned with NTPVT\n"); + + return params; +} + template BlackOilBrineParams setupBrineParams(bool, const EclipseState&); template BlackOilBrineParams setupBrineParams(bool, const EclipseState&); +template BlackOilExtboParams +setupExtboParams(bool, const EclipseState&); + } diff --git a/ebos/eclblackoilmoduleinit.hh b/ebos/eclblackoilmoduleinit.hh index 08c8b063c84..10482761b29 100644 --- a/ebos/eclblackoilmoduleinit.hh +++ b/ebos/eclblackoilmoduleinit.hh @@ -26,11 +26,17 @@ namespace Opm { class EclipseState; template struct BlackOilBrineParams; +template struct BlackOilExtboParams; //! \brief Setup parameters for brine module from an EclipseState. template BlackOilBrineParams setupBrineParams(bool enableBrine, const EclipseState& eclState); + +//! \brief Setup parameters for extbo module from an EclipseState. +template +BlackOilExtboParams setupExtboParams(bool enableExtbo, + const EclipseState& eclState); } #endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 794182bd5d4..e91d75d0369 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -786,10 +786,10 @@ public: const auto& vanguard = simulator.vanguard(); BrineModule::setParams(setupBrineParams(enableBrine, vanguard.eclState())); + ExtboModule::setParams(setupExtboParams(enableExtbo, vanguard.eclState())); SolventModule::initFromState(vanguard.eclState(), vanguard.schedule()); PolymerModule::initFromState(vanguard.eclState()); FoamModule::initFromState(vanguard.eclState()); - ExtboModule::initFromState(vanguard.eclState()); MICPModule::initFromState(vanguard.eclState()); // create the ECL writer From 49cc85626eed27ed09129263f0be83e7ab538ab1 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Sep 2022 11:41:39 +0200 Subject: [PATCH 3/6] import code for initializing foam module from eclipsestate this allows us to keep this code in a separate compile unit, limiting the amount of times it has to be built. --- ebos/eclblackoilmoduleinit.cc | 97 +++++++++++++++++++++++++++++++++++ ebos/eclblackoilmoduleinit.hh | 7 +++ ebos/eclproblem.hh | 2 +- 3 files changed, 105 insertions(+), 1 deletion(-) diff --git a/ebos/eclblackoilmoduleinit.cc b/ebos/eclblackoilmoduleinit.cc index 58c466b35fd..c0341299ed4 100644 --- a/ebos/eclblackoilmoduleinit.cc +++ b/ebos/eclblackoilmoduleinit.cc @@ -23,6 +23,8 @@ #include #include +#include +#include #include #include #include @@ -39,6 +41,7 @@ #include #include +#include #include #include @@ -268,6 +271,97 @@ BlackOilExtboParams setupExtboParams(bool enableExtbo, return params; } +template +BlackOilFoamParams setupFoamParams(bool enableFoam, + const EclipseState& eclState) +{ + BlackOilFoamParams params; + // some sanity checks: if foam is enabled, the FOAM keyword must be + // present, if foam is disabled the keyword must not be present. + if (enableFoam && !eclState.runspec().phases().active(Phase::FOAM)) { + throw std::runtime_error("Non-trivial foam treatment requested at compile time, but " + "the deck does not contain the FOAM keyword"); + } + else if (!enableFoam && eclState.runspec().phases().active(Phase::FOAM)) { + throw std::runtime_error("Foam treatment disabled at compile time, but the deck " + "contains the FOAM keyword"); + } + + if (!eclState.runspec().phases().active(Phase::FOAM)) { + return params; // foam treatment is supposed to be disabled + } + + // Check that only implemented options are used. + // We only support the default values of FOAMOPTS (GAS, TAB). + if (eclState.getInitConfig().getFoamConfig().getTransportPhase() != Phase::GAS) { + throw std::runtime_error("In FOAMOPTS, only GAS is allowed for the transport phase."); + } + if (eclState.getInitConfig().getFoamConfig().getMobilityModel() != FoamConfig::MobilityModel::TAB) { + throw std::runtime_error("In FOAGMOPTS, only TAB is allowed for the gas mobility factor reduction model."); + } + + const auto& tableManager = eclState.getTableManager(); + const unsigned int numSatRegions = tableManager.getTabdims().getNumSatTables(); + params.setNumSatRegions(numSatRegions); + const unsigned int numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + params.gasMobilityMultiplierTable_.resize(numPvtRegions); + + // Get and check FOAMROCK data. + const FoamConfig& foamConf = eclState.getInitConfig().getFoamConfig(); + if (numSatRegions != foamConf.size()) { + throw std::runtime_error("Inconsistent sizes, number of saturation regions differ from the number of elements " + "in FoamConfig, which typically corresponds to the number of records in FOAMROCK."); + } + + // Get and check FOAMADS data. + const auto& foamadsTables = tableManager.getFoamadsTables(); + if (foamadsTables.empty()) { + throw std::runtime_error("FOAMADS must be specified in FOAM runs"); + } + if (numSatRegions != foamadsTables.size()) { + throw std::runtime_error("Inconsistent sizes, number of saturation regions differ from the " + "number of FOAMADS tables."); + } + + // Set data that vary with saturation region. + for (std::size_t satReg = 0; satReg < numSatRegions; ++satReg) { + const auto& rec = foamConf.getRecord(satReg); + params.foamCoefficients_[satReg] = typename BlackOilFoamParams::FoamCoefficients(); + params.foamCoefficients_[satReg].fm_min = rec.minimumSurfactantConcentration(); + params.foamCoefficients_[satReg].fm_surf = rec.referenceSurfactantConcentration(); + params.foamCoefficients_[satReg].ep_surf = rec.exponent(); + params.foamRockDensity_[satReg] = rec.rockDensity(); + params.foamAllowDesorption_[satReg] = rec.allowDesorption(); + const auto& foamadsTable = foamadsTables.template getTable(satReg); + const auto& conc = foamadsTable.getFoamConcentrationColumn(); + const auto& ads = foamadsTable.getAdsorbedFoamColumn(); + params.adsorbedFoamTable_[satReg].setXYContainers(conc, ads); + } + + // Get and check FOAMMOB data. + const auto& foammobTables = tableManager.getFoammobTables(); + if (foammobTables.empty()) { + // When in the future adding support for the functional + // model, FOAMMOB will not be required anymore (functional + // family of keywords can be used instead, FOAMFSC etc.). + throw std::runtime_error("FOAMMOB must be specified in FOAM runs"); + } + if (numPvtRegions != foammobTables.size()) { + throw std::runtime_error("Inconsistent sizes, number of PVT regions differ from the " + "number of FOAMMOB tables."); + } + + // Set data that vary with PVT region. + for (std::size_t pvtReg = 0; pvtReg < numPvtRegions; ++pvtReg) { + const auto& foammobTable = foammobTables.template getTable(pvtReg); + const auto& conc = foammobTable.getFoamConcentrationColumn(); + const auto& mobMult = foammobTable.getMobilityMultiplierColumn(); + params.gasMobilityMultiplierTable_[pvtReg].setXYContainers(conc, mobMult); + } + + return params; +} + template BlackOilBrineParams setupBrineParams(bool, const EclipseState&); template BlackOilBrineParams @@ -276,4 +370,7 @@ setupBrineParams(bool, const EclipseState&); template BlackOilExtboParams setupExtboParams(bool, const EclipseState&); +template BlackOilFoamParams +setupFoamParams(bool, const EclipseState&); + } diff --git a/ebos/eclblackoilmoduleinit.hh b/ebos/eclblackoilmoduleinit.hh index 10482761b29..58da6ea0868 100644 --- a/ebos/eclblackoilmoduleinit.hh +++ b/ebos/eclblackoilmoduleinit.hh @@ -27,6 +27,7 @@ namespace Opm { class EclipseState; template struct BlackOilBrineParams; template struct BlackOilExtboParams; +template struct BlackOilFoamParams; //! \brief Setup parameters for brine module from an EclipseState. template @@ -37,6 +38,12 @@ BlackOilBrineParams setupBrineParams(bool enableBrine, template BlackOilExtboParams setupExtboParams(bool enableExtbo, const EclipseState& eclState); + +//! \brief Setup parameters for foam module from an EclipseState. +template +BlackOilFoamParams setupFoamParams(bool enableFoam, + const EclipseState& eclState); + } #endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index e91d75d0369..db6484257b3 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -787,9 +787,9 @@ public: BrineModule::setParams(setupBrineParams(enableBrine, vanguard.eclState())); ExtboModule::setParams(setupExtboParams(enableExtbo, vanguard.eclState())); + FoamModule::setParams(setupFoamParams(enableFoam, vanguard.eclState())); SolventModule::initFromState(vanguard.eclState(), vanguard.schedule()); PolymerModule::initFromState(vanguard.eclState()); - FoamModule::initFromState(vanguard.eclState()); MICPModule::initFromState(vanguard.eclState()); // create the ECL writer From 8bdeb026205a9fc1319beaa723767e981f9b6ceb Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Sep 2022 11:41:39 +0200 Subject: [PATCH 4/6] import code for initializing micp module from eclipsestate this allows us to keep this code in a separate compile unit, limiting the amount of times it has to be built. --- ebos/eclblackoilmoduleinit.cc | 48 +++++++++++++++++++++++++++++++++++ ebos/eclblackoilmoduleinit.hh | 6 +++++ ebos/eclproblem.hh | 2 +- 3 files changed, 55 insertions(+), 1 deletion(-) diff --git a/ebos/eclblackoilmoduleinit.cc b/ebos/eclblackoilmoduleinit.cc index c0341299ed4..cfb0c53e442 100644 --- a/ebos/eclblackoilmoduleinit.cc +++ b/ebos/eclblackoilmoduleinit.cc @@ -42,6 +42,7 @@ #include #include #include +#include #include #include @@ -362,6 +363,50 @@ BlackOilFoamParams setupFoamParams(bool enableFoam, return params; } +template +BlackOilMICPParams setupMICPParams(bool enableMICP, + const EclipseState& eclState) +{ + BlackOilMICPParams params; + // some sanity checks: if MICP is enabled, the MICP keyword must be + // present, if MICP is disabled the keyword must not be present. + if (enableMICP && !eclState.runspec().micp()) { + throw std::runtime_error("Non-trivial MICP treatment requested at compile time, but " + "the deck does not contain the MICP keyword"); + } + else if (!enableMICP && eclState.runspec().micp()) { + throw std::runtime_error("MICP treatment disabled at compile time, but the deck " + "contains the MICP keyword"); + } + + if (!eclState.runspec().micp()) + return params; // MICP treatment is supposed to be disabled*/ + + // initialize the objects which deal with the MICPpara keyword + const auto& mp = eclState.getMICPpara(); + params.densityBiofilm_ = mp.getDensityBiofilm(); + params.densityCalcite_ = mp.getDensityCalcite(); + params.detachmentRate_ = mp.getDetachmentRate(); + params.criticalPorosity_ = mp.getCriticalPorosity(); + params.fittingFactor_ = mp.getFittingFactor(); + params.halfVelocityOxygen_ = mp.getHalfVelocityOxygen(); + params.halfVelocityUrea_ = mp.getHalfVelocityUrea(); + params.maximumGrowthRate_ = mp.getMaximumGrowthRate(); + params.maximumUreaUtilization_ = mp.getMaximumUreaUtilization(); + params.microbialAttachmentRate_ = mp.getMicrobialAttachmentRate(); + params.microbialDeathRate_ = mp.getMicrobialDeathRate(); + params.minimumPermeability_ = mp.getMinimumPermeability(); + params.oxygenConsumptionFactor_ = mp.getOxygenConsumptionFactor(); + params.yieldGrowthCoefficient_ = mp.getYieldGrowthCoefficient(); + params.maximumOxygenConcentration_ = mp.getMaximumOxygenConcentration(); + params.maximumUreaConcentration_ = mp.getMaximumUreaConcentration(); + params.toleranceBeforeClogging_ = mp.getToleranceBeforeClogging(); + // obtain the porosity for the clamp in the blackoilnewtonmethod + params.phi_ = eclState.fieldProps().get_double("PORO"); + + return params; +} + template BlackOilBrineParams setupBrineParams(bool, const EclipseState&); template BlackOilBrineParams @@ -373,4 +418,7 @@ setupExtboParams(bool, const EclipseState&); template BlackOilFoamParams setupFoamParams(bool, const EclipseState&); +template BlackOilMICPParams +setupMICPParams(bool, const EclipseState&); + } diff --git a/ebos/eclblackoilmoduleinit.hh b/ebos/eclblackoilmoduleinit.hh index 58da6ea0868..b69cba7045c 100644 --- a/ebos/eclblackoilmoduleinit.hh +++ b/ebos/eclblackoilmoduleinit.hh @@ -28,6 +28,7 @@ class EclipseState; template struct BlackOilBrineParams; template struct BlackOilExtboParams; template struct BlackOilFoamParams; +template struct BlackOilMICPParams; //! \brief Setup parameters for brine module from an EclipseState. template @@ -44,6 +45,11 @@ template BlackOilFoamParams setupFoamParams(bool enableFoam, const EclipseState& eclState); +//! \brief Setup parameters for foam module from an EclipseState. +template +BlackOilMICPParams setupMICPParams(bool enableMICP, + const EclipseState& eclState); + } #endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index db6484257b3..39b027ccd1b 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -788,9 +788,9 @@ public: vanguard.eclState())); ExtboModule::setParams(setupExtboParams(enableExtbo, vanguard.eclState())); FoamModule::setParams(setupFoamParams(enableFoam, vanguard.eclState())); + MICPModule::setParams(setupMICPParams(enableMICP, vanguard.eclState())); SolventModule::initFromState(vanguard.eclState(), vanguard.schedule()); PolymerModule::initFromState(vanguard.eclState()); - MICPModule::initFromState(vanguard.eclState()); // create the ECL writer eclWriter_.reset(new EclWriterType(simulator)); From 64fe730126bba304a638ea64ab4aa8e3f3db9dea Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Sep 2022 11:41:39 +0200 Subject: [PATCH 5/6] import code for initializing polymer module from eclipsestate this allows us to keep this code in a separate compile unit, limiting the amount of times it has to be built. --- ebos/eclblackoilmoduleinit.cc | 268 ++++++++++++++++++++++++++++++++++ ebos/eclblackoilmoduleinit.hh | 6 + ebos/eclproblem.hh | 3 +- 3 files changed, 276 insertions(+), 1 deletion(-) diff --git a/ebos/eclblackoilmoduleinit.cc b/ebos/eclblackoilmoduleinit.cc index cfb0c53e442..bac1c784d19 100644 --- a/ebos/eclblackoilmoduleinit.cc +++ b/ebos/eclblackoilmoduleinit.cc @@ -22,6 +22,8 @@ #include #include +#include + #include #include #include @@ -29,6 +31,11 @@ #include #include #include +#include +#include +#include +#include +#include #include #include #include @@ -43,6 +50,7 @@ #include #include #include +#include #include #include @@ -407,6 +415,261 @@ BlackOilMICPParams setupMICPParams(bool enableMICP, return params; } +template +BlackOilPolymerParams setupPolymerParams(bool enablePolymer, + const EclipseState& eclState) +{ + BlackOilPolymerParams params; + // some sanity checks: if polymers are enabled, the POLYMER keyword must be + // present, if polymers are disabled the keyword must not be present. + if (enablePolymer && !eclState.runspec().phases().active(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 (!enablePolymer && eclState.runspec().phases().active(Phase::POLYMER)) { + throw std::runtime_error("Polymer treatment disabled at compile time, but the deck " + "contains the POLYMER keyword"); + } + + if (enablePolymerMolarWeight && !eclState.runspec().phases().active(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 (!enablePolymerMolarWeight && eclState.runspec().phases().active(Phase::POLYMW)) { + throw std::runtime_error("Polymer molecular weight tracking is disabled at compile time, but the deck " + "contains the POLYMW keyword"); + } + + if (enablePolymerMolarWeight && !enablePolymer) { + throw std::runtime_error("Polymer molecular weight tracking is enabled while polymer treatment " + "is disabled at compile time"); + } + + if (!eclState.runspec().phases().active(Phase::POLYMER)) + return params; // polymer treatment is supposed to be disabled + + const auto& tableManager = eclState.getTableManager(); + + unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + params.setNumSatRegions(numSatRegions); + + // initialize the objects which deal with the PLYROCK keyword + const auto& plyrockTables = tableManager.getPlyrockTables(); + if (!plyrockTables.empty()) { + assert(numSatRegions == plyrockTables.size()); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) { + const auto& plyrockTable = plyrockTables.template getTable(satRegionIdx); + params.setPlyrock(satRegionIdx, + plyrockTable.getDeadPoreVolumeColumn()[0], + plyrockTable.getResidualResistanceFactorColumn()[0], + plyrockTable.getRockDensityFactorColumn()[0], + static_cast::AdsorptionBehaviour>(plyrockTable.getAdsorbtionIndexColumn()[0]), + plyrockTable.getMaxAdsorbtionColumn()[0]); + } + } + else { + throw std::runtime_error("PLYROCK must be specified in POLYMER runs\n"); + } + + // initialize the objects which deal with the PLYADS keyword + const auto& plyadsTables = tableManager.getPlyadsTables(); + if (!plyadsTables.empty()) { + assert(numSatRegions == plyadsTables.size()); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) { + const auto& plyadsTable = plyadsTables.template getTable(satRegionIdx); + // Copy data + const auto& c = plyadsTable.getPolymerConcentrationColumn(); + const auto& ads = plyadsTable.getAdsorbedPolymerColumn(); + params.plyadsAdsorbedPolymer_[satRegionIdx].setXYContainers(c, ads); + } + } + else { + throw std::runtime_error("PLYADS must be specified in POLYMER runs\n"); + } + + + unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + params.plyviscViscosityMultiplierTable_.resize(numPvtRegions); + + // initialize the objects which deal with the PLYVISC keyword + const auto& plyviscTables = tableManager.getPlyviscTables(); + if (!plyviscTables.empty()) { + // different viscosity model is used for POLYMW + if (enablePolymerMolarWeight) { + OpmLog::warning("PLYVISC should not be used in POLYMW runs, " + "it will have no effect. A viscosity model based on PLYVMH is used instead.\n"); + } + else { + assert(numPvtRegions == plyviscTables.size()); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) { + const auto& plyadsTable = plyviscTables.template getTable(pvtRegionIdx); + // Copy data + const auto& c = plyadsTable.getPolymerConcentrationColumn(); + const auto& visc = plyadsTable.getViscosityMultiplierColumn(); + params.plyviscViscosityMultiplierTable_[pvtRegionIdx].setXYContainers(c, visc); + } + } + } + else if (!enablePolymerMolarWeight) { + throw std::runtime_error("PLYVISC must be specified in POLYMER runs\n"); + } + + // initialize the objects which deal with the PLYMAX keyword + const auto& plymaxTables = tableManager.getPlymaxTables(); + const unsigned numMixRegions = plymaxTables.size(); + params.setNumMixRegions(numMixRegions, enablePolymerMolarWeight); + if (!plymaxTables.empty()) { + for (unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++ mixRegionIdx) { + const auto& plymaxTable = plymaxTables.template getTable(mixRegionIdx); + params.plymaxMaxConcentration_[mixRegionIdx] = plymaxTable.getPolymerConcentrationColumn()[0]; + } + } + else { + throw std::runtime_error("PLYMAX must be specified in POLYMER runs\n"); + } + + if (!eclState.getTableManager().getPlmixparTable().empty()) { + if (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) { + params.plymixparToddLongstaff_[mixRegionIdx] = plmixparTable[mixRegionIdx].todd_langstaff; + } + } + } + else if (!enablePolymerMolarWeight) { + throw std::runtime_error("PLMIXPAR must be specified in POLYMER runs\n"); + } + + params.hasPlyshlog_ = eclState.getTableManager().hasTables("PLYSHLOG"); + params.hasShrate_ = eclState.getTableManager().useShrate(); + + if ((params.hasPlyshlog_ || params.hasShrate_) && enablePolymerMolarWeight) { + OpmLog::warning("PLYSHLOG and SHRATE should not be used in POLYMW runs, they will have no effect.\n"); + } + + if (params.hasPlyshlog_ && !enablePolymerMolarWeight) { + const auto& plyshlogTables = tableManager.getPlyshlogTables(); + assert(numPvtRegions == plyshlogTables.size()); + params.plyshlogShearEffectRefMultiplier_.resize(numPvtRegions); + params.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 = params.hasShrate_? unitSystem.parse("1/Time").getSIScaling() : unitSystem.parse("Length/Time").getSIScaling(); + for (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 = params.plyviscViscosityMultiplierTable_[pvtRegionIdx].eval(plyshlogRefPolymerConcentration, /*extrapolate=*/true); + // convert the table using referece conditions + for (size_t i = 0; i < waterVelocity.size(); ++i) { + shearMultiplier[i] *= refViscMult; + shearMultiplier[i] -= 1; + shearMultiplier[i] /= (refViscMult - 1); + shearMultiplier[i] = shearMultiplier[i]; + } + params.plyshlogShearEffectRefMultiplier_[pvtRegionIdx].resize(waterVelocity.size()); + params.plyshlogShearEffectRefLogVelocity_[pvtRegionIdx].resize(waterVelocity.size()); + + for (size_t i = 0; i < waterVelocity.size(); ++i) { + params.plyshlogShearEffectRefMultiplier_[pvtRegionIdx][i] = shearMultiplier[i]; + params.plyshlogShearEffectRefLogVelocity_[pvtRegionIdx][i] = waterVelocity[i]; + } + } + } + + if (params.hasShrate_ && !enablePolymerMolarWeight) { + if (!params.hasPlyshlog_) { + throw std::runtime_error("PLYSHLOG must be specified if SHRATE is used in POLYMER runs\n"); + } + const auto& shrateTable = eclState.getTableManager().getShrateTable(); + params.shrate_.resize(numPvtRegions); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) { + if (shrateTable.empty()) { + params.shrate_[pvtRegionIdx] = 4.8; //default; + } + else if (shrateTable.size() == numPvtRegions) { + params.shrate_[pvtRegionIdx] = shrateTable[pvtRegionIdx].rate; + } + else { + throw std::runtime_error("SHRATE must either have 0 or number of NUMPVT entries\n"); + } + } + } + + if constexpr (enablePolymerMolarWeight) { + const auto& plyvmhTable = eclState.getTableManager().getPlyvmhTable(); + if (!plyvmhTable.empty()) { + assert(plyvmhTable.size() == numMixRegions); + for (size_t regionIdx = 0; regionIdx < numMixRegions; ++regionIdx) { + params.plyvmhCoefficients_[regionIdx].k_mh = plyvmhTable[regionIdx].k_mh; + params.plyvmhCoefficients_[regionIdx].a_mh = plyvmhTable[regionIdx].a_mh; + params.plyvmhCoefficients_[regionIdx].gamma = plyvmhTable[regionIdx].gamma; + params.plyvmhCoefficients_[regionIdx].kappa = plyvmhTable[regionIdx].kappa; + } + } + else { + throw std::runtime_error("PLYVMH keyword must be specified in POLYMW rus \n"); + } + + using TabulatedTwoDFunction = typename BlackOilPolymerParams::TabulatedTwoDFunction; + + // 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(); + TabulatedTwoDFunction tablefunc(throughput, watervelocity, molecularweight, true, false); + params.plymwinjTables_[tableNumber] = std::move(tablefunc); + } + + // 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(); + TabulatedTwoDFunction tablefunc(throughput, watervelocity, skinpressure, true, false); + params.skprwatTables_[tableNumber] = std::move(tablefunc); + } + + // 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(); + typename BlackOilPolymerParams::SkprpolyTable tablefunc = + {refPolymerConcentration, + TabulatedTwoDFunction(throughput, watervelocity, skinpressure, true, false)}; + params.skprpolyTables_[tableNumber] = std::move(tablefunc); + } + } + + return params; +} + template BlackOilBrineParams setupBrineParams(bool, const EclipseState&); template BlackOilBrineParams @@ -421,4 +684,9 @@ setupFoamParams(bool, const EclipseState&); template BlackOilMICPParams setupMICPParams(bool, const EclipseState&); +template BlackOilPolymerParams +setupPolymerParams(bool, const EclipseState&); +template BlackOilPolymerParams +setupPolymerParams(bool, const EclipseState&); + } diff --git a/ebos/eclblackoilmoduleinit.hh b/ebos/eclblackoilmoduleinit.hh index b69cba7045c..9d67be31d42 100644 --- a/ebos/eclblackoilmoduleinit.hh +++ b/ebos/eclblackoilmoduleinit.hh @@ -29,6 +29,7 @@ template struct BlackOilBrineParams; template struct BlackOilExtboParams; template struct BlackOilFoamParams; template struct BlackOilMICPParams; +template struct BlackOilPolymerParams; //! \brief Setup parameters for brine module from an EclipseState. template @@ -50,6 +51,11 @@ template BlackOilMICPParams setupMICPParams(bool enableMICP, const EclipseState& eclState); +//! \brief Setup parameters for polymer module from an EclipseState. +template +BlackOilPolymerParams setupPolymerParams(bool enablePolymer, + const EclipseState& eclState); + } #endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 39b027ccd1b..1337a3c74d3 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -789,8 +789,9 @@ public: ExtboModule::setParams(setupExtboParams(enableExtbo, vanguard.eclState())); FoamModule::setParams(setupFoamParams(enableFoam, vanguard.eclState())); MICPModule::setParams(setupMICPParams(enableMICP, vanguard.eclState())); + PolymerModule::setParams(setupPolymerParams(enablePolymer, + vanguard.eclState())); SolventModule::initFromState(vanguard.eclState(), vanguard.schedule()); - PolymerModule::initFromState(vanguard.eclState()); // create the ECL writer eclWriter_.reset(new EclWriterType(simulator)); From 94492e181030177c536ae1c0e60e4d8795a80e7d Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Sep 2022 11:41:39 +0200 Subject: [PATCH 6/6] import code for initializing solvent module from eclipsestate this allows us to keep this code in a separate compile unit, limiting the amount of times it has to be built. --- ebos/eclblackoilmoduleinit.cc | 242 ++++++++++++++++++++++++++++++++++ ebos/eclblackoilmoduleinit.hh | 8 ++ ebos/eclproblem.hh | 4 +- 3 files changed, 253 insertions(+), 1 deletion(-) diff --git a/ebos/eclblackoilmoduleinit.cc b/ebos/eclblackoilmoduleinit.cc index bac1c784d19..a06e1c62836 100644 --- a/ebos/eclblackoilmoduleinit.cc +++ b/ebos/eclblackoilmoduleinit.cc @@ -51,6 +51,7 @@ #include #include #include +#include #include #include @@ -670,6 +671,244 @@ BlackOilPolymerParams setupPolymerParams(bool enablePolymer, return params; } +template +BlackOilSolventParams setupSolventParams(bool enableSolvent, + const EclipseState& eclState, + const Schedule& schedule) +{ + BlackOilSolventParams params; + // some sanity checks: if solvents are enabled, the SOLVENT keyword must be + // present, if solvents are disabled the keyword must not be present. + if (enableSolvent && !eclState.runspec().phases().active(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 (!enableSolvent && eclState.runspec().phases().active(Phase::SOLVENT)) + throw std::runtime_error("Solvent treatment disabled at compile time, but the deck " + "contains the SOLVENT keyword"); + + if (!eclState.runspec().phases().active(Phase::SOLVENT)) + return params; // solvent treatment is supposed to be disabled + + params.solventPvt_.initFromState(eclState, schedule); + + const auto& tableManager = eclState.getTableManager(); + // initialize the objects which deal with the SSFN keyword + const auto& ssfnTables = tableManager.getSsfnTables(); + unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + params.setNumSatRegions(numSatRegions); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) { + const auto& ssfnTable = ssfnTables.template getTable(satRegionIdx); + params.ssfnKrg_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(), + ssfnTable.getGasRelPermMultiplierColumn(), + /*sortInput=*/true); + params.ssfnKrs_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(), + ssfnTable.getSolventRelPermMultiplierColumn(), + /*sortInput=*/true); + } + + // initialize the objects needed for miscible solvent and oil simulations + params.isMiscible_ = false; + if (!eclState.getTableManager().getMiscTables().empty()) { + params.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 + params.sof2Krn_.resize(numSatRegions); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) { + const auto& sof2Table = sof2Tables.template getTable(satRegionIdx); + params.sof2Krn_[satRegionIdx].setXYContainers(sof2Table.getSoColumn(), + sof2Table.getKroColumn(), + /*sortInput=*/true); + } + } + else + throw std::runtime_error("SOF2 must be specified in MISCIBLE (SOLVENT) runs\n"); + + const auto& miscTables = tableManager.getMiscTables(); + if (!miscTables.empty()) { + assert(numMiscRegions == miscTables.size()); + + // resize the attributes of the object + params.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(); + params.misc_[miscRegionIdx].setXYContainers(solventFraction, misc); + } + } + else + throw std::runtime_error("MISC must be specified in MISCIBLE (SOLVENT) runs\n"); + + using TabulatedFunction = typename BlackOilSolventParams::TabulatedFunction; + + // resize the attributes of the object + params.pmisc_.resize(numMiscRegions); + const auto& pmiscTables = tableManager.getPmiscTables(); + if (!pmiscTables.empty()) { + assert(numMiscRegions == pmiscTables.size()); + + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + const auto& pmiscTable = pmiscTables.template getTable(regionIdx); + + // Copy data + const auto& po = pmiscTable.getOilPhasePressureColumn(); + const auto& pmisc = pmiscTable.getMiscibilityColumn(); + + params.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) { + params.pmisc_[regionIdx] = constant; + } + } + + // miscible relative permeability multipleiers + params.msfnKrsg_.resize(numSatRegions); + params.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(); + + params.msfnKrsg_[regionIdx].setXYContainers(Ssg, krsg); + params.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) { + params.setMsfn(regionIdx, unit, invUnit); + } + } + // resize the attributes of the object + params.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(); + + params.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) { + params.sorwmis_[regionIdx] = zero; + } + } + + // resize the attributes of the object + params.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); + + // Copy data + const auto& sw = sgcwmisTable.getWaterSaturationColumn(); + const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn(); + + params.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) + params.sgcwmis_[regionIdx] = zero; + } + + const auto& tlmixpar = eclState.getTableManager().getTLMixpar(); + if (!tlmixpar.empty()) { + // resize the attributes of the object + params.tlMixParamViscosity_.resize(numMiscRegions); + params.tlMixParamDensity_.resize(numMiscRegions); + + assert(numMiscRegions == tlmixpar.size()); + for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) { + const auto& tlp = tlmixpar[regionIdx]; + params.tlMixParamViscosity_[regionIdx] = tlp.viscosity_parameter; + params.tlMixParamDensity_[regionIdx] = tlp.density_parameter; + } + } + else + throw std::runtime_error("TLMIXPAR must be specified in MISCIBLE (SOLVENT) runs\n"); + + // resize the attributes of the object + params.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(); + + params.tlPMixTable_[regionIdx].setXYContainers(po, tlpmixpa); + } + } + else { + // if empty keyword. Try to use the pmisc table as default. + if (params.pmisc_.size() > 0) + params.tlPMixTable_ = params.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."); + } + } + 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) + params.tlPMixTable_[regionIdx] = ones; + } + } + + return params; +} + template BlackOilBrineParams setupBrineParams(bool, const EclipseState&); template BlackOilBrineParams @@ -689,4 +928,7 @@ setupPolymerParams(bool, const EclipseState&); template BlackOilPolymerParams setupPolymerParams(bool, const EclipseState&); +template BlackOilSolventParams +setupSolventParams(bool, const EclipseState&, const Schedule&); + } diff --git a/ebos/eclblackoilmoduleinit.hh b/ebos/eclblackoilmoduleinit.hh index 9d67be31d42..49d51583708 100644 --- a/ebos/eclblackoilmoduleinit.hh +++ b/ebos/eclblackoilmoduleinit.hh @@ -25,11 +25,13 @@ namespace Opm { class EclipseState; +class Schedule; template struct BlackOilBrineParams; template struct BlackOilExtboParams; template struct BlackOilFoamParams; template struct BlackOilMICPParams; template struct BlackOilPolymerParams; +template struct BlackOilSolventParams; //! \brief Setup parameters for brine module from an EclipseState. template @@ -56,6 +58,12 @@ template BlackOilPolymerParams setupPolymerParams(bool enablePolymer, const EclipseState& eclState); +//! \brief Setup parameters for solvent module from an EclipseState. +template +BlackOilSolventParams setupSolventParams(bool enableSolvent, + const EclipseState& eclState, + const Schedule& schedule); + } #endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 1337a3c74d3..a94aae1d5dc 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -791,7 +791,9 @@ public: MICPModule::setParams(setupMICPParams(enableMICP, vanguard.eclState())); PolymerModule::setParams(setupPolymerParams(enablePolymer, vanguard.eclState())); - SolventModule::initFromState(vanguard.eclState(), vanguard.schedule()); + SolventModule::setParams(setupSolventParams(enableSolvent, + vanguard.eclState(), + vanguard.schedule())); // create the ECL writer eclWriter_.reset(new EclWriterType(simulator));