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_LINEAR_MATERIAL_HPP
00028 #define OPM_LINEAR_MATERIAL_HPP
00029
00030 #include "LinearMaterialParams.hpp"
00031
00032 #include <opm/material/common/MathToolbox.hpp>
00033 #include <opm/common/Valgrind.hpp>
00034 #include <opm/common/Exceptions.hpp>
00035 #include <opm/common/ErrorMacros.hpp>
00036
00037 #include <algorithm>
00038 #include <type_traits>
00039
00040 namespace Opm {
00041
00052 template <class TraitsT, class ParamsT = LinearMaterialParams<TraitsT> >
00053 class LinearMaterial : public TraitsT
00054 {
00055 public:
00056 typedef TraitsT Traits;
00057 typedef ParamsT Params;
00058 typedef typename Traits::Scalar Scalar;
00059
00061 static const int numPhases = Traits::numPhases;
00062
00065 static const bool implementsTwoPhaseApi = (numPhases == 2);
00066
00069 static const bool implementsTwoPhaseSatApi = (numPhases == 2);
00070
00073 static const bool isSaturationDependent = true;
00074
00077 static const bool isPressureDependent = false;
00078
00081 static const bool isTemperatureDependent = false;
00082
00085 static const bool isCompositionDependent = false;
00086
00099 template <class ContainerT, class FluidState>
00100 static void capillaryPressures(ContainerT& values,
00101 const Params& params,
00102 const FluidState& state)
00103 {
00104 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00105
00106 for (unsigned phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
00107 const Evaluation& S =
00108 Opm::decay<Evaluation>(state.saturation(phaseIdx));
00109 Valgrind::CheckDefined(S);
00110
00111 values[phaseIdx] =
00112 S*params.pcMaxSat(phaseIdx) +
00113 (1.0 - S)*params.pcMinSat(phaseIdx);
00114 }
00115 }
00116
00120 template <class ContainerT, class FluidState>
00121 static void saturations(ContainerT& ,
00122 const Params& ,
00123 const FluidState& )
00124 {
00125 OPM_THROW(std::runtime_error, "Not implemented: LinearMaterial::saturations()");
00126 }
00127
00131 template <class ContainerT, class FluidState>
00132 static void relativePermeabilities(ContainerT& values,
00133 const Params& ,
00134 const FluidState& state)
00135 {
00136 typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
00137
00138 for (unsigned phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
00139 const Evaluation& S =
00140 Opm::decay<Evaluation>(state.saturation(phaseIdx));
00141 Valgrind::CheckDefined(S);
00142
00143 values[phaseIdx] = Opm::max(Opm::min(S,1.0),0.0);
00144 }
00145 }
00146
00150 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00151 static Evaluation pcnw(const Params& params, const FluidState& fs)
00152 {
00153 const Evaluation& Sw =
00154 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00155 Valgrind::CheckDefined(Sw);
00156
00157 const Evaluation& wPhasePressure =
00158 Sw*params.pcMaxSat(Traits::wettingPhaseIdx) +
00159 (1.0 - Sw)*params.pcMinSat(Traits::wettingPhaseIdx);
00160
00161 const Evaluation& Sn =
00162 Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
00163 Valgrind::CheckDefined(Sn);
00164
00165 const Evaluation& nPhasePressure =
00166 Sn*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
00167 (1.0 - Sn)*params.pcMinSat(Traits::nonWettingPhaseIdx);
00168
00169 return nPhasePressure - wPhasePressure;
00170 }
00171
00172
00173 template <class Evaluation = Scalar>
00174 static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
00175 twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
00176 {
00177 const Evaluation& wPhasePressure =
00178 Sw*params.pcMaxSat(Traits::wettingPhaseIdx) +
00179 (1.0 - Sw)*params.pcMinSat(Traits::wettingPhaseIdx);
00180
00181 const Evaluation& nPhasePressure =
00182 (1.0 - Sw)*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
00183 Sw*params.pcMinSat(Traits::nonWettingPhaseIdx);
00184
00185 return nPhasePressure - wPhasePressure;
00186 }
00187
00192 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00193 static Evaluation Sw(const Params& , const FluidState& )
00194 { OPM_THROW(std::runtime_error, "Not implemented: Sw()"); }
00195
00196 template <class Evaluation = Scalar>
00197 static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
00198 twoPhaseSatSw(const Params& , const Evaluation& )
00199 { OPM_THROW(std::runtime_error, "Not implemented: twoPhaseSatSw()"); }
00200
00205 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00206 static Evaluation Sn(const Params& , const FluidState& )
00207 { OPM_THROW(std::runtime_error, "Not implemented: Sn()"); }
00208
00209 template <class Evaluation = Scalar>
00210 static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
00211 twoPhaseSatSn(const Params& , const Evaluation& )
00212 { OPM_THROW(std::runtime_error, "Not implemented: twoPhaseSatSn()"); }
00213
00220 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00221 static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
00222 Sg(const Params& , const FluidState& )
00223 { OPM_THROW(std::runtime_error, "Not implemented: Sg()"); }
00224
00228 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00229 static Evaluation krw(const Params& , const FluidState& fs)
00230 {
00231 const Evaluation& Sw =
00232 Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
00233 return Opm::max(0.0, Opm::min(1.0, Sw));
00234 }
00235
00236 template <class Evaluation = Scalar>
00237 static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
00238 twoPhaseSatKrw(const Params& , const Evaluation& Sw)
00239 { return Opm::max(0.0, Opm::min(1.0, Sw)); }
00240
00244 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00245 static Evaluation krn(const Params& , const FluidState& fs)
00246 {
00247 const Evaluation& Sn =
00248 Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
00249 return Opm::max(0.0, Opm::min(1.0, Sn));
00250 }
00251
00252 template <class Evaluation = Scalar>
00253 static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
00254 twoPhaseSatKrn(const Params& , const Evaluation& Sw)
00255 {
00256 return Opm::max(0.0, Opm::min(1.0, Sw));
00257 }
00258
00264 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00265 static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
00266 krg(const Params& , const FluidState& fs)
00267 {
00268 const Evaluation& Sg =
00269 Opm::decay<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
00270 return Opm::max(0.0, Opm::min(1.0, Sg));
00271 }
00272
00278 template <class FluidState, class Evaluation = Scalar>
00279 static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
00280 pcgn(const Params& params, const FluidState& fs)
00281 {
00282 const Evaluation& Sn =
00283 Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
00284 Valgrind::CheckDefined(Sn);
00285
00286 const Evaluation& nPhasePressure =
00287 Sn*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
00288 (1.0 - Sn)*params.pcMinSat(Traits::nonWettingPhaseIdx);
00289
00290 const Evaluation& Sg =
00291 Opm::decay<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
00292 Valgrind::CheckDefined(Sg);
00293
00294 const Evaluation& gPhasePressure =
00295 Sg*params.pcMaxSat(Traits::gasPhaseIdx) +
00296 (1.0 - Sg)*params.pcMinSat(Traits::gasPhaseIdx);
00297
00298 return gPhasePressure - nPhasePressure;
00299 }
00300 };
00301 }
00302
00303 #endif