All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
NcpFlash.hpp
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  This file is part of the Open Porous Media project (OPM).
5 
6  OPM 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 2 of the License, or
9  (at your option) any later version.
10 
11  OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
27 #ifndef OPM_NCP_FLASH_HPP
28 #define OPM_NCP_FLASH_HPP
29 
36 #include <opm/common/Valgrind.hpp>
37 
38 #include <opm/common/ErrorMacros.hpp>
39 #include <opm/common/Exceptions.hpp>
40 
41 #include <opm/common/utility/platform_dependent/disable_warnings.h>
42 #include <dune/common/fvector.hh>
43 #include <dune/common/fmatrix.hh>
44 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
45 
46 #include <limits>
47 #include <iostream>
48 
49 namespace Opm {
50 
90 template <class Scalar, class FluidSystem>
91 class NcpFlash
92 {
93  enum { numPhases = FluidSystem::numPhases };
94  enum { numComponents = FluidSystem::numComponents };
95 
96  enum {
97  p0PvIdx = 0,
98  S0PvIdx = 1,
99  x00PvIdx = S0PvIdx + numPhases - 1
100  };
101 
102  static const int numEq = numPhases*(numComponents + 1);
103 
104 public:
108  template <class FluidState, class Evaluation = typename FluidState::Scalar>
109  static void guessInitial(FluidState& fluidState,
110  const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
111  {
112  // the sum of all molarities
113  Evaluation sumMoles = 0;
114  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
115  sumMoles += globalMolarities[compIdx];
116 
117  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
118  // composition
119  for (unsigned 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  typename FluidSystem::template ParameterCache<Evaluation> paramCache;
133  paramCache.updateAll(fluidState);
134  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
135  for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
136  const typename FluidState::Scalar phi =
137  FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
138  fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
139  }
140  }
141  }
142 
149  template <class MaterialLaw, class FluidState>
150  static void solve(FluidState& fluidState,
151  const typename MaterialLaw::Params& matParams,
152  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
153  const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
154  Scalar tolerance = -1.0)
155  {
156  typedef typename FluidState::Scalar InputEval;
157 
158  typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
159  typedef Dune::FieldVector<InputEval, numEq> Vector;
160 
161  typedef Opm::DenseAd::Evaluation</*Scalar=*/InputEval,
162  /*numDerivs=*/numEq> FlashEval;
163 
164  typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
165  typedef Opm::CompositionalFluidState<FlashEval, FluidSystem, /*energy=*/false> FlashFluidState;
166 
167  Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
168 
169  if (tolerance <= 0)
170  tolerance = std::min<Scalar>(1e-3,
171  1e8*std::numeric_limits<Scalar>::epsilon());
172 
173  typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
174  flashParamCache.assignPersistentData(paramCache);
175 
177  // Newton method
179 
180  // Jacobian matrix
181  Matrix J;
182  // solution, i.e. phase composition
183  Vector deltaX;
184  // right hand side
185  Vector b;
186 
187  Valgrind::SetUndefined(J);
188  Valgrind::SetUndefined(deltaX);
189  Valgrind::SetUndefined(b);
190 
191  FlashFluidState flashFluidState;
192  assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
193 
194  // copy the global molarities to a vector of evaluations. Remember that the
195  // global molarities are constants. (but we need to copy them to a vector of
196  // FlashEvals anyway in order to avoid getting into hell's kitchen.)
197  Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
198  for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
199  flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
200 
201  FlashDefectVector defect;
202  const unsigned nMax = 50; // <- maximum number of newton iterations
203  for (unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
204  // calculate the defect of the flash equations and their derivatives
205  evalDefect_(defect, flashFluidState, flashGlobalMolarities);
206  Valgrind::CheckDefined(defect);
207 
208  // create field matrices and vectors out of the evaluation vector to solve
209  // the linear system of equations.
210  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
211  for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
212  J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
213 
214  b[eqIdx] = defect[eqIdx].value();
215  }
216  Valgrind::CheckDefined(J);
217  Valgrind::CheckDefined(b);
218 
219  // Solve J*x = b
220  deltaX = 0.0;
221  try { J.solve(deltaX, b); }
222  catch (const Dune::FMatrixError& e) {
223  throw Opm::NumericalProblem(e.what());
224  }
225  Valgrind::CheckDefined(deltaX);
226 
227  // update the fluid quantities.
228  Scalar relError = update_<MaterialLaw>(flashFluidState, matParams, flashParamCache, deltaX);
229 
230  if (relError < tolerance) {
231  assignOutputFluidState_(flashFluidState, fluidState);
232  return;
233  }
234  }
235 
236  OPM_THROW(NumericalProblem,
237  "NcpFlash solver failed: "
238  "{c_alpha^kappa} = {" << globalMolarities << "}, "
239  << "T = " << fluidState.temperature(/*phaseIdx=*/0));
240  }
241 
249  template <class FluidState, class ComponentVector>
250  static void solve(FluidState& fluidState,
251  const ComponentVector& globalMolarities,
252  Scalar tolerance = 0.0)
253  {
254  typedef NullMaterialTraits<Scalar, numPhases> MaterialTraits;
255  typedef NullMaterial<MaterialTraits> MaterialLaw;
256  typedef typename MaterialLaw::Params MaterialLawParams;
257 
258  MaterialLawParams matParams;
259  solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
260  }
261 
262 
263 protected:
264  template <class FluidState>
265  static void printFluidState_(const FluidState& fluidState)
266  {
267  typedef typename FluidState::Scalar FsScalar;
268 
269  std::cout << "saturations: ";
270  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
271  std::cout << fluidState.saturation(phaseIdx) << " ";
272  std::cout << "\n";
273 
274  std::cout << "pressures: ";
275  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
276  std::cout << fluidState.pressure(phaseIdx) << " ";
277  std::cout << "\n";
278 
279  std::cout << "densities: ";
280  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
281  std::cout << fluidState.density(phaseIdx) << " ";
282  std::cout << "\n";
283 
284  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
285  std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
286  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
287  std::cout << fluidState.moleFraction(phaseIdx, compIdx) << " ";
288  }
289  std::cout << "\n";
290  }
291 
292  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
293  std::cout << "fugacities " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
294  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
295  std::cout << fluidState.fugacity(phaseIdx, compIdx) << " ";
296  }
297  std::cout << "\n";
298  }
299 
300  std::cout << "global component molarities: ";
301  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
302  FsScalar sum = 0;
303  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
304  sum += fluidState.saturation(phaseIdx)*fluidState.molarity(phaseIdx, compIdx);
305  }
306  std::cout << sum << " ";
307  }
308  std::cout << "\n";
309  }
310 
311  template <class MaterialLaw, class InputFluidState, class FlashFluidState>
312  static void assignFlashFluidState_(const InputFluidState& inputFluidState,
313  FlashFluidState& flashFluidState,
314  const typename MaterialLaw::Params& matParams,
315  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
316  {
317  typedef typename FlashFluidState::Scalar FlashEval;
318 
319  // copy the temperature: even though the model which uses the flash solver might
320  // be non-isothermal, the flash solver does not consider energy. (it could be
321  // modified to do so relatively easily, but it would come at increased
322  // computational cost and normally temperature instead of "total internal energy
323  // of the fluids" is specified.)
324  flashFluidState.setTemperature(inputFluidState.temperature(/*phaseIdx=*/0));
325 
326  // copy the saturations: the first N-1 phases are primary variables, the last one
327  // is one minus the sum of the former.
328  FlashEval Slast = 1.0;
329  for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
330  FlashEval S = inputFluidState.saturation(phaseIdx);
331  S.setDerivative(S0PvIdx + phaseIdx, 1.0);
332 
333  Slast -= S;
334 
335  flashFluidState.setSaturation(phaseIdx, S);
336  }
337  flashFluidState.setSaturation(numPhases - 1, Slast);
338 
339  // copy the pressures: the first pressure is the first primary variable, the
340  // remaining ones are given as p_beta = p_alpha + p_calpha,beta
341  FlashEval p0 = inputFluidState.pressure(0);
342  p0.setDerivative(p0PvIdx, 1.0);
343 
344  std::array<FlashEval, numPhases> pc;
345  MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
346  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
347  flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
348 
349  // copy the mole fractions: all of them are primary variables
350  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
351  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
352  FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
353  x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
354  flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
355  }
356  }
357 
358  flashParamCache.updateAll(flashFluidState);
359 
360  // compute the density of each phase and the fugacity coefficient of each
361  // component in each phase.
362  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
363  const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
364  flashFluidState.setDensity(phaseIdx, rho);
365 
366  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
367  const FlashEval& fugCoeff = FluidSystem::fugacityCoefficient(flashFluidState, flashParamCache, phaseIdx, compIdx);
368  flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
369  }
370  }
371  }
372 
373  template <class FlashFluidState, class OutputFluidState>
374  static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
375  OutputFluidState& outputFluidState)
376  {
377  outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value());
378 
379  // copy the saturations, pressures and densities
380  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
381  const auto& S = flashFluidState.saturation(phaseIdx).value();
382  outputFluidState.setSaturation(phaseIdx, S);
383 
384  const auto& p = flashFluidState.pressure(phaseIdx).value();
385  outputFluidState.setPressure(phaseIdx, p);
386 
387  const auto& rho = flashFluidState.density(phaseIdx).value();
388  outputFluidState.setDensity(phaseIdx, rho);
389  }
390 
391  // copy the mole fractions and fugacity coefficients
392  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
393  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
394  const auto& moleFrac =
395  flashFluidState.moleFraction(phaseIdx, compIdx).value();
396  outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
397 
398  const auto& fugCoeff =
399  flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
400  outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
401  }
402  }
403  }
404 
405  template <class FlashFluidState, class FlashDefectVector, class FlashComponentVector>
406  static void evalDefect_(FlashDefectVector& b,
407  const FlashFluidState& fluidState,
408  const FlashComponentVector& globalMolarities)
409  {
410  typedef typename FlashFluidState::Scalar FlashEval;
411 
412  unsigned eqIdx = 0;
413 
414  // fugacity of any component must be equal in all phases
415  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
416  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
417  b[eqIdx] =
418  fluidState.fugacity(/*phaseIdx=*/0, compIdx)
419  - fluidState.fugacity(phaseIdx, compIdx);
420  ++eqIdx;
421  }
422  }
423  assert(eqIdx == numComponents*(numPhases - 1));
424 
425  // the fact saturations must sum up to 1 is included implicitly and also,
426  // capillary pressures are treated implicitly!
427 
428  // global molarities are given
429  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
430  b[eqIdx] = 0.0;
431  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
432  b[eqIdx] +=
433  fluidState.saturation(phaseIdx)
434  * fluidState.molarity(phaseIdx, compIdx);
435  }
436 
437  b[eqIdx] -= globalMolarities[compIdx];
438  ++eqIdx;
439  }
440 
441  // model assumptions (-> non-linear complementarity functions) must be adhered
442  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
443  FlashEval oneMinusSumMoleFrac = 1.0;
444  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
445  oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
446 
447  if (oneMinusSumMoleFrac > fluidState.saturation(phaseIdx))
448  b[eqIdx] = fluidState.saturation(phaseIdx);
449  else
450  b[eqIdx] = oneMinusSumMoleFrac;
451 
452  ++eqIdx;
453  }
454  }
455 
456  template <class MaterialLaw, class FlashFluidState, class EvalVector>
457  static Scalar update_(FlashFluidState& fluidState,
458  const typename MaterialLaw::Params& matParams,
459  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
460  const EvalVector& deltaX)
461  {
462  // note that it is possible that FlashEval::Scalar is an Evaluation itself
463  typedef typename FlashFluidState::Scalar FlashEval;
464  typedef typename FlashEval::ValueType InnerEval;
465 
466 #ifndef NDEBUG
467  // make sure we don't swallow non-finite update vectors
468  assert(deltaX.dimension == numEq);
469  for (unsigned i = 0; i < numEq; ++i)
470  assert(std::isfinite(Opm::scalarValue(deltaX[i])));
471 #endif
472 
473  Scalar relError = 0;
474  for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
475  FlashEval tmp = getQuantity_(fluidState, pvIdx);
476  InnerEval delta = deltaX[pvIdx];
477 
478  relError = std::max(relError,
479  std::abs(Opm::scalarValue(delta))
480  * quantityWeight_(fluidState, pvIdx));
481 
482  if (isSaturationIdx_(pvIdx)) {
483  // dampen to at most 25% change in saturation per iteration
484  delta = Opm::min(0.25, Opm::max(-0.25, delta));
485  }
486  else if (isMoleFracIdx_(pvIdx)) {
487  // dampen to at most 20% change in mole fraction per iteration
488  delta = Opm::min(0.20, Opm::max(-0.20, delta));
489  }
490  else if (isPressureIdx_(pvIdx)) {
491  // dampen to at most 50% change in pressure per iteration
492  delta = Opm::min(0.5*fluidState.pressure(0).value(),
493  Opm::max(-0.5*fluidState.pressure(0).value(),
494  delta));
495  }
496 
497  tmp -= delta;
498  setQuantity_(fluidState, pvIdx, tmp);
499  }
500 
501  completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
502 
503  return relError;
504  }
505 
506  template <class MaterialLaw, class FlashFluidState>
507  static void completeFluidState_(FlashFluidState& flashFluidState,
508  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
509  const typename MaterialLaw::Params& matParams)
510  {
511  typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
512 
513  typedef typename FlashFluidState::Scalar FlashEval;
514 
515  // calculate the saturation of the last phase as a function of
516  // the other saturations
517  FlashEval sumSat = 0.0;
518  for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
519  sumSat += flashFluidState.saturation(phaseIdx);
520  flashFluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
521 
522  // update the pressures using the material law (saturations
523  // and first pressure are already set because it is implicitly
524  // solved for.)
525  Dune::FieldVector<FlashEval, numPhases> pC;
526  MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
527  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
528  flashFluidState.setPressure(phaseIdx,
529  flashFluidState.pressure(0)
530  + (pC[phaseIdx] - pC[0]));
531 
532  // update the parameter cache
533  paramCache.updateAll(flashFluidState, /*except=*/ParamCache::Temperature);
534 
535  // update all densities and fugacity coefficients
536  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
537  const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
538  flashFluidState.setDensity(phaseIdx, rho);
539 
540  for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
541  const FlashEval& phi = FluidSystem::fugacityCoefficient(flashFluidState, paramCache, phaseIdx, compIdx);
542  flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
543  }
544  }
545  }
546 
547  static bool isPressureIdx_(unsigned pvIdx)
548  { return pvIdx == 0; }
549 
550  static bool isSaturationIdx_(unsigned pvIdx)
551  { return 1 <= pvIdx && pvIdx < numPhases; }
552 
553  static bool isMoleFracIdx_(unsigned pvIdx)
554  { return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
555 
556  // retrieves a quantity from the fluid state
557  template <class FluidState>
558  static const typename FluidState::Scalar& getQuantity_(const FluidState& fluidState, unsigned pvIdx)
559  {
560  assert(pvIdx < numEq);
561 
562  // first pressure
563  if (pvIdx < 1) {
564  unsigned phaseIdx = 0;
565  return fluidState.pressure(phaseIdx);
566  }
567  // first M - 1 saturations
568  else if (pvIdx < numPhases) {
569  unsigned phaseIdx = pvIdx - 1;
570  return fluidState.saturation(phaseIdx);
571  }
572  // mole fractions
573  else // if (pvIdx < numPhases + numPhases*numComponents)
574  {
575  unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
576  unsigned compIdx = (pvIdx - numPhases)%numComponents;
577  return fluidState.moleFraction(phaseIdx, compIdx);
578  }
579  }
580 
581  // set a quantity in the fluid state
582  template <class FluidState>
583  static void setQuantity_(FluidState& fluidState,
584  unsigned pvIdx,
585  const typename FluidState::Scalar& value)
586  {
587  assert(pvIdx < numEq);
588 
589  Valgrind::CheckDefined(value);
590  // first pressure
591  if (pvIdx < 1) {
592  unsigned phaseIdx = 0;
593  fluidState.setPressure(phaseIdx, value);
594  }
595  // first M - 1 saturations
596  else if (pvIdx < numPhases) {
597  unsigned phaseIdx = pvIdx - 1;
598  fluidState.setSaturation(phaseIdx, value);
599  }
600  // mole fractions
601  else {
602  assert(pvIdx < numPhases + numPhases*numComponents);
603  unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
604  unsigned compIdx = (pvIdx - numPhases)%numComponents;
605  fluidState.setMoleFraction(phaseIdx, compIdx, value);
606  }
607  }
608 
609  template <class FluidState>
610  static Scalar quantityWeight_(const FluidState& /*fluidState*/, unsigned pvIdx)
611  {
612  // first pressure
613  if (pvIdx < 1)
614  return 1e-6;
615  // first M - 1 saturations
616  else if (pvIdx < numPhases)
617  return 1.0;
618  // mole fractions
619  else // if (pvIdx < numPhases + numPhases*numComponents)
620  return 1.0;
621  }
622 };
623 
624 } // namespace Opm
625 
626 #endif
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &globalMolarities)
Guess initial values for all quantities.
Definition: NcpFlash.hpp:109
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
This file contains helper classes for the material laws.
static void solve(FluidState &fluidState, const typename MaterialLaw::Params &matParams, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache, const Dune::FieldVector< typename FluidState::Scalar, numComponents > &globalMolarities, Scalar tolerance=-1.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: NcpFlash.hpp:150
A generic traits class which does not provide any indices.
Definition: MaterialTraits.hpp:43
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Definition: NullMaterial.hpp:47
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
static void solve(FluidState &fluidState, const ComponentVector &globalMolarities, Scalar tolerance=0.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: NcpFlash.hpp:250
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:57
Determines the phase compositions, pressures and saturations given the total mass of all components...
Definition: NcpFlash.hpp:91
Representation of an evaluation of a function and its derivatives w.r.t.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: CompositionalFluidState.hpp:46