implement a two-phase fluid-matrix interaction for ECL end-point scaling

This commit is contained in:
Andreas Lauser 2015-07-28 17:24:09 +02:00
parent 9343d3d720
commit bb100ced1c
5 changed files with 1381 additions and 0 deletions

View File

@ -0,0 +1,224 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
Copyright (C) 2015 by Andreas Lauser
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/>.
*/
/*!
* \file
* \copydoc Opm::EclEpsConfig
*/
#ifndef OPM_ECL_EPS_CONFIG_HPP
#define OPM_ECL_EPS_CONFIG_HPP
#if HAVE_OPM_PARSER
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#endif
#include <string>
#include <cassert>
#include <algorithm>
namespace Opm {
/*!
* \brief Specified which fluids are involved in a given twophase material law for
* endpoint scaling.
*/
enum EclTwoPhaseSystemType {
EclGasOilSystem,
EclOilWaterSystem,
};
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Specifies the configuration used by the endpoint scaling code
*
* This means which quantities are scaled and how this is to be done.
*/
class EclEpsConfig
{
public:
EclEpsConfig()
{
// by default, do not scale anything
enableSatScaling_ = false;
enablePcScaling_ = false;
enableKrwScaling_ = false;
enableKrnScaling_ = false;
enableThreePointKrSatScaling_ = false;
}
/*!
* \brief Specify whether saturation scaling is enabled.
*/
void setEnableSatScaling(bool yesno)
{ enableSatScaling_ = yesno; };
/*!
* \brief Returns whether saturation scaling is enabled.
*/
bool enableSatScaling() const
{ return enableSatScaling_; };
/*!
* \brief Specify whether three point saturation scaling is enabled for the relative
* permeabilities.
*/
void setEnableThreePointKrSatScaling(bool yesno)
{ enableThreePointKrSatScaling_ = yesno; };
/*!
* \brief Returns whether three point saturation scaling is enabled for the relative
* permeabilities.
*/
bool enableThreePointKrSatScaling() const
{ return enableThreePointKrSatScaling_; };
/*!
* \brief Specify whether relative permeability scaling is enabled for the wetting phase.
*/
void setEnableKrwScaling(bool yesno)
{ enableKrwScaling_ = yesno; };
/*!
* \brief Returns whether relative permeability scaling is enabled for the wetting phase.
*/
bool enableKrwScaling() const
{ return enableKrwScaling_; };
/*!
* \brief Specify whether relative permeability scaling is enabled for the non-wetting phase.
*/
void setEnableKrnScaling(bool yesno)
{ enableKrnScaling_ = yesno; };
/*!
* \brief Returns whether relative permeability scaling is enabled for the non-wetting phase.
*/
bool enableKrnScaling() const
{ return enableKrnScaling_; };
/*!
* \brief Specify whether capillary pressure scaling is enabled.
*/
void setEnablePcScaling(bool yesno)
{ enablePcScaling_ = yesno; };
/*!
* \brief Returns whether capillary pressure scaling is enabled.
*/
bool enablePcScaling() const
{ return enablePcScaling_; };
#if HAVE_OPM_PARSER
/*!
* \brief Reads all relevant material parameters form a cell of a parsed ECL deck.
*
* This requires that the opm-parser module is available.
*/
void initFromDeck(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
Opm::EclTwoPhaseSystemType twoPhaseSystemType)
{
// find out if endpoint scaling is used in the first place
if (!deck->hasKeyword("ENDSCALE")) {
// it is not used, i.e., just set all enable$Foo attributes to 0 and be done
// with it.
enableSatScaling_ = false;
enableThreePointKrSatScaling_ = false;
enablePcScaling_ = false;
enableKrwScaling_ = false;
enableKrnScaling_ = false;
return;
}
// endpoint scaling is used, i.e., at least saturation scaling needs to be enabled
enableSatScaling_ = true;
// check if three-point scaling is to be used for the saturations of the relative
// permeabilities
if (deck->hasKeyword("SCALECRS")) {
// if the deck features the SCALECRS keyword, it must be set to 'YES'
Opm::DeckKeywordConstPtr scalecrsKeyword = deck->getKeyword("SCALECRS");
std::string scalecrsValue =
scalecrsKeyword->getRecord(0)->getItem("VALUE")->getString(0);
// convert the value of the SCALECRS keyword to upper case, just to be sure
std::transform(scalecrsValue.begin(),
scalecrsValue.end(),
scalecrsValue.begin(),
::toupper);
if (scalecrsValue == "YES" || scalecrsValue == "Y")
enableThreePointKrSatScaling_ = true;
else
enableThreePointKrSatScaling_ = false;
}
else
enableThreePointKrSatScaling_ = false;
// check if we are supposed to scale the Y axis of the capillary pressure
if (twoPhaseSystemType == EclOilWaterSystem)
enablePcScaling_ =
eclState->hasDoubleGridProperty("PCW")
|| eclState->hasDoubleGridProperty("SWATINIT");
else {
assert(twoPhaseSystemType == EclGasOilSystem);
enablePcScaling_ = eclState->hasDoubleGridProperty("PCG");
}
// check if we are supposed to scale the Y axis of the wetting phase relperm
if (twoPhaseSystemType == EclOilWaterSystem)
enableKrwScaling_ = eclState->hasDoubleGridProperty("KRW");
else {
assert(twoPhaseSystemType == EclGasOilSystem);
enableKrwScaling_ = eclState->hasDoubleGridProperty("KRO");
}
// check if we are supposed to scale the Y axis of the non-wetting phase relperm
if (twoPhaseSystemType == EclOilWaterSystem)
enableKrnScaling_ = eclState->hasDoubleGridProperty("KRO");
else {
assert(twoPhaseSystemType == EclGasOilSystem);
enableKrnScaling_ = eclState->hasDoubleGridProperty("KRG");
}
}
#endif
private:
// enable scaling of the input saturations (i.e., rescale the x-Axis)
bool enableSatScaling_;
// use three (instead of two) points to scale the saturations for the relative
// permeabilities.
//
// This means that two piecewise linear functions are used for saturation scaling
// instead of a single linear one
bool enableThreePointKrSatScaling_;
// enable the scaling of the capillary pressure and relative permeability outputs
// (i.e., rescale the y-Axis)
bool enablePcScaling_;
bool enableKrwScaling_;
bool enableKrnScaling_;
};
} // namespace Opm
#endif

