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_NCP_FLASH_HPP
00028 #define OPM_NCP_FLASH_HPP
00029
00030 #include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
00031 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
00032 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
00033 #include <opm/material/densead/Evaluation.hpp>
00034 #include <opm/material/densead/Math.hpp>
00035 #include <opm/material/common/MathToolbox.hpp>
00036 #include <opm/common/Valgrind.hpp>
00037
00038 #include <opm/common/ErrorMacros.hpp>
00039 #include <opm/common/Exceptions.hpp>
00040
00041 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00042 #include <dune/common/fvector.hh>
00043 #include <dune/common/fmatrix.hh>
00044 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00045
00046 #include <limits>
00047 #include <iostream>
00048
00049 namespace Opm {
00050
00090 template <class Scalar, class FluidSystem>
00091 class NcpFlash
00092 {
00093 enum { numPhases = FluidSystem::numPhases };
00094 enum { numComponents = FluidSystem::numComponents };
00095
00096 enum {
00097 p0PvIdx = 0,
00098 S0PvIdx = 1,
00099 x00PvIdx = S0PvIdx + numPhases - 1
00100 };
00101
00102 static const int numEq = numPhases*(numComponents + 1);
00103
00104 public:
00108 template <class FluidState, class Evaluation = typename FluidState::Scalar>
00109 static void guessInitial(FluidState& fluidState,
00110 const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
00111 {
00112
00113 Evaluation sumMoles = 0;
00114 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
00115 sumMoles += globalMolarities[compIdx];
00116
00117 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
00118
00119 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
00120 fluidState.setMoleFraction(phaseIdx,
00121 compIdx,
00122 globalMolarities[compIdx]/sumMoles);
00123
00124
00125 fluidState.setPressure(phaseIdx, 1.0135e5);
00126
00127
00128 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
00129 }
00130
00131
00132 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
00133 paramCache.updateAll(fluidState);
00134 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
00135 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
00136 const typename FluidState::Scalar phi =
00137 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
00138 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
00139 }
00140 }
00141 }
00142
00149 template <class MaterialLaw, class FluidState>
00150 static void solve(FluidState& fluidState,
00151 const typename MaterialLaw::Params& matParams,
00152 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
00153 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
00154 Scalar tolerance = -1.0)
00155 {
00156 typedef typename FluidState::Scalar InputEval;
00157
00158 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
00159 typedef Dune::FieldVector<InputEval, numEq> Vector;
00160
00161 typedef Opm::DenseAd::Evaluation<InputEval,
00162 numEq> FlashEval;
00163
00164 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
00165 typedef Opm::CompositionalFluidState<FlashEval, FluidSystem, false> FlashFluidState;
00166
00167 Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
00168
00169 if (tolerance <= 0)
00170 tolerance = std::min<Scalar>(1e-3,
00171 1e8*std::numeric_limits<Scalar>::epsilon());
00172
00173 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
00174 flashParamCache.assignPersistentData(paramCache);
00175
00177
00179
00180
00181 Matrix J;
00182
00183 Vector deltaX;
00184
00185 Vector b;
00186
00187 Valgrind::SetUndefined(J);
00188 Valgrind::SetUndefined(deltaX);
00189 Valgrind::SetUndefined(b);
00190
00191 FlashFluidState flashFluidState;
00192 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
00193
00194
00195
00196
00197 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
00198 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
00199 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
00200
00201 FlashDefectVector defect;
00202 const unsigned nMax = 50;
00203 for (unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
00204
00205 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
00206 Valgrind::CheckDefined(defect);
00207
00208
00209
00210 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
00211 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
00212 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
00213
00214 b[eqIdx] = defect[eqIdx].value();
00215 }
00216 Valgrind::CheckDefined(J);
00217 Valgrind::CheckDefined(b);
00218
00219
00220 deltaX = 0.0;
00221 try { J.solve(deltaX, b); }
00222 catch (const Dune::FMatrixError& e) {
00223 throw Opm::NumericalProblem(e.what());
00224 }
00225 Valgrind::CheckDefined(deltaX);
00226
00227
00228 Scalar relError = update_<MaterialLaw>(flashFluidState, matParams, flashParamCache, deltaX);
00229
00230 if (relError < tolerance) {
00231 assignOutputFluidState_(flashFluidState, fluidState);
00232 return;
00233 }
00234 }
00235
00236 OPM_THROW(NumericalProblem,
00237 "NcpFlash solver failed: "
00238 "{c_alpha^kappa} = {" << globalMolarities << "}, "
00239 << "T = " << fluidState.temperature(0));
00240 }
00241
00249 template <class FluidState, class ComponentVector>
00250 static void solve(FluidState& fluidState,
00251 const ComponentVector& globalMolarities,
00252 Scalar tolerance = 0.0)
00253 {
00254 typedef NullMaterialTraits<Scalar, numPhases> MaterialTraits;
00255 typedef NullMaterial<MaterialTraits> MaterialLaw;
00256 typedef typename MaterialLaw::Params MaterialLawParams;
00257
00258 MaterialLawParams matParams;
00259 solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
00260 }
00261
00262
00263 protected:
00264 template <class FluidState>
00265 static void printFluidState_(const FluidState& fluidState)
00266 {
00267 typedef typename FluidState::Scalar FsScalar;
00268
00269 std::cout << "saturations: ";
00270 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00271 std::cout << fluidState.saturation(phaseIdx) << " ";
00272 std::cout << "\n";
00273
00274 std::cout << "pressures: ";
00275 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00276 std::cout << fluidState.pressure(phaseIdx) << " ";
00277 std::cout << "\n";
00278
00279 std::cout << "densities: ";
00280 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00281 std::cout << fluidState.density(phaseIdx) << " ";
00282 std::cout << "\n";
00283
00284 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00285 std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
00286 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00287 std::cout << fluidState.moleFraction(phaseIdx, compIdx) << " ";
00288 }
00289 std::cout << "\n";
00290 }
00291
00292 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00293 std::cout << "fugacities " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
00294 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00295 std::cout << fluidState.fugacity(phaseIdx, compIdx) << " ";
00296 }
00297 std::cout << "\n";
00298 }
00299
00300 std::cout << "global component molarities: ";
00301 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00302 FsScalar sum = 0;
00303 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00304 sum += fluidState.saturation(phaseIdx)*fluidState.molarity(phaseIdx, compIdx);
00305 }
00306 std::cout << sum << " ";
00307 }
00308 std::cout << "\n";
00309 }
00310
00311 template <class MaterialLaw, class InputFluidState, class FlashFluidState>
00312 static void assignFlashFluidState_(const InputFluidState& inputFluidState,
00313 FlashFluidState& flashFluidState,
00314 const typename MaterialLaw::Params& matParams,
00315 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
00316 {
00317 typedef typename FlashFluidState::Scalar FlashEval;
00318
00319
00320
00321
00322
00323
00324 flashFluidState.setTemperature(inputFluidState.temperature(0));
00325
00326
00327
00328 FlashEval Slast = 1.0;
00329 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
00330 FlashEval S = inputFluidState.saturation(phaseIdx);
00331 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
00332
00333 Slast -= S;
00334
00335 flashFluidState.setSaturation(phaseIdx, S);
00336 }
00337 flashFluidState.setSaturation(numPhases - 1, Slast);
00338
00339
00340
00341 FlashEval p0 = inputFluidState.pressure(0);
00342 p0.setDerivative(p0PvIdx, 1.0);
00343
00344 std::array<FlashEval, numPhases> pc;
00345 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
00346 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
00347 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
00348
00349
00350 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00351 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00352 FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
00353 x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
00354 flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
00355 }
00356 }
00357
00358 flashParamCache.updateAll(flashFluidState);
00359
00360
00361
00362 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00363 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
00364 flashFluidState.setDensity(phaseIdx, rho);
00365
00366 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00367 const FlashEval& fugCoeff = FluidSystem::fugacityCoefficient(flashFluidState, flashParamCache, phaseIdx, compIdx);
00368 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
00369 }
00370 }
00371 }
00372
00373 template <class FlashFluidState, class OutputFluidState>
00374 static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
00375 OutputFluidState& outputFluidState)
00376 {
00377 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
00378
00379
00380 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00381 const auto& S = flashFluidState.saturation(phaseIdx).value();
00382 outputFluidState.setSaturation(phaseIdx, S);
00383
00384 const auto& p = flashFluidState.pressure(phaseIdx).value();
00385 outputFluidState.setPressure(phaseIdx, p);
00386
00387 const auto& rho = flashFluidState.density(phaseIdx).value();
00388 outputFluidState.setDensity(phaseIdx, rho);
00389 }
00390
00391
00392 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00393 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00394 const auto& moleFrac =
00395 flashFluidState.moleFraction(phaseIdx, compIdx).value();
00396 outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
00397
00398 const auto& fugCoeff =
00399 flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
00400 outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
00401 }
00402 }
00403 }
00404
00405 template <class FlashFluidState, class FlashDefectVector, class FlashComponentVector>
00406 static void evalDefect_(FlashDefectVector& b,
00407 const FlashFluidState& fluidState,
00408 const FlashComponentVector& globalMolarities)
00409 {
00410 typedef typename FlashFluidState::Scalar FlashEval;
00411
00412 unsigned eqIdx = 0;
00413
00414
00415 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00416 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
00417 b[eqIdx] =
00418 fluidState.fugacity(0, compIdx)
00419 - fluidState.fugacity(phaseIdx, compIdx);
00420 ++eqIdx;
00421 }
00422 }
00423 assert(eqIdx == numComponents*(numPhases - 1));
00424
00425
00426
00427
00428
00429 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
00430 b[eqIdx] = 0.0;
00431 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00432 b[eqIdx] +=
00433 fluidState.saturation(phaseIdx)
00434 * fluidState.molarity(phaseIdx, compIdx);
00435 }
00436
00437 b[eqIdx] -= globalMolarities[compIdx];
00438 ++eqIdx;
00439 }
00440
00441
00442 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00443 FlashEval oneMinusSumMoleFrac = 1.0;
00444 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
00445 oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
00446
00447 if (oneMinusSumMoleFrac > fluidState.saturation(phaseIdx))
00448 b[eqIdx] = fluidState.saturation(phaseIdx);
00449 else
00450 b[eqIdx] = oneMinusSumMoleFrac;
00451
00452 ++eqIdx;
00453 }
00454 }
00455
00456 template <class MaterialLaw, class FlashFluidState, class EvalVector>
00457 static Scalar update_(FlashFluidState& fluidState,
00458 const typename MaterialLaw::Params& matParams,
00459 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
00460 const EvalVector& deltaX)
00461 {
00462
00463 typedef typename FlashFluidState::Scalar FlashEval;
00464 typedef typename FlashEval::ValueType InnerEval;
00465
00466 #ifndef NDEBUG
00467
00468 assert(deltaX.dimension == numEq);
00469 for (unsigned i = 0; i < numEq; ++i)
00470 assert(std::isfinite(Opm::scalarValue(deltaX[i])));
00471 #endif
00472
00473 Scalar relError = 0;
00474 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
00475 FlashEval tmp = getQuantity_(fluidState, pvIdx);
00476 InnerEval delta = deltaX[pvIdx];
00477
00478 relError = std::max(relError,
00479 std::abs(Opm::scalarValue(delta))
00480 * quantityWeight_(fluidState, pvIdx));
00481
00482 if (isSaturationIdx_(pvIdx)) {
00483
00484 delta = Opm::min(0.25, Opm::max(-0.25, delta));
00485 }
00486 else if (isMoleFracIdx_(pvIdx)) {
00487
00488 delta = Opm::min(0.20, Opm::max(-0.20, delta));
00489 }
00490 else if (isPressureIdx_(pvIdx)) {
00491
00492 delta = Opm::min(0.5*fluidState.pressure(0).value(),
00493 Opm::max(-0.5*fluidState.pressure(0).value(),
00494 delta));
00495 }
00496
00497 tmp -= delta;
00498 setQuantity_(fluidState, pvIdx, tmp);
00499 }
00500
00501 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
00502
00503 return relError;
00504 }
00505
00506 template <class MaterialLaw, class FlashFluidState>
00507 static void completeFluidState_(FlashFluidState& flashFluidState,
00508 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
00509 const typename MaterialLaw::Params& matParams)
00510 {
00511 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
00512
00513 typedef typename FlashFluidState::Scalar FlashEval;
00514
00515
00516
00517 FlashEval sumSat = 0.0;
00518 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
00519 sumSat += flashFluidState.saturation(phaseIdx);
00520 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
00521
00522
00523
00524
00525 Dune::FieldVector<FlashEval, numPhases> pC;
00526 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
00527 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
00528 flashFluidState.setPressure(phaseIdx,
00529 flashFluidState.pressure(0)
00530 + (pC[phaseIdx] - pC[0]));
00531
00532
00533 paramCache.updateAll(flashFluidState, ParamCache::Temperature);
00534
00535
00536 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
00537 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
00538 flashFluidState.setDensity(phaseIdx, rho);
00539
00540 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
00541 const FlashEval& phi = FluidSystem::fugacityCoefficient(flashFluidState, paramCache, phaseIdx, compIdx);
00542 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
00543 }
00544 }
00545 }
00546
00547 static bool isPressureIdx_(unsigned pvIdx)
00548 { return pvIdx == 0; }
00549
00550 static bool isSaturationIdx_(unsigned pvIdx)
00551 { return 1 <= pvIdx && pvIdx < numPhases; }
00552
00553 static bool isMoleFracIdx_(unsigned pvIdx)
00554 { return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
00555
00556
00557 template <class FluidState>
00558 static const typename FluidState::Scalar& getQuantity_(const FluidState& fluidState, unsigned pvIdx)
00559 {
00560 assert(pvIdx < numEq);
00561
00562
00563 if (pvIdx < 1) {
00564 unsigned phaseIdx = 0;
00565 return fluidState.pressure(phaseIdx);
00566 }
00567
00568 else if (pvIdx < numPhases) {
00569 unsigned phaseIdx = pvIdx - 1;
00570 return fluidState.saturation(phaseIdx);
00571 }
00572
00573 else
00574 {
00575 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
00576 unsigned compIdx = (pvIdx - numPhases)%numComponents;
00577 return fluidState.moleFraction(phaseIdx, compIdx);
00578 }
00579 }
00580
00581
00582 template <class FluidState>
00583 static void setQuantity_(FluidState& fluidState,
00584 unsigned pvIdx,
00585 const typename FluidState::Scalar& value)
00586 {
00587 assert(pvIdx < numEq);
00588
00589 Valgrind::CheckDefined(value);
00590
00591 if (pvIdx < 1) {
00592 unsigned phaseIdx = 0;
00593 fluidState.setPressure(phaseIdx, value);
00594 }
00595
00596 else if (pvIdx < numPhases) {
00597 unsigned phaseIdx = pvIdx - 1;
00598 fluidState.setSaturation(phaseIdx, value);
00599 }
00600
00601 else {
00602 assert(pvIdx < numPhases + numPhases*numComponents);
00603 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
00604 unsigned compIdx = (pvIdx - numPhases)%numComponents;
00605 fluidState.setMoleFraction(phaseIdx, compIdx, value);
00606 }
00607 }
00608
00609 template <class FluidState>
00610 static Scalar quantityWeight_(const FluidState& , unsigned pvIdx)
00611 {
00612
00613 if (pvIdx < 1)
00614 return 1e-6;
00615
00616 else if (pvIdx < numPhases)
00617 return 1.0;
00618
00619 else
00620 return 1.0;
00621 }
00622 };
00623
00624 }
00625
00626 #endif