74 static const int numPhases = FluidSystem::numPhases;
75 static const int numComponents = FluidSystem::numComponents;
76 static_assert(numPhases == numComponents,
77 "Immiscibility assumes that the number of phases"
78 " is equal to the number of components");
85 static const int numEq = numPhases;
91 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
93 const Dune::FieldVector<Evaluation, numComponents>& )
95 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
97 fluidState.setPressure(phaseIdx, 1e5);
100 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
110 template <
class MaterialLaw,
class Flu
idState>
111 static void solve(FluidState& fluidState,
112 const typename MaterialLaw::Params& matParams,
113 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
114 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
115 Scalar tolerance = -1)
117 typedef typename FluidState::Scalar InputEval;
122 bool allIncompressible =
true;
123 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
124 if (FluidSystem::isCompressible(phaseIdx)) {
125 allIncompressible =
false;
130 if (allIncompressible) {
134 paramCache.updateAll(fluidState);
135 solveAllIncompressible_(fluidState, paramCache, globalMolarities);
139 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
140 typedef Dune::FieldVector<InputEval, numEq> Vector;
143 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
147 tolerance = std::min<Scalar>(1e-5,
148 1e8*std::numeric_limits<Scalar>::epsilon());
150 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
151 flashParamCache.assignPersistentData(paramCache);
164 Valgrind::SetUndefined(J);
165 Valgrind::SetUndefined(deltaX);
166 Valgrind::SetUndefined(b);
168 FlashFluidState flashFluidState;
169 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
174 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
175 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
176 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
178 FlashDefectVector defect;
179 const unsigned nMax = 50;
180 for (
unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
182 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
183 Valgrind::CheckDefined(defect);
187 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
188 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
189 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
191 b[eqIdx] = defect[eqIdx].value();
193 Valgrind::CheckDefined(J);
194 Valgrind::CheckDefined(b);
199 try { J.solve(deltaX, b); }
200 catch (
const Dune::FMatrixError& e) {
203 Valgrind::CheckDefined(deltaX);
206 Scalar relError = update_<MaterialLaw>(flashFluidState, flashParamCache, matParams, deltaX);
208 if (relError < tolerance) {
209 assignOutputFluidState_(flashFluidState, fluidState);
214 std::string msg =
"ImmiscibleFlash solver failed: "
215 "{c_alpha^kappa} = {";
216 for (
const auto& v : globalMolarities)
217 msg +=
" " + std::to_string(v);
219 msg += std::to_string(fluidState.temperature(0));
225 template <
class MaterialLaw,
class InputFlu
idState,
class FlashFlu
idState>
226 static void assignFlashFluidState_(
const InputFluidState& inputFluidState,
227 FlashFluidState& flashFluidState,
228 const typename MaterialLaw::Params& matParams,
229 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
231 typedef typename FlashFluidState::Scalar FlashEval;
238 flashFluidState.setTemperature(inputFluidState.temperature(0));
242 FlashEval Slast = 1.0;
243 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
244 FlashEval S = inputFluidState.saturation(phaseIdx);
245 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
249 flashFluidState.setSaturation(phaseIdx, S);
251 flashFluidState.setSaturation(numPhases - 1, Slast);
255 FlashEval p0 = inputFluidState.pressure(0);
256 p0.setDerivative(p0PvIdx, 1.0);
258 std::array<FlashEval, numPhases> pc;
259 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
260 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
261 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
263 flashParamCache.updateAll(flashFluidState);
266 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
267 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
268 flashFluidState.setDensity(phaseIdx, rho);
272 template <
class FlashFlu
idState,
class OutputFlu
idState>
273 static void assignOutputFluidState_(
const FlashFluidState& flashFluidState,
274 OutputFluidState& outputFluidState)
276 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
279 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
280 const auto& S = flashFluidState.saturation(phaseIdx).value();
281 outputFluidState.setSaturation(phaseIdx, S);
283 const auto& p = flashFluidState.pressure(phaseIdx).value();
284 outputFluidState.setPressure(phaseIdx, p);
286 const auto& rho = flashFluidState.density(phaseIdx).value();
287 outputFluidState.setDensity(phaseIdx, rho);
291 template <
class Flu
idState>
292 static void solveAllIncompressible_(FluidState& fluidState,
293 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
294 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities)
296 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
297 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
298 fluidState.setDensity(phaseIdx, rho);
301 globalMolarities[phaseIdx]
302 / fluidState.molarDensity(phaseIdx);
303 fluidState.setSaturation(phaseIdx, saturation);
307 template <
class Flu
idState,
class FlashDefectVector,
class FlashComponentVector>
308 static void evalDefect_(FlashDefectVector& b,
309 const FluidState& fluidState,
310 const FlashComponentVector& globalMolarities)
313 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
315 fluidState.saturation(phaseIdx)
316 * fluidState.molarity(phaseIdx, phaseIdx);
317 b[phaseIdx] -= globalMolarities[phaseIdx];
321 template <
class MaterialLaw,
class FlashFlu
idState,
class EvalVector>
322 static Scalar update_(FlashFluidState& fluidState,
323 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
324 const typename MaterialLaw::Params& matParams,
325 const EvalVector& deltaX)
328 typedef typename FlashFluidState::Scalar FlashEval;
329 typedef MathToolbox<FlashEval> FlashEvalToolbox;
331 typedef typename FlashEvalToolbox::ValueType InnerEval;
335 assert(deltaX.dimension == numEq);
336 for (
unsigned i = 0; i < numEq; ++i)
337 assert(std::isfinite(scalarValue(deltaX[i])));
341 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
342 FlashEval tmp = getQuantity_(fluidState, pvIdx);
343 InnerEval delta = deltaX[pvIdx];
345 relError = std::max(relError,
346 std::abs(scalarValue(delta))
347 * quantityWeight_(fluidState, pvIdx));
349 if (isSaturationIdx_(pvIdx)) {
352 delta = min(0.25, max(-0.25, delta));
354 else if (isPressureIdx_(pvIdx)) {
357 delta = min(0.5*fluidState.pressure(0).value(),
358 max(-0.50*fluidState.pressure(0).value(), delta));
362 setQuantity_(fluidState, pvIdx, tmp);
365 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
370 template <
class MaterialLaw,
class FlashFlu
idState>
371 static void completeFluidState_(FlashFluidState& flashFluidState,
372 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
373 const typename MaterialLaw::Params& matParams)
375 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
377 typedef typename FlashFluidState::Scalar FlashEval;
381 FlashEval sumSat = 0.0;
382 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
383 sumSat += flashFluidState.saturation(phaseIdx);
384 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
389 std::array<FlashEval, numPhases> pC;
390 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
391 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
392 flashFluidState.setPressure(phaseIdx,
393 flashFluidState.pressure(0)
394 + (pC[phaseIdx] - pC[0]));
397 paramCache.updateAll(flashFluidState, ParamCache::Temperature|ParamCache::Composition);
400 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
401 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
402 flashFluidState.setDensity(phaseIdx, rho);
406 static bool isPressureIdx_(
unsigned pvIdx)
407 {
return pvIdx == 0; }
409 static bool isSaturationIdx_(
unsigned pvIdx)
410 {
return 1 <= pvIdx && pvIdx < numPhases; }
413 template <
class Flu
idState>
414 static const typename FluidState::Scalar& getQuantity_(
const FluidState& fluidState,
unsigned pvIdx)
416 assert(pvIdx < numEq);
420 unsigned phaseIdx = 0;
421 return fluidState.pressure(phaseIdx);
425 assert(pvIdx < numPhases);
426 unsigned phaseIdx = pvIdx - 1;
427 return fluidState.saturation(phaseIdx);
432 template <
class Flu
idState>
433 static void setQuantity_(FluidState& fluidState,
435 const typename FluidState::Scalar& value)
437 assert(pvIdx < numEq);
441 unsigned phaseIdx = 0;
442 fluidState.setPressure(phaseIdx, value);
446 assert(pvIdx < numPhases);
447 unsigned phaseIdx = pvIdx - 1;
448 fluidState.setSaturation(phaseIdx, value);
452 template <
class Flu
idState>
453 static Scalar quantityWeight_(
const FluidState& ,
unsigned pvIdx)
460 assert(pvIdx < numPhases);