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_IMMISCIBLE_FLASH_HPP
00028 #define OPM_IMMISCIBLE_FLASH_HPP
00029
00030 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
00031 #include <opm/material/densead/Evaluation.hpp>
00032 #include <opm/material/densead/Math.hpp>
00033 #include <opm/material/common/MathToolbox.hpp>
00034 #include <opm/common/Valgrind.hpp>
00035
00036 #include <opm/common/ErrorMacros.hpp>
00037 #include <opm/common/Exceptions.hpp>
00038
00039 #include <dune/common/fvector.hh>
00040 #include <dune/common/fmatrix.hh>
00041
00042 #include <limits>
00043 #include <iostream>
00044
00045 namespace Opm {
00046
00073 template <class Scalar, class FluidSystem>
00074 class ImmiscibleFlash
00075 {
00076 static const int numPhases = FluidSystem::numPhases;
00077 static const int numComponents = FluidSystem::numComponents;
00078 static_assert(numPhases == numComponents,
00079 "Immiscibility assumes that the number of phases"
00080 " is equal to the number of components");
00081
00082 enum {
00083 p0PvIdx = 0,
00084 S0PvIdx = 1
00085 };
00086
00087 static const int numEq = numPhases;
00088
00089 public:
00093 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00094 static void guessInitial(FluidState& fluidState,
00095 const Dune::FieldVector<Evaluation, numComponents>& )
00096 {
00097 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
00098
00099 fluidState.setPressure(phaseIdx, 1e5);
00100
00101
00102 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
00103 }
00104 }
00105
00112 template <class MaterialLaw, class FluidState>
00113 static void solve(FluidState& fluidState,
00114 const typename MaterialLaw::Params& matParams,
00115 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
00116 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
00117 Scalar tolerance = -1)
00118 {
00119 typedef typename FluidState::Scalar InputEval;
00120
00122
00124 bool allIncompressible = true;
00125 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00126 if (FluidSystem::isCompressible(phaseIdx)) {
00127 allIncompressible = false;
00128 break;
00129 }
00130 }
00131
00132 if (allIncompressible) {
00133
00134
00135
00136 paramCache.updateAll(fluidState);
00137 solveAllIncompressible_(fluidState, paramCache, globalMolarities);
00138 return;
00139 }
00140
00141 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
00142 typedef Dune::FieldVector<InputEval, numEq> Vector;
00143
00144 typedef Opm::DenseAd::Evaluation<InputEval, numEq> FlashEval;
00145 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
00146 typedef Opm::ImmiscibleFluidState<FlashEval, FluidSystem> FlashFluidState;
00147
00148 Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
00149
00150 if (tolerance <= 0)
00151 tolerance = std::min<Scalar>(1e-5,
00152 1e8*std::numeric_limits<Scalar>::epsilon());
00153
00154 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
00155 flashParamCache.assignPersistentData(paramCache);
00156
00158
00160
00161
00162 Matrix J;
00163
00164 Vector deltaX;
00165
00166 Vector b;
00167
00168 Valgrind::SetUndefined(J);
00169 Valgrind::SetUndefined(deltaX);
00170 Valgrind::SetUndefined(b);
00171
00172 FlashFluidState flashFluidState;
00173 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
00174
00175
00176
00177
00178 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
00179 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
00180 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
00181
00182 FlashDefectVector defect;
00183 const unsigned nMax = 50;
00184 for (unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
00185
00186 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
00187 Valgrind::CheckDefined(defect);
00188
00189
00190
00191 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
00192 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
00193 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
00194
00195 b[eqIdx] = defect[eqIdx].value();
00196 }
00197 Valgrind::CheckDefined(J);
00198 Valgrind::CheckDefined(b);
00199
00200
00201 deltaX = 0;
00202
00203 try { J.solve(deltaX, b); }
00204 catch (Dune::FMatrixError e) {
00205 throw Opm::NumericalProblem(e.what());
00206 }
00207 Valgrind::CheckDefined(deltaX);
00208
00209
00210 Scalar relError = update_<MaterialLaw>(flashFluidState, flashParamCache, matParams, deltaX);
00211
00212 if (relError < tolerance) {
00213 assignOutputFluidState_(flashFluidState, fluidState);
00214 return;
00215 }
00216 }
00217
00218 OPM_THROW(Opm::NumericalProblem,
00219 "ImmiscibleFlash solver failed: "
00220 "{c_alpha^kappa} = {" << globalMolarities << "}, "
00221 << "T = " << fluidState.temperature(0));
00222 }
00223
00224
00225 protected:
00226 template <class FluidState>
00227 static void printFluidState_(const FluidState& fs)
00228 {
00229 std::cout << "saturations: ";
00230 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00231 std::cout << fs.saturation(phaseIdx) << " ";
00232 std::cout << "\n";
00233
00234 std::cout << "pressures: ";
00235 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00236 std::cout << fs.pressure(phaseIdx) << " ";
00237 std::cout << "\n";
00238
00239 std::cout << "densities: ";
00240 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00241 std::cout << fs.density(phaseIdx) << " ";
00242 std::cout << "\n";
00243
00244 std::cout << "global component molarities: ";
00245 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00246 Scalar sum = 0;
00247 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00248 sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
00249 }
00250 std::cout << sum << " ";
00251 }
00252 std::cout << "\n";
00253 }
00254
00255 template <class MaterialLaw, class InputFluidState, class FlashFluidState>
00256 static void assignFlashFluidState_(const InputFluidState& inputFluidState,
00257 FlashFluidState& flashFluidState,
00258 const typename MaterialLaw::Params& matParams,
00259 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
00260 {
00261 typedef typename FlashFluidState::Scalar FlashEval;
00262
00263
00264
00265
00266
00267
00268 flashFluidState.setTemperature(inputFluidState.temperature(0));
00269
00270
00271
00272 FlashEval Slast = 1.0;
00273 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
00274 FlashEval S = inputFluidState.saturation(phaseIdx);
00275 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
00276
00277 Slast -= S;
00278
00279 flashFluidState.setSaturation(phaseIdx, S);
00280 }
00281 flashFluidState.setSaturation(numPhases - 1, Slast);
00282
00283
00284
00285 FlashEval p0 = inputFluidState.pressure(0);
00286 p0.setDerivative(p0PvIdx, 1.0);
00287
00288 std::array<FlashEval, numPhases> pc;
00289 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
00290 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00291 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
00292
00293 flashParamCache.updateAll(flashFluidState);
00294
00295
00296 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00297 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
00298 flashFluidState.setDensity(phaseIdx, rho);
00299 }
00300 }
00301
00302 template <class FlashFluidState, class OutputFluidState>
00303 static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
00304 OutputFluidState& outputFluidState)
00305 {
00306 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
00307
00308
00309 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00310 const auto& S = flashFluidState.saturation(phaseIdx).value();
00311 outputFluidState.setSaturation(phaseIdx, S);
00312
00313 const auto& p = flashFluidState.pressure(phaseIdx).value();
00314 outputFluidState.setPressure(phaseIdx, p);
00315
00316 const auto& rho = flashFluidState.density(phaseIdx).value();
00317 outputFluidState.setDensity(phaseIdx, rho);
00318 }
00319 }
00320
00321 template <class FluidState>
00322 static void solveAllIncompressible_(FluidState& fluidState,
00323 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
00324 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities)
00325 {
00326 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00327 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
00328 fluidState.setDensity(phaseIdx, rho);
00329
00330 Scalar saturation =
00331 globalMolarities[phaseIdx]
00332 / fluidState.molarDensity(phaseIdx);
00333 fluidState.setSaturation(phaseIdx, saturation);
00334 }
00335 }
00336
00337 template <class FluidState, class FlashDefectVector, class FlashComponentVector>
00338 static void evalDefect_(FlashDefectVector& b,
00339 const FluidState& fluidState,
00340 const FlashComponentVector& globalMolarities)
00341 {
00342
00343 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00344 b[phaseIdx] =
00345 fluidState.saturation(phaseIdx)
00346 * fluidState.molarity(phaseIdx, phaseIdx);
00347 b[phaseIdx] -= globalMolarities[phaseIdx];
00348 }
00349 }
00350
00351 template <class MaterialLaw, class FlashFluidState, class EvalVector>
00352 static Scalar update_(FlashFluidState& fluidState,
00353 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
00354 const typename MaterialLaw::Params& matParams,
00355 const EvalVector& deltaX)
00356 {
00357
00358 typedef typename FlashFluidState::Scalar FlashEval;
00359 typedef Opm::MathToolbox<FlashEval> FlashEvalToolbox;
00360
00361 typedef typename FlashEvalToolbox::ValueType InnerEval;
00362
00363 #ifndef NDEBUG
00364
00365 assert(deltaX.dimension == numEq);
00366 for (unsigned i = 0; i < numEq; ++i)
00367 assert(std::isfinite(Opm::scalarValue(deltaX[i])));
00368 #endif
00369
00370 Scalar relError = 0;
00371 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
00372 FlashEval tmp = getQuantity_(fluidState, pvIdx);
00373 InnerEval delta = deltaX[pvIdx];
00374
00375 relError = std::max(relError,
00376 std::abs(Opm::scalarValue(delta))
00377 * quantityWeight_(fluidState, pvIdx));
00378
00379 if (isSaturationIdx_(pvIdx)) {
00380
00381
00382 delta = Opm::min(0.25, Opm::max(-0.25, delta));
00383 }
00384 else if (isPressureIdx_(pvIdx)) {
00385
00386
00387 delta = Opm::min(0.5*fluidState.pressure(0).value(),
00388 Opm::max(-0.50*fluidState.pressure(0).value(), delta));
00389 }
00390
00391 tmp -= delta;
00392 setQuantity_(fluidState, pvIdx, tmp);
00393 }
00394
00395 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
00396
00397 return relError;
00398 }
00399
00400 template <class MaterialLaw, class FlashFluidState>
00401 static void completeFluidState_(FlashFluidState& flashFluidState,
00402 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
00403 const typename MaterialLaw::Params& matParams)
00404 {
00405 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
00406
00407 typedef typename FlashFluidState::Scalar FlashEval;
00408
00409
00410
00411 FlashEval sumSat = 0.0;
00412 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
00413 sumSat += flashFluidState.saturation(phaseIdx);
00414 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
00415
00416
00417
00418
00419 Dune::FieldVector<FlashEval, numPhases> pC;
00420 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
00421 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
00422 flashFluidState.setPressure(phaseIdx,
00423 flashFluidState.pressure(0)
00424 + (pC[phaseIdx] - pC[0]));
00425
00426
00427 paramCache.updateAll(flashFluidState, ParamCache::Temperature|ParamCache::Composition);
00428
00429
00430 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
00431 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
00432 flashFluidState.setDensity(phaseIdx, rho);
00433 }
00434 }
00435
00436 static bool isPressureIdx_(unsigned pvIdx)
00437 { return pvIdx == 0; }
00438
00439 static bool isSaturationIdx_(unsigned pvIdx)
00440 { return 1 <= pvIdx && pvIdx < numPhases; }
00441
00442
00443 template <class FluidState>
00444 static const typename FluidState::Scalar& getQuantity_(const FluidState& fluidState, unsigned pvIdx)
00445 {
00446 assert(pvIdx < numEq);
00447
00448
00449 if (pvIdx < 1) {
00450 unsigned phaseIdx = 0;
00451 return fluidState.pressure(phaseIdx);
00452 }
00453
00454 else {
00455 assert(pvIdx < numPhases);
00456 unsigned phaseIdx = pvIdx - 1;
00457 return fluidState.saturation(phaseIdx);
00458 }
00459 }
00460
00461
00462 template <class FluidState>
00463 static void setQuantity_(FluidState& fluidState,
00464 unsigned pvIdx,
00465 const typename FluidState::Scalar& value)
00466 {
00467 assert(pvIdx < numEq);
00468
00469
00470 if (pvIdx < 1) {
00471 unsigned phaseIdx = 0;
00472 fluidState.setPressure(phaseIdx, value);
00473 }
00474
00475 else {
00476 assert(pvIdx < numPhases);
00477 unsigned phaseIdx = pvIdx - 1;
00478 fluidState.setSaturation(phaseIdx, value);
00479 }
00480 }
00481
00482 template <class FluidState>
00483 static Scalar quantityWeight_(const FluidState& , unsigned pvIdx)
00484 {
00485
00486 if (pvIdx < 1)
00487 return 1e-8;
00488
00489 else {
00490 assert(pvIdx < numPhases);
00491 return 1.0;
00492 }
00493 }
00494 };
00495
00496 }
00497
00498 #endif