00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_SPLINE_TWO_PHASE_MATERIAL_HPP
00028 #define OPM_SPLINE_TWO_PHASE_MATERIAL_HPP
00029
00030 #include "SplineTwoPhaseMaterialParams.hpp"
00031
00032 #include <opm/common/ErrorMacros.hpp>
00033 #include <opm/common/Exceptions.hpp>
00034
00035 #include <algorithm>
00036 #include <cmath>
00037 #include <cassert>
00038
00039 namespace Opm {
00046 template <class TraitsT, class ParamsT = SplineTwoPhaseMaterialParams<TraitsT> >
00047 class SplineTwoPhaseMaterial : public TraitsT
00048 {
00049 typedef typename ParamsT::SamplePoints SamplePoints;
00050
00051 public:
00053 typedef TraitsT Traits;
00054
00056 typedef ParamsT Params;
00057
00059 typedef typename Traits::Scalar Scalar;
00060
00062 static const int numPhases = Traits::numPhases;
00063 static_assert(numPhases == 2,
00064 "The piecewise linear two-phase capillary pressure law only"
00065 "applies to the case of two fluid phases");
00066
00069 static const bool implementsTwoPhaseApi = true;
00070
00073 static const bool implementsTwoPhaseSatApi = true;
00074
00077 static const bool isSaturationDependent = true;
00078
00081 static const bool isPressureDependent = false;
00082
00085 static const bool isTemperatureDependent = false;
00086
00089 static const bool isCompositionDependent = false;
00090
00094 template <class Container, class FluidState>
00095 static void capillaryPressures(Container& values, const Params& params, const FluidState& fluidState)
00096 {
00097 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00098
00099 values[Traits::wettingPhaseIdx] = 0.0;
00100 values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fluidState);
00101 }
00102
00107 template <class Container, class FluidState>
00108 static void saturations(Container& , const Params& , const FluidState& )
00109 { OPM_THROW(std::logic_error, "Not implemented: saturations()"); }
00110
00114 template <class Container, class FluidState>
00115 static void relativePermeabilities(Container& values, const Params& params, const FluidState& fluidState)
00116 {
00117 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00118
00119 values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
00120 values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
00121 }
00122
00126 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00127 static Evaluation pcnw(const Params& params, const FluidState& fluidState)
00128 {
00129 const Evaluation& Sw =
00130 Opm::decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
00131
00132 return twoPhaseSatPcnw(params, Sw);
00133 }
00134
00138 template <class Evaluation>
00139 static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00140 {
00141
00142 const auto& pcnwSpline = params.pcnwSpline();
00143 if (Sw <= pcnwSpline.xAt(0))
00144 return Evaluation(pcnwSpline.valueAt(0));
00145 if (Sw >= pcnwSpline.xAt(pcnwSpline.numSamples() - 1))
00146 return Evaluation(pcnwSpline.valueAt(pcnwSpline.numSamples() - 1));
00147
00148 return pcnwSpline.eval(Sw);
00149 }
00150
00151 template <class Evaluation>
00152 static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnw)
00153 {
00154 static const Evaluation nil(0.0);
00155
00156
00157 const auto& pcnwSpline = params.pcnwSpline();
00158 if (pcnw >= pcnwSpline.valueAt(0))
00159 return Evaluation(pcnwSpline.xAt(0));
00160 if (pcnw <= pcnwSpline.y(pcnwSpline.numSamples() - 1))
00161 return Evaluation(pcnwSpline.xAt(pcnwSpline.numSamples() - 1));
00162
00163
00164
00165 return pcnwSpline.intersect(nil, nil, nil, pcnw);
00166 }
00167
00171 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00172 static Evaluation Sw(const Params& , const FluidState& )
00173 { OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
00174
00175 template <class Evaluation>
00176 static Evaluation twoPhaseSatSw(const Params& , const Evaluation& )
00177 { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); }
00178
00183 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00184 static Evaluation Sn(const Params& params, const FluidState& fluidState)
00185 { return 1 - Sw<FluidState, Evaluation>(params, fluidState); }
00186
00187 template <class Evaluation>
00188 static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pC)
00189 { return 1 - twoPhaseSatSw(params, pC); }
00190
00195 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00196 static Evaluation krw(const Params& params, const FluidState& fluidState)
00197 {
00198 const Evaluation& Sw =
00199 Opm::decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
00200
00201 return twoPhaseSatKrw(params, Sw);
00202 }
00203
00204 template <class Evaluation>
00205 static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
00206 {
00207 const auto& krwSpline = params.krwSpline();
00208 if (Sw <= krwSpline.xAt(0))
00209 return Evaluation(krwSpline.valueAt(0));
00210 if (Sw >= krwSpline.xAt(krwSpline.numSamples() - 1))
00211 return Evaluation(krwSpline.valueAt(krwSpline.numSamples() - 1));
00212
00213 return krwSpline.eval(Sw);
00214 }
00215
00216 template <class Evaluation>
00217 static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
00218 {
00219 static const Evaluation nil(0.0);
00220
00221 const auto& krwSpline = params.krwSpline();
00222 if (krw <= krwSpline.valueAt(0))
00223 return Evaluation(krwSpline.xAt(0));
00224 if (krw >= krwSpline.valueAt(krwSpline.numSamples() - 1))
00225 return Evaluation(krwSpline.xAt(krwSpline.numSamples() - 1));
00226
00227 return krwSpline.intersect(nil, nil, nil, krw);
00228 }
00229
00234 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00235 static Evaluation krn(const Params& params, const FluidState& fluidState)
00236 {
00237 const Evaluation& Sn =
00238 Opm::decay<Evaluation>(fluidState.saturation(Traits::nonWettingPhaseIdx));
00239
00240 return twoPhaseSatKrn(params, 1.0 - Sn);
00241 }
00242
00243 template <class Evaluation>
00244 static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
00245 {
00246 const auto& krnSpline = params.krnSpline();
00247 if (Sw <= krnSpline.xAt(0))
00248 return Evaluation(krnSpline.valueAt(0));
00249 if (Sw >= krnSpline.xAt(krnSpline.numSamples() - 1))
00250 return Evaluation(krnSpline.valueAt(krnSpline.numSamples() - 1));
00251
00252 return krnSpline.eval(Sw);
00253 }
00254
00255 template <class Evaluation>
00256 static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
00257 {
00258 static const Evaluation nil(0.0);
00259
00260 const auto& krnSpline = params.krnSpline();
00261 if (krn >= krnSpline.valueAt(0))
00262 return Evaluation(krnSpline.xAt(0));
00263 if (krn <= krnSpline.valueAt(krnSpline.numSamples() - 1))
00264 return Evaluation(krnSpline.xAt(krnSpline.numSamples() - 1));
00265
00266 return krnSpline.intersect(nil, nil, nil, krn);
00267 }
00268 };
00269 }
00270
00271 #endif