diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 173893532..01ace7fcf 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -55,6 +55,12 @@ endmacro() # originally generated with the command: # find opm -name '*.c*' -printf '\t%p\n' | sort list (APPEND MAIN_SOURCE_FILES + opm/models/blackoil/blackoilbrineparams.cpp + opm/models/blackoil/blackoilextboparams.cpp + opm/models/blackoil/blackoilfoamparams.cpp + opm/models/blackoil/blackoilmicpparams.cpp + opm/models/blackoil/blackoilpolymerparams.cpp + opm/models/blackoil/blackoilsolventparams.cpp opm/simulators/flow/ActionHandler.cpp opm/simulators/flow/Banners.cpp opm/simulators/flow/BlackoilModelParameters.cpp @@ -516,34 +522,34 @@ list (APPEND TEST_DATA_FILES list (APPEND PUBLIC_HEADER_FILES opm/models/blackoil/blackoilboundaryratevector.hh opm/models/blackoil/blackoilbrinemodules.hh - opm/models/blackoil/blackoilbrineparams.hh + opm/models/blackoil/blackoilbrineparams.hpp opm/models/blackoil/blackoildarcyfluxmodule.hh opm/models/blackoil/blackoildiffusionmodule.hh opm/models/blackoil/blackoildispersionmodule.hh opm/models/blackoil/blackoilenergymodules.hh opm/models/blackoil/blackoilextbomodules.hh - opm/models/blackoil/blackoilextboparams.hh + opm/models/blackoil/blackoilextboparams.hpp opm/models/blackoil/blackoilextensivequantities.hh opm/models/blackoil/blackoilfoammodules.hh - opm/models/blackoil/blackoilfoamparams.hh + opm/models/blackoil/blackoilfoamparams.hpp opm/models/blackoil/blackoilindices.hh opm/models/blackoil/blackoilintensivequantities.hh opm/models/blackoil/blackoillocalresidual.hh opm/models/blackoil/blackoillocalresidualtpfa.hh opm/models/blackoil/blackoilmicpmodules.hh - opm/models/blackoil/blackoilmicpparams.hh + opm/models/blackoil/blackoilmicpparams.hpp opm/models/blackoil/blackoilmodel.hh opm/models/blackoil/blackoilnewtonmethod.hh opm/models/blackoil/blackoilnewtonmethodparameters.hh opm/models/blackoil/blackoilonephaseindices.hh opm/models/blackoil/blackoilpolymermodules.hh - opm/models/blackoil/blackoilpolymerparams.hh + opm/models/blackoil/blackoilpolymerparams.hpp opm/models/blackoil/blackoilprimaryvariables.hh opm/models/blackoil/blackoilproblem.hh opm/models/blackoil/blackoilproperties.hh opm/models/blackoil/blackoilratevector.hh opm/models/blackoil/blackoilsolventmodules.hh - opm/models/blackoil/blackoilsolventparams.hh + opm/models/blackoil/blackoilsolventparams.hpp opm/models/blackoil/blackoiltwophaseindices.hh opm/models/common/darcyfluxmodule.hh opm/models/common/diffusionmodule.hh diff --git a/opm/models/blackoil/blackoilbrinemodules.hh b/opm/models/blackoil/blackoilbrinemodules.hh index 814dc0784..fda6f814e 100644 --- a/opm/models/blackoil/blackoilbrinemodules.hh +++ b/opm/models/blackoil/blackoilbrinemodules.hh @@ -28,24 +28,16 @@ #ifndef EWOMS_BLACK_OIL_BRINE_MODULE_HH #define EWOMS_BLACK_OIL_BRINE_MODULE_HH -#include "blackoilproperties.hh" +#include +#include -#include +#include -#if HAVE_ECL_INPUT -#include -#include -#include -#include -#include -#include -#include -#endif +#include #include #include -#include namespace Opm { /*! @@ -85,79 +77,11 @@ class BlackOilBrineModule static constexpr unsigned numPhases = FluidSystem::numPhases; public: - -#if HAVE_ECL_INPUT - /*! - * \brief Initialize all internal data structures needed by the brine module - */ - static void initFromState(const EclipseState& eclState) + //! \brief Set parameters. + static void setParams(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; // 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); - params_.saltdenTable_.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(); - params_.saltdenTable_[pvtRegionIdx] = saltsolTable.getSaltdenColumn().front(); - } - } - - const TableContainer& pcfactTables = tableManager.getPcfactTables(); - if (!pcfactTables.empty()) { - unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); - params_.pcfactTable_.resize(numSatRegions); - for (size_t i = 0; i < pcfactTables.size(); ++i) { - const PcfactTable& pcfactTable = pcfactTables.getTable(i); - params_.pcfactTable_[i].setXYContainers(pcfactTable.getPorosityChangeColumn(), pcfactTable.getPcMultiplierColumn()); - } - } - } + params_ = params; } -#endif /*! * \brief Register all run-time parameters for the black-oil brine module. diff --git a/opm/models/blackoil/blackoilbrineparams.cpp b/opm/models/blackoil/blackoilbrineparams.cpp new file mode 100644 index 000000000..4173d5f3f --- /dev/null +++ b/opm/models/blackoil/blackoilbrineparams.cpp @@ -0,0 +1,134 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + +#if HAVE_ECL_INPUT +#include +#include +#include +#include +#include +#include +#include +#endif + +#include +#include +#include + +namespace Opm { + +#if HAVE_ECL_INPUT +template +template +void BlackOilBrineParams:: +initFromState(const 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)) { + 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)) { + 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; // brine treatment is supposed to be disabled + } + + const auto& tableManager = eclState.getTableManager(); + + unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + 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()) { + 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(); + bdensityTable_[pvtRegionIdx].setXYContainers(c, bdensityTable); + } + } + + if constexpr (enableSaltPrecipitation) { + const TableContainer& permfactTables = tableManager.getPermfactTables(); + permfactTable_.resize(numPvtRegions); + for (std::size_t i = 0; i < permfactTables.size(); ++i) { + const PermfactTable& permfactTable = permfactTables.getTable(i); + permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn()); + } + + const TableContainer& saltsolTables = tableManager.getSaltsolTables(); + if (!saltsolTables.empty()) { + saltsolTable_.resize(numPvtRegions); + saltdenTable_.resize(numPvtRegions); + assert(numPvtRegions == saltsolTables.size()); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { + const SaltsolTable& saltsolTable = saltsolTables.getTable(pvtRegionIdx ); + saltsolTable_[pvtRegionIdx] = saltsolTable.getSaltsolColumn().front(); + saltdenTable_[pvtRegionIdx] = saltsolTable.getSaltdenColumn().front(); + } + } + + const TableContainer& pcfactTables = tableManager.getPcfactTables(); + if (!pcfactTables.empty()) { + unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + pcfactTable_.resize(numSatRegions); + for (std::size_t i = 0; i < pcfactTables.size(); ++i) { + const PcfactTable& pcfactTable = pcfactTables.getTable(i); + pcfactTable_[i].setXYContainers(pcfactTable.getPorosityChangeColumn(), pcfactTable.getPcMultiplierColumn()); + } + } + } +} +#endif + +#define INSTANTIATE_TYPE(T) \ + template struct BlackOilBrineParams; \ + template void BlackOilBrineParams::initFromState(const EclipseState&); \ + template void BlackOilBrineParams::initFromState(const EclipseState&); \ + template void BlackOilBrineParams::initFromState(const EclipseState&); \ + template void BlackOilBrineParams::initFromState(const EclipseState&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/models/blackoil/blackoilbrineparams.hh b/opm/models/blackoil/blackoilbrineparams.hpp similarity index 82% rename from opm/models/blackoil/blackoilbrineparams.hh rename to opm/models/blackoil/blackoilbrineparams.hpp index 67bf50f81..dafa838d2 100644 --- a/opm/models/blackoil/blackoilbrineparams.hh +++ b/opm/models/blackoil/blackoilbrineparams.hpp @@ -25,8 +25,8 @@ * * \brief Contains the parameters required to extend the black-oil model by brine. */ -#ifndef EWOMS_BLACK_OIL_BRINE_PARAMS_HH -#define EWOMS_BLACK_OIL_BRINE_PARAMS_HH +#ifndef OPM_BLACK_OIL_BRINE_PARAMS_HPP +#define OPM_BLACK_OIL_BRINE_PARAMS_HPP #include @@ -34,9 +34,19 @@ namespace Opm { +#if HAVE_ECL_INPUT +class EclipseState; +#endif + //! \brief Struct holding the parameters for the BlackoilBrineModule class. template -struct BlackOilBrineParams { +struct BlackOilBrineParams +{ +#if HAVE_ECL_INPUT + template + void initFromState(const EclipseState& eclState); +#endif + using TabulatedFunction = Tabulated1DFunction; std::vector bdensityTable_; @@ -49,4 +59,4 @@ struct BlackOilBrineParams { } // namespace Opm -#endif +#endif // OPM_BLACK_OIL_BRINE_PARAMS_HPP diff --git a/opm/models/blackoil/blackoilextbomodules.hh b/opm/models/blackoil/blackoilextbomodules.hh index c8cc28dbf..188c69c31 100644 --- a/opm/models/blackoil/blackoilextbomodules.hh +++ b/opm/models/blackoil/blackoilextbomodules.hh @@ -32,30 +32,20 @@ #ifndef EWOMS_BLACK_OIL_EXTBO_MODULE_HH #define EWOMS_BLACK_OIL_EXTBO_MODULE_HH -#include "blackoilproperties.hh" +#include +#include -#include +#include + +#include //#include //TODO: Missing ... -#if HAVE_ECL_INPUT -#include -#include -#include -#include -#include -#include -#include -#include -#include -#endif - #include #include #include #include -#include namespace Opm { @@ -96,164 +86,11 @@ class BlackOilExtboModule static constexpr bool blackoilConserveSurfaceVolume = getPropValue(); public: -#if HAVE_ECL_INPUT - /*! - * \brief Initialize all internal data structures needed by the solvent module - */ - static void initFromState(const EclipseState& eclState) + //! \brief Set parameters. + static void setParams(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; // 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(); - - 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"); + params_ = params; } -#endif /*! * \brief Register all run-time parameters for the black-oil solvent module. diff --git a/opm/models/blackoil/blackoilextboparams.cpp b/opm/models/blackoil/blackoilextboparams.cpp new file mode 100644 index 000000000..efd27d0e4 --- /dev/null +++ b/opm/models/blackoil/blackoilextboparams.cpp @@ -0,0 +1,220 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + +#if HAVE_ECL_INPUT +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +#include +#include + +namespace Opm { + +#if HAVE_ECL_INPUT +template +template +void BlackOilExtboParams:: +initFromState(const 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)) { + 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; // 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(); + + BO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + BG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + X_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + Y_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + VISCO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + VISCG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + + PBUB_RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + PBUB_RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme}); + + zLim_.resize(numPvtRegions); + + const bool extractCmpFromPvt = true; //: Default values used in [*] + oilCmp_.resize(numPvtRegions); + 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); //-------------"------------- + 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; + + BO_[regionIdx].appendXPos(ZCO2); + BG_[regionIdx].appendXPos(ZCO2); + + RS_[regionIdx].appendXPos(ZCO2); + RV_[regionIdx].appendXPos(ZCO2); + + X_[regionIdx].appendXPos(ZCO2); + Y_[regionIdx].appendXPos(ZCO2); + + VISCO_[regionIdx].appendXPos(ZCO2); + VISCG_[regionIdx].appendXPos(ZCO2); + + PBUB_RS_[regionIdx].appendXPos(ZCO2); + PBUB_RV_[regionIdx].appendXPos(ZCO2); + + const auto& underSaturatedTable = pvtsolTable.getUnderSaturatedTable(outerIdx); + std::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; + 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; + } + + BO_[regionIdx].appendSamplePoint(outerIdx, po, bo); + BG_[regionIdx].appendSamplePoint(outerIdx, po, bg); + RS_[regionIdx].appendSamplePoint(outerIdx, po, rs); + RV_[regionIdx].appendSamplePoint(outerIdx, po, rv); + X_[regionIdx].appendSamplePoint(outerIdx, po, xv); + Y_[regionIdx].appendSamplePoint(outerIdx, po, yv); + VISCO_[regionIdx].appendSamplePoint(outerIdx, po, mo); + VISCG_[regionIdx].appendSamplePoint(outerIdx, po, mg); + break; + } + + bo0 = bo; + po0 = po; + + BO_[regionIdx].appendSamplePoint(outerIdx, po, bo); + BG_[regionIdx].appendSamplePoint(outerIdx, po, bg); + + RS_[regionIdx].appendSamplePoint(outerIdx, po, rs); + RV_[regionIdx].appendSamplePoint(outerIdx, po, rv); + + X_[regionIdx].appendSamplePoint(outerIdx, po, xv); + Y_[regionIdx].appendSamplePoint(outerIdx, po, yv); + + VISCO_[regionIdx].appendSamplePoint(outerIdx, po, mo); + VISCG_[regionIdx].appendSamplePoint(outerIdx, po, mg); + + // rs,rv -> pressure + PBUB_RS_[regionIdx].appendSamplePoint(outerIdx, rs, po); + PBUB_RV_[regionIdx].appendSamplePoint(outerIdx, rv, po); + } + } + oilCmp_[regionIdx].setXYContainers(zArg, oilCmp, /*sortInput=*/false); + 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"); +} +#endif + +#define INSTANTIATE_TYPE(T) \ + template struct BlackOilExtboParams; \ + template void BlackOilExtboParams::initFromState(const EclipseState&); \ + template void BlackOilExtboParams::initFromState(const EclipseState&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/models/blackoil/blackoilextboparams.hh b/opm/models/blackoil/blackoilextboparams.hpp similarity index 88% rename from opm/models/blackoil/blackoilextboparams.hh rename to opm/models/blackoil/blackoilextboparams.hpp index 013c24665..656596d32 100644 --- a/opm/models/blackoil/blackoilextboparams.hh +++ b/opm/models/blackoil/blackoilextboparams.hpp @@ -29,8 +29,8 @@ * for CO2 EOR Simulations.” in ECMOR XVII – The 17th European Conference on the * Mathematics of Oil Recovery, September 2020. */ -#ifndef EWOMS_BLACK_OIL_EXTBO_PARAMS_HH -#define EWOMS_BLACK_OIL_EXTBO_PARAMS_HH +#ifndef OPM_BLACK_OIL_EXTBO_PARAMS_HPP +#define OPM_BLACK_OIL_EXTBO_PARAMS_HPP #include #include @@ -39,12 +39,23 @@ namespace Opm { +#if HAVE_ECL_INPUT +class EclipseState; +#endif + //! \brief Struct holding the parameters for the BlackoilExtboModule class. template -struct BlackOilExtboParams { +struct BlackOilExtboParams +{ using TabulatedFunction = Tabulated1DFunction; using Tabulated2DFunction = UniformXTabulated2DFunction; +#if HAVE_ECL_INPUT + template + void initFromState(const EclipseState& eclState); +#endif + + std::vector X_; std::vector Y_; std::vector PBUB_RS_; @@ -65,4 +76,4 @@ struct BlackOilExtboParams { } // namespace Opm -#endif +#endif // OPM_BLACK_OIL_EXTBO_PARAMS_HPP diff --git a/opm/models/blackoil/blackoilfoammodules.hh b/opm/models/blackoil/blackoilfoammodules.hh index 737d554b8..17ca0d028 100644 --- a/opm/models/blackoil/blackoilfoammodules.hh +++ b/opm/models/blackoil/blackoilfoammodules.hh @@ -28,23 +28,24 @@ #ifndef EWOMS_BLACK_OIL_FOAM_MODULE_HH #define EWOMS_BLACK_OIL_FOAM_MODULE_HH -#include "blackoilproperties.hh" - #include #include +#include + +#include +#include + +#include +#include + #if HAVE_ECL_INPUT #include #include #include #endif -#include - -#include -#include - #include namespace Opm { @@ -87,93 +88,11 @@ class BlackOilFoamModule enum { enableSolvent = getPropValue() }; public: -#if HAVE_ECL_INPUT - /*! - * \brief Initialize all internal data structures needed by the foam module - */ - static void initFromState(const EclipseState& eclState) + //! \brief Set parameters. + static void setParams(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; // foam treatment is supposed to be disabled - } - - params_.transport_phase_ = eclState.getInitConfig().getFoamConfig().getTransportPhase(); - - if (eclState.getInitConfig().getFoamConfig().getMobilityModel() != FoamConfig::MobilityModel::TAB) { - throw std::runtime_error("In FOAMOPTS, 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); - } + params_ = params; } -#endif /*! * \brief Register all run-time parameters for the black-oil foam module. diff --git a/opm/models/blackoil/blackoilfoamparams.cpp b/opm/models/blackoil/blackoilfoamparams.cpp new file mode 100644 index 000000000..71f3d779b --- /dev/null +++ b/opm/models/blackoil/blackoilfoamparams.cpp @@ -0,0 +1,150 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + +#if HAVE_ECL_INPUT +#include +#include +#include +#endif + +#include +#include + +namespace Opm { + +#if HAVE_ECL_INPUT +template +template +void BlackOilFoamParams:: +initFromState(const EclipseState& eclState) +{ + // some sanity checks: if foam is enabled, the FOAM keyword must be + // present, if foam is disabled the keyword must not be present. + if constexpr (enableFoam) { + if (!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 (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; // foam treatment is supposed to be disabled + } + + transport_phase_ = eclState.getInitConfig().getFoamConfig().getTransportPhase(); + + if (eclState.getInitConfig().getFoamConfig().getMobilityModel() != FoamConfig::MobilityModel::TAB) { + throw std::runtime_error("In FOAMOPTS, only TAB is allowed for the gas mobility factor reduction model."); + } + + const auto& tableManager = eclState.getTableManager(); + const unsigned int numSatRegions = tableManager.getTabdims().getNumSatTables(); + setNumSatRegions(numSatRegions); + const unsigned int numPvtRegions = tableManager.getTabdims().getNumPVTTables(); + 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); + foamCoefficients_[satReg] = FoamCoefficients(); + foamCoefficients_[satReg].fm_min = rec.minimumSurfactantConcentration(); + foamCoefficients_[satReg].fm_surf = rec.referenceSurfactantConcentration(); + foamCoefficients_[satReg].ep_surf = rec.exponent(); + foamRockDensity_[satReg] = rec.rockDensity(); + foamAllowDesorption_[satReg] = rec.allowDesorption(); + const auto& foamadsTable = foamadsTables.template getTable(satReg); + const auto& conc = foamadsTable.getFoamConcentrationColumn(); + const auto& ads = foamadsTable.getAdsorbedFoamColumn(); + 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(); + gasMobilityMultiplierTable_[pvtReg].setXYContainers(conc, mobMult); + } +} +#endif + +template +void BlackOilFoamParams::setNumSatRegions(unsigned numRegions) +{ + foamCoefficients_.resize(numRegions); + foamRockDensity_.resize(numRegions); + foamAllowDesorption_.resize(numRegions); + adsorbedFoamTable_.resize(numRegions); +} + +#define INSTANTIATE_TYPE(T) \ + template struct BlackOilFoamParams; \ + template void BlackOilFoamParams::initFromState(const EclipseState&); \ + template void BlackOilFoamParams::initFromState(const EclipseState&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/models/blackoil/blackoilfoamparams.hh b/opm/models/blackoil/blackoilfoamparams.hpp similarity index 86% rename from opm/models/blackoil/blackoilfoamparams.hh rename to opm/models/blackoil/blackoilfoamparams.hpp index 7f9e0b4f8..99275d383 100644 --- a/opm/models/blackoil/blackoilfoamparams.hh +++ b/opm/models/blackoil/blackoilfoamparams.hpp @@ -25,31 +25,36 @@ * * \brief Contains the parameters to extend the black-oil model to include the effects of foam. */ -#ifndef EWOMS_BLACK_OIL_FOAM_PARAMS_HH -#define EWOMS_BLACK_OIL_FOAM_PARAMS_HH +#ifndef OPM_BLACK_OIL_FOAM_PARAMS_HPP +#define OPM_BLACK_OIL_FOAM_PARAMS_HPP -#include #include #include namespace Opm { +enum class Phase; + +#if HAVE_ECL_INPUT +class EclipseState; +#endif + //! \brief Struct holding the parameters for the BlackoilFoamModule class. template -struct BlackOilFoamParams { +struct BlackOilFoamParams +{ using TabulatedFunction = Tabulated1DFunction; +#if HAVE_ECL_INPUT + template + void initFromState(const EclipseState& eclState); +#endif + /*! * \brief Specify the number of saturation regions. */ - void setNumSatRegions(unsigned numRegions) - { - foamCoefficients_.resize(numRegions); - foamRockDensity_.resize(numRegions); - foamAllowDesorption_.resize(numRegions); - adsorbedFoamTable_.resize(numRegions); - } + void setNumSatRegions(unsigned numRegions); // a struct containing constants to calculate change to relative permeability, // based on model (1-9) in Table 1 of @@ -84,4 +89,4 @@ struct BlackOilFoamParams { } // namespace Opm -#endif +#endif // OPM_BLACK_OIL_FOAM_PARAMS_HPP diff --git a/opm/models/blackoil/blackoilmicpmodules.hh b/opm/models/blackoil/blackoilmicpmodules.hh index da49371a1..f027237f6 100644 --- a/opm/models/blackoil/blackoilmicpmodules.hh +++ b/opm/models/blackoil/blackoilmicpmodules.hh @@ -28,22 +28,15 @@ #ifndef EWOMS_BLACK_OIL_MICP_MODULE_HH #define EWOMS_BLACK_OIL_MICP_MODULE_HH -#include "blackoilproperties.hh" - -#include -#include - -#if HAVE_ECL_INPUT -#include -#include -#endif - #include -#include +#include +#include + +#include + #include #include -#include namespace Opm { /*! @@ -86,56 +79,11 @@ class BlackOilMICPModule static constexpr unsigned numEq = getPropValue(); public: - -#if HAVE_ECL_INPUT - // - //* \brief Initialize all internal data structures needed by the MICP module - // - static void initFromState(const EclipseState& eclState) + //! \brief Set parameters. + static void setParams(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; // MICP treatment is supposed to be disabled*/ - - // initialize the objects which deal with the MICPpara keyword - const auto& MICPpara = eclState.getMICPpara(); - setMICPpara(MICPpara.getDensityBiofilm(), - MICPpara.getDensityCalcite(), - MICPpara.getDetachmentRate(), - MICPpara.getCriticalPorosity(), - MICPpara.getFittingFactor(), - MICPpara.getHalfVelocityOxygen(), - MICPpara.getHalfVelocityUrea(), - MICPpara.getMaximumGrowthRate(), - MICPpara.getMaximumUreaUtilization(), - MICPpara.getMicrobialAttachmentRate(), - MICPpara.getMicrobialDeathRate(), - MICPpara.getMinimumPermeability(), - MICPpara.getOxygenConsumptionFactor(), - MICPpara.getYieldGrowthCoefficient(), - MICPpara.getMaximumOxygenConcentration(), - MICPpara.getMaximumUreaConcentration(), - MICPpara.getToleranceBeforeClogging()); - // obtain the porosity for the clamp in the blackoilnewtonmethod - if constexpr (std::is_same_v) { - const auto phi = eclState.fieldProps().get_double("PORO"); - params_.phi_.resize(phi.size()); - std::copy(phi.begin(), phi.end(), params_.phi_.begin()); - } else { - params_.phi_ = eclState.fieldProps().get_double("PORO"); - } + params_ = params; } -#endif /*! * \brief The simulator stops if "clogging" has been (almost) reached in any of the cells. @@ -151,48 +99,6 @@ public: throw std::logic_error("Clogging has been (almost) reached in at least one cell\n"); } - /*! - * \brief Specify the MICP properties a single region. - * - * The index of specified here must be in range [0, numSatRegions) - */ - static void setMICPpara(const Scalar& densityBiofilm, - const Scalar& densityCalcite, - const Scalar& detachmentRate, - const Scalar& criticalPorosity, - const Scalar& fittingFactor, - const Scalar& halfVelocityOxygen, - const Scalar& halfVelocityUrea, - const Scalar& maximumGrowthRate, - const Scalar& maximumUreaUtilization, - const Scalar& microbialAttachmentRate, - const Scalar& microbialDeathRate, - const Scalar& minimumPermeability, - const Scalar& oxygenConsumptionFactor, - const Scalar& yieldGrowthCoefficient, - const Scalar& maximumOxygenConcentration, - const Scalar& maximumUreaConcentration, - const Scalar& toleranceBeforeClogging) - { - params_.densityBiofilm_ = densityBiofilm; - params_.densityCalcite_ = densityCalcite; - params_.detachmentRate_ = detachmentRate; - params_.criticalPorosity_ = criticalPorosity; - params_.fittingFactor_ = fittingFactor; - params_.halfVelocityOxygen_ = halfVelocityOxygen; - params_.halfVelocityUrea_ = halfVelocityUrea; - params_.maximumGrowthRate_ = maximumGrowthRate; - params_.maximumUreaUtilization_ = maximumUreaUtilization; - params_.microbialAttachmentRate_ = microbialAttachmentRate; - params_.microbialDeathRate_ = microbialDeathRate; - params_.minimumPermeability_ = minimumPermeability; - params_.oxygenConsumptionFactor_ = oxygenConsumptionFactor; - params_.yieldGrowthCoefficient_ = yieldGrowthCoefficient; - params_.maximumOxygenConcentration_ = maximumOxygenConcentration; - params_.maximumUreaConcentration_ = maximumUreaConcentration; - params_.toleranceBeforeClogging_ = toleranceBeforeClogging; - } - /*! * \brief Register all run-time parameters for the black-oil MICP module. */ diff --git a/opm/models/blackoil/blackoilmicpparams.cpp b/opm/models/blackoil/blackoilmicpparams.cpp new file mode 100644 index 000000000..9afaa55ed --- /dev/null +++ b/opm/models/blackoil/blackoilmicpparams.cpp @@ -0,0 +1,104 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + +#if HAVE_ECL_INPUT +#include +#include +#endif + +#include +#include +#include + +namespace Opm { + +#if HAVE_ECL_INPUT +template +template +void BlackOilMICPParams:: +initFromState(const EclipseState& eclState) +{ + // some sanity checks: if MICP is enabled, the MICP keyword must be + // present, if MICP is disabled the keyword must not be present. + if constexpr (enableMICP) { + if (!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 (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; // MICP treatment is supposed to be disabled*/ + + // initialize the objects which deal with the MICPpara keyword + const auto& MICPpara = eclState.getMICPpara(); + densityBiofilm_ = MICPpara.getDensityBiofilm(); + densityCalcite_ = MICPpara.getDensityCalcite(); + detachmentRate_ = MICPpara.getDetachmentRate(); + criticalPorosity_ = MICPpara.getCriticalPorosity(); + fittingFactor_ = MICPpara.getFittingFactor(); + halfVelocityOxygen_ = MICPpara.getHalfVelocityOxygen(); + halfVelocityUrea_ = MICPpara.getHalfVelocityUrea(); + maximumGrowthRate_ = MICPpara.getMaximumGrowthRate(); + maximumUreaUtilization_ = MICPpara.getMaximumUreaUtilization(); + microbialAttachmentRate_ = MICPpara.getMicrobialAttachmentRate(); + microbialDeathRate_ = MICPpara.getMicrobialDeathRate(); + minimumPermeability_ = MICPpara.getMinimumPermeability(); + oxygenConsumptionFactor_ = MICPpara.getOxygenConsumptionFactor(); + yieldGrowthCoefficient_ = MICPpara.getYieldGrowthCoefficient(); + maximumOxygenConcentration_ = MICPpara.getMaximumOxygenConcentration(); + maximumUreaConcentration_ = MICPpara.getMaximumUreaConcentration(); + toleranceBeforeClogging_ = MICPpara.getToleranceBeforeClogging(); + + // obtain the porosity for the clamp in the blackoilnewtonmethod + if constexpr (std::is_same_v) { + const auto phi = eclState.fieldProps().get_double("PORO"); + phi_.resize(phi.size()); + std::copy(phi.begin(), phi.end(), phi_.begin()); + } else { + phi_ = eclState.fieldProps().get_double("PORO"); + } +} +#endif + +#define INSTANTIATE_TYPE(T) \ + template struct BlackOilMICPParams; \ + template void BlackOilMICPParams::initFromState(const EclipseState&); \ + template void BlackOilMICPParams::initFromState(const EclipseState&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/models/blackoil/blackoilmicpparams.hh b/opm/models/blackoil/blackoilmicpparams.hpp similarity index 85% rename from opm/models/blackoil/blackoilmicpparams.hh rename to opm/models/blackoil/blackoilmicpparams.hpp index ef25010e0..13dcc67b5 100644 --- a/opm/models/blackoil/blackoilmicpparams.hh +++ b/opm/models/blackoil/blackoilmicpparams.hpp @@ -25,16 +25,26 @@ * * \brief Contains the parameters required to extend the black-oil model by MICP. */ -#ifndef EWOMS_BLACK_OIL_MICP_PARAMS_HH -#define EWOMS_BLACK_OIL_MICP_PARAMS_HH +#ifndef OPM_BLACK_OIL_MICP_PARAMS_HPP +#define OPM_BLACK_OIL_MICP_PARAMS_HPP #include namespace Opm { +#if HAVE_ECL_INPUT +class EclipseState; +#endif + //! \brief Struct holding the parameters for the BlackOilMICPModule class. template -struct BlackOilMICPParams { +struct BlackOilMICPParams +{ +#if HAVE_ECL_INPUT + template + void initFromState(const EclipseState& eclState); +#endif + Scalar densityBiofilm_; Scalar densityCalcite_; Scalar detachmentRate_; @@ -57,4 +67,4 @@ struct BlackOilMICPParams { } // namespace Opm -#endif +#endif // OPM_BLACK_OIL_MICP_PARAMS_HPP diff --git a/opm/models/blackoil/blackoilpolymermodules.hh b/opm/models/blackoil/blackoilpolymermodules.hh index e5d6ef21e..bf21eb4f8 100644 --- a/opm/models/blackoil/blackoilpolymermodules.hh +++ b/opm/models/blackoil/blackoilpolymermodules.hh @@ -28,23 +28,16 @@ #ifndef EWOMS_BLACK_OIL_POLYMER_MODULE_HH #define EWOMS_BLACK_OIL_POLYMER_MODULE_HH -#include "blackoilproperties.hh" - -#include -#include +#include #include -#if HAVE_ECL_INPUT -#include -#include -#include -#include -#include -#include -#endif +#include +#include -#include +#include + +#include #include #include @@ -91,259 +84,11 @@ class BlackOilPolymerModule static constexpr unsigned numPhases = FluidSystem::numPhases; public: -#if HAVE_ECL_INPUT - /*! - * \brief Initialize all internal data structures needed by the polymer module - */ - static void initFromState(const EclipseState& eclState) + //! \brief Set parameters. + static void setParams(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; // 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"); - } - - // 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); - } - } + params_ = params; } -#endif - /*! * \brief get the PLYMWINJ table */ @@ -1133,7 +878,6 @@ public: { throw std::runtime_error("waterShearFactor() called but polymers are disabled"); } }; - } // namespace Opm #endif diff --git a/opm/models/blackoil/blackoilpolymerparams.cpp b/opm/models/blackoil/blackoilpolymerparams.cpp new file mode 100644 index 000000000..219bfff89 --- /dev/null +++ b/opm/models/blackoil/blackoilpolymerparams.cpp @@ -0,0 +1,403 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + +#if HAVE_ECL_INPUT +#include +#include +#include +#include +#include +#include +#endif + +#include +#include +#include +#include + +namespace { + +#if FLOW_INSTANTIATE_FLOAT && HAVE_ECL_INPUT +std::vector> +doubleVecsToFloat(const std::vector>& input) +{ + std::vector> output; + output.reserve(input.size()); + for (std::size_t i = 0; i < input.size(); ++i) { + output.emplace_back(input[i].begin(), input[i].end()); + } + + return output; +} +#endif + +} + +namespace Opm { + +#if HAVE_ECL_INPUT +template +template +void BlackOilPolymerParams:: +initFromState(const 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)) { + 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)) { + 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)) { + 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)) { + throw std::runtime_error("Polymer molecular weight tracking is disabled at compile time, but the deck " + "contains the POLYMW keyword"); + } + } + + if constexpr (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; // polymer treatment is supposed to be disabled + } + + const auto& tableManager = eclState.getTableManager(); + + unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + 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); + setPlyrock(satRegionIdx, + plyrockTable.getDeadPoreVolumeColumn()[0], + plyrockTable.getResidualResistanceFactorColumn()[0], + plyrockTable.getRockDensityFactorColumn()[0], + static_cast(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(); + plyadsAdsorbedPolymer_[satRegionIdx].setXYContainers(c, ads); + } + } + 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 + const auto& plyviscTables = tableManager.getPlyviscTables(); + if (!plyviscTables.empty()) { + // different viscosity model is used for POLYMW + if constexpr (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(); + plyviscViscosityMultiplierTable_[pvtRegionIdx].setXYContainers(c, visc); + } + } + } + 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 + const auto& plymaxTables = tableManager.getPlymaxTables(); + const unsigned numMixRegions = plymaxTables.size(); + setNumMixRegions(numMixRegions, enablePolymerMolarWeight); + if (!plymaxTables.empty()) { + for (unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++mixRegionIdx) { + const auto& plymaxTable = plymaxTables.template getTable(mixRegionIdx); + plymaxMaxConcentration_[mixRegionIdx] = plymaxTable.getPolymerConcentrationColumn()[0]; + } + } + else { + throw std::runtime_error("PLYMAX must be specified in POLYMER runs\n"); + } + + if (!eclState.getTableManager().getPlmixparTable().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); + shearMultiplier[i] = shearMultiplier[i]; + } + 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 (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"); + } + } + } + + 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; + } + } + else { + throw std::runtime_error("PLYVMH keyword must be specified in POLYMW rus \n"); + } + + // 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 = doubleVecsToFloat(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); + } + } + + // 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 = doubleVecsToFloat(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); + } + } + + // 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 = doubleVecsToFloat(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); + } + } + } +} +#endif // HAVE_ECL_INPUT + +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:: +setPlyrock(unsigned satRegionIdx, + const Scalar& plyrockDeadPoreVolume, + const Scalar& plyrockResidualResistanceFactor, + const Scalar& plyrockRockDensityFactor, + const Scalar& plyrockAdsorbtionIndex, + const Scalar& plyrockMaxAdsorbtion) +{ + plyrockDeadPoreVolume_[satRegionIdx] = plyrockDeadPoreVolume; + plyrockResidualResistanceFactor_[satRegionIdx] = plyrockResidualResistanceFactor; + plyrockRockDensityFactor_[satRegionIdx] = plyrockRockDensityFactor; + plyrockAdsorbtionIndex_[satRegionIdx] = plyrockAdsorbtionIndex; + plyrockMaxAdsorbtion_[satRegionIdx] = plyrockMaxAdsorbtion; +} + +#define INSTANTIATE_TYPE(T) \ + template struct BlackOilPolymerParams; \ + template void BlackOilPolymerParams::initFromState(const EclipseState&); \ + template void BlackOilPolymerParams::initFromState(const EclipseState&); \ + template void BlackOilPolymerParams::initFromState(const EclipseState&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/models/blackoil/blackoilpolymerparams.hh b/opm/models/blackoil/blackoilpolymerparams.hpp similarity index 76% rename from opm/models/blackoil/blackoilpolymerparams.hh rename to opm/models/blackoil/blackoilpolymerparams.hpp index 9fdb63608..c67a3b133 100644 --- a/opm/models/blackoil/blackoilpolymerparams.hh +++ b/opm/models/blackoil/blackoilpolymerparams.hpp @@ -25,8 +25,8 @@ * * \brief Contains the parameters required to extend the black-oil model by polymer. */ -#ifndef EWOMS_BLACK_OIL_POLYMER_PARAMS_HH -#define EWOMS_BLACK_OIL_POLYMER_PARAMS_HH +#ifndef OPM_BLACK_OIL_POLYMER_PARAMS_HPP +#define OPM_BLACK_OIL_POLYMER_PARAMS_HPP #include #include @@ -36,6 +36,10 @@ namespace Opm { +#if HAVE_ECL_INPUT +class EclipseState; +#endif + //! \brief Struct holding the parameters for the BlackOilPolymerModule class. template struct BlackOilPolymerParams { @@ -44,35 +48,24 @@ struct BlackOilPolymerParams { enum AdsorptionBehaviour { Desorption = 1, NoDesorption = 2 }; +#if HAVE_ECL_INPUT + template + void initFromState(const EclipseState& eclState); +#endif + /*! * \brief Specify the number of satuation regions. * * This must be called before setting the PLYROCK and PLYADS of any region. */ - void setNumSatRegions(unsigned numRegions) - { - plyrockDeadPoreVolume_.resize(numRegions); - plyrockResidualResistanceFactor_.resize(numRegions); - plyrockRockDensityFactor_.resize(numRegions); - plyrockAdsorbtionIndex_.resize(numRegions); - plyrockMaxAdsorbtion_.resize(numRegions); - plyadsAdsorbedPolymer_.resize(numRegions); - } + 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) - { - plymaxMaxConcentration_.resize(numRegions); - plymixparToddLongstaff_.resize(numRegions); - - if (enablePolymerMolarWeight) { - plyvmhCoefficients_.resize(numRegions); - } - } + void setNumMixRegions(unsigned numRegions, bool enablePolymerMolarWeight); /*! * \brief Specify the polymer rock properties a single region. @@ -84,14 +77,7 @@ struct BlackOilPolymerParams { const Scalar& plyrockResidualResistanceFactor, const Scalar& plyrockRockDensityFactor, const Scalar& plyrockAdsorbtionIndex, - const Scalar& plyrockMaxAdsorbtion) - { - plyrockDeadPoreVolume_[satRegionIdx] = plyrockDeadPoreVolume; - plyrockResidualResistanceFactor_[satRegionIdx] = plyrockResidualResistanceFactor; - plyrockRockDensityFactor_[satRegionIdx] = plyrockRockDensityFactor; - plyrockAdsorbtionIndex_[satRegionIdx] = plyrockAdsorbtionIndex; - plyrockMaxAdsorbtion_[satRegionIdx] = plyrockMaxAdsorbtion; - } + const Scalar& plyrockMaxAdsorbtion); // a struct containing the constants to calculate polymer viscosity // based on Mark-Houwink equation and Huggins equation, the constants are provided @@ -132,4 +118,4 @@ struct BlackOilPolymerParams { } // namespace Opm -#endif +#endif // OPM_BLACK_OIL_POLYMER_PARAMS_HPP diff --git a/opm/models/blackoil/blackoilsolventmodules.hh b/opm/models/blackoil/blackoilsolventmodules.hh index b0052e506..4964f2043 100644 --- a/opm/models/blackoil/blackoilsolventmodules.hh +++ b/opm/models/blackoil/blackoilsolventmodules.hh @@ -28,28 +28,17 @@ #ifndef EWOMS_BLACK_OIL_SOLVENT_MODULE_HH #define EWOMS_BLACK_OIL_SOLVENT_MODULE_HH -#include "blackoilproperties.hh" - #include -#include -#include -#include -#include - #include -#if HAVE_ECL_INPUT -#include -#include -#include -#include -#include -#include -#include -#include -#include -#endif +#include +#include + +#include +#include + +#include #include @@ -96,261 +85,12 @@ class BlackOilSolventModule static constexpr bool blackoilConserveSurfaceVolume = getPropValue(); static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx; - public: -#if HAVE_ECL_INPUT - /*! - * \brief Initialize all internal data structures needed by the solvent module - */ - static void initFromState(const EclipseState& eclState, const Schedule& schedule) + //! \brief Set parameters. + static void setParams(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; // solvent treatment is supposed to be disabled - - params_.co2sol_ = eclState.runspec().co2Sol(); - params_.h2sol_ = eclState.runspec().h2Sol(); - - if (isCO2Sol() && isH2Sol()) { - throw std::runtime_error("CO2SOL and H2SOL can not be used together"); - } - - if (isCO2Sol() || isH2Sol() ) { - if (isCO2Sol()) { - params_.co2GasPvt_.initFromState(eclState, schedule); - params_.brineCo2Pvt_.initFromState(eclState, schedule); - } else { - params_.h2GasPvt_.initFromState(eclState, schedule); - params_.brineH2Pvt_.initFromState(eclState, schedule); - } - if (eclState.getSimulationConfig().hasDISGASW()) { - params_.rsSolw_active_ = true; - } - } else - 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 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 - 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"); - - // 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; - } - } + params_ = params; } -#endif /*! * \brief Specify the solvent PVT of a all PVT regions. diff --git a/opm/models/blackoil/blackoilsolventparams.cpp b/opm/models/blackoil/blackoilsolventparams.cpp new file mode 100644 index 000000000..606bddce7 --- /dev/null +++ b/opm/models/blackoil/blackoilsolventparams.cpp @@ -0,0 +1,334 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 + +#if HAVE_ECL_INPUT +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +namespace Opm { + +#if HAVE_ECL_INPUT +template +template +void BlackOilSolventParams:: +initFromState(const EclipseState& eclState, const Schedule& schedule) +{ + // 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)) { + 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)) { + 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; // solvent treatment is supposed to be disabled + } + + co2sol_ = eclState.runspec().co2Sol(); + h2sol_ = eclState.runspec().h2Sol(); + + if (co2sol_ && h2sol_) { + throw std::runtime_error("CO2SOL and H2SOL can not be used together"); + } + + if (co2sol_ || h2sol_) { + if (co2sol_) { + co2GasPvt_.initFromState(eclState, schedule); + brineCo2Pvt_.initFromState(eclState, schedule); + } else { + h2GasPvt_.initFromState(eclState, schedule); + brineH2Pvt_.initFromState(eclState, schedule); + } + if (eclState.getSimulationConfig().hasDISGASW()) { + rsSolw_active_ = true; + } + } else + 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(); + setNumSatRegions(numSatRegions); + for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++satRegionIdx) { + const auto& ssfnTable = ssfnTables.template getTable(satRegionIdx); + ssfnKrg_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(), + ssfnTable.getGasRelPermMultiplierColumn(), + /*sortInput=*/true); + ssfnKrs_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(), + 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"); + } + + 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); + + // 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()); + + 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(); + + 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; + } + } + + // 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); + } + } + // 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); + + // Copy data + const auto& sw = sorwmisTable.getWaterSaturationColumn(); + const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn(); + + 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; + } + } + + // 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); + + // Copy data + const auto& sw = sgcwmisTable.getWaterSaturationColumn(); + const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn(); + + 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; + } + + 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"); + + // 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."); + } + } + } + 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; + } + } +} +#endif + +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; +} + +#define INSTANTIATE_TYPE(T) \ + template struct BlackOilSolventParams; \ + template void BlackOilSolventParams::initFromState(const EclipseState&, const Schedule&); \ + template void BlackOilSolventParams::initFromState(const EclipseState&, const Schedule&); + +INSTANTIATE_TYPE(double) + +#if FLOW_INSTANTIATE_FLOAT +INSTANTIATE_TYPE(float) +#endif + +} // namespace Opm diff --git a/opm/models/blackoil/blackoilsolventparams.hh b/opm/models/blackoil/blackoilsolventparams.hpp similarity index 88% rename from opm/models/blackoil/blackoilsolventparams.hh rename to opm/models/blackoil/blackoilsolventparams.hpp index c96c91d38..3962bf63c 100644 --- a/opm/models/blackoil/blackoilsolventparams.hh +++ b/opm/models/blackoil/blackoilsolventparams.hpp @@ -25,8 +25,8 @@ * * \brief Contains the parameters required to extend the black-oil model by solvents. */ -#ifndef EWOMS_BLACK_OIL_SOLVENT_PARAMS_HH -#define EWOMS_BLACK_OIL_SOLVENT_PARAMS_HH +#ifndef OPM_BLACK_OIL_SOLVENT_PARAMS_HPP +#define OPM_BLACK_OIL_SOLVENT_PARAMS_HPP #include #include @@ -38,25 +38,32 @@ namespace Opm { +#if HAVE_ECL_INPUT +class EclipseState; +class Schedule; +#endif + //! \brief Struct holding the parameters for the BlackOilSolventModule class. template -struct BlackOilSolventParams { +struct BlackOilSolventParams +{ + using BrineCo2Pvt = ::Opm::BrineCo2Pvt; + using BrineH2Pvt = ::Opm::BrineH2Pvt; + using Co2GasPvt = ::Opm::Co2GasPvt; + using H2GasPvt = ::Opm::H2GasPvt; + using SolventPvt = ::Opm::SolventPvt; using TabulatedFunction = Tabulated1DFunction; - using SolventPvt = ::Opm::SolventPvt; - SolventPvt solventPvt_; +#if HAVE_ECL_INPUT + template + void initFromState(const EclipseState& eclState, const Schedule& schedule); +#endif - using Co2GasPvt = ::Opm::Co2GasPvt; - Co2GasPvt co2GasPvt_; - - using H2GasPvt = ::Opm::H2GasPvt; - H2GasPvt h2GasPvt_; - - using BrineCo2Pvt = ::Opm::BrineCo2Pvt; BrineCo2Pvt brineCo2Pvt_; - - using BrineH2Pvt = ::Opm::BrineH2Pvt; BrineH2Pvt brineH2Pvt_; + Co2GasPvt co2GasPvt_; + H2GasPvt h2GasPvt_; + SolventPvt solventPvt_; std::vector ssfnKrg_; // the krg(Fs) column of the SSFN table std::vector ssfnKrs_; // the krs(Fs) column of the SSFN table @@ -82,11 +89,7 @@ struct BlackOilSolventParams { * * This must be called before setting the SSFN of any region. */ - void setNumSatRegions(unsigned numRegions) - { - ssfnKrg_.resize(numRegions); - ssfnKrs_.resize(numRegions); - } + void setNumSatRegions(unsigned numRegions); /*! * \brief Specify miscible relative permeability multipliers of a single region. @@ -95,13 +98,9 @@ struct BlackOilSolventParams { */ void setMsfn(unsigned satRegionIdx, const TabulatedFunction& msfnKrsg, - const TabulatedFunction& msfnKro) - { - msfnKrsg_[satRegionIdx] = msfnKrsg; - msfnKro_[satRegionIdx] = msfnKro; - } + const TabulatedFunction& msfnKro); }; } // namespace Opm -#endif +#endif // OPM_BLACK_OIL_SOLVENT_PARAMS_HPP diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index d750edb0d..e64120946 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -274,14 +274,34 @@ public: this->model().addOutputModule(new VtkTracerModule(simulator)); // Tell the black-oil extensions to initialize their internal data structures const auto& vanguard = simulator.vanguard(); - 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()); - DispersionModule::initFromState(vanguard.eclState()); + + BlackOilBrineParams brineParams; + brineParams.template initFromState(vanguard.eclState()); + BrineModule::setParams(std::move(brineParams)); + DiffusionModule::initFromState(vanguard.eclState()); + DispersionModule::initFromState(vanguard.eclState()); + + BlackOilExtboParams extboParams; + extboParams.template initFromState(vanguard.eclState()); + ExtboModule::setParams(std::move(extboParams)); + + BlackOilFoamParams foamParams; + foamParams.template initFromState(vanguard.eclState()); + FoamModule::setParams(std::move(foamParams)); + + BlackOilMICPParams micpParams; + micpParams.template initFromState(vanguard.eclState()); + MICPModule::setParams(std::move(micpParams)); + + BlackOilPolymerParams polymerParams; + polymerParams.template initFromState(vanguard.eclState()); + PolymerModule::setParams(std::move(polymerParams)); + + BlackOilSolventParams solventParams; + solventParams.template initFromState(vanguard.eclState(), vanguard.schedule()); + SolventModule::setParams(std::move(solventParams)); // create the ECL writer eclWriter_ = std::make_unique(simulator);