OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5574 from akva2/blackoilparams_tu

Browse Source 
BlackoilXXParams: introduce translation units
This commit is contained in:
Bård Skaflestad 2024-09-04 11:31:19 +02:00 committed by GitHub
commit 44d22a05ba
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
20 changed files with 1536 additions and 1074 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
CMakeLists_files.cmake
View File

@ -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

90
opm/models/blackoil/blackoilbrinemodules.hh
View File

@ -28,24 +28,16 @@
#ifndef EWOMS_BLACK_OIL_BRINE_MODULE_HH
#define EWOMS_BLACK_OIL_BRINE_MODULE_HH
#include "blackoilproperties.hh"
#include <opm/models/blackoil/blackoilbrineparams.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/blackoil/blackoilbrineparams.hh>
#include <opm/models/discretization/common/fvbaseproperties.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/PcfactTable.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 <opm/models/utils/basicproperties.hh>
#include <dune/common/fvector.hh>
#include <string>
#include <math.h>
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<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);
params_.saltdenTable_.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_.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<PcfactTable>(i);
params_.pcfactTable_[i].setXYContainers(pcfactTable.getPorosityChangeColumn(), pcfactTable.getPcMultiplierColumn());
}
}
}
params_ = params;
}
#endif
/*!
* \brief Register all run-time parameters for the black-oil brine module.

134
opm/models/blackoil/blackoilbrineparams.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <config.h>
#include <opm/models/blackoil/blackoilbrineparams.hpp>
#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/PcfactTable.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 <cassert>
#include <cstddef>
#include <stdexcept>
namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar>
template<bool enableBrine, bool enableSaltPrecipitation>
void BlackOilBrineParams<Scalar>::
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<PermfactTable>(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<SaltsolTable>(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<PcfactTable>(i);
pcfactTable_[i].setXYContainers(pcfactTable.getPorosityChangeColumn(), pcfactTable.getPcMultiplierColumn());
}
}
}
}
#endif
#define INSTANTIATE_TYPE(T) \
template struct BlackOilBrineParams<T>; \
template void BlackOilBrineParams<T>::initFromState<false,false>(const EclipseState&); \
template void BlackOilBrineParams<T>::initFromState<true,false>(const EclipseState&); \
template void BlackOilBrineParams<T>::initFromState<false,true>(const EclipseState&); \
template void BlackOilBrineParams<T>::initFromState<true,true>(const EclipseState&);
INSTANTIATE_TYPE(double)
#if FLOW_INSTANTIATE_FLOAT
INSTANTIATE_TYPE(float)
#endif
} // namespace Opm

18
opm/models/blackoil/blackoilbrineparams.hh → opm/models/blackoil/blackoilbrineparams.hpp
View File

@ -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 <opm/material/common/Tabulated1DFunction.hpp>
@ -34,9 +34,19 @@
namespace Opm {
#if HAVE_ECL_INPUT
class EclipseState;
#endif
//! \brief Struct holding the parameters for the BlackoilBrineModule class.
template<class Scalar>
struct BlackOilBrineParams {
struct BlackOilBrineParams
{
#if HAVE_ECL_INPUT
template<bool enableBrine, bool enableSaltPrecipitation>
void initFromState(const EclipseState& eclState);
#endif
using TabulatedFunction = Tabulated1DFunction<Scalar>;
std::vector<TabulatedFunction> bdensityTable_;
@ -49,4 +59,4 @@ struct BlackOilBrineParams {
} // namespace Opm
#endif
#endif // OPM_BLACK_OIL_BRINE_PARAMS_HPP

179
opm/models/blackoil/blackoilextbomodules.hh
View File

@ -32,30 +32,20 @@
#ifndef EWOMS_BLACK_OIL_EXTBO_MODULE_HH
#define EWOMS_BLACK_OIL_EXTBO_MODULE_HH
#include "blackoilproperties.hh"
#include <opm/models/blackoil/blackoilextboparams.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/blackoil/blackoilextboparams.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/utils/basicproperties.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>
#include <stdexcept>
#include <string>
#include <vector>
namespace Opm {
@ -96,164 +86,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 parameters.
static void setParams(BlackOilExtboParams<Scalar>&& 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; //<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_ = params;
}
#endif
/*!
* \brief Register all run-time parameters for the black-oil solvent module.

220
opm/models/blackoil/blackoilextboparams.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <config.h>
#include <opm/models/blackoil/blackoilextboparams.hpp>
#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 <cstddef>
#include <stdexcept>
namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar>
template<bool enableExtbo>
void BlackOilExtboParams<Scalar>::
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; //<false>: 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<Scalar> oilCmp(saturatedTable.numRows(), -4.0e-9); //Default values used in [*]
std::vector<Scalar> gasCmp(saturatedTable.numRows(), -0.08); //-------------"-------------
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;
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<T>; \
template void BlackOilExtboParams<T>::initFromState<false>(const EclipseState&); \
template void BlackOilExtboParams<T>::initFromState<true>(const EclipseState&);
INSTANTIATE_TYPE(double)
#if FLOW_INSTANTIATE_FLOAT
INSTANTIATE_TYPE(float)
#endif
} // namespace Opm

19
opm/models/blackoil/blackoilextboparams.hh → opm/models/blackoil/blackoilextboparams.hpp
View File

@ -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 <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/common/UniformXTabulated2DFunction.hpp>
@ -39,12 +39,23 @@
namespace Opm {
#if HAVE_ECL_INPUT
class EclipseState;
#endif
//! \brief Struct holding the parameters for the BlackoilExtboModule class.
template<class Scalar>
struct BlackOilExtboParams {
struct BlackOilExtboParams
{
using TabulatedFunction = Tabulated1DFunction<Scalar>;
using Tabulated2DFunction = UniformXTabulated2DFunction<Scalar>;
#if HAVE_ECL_INPUT
template<bool enableExtbo>
void initFromState(const EclipseState& eclState);
#endif
std::vector<Tabulated2DFunction> X_;
std::vector<Tabulated2DFunction> Y_;
std::vector<Tabulated2DFunction> PBUB_RS_;
@ -65,4 +76,4 @@ struct BlackOilExtboParams {
} // namespace Opm
#endif
#endif // OPM_BLACK_OIL_EXTBO_PARAMS_HPP

103
opm/models/blackoil/blackoilfoammodules.hh
View File

@ -28,23 +28,24 @@
#ifndef EWOMS_BLACK_OIL_FOAM_MODULE_HH
#define EWOMS_BLACK_OIL_FOAM_MODULE_HH
#include "blackoilproperties.hh"
#include <dune/common/fvector.hh>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/input/eclipse/EclipseState/Phase.hpp>
#include <opm/models/blackoil/blackoilfoamparams.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/discretization/common/fvbaseparameters.hh>
#include <opm/models/discretization/common/fvbaseproperties.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 <opm/models/blackoil/blackoilfoamparams.hh>
#include <opm/models/discretization/common/fvbaseparameters.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <string>
namespace Opm {
@ -87,93 +88,11 @@ class BlackOilFoamModule
enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
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<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
}
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<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_ = params;
}
#endif
/*!
* \brief Register all run-time parameters for the black-oil foam module.

150
opm/models/blackoil/blackoilfoamparams.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <config.h>
#include <opm/models/blackoil/blackoilfoamparams.hpp>
#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 <cstddef>
#include <stdexcept>
namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar>
template<bool enableFoam>
void BlackOilFoamParams<Scalar>::
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<FoamadsTable>(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<FoammobTable>(pvtReg);
const auto& conc = foammobTable.getFoamConcentrationColumn();
const auto& mobMult = foammobTable.getMobilityMultiplierColumn();
gasMobilityMultiplierTable_[pvtReg].setXYContainers(conc, mobMult);
}
}
#endif
template<class Scalar>
void BlackOilFoamParams<Scalar>::setNumSatRegions(unsigned numRegions)
{
foamCoefficients_.resize(numRegions);
foamRockDensity_.resize(numRegions);
foamAllowDesorption_.resize(numRegions);
adsorbedFoamTable_.resize(numRegions);
}
#define INSTANTIATE_TYPE(T) \
template struct BlackOilFoamParams<T>; \
template void BlackOilFoamParams<T>::initFromState<false>(const EclipseState&); \
template void BlackOilFoamParams<T>::initFromState<true>(const EclipseState&);
INSTANTIATE_TYPE(double)
#if FLOW_INSTANTIATE_FLOAT
INSTANTIATE_TYPE(float)
#endif
} // namespace Opm

29
opm/models/blackoil/blackoilfoamparams.hh → opm/models/blackoil/blackoilfoamparams.hpp
View File

@ -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 <opm/input/eclipse/EclipseState/Phase.hpp>
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <vector>
namespace Opm {
enum class Phase;
#if HAVE_ECL_INPUT
class EclipseState;
#endif
//! \brief Struct holding the parameters for the BlackoilFoamModule class.
template<class Scalar>
struct BlackOilFoamParams {
struct BlackOilFoamParams
{
using TabulatedFunction = Tabulated1DFunction<Scalar>;
#if HAVE_ECL_INPUT
template<bool enableFoam>
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

110
opm/models/blackoil/blackoilmicpmodules.hh
View File

@ -28,22 +28,15 @@
#ifndef EWOMS_BLACK_OIL_MICP_MODULE_HH
#define EWOMS_BLACK_OIL_MICP_MODULE_HH
#include "blackoilproperties.hh"
#include <opm/models/blackoil/blackoilmicpparams.hh>
#include <opm/models/io/vtkblackoilmicpmodule.hh>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/MICPpara.hpp>
#endif
#include <dune/common/fvector.hh>
#include <cmath>
#include <opm/models/blackoil/blackoilmicpparams.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/io/vtkblackoilmicpmodule.hh>
#include <cstddef>
#include <stdexcept>
#include <string>
namespace Opm {
/*!
@ -86,56 +79,11 @@ class BlackOilMICPModule
static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
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<Scalar>&& 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<Scalar, float>) {
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.
*/