View File

@ -0,0 +1,480 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
Copyright (C) 2015 by Andreas Lauser
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/>.
*/
/*!
* \file
* \copydoc Opm::EclEpsTwoPhaseLawPoints
*/
#ifndef OPM_ECL_EPS_SCALING_POINTS_HPP
#define OPM_ECL_EPS_SCALING_POINTS_HPP
#include "EclEpsConfig.hpp"
#if HAVE_OPM_PARSER
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#endif
#include <string>
#include <iostream>
#include <cassert>
#include <algorithm>
namespace Opm {
/*!
* \brief Collects all grid properties which are relevant for end point scaling.
*
* This class is used for both, the drainage and the imbibition variants of the ECL
* keywords.
*/
class EclEpsGridProperties
{
typedef std::vector<int> IntData;
typedef std::vector<double> DoubleData;
public:
#if HAVE_OPM_PARSER
void initFromDeck(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
bool useImbibition)
{
std::string kwPrefix = useImbibition?"I":"";
if (useImbibition)
satnum = &eclState->getIntGridProperty("IMBNUM")->getData();
else
satnum = &eclState->getIntGridProperty("SATNUM")->getData();
retrieveGridPropertyData_(&swl, deck, eclState, kwPrefix+"SWL");
retrieveGridPropertyData_(&sgl, deck, eclState, kwPrefix+"SGL");
retrieveGridPropertyData_(&swcr, deck, eclState, kwPrefix+"SWCR");
retrieveGridPropertyData_(&sgcr, deck, eclState, kwPrefix+"SGCR");
retrieveGridPropertyData_(&sowcr, deck, eclState, kwPrefix+"SOWCR");
retrieveGridPropertyData_(&sogcr, deck, eclState, kwPrefix+"SOGCR");
retrieveGridPropertyData_(&swu, deck, eclState, kwPrefix+"SWU");
retrieveGridPropertyData_(&sgu, deck, eclState, kwPrefix+"SGU");
retrieveGridPropertyData_(&pcw, deck, eclState, kwPrefix+"PCW");
retrieveGridPropertyData_(&pcg, deck, eclState, kwPrefix+"PCG");
retrieveGridPropertyData_(&krw, deck, eclState, kwPrefix+"KRW");
retrieveGridPropertyData_(&kro, deck, eclState, kwPrefix+"KRO");
retrieveGridPropertyData_(&krg, deck, eclState, kwPrefix+"KRG");
}
#endif
const IntData* satnum;
const DoubleData* swl;
const DoubleData* sgl;
const DoubleData* swcr;
const DoubleData* sgcr;
const DoubleData* sowcr;
const DoubleData* sogcr;
const DoubleData* swu;
const DoubleData* sgu;
const DoubleData* pcw;
const DoubleData* pcg;
const DoubleData* krw;
const DoubleData* kro;
const DoubleData* krg;
private:
#if HAVE_OPM_PARSER
// this method makes sure that a grid property is not created if it is not explicitly
// mentioned in the deck. (saves memory.)
void retrieveGridPropertyData_(const DoubleData **data,
Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
const std::string& properyName)
{
(*data) = 0;
if (eclState->hasDoubleGridProperty(properyName))
(*data) = &eclState->getDoubleGridProperty(properyName)->getData();
}
#endif
};
/*!
* \brief This structure represents all values which can be possibly used as scaling
* points in the endpoint scaling code.
*
* Depending on the exact configuration, some of these quantities are not used as actual
* scaling points. It is easier to extract all of them at once, though.
*/
template <class Scalar>
struct EclEpsScalingPointsInfo
{
// connate saturations
Scalar Swl; // oil
Scalar Sgl; // gas
Scalar Sowl; // oil for the oil-water system
Scalar Sogl; // oil for the gas-oil system
// critical water and gas saturations
Scalar Swcr; // oil
Scalar Sgcr; // gas
Scalar Sowcr; // oil for the oil-water system
Scalar Sogcr; // oil for the gas-oil system
// maximum saturations
Scalar Swu; // oil
Scalar Sgu; // gas
Scalar Sowu; // oil for the oil-water system
Scalar Sogu; // oil for the gas-oil system
// maximum capillary pressures
Scalar maxPcow; // maximum capillary pressure of the oil-water system
Scalar maxPcgo; // maximum capillary pressure of the gas-oil system
// maximum relative permabilities
Scalar maxKrw; // maximum relative permability of water
Scalar maxKrow; // maximum relative permability of oil in the oil-water system
Scalar maxKrog; // maximum relative permability of oil in the gas-oil system
Scalar maxKrg; // maximum relative permability of gas
void print() const
{
std::cout << " Swl: " << Swl << "\n"
<< " Sgl: " << Sgl << "\n"
<< " Sowl: " << Sowl << "\n"
<< " Sogl: " << Sogl << "\n"
<< " Swcr: " << Swcr << "\n"
<< " Sgcr: " << Sgcr << "\n"
<< " Sowcr: " << Sowcr << "\n"
<< " Sogcr: " << Sogcr << "\n"
<< " Swu: " << Swu << "\n"
<< " Sgu: " << Sgu << "\n"
<< " Sowu: " << Sowu << "\n"
<< " Sogu: " << Sogu << "\n"
<< " maxPcow: " << maxPcow << "\n"
<< " maxPcgo: " << maxPcgo << "\n"
<< " maxKrw: " << maxKrw << "\n"
<< " maxKrg: " << maxKrg << "\n"
<< " maxKrow: " << maxKrow << "\n"
<< " maxKrog: " << maxKrog << "\n";
}
#if HAVE_OPM_PARSER
/*!
* \brief Extract the values of the unscaled scaling parameters.
*
* I.e., the values which are used for the nested Fluid-Matrix interactions and which
* are produced by them.
*/
void extractUnscaled(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
int satRegionIdx)
{
// TODO: support for the SOF2/SOF3 keyword family
const auto& swofTable = eclState->getSwofTables()[satRegionIdx];
const auto& sgofTable = eclState->getSgofTables()[satRegionIdx];
// connate saturations
Swl = swofTable.getSwColumn().front();
Sowl = 1.0 - swofTable.getSwColumn().back();
Sgl = sgofTable.getSgColumn().front();
Sogl = 1.0 - sgofTable.getSgColumn().back();
// critical water saturation
for (unsigned rowIdx = 0; rowIdx < swofTable.numRows(); ++ rowIdx) {
if (swofTable.getKrwColumn()[rowIdx] > 0) {
assert(rowIdx > 0);
Swcr = swofTable.getSwColumn()[rowIdx - 1];
break;
};
}
// critical oil saturation of oil-water system
for (int rowIdx = swofTable.numRows() - 1; rowIdx >= 0; -- rowIdx) {
if (swofTable.getKrowColumn()[rowIdx] > 0) {
assert(rowIdx < (int) swofTable.numRows() - 1);
Sowcr = 1.0 - swofTable.getSwColumn()[rowIdx + 1];
break;
};
}
// critical gas saturation
for (unsigned rowIdx = 0; rowIdx < sgofTable.numRows(); ++ rowIdx) {
if (sgofTable.getKrgColumn()[rowIdx] > 0) {
assert(rowIdx > 0);
Sgcr = sgofTable.getSgColumn()[rowIdx - 1];
break;
};
}
// critical oil saturation of gas-oil system
for (int rowIdx = sgofTable.numRows() - 1; rowIdx >= 0; -- rowIdx) {
if (sgofTable.getKrogColumn()[rowIdx] > 0) {
assert(rowIdx < (int) sgofTable.numRows() - 1);
Sogcr = 1.0 - sgofTable.getSgColumn()[rowIdx + 1];
break;
};
}
// maximum saturations
Swu = swofTable.getSwColumn().back();
Sowu = 1.0 - swofTable.getSwColumn().front();
Sgu = sgofTable.getSgColumn().back();
Sogu = 1.0 - sgofTable.getSgColumn().front();
// maximum capillary pressures
maxPcow = swofTable.getPcowColumn().front();
maxPcgo = sgofTable.getPcogColumn().back();
// maximum relative permeabilities
maxKrw = swofTable.getKrwColumn().back();
maxKrow = swofTable.getKrowColumn().front();
maxKrg = sgofTable.getKrgColumn().back();
maxKrog = sgofTable.getKrogColumn().front();
}
#endif
/*!
* \brief Extract the values of the scaled scaling parameters.
*
* I.e., the values which are "seen" by the physical model.
*/
void extractScaled(const EclEpsGridProperties& epsProperties, int cartesianCellIdx)
{
// overwrite the unscaled values with the values for the cell if it is
// explicitly specified by the corresponding keyword.
extractGridPropertyValue_(Swl, epsProperties.swl, cartesianCellIdx);
extractGridPropertyValue_(Sgl, epsProperties.sgl, cartesianCellIdx);
extractGridPropertyValue_(Swcr, epsProperties.swcr, cartesianCellIdx);
extractGridPropertyValue_(Sgcr, epsProperties.sgcr, cartesianCellIdx);
extractGridPropertyValue_(Sowcr, epsProperties.sowcr, cartesianCellIdx);
extractGridPropertyValue_(Sogcr, epsProperties.sogcr, cartesianCellIdx);
extractGridPropertyValue_(Swu, epsProperties.swu, cartesianCellIdx);
extractGridPropertyValue_(Sgu, epsProperties.sgu, cartesianCellIdx);
extractGridPropertyValue_(maxPcow, epsProperties.pcw, cartesianCellIdx);
extractGridPropertyValue_(maxPcgo, epsProperties.pcg, cartesianCellIdx);
extractGridPropertyValue_(maxKrw, epsProperties.krw, cartesianCellIdx);
extractGridPropertyValue_(maxKrg, epsProperties.krg, cartesianCellIdx);
// quite likely that's wrong!
extractGridPropertyValue_(maxKrow, epsProperties.kro, cartesianCellIdx);
extractGridPropertyValue_(maxKrog, epsProperties.kro, cartesianCellIdx);
}
private:
void extractGridPropertyValue_(Scalar& targetValue,
const std::vector<double>* propData,
int cartesianCellIdx)
{
if (!propData)
return;
targetValue = (*propData)[cartesianCellIdx];
}
};
/*!
* \ingroup FluidMatrixInteractions
*
* \brief Represents the points on the X and Y axis to be scaled if endpoint scaling is
* used.
*/
template <class Scalar>
class EclEpsScalingPoints
{
public:
/*!
* \brief Assigns the scaling points which actually ought to be used.
*/
void init(const EclEpsScalingPointsInfo<Scalar>& epsInfo,
const EclEpsConfig& config,
EclTwoPhaseSystemType epsSystemType)
{
if (epsSystemType == EclOilWaterSystem) {
// saturation scaling for capillary pressure
saturationPcPoints_[0] = epsInfo.Swl;
saturationPcPoints_[1] = epsInfo.Swu;
// krw saturation scaling endpoints
if (config.enableThreePointKrSatScaling()) {
saturationKrwPoints_[0] = epsInfo.Swcr;
saturationKrwPoints_[1] = 1.0 - epsInfo.Sowcr - epsInfo.Sgl;
saturationKrwPoints_[2] = epsInfo.Swu;
}
else {
saturationKrwPoints_[0] = epsInfo.Swcr;
saturationKrwPoints_[1] = epsInfo.Swu;
}
// krn saturation scaling endpoints (with the non-wetting phase being oil).
// because opm-material specifies non-wetting phase relperms in terms of the
// wetting phase saturations, the code here uses 1 minus the values specified
// by the Eclipse TD and the order of the scaling points is reversed
if (config.enableThreePointKrSatScaling()) {
saturationKrnPoints_[2] = 1.0 - epsInfo.Sowcr;
saturationKrnPoints_[1] = epsInfo.Swcr + epsInfo.Sgl;
saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
}
else {
saturationKrnPoints_[1] = 1 - epsInfo.Sowcr;
saturationKrnPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
}
maxPcnw_ = epsInfo.maxPcow;
maxKrw_ = epsInfo.maxKrw;
maxKrn_ = epsInfo.maxKrow;
}
else {
assert(epsSystemType == EclGasOilSystem);
// saturation scaling for capillary pressure
saturationPcPoints_[0] = 1.0 - epsInfo.Sgu;
saturationPcPoints_[1] = 1.0 - epsInfo.Sgl;
// krw saturation scaling endpoints
if (config.enableThreePointKrSatScaling()) {
saturationKrwPoints_[0] = epsInfo.Sogcr;
saturationKrwPoints_[1] = 1 - epsInfo.Sgcr - epsInfo.Swl;
saturationKrwPoints_[2] = 1 - epsInfo.Swl - epsInfo.Sgl;
}
else {
saturationKrwPoints_[0] = epsInfo.Swl + epsInfo.Sgl;
saturationKrwPoints_[1] = 1.0 - epsInfo.Sogcr;
}
// krn saturation scaling endpoints (with the non-wetting phase being gas).
// because opm-material specifies non-wetting phase relperms in terms of the
// wetting phase saturations, the code here uses 1 minus the values specified
// by the Eclipse TD and the order of the scaling points is reversed
if (config.enableThreePointKrSatScaling()) {
saturationKrnPoints_[2] = 1.0 - epsInfo.Sgcr;
saturationKrnPoints_[1] = epsInfo.Sogcr + epsInfo.Swl;
saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
}
else {
saturationKrnPoints_[1] = 1.0 - epsInfo.Sgcr;
saturationKrnPoints_[0] = 1.0 - epsInfo.Sgu;
}
maxPcnw_ = epsInfo.maxPcgo;
maxKrw_ = epsInfo.maxKrog;
maxKrn_ = epsInfo.maxKrg;
}
}
/*!
* \brief Sets an saturation value for capillary pressure saturation scaling
*/
void setSaturationPcPoint(int pointIdx, Scalar value)
{ saturationPcPoints_[pointIdx] = value; }
/*!
* \brief Returns the points used for capillary pressure saturation scaling
*/
const std::array<Scalar, 2>& saturationPcPoints() const
{ return saturationPcPoints_; }
/*!
* \brief Sets an saturation value for wetting-phase relperm saturation scaling
*/
void setSaturationKrwPoint(int pointIdx, Scalar value)
{ saturationKrwPoints_[pointIdx] = value; }
/*!
* \brief Returns the points used for wetting phase relperm saturation scaling
*/
const std::array<Scalar, 3>& saturationKrwPoints() const
{ return saturationKrwPoints_; }
/*!
* \brief Sets an saturation value for non-wetting phase relperm saturation scaling
*/
void setSaturationKrnPoint(int pointIdx, Scalar value)
{ saturationKrnPoints_[pointIdx] = value; }
/*!
* \brief Returns the points used for non-wetting phase relperm saturation scaling
*/
const std::array<Scalar, 3>& saturationKrnPoints() const
{ return saturationKrnPoints_; }
/*!
* \brief Sets the maximum capillary pressure
*/
void setMaxPcnw(Scalar value)
{ maxPcnw_ = value; }
/*!
* \brief Returns the maximum capillary pressure
*/
Scalar maxPcnw() const
{ return maxPcnw_; }
/*!
* \brief Sets the maximum wetting phase relative permeability
*/
void setMaxKrw(Scalar value)
{ maxKrw_ = value; }
/*!
* \brief Returns the maximum wetting phase relative permeability
*/
Scalar maxKrw() const
{ return maxKrw_; }
/*!
* \brief Sets the maximum wetting phase relative permeability
*/
void setMaxKrn(Scalar value)
{ maxKrn_ = value; }
/*!
* \brief Returns the maximum wetting phase relative permeability
*/
Scalar maxKrn() const
{ return maxKrn_; }
void print() const
{
std::cout << " saturationKrnPoints_[0]: " << saturationKrnPoints_[0] << "\n"
<< " saturationKrnPoints_[1]: " << saturationKrnPoints_[1] << "\n"
<< " saturationKrnPoints_[2]: " << saturationKrnPoints_[2] << "\n";
}
private:
// The the points used for the "y-axis" scaling of capillary pressure
Scalar maxPcnw_;
// The the points used for the "y-axis" scaling of wetting phase relative permability
Scalar maxKrw_;
// The the points used for the "y-axis" scaling of non-wetting phase relative permability
Scalar maxKrn_;
// The the points used for saturation ("x-axis") scaling of capillary pressure
std::array<Scalar, 2> saturationPcPoints_;
// The the points used for saturation ("x-axis") scaling of wetting phase relative permeability
std::array<Scalar, 3> saturationKrwPoints_;
// The the points used for saturation ("x-axis") scaling of non-wetting phase relative permeability
std::array<Scalar, 3> saturationKrnPoints_;
};
} // namespace Opm
#endif

