Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change how the blackoilmodules are setup #725

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
72 changes: 4 additions & 68 deletions opm/models/blackoil/blackoilbrinemodules.hh
Original file line number Diff line number Diff line change
Expand Up @@ -32,21 +32,13 @@

#include <opm/models/blackoil/blackoilbrineparams.hh>

#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PvtwsaltTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PermfactTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SaltSolubilityTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
#endif

#include <dune/common/fvector.hh>

#include <string>
#include <math.h>

namespace Opm {

/*!
* \ingroup BlackOil
* \brief Contains the high level supplements required to extend the black oil
Expand Down Expand Up @@ -84,67 +76,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 the parameters.
static void setParams(BlackOilBrineParams<Scalar> 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<PermfactTable>(i);
params_.permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn());
}

const TableContainer& saltsolTables = tableManager.getSaltsolTables();
if (!saltsolTables.empty()) {
params_.saltsolTable_.resize(numPvtRegions);
assert(numPvtRegions == saltsolTables.size());
for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
const SaltsolTable& saltsolTable = saltsolTables.getTable<SaltsolTable>(pvtRegionIdx );
params_.saltsolTable_[pvtRegionIdx] = saltsolTable.getSaltsolColumn().front();
}
}
}
params_ = std::move(params);
}
#endif

/*!
* \brief Register all run-time parameters for the black-oil brine module.
Expand Down
173 changes: 3 additions & 170 deletions opm/models/blackoil/blackoilextbomodules.hh
Original file line number Diff line number Diff line change
Expand Up @@ -36,20 +36,6 @@

#include <opm/models/blackoil/blackoilextboparams.hh>

//#include <opm/models/io/vtkBlackOilExtboModule.hh> //TODO: Missing ...

#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SsfnTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/Sof2Table.hpp>
#include <opm/input/eclipse/EclipseState/Tables/MsfnTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PmiscTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/MiscTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SorwmisTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SgcwmisTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
#endif

#include <dune/common/fvector.hh>

#include <cmath>
Expand Down Expand Up @@ -96,164 +82,11 @@ class BlackOilExtboModule
static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();

public:
#if HAVE_ECL_INPUT
/*!
* \brief Initialize all internal data structures needed by the solvent module
*/
static void initFromState(const EclipseState& eclState)
//! \brief Set the parameters.
static void setParams(BlackOilExtboParams<Scalar> params)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider taking rvalue ref argument here. As written, we will have minimal overhead (one extra move) when called with a temporary, as it is used downstream, but if someone calls this with an lvalue it will incur a copy. The copy can be avoided if the user remembers to std::move() the argument in the call to this function, but making the argument an rvalue ref will force it, avoiding any chance of copying. (If you actually need to keep the original at the call site, you would need to make a temporary copy there.)

{
// 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; //<false>: 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<Scalar> oilCmp(saturatedTable.numRows(), -4.0e-9); //Default values used in [*]
std::vector<Scalar> gasCmp(saturatedTable.numRows(), -0.08); //-------------"-------------
params_.zLim_[regionIdx] = 0.7; //-------------"-------------
std::vector<Scalar> 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_ = std::move(params);
}
#endif

/*!
* \brief Register all run-time parameters for the black-oil solvent module.
Expand Down
97 changes: 3 additions & 94 deletions opm/models/blackoil/blackoilfoammodules.hh
Original file line number Diff line number Diff line change
Expand Up @@ -34,12 +34,6 @@

#include <opm/models/blackoil/blackoilfoamparams.hh>

#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/FoamadsTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/FoammobTable.hpp>
#endif

#include <dune/common/fvector.hh>

#include <string>
Expand Down Expand Up @@ -83,96 +77,11 @@ class BlackOilFoamModule
static constexpr unsigned numPhases = FluidSystem::numPhases;

public:
#if HAVE_ECL_INPUT
/*!
* \brief Initialize all internal data structures needed by the foam module
*/
static void initFromState(const EclipseState& eclState)
//! \brief Set the parameters.
static void setParams(BlackOilFoamParams<Scalar> 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
}

// Check that only implemented options are used.
// We only support the default values of FOAMOPTS (GAS, TAB).
if (eclState.getInitConfig().getFoamConfig().getTransportPhase() != Phase::GAS) {
throw std::runtime_error("In FOAMOPTS, only GAS is allowed for the transport phase.");
}
if (eclState.getInitConfig().getFoamConfig().getMobilityModel() != FoamConfig::MobilityModel::TAB) {
throw std::runtime_error("In 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<Scalar>::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<FoamadsTable>(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<FoammobTable>(pvtReg);
const auto& conc = foammobTable.getFoamConcentrationColumn();
const auto& mobMult = foammobTable.getMobilityMultiplierColumn();
params_.gasMobilityMultiplierTable_[pvtReg].setXYContainers(conc, mobMult);
}
params_ = std::move(params);
}
#endif

/*!
* \brief Register all run-time parameters for the black-oil foam module.
Expand Down