mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-25 16:51:00 -06:00
1121 lines
45 KiB
C++
1121 lines
45 KiB
C++
/*
|
|
Copyright 2013 SINTEF ICT, Applied Mathematics.
|
|
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services.
|
|
Copyright 2015 NTNU.
|
|
|
|
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 3 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/>.
|
|
*/
|
|
|
|
#include <config.h>
|
|
|
|
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
|
|
#include <opm/autodiff/BlackoilModelEnums.hpp>
|
|
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
|
|
|
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
|
#include <opm/core/props/BlackoilPhases.hpp>
|
|
#include <opm/parser/eclipse/Units/Units.hpp>
|
|
#include <opm/grid/utility/extractPvtTableIndex.hpp>
|
|
|
|
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
|
|
#include <opm/material/common/Exceptions.hpp>
|
|
|
|
#include <opm/parser/eclipse/Deck/Deck.hpp>
|
|
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
|
|
|
#include <opm/common/ErrorMacros.hpp>
|
|
|
|
namespace Opm
|
|
{
|
|
// Making these typedef to make the code more readable.
|
|
typedef BlackoilPropsAdFromDeck::ADB ADB;
|
|
typedef BlackoilPropsAdFromDeck::V V;
|
|
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
|
|
|
|
/// Constructor wrapping an opm-core black oil interface.
|
|
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const Opm::Deck& deck,
|
|
const Opm::EclipseState& eclState,
|
|
std::shared_ptr<MaterialLawManager> materialLawManager,
|
|
const UnstructuredGrid& grid,
|
|
const bool init_rock)
|
|
{
|
|
init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims,
|
|
init_rock);
|
|
}
|
|
|
|
#ifdef HAVE_OPM_GRID
|
|
/// Constructor wrapping an opm-core black oil interface.
|
|
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const Opm::Deck& deck,
|
|
const Opm::EclipseState& eclState,
|
|
const Dune::CpGrid& grid,
|
|
const bool init_rock )
|
|
{
|
|
auto materialLawManager = std::make_shared<MaterialLawManager>();
|
|
unsigned number_of_cells = grid.size(0);
|
|
std::vector<int> compressedToCartesianIdx(number_of_cells);
|
|
for (unsigned cellIdx = 0; cellIdx < number_of_cells; ++cellIdx) {
|
|
compressedToCartesianIdx[cellIdx] = grid.globalCell()[cellIdx];
|
|
}
|
|
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
|
|
init(deck, eclState, materialLawManager, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
|
|
static_cast<const int*>(&grid.logicalCartesianSize()[0]),
|
|
init_rock);
|
|
}
|
|
#endif
|
|
|
|
/// Constructor wrapping an opm-core black oil interface.
|
|
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const Opm::Deck& deck,
|
|
const Opm::EclipseState& eclState,
|
|
const UnstructuredGrid& grid,
|
|
const bool init_rock)
|
|
{
|
|
auto materialLawManager = std::make_shared<MaterialLawManager>();
|
|
std::vector<int> compressedToCartesianIdx(grid.number_of_cells);
|
|
for (int cellIdx = 0; cellIdx < grid.number_of_cells; ++cellIdx) {
|
|
if (grid.global_cell) {
|
|
compressedToCartesianIdx[cellIdx] = grid.global_cell[cellIdx];
|
|
}
|
|
else {
|
|
compressedToCartesianIdx[cellIdx] = cellIdx;
|
|
}
|
|
}
|
|
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
|
|
init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims,
|
|
init_rock);
|
|
}
|
|
|
|
#ifdef HAVE_OPM_GRID
|
|
/// Constructor wrapping an opm-core black oil interface.
|
|
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const Opm::Deck& deck,
|
|
const Opm::EclipseState& eclState,
|
|
std::shared_ptr<MaterialLawManager> materialLawManager,
|
|
const Dune::CpGrid& grid,
|
|
const bool init_rock )
|
|
{
|
|
init(deck, eclState, materialLawManager, grid.numCells(), static_cast<const int*>(&grid.globalCell()[0]),
|
|
static_cast<const int*>(&grid.logicalCartesianSize()[0]),
|
|
init_rock);
|
|
}
|
|
#endif
|
|
|
|
/// Constructor for properties on a subgrid
|
|
BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& props,
|
|
std::shared_ptr<MaterialLawManager> materialLawManager,
|
|
const int number_of_cells)
|
|
: rock_(number_of_cells), satprops_(new SaturationPropsFromDeck())
|
|
{
|
|
const int original_size = props.cellPvtRegionIdx_.size();
|
|
if (number_of_cells > original_size) {
|
|
OPM_THROW(std::runtime_error, "The number of cells is larger than the one of the original grid!");
|
|
}
|
|
if (number_of_cells < 0) {
|
|
OPM_THROW(std::runtime_error, "The number of cells is has to be larger than 0.");
|
|
}
|
|
|
|
materialLawManager_ = materialLawManager;
|
|
|
|
// Copy properties that do not depend on the postion within the grid.
|
|
phase_usage_ = props.phase_usage_;
|
|
vap1_ = props.vap1_;
|
|
vap2_ = props.vap2_;
|
|
vap_satmax_guard_ = props.vap_satmax_guard_;
|
|
// For data that is dependant on the subgrid we simply allocate space
|
|
// and initialize with obviously bogus numbers.
|
|
cellPvtRegionIdx_.resize(number_of_cells, std::numeric_limits<int>::min());
|
|
satprops_->init(phase_usage_, materialLawManager_);
|
|
}
|
|
|
|
/// Initializes the properties.
|
|
void BlackoilPropsAdFromDeck::init(const Opm::Deck& deck,
|
|
const Opm::EclipseState& eclState,
|
|
std::shared_ptr<MaterialLawManager> materialLawManager,
|
|
int number_of_cells,
|
|
const int* global_cell,
|
|
const int* cart_dims,
|
|
const bool init_rock)
|
|
{
|
|
materialLawManager_ = materialLawManager;
|
|
|
|
// retrieve the cell specific PVT table index from the deck
|
|
// and using the grid...
|
|
extractPvtTableIndex(cellPvtRegionIdx_, eclState, number_of_cells, global_cell);
|
|
|
|
if (init_rock){
|
|
rock_.init(eclState, number_of_cells, global_cell, cart_dims);
|
|
}
|
|
|
|
phase_usage_ = phaseUsageFromDeck(deck);
|
|
|
|
if (!FluidSystem::isInitialized()) {
|
|
// make sure that we don't initialize the fluid system twice
|
|
FluidSystem::initFromDeck(deck, eclState);
|
|
}
|
|
|
|
// Oil vaporization controls (kw VAPPARS)
|
|
vap1_ = vap2_ = 0.0;
|
|
if (deck.hasKeyword("VAPPARS") && deck.hasKeyword("VAPOIL") && deck.hasKeyword("DISGAS")) {
|
|
vap1_ = deck.getKeyword("VAPPARS").getRecord(0).getItem(0).get< double >(0);
|
|
vap2_ = deck.getKeyword("VAPPARS").getRecord(0).getItem(1).get< double >(0);
|
|
} else if (deck.hasKeyword("VAPPARS")) {
|
|
OPM_THROW(std::runtime_error, "Input has VAPPARS, but missing VAPOIL and/or DISGAS\n");
|
|
}
|
|
satOilMax_.resize(number_of_cells, 0.0);
|
|
|
|
SaturationPropsFromDeck* ptr
|
|
= new SaturationPropsFromDeck();
|
|
satprops_.reset(ptr);
|
|
ptr->init(deck, materialLawManager_);
|
|
|
|
if (phase_usage_.num_phases != satprops_->numPhases()) {
|
|
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck() - "
|
|
"Inconsistent number of phases in pvt data (" << phase_usage_.num_phases
|
|
<< ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
|
|
}
|
|
vap_satmax_guard_ = 0.01;
|
|
}
|
|
|
|
////////////////////////////
|
|
// Rock interface //
|
|
////////////////////////////
|
|
|
|
/// \return D, the number of spatial dimensions.
|
|
int BlackoilPropsAdFromDeck::numDimensions() const
|
|
{
|
|
return rock_.numDimensions();
|
|
}
|
|
|
|
/// \return N, the number of cells.
|
|
int BlackoilPropsAdFromDeck::numCells() const
|
|
{
|
|
return rock_.numCells();
|
|
}
|
|
|
|
/// \return Array of N porosity values.
|
|
const double* BlackoilPropsAdFromDeck::porosity() const
|
|
{
|
|
return rock_.porosity();
|
|
}
|
|
|
|
/// \return Array of ND^2 permeability values.
|
|
/// The D^2 permeability values for a cell are organized as a matrix,
|
|
/// which is symmetric (so ordering does not matter).
|
|
const double* BlackoilPropsAdFromDeck::permeability() const
|
|
{
|
|
return rock_.permeability();
|
|
}
|
|
|
|
|
|
////////////////////////////
|
|
// Fluid interface //
|
|
////////////////////////////
|
|
|
|
/// \return Number of active phases (also the number of components).
|
|
int BlackoilPropsAdFromDeck::numPhases() const
|
|
{
|
|
return phase_usage_.num_phases;
|
|
}
|
|
|
|
/// \return Object describing the active phases.
|
|
PhaseUsage BlackoilPropsAdFromDeck::phaseUsage() const
|
|
{
|
|
return phase_usage_;
|
|
}
|
|
|
|
// ------ Density ------
|
|
|
|
/// Densities of stock components at surface conditions.
|
|
/// \param[in] phaseIdx
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n density values for phase given by phaseIdx.
|
|
V BlackoilPropsAdFromDeck::surfaceDensity(const int phaseIdx, const Cells& cells) const
|
|
{
|
|
assert( !(phaseIdx > numPhases()));
|
|
const int n = cells.size();
|
|
V rhos = V::Zero(n);
|
|
for (int cellIdx = 0; cellIdx < n; ++cellIdx) {
|
|
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
|
|
rhos[cellIdx] = FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
|
|
}
|
|
return rhos;
|
|
}
|
|
|
|
|
|
// ------ Viscosity ------
|
|
|
|
|
|
/// Water viscosity.
|
|
/// \param[in] pw Array of n water pressure values.
|
|
/// \param[in] T Array of n temperature values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n viscosity values.
|
|
ADB BlackoilPropsAdFromDeck::muWat(const ADB& pw,
|
|
const ADB& T,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Water]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(pw.size() == n);
|
|
|
|
V mu(n);
|
|
V dmudp(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/1> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 0.0;
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(pw.value()[i]);
|
|
TEval.setValue(T.value()[i]);
|
|
|
|
const Eval& muEval = FluidSystem::waterPvt().viscosity(pvtRegionIdx, TEval, pEval);
|
|
|
|
mu[i] = muEval.value();
|
|
dmudp[i] = muEval.derivative(0);
|
|
}
|
|
|
|
if (pw.derivative().empty()) {
|
|
return ADB::constant(std::move(mu));
|
|
} else {
|
|
ADB::M dmudp_diag(dmudp.matrix().asDiagonal());
|
|
const int num_blocks = pw.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(dmudp_diag, pw.derivative()[block], jacs[block]);
|
|
}
|
|
return ADB::function(std::move(mu), std::move(jacs));
|
|
}
|
|
}
|
|
|
|
/// Oil viscosity.
|
|
/// \param[in] po Array of n oil pressure values.
|
|
/// \param[in] T Array of n temperature values.
|
|
/// \param[in] rs Array of n gas solution factor values.
|
|
/// \param[in] cond Array of n taxonomies classifying fluid condition.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n viscosity values.
|
|
ADB BlackoilPropsAdFromDeck::muOil(const ADB& po,
|
|
const ADB& T,
|
|
const ADB& rs,
|
|
const std::vector<PhasePresence>& cond,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Oil]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(po.size() == n);
|
|
V mu(n);
|
|
V dmudp(n);
|
|
V dmudr(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/2> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 0.0;
|
|
Eval RsEval = 0.0;
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
RsEval.setDerivative(1, 1.0);
|
|
|
|
Eval muEval;
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(po.value()[i]);
|
|
TEval.setValue(T.value()[i]);
|
|
|
|
if (cond[i].hasFreeGas()) {
|
|
muEval = FluidSystem::oilPvt().saturatedViscosity(pvtRegionIdx, TEval, pEval);
|
|
}
|
|
else {
|
|
if (phase_usage_.phase_used[Gas]) {
|
|
RsEval.setValue(rs.value()[i]);
|
|
}
|
|
muEval = FluidSystem::oilPvt().viscosity(pvtRegionIdx, TEval, pEval, RsEval);
|
|
}
|
|
|
|
mu[i] = muEval.value();
|
|
dmudp[i] = muEval.derivative(0);
|
|
dmudr[i] = muEval.derivative(1);
|
|
}
|
|
|
|
ADB::M dmudp_diag(dmudp.matrix().asDiagonal());
|
|
ADB::M dmudr_diag(dmudr.matrix().asDiagonal());
|
|
const int num_blocks = po.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(dmudp_diag, po.derivative()[block], jacs[block]);
|
|
if (phase_usage_.phase_used[Gas]) {
|
|
ADB::M temp;
|
|
fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
|
|
jacs[block] += temp;
|
|
}
|
|
}
|
|
return ADB::function(std::move(mu), std::move(jacs));
|
|
}
|
|
|
|
/// Gas viscosity.
|
|
/// \param[in] pg Array of n gas pressure values.
|
|
/// \param[in] T Array of n temperature values.
|
|
/// \param[in] rv Array of n vapor oil/gas ratio
|
|
/// \param[in] cond Array of n taxonomies classifying fluid condition.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n viscosity values.
|
|
ADB BlackoilPropsAdFromDeck::muGas(const ADB& pg,
|
|
const ADB& T,
|
|
const ADB& rv,
|
|
const std::vector<PhasePresence>& cond,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Gas]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(pg.value().size() == n);
|
|
V mu(n);
|
|
V dmudp(n);
|
|
V dmudr(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/2> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 0.0;
|
|
Eval RvEval = 0.0;
|
|
Eval muEval;
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
RvEval.setDerivative(1, 1.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(pg.value()[i]);
|
|
TEval.setValue(T.value()[i]);
|
|
|
|
if (cond[i].hasFreeOil()) {
|
|
muEval = FluidSystem::gasPvt().saturatedViscosity(pvtRegionIdx, TEval, pEval);
|
|
}
|
|
else {
|
|
RvEval.setValue(rv.value()[i]);
|
|
muEval = FluidSystem::gasPvt().viscosity(pvtRegionIdx, TEval, pEval, RvEval);
|
|
}
|
|
|
|
mu[i] = muEval.value();
|
|
dmudp[i] = muEval.derivative(0);
|
|
dmudr[i] = muEval.derivative(1);
|
|
}
|
|
|
|
ADB::M dmudp_diag(dmudp.matrix().asDiagonal());
|
|
ADB::M dmudr_diag(dmudr.matrix().asDiagonal());
|
|
const int num_blocks = pg.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(dmudp_diag, pg.derivative()[block], jacs[block]);
|
|
ADB::M temp;
|
|
fastSparseProduct(dmudr_diag, rv.derivative()[block], temp);
|
|
jacs[block] += temp;
|
|
}
|
|
return ADB::function(std::move(mu), std::move(jacs));
|
|
}
|
|
|
|
|
|
// ------ Formation volume factor (b) ------
|
|
|
|
|
|
/// Water formation volume factor.
|
|
/// \param[in] pw Array of n water pressure values.
|
|
/// \param[in] T Array of n temperature values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n formation volume factor values.
|
|
ADB BlackoilPropsAdFromDeck::bWat(const ADB& pw,
|
|
const ADB& T,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Water]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call bWat(): water phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(pw.size() == n);
|
|
|
|
V b(n);
|
|
V dbdp(n);
|
|
V dbdr(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/1> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 0.0;
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(pw.value()[i]);
|
|
TEval.setValue(T.value()[i]);
|
|
|
|
const Eval& bEval = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
|
|
|
|
b[i] = bEval.value();
|
|
dbdp[i] = bEval.derivative(0);
|
|
}
|
|
|
|
ADB::M dbdp_diag(dbdp.matrix().asDiagonal());
|
|
const int num_blocks = pw.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(dbdp_diag, pw.derivative()[block], jacs[block]);
|
|
}
|
|
return ADB::function(std::move(b), std::move(jacs));
|
|
}
|
|
|
|
/// Oil formation volume factor.
|
|
/// \param[in] po Array of n oil pressure values.
|
|
/// \param[in] T Array of n temperature values.
|
|
/// \param[in] rs Array of n gas solution factor values.
|
|
/// \param[in] cond Array of n taxonomies classifying fluid condition.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n formation volume factor values.
|
|
ADB BlackoilPropsAdFromDeck::bOil(const ADB& po,
|
|
const ADB& T,
|
|
const ADB& rs,
|
|
const std::vector<PhasePresence>& cond,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Oil]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call bOil(): oil phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(po.size() == n);
|
|
|
|
V b(n);
|
|
V dbdp(n);
|
|
V dbdr(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/2> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 0.0;
|
|
Eval RsEval = 0.0;
|
|
Eval bEval;
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
RsEval.setDerivative(1, 1.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(po.value()[i]);
|
|
TEval.setValue(T.value()[i]);
|
|
|
|
//RS/RV only makes sense when gas phase is active
|
|
if (cond[i].hasFreeGas()) {
|
|
bEval = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
|
|
}
|
|
else {
|
|
if (rs.size() == 0) {
|
|
RsEval.setValue(0.0);
|
|
}
|
|
else {
|
|
RsEval.setValue(rs.value()[i]);
|
|
}
|
|
bEval = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval, RsEval);
|
|
}
|
|
|
|
b[i] = bEval.value();
|
|
dbdp[i] = bEval.derivative(0);
|
|
dbdr[i] = bEval.derivative(1);
|
|
}
|
|
|
|
ADB::M dbdp_diag(dbdp.matrix().asDiagonal());
|
|
ADB::M dbdr_diag(dbdr.matrix().asDiagonal());
|
|
const int num_blocks = po.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(dbdp_diag, po.derivative()[block], jacs[block]);
|
|
if (phase_usage_.phase_used[Gas]) {
|
|
ADB::M temp;
|
|
fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
|
|
jacs[block] += temp;
|
|
}
|
|
}
|
|
return ADB::function(std::move(b), std::move(jacs));
|
|
}
|
|
|
|
/// Gas formation volume factor.
|
|
/// \param[in] pg Array of n gas pressure values.
|
|
/// \param[in] T Array of n temperature values.
|
|
/// \param[in] rv Array of n vapor oil/gas ratio
|
|
/// \param[in] cond Array of n objects, each specifying which phases are present with non-zero saturation in a cell.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n formation volume factor values.
|
|
ADB BlackoilPropsAdFromDeck::bGas(const ADB& pg,
|
|
const ADB& T,
|
|
const ADB& rv,
|
|
const std::vector<PhasePresence>& cond,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Gas]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call bGas(): gas phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(pg.size() == n);
|
|
|
|
V b(n);
|
|
V dbdp(n);
|
|
V dbdr(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/2> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 0.0;
|
|
Eval RvEval = 0.0;
|
|
Eval bEval;
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
RvEval.setDerivative(1, 1.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(pg.value()[i]);
|
|
TEval.setValue(T.value()[i]);
|
|
|
|
if (cond[i].hasFreeOil()) {
|
|
bEval = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval);
|
|
}
|
|
else {
|
|
RvEval.setValue(rv.value()[i]);
|
|
bEval = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, TEval, pEval, RvEval);
|
|
}
|
|
|
|
b[i] = bEval.value();
|
|
dbdp[i] = bEval.derivative(0);
|
|
dbdr[i] = bEval.derivative(1);
|
|
}
|
|
|
|
ADB::M dbdp_diag(dbdp.matrix().asDiagonal());
|
|
ADB::M dbdr_diag(dbdr.matrix().asDiagonal());
|
|
const int num_blocks = pg.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]);
|
|
ADB::M temp;
|
|
fastSparseProduct(dbdr_diag, rv.derivative()[block], temp);
|
|
jacs[block] += temp;
|
|
}
|
|
return ADB::function(std::move(b), std::move(jacs));
|
|
}
|
|
|
|
|
|
|
|
// ------ Rs bubble point curve ------
|
|
|
|
/// Bubble point curve for Rs as function of oil pressure.
|
|
/// \param[in] po Array of n oil pressure values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n bubble point values for Rs.
|
|
ADB BlackoilPropsAdFromDeck::rsSat(const ADB& po,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Oil]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call rsSat(): oil phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(po.size() == n);
|
|
V rbub(n);
|
|
V drbubdp(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/1> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 293.15; // temperature is not supported by this API!
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(po.value()[i]);
|
|
|
|
const Eval& RsEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx, TEval, pEval);
|
|
|
|
rbub[i] = RsEval.value();
|
|
drbubdp[i] = RsEval.derivative(0);
|
|
}
|
|
|
|
ADB::M drbubdp_diag(drbubdp.matrix().asDiagonal());
|
|
const int num_blocks = po.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(drbubdp_diag, po.derivative()[block], jacs[block]);
|
|
}
|
|
return ADB::function(std::move(rbub), std::move(jacs));
|
|
}
|
|
|
|
/// Bubble point curve for Rs as function of oil pressure.
|
|
/// \param[in] po Array of n oil pressure values.
|
|
/// \param[in] so Array of n oil saturation values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n bubble point values for Rs.
|
|
ADB BlackoilPropsAdFromDeck::rsSat(const ADB& po,
|
|
const ADB& so,
|
|
const Cells& cells) const
|
|
{
|
|
ADB rs = rsSat(po, cells);
|
|
applyVap(rs, so, cells, vap2_);
|
|
return rs;
|
|
}
|
|
|
|
// ------ Rv condensation curve ------
|
|
|
|
/// Condensation curve for Rv as function of oil pressure.
|
|
/// \param[in] pg Array of n gas pressure values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n condensation point values for Rv.
|
|
ADB BlackoilPropsAdFromDeck::rvSat(const ADB& pg,
|
|
const Cells& cells) const
|
|
{
|
|
if (!phase_usage_.phase_used[Gas]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call rvSat(): gas phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
assert(pg.size() == n);
|
|
V rv(n);
|
|
V drvdp(n);
|
|
|
|
typedef Opm::DenseAd::Evaluation<double, /*size=*/1> Eval;
|
|
|
|
Eval pEval = 0.0;
|
|
Eval TEval = 293.15; // temperature is not supported by this API!
|
|
|
|
pEval.setDerivative(0, 1.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
pEval.setValue(pg.value()[i]);
|
|
|
|
const Eval& RvEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx, TEval, pEval);
|
|
|
|
rv[i] = RvEval.value();
|
|
drvdp[i] = RvEval.derivative(0);
|
|
}
|
|
|
|
ADB::M drvdp_diag(drvdp.matrix().asDiagonal());
|
|
const int num_blocks = pg.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
fastSparseProduct(drvdp_diag, pg.derivative()[block], jacs[block]);
|
|
}
|
|
return ADB::function(std::move(rv), std::move(jacs));
|
|
}
|
|
|
|
/// Condensation curve for Rv as function of oil pressure.
|
|
/// \param[in] po Array of n oil pressure values.
|
|
/// \param[in] so Array of n oil saturation values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
|
|
/// \return Array of n condensation point values for Rv.
|
|
ADB BlackoilPropsAdFromDeck::rvSat(const ADB& po,
|
|
const ADB& so,
|
|
const Cells& cells) const
|
|
{
|
|
ADB rv = rvSat(po, cells);
|
|
applyVap(rv, so, cells, vap1_);
|
|
return rv;
|
|
}
|
|
|
|
// ------ Relative permeability ------
|
|
|
|
/// Relative permeabilities for all phases.
|
|
/// \param[in] sw Array of n water saturation values.
|
|
/// \param[in] so Array of n oil saturation values.
|
|
/// \param[in] sg Array of n gas saturation values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
|
/// \return An std::vector with 3 elements, each an array of n relperm values,
|
|
/// containing krw, kro, krg. Use PhaseIndex for indexing into the result.
|
|
std::vector<ADB> BlackoilPropsAdFromDeck::relperm(const ADB& sw,
|
|
const ADB& so,
|
|
const ADB& sg,
|
|
const Cells& cells) const
|
|
{
|
|
const int n = cells.size();
|
|
const int np = numPhases();
|
|
Block s_all(n, np);
|
|
if (phase_usage_.phase_used[Water]) {
|
|
assert(sw.value().size() == n);
|
|
s_all.col(phase_usage_.phase_pos[Water]) = sw.value();
|
|
}
|
|
if (phase_usage_.phase_used[Oil]) {
|
|
assert(so.value().size() == n);
|
|
s_all.col(phase_usage_.phase_pos[Oil]) = so.value();
|
|
} else {
|
|
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
|
|
}
|
|
if (phase_usage_.phase_used[Gas]) {
|
|
assert(sg.value().size() == n);
|
|
s_all.col(phase_usage_.phase_pos[Gas]) = sg.value();
|
|
}
|
|
Block kr(n, np);
|
|
Block dkr(n, np*np);
|
|
satprops_->relperm(n, s_all.data(), cells.data(), kr.data(), dkr.data());
|
|
const int num_blocks = so.numBlocks();
|
|
std::vector<ADB> relperms;
|
|
relperms.reserve(3);
|
|
typedef const ADB* ADBPtr;
|
|
ADBPtr s[3] = { &sw, &so, &sg };
|
|
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
|
if (phase_usage_.phase_used[phase1]) {
|
|
const int phase1_pos = phase_usage_.phase_pos[phase1];
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols());
|
|
}
|
|
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
|
if (!phase_usage_.phase_used[phase2]) {
|
|
continue;
|
|
}
|
|
const int phase2_pos = phase_usage_.phase_pos[phase2];
|
|
// Assemble dkr1/ds2.
|
|
const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
|
|
|
ADB::M dkr1_ds2_diag(dkr.col(column).matrix().asDiagonal());
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
ADB::M temp;
|
|
fastSparseProduct(dkr1_ds2_diag, s[phase2]->derivative()[block], temp);
|
|
jacs[block] += temp;
|
|
}
|
|
}
|
|
ADB::V val = kr.col(phase1_pos);
|
|
relperms.emplace_back(ADB::function(std::move(val), std::move(jacs)));
|
|
} else {
|
|
relperms.emplace_back(ADB::null());
|
|
}
|
|
}
|
|
return relperms;
|
|
}
|
|
|
|
std::vector<ADB> BlackoilPropsAdFromDeck::capPress(const ADB& sw,
|
|
const ADB& so,
|
|
const ADB& sg,
|
|
const Cells& cells) const
|
|
{
|
|
const int nCells = cells.size();
|
|
const int nActivePhases = numPhases();
|
|
const int nBlocks = so.numBlocks();
|
|
|
|
Block activeSat(nCells, nActivePhases);
|
|
if (phase_usage_.phase_used[Water]) {
|
|
assert(sw.value().size() == nCells);
|
|
activeSat.col(phase_usage_.phase_pos[Water]) = sw.value();
|
|
}
|
|
if (phase_usage_.phase_used[Oil]) {
|
|
assert(so.value().size() == nCells);
|
|
activeSat.col(phase_usage_.phase_pos[Oil]) = so.value();
|
|
} else {
|
|
OPM_THROW(std::runtime_error, "BlackoilPropsAdFromDeck::relperm() assumes oil phase is active.");
|
|
}
|
|
if (phase_usage_.phase_used[Gas]) {
|
|
assert(sg.value().size() == nCells);
|
|
activeSat.col(phase_usage_.phase_pos[Gas]) = sg.value();
|
|
}
|
|
|
|
Block pc(nCells, nActivePhases);
|
|
Block dpc(nCells, nActivePhases*nActivePhases);
|
|
satprops_->capPress(nCells, activeSat.data(), cells.data(), pc.data(), dpc.data());
|
|
|
|
std::vector<ADB> adbCapPressures;
|
|
adbCapPressures.reserve(3);
|
|
const ADB* s[3] = { &sw, &so, &sg };
|
|
for (int phase1 = 0; phase1 < 3; ++phase1) {
|
|
if (phase_usage_.phase_used[phase1]) {
|
|
const int phase1_pos = phase_usage_.phase_pos[phase1];
|
|
std::vector<ADB::M> jacs(nBlocks);
|
|
for (int block = 0; block < nBlocks; ++block) {
|
|
jacs[block] = ADB::M(nCells, s[phase1]->derivative()[block].cols());
|
|
}
|
|
for (int phase2 = 0; phase2 < 3; ++phase2) {
|
|
if (!phase_usage_.phase_used[phase2])
|
|
continue;
|
|
const int phase2_pos = phase_usage_.phase_pos[phase2];
|
|
// Assemble dpc1/ds2.
|
|
const int column = phase1_pos + nActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
|
|
ADB::M dpc1_ds2_diag(dpc.col(column).matrix().asDiagonal());
|
|
for (int block = 0; block < nBlocks; ++block) {
|
|
ADB::M temp;
|
|
fastSparseProduct(dpc1_ds2_diag, s[phase2]->derivative()[block], temp);
|
|
jacs[block] += temp;
|
|
}
|
|
}
|
|
ADB::V val = pc.col(phase1_pos);
|
|
adbCapPressures.emplace_back(ADB::function(std::move(val), std::move(jacs)));
|
|
} else {
|
|
adbCapPressures.emplace_back(ADB::null());
|
|
}
|
|
}
|
|
return adbCapPressures;
|
|
}
|
|
|
|
/// Saturation update for hysteresis behavior.
|
|
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
|
void BlackoilPropsAdFromDeck::updateSatHyst(const std::vector<double>& saturation,
|
|
const std::vector<int>& cells)
|
|
{
|
|
const int n = cells.size();
|
|
satprops_->updateSatHyst(n, cells.data(), saturation.data());
|
|
}
|
|
|
|
/// Set gas-oil hysteresis parameters
|
|
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
|
|
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
|
|
void BlackoilPropsAdFromDeck::setGasOilHystParams(const std::vector<double>& pcswmdc,
|
|
const std::vector<double>& krnswdc,
|
|
const std::vector<int>& cells)
|
|
{
|
|
const size_t n = cells.size();
|
|
assert(pcswmdc.size() == n);
|
|
assert(krnswdc.size() == n);
|
|
satprops_->setGasOilHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
|
|
}
|
|
|
|
/// Get gas-oil hysteresis parameters
|
|
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
|
|
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
|
|
void BlackoilPropsAdFromDeck::getGasOilHystParams(std::vector<double>& pcswmdc,
|
|
std::vector<double>& krnswdc,
|
|
const std::vector<int>& cells) const
|
|
{
|
|
const size_t n = cells.size();
|
|
pcswmdc.resize(n);
|
|
krnswdc.resize(n);
|
|
satprops_->getGasOilHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
|
|
}
|
|
|
|
/// Set oil-water hysteresis parameters
|
|
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
|
|
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
|
|
void BlackoilPropsAdFromDeck::setOilWaterHystParams(const std::vector<double>& pcswmdc,
|
|
const std::vector<double>& krnswdc,
|
|
const std::vector<int>& cells)
|
|
{
|
|
const size_t n = cells.size();
|
|
assert(pcswmdc.size() == n);
|
|
assert(krnswdc.size() == n);
|
|
satprops_->setOilWaterHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
|
|
}
|
|
|
|
/// Get oil-water hysteresis parameters
|
|
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
|
|
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
|
|
void BlackoilPropsAdFromDeck::getOilWaterHystParams(std::vector<double>& pcswmdc,
|
|
std::vector<double>& krnswdc,
|
|
const std::vector<int>& cells) const
|
|
{
|
|
const size_t n = cells.size();
|
|
pcswmdc.resize(n);
|
|
krnswdc.resize(n);
|
|
satprops_->getOilWaterHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
|
|
}
|
|
|
|
/// Update for max oil saturation.
|
|
void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
|
|
{
|
|
const int n = satOilMax_.size();
|
|
const int np = phase_usage_.num_phases;
|
|
const int posOil = phase_usage_.phase_pos[Oil];
|
|
const double* s = saturation.data();
|
|
for (int i=0; i<n; ++i) {
|
|
if (satOilMax_[i] < s[np*i+posOil]) {
|
|
satOilMax_[i] = s[np*i+posOil];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Get max oil saturation vector
|
|
const std::vector<double>& BlackoilPropsAdFromDeck::satOilMax() const
|
|
{
|
|
return satOilMax_;
|
|
}
|
|
|
|
// Set max oil saturation vector
|
|
void BlackoilPropsAdFromDeck::setSatOilMax(const std::vector<double>& max_sat) {
|
|
assert(satOilMax_.size() == max_sat.size());
|
|
satOilMax_ = max_sat;
|
|
}
|
|
|
|
/// Bubble point pressures
|
|
std::vector<double> BlackoilPropsAdFromDeck::bubblePointPressure(const Cells& cells,
|
|
const V& T,
|
|
const V& rs) const
|
|
{
|
|
if (!phase_usage_.phase_used[Gas]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call bubblePointPressure(): gas phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
std::vector<double> Pb(n, 0.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
try {
|
|
Pb[i] = FluidSystem::oilPvt().saturationPressure(pvtRegionIdx, T[i], rs[i]);
|
|
}
|
|
catch (const NumericalIssue&) {
|
|
// Ignore
|
|
}
|
|
}
|
|
|
|
return Pb;
|
|
}
|
|
|
|
/// Dew point pressures
|
|
std::vector<double> BlackoilPropsAdFromDeck::dewPointPressure(const Cells& cells,
|
|
const V& T,
|
|
const V& rv) const
|
|
{
|
|
if (!phase_usage_.phase_used[Gas]) {
|
|
OPM_THROW(std::runtime_error, "Cannot call dewPointPressure(): gas phase not active.");
|
|
}
|
|
const int n = cells.size();
|
|
std::vector<double> Pd(n, 0.0);
|
|
|
|
for (int i = 0; i < n; ++i) {
|
|
unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
|
|
try {
|
|
Pd[i] = FluidSystem::gasPvt().saturationPressure(pvtRegionIdx, T[i], rv[i]);
|
|
}
|
|
catch (const NumericalIssue&) {
|
|
// Ignore
|
|
}
|
|
}
|
|
|
|
return Pd;
|
|
}
|
|
|
|
/// Set capillary pressure scaling according to pressure diff. and initial water saturation.
|
|
/// \param[in] saturation Array of n*numPhases saturation values.
|
|
/// \param[in] pc Array of n*numPhases capillary pressure values.
|
|
void BlackoilPropsAdFromDeck::setSwatInitScaling(const std::vector<double>& saturation,
|
|
const std::vector<double>& pc)
|
|
{
|
|
const int nc = rock_.numCells();
|
|
const int numActivePhases = numPhases();
|
|
for (int i = 0; i < nc; ++i) {
|
|
double pcow = pc[numActivePhases*i + phase_usage_.phase_pos[Water]];
|
|
double swat = saturation[numActivePhases*i + phase_usage_.phase_pos[Water]];
|
|
satprops_->swatInitScaling(i, pcow, swat);
|
|
}
|
|
}
|
|
|
|
|
|
/// Apply correction to rs/rv according to kw VAPPARS
|
|
/// \param[in/out] r Array of n rs/rv values.
|
|
/// \param[in] so Array of n oil saturation values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the r and so values.
|
|
/// \param[in] vap Correction parameter.
|
|
void BlackoilPropsAdFromDeck::applyVap(V& r,
|
|
const V& so,
|
|
const std::vector<int>& cells,
|
|
const double vap) const
|
|
{
|
|
if (vap > 0.0) {
|
|
const int n = cells.size();
|
|
V factor = V::Ones(n, 1);
|
|
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
|
for (int i=0; i<n; ++i) {
|
|
if (satOilMax_[cells[i]] > vap_satmax_guard_ && so[i] < satOilMax_[cells[i]]) {
|
|
// guard against too small saturation values.
|
|
const double so_i= std::max(so[i],eps_sqrt);
|
|
factor[i] = std::pow(so_i/satOilMax_[cells[i]], vap);
|
|
}
|
|
}
|
|
r = factor*r;
|
|
}
|
|
}
|
|
|
|
/// Apply correction to rs/rv according to kw VAPPARS
|
|
/// \param[in/out] r Array of n rs/rv values.
|
|
/// \param[in] so Array of n oil saturation values.
|
|
/// \param[in] cells Array of n cell indices to be associated with the r and so values.
|
|
/// \param[in] vap Correction parameter.
|
|
void BlackoilPropsAdFromDeck::applyVap(ADB& r,
|
|
const ADB& so,
|
|
const std::vector<int>& cells,
|
|
const double vap) const
|
|
{
|
|
if (vap > 0.0) {
|
|
const int n = cells.size();
|
|
V factor = V::Ones(n, 1);
|
|
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
|
V dfactor_dso = V::Zero(n, 1);
|
|
for (int i=0; i<n; ++i) {
|
|
if (satOilMax_[cells[i]] > vap_satmax_guard_ && so.value()[i] < satOilMax_[cells[i]]) {
|
|
// guard against too small saturation values.
|
|
const double so_i= std::max(so.value()[i],eps_sqrt);
|
|
factor[i] = std::pow(so_i/satOilMax_[cells[i]], vap);
|
|
dfactor_dso[i] = vap*std::pow(so_i/satOilMax_[cells[i]], vap-1.0)/satOilMax_[cells[i]];
|
|
}
|
|
}
|
|
ADB::M dfactor_dso_diag(dfactor_dso.matrix().asDiagonal());
|
|
const int num_blocks = so.numBlocks();
|
|
std::vector<ADB::M> jacs(num_blocks);
|
|
for (int block = 0; block < num_blocks; ++block) {
|
|
jacs[block] = dfactor_dso_diag * so.derivative()[block];
|
|
}
|
|
r = ADB::function(std::move(factor), std::move(jacs))*r;
|
|
}
|
|
}
|
|
|
|
/// Obtain the scaled critical oil in gas saturation values.
|
|
/// \param[in] cells Array of cell indices.
|
|
/// \return Array of critical oil in gas saturaion values.
|
|
V BlackoilPropsAdFromDeck::scaledCriticalOilinGasSaturations(const Cells& cells) const {
|
|
|
|
|
|
assert(phase_usage_.phase_used[Gas]);
|
|
assert(phase_usage_.phase_used[Oil]);
|
|
|
|
const int n = cells.size();
|
|
V sogcr = V::Zero(n);
|
|
const MaterialLawManager& materialLawManager = satprops_->materialLawManager();
|
|
for (int i = 0; i < n; ++i) {
|
|
const auto& scaledDrainageInfo =
|
|
materialLawManager.oilWaterScaledEpsInfoDrainage(cells[i]);
|
|
sogcr[i] = scaledDrainageInfo.Sogcr;
|
|
}
|
|
return sogcr;
|
|
}
|
|
|
|
/// Obtain the scaled critical gas saturation values.
|
|
/// \param[in] cells Array of cell indices.
|
|
/// \return Array of scaled critical gas saturaion values.
|
|
V BlackoilPropsAdFromDeck::scaledCriticalGasSaturations(const Cells& cells) const {
|
|
|
|
assert(phase_usage_.phase_used[Gas]);
|
|
|
|
const int n = cells.size();
|
|
V sgcr = V::Zero(n);
|
|
const MaterialLawManager& materialLawManager = satprops_->materialLawManager();
|
|
for (int i = 0; i < n; ++i) {
|
|
const auto& scaledDrainageInfo =
|
|
materialLawManager.oilWaterScaledEpsInfoDrainage(cells[i]);
|
|
sgcr[i] = scaledDrainageInfo.Sgcr;
|
|
}
|
|
return sgcr;
|
|
}
|
|
|
|
|
|
} // namespace Opm
|
|
|