My Project
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 
37 
39 
40 #include <dune/common/fvector.hh>
41 #include <dune/common/fmatrix.hh>
42 #include <dune/common/version.hh>
43 
44 #include <limits>
45 #include <iostream>
46 
47 namespace Opm {
48 
88 template <class Scalar, class FluidSystem>
89 class NcpFlash
90 {
91  enum { numPhases = FluidSystem::numPhases };
92  enum { numComponents = FluidSystem::numComponents };
93 
94  enum {
95  p0PvIdx = 0,
96  S0PvIdx = 1,
97  x00PvIdx = S0PvIdx + numPhases - 1
98  };
99 
100  static const int numEq = numPhases*(numComponents + 1);
101 
102 public:
106  template <class FluidState, class Evaluation = typename FluidState::Scalar>
107  static void guessInitial(FluidState& fluidState,
108  const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
109  {
110  // the sum of all molarities
111  Evaluation sumMoles = 0;
112  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
113  sumMoles += globalMolarities[compIdx];
114 
115  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
116  // composition
117  for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
118  fluidState.setMoleFraction(phaseIdx,
119  compIdx,
120  globalMolarities[compIdx]/sumMoles);
121 
122  // pressure. use atmospheric pressure as initial guess
123  fluidState.setPressure(phaseIdx, 1.0135e5);
124 
125  // saturation. assume all fluids to be equally distributed
126  fluidState.setSaturation(phaseIdx, 1.0/numPhases);
127  }
128 
129  // set the fugacity coefficients of all components in all phases
130  typename FluidSystem::template ParameterCache<Evaluation> paramCache;
131  paramCache.updateAll(fluidState);
132  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
133  for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
134  const typename FluidState::Scalar phi =
135  FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
136  fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
137  }
138  }
139  }
140 
147  template <class MaterialLaw, class FluidState>
148  static void solve(FluidState& fluidState,
149  const typename MaterialLaw::Params& matParams,
150  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
151  const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
152  Scalar tolerance = -1.0)
153  {
154  typedef typename FluidState::Scalar InputEval;
155 
156  typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
157  typedef Dune::FieldVector<InputEval, numEq> Vector;
158 
159  typedef DenseAd::Evaluation</*Scalar=*/InputEval,
160  /*numDerivs=*/numEq> FlashEval;
161 
162  typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
163  typedef CompositionalFluidState<FlashEval, FluidSystem, /*energy=*/false> FlashFluidState;
164 
165 #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
166  Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
167 #endif
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 NumericalIssue(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  std::ostringstream oss;
237  oss << "NcpFlash solver failed:"
238  << " {c_alpha^kappa} = {" << globalMolarities << "}, "
239  << " T = " << fluidState.temperature(/*phaseIdx=*/0);
240  throw NumericalIssue(oss.str());
241  }
242 
250  template <class FluidState, class ComponentVector>
251  static void solve(FluidState& fluidState,
252  const ComponentVector& globalMolarities,
253  Scalar tolerance = 0.0)
254  {
255  typedef NullMaterialTraits<Scalar, numPhases> MaterialTraits;
256  typedef NullMaterial<MaterialTraits> MaterialLaw;
257  typedef typename MaterialLaw::Params MaterialLawParams;
258 
259  MaterialLawParams matParams;
260  solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
261  }
262 
263 
264 protected:
265  template <class FluidState>
266  static void printFluidState_(const FluidState& fluidState)
267  {
268  typedef typename FluidState::Scalar FsScalar;
269 
270  std::cout << "saturations: ";
271  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
272  std::cout << fluidState.saturation(phaseIdx) << " ";
273  std::cout << "\n";
274 
275  std::cout << "pressures: ";
276  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
277  std::cout << fluidState.pressure(phaseIdx) << " ";
278  std::cout << "\n";
279 
280  std::cout << "densities: ";
281  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
282  std::cout << fluidState.density(phaseIdx) << " ";
283  std::cout << "\n";
284 
285  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
286  std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
287  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
288  std::cout << fluidState.moleFraction(phaseIdx, compIdx) << " ";
289  }
290  std::cout << "\n";
291  }
292 
293  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
294  std::cout << "fugacities " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
295  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
296  std::cout << fluidState.fugacity(phaseIdx, compIdx) << " ";
297  }
298  std::cout << "\n";
299  }
300 
301  std::cout << "global component molarities: ";
302  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
303  FsScalar sum = 0;
304  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
305  sum += fluidState.saturation(phaseIdx)*fluidState.molarity(phaseIdx, compIdx);
306  }
307  std::cout << sum << " ";
308  }
309  std::cout << "\n";
310  }
311 
312  template <class MaterialLaw, class InputFluidState, class FlashFluidState>
313  static void assignFlashFluidState_(const InputFluidState& inputFluidState,
314  FlashFluidState& flashFluidState,
315  const typename MaterialLaw::Params& matParams,
316  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
317  {
318  typedef typename FlashFluidState::Scalar FlashEval;
319 
320  // copy the temperature: even though the model which uses the flash solver might
321  // be non-isothermal, the flash solver does not consider energy. (it could be
322  // modified to do so relatively easily, but it would come at increased
323  // computational cost and normally temperature instead of "total internal energy
324  // of the fluids" is specified.)
325  flashFluidState.setTemperature(inputFluidState.temperature(/*phaseIdx=*/0));
326 
327  // copy the saturations: the first N-1 phases are primary variables, the last one
328  // is one minus the sum of the former.
329  FlashEval Slast = 1.0;
330  for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
331  FlashEval S = inputFluidState.saturation(phaseIdx);
332  S.setDerivative(S0PvIdx + phaseIdx, 1.0);
333 
334  Slast -= S;
335 
336  flashFluidState.setSaturation(phaseIdx, S);
337  }
338  flashFluidState.setSaturation(numPhases - 1, Slast);
339 
340  // copy the pressures: the first pressure is the first primary variable, the
341  // remaining ones are given as p_beta = p_alpha + p_calpha,beta
342  FlashEval p0 = inputFluidState.pressure(0);
343  p0.setDerivative(p0PvIdx, 1.0);
344 
345  std::array<FlashEval, numPhases> pc;
346  MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
347  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
348  flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
349 
350  // copy the mole fractions: all of them are primary variables
351  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
352  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
353  FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
354  x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
355  flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
356  }
357  }
358 
359  flashParamCache.updateAll(flashFluidState);
360 
361  // compute the density of each phase and the fugacity coefficient of each
362  // component in each phase.
363  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
364  const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
365  flashFluidState.setDensity(phaseIdx, rho);
366 
367  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
368  const FlashEval& fugCoeff = FluidSystem::fugacityCoefficient(flashFluidState, flashParamCache, phaseIdx, compIdx);
369  flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
370  }
371  }
372  }
373 
374  template <class FlashFluidState, class OutputFluidState>
375  static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
376  OutputFluidState& outputFluidState)
377  {
378  outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value());
379 
380  // copy the saturations, pressures and densities
381  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
382  const auto& S = flashFluidState.saturation(phaseIdx).value();
383  outputFluidState.setSaturation(phaseIdx, S);
384 
385  const auto& p = flashFluidState.pressure(phaseIdx).value();
386  outputFluidState.setPressure(phaseIdx, p);
387 
388  const auto& rho = flashFluidState.density(phaseIdx).value();
389  outputFluidState.setDensity(phaseIdx, rho);
390  }
391 
392  // copy the mole fractions and fugacity coefficients
393  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
394  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
395  const auto& moleFrac =
396  flashFluidState.moleFraction(phaseIdx, compIdx).value();
397  outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
398 
399  const auto& fugCoeff =
400  flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
401  outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
402  }
403  }
404  }
405 
406  template <class FlashFluidState, class FlashDefectVector, class FlashComponentVector>
407  static void evalDefect_(FlashDefectVector& b,
408  const FlashFluidState& fluidState,
409  const FlashComponentVector& globalMolarities)
410  {
411  typedef typename FlashFluidState::Scalar FlashEval;
412 
413  unsigned eqIdx = 0;
414 
415  // fugacity of any component must be equal in all phases
416  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
417  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
418  b[eqIdx] =
419  fluidState.fugacity(/*phaseIdx=*/0, compIdx)
420  - fluidState.fugacity(phaseIdx, compIdx);
421  ++eqIdx;
422  }
423  }
424  assert(eqIdx == numComponents*(numPhases - 1));
425 
426  // the fact saturations must sum up to 1 is included implicitly and also,
427  // capillary pressures are treated implicitly!
428 
429  // global molarities are given
430  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
431  b[eqIdx] = 0.0;
432  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
433  b[eqIdx] +=
434  fluidState.saturation(phaseIdx)
435  * fluidState.molarity(phaseIdx, compIdx);
436  }
437 
438  b[eqIdx] -= globalMolarities[compIdx];
439  ++eqIdx;
440  }
441 
442  // model assumptions (-> non-linear complementarity functions) must be adhered
443  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
444  FlashEval oneMinusSumMoleFrac = 1.0;
445  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
446  oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
447 
448  if (oneMinusSumMoleFrac > fluidState.saturation(phaseIdx))
449  b[eqIdx] = fluidState.saturation(phaseIdx);
450  else
451  b[eqIdx] = oneMinusSumMoleFrac;
452 
453  ++eqIdx;
454  }
455  }
456 
457  template <class MaterialLaw, class FlashFluidState, class EvalVector>
458  static Scalar update_(FlashFluidState& fluidState,
459  const typename MaterialLaw::Params& matParams,
460  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
461  const EvalVector& deltaX)
462  {
463  // note that it is possible that FlashEval::Scalar is an Evaluation itself
464  typedef typename FlashFluidState::Scalar FlashEval;
465  typedef typename FlashEval::ValueType InnerEval;
466 
467 #ifndef NDEBUG
468  // make sure we don't swallow non-finite update vectors
469  assert(deltaX.dimension == numEq);
470  for (unsigned i = 0; i < numEq; ++i)
471  assert(std::isfinite(scalarValue(deltaX[i])));
472 #endif
473 
474  Scalar relError = 0;
475  for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
476  FlashEval tmp = getQuantity_(fluidState, pvIdx);
477  InnerEval delta = deltaX[pvIdx];
478 
479  relError = std::max(relError,
480  std::abs(scalarValue(delta))
481  * quantityWeight_(fluidState, pvIdx));
482 
483  if (isSaturationIdx_(pvIdx)) {
484  // dampen to at most 25% change in saturation per iteration
485  delta = min(0.25, max(-0.25, delta));
486  }
487  else if (isMoleFracIdx_(pvIdx)) {
488  // dampen to at most 20% change in mole fraction per iteration
489  delta = min(0.20, max(-0.20, delta));
490  }
491  else if (isPressureIdx_(pvIdx)) {
492  // dampen to at most 50% change in pressure per iteration
493  delta = min(0.5*fluidState.pressure(0).value(),
494  max(-0.5*fluidState.pressure(0).value(),
495  delta));
496  }
497 
498  tmp -= delta;
499  setQuantity_(fluidState, pvIdx, tmp);
500  }
501 
502  completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
503 
504  return relError;
505  }
506 
507  template <class MaterialLaw, class FlashFluidState>
508  static void completeFluidState_(FlashFluidState& flashFluidState,
509  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
510  const typename MaterialLaw::Params& matParams)
511  {
512  typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
513 
514  typedef typename FlashFluidState::Scalar FlashEval;
515 
516  // calculate the saturation of the last phase as a function of
517  // the other saturations
518  FlashEval sumSat = 0.0;
519  for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
520  sumSat += flashFluidState.saturation(phaseIdx);
521  flashFluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
522 
523  // update the pressures using the material law (saturations
524  // and first pressure are already set because it is implicitly
525  // solved for.)
526  std::array<FlashEval, numPhases> pC;
527  MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
528  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
529  flashFluidState.setPressure(phaseIdx,
530  flashFluidState.pressure(0)
531  + (pC[phaseIdx] - pC[0]));
532 
533  // update the parameter cache
534  paramCache.updateAll(flashFluidState, /*except=*/ParamCache::Temperature);
535 
536  // update all densities and fugacity coefficients
537  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
538  const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
539  flashFluidState.setDensity(phaseIdx, rho);
540 
541  for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
542  const FlashEval& phi = FluidSystem::fugacityCoefficient(flashFluidState, paramCache, phaseIdx, compIdx);
543  flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
544  }
545  }
546  }
547 
548  static bool isPressureIdx_(unsigned pvIdx)
549  { return pvIdx == 0; }
550 
551  static bool isSaturationIdx_(unsigned pvIdx)
552  { return 1 <= pvIdx && pvIdx < numPhases; }
553 
554  static bool isMoleFracIdx_(unsigned pvIdx)
555  { return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
556 
557  // retrieves a quantity from the fluid state
558  template <class FluidState>
559  static const typename FluidState::Scalar& getQuantity_(const FluidState& fluidState, unsigned pvIdx)
560  {
561  assert(pvIdx < numEq);
562 
563  // first pressure
564  if (pvIdx < 1) {
565  unsigned phaseIdx = 0;
566  return fluidState.pressure(phaseIdx);
567  }
568  // first M - 1 saturations
569  else if (pvIdx < numPhases) {
570  unsigned phaseIdx = pvIdx - 1;
571  return fluidState.saturation(phaseIdx);
572  }
573  // mole fractions
574  else // if (pvIdx < numPhases + numPhases*numComponents)
575  {
576  unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
577  unsigned compIdx = (pvIdx - numPhases)%numComponents;
578  return fluidState.moleFraction(phaseIdx, compIdx);
579  }
580  }
581 
582  // set a quantity in the fluid state
583  template <class FluidState>
584  static void setQuantity_(FluidState& fluidState,
585  unsigned pvIdx,
586  const typename FluidState::Scalar& value)
587  {
588  assert(pvIdx < numEq);
589 
590  Valgrind::CheckDefined(value);
591  // first pressure
592  if (pvIdx < 1) {
593  unsigned phaseIdx = 0;
594  fluidState.setPressure(phaseIdx, value);
595  }
596  // first M - 1 saturations
597  else if (pvIdx < numPhases) {
598  unsigned phaseIdx = pvIdx - 1;
599  fluidState.setSaturation(phaseIdx, value);
600  }
601  // mole fractions
602  else {
603  assert(pvIdx < numPhases + numPhases*numComponents);
604  unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
605  unsigned compIdx = (pvIdx - numPhases)%numComponents;
606  fluidState.setMoleFraction(phaseIdx, compIdx, value);
607  }
608  }
609 
610  template <class FluidState>
611  static Scalar quantityWeight_(const FluidState& /*fluidState*/, unsigned pvIdx)
612  {
613  // first pressure
614  if (pvIdx < 1)
615  return 1e-6;
616  // first M - 1 saturations
617  else if (pvIdx < numPhases)
618  return 1.0;
619  // mole fractions
620  else // if (pvIdx < numPhases + numPhases*numComponents)
621  return 1.0;
622  }
623 };
624 
625 } // namespace Opm
626 
627 #endif
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Representation of an evaluation of a function and its derivatives w.r.t.
Provides the opm-material specific exception classes.
This file contains helper classes for the material laws.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Some templates to wrap the valgrind client request macros.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: CompositionalFluidState.hpp:46
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:59
Determines the phase compositions, pressures and saturations given the total mass of all components.
Definition: NcpFlash.hpp:90
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &globalMolarities)
Guess initial values for all quantities.
Definition: NcpFlash.hpp:107
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:148
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:251
A generic traits class which does not provide any indices.
Definition: MaterialTraits.hpp:45
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Definition: NullMaterial.hpp:46
Definition: Exceptions.hpp:46