3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
ncpflash.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_NCP_FLASH_HH
26#define DUMUX_NCP_FLASH_HH
27
28#include <dune/common/fvector.hh>
29#include <dune/common/fmatrix.hh>
30
33
34namespace Dumux {
35
77template <class Scalar, class FluidSystem>
79{
80 enum { numPhases = FluidSystem::numPhases };
81 enum { numComponents = FluidSystem::numComponents };
82
83 using ParameterCache = typename FluidSystem::ParameterCache;
84
85 static constexpr int numEq = numPhases*(numComponents + 1);
86
87 using Matrix = Dune::FieldMatrix<Scalar, numEq, numEq>;
88 using Vector = Dune::FieldVector<Scalar, numEq>;
89
90public:
91 using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
92
97 explicit NcpFlash(int wettingPhaseIdx = 0)
98 : wettingPhaseIdx_(wettingPhaseIdx)
99 {}
100
107 template <class FluidState>
108 void guessInitial(FluidState &fluidState,
109 ParameterCache &paramCache,
110 const ComponentVector &globalMolarities)
111 {
112 // the sum of all molarities
113 Scalar sumMoles = 0;
114 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
115 sumMoles += globalMolarities[compIdx];
116
117 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
118 // composition
119 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
120 fluidState.setMoleFraction(phaseIdx,
121 compIdx,
122 globalMolarities[compIdx]/sumMoles);
123
124 // pressure. use atmospheric pressure as initial guess
125 fluidState.setPressure(phaseIdx, 1.0135e5);
126
127 // saturation. assume all fluids to be equally distributed
128 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
129 }
130
131 // set the fugacity coefficients of all components in all phases
132 paramCache.updateAll(fluidState);
133 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
134 for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
135 Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
136 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
137 }
138 }
139 }
140
151 template <class MaterialLaw, class FluidState>
152 void solve(FluidState &fluidState,
153 ParameterCache &paramCache,
154 const MaterialLaw& material,
155 const ComponentVector &globalMolarities)
156 {
158 // Newton method
160
161 // Jacobian matrix
162 Matrix J;
163 // solution, i.e. phase composition
164 Vector deltaX;
165 // right hand side
166 Vector b;
167
168 // make the fluid state consistent with the fluid system.
169 completeFluidState_(fluidState, paramCache, material);
170
171 /*
172 std::cout << "--------------------\n";
173 std::cout << "globalMolarities: ";
174 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
175 std::cout << globalMolarities[compIdx] << " ";
176 std::cout << "\n";
177 */
178
179 const int nMax = 50; // <- maximum number of newton iterations
180 for (int nIdx = 0; nIdx < nMax; ++nIdx) {
181 // calculate Jacobian matrix and right hand side
182 linearize_(J, b, fluidState, paramCache, material, globalMolarities);
183
184 // Solve J*x = b
185 deltaX = 0;
186
187 try {
188
189 // simple preconditioning of the linear system to reduce condition number
190 int eqIdx = 0;
191 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
192 {
193 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
194 {
195 J[eqIdx] *= 10e-7; // roughly the magnitude of the fugacities
196 b[eqIdx] *= 10e-7;
197 ++eqIdx;
198 }
199 }
200
201 J.solve(deltaX, b);
202
203 }
204 catch (Dune::FMatrixError &e)
205 {
206 /*
207 printFluidState_(fluidState);
208 std::cout << "error: " << e << "\n";
209 std::cout << "b: " << b << "\n";
210 std::cout << "J: " << J << "\n";
211 */
212
213 throw NumericalProblem(e.what());
214 }
215
216 /*
217 printFluidState_(fluidState);
218 std::cout << "J:\n";
219 for (int i = 0; i < numEq; ++i) {
220 for (int j = 0; j < numEq; ++j) {
221 std::ostringstream os;
222 os << J[i][j];
223
224 std::string s(os.str());
225 do {
226 s += " ";
227 } while (s.size() < 20);
228 std::cout << s;
229 }
230 std::cout << "\n";
231 }
232
233 std::cout << "deltaX: " << deltaX << "\n";
234 std::cout << "---------------\n";
235 */
236
237 // update the fluid quantities.
238 Scalar relError = update_(fluidState, paramCache, material, deltaX);
239
240 if (relError < 1e-9)
241 return;
242 }
243
244 /*
245 printFluidState_(fluidState);
246 std::cout << "globalMolarities: ";
247 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
248 std::cout << globalMolarities[compIdx] << " ";
249 std::cout << "\n";
250 */
251
252 DUNE_THROW(NumericalProblem,
253 "Flash calculation failed."
254 " {c_alpha^kappa} = {" << globalMolarities << "}, T = "
255 << fluidState.temperature(/*phaseIdx=*/0));
256 }
257
258
259protected:
260 template <class FluidState>
261 void printFluidState_(const FluidState &fs)
262 {
263 std::cout << "saturations: ";
264 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
265 std::cout << fs.saturation(phaseIdx) << " ";
266 std::cout << "\n";
267
268 std::cout << "pressures: ";
269 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
270 std::cout << fs.pressure(phaseIdx) << " ";
271 std::cout << "\n";
272
273 std::cout << "densities: ";
274 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
275 std::cout << fs.density(phaseIdx) << " ";
276 std::cout << "\n";
277
278 std::cout << "molar densities: ";
279 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
280 std::cout << fs.molarDensity(phaseIdx) << " ";
281 std::cout << "\n";
282
283 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
284 std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
285 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
286 std::cout << fs.moleFraction(phaseIdx, compIdx) << " ";
287 }
288 std::cout << "\n";
289 }
290
291 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
292 std::cout << "fugacities " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
293 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
294 std::cout << fs.fugacity(phaseIdx, compIdx) << " ";
295 }
296 std::cout << "\n";
297 }
298
299 std::cout << "global component molarities: ";
300 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
301 Scalar sum = 0;
302 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
303 sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
304 }
305 std::cout << sum << " ";
306 }
307 std::cout << "\n";
308 }
309
310 template <class MaterialLaw, class FluidState>
311 void linearize_(Matrix &J,
312 Vector &b,
313 FluidState &fluidState,
314 ParameterCache &paramCache,
315 const MaterialLaw& material,
316 const ComponentVector &globalMolarities)
317 {
318 FluidState origFluidState(fluidState);
319 ParameterCache origParamCache(paramCache);
320
321 Vector tmp;
322
323 // reset jacobian
324 J = 0;
325
326 calculateDefect_(b, fluidState, fluidState, globalMolarities);
327
328 // assemble jacobian matrix
329 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
331 // approximately calculate partial derivatives of the
332 // fugacity defect of all components in regard to the mole
333 // fraction of the i-th component. This is done via
334 // forward differences
335
336 // deviate the mole fraction of the i-th component
337 Scalar x_i = getQuantity_(fluidState, pvIdx);
338 const Scalar eps = 1e-8/quantityWeight_(fluidState, pvIdx);
339 setQuantity_(fluidState, paramCache, material, pvIdx, x_i + eps);
340
341 // compute derivative of the defect
342 calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
343 tmp -= b;
344 tmp /= eps;
345
346 // store derivative in jacobian matrix
347 for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
348 J[eqIdx][pvIdx] = tmp[eqIdx];
349
350 // fluid state and parameter cache to their original values
351 fluidState = origFluidState;
352 paramCache = origParamCache;
353
354 // end forward differences
356 }
357 }
358
359 template <class FluidState>
360 void calculateDefect_(Vector &b,
361 const FluidState &fluidStateEval,
362 const FluidState &fluidState,
363 const ComponentVector &globalMolarities)
364 {
365 int eqIdx = 0;
366
367 // fugacity of any component must be equal in all phases
368 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
369 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
370 b[eqIdx] =
371 fluidState.fugacity(/*phaseIdx=*/0, compIdx) -
372 fluidState.fugacity(phaseIdx, compIdx);
373 ++eqIdx;
374 }
375 }
376
377 assert(eqIdx == numComponents*(numPhases - 1));
378
379 // the fact saturations must sum up to 1 is included explicitly!
380
381 // capillary pressures are explicitly included!
382
383 // global molarities are given
384 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
385 b[eqIdx] = 0.0;
386 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
387 b[eqIdx] +=
388 fluidState.saturation(phaseIdx)
389 * fluidState.molarity(phaseIdx, compIdx);
390 }
391
392 b[eqIdx] -= globalMolarities[compIdx];
393 ++eqIdx;
394 }
395
396 // model assumptions (-> non-linear complementarity functions)
397 // must be adhered
398 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
399 Scalar sumMoleFracEval = 0.0;
400 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
401 sumMoleFracEval += fluidStateEval.moleFraction(phaseIdx, compIdx);
402
403 if (1.0 - sumMoleFracEval > fluidStateEval.saturation(phaseIdx)) {
404 b[eqIdx] = fluidState.saturation(phaseIdx);
405 }
406 else {
407 Scalar sumMoleFrac = 0.0;
408 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
409 sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
410 b[eqIdx] = 1.0 - sumMoleFrac;
411 }
412
413 ++eqIdx;
414 }
415 }
416
417 template <class MaterialLaw, class FluidState>
418 Scalar update_(FluidState &fluidState,
419 ParameterCache &paramCache,
420 const MaterialLaw& material,
421 const Vector &deltaX)
422 {
423 Scalar relError = 0;
424 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
425 Scalar tmp = getQuantity_(fluidState, pvIdx);
426 Scalar delta = deltaX[pvIdx];
427 using std::max;
428 using std::abs;
429 using std::min;
430 relError = max(relError, abs(delta)*quantityWeight_(fluidState, pvIdx));
431
432 if (isSaturationIdx_(pvIdx)) {
433 // dampen to at most 20% change in saturation per
434 // iteration
435 delta = min(0.2, max(-0.2, delta));
436 }
437 else if (isMoleFracIdx_(pvIdx)) {
438 // dampen to at most 15% change in mole fraction per
439 // iteration
440 delta = min(0.15, max(-0.15, delta));
441 }
442 else if (isPressureIdx_(pvIdx)) {
443 // dampen to at most 15% change in pressure per
444 // iteration
445 delta = min(0.15*fluidState.pressure(0), max(-0.15*fluidState.pressure(0), delta));
446 }
447
448 setQuantityRaw_(fluidState, pvIdx, tmp - delta);
449 }
450
451 // make sure all saturations, pressures and mole fractions are non-negative
452 Scalar sumSat = 0;
453 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
454 Scalar value = fluidState.saturation(phaseIdx);
455 if (value < -0.05) {
456 value = -0.05;
457 fluidState.setSaturation(phaseIdx, value);
458 }
459 sumSat += value;
460
461 value = fluidState.pressure(phaseIdx);
462 if (value < 0)
463 fluidState.setPressure(phaseIdx, 0.0);
464
465 bool allMoleFractionsAreZero = true;
466 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
467 value = fluidState.moleFraction(phaseIdx, compIdx);
468 if (value > 0)
469 allMoleFractionsAreZero = false;
470 else if (value < 0)
471 fluidState.setMoleFraction(phaseIdx, compIdx, 0.0);
472 }
473
474 // set at least one mole fraction to a positive value
475 // to prevent breakdown of the linear solver in completeFluidState_
476 if (allMoleFractionsAreZero)
477 fluidState.setMoleFraction(phaseIdx, /*compIdx=*/0, 1e-10);
478 }
479
480 // last saturation
481 if (sumSat > 1.05) {
482 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
483 Scalar value = fluidState.saturation(phaseIdx)/(0.95*sumSat);
484 fluidState.setSaturation(phaseIdx, value);
485 }
486 }
487
488 completeFluidState_(fluidState, paramCache, material);
489
490 return relError;
491 }
492
493 template <class MaterialLaw, class FluidState>
494 void completeFluidState_(FluidState &fluidState,
495 ParameterCache &paramCache,
496 const MaterialLaw& material)
497 {
498 // calculate the saturation of the last phase as a function of
499 // the other saturations
500 Scalar sumSat = 0.0;
501 for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
502 sumSat += fluidState.saturation(phaseIdx);
503 fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
504
505 // update the pressures using the material law (saturations
506 // and first pressure are already set because it is implicitly
507 // solved for.)
508 const auto pc = material.capillaryPressures(fluidState, wettingPhaseIdx_);
509 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
510 fluidState.setPressure(phaseIdx,
511 fluidState.pressure(0)
512 + (pc[phaseIdx] - pc[0]));
513
514 // update the parameter cache
515 paramCache.updateAll(fluidState, /*except=*/ParameterCache::Temperature);
516
517 // update all densities and fugacity coefficients
518 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
519
520 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
521 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
522 fluidState.setDensity(phaseIdx, rho);
523 fluidState.setMolarDensity(phaseIdx, rhoMolar);
524
525 for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
526 Scalar phi = FluidSystem::fugacityCoefficient( fluidState, paramCache, phaseIdx, compIdx);
527 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
528 }
529 }
530 }
531
532 bool isPressureIdx_(int pvIdx)
533 { return pvIdx == 0; }
534
535 bool isSaturationIdx_(int pvIdx)
536 { return 1 <= pvIdx && pvIdx < numPhases; }
537
538 bool isMoleFracIdx_(int pvIdx)
539 { return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
540
541 // retrieves a quantity from the fluid state
542 template <class FluidState>
543 Scalar getQuantity_(const FluidState &fs, int pvIdx)
544 {
545 assert(pvIdx < numEq);
546
547 // first pressure
548 if (pvIdx < 1) {
549 int phaseIdx = 0;
550 return fs.pressure(phaseIdx);
551 }
552 // first M - 1 saturations
553 else if (pvIdx < numPhases) {
554 int phaseIdx = pvIdx - 1;
555 return fs.saturation(phaseIdx);
556 }
557 // mole fractions
558 else // if (pvIdx < numPhases + numPhases*numComponents)
559 {
560 int phaseIdx = (pvIdx - numPhases)/numComponents;
561 int compIdx = (pvIdx - numPhases)%numComponents;
562 return fs.moleFraction(phaseIdx, compIdx);
563 }
564 }
565
566 // set a quantity in the fluid state
567 template <class MaterialLaw, class FluidState>
568 void setQuantity_(FluidState &fs,
569 ParameterCache &paramCache,
570 const MaterialLaw& material,
571 int pvIdx,
572 Scalar value)
573 {
574 assert(0 <= pvIdx && pvIdx < numEq);
575
576 if (pvIdx < 1) { // <- first pressure
577 Scalar delta = value - fs.pressure(0);
578
579 // set all pressures. here we assume that the capillary
580 // pressure does not depend on absolute pressure.
581 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
582 fs.setPressure(phaseIdx, fs.pressure(phaseIdx) + delta);
583 paramCache.updateAllPressures(fs);
584
585 // update all densities and fugacity coefficients
586 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
587 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
588 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
589 fs.setDensity(phaseIdx, rho);
590 fs.setMolarDensity(phaseIdx, rhoMolar);
591
592 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
593 Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
594 fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
595 }
596 }
597 }
598 else if (pvIdx < numPhases) { // <- first M - 1 saturations
599 Scalar delta = value - fs.saturation(/*phaseIdx=*/pvIdx - 1);
600 fs.setSaturation(/*phaseIdx=*/pvIdx - 1, value);
601
602 // set last saturation (-> minus the change of the
603 // satuation of the other phase)
604 fs.setSaturation(/*phaseIdx=*/numPhases - 1,
605 fs.saturation(numPhases - 1) - delta);
606
607 // update all fluid pressures using the capillary pressure
608 // law
609 const auto pc = material.capillaryPressures(fs, wettingPhaseIdx_);
610 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
611 fs.setPressure(phaseIdx,
612 fs.pressure(0)
613 + (pc[phaseIdx] - pc[0]));
614 paramCache.updateAllPressures(fs);
615
616 // update all densities and fugacity coefficients
617 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
618 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
619 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
620 fs.setDensity(phaseIdx, rho);
621 fs.setMolarDensity(phaseIdx, rhoMolar);
622
623 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
624 Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
625 fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
626 }
627 }
628 }
629 else if (pvIdx < numPhases + numPhases*numComponents) // <- mole fractions
630 {
631 int phaseIdx = (pvIdx - numPhases)/numComponents;
632 int compIdx = (pvIdx - numPhases)%numComponents;
633
634 fs.setMoleFraction(phaseIdx, compIdx, value);
635 paramCache.updateSingleMoleFraction(fs, phaseIdx, compIdx);
636
637 // update the density of the phase
638 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
639 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
640 fs.setDensity(phaseIdx, rho);
641 fs.setMolarDensity(phaseIdx, rhoMolar);
642
643 // if the phase's fugacity coefficients are composition
644 // dependent, update them as well.
645 if (!FluidSystem::isIdealMixture(phaseIdx)) {
646 for (int compIdx2 = 0; compIdx2 < numComponents; ++compIdx2) {
647 Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx2);
648 fs.setFugacityCoefficient(phaseIdx, compIdx2, phi);
649 }
650 }
651 }
652 else {
653 assert(false);
654 }
655 }
656
657 // set a quantity in the fluid state
658 template <class FluidState>
659 void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
660 {
661 assert(pvIdx < numEq);
662
663 // first pressure
664 if (pvIdx < 1) {
665 int phaseIdx = 0;
666 fs.setPressure(phaseIdx, value);
667 }
668 // first M - 1 saturations
669 else if (pvIdx < numPhases) {
670 int phaseIdx = pvIdx - 1;
671 fs.setSaturation(phaseIdx, value);
672 }
673 // mole fractions
674 else // if (pvIdx < numPhases + numPhases*numComponents)
675 {
676 int phaseIdx = (pvIdx - numPhases)/numComponents;
677 int compIdx = (pvIdx - numPhases)%numComponents;
678 fs.setMoleFraction(phaseIdx, compIdx, value);
679 }
680 }
681
682 template <class FluidState>
683 Scalar quantityWeight_(const FluidState &fs, int pvIdx)
684 {
685 // first pressure
686 if (pvIdx < 1)
687 return 1.0/1e5;
688 // first M - 1 saturations
689 else if (pvIdx < numPhases)
690 return 1.0;
691 // mole fractions
692 else // if (pvIdx < numPhases + numPhases*numComponents)
693 return 1.0;
694 }
695
696private:
697 int wettingPhaseIdx_ = 0;
698};
699
700} // end namespace Dumux
701
702#endif
Some exceptions thrown in DuMux
Makes the capillary pressure-saturation relations available under the M-phase API for material laws.
Definition: adapt.hh:29
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
Determines the phase compositions, pressures and saturations given the total mass of all components.
Definition: ncpflash.hh:79
Dune::FieldVector< Scalar, numComponents > ComponentVector
Definition: ncpflash.hh:91
NcpFlash(int wettingPhaseIdx=0)
Contruct a new flash.
Definition: ncpflash.hh:97
bool isSaturationIdx_(int pvIdx)
Definition: ncpflash.hh:535
Scalar update_(FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material, const Vector &deltaX)
Definition: ncpflash.hh:418
void completeFluidState_(FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material)
Definition: ncpflash.hh:494
bool isPressureIdx_(int pvIdx)
Definition: ncpflash.hh:532
Scalar getQuantity_(const FluidState &fs, int pvIdx)
Definition: ncpflash.hh:543
void linearize_(Matrix &J, Vector &b, FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material, const ComponentVector &globalMolarities)
Definition: ncpflash.hh:311
void guessInitial(FluidState &fluidState, ParameterCache &paramCache, const ComponentVector &globalMolarities)
Guess initial values for all quantities.
Definition: ncpflash.hh:108
void solve(FluidState &fluidState, ParameterCache &paramCache, const MaterialLaw &material, const ComponentVector &globalMolarities)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: ncpflash.hh:152
void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
Definition: ncpflash.hh:659
Scalar quantityWeight_(const FluidState &fs, int pvIdx)
Definition: ncpflash.hh:683
bool isMoleFracIdx_(int pvIdx)
Definition: ncpflash.hh:538
void setQuantity_(FluidState &fs, ParameterCache &paramCache, const MaterialLaw &material, int pvIdx, Scalar value)
Definition: ncpflash.hh:568
void printFluidState_(const FluidState &fs)
Definition: ncpflash.hh:261
void calculateDefect_(Vector &b, const FluidState &fluidStateEval, const FluidState &fluidState, const ComponentVector &globalMolarities)
Definition: ncpflash.hh:360