diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index 6793fb373..2a9bddff5 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -48,32 +48,36 @@ namespace Opm { * that these only depend on saturation.) */ template > + 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::value, + static_assert(std::is_same::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 diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp index c46a51602..d3be5bf7a 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp @@ -37,12 +37,15 @@ namespace Opm { * Essentially, this class just stores the two parameter objects for * the twophase capillary pressure laws. */ -template +template 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: diff --git a/opm/material/fluidmatrixinteractions/MaterialTraits.hpp b/opm/material/fluidmatrixinteractions/MaterialTraits.hpp index 439391a87..0d5f73d87 100644 --- a/opm/material/fluidmatrixinteractions/MaterialTraits.hpp +++ b/opm/material/fluidmatrixinteractions/MaterialTraits.hpp @@ -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"); }; diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp index ab56a6c33..ef5f4de89 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp @@ -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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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 ¶ms, 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; diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp index 172bad0ee..55d2536a4 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp @@ -43,8 +43,7 @@ class PiecewiseLinearTwoPhaseMaterialParams typedef typename TraitsT::Scalar Scalar; public: - typedef std::pair SamplePoint; - typedef std::vector SamplePoints; + typedef std::vector 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 &SwColumn = swofTable.getSwColumn(); - const std::vector &krwColumn = swofTable.getKrwColumn(); - const std::vector &krowColumn = swofTable.getKrowColumn(); - const std::vector &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 &SgColumn = sgofTable.getSgColumn(); - const std::vector &krgColumn = sgofTable.getKrgColumn(); - const std::vector &krogColumn = sgofTable.getKrogColumn(); - const std::vector &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_;