View File

@ -0,0 +1,516 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
Copyright (C) 2015 by Andreas Lauser
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/>.
*/
/*!
* \file
* \copydoc Opm::EclEpsTwoPhaseLaw
*/
#ifndef OPM_ECL_EPS_TWO_PHASE_LAW_HPP
#define OPM_ECL_EPS_TWO_PHASE_LAW_HPP
#include "EclEpsTwoPhaseLawParams.hpp"
#include <opm/material/fluidstates/SaturationOverlayFluidState.hpp>
#include <opm/material/common/ErrorMacros.hpp>
#include <opm/material/common/Exceptions.hpp>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief This material law takes a material law defined for effective
* saturations and converts it to a material law defined on absolute
* saturations.
*
* This class provides the endpoint scaling functionality used by the ECL reservoir
* simulator. The basic purpose of this class is the same as the one of \a EffToAbsLaw,
* but it is quite a bit more complex.
*/
template <class EffLawT,
class ParamsT = EclEpsTwoPhaseLawParams<EffLawT> >
class EclEpsTwoPhaseLaw : public EffLawT::Traits
{
typedef EffLawT EffLaw;
public:
typedef typename EffLaw::Traits Traits;
typedef ParamsT Params;
typedef typename EffLaw::Scalar Scalar;
enum { wettingPhaseIdx = Traits::wettingPhaseIdx };
enum { nonWettingPhaseIdx = Traits::nonWettingPhaseIdx };
//! The number of fluid phases
static const int numPhases = EffLaw::numPhases;
static_assert(numPhases == 2,
"The endpoint scaling applies to the nested twophase laws, not to "
"the threephase one!");
//! Specify whether this material law implements the two-phase
//! convenience API
static const bool implementsTwoPhaseApi = true;
static_assert(EffLaw::implementsTwoPhaseApi,
"The material laws put into EclEpsTwoPhaseLaw must implement the "
"two-phase material law API!");
//! Specify whether this material law implements the two-phase
//! convenience API which only depends on the phase saturations
static const bool implementsTwoPhaseSatApi = true;
static_assert(EffLaw::implementsTwoPhaseSatApi,
"The material laws put into EclEpsTwoPhaseLaw must implement the "
"two-phase material law saturation API!");
//! Specify whether the quantities defined by this material law
//! are saturation dependent
static const bool isSaturationDependent = true;
//! Specify whether the quantities defined by this material law
//! are dependent on the absolute pressure
static const bool isPressureDependent = false;
//! Specify whether the quantities defined by this material law
//! are temperature dependent
static const bool isTemperatureDependent = false;
//! Specify whether the quantities defined by this material law
//! are dependent on the phase composition
static const bool isCompositionDependent = false;
/*!
* \brief The capillary pressure-saturation curves depending on absolute saturations.
*
* \param values A random access container which stores the
* relative pressure of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the capillary pressure
* ought to be calculated
*/
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The capillaryPressures(fs) method is not yet implemented");
}
/*!
* \brief The relative permeability-saturation curves depending on absolute saturations.
*
* \param values A random access container which stores the
* relative permeability of each fluid phase.
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the relative permeabilities
* ought to be calculated
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The pcnw(fs) method is not yet implemented");
}
/*!
* \brief The capillary pressure-saturation curve.
*
*
* \param params A object that stores the appropriate coefficients
* for the respective law.
*
* \return Capillary pressure [Pa] calculated by specific
* constitutive relation (e.g. Brooks & Corey, van
* Genuchten, linear...)
*/
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The pcnw(fs) method is not yet implemented");
}
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& SwAbs)
{
const Evaluation& SwEff = effectiveSaturationPc(params, SwAbs);
const Evaluation& pcUnscaled = EffLaw::twoPhaseSatPcnw(params.effectiveLawParams(), SwEff);
return scaleCapillaryPressure_(params, pcUnscaled);
}
template <class Evaluation>
static Evaluation twoPhaseSatPcnwInv(const Params &params, const Evaluation& pcnw)
{
Evaluation pcnwUnscaled = scaleCapillaryPressureInv_(params, pcnw);
Evaluation SwUnscaled = EffLaw::twoPhaseSatPcnwInv(params.effectiveLawParams(), pcnwUnscaled);
return effectiveSaturationPcInv(params, SwUnscaled);
}
/*!
* \brief The saturation-capillary pressure curves.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The saturations(fs) method is not yet implemented");
}
/*!
* \brief Calculate wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The Sw(fs) method is not yet implemented");
}
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pc)
{
OPM_THROW(NotAvailable,
"The twoPhaseSatSw(pc) method is not yet implemented");
}
/*!
* \brief Calculate non-wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The Sn(pc) method is not yet implemented");
}
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pc)
{
OPM_THROW(NotAvailable,
"The twoPhaseSatSn(pc) method is not yet implemented");
}
/*!
* \brief The relative permeability for the wetting phase.
*
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Relative permeability of the wetting phase calculated as implied by EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*
*/
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The krw(fs) method is not yet implemented");
}
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
const Evaluation& Swe = effectiveSaturationKrw(params, Sw);
const Evaluation& rawKrw = EffLaw::twoPhaseSatKrw(params.effectiveLawParams(), Swe);
return scaleKrw_(params, rawKrw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrwInv(const Params &params, const Evaluation& krw)
{
Evaluation krwUnscaled = scaleKrwInv_(params, krw);
Evaluation SwUnscaled = EffLaw::twoPhaseSatKrwInv(params.effectiveLawParams(), krwUnscaled);
return effectiveSaturationKrwInv(params, SwUnscaled);
}
/*!
* \brief The relative permeability of the non-wetting phase.
*/
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
OPM_THROW(NotAvailable,
"The krn(fs) method is not yet implemented");
}
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
const Evaluation& Swe = effectiveSaturationKrn(params, Sw);
const Evaluation& rawKrn = EffLaw::twoPhaseSatKrn(params.effectiveLawParams(), Swe);
return scaleKrn_(params, rawKrn);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrnInv(const Params &params, const Evaluation& krn)
{
Evaluation krnUnscaled = scaleKrnInv_(params, krn);
Evaluation SwUnscaled = EffLaw::twoPhaseSatKrnInv(params.effectiveLawParams(), krnUnscaled);
return effectiveSaturationKrnInv(params, SwUnscaled);
}
/*!
* \brief Convert an absolute saturation to an effective one for capillary pressure.
*
* The effective saturation is then feed into the "raw" capillary pressure law.
*/
template <class Evaluation>
static Evaluation effectiveSaturationPc(const Params &params, const Evaluation& Sw)
{
if (!params.config().enableSatScaling())
return Sw;
// the saturations of capillary pressure are always scaled using two-point
// scaling
return scaleSatTwoPoint_(Sw,
params.unscaledPoints().saturationPcPoints(),
params.scaledPoints().saturationPcPoints());
}
template <class Evaluation>
static Evaluation effectiveSaturationPcInv(const Params &params, const Evaluation& Sw)
{
if (!params.config().enableSatScaling())
return Sw;
// the saturations of capillary pressure are always scaled using two-point
// scaling
return scaleSatTwoPointInv_(Sw,
params.unscaledPoints().saturationPcPoints(),
params.scaledPoints().saturationPcPoints());
}
/*!
* \brief Convert an absolute saturation to an effective one for the scaling of the
* relperm of the wetting phase.
*/
template <class Evaluation>
static Evaluation effectiveSaturationKrw(const Params &params, const Evaluation& Sw)
{
if (!params.config().enableSatScaling())
return Sw;
if (params.config().enableThreePointKrSatScaling()) {
return scaleSatThreePoint_(Sw,
params.unscaledPoints().saturationKrwPoints(),
params.scaledPoints().saturationKrwPoints());
}
else { // two-point relperm saturation scaling
return scaleSatTwoPoint_(Sw,
params.unscaledPoints().saturationKrwPoints(),
params.scaledPoints().saturationKrwPoints());
}
}
template <class Evaluation>
static Evaluation effectiveSaturationKrwInv(const Params &params, const Evaluation& Sw)
{
if (!params.config().enableSatScaling())
return Sw;
if (params.config().enableThreePointKrSatScaling()) {
return scaleSatThreePointInv_(Sw,
params.unscaledPoints().saturationKrwPoints(),
params.scaledPoints().saturationKrwPoints());
}
else { // two-point relperm saturation scaling
return scaleSatTwoPointInv_(Sw,
params.unscaledPoints().saturationKrwPoints(),
params.scaledPoints().saturationKrwPoints());
}
}
/*!
* \brief Convert an absolute saturation to an effective one for the scaling of the
* relperm of the non-wetting phase.
*/
template <class Evaluation>
static Evaluation effectiveSaturationKrn(const Params &params, const Evaluation& Sw)
{
if (!params.config().enableSatScaling())
return Sw;
if (params.config().enableThreePointKrSatScaling())
return scaleSatThreePoint_(Sw,
params.unscaledPoints().saturationKrnPoints(),
params.scaledPoints().saturationKrnPoints());
else // two-point relperm saturation scaling
return scaleSatTwoPoint_(Sw,
params.unscaledPoints().saturationKrnPoints(),
params.scaledPoints().saturationKrnPoints());
}
template <class Evaluation>
static Evaluation effectiveSaturationKrnInv(const Params &params, const Evaluation& Sw)
{
if (!params.config().enableSatScaling())
return Sw;
if (params.config().enableThreePointKrSatScaling()) {
return scaleSatThreePointInv_(Sw,
params.unscaledPoints().saturationKrnPoints(),
params.scaledPoints().saturationKrnPoints());
}
else { // two-point relperm saturation scaling
return scaleSatTwoPointInv_(Sw,
params.unscaledPoints().saturationKrnPoints(),
params.scaledPoints().saturationKrnPoints());
}
}
private:
template <class Evaluation, class PointsContainer>
static Evaluation scaleSatTwoPoint_(const Evaluation& S,
const PointsContainer& unscaledSats,
const PointsContainer& scaledSats)
{
return
unscaledSats[0]
+
(S - scaledSats[0])*((unscaledSats[1] - unscaledSats[0])/(scaledSats[1] - scaledSats[0]));
}
template <class Evaluation, class PointsContainer>
static Evaluation scaleSatTwoPointInv_(const Evaluation& S,
const PointsContainer& unscaledSats,
const PointsContainer& scaledSats)
{
return
scaledSats[0]
+
(S - unscaledSats[0])*((scaledSats[1] - scaledSats[0])/(unscaledSats[1] - unscaledSats[0]));
}
template <class Evaluation, class PointsContainer>
static Evaluation scaleSatThreePoint_(const Evaluation& S,
const PointsContainer& unscaledSats,
const PointsContainer& scaledSats)
{
if (unscaledSats[1] >= unscaledSats[2])
return scaleSatTwoPoint_(S, unscaledSats, scaledSats);
if (S < scaledSats[1])
return
unscaledSats[0]
+
(S - scaledSats[0])*((unscaledSats[1] - unscaledSats[0])/(scaledSats[1] - scaledSats[0]));
else
return
unscaledSats[1]
+
(S - scaledSats[1])*((unscaledSats[2] - unscaledSats[1])/(scaledSats[2] - scaledSats[1]));
}
template <class Evaluation, class PointsContainer>
static Evaluation scaleSatThreePointInv_(const Evaluation& S,
const PointsContainer& unscaledSats,
const PointsContainer& scaledSats)
{
if (unscaledSats[1] >= unscaledSats[2])
return scaleSatTwoPointInv_(S, unscaledSats, scaledSats);
if (S < unscaledSats[1])
return
scaledSats[0]
+
(S - unscaledSats[0])*((scaledSats[1] - scaledSats[0])/(unscaledSats[1] - unscaledSats[0]));
else
return
scaledSats[1]
+
(S - unscaledSats[1])*((scaledSats[2] - scaledSats[1])/(unscaledSats[2] - unscaledSats[1]));
}
/*!
* \brief Scale the capillary pressure according to the given parameters
*/
template <class Evaluation>
static Evaluation scaleCapillaryPressure_(const Params &params, const Evaluation& pc)
{
if (!params.config().enablePcScaling())
return pc;
Scalar alpha = params.scaledPoints().maxPcnw()/params.unscaledPoints().maxPcnw();
return pc*alpha;
}
template <class Evaluation>
static Evaluation scaleCapillaryPressureInv_(const Params &params, const Evaluation& pc)
{
if (!params.config().enablePcScaling())
return pc;
Scalar alpha = params.unscaledPoints().maxPcnw()/params.scaledPoints().maxPcnw();
return pc/alpha;
}
/*!
* \brief Scale the wetting phase relative permeability of a phase according to the given parameters
*/
template <class Evaluation>
static Evaluation scaleKrw_(const Params &params, const Evaluation& krw)
{
if (!params.config().enableKrwScaling())
return krw;
// TODO: three point krw y-scaling
Scalar alpha = params.scaledPoints().maxKrw()/params.unscaledPoints().maxKrw();
return krw*alpha;
}
template <class Evaluation>
static Evaluation scaleKrwInv_(const Params &params, const Evaluation& krw)
{
if (!params.config().enableKrwScaling())
return krw;
Scalar alpha = params.unscaledPoints().maxKrw()/params.scaledPoints().maxKrw();
return krw*alpha;
}
/*!
* \brief Scale the non-wetting phase relative permeability of a phase according to the given parameters
*/
template <class Evaluation>
static Evaluation scaleKrn_(const Params &params, const Evaluation& krn)
{
if (!params.config().enableKrnScaling())
return krn;
//TODO: three point krn y-scaling
Scalar alpha = params.scaledPoints().maxKrn()/params.unscaledPoints().maxKrn();
return krn*alpha;
}
template <class Evaluation>
static Evaluation scaleKrnInv_(const Params &params, const Evaluation& krn)
{
if (!params.config().enableKrnScaling())
return krn;
Scalar alpha = params.unscaledPoints().maxKrn()/params.scaledPoints().maxKrn();
return krn*alpha;
}
};
} // namespace Opm
#endif