104
opm/models/blackoil/blackoilmicpparams.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <config.h>
#include <opm/models/blackoil/blackoilmicpparams.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/MICPpara.hpp>
#endif
#include <algorithm>
#include <stdexcept>
#include <type_traits>
namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar>
template<bool enableMICP>
void BlackOilMICPParams<Scalar>::
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<Scalar, float>) {
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<T>; \
template void BlackOilMICPParams<T>::initFromState<false>(const EclipseState&); \
template void BlackOilMICPParams<T>::initFromState<true>(const EclipseState&);
INSTANTIATE_TYPE(double)
#if FLOW_INSTANTIATE_FLOAT
INSTANTIATE_TYPE(float)
#endif
} // namespace Opm

18
opm/models/blackoil/blackoilmicpparams.hh → opm/models/blackoil/blackoilmicpparams.hpp
View File

@ -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 <vector>
namespace Opm {
#if HAVE_ECL_INPUT
class EclipseState;
#endif
//! \brief Struct holding the parameters for the BlackOilMICPModule class.
template<class Scalar>
struct BlackOilMICPParams {
struct BlackOilMICPParams
{
#if HAVE_ECL_INPUT
template<bool enableMICP>
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

274
opm/models/blackoil/blackoilpolymermodules.hh
View File

@ -28,23 +28,16 @@
#ifndef EWOMS_BLACK_OIL_POLYMER_MODULE_HH
#define EWOMS_BLACK_OIL_POLYMER_MODULE_HH
#include "blackoilproperties.hh"
#include <opm/models/blackoil/blackoilpolymerparams.hh>
#include <opm/models/io/vtkblackoilpolymermodule.hh>
#include <dune/common/fvector.hh>
#include <opm/common/OpmLog/OpmLog.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyadsTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlymaxTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyrockTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyshlogTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyviscTable.hpp>
#endif
#include <opm/models/blackoil/blackoilpolymerparams.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <dune/common/fvector.hh>
#include <opm/models/io/vtkblackoilpolymermodule.hh>
#include <opm/models/utils/propertysystem.hh>
#include <cmath>
#include <stdexcept>
@ -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<Scalar>&& 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<PlyrockTable>(satRegionIdx);
params_.setPlyrock(satRegionIdx,
plyrockTable.getDeadPoreVolumeColumn()[0],
plyrockTable.getResidualResistanceFactorColumn()[0],
plyrockTable.getRockDensityFactorColumn()[0],
static_cast<typename BlackOilPolymerParams<Scalar>::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<PlyadsTable>(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<PlyviscTable>(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<PlymaxTable>(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<PlyshlogTable>(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<double>& throughput = plymwinjtable.getThroughputs();
const std::vector<double>& watervelocity = plymwinjtable.getVelocities();
const std::vector<std::vector<double>>& 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<double>& throughput = skprwattable.getThroughputs();
const std::vector<double>& watervelocity = skprwattable.getVelocities();
const std::vector<std::vector<double>>& 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<double>& throughput = skprpolytable.getThroughputs();
const std::vector<double>& watervelocity = skprpolytable.getVelocities();
const std::vector<std::vector<double>>& skinpressure = skprpolytable.getSkinPressures();
const double refPolymerConcentration = skprpolytable.referenceConcentration();
typename BlackOilPolymerParams<Scalar>::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

403
opm/models/blackoil/blackoilpolymerparams.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <config.h>
#include <opm/models/blackoil/blackoilpolymerparams.hpp>
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyadsTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlymaxTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyrockTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyshlogTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PlyviscTable.hpp>
#endif
#include <cassert>
#include <cstddef>
#include <stdexcept>
#include <type_traits>
namespace {
#if FLOW_INSTANTIATE_FLOAT && HAVE_ECL_INPUT
std::vector<std::vector<float>>
doubleVecsToFloat(const std::vector<std::vector<double>>& input)
{
std::vector<std::vector<float>> 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<class Scalar>
template<bool enablePolymer, bool enablePolymerMolarWeight>
void BlackOilPolymerParams<Scalar>::
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<PlyrockTable>(satRegionIdx);
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<PlyadsTable>(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<PlyviscTable>(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<PlymaxTable>(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<PlyshlogTable>(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<double>& throughput = plymwinjtable.getThroughputs();
const std::vector<double>& watervelocity = plymwinjtable.getVelocities();
const std::vector<std::vector<double>>& molecularweight = plymwinjtable.getMoleWeights();
if constexpr (std::is_same_v<Scalar, float>) {
const std::vector<Scalar> tp(throughput.begin(), throughput.end());
const std::vector<Scalar> 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<double>& throughput = skprwattable.getThroughputs();
const std::vector<double>& watervelocity = skprwattable.getVelocities();
const std::vector<std::vector<double>>& skinpressure = skprwattable.getSkinPressures();
if constexpr (std::is_same_v<Scalar, float>) {
const std::vector<Scalar> tp(throughput.begin(), throughput.end());
const std::vector<Scalar> 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<double>& throughput = skprpolytable.getThroughputs();
const std::vector<double>& watervelocity = skprpolytable.getVelocities();
const std::vector<std::vector<double>>& skinpressure = skprpolytable.getSkinPressures();
const double refPolymerConcentration = skprpolytable.referenceConcentration();
if constexpr (std::is_same_v<Scalar, float>) {
const std::vector<Scalar> tp(throughput.begin(), throughput.end());
const std::vector<Scalar> 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<class Scalar>
void BlackOilPolymerParams<Scalar>::
setNumSatRegions(unsigned numRegions)
{
plyrockDeadPoreVolume_.resize(numRegions);
plyrockResidualResistanceFactor_.resize(numRegions);
plyrockRockDensityFactor_.resize(numRegions);
plyrockAdsorbtionIndex_.resize(numRegions);
plyrockMaxAdsorbtion_.resize(numRegions);
plyadsAdsorbedPolymer_.resize(numRegions);
}
template<class Scalar>
void BlackOilPolymerParams<Scalar>::
setNumMixRegions(unsigned numRegions, bool enablePolymerMolarWeight)
{
plymaxMaxConcentration_.resize(numRegions);
plymixparToddLongstaff_.resize(numRegions);
if (enablePolymerMolarWeight) {
plyvmhCoefficients_.resize(numRegions);
}
}
template<class Scalar>
void BlackOilPolymerParams<Scalar>::
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<T>; \
template void BlackOilPolymerParams<T>::initFromState<false,false>(const EclipseState&); \
template void BlackOilPolymerParams<T>::initFromState<true,false>(const EclipseState&); \
template void BlackOilPolymerParams<T>::initFromState<true,true>(const EclipseState&);
INSTANTIATE_TYPE(double)
#if FLOW_INSTANTIATE_FLOAT
INSTANTIATE_TYPE(float)
#endif
} // namespace Opm

44
opm/models/blackoil/blackoilpolymerparams.hh → opm/models/blackoil/blackoilpolymerparams.hpp
View File

@ -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 <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/common/IntervalTabulated2DFunction.hpp>
@ -36,6 +36,10 @@
namespace Opm {
#if HAVE_ECL_INPUT
class EclipseState;
#endif
//! \brief Struct holding the parameters for the BlackOilPolymerModule class.
template<class Scalar>
struct BlackOilPolymerParams {
@ -44,35 +48,24 @@ struct BlackOilPolymerParams {
enum AdsorptionBehaviour { Desorption = 1, NoDesorption = 2 };
#if HAVE_ECL_INPUT
template<bool enablePolymer, bool enablePolymerMolarWeight>
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

280
opm/models/blackoil/blackoilsolventmodules.hh
View File

@ -28,28 +28,17 @@
#ifndef EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
#define EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
#include "blackoilproperties.hh"
#include <opm/common/Exceptions.hpp>
#include <opm/models/blackoil/blackoilsolventparams.hh>
#include <opm/models/io/vtkblackoilsolventmodule.hh>
#include <opm/models/common/multiphasebaseparameters.hh>
#include <opm/models/common/quantitycallbacks.hh>
#include <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
#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 <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/blackoil/blackoilsolventparams.hpp>
#include <opm/models/common/multiphasebaseparameters.hh>
#include <opm/models/common/quantitycallbacks.hh>
#include <opm/models/io/vtkblackoilsolventmodule.hh>
#include <opm/material/common/Valgrind.hpp>
@ -96,261 +85,12 @@ class BlackOilSolventModule
static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
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<Scalar>&& 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<SsfnTable>(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<Sof2Table>(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<MiscTable>(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<PmiscTable>(regionIdx);
// Copy data
const auto& po = pmiscTable.getOilPhasePressureColumn();
const auto& pmisc = pmiscTable.getMiscibilityColumn();
params_.pmisc_[regionIdx].setXYContainers(po, pmisc);
}
}
else {
std::vector<double> x = {0.0,1.0e20};
std::vector<double> 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<MsfnTable>(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<double> x = {0.0,1.0};
std::vector<double> 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<SorwmisTable>(regionIdx);
// Copy data
const auto& sw = sorwmisTable.getWaterSaturationColumn();
const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn();
params_.sorwmis_[regionIdx].setXYContainers(sw, sorwmis);
}
}
else {
// default
std::vector<double> x = {0.0,1.0};
std::vector<double> 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<SgcwmisTable>(regionIdx);
// Copy data
const auto& sw = sgcwmisTable.getWaterSaturationColumn();
const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn();
params_.sgcwmis_[regionIdx].setXYContainers(sw, sgcwmis);
}
}
else {
// default
std::vector<double> x = {0.0,1.0};
std::vector<double> 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<TlpmixpaTable>(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<double> x = {0.0,1.0e20};
std::vector<double> 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.

334
opm/models/blackoil/blackoilsolventparams.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <config.h>
#include <opm/models/blackoil/blackoilsolventparams.hpp>
#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
namespace Opm {
#if HAVE_ECL_INPUT
template<class Scalar>
template<bool enableSolvent>
void BlackOilSolventParams<Scalar>::
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<SsfnTable>(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<Sof2Table>(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<MiscTable>(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<PmiscTable>(regionIdx);
// Copy data
const auto& po = pmiscTable.getOilPhasePressureColumn();
const auto& pmisc = pmiscTable.getMiscibilityColumn();
pmisc_[regionIdx].setXYContainers(po, pmisc);
}
}
else {
std::vector<double> x = {0.0,1.0e20};
std::vector<double> 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<MsfnTable>(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<double> x = {0.0,1.0};
std::vector<double> 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<SorwmisTable>(regionIdx);
// Copy data
const auto& sw = sorwmisTable.getWaterSaturationColumn();
const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn();
sorwmis_[regionIdx].setXYContainers(sw, sorwmis);
}
}
else {
// default
std::vector<double> x = {0.0,1.0};
std::vector<double> 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<SgcwmisTable>(regionIdx);
// Copy data
const auto& sw = sgcwmisTable.getWaterSaturationColumn();
const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn();
sgcwmis_[regionIdx].setXYContainers(sw, sgcwmis);
}
}
else {
// default
std::vector<double> x = {0.0,1.0};
std::vector<double> 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<TlpmixpaTable>(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<double> x = {0.0,1.0e20};
std::vector<double> y = {1.0,1.0};
TabulatedFunction ones = TabulatedFunction(2, x, y);
for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx)
tlPMixTable_[regionIdx] = ones;
}
}
}
#endif
template<class Scalar>
void BlackOilSolventParams<Scalar>::
setNumSatRegions(unsigned numRegions)
{
ssfnKrg_.resize(numRegions);
ssfnKrs_.resize(numRegions);
}
template<class Scalar>
void BlackOilSolventParams<Scalar>::
setMsfn(unsigned satRegionIdx,
const TabulatedFunction& msfnKrsg,
const TabulatedFunction& msfnKro)
{
msfnKrsg_[satRegionIdx] = msfnKrsg;
msfnKro_[satRegionIdx] = msfnKro;
}
#define INSTANTIATE_TYPE(T) \
template struct BlackOilSolventParams<T>; \
template void BlackOilSolventParams<T>::initFromState<false>(const EclipseState&, const Schedule&); \
template void BlackOilSolventParams<T>::initFromState<true>(const EclipseState&, const Schedule&);
INSTANTIATE_TYPE(double)
#if FLOW_INSTANTIATE_FLOAT
INSTANTIATE_TYPE(float)
#endif
} // namespace Opm

49
opm/models/blackoil/blackoilsolventparams.hh → opm/models/blackoil/blackoilsolventparams.hpp
View File

@ -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 <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
#include <opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp>
@ -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<class Scalar>
struct BlackOilSolventParams {
struct BlackOilSolventParams
{
using BrineCo2Pvt = ::Opm::BrineCo2Pvt<Scalar>;
using BrineH2Pvt = ::Opm::BrineH2Pvt<Scalar>;
using Co2GasPvt = ::Opm::Co2GasPvt<Scalar>;
using H2GasPvt = ::Opm::H2GasPvt<Scalar>;
using SolventPvt = ::Opm::SolventPvt<Scalar>;
using TabulatedFunction = Tabulated1DFunction<Scalar>;
using SolventPvt = ::Opm::SolventPvt<Scalar>;
SolventPvt solventPvt_;
#if HAVE_ECL_INPUT
template<bool enableSolvent>
void initFromState(const EclipseState& eclState, const Schedule& schedule);
#endif
using Co2GasPvt = ::Opm::Co2GasPvt<Scalar>;
Co2GasPvt co2GasPvt_;
using H2GasPvt = ::Opm::H2GasPvt<Scalar>;
H2GasPvt h2GasPvt_;
using BrineCo2Pvt = ::Opm::BrineCo2Pvt<Scalar>;
BrineCo2Pvt brineCo2Pvt_;
using BrineH2Pvt = ::Opm::BrineH2Pvt<Scalar>;
BrineH2Pvt brineH2Pvt_;
Co2GasPvt co2GasPvt_;
H2GasPvt h2GasPvt_;
SolventPvt solventPvt_;
std::vector<TabulatedFunction> ssfnKrg_; // the krg(Fs) column of the SSFN table
std::vector<TabulatedFunction> 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

34
opm/simulators/flow/FlowProblem.hpp
View File

@ -274,14 +274,34 @@ public:
this->model().addOutputModule(new VtkTracerModule<TypeTag>(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<Scalar> brineParams;
brineParams.template initFromState<enableBrine,
enableSaltPrecipitation>(vanguard.eclState());
BrineModule::setParams(std::move(brineParams));
DiffusionModule::initFromState(vanguard.eclState());
DispersionModule::initFromState(vanguard.eclState());
BlackOilExtboParams<Scalar> extboParams;
extboParams.template initFromState<enableExtbo>(vanguard.eclState());
ExtboModule::setParams(std::move(extboParams));
BlackOilFoamParams<Scalar> foamParams;
foamParams.template initFromState<enableFoam>(vanguard.eclState());
FoamModule::setParams(std::move(foamParams));
BlackOilMICPParams<Scalar> micpParams;
micpParams.template initFromState<enableMICP>(vanguard.eclState());
MICPModule::setParams(std::move(micpParams));
BlackOilPolymerParams<Scalar> polymerParams;
polymerParams.template initFromState<enablePolymer, enablePolymerMolarWeight>(vanguard.eclState());
PolymerModule::setParams(std::move(polymerParams));
BlackOilSolventParams<Scalar> solventParams;
solventParams.template initFromState<enableSolvent>(vanguard.eclState(), vanguard.schedule());
SolventModule::setParams(std::move(solventParams));
// create the ECL writer
eclWriter_ = std::make_unique<EclWriterType>(simulator);