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_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
00028 #define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
00029
00030 #include <opm/material/common/MathToolbox.hpp>
00031
00032
00033 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00034
00035 #include <dune/common/fvector.hh>
00036 #include <dune/common/fmatrix.hh>
00037
00038 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00039
00040
00041 #include <opm/common/Exceptions.hpp>
00042 #include <opm/common/ErrorMacros.hpp>
00043 #include <opm/common/Valgrind.hpp>
00044
00045 namespace Opm {
00046
00054 template <class Scalar>
00055 class MMPCAuxConstraint
00056 {
00057 public:
00058 MMPCAuxConstraint()
00059 {}
00060
00061 MMPCAuxConstraint(unsigned phaseIndex, unsigned compIndex, Scalar val)
00062 : phaseIdx_(phaseIndex)
00063 , compIdx_(compIndex)
00064 , value_(val)
00065 {}
00066
00070 void set(unsigned phaseIndex, unsigned compIndex, Scalar val)
00071 {
00072 phaseIdx_ = phaseIndex;
00073 compIdx_ = compIndex;
00074 value_ = val;
00075 }
00076
00081 unsigned phaseIdx() const
00082 { return phaseIdx_; }
00083
00088 unsigned compIdx() const
00089 { return compIdx_; }
00090
00095 Scalar value() const
00096 { return value_; }
00097
00098 private:
00099 unsigned phaseIdx_;
00100 unsigned compIdx_;
00101 Scalar value_;
00102 };
00103
00127 template <class Scalar, class FluidSystem, class Evaluation = Scalar>
00128 class MiscibleMultiPhaseComposition
00129 {
00130 static const int numPhases = FluidSystem::numPhases;
00131 static const int numComponents = FluidSystem::numComponents;
00132
00133 typedef MathToolbox<Evaluation> Toolbox;
00134
00135 static_assert(numPhases <= numComponents,
00136 "This solver requires that the number fluid phases is smaller or equal "
00137 "to the number of components");
00138
00139
00140 public:
00164 template <class FluidState, class ParameterCache>
00165 static void solve(FluidState& fluidState,
00166 ParameterCache& paramCache,
00167 int phasePresence,
00168 const MMPCAuxConstraint<Evaluation>* auxConstraints,
00169 unsigned numAuxConstraints,
00170 bool setViscosity,
00171 bool setInternalEnergy)
00172 {
00173 static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
00174 "The scalar type of the fluid state must be 'Evaluation'");
00175
00176 #ifndef NDEBUG
00177
00178
00179
00180
00181 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00182 assert(FluidSystem::isIdealMixture(phaseIdx));
00183 }
00184 #endif
00185
00186
00187 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00188 paramCache.updatePhase(fluidState, phaseIdx);
00189
00190
00191
00192
00193 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00194 Evaluation fugCoeff = Opm::decay<Evaluation>(
00195 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
00196 fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
00197 }
00198 }
00199
00200
00201
00202 static const int numEq = numComponents*numPhases;
00203 Dune::FieldMatrix<Evaluation, numEq, numEq> M(0.0);
00204 Dune::FieldVector<Evaluation, numEq> x(0.0);
00205 Dune::FieldVector<Evaluation, numEq> b(0.0);
00206
00207
00208
00209 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00210 const Evaluation& entryCol1 =
00211 fluidState.fugacityCoefficient(0, compIdx)
00212 *fluidState.pressure(0);
00213 unsigned col1Idx = compIdx;
00214
00215 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
00216 unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
00217 unsigned col2Idx = phaseIdx*numComponents + compIdx;
00218
00219 const Evaluation& entryCol2 =
00220 fluidState.fugacityCoefficient(phaseIdx, compIdx)
00221 *fluidState.pressure(phaseIdx);
00222
00223 M[rowIdx][col1Idx] = entryCol1;
00224 M[rowIdx][col2Idx] = -entryCol2;
00225 }
00226 }
00227
00228
00229
00230
00231 unsigned presentPhases = 0;
00232 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00233 if (!(phasePresence& (1 << phaseIdx)))
00234 continue;
00235
00236 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
00237 presentPhases += 1;
00238
00239 b[rowIdx] = 1.0;
00240 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00241 unsigned colIdx = phaseIdx*numComponents + compIdx;
00242
00243 M[rowIdx][colIdx] = 1.0;
00244 }
00245 }
00246
00247 assert(presentPhases + numAuxConstraints == numComponents);
00248
00249
00250 for (unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
00251 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
00252 b[rowIdx] = auxConstraints[auxEqIdx].value();
00253
00254 unsigned colIdx = auxConstraints[auxEqIdx].phaseIdx()*numComponents + auxConstraints[auxEqIdx].compIdx();
00255 M[rowIdx][colIdx] = 1.0;
00256 }
00257
00258
00259 try {
00260 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-50);
00261 M.solve(x, b);
00262 }
00263 catch (const Dune::FMatrixError& e) {
00264 OPM_THROW(NumericalProblem,
00265 "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() << "; M="<<M);
00266 }
00267 catch (...) {
00268 throw;
00269 }
00270
00271
00272
00273
00274 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00275 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00276 unsigned rowIdx = phaseIdx*numComponents + compIdx;
00277 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
00278 }
00279 paramCache.updateComposition(fluidState, phaseIdx);
00280
00281 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
00282 fluidState.setDensity(phaseIdx, rho);
00283
00284 if (setViscosity) {
00285 const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
00286 fluidState.setViscosity(phaseIdx, mu);
00287 }
00288
00289 if (setInternalEnergy) {
00290 const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
00291 fluidState.setEnthalpy(phaseIdx, h);
00292 }
00293 }
00294 }
00295
00303 template <class FluidState, class ParameterCache>
00304 static void solve(FluidState& fluidState,
00305 ParameterCache& paramCache,
00306 bool setViscosity,
00307 bool setInternalEnergy)
00308 {
00309 solve(fluidState,
00310 paramCache,
00311 0xffffff,
00312 0,
00313 0,
00314 setViscosity,
00315 setInternalEnergy);
00316 }
00317 };
00318
00319 }
00320
00321 #endif