View File

@ -0,0 +1,152 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
Copyright (C) 2015 by Andreas Lauser
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/>.
*/
/*!
* \file
* \copydoc Opm::EclEpsTwoPhaseLawParams
*/
#ifndef OPM_ECL_EPS_TWO_PHASE_LAW_PARAMS_HPP
#define OPM_ECL_EPS_TWO_PHASE_LAW_PARAMS_HPP
#include "EclEpsConfig.hpp"
#include "EclEpsScalingPoints.hpp"
#if HAVE_OPM_PARSER
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#endif
#include <string>
#include <memory>
#include <cassert>
#include <algorithm>
namespace Opm {
/*!
* \ingroup FluidMatrixInteractions
*
* \brief A default implementation of the parameters for the material law adapter class
* which implements ECL endpoint scaleing .
*/
template <class EffLawT>
class EclEpsTwoPhaseLawParams
{
typedef typename EffLawT::Params EffLawParams;
typedef typename EffLawParams::Traits::Scalar Scalar;
public:
typedef typename EffLawParams::Traits Traits;
typedef Opm::EclEpsScalingPoints<Scalar> ScalingPoints;
EclEpsTwoPhaseLawParams()
{
#ifndef NDEBUG
finalized_ = false;
#endif
}
/*!
* \brief Calculate all dependent quantities once the independent
* quantities of the parameter object have been set.
*/
void finalize()
{
#ifndef NDEBUG
assert(config_);
if (config_->enableSatScaling()) {
assert(unscaledPoints_);
assert(scaledPoints_);
}
assert(effectiveLawParams_);
finalized_ = true;
#endif
}
/*!
* \brief Set the endpoint scaling configuration object.
*/
void setConfig(std::shared_ptr<EclEpsConfig> value)
{ config_ = value; }
/*!
* \brief Returns the endpoint scaling configuration object.
*/
const EclEpsConfig& config() const
{ return *config_; }
/*!
* \brief Set the scaling points which are seen by the nested material law
*/
void setUnscaledPoints(std::shared_ptr<ScalingPoints> value)
{ unscaledPoints_ = value; }
/*!
* \brief Returns the scaling points which are seen by the nested material law
*/
const ScalingPoints& unscaledPoints() const
{ return *unscaledPoints_; }
/*!
* \brief Set the scaling points which are seen by the physical model
*/
void setScaledPoints(std::shared_ptr<ScalingPoints> value)
{ scaledPoints_ = value; }
/*!
* \brief Returns the scaling points which are seen by the physical model
*/
const ScalingPoints& scaledPoints() const
{ return *scaledPoints_; }
/*!
* \brief Sets the parameter object for the effective/nested material law.
*/
void setEffectiveLawParams(std::shared_ptr<EffLawParams> value)
{ effectiveLawParams_ = value; }
/*!
* \brief Returns the parameter object for the effective/nested material law.
*/
const EffLawParams& effectiveLawParams() const
{ return *effectiveLawParams_; }
private:
#ifndef NDEBUG
void assertFinalized_() const
{ assert(finalized_); }
bool finalized_;
#else
void assertFinalized_() const
{ }
#endif
std::shared_ptr<EffLawParams> effectiveLawParams_;
std::shared_ptr<EclEpsConfig> config_;
std::shared_ptr<ScalingPoints> unscaledPoints_;
std::shared_ptr<ScalingPoints> scaledPoints_;
};
} // namespace Opm
#endif

View File

@ -43,6 +43,7 @@
#include <opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp>
#include <opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp>
// include the helper classes to construct traits
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
@ -375,5 +376,13 @@ int main(int argc, char **argv)
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
{
typedef Opm::BrooksCorey<TwoPhaseTraits> RawMaterialLaw;
typedef Opm::EclEpsTwoPhaseLaw<RawMaterialLaw> MaterialLaw;
testGenericApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseApi<MaterialLaw, TwoPhaseFluidState>();
testTwoPhaseSatApi<MaterialLaw, TwoPhaseFluidState>();
}
return 0;
}