OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make the EclDefaultMaterial actually work

Browse Source 
in the sense that it is now used by an eWoms simulator
This commit is contained in:
Andreas Lauser
2014-04-28 19:11:53 +02:00
parent 90e1375e0b
commit 5eaedaddf6
5 changed files with 109 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files

36
opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
View File

@@ -48,32 +48,36 @@ namespace Opm {
* that these only depend on saturation.)
*/
template <class TraitsT,
class GasOilMaterial,
class OilWaterMaterial,
class GasOilMaterialLawT,
class OilWaterMaterialLawT,
class ParamsT = EclDefaultMaterialParams<TraitsT,
typename GasOilMaterial::Params,
typename OilWaterMaterial::Params> >
typename GasOilMaterialLawT::Params,
typename OilWaterMaterialLawT::Params> >
class EclDefaultMaterial : public TraitsT
{
public:
typedef GasOilMaterialLawT GasOilMaterialLaw;
typedef OilWaterMaterialLawT OilWaterMaterialLaw;
// some safety checks
static_assert(TraitsT::numPhases == 3,
"The number of phases considered by this capillary pressure "
"law is always three!");
static_assert(GasOilMaterial::numPhases == 2,
static_assert(GasOilMaterialLaw::numPhases == 2,
"The number of phases considered by the gas-oil capillary "
"pressure law must be two!");
static_assert(OilWaterMaterial::numPhases == 2,
static_assert(OilWaterMaterialLaw::numPhases == 2,
"The number of phases considered by the oil-water capillary "
"pressure law must be two!");
static_assert(std::is_same<typename GasOilMaterial::Scalar,
typename OilWaterMaterial::Scalar>::value,
static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
typename OilWaterMaterialLaw::Scalar>::value,
"The two two-phase capillary pressure laws must use the same "
"type of floating point values.");
static_assert(GasOilMaterial::implementsTwoPhaseSatApi,
static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
"The gas-oil material law must implement the two-phase saturation "
"only API to for the default Ecl capillary pressure law!");
static_assert(OilWaterMaterial::implementsTwoPhaseSatApi,
static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
"The oil-water material law must implement the two-phase saturation "
"only API to for the default Ecl capillary pressure law!");
@@ -149,7 +153,7 @@ public:
const FluidState &fs)
{
Scalar Sw = 1 - fs.saturation(gasPhaseIdx);
return GasOilMaterial::twoPhaseSatPcnw(params.gasOilParams(), Sw);
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
}
/*!
@@ -166,7 +170,7 @@ public:
const FluidState &fs)
{
Scalar Sw = fs.saturation(waterPhaseIdx);
return OilWaterMaterial::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
}
/*!
@@ -243,7 +247,7 @@ public:
const FluidState &fluidState)
{
Scalar Sw = 1 - fluidState.saturation(gasPhaseIdx);
return GasOilMaterial::twoPhaseSatKrn(params.oilWaterParams(), Sw);
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
}
/*!
@@ -254,7 +258,7 @@ public:
const FluidState &fluidState)
{
Scalar Sw = fluidState.saturation(waterPhaseIdx);
return OilWaterMaterial::twoPhaseSatKrw(params.oilWaterParams(), Sw);
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
}
/*!
@@ -272,8 +276,8 @@ public:
// probably only relevant if hysteresis is enabled...
Scalar Swco = 0; // todo!
Scalar krog = GasOilMaterial::twoPhaseSatKrw(params.oilWaterParams(), So + Swco);
Scalar krow = OilWaterMaterial::twoPhaseSatKrn(params.oilWaterParams(), 1 - So);
Scalar krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So + Swco);
Scalar krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), 1 - So);
if (Sg + Sw - Swco < 1e-30)
return 1.0; // avoid division by zero

7
opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp
View File

@@ -37,12 +37,15 @@ namespace Opm {
* Essentially, this class just stores the two parameter objects for
* the twophase capillary pressure laws.
*/
template<class Traits, class GasOilParams, class OilWaterParams>
template<class Traits, class GasOilParamsT, class OilWaterParamsT>
class EclDefaultMaterialParams
{
typedef typename Traits::Scalar Scalar;
enum { numPhases = 3 };
public:
typedef GasOilParamsT GasOilParams;
typedef OilWaterParamsT OilWaterParams;
/*!
* \brief The default constructor.
*/
@@ -84,7 +87,7 @@ public:
/*!
* \brief The parameter object for the oil-water twophase law.
*/
void oilWaterParams(const OilWaterParams& val)
void setOilWaterParams(const OilWaterParams& val)
{ oilWaterParams_ = val; }
private:

5
opm/material/fluidmatrixinteractions/MaterialTraits.hpp
View File

@@ -68,11 +68,6 @@ public:
static const int nonWettingPhaseIdx = nonWettingasPhaseIdxV;
// some safety checks...
static_assert(0 <= wettingPhaseIdx && wettingPhaseIdx < numPhases,
"wettingPhaseIdx is out of range");
static_assert(0 <= nonWettingPhaseIdx && nonWettingPhaseIdx < numPhases,
"nonWettingPhaseIdx is out of range");
static_assert(wettingPhaseIdx != nonWettingPhaseIdx,
"wettingPhaseIdx and nonWettingPhaseIdx must be different");
};

104
opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp
View File

@@ -130,7 +130,10 @@ public:
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = 0;
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::nonWettingPhaseIdx] = evalDeriv_(params.pcnwSamples(), state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] =
evalDeriv_(params.SwSamples(),
params.pcnwSamples(),
state.saturation(Traits::wettingPhaseIdx));
}
}
@@ -190,12 +193,12 @@ public:
int satPhaseIdx)
{
if (satPhaseIdx == Traits::wettingPhaseIdx) {
values[Traits::wettingPhaseIdx] = evalDeriv_(params.krwSamples(), state.saturation(Traits::wettingPhaseIdx));
values[Traits::wettingPhaseIdx] = evalDeriv_(params.SwSamples(), params.krwSamples(), state.saturation(Traits::wettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = 0;
}
else {
values[Traits::wettingPhaseIdx] = 0;
values[Traits::nonWettingPhaseIdx] = - evalDeriv_(params.krwSamples(), 1 - state.saturation(Traits::nonWettingPhaseIdx));
values[Traits::nonWettingPhaseIdx] = - evalDeriv_(params.SwSamples(), params.krwSamples(), 1 - state.saturation(Traits::nonWettingPhaseIdx));
}
}
@@ -259,7 +262,7 @@ public:
* \brief The saturation-capillary pressure curve
*/
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
{ return evalDeriv_(params.pcnwSamples(), Sw); }
{ return eval_(params.SwSamples(), params.pcnwSamples(), Sw); }
/*!
* \brief The saturation-capillary pressure curve
@@ -293,7 +296,9 @@ public:
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
assert(0 < Sw && Sw < 1);
return evalDeriv_(params.pcnwSamples(), Sw);
return evalDeriv_(params.SwSamples(),
params.pcnwSamples(),
Sw);
}
/*!
@@ -306,12 +311,12 @@ public:
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{
if (Sw < params.krwSamples().front().first)
return params.krwSamples().front().second;
else if (Sw > params.krwSamples().back().first)
return params.krwSamples().back().second;
if (Sw < params.SwSamples().front())
return params.krwSamples().front();
else if (Sw > params.SwSamples().back())
return params.krwSamples().back();
return eval_(params.krwSamples(), Sw);
return eval_(params.SwSamples(), params.krwSamples(), Sw);
}
/*!
@@ -325,12 +330,14 @@ public:
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
{
if (Sw < params.krwSamples().front().first)
if (Sw < params.SwSamples().front())
return 0;
else if (Sw > params.krwSamples().back().first)
else if (Sw > params.SwSamples().back())
return 0;
return evalDeriv_(params.krwSamples(), Sw);
return evalDeriv_(params.SwSamples(),
params.krwSamples(),
Sw);
}
/*!
@@ -343,12 +350,14 @@ public:
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{
if (Sw < params.krnSamples().front().first)
return params.krnSamples().front().second;
else if (Sw > params.krnSamples().back().first)
return params.krnSamples().back().second;
if (Sw < params.SwSamples().front())
return params.krnSamples().front();
else if (Sw > params.SwSamples().back())
return params.krnSamples().back();
return eval_(params.krnSamples(), Sw);
return eval_(params.SwSamples(),
params.krnSamples(),
Sw);
}
/*!
@@ -362,52 +371,61 @@ public:
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{
if (Sw < params.krwSamples().front().first)
if (Sw < params.SwSamples().front())
return 0;
else if (Sw > params.krwSamples().back().first)
else if (Sw > params.SwSamples().back())
return 0;
return evalDeriv_(params.krnSamples(), Sw);
return evalDeriv_(params.SwSamples(),
params.krnSamples(),
Sw);
}
private:
static Scalar eval_(const SamplePoints &samples, Scalar x)
static Scalar eval_(const SamplePoints &xSamples,
const SamplePoints &ySamples,
Scalar x)
{
int segIdx = findSegmentIndex_(samples, x);
if (x >= samples.front().first && x <= samples.back().first) {
assert(samples[segIdx].first <= x);
assert(samples[segIdx + 1].first >= x);
}
Scalar alpha =
(x - samples[segIdx].first)
/ (samples[segIdx + 1].first - samples[segIdx].first);
return
samples[segIdx].second +
(samples[segIdx + 1].second - samples[segIdx].second)*alpha;
int segIdx = findSegmentIndex_(xSamples, x);
Scalar x0 = xSamples[segIdx];
Scalar x1 = xSamples[segIdx + 1];
Scalar y0 = ySamples[segIdx];
Scalar y1 = ySamples[segIdx + 1];
return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
}
static Scalar evalDeriv_(const SamplePoints &samples, Scalar x)
static Scalar evalDeriv_(const SamplePoints &xSamples,
const SamplePoints &ySamples,
Scalar x)
{
int segIdx = findSegmentIndex_(samples, x);
return
(samples[segIdx + 1].second - samples[segIdx].second)
/(samples[segIdx + 1].first - samples[segIdx].first);
int segIdx = findSegmentIndex_(xSamples, x);
Scalar x0 = xSamples[segIdx];
Scalar x1 = xSamples[segIdx + 1];
Scalar y0 = ySamples[segIdx];
Scalar y1 = ySamples[segIdx + 1];
return (y1 - y0)/(x1 - x0);
}
static int findSegmentIndex_(const SamplePoints &samples, Scalar x)
static int findSegmentIndex_(const SamplePoints &xSamples, Scalar x)
{
int n = samples.size() - 1;
int n = xSamples.size() - 1;
assert(n >= 1); // we need at least two sampling points!
if (samples[n].first < x)
if (xSamples[n] < x)
return n - 1;
else if (samples[0].first > x)
else if (xSamples[0] > x)
return 0;
// bisection
int lowIdx = 0, highIdx = n;
while (lowIdx + 1 < highIdx) {
int curIdx = (lowIdx + highIdx)/2;
if (samples[curIdx].first < x)
if (xSamples[curIdx] < x)
lowIdx = curIdx;
else
highIdx = curIdx;

64
opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp
View File

@@ -43,8 +43,7 @@ class PiecewiseLinearTwoPhaseMaterialParams
typedef typename TraitsT::Scalar Scalar;
public:
typedef std::pair<Scalar, Scalar> SamplePoint;
typedef std::vector<SamplePoint> SamplePoints;
typedef std::vector<Scalar> SamplePoints;
typedef TraitsT Traits;
@@ -64,53 +63,35 @@ public:
#ifndef NDEBUG
finalized_ = true;
#endif
}
/*!
* \brief Read the relevant curves from the table specified by the
* SWOF keyword of an ECLIPSE input file
*/
void readFromSwof(const Opm::SwofTable &swofTable)
{
const std::vector<double> &SwColumn = swofTable.getSwColumn();
const std::vector<double> &krwColumn = swofTable.getKrwColumn();
const std::vector<double> &krowColumn = swofTable.getKrowColumn();
const std::vector<double> &pcowColumn = swofTable.getPcowColumn();
int numRows = swofTable.numRows();
for (int rowIdx = 0; rowIdx < numRows; ++rowIdx) {
pcwnSamples_[rowIdx].first = SwColumn[rowIdx];
pcwnSamples_[rowIdx].second = - pcowColumn[rowIdx];
// revert the order of the sampling points if they were given
// in reverse direction.
if (SwSamples_.front() > SwSamples_.back()) {
for (size_t origSampleIdx = 0;
origSampleIdx < SwSamples_.size() / 2;
++ origSampleIdx)
{
size_t newSampleIdx = SwSamples_.size() - origSampleIdx;
krwSamples_[rowIdx].first = SwColumn[rowIdx];
krwSamples_[rowIdx].second = krwColumn[rowIdx];
krnSamples_[rowIdx].first = SwColumn[rowIdx];
krnSamples_[rowIdx].second = krowColumn[rowIdx];
std::swap(SwSamples_[origSampleIdx], SwSamples_[newSampleIdx]);
std::swap(pcwnSamples_[origSampleIdx], pcwnSamples_[newSampleIdx]);
std::swap(krwSamples_[origSampleIdx], krwSamples_[newSampleIdx]);
std::swap(krnSamples_[origSampleIdx], krnSamples_[newSampleIdx]);
}
}
}
/*!
* \brief Read the relevant curves from the table specified by the
* SGOF keyword of an ECLIPSE input file
* \brief Return the wetting-phase saturation values of all sampling points.
*/
void readFromSgof(const Opm::SgofTable &sgofTable)
{
const std::vector<double> &SgColumn = sgofTable.getSgColumn();
const std::vector<double> &krgColumn = sgofTable.getKrgColumn();
const std::vector<double> &krogColumn = sgofTable.getKrogColumn();
const std::vector<double> &pcogColumn = sgofTable.getPcogColumn();
int numRows = sgofTable.numRows();
for (int rowIdx = 0; rowIdx < numRows; ++rowIdx) {
pcwnSamples_[rowIdx].first = 1 - SgColumn[rowIdx];
pcwnSamples_[rowIdx].second = pcogColumn[rowIdx];
const SamplePoints& SwSamples() const
{ assertFinalized_(); return SwSamples_; }
krwSamples_[rowIdx].first = 1 - SgColumn[rowIdx];
krwSamples_[rowIdx].second = krogColumn[rowIdx];
krnSamples_[rowIdx].first = 1 - SgColumn[rowIdx];
krnSamples_[rowIdx].second = krgColumn[rowIdx];
}
}
/*!
* \brief Set the wetting-phase saturation values of all sampling points.
*/
void setSwSamples(const SamplePoints& samples)
{ SwSamples_ = samples; }
/*!
* \brief Return the sampling points for the capillary pressure curve.
@@ -175,6 +156,7 @@ private:
{ }
#endif
SamplePoints SwSamples_;
SamplePoints pcwnSamples_;
SamplePoints krwSamples_;
SamplePoints krnSamples_;