All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ImmiscibleFlash.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_IMMISCIBLE_FLASH_HPP
28 #define OPM_IMMISCIBLE_FLASH_HPP
29 
34 #include <opm/common/Valgrind.hpp>
35 
36 #include <opm/common/ErrorMacros.hpp>
37 #include <opm/common/Exceptions.hpp>
38 
39 #include <dune/common/fvector.hh>
40 #include <dune/common/fmatrix.hh>
41 
42 #include <limits>
43 #include <iostream>
44 
45 namespace Opm {
46 
73 template <class Scalar, class FluidSystem>
75 {
76  static const int numPhases = FluidSystem::numPhases;
77  static const int numComponents = FluidSystem::numComponents;
78  static_assert(numPhases == numComponents,
79  "Immiscibility assumes that the number of phases"
80  " is equal to the number of components");
81 
82  enum {
83  p0PvIdx = 0,
84  S0PvIdx = 1
85  };
86 
87  static const int numEq = numPhases;
88 
89 public:
93  template <class FluidState, class Evaluation = typename FluidState::Scalar>
94  static void guessInitial(FluidState& fluidState,
95  const Dune::FieldVector<Evaluation, numComponents>& /*globalMolarities*/)
96  {
97  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
98  // pressure. use 1 bar as initial guess
99  fluidState.setPressure(phaseIdx, 1e5);
100 
101  // saturation. assume all fluids to be equally distributed
102  fluidState.setSaturation(phaseIdx, 1.0/numPhases);
103  }
104  }
105 
112  template <class MaterialLaw, class FluidState>
113  static void solve(FluidState& fluidState,
114  const typename MaterialLaw::Params& matParams,
115  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
116  const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
117  Scalar tolerance = -1)
118  {
119  typedef typename FluidState::Scalar InputEval;
120 
122  // Check if all fluid phases are incompressible
124  bool allIncompressible = true;
125  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
126  if (FluidSystem::isCompressible(phaseIdx)) {
127  allIncompressible = false;
128  break;
129  }
130  }
131 
132  if (allIncompressible) {
133  // yes, all fluid phases are incompressible. in this case the flash solver
134  // can only determine the saturations, not the pressures. (but this
135  // determination is much simpler than a full flash calculation.)
136  paramCache.updateAll(fluidState);
137  solveAllIncompressible_(fluidState, paramCache, globalMolarities);
138  return;
139  }
140 
141  typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
142  typedef Dune::FieldVector<InputEval, numEq> Vector;
143 
145  typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
147 
148  Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
149 
150  if (tolerance <= 0)
151  tolerance = std::min<Scalar>(1e-5,
152  1e8*std::numeric_limits<Scalar>::epsilon());
153 
154  typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
155  flashParamCache.assignPersistentData(paramCache);
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  Valgrind::SetUndefined(J);
169  Valgrind::SetUndefined(deltaX);
170  Valgrind::SetUndefined(b);
171 
172  FlashFluidState flashFluidState;
173  assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
174 
175  // copy the global molarities to a vector of evaluations. Remember that the
176  // global molarities are constants. (but we need to copy them to a vector of
177  // FlashEvals anyway in order to avoid getting into hell's kitchen.)
178  Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
179  for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
180  flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
181 
182  FlashDefectVector defect;
183  const unsigned nMax = 50; // <- maximum number of newton iterations
184  for (unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
185  // calculate Jacobian matrix and right hand side
186  evalDefect_(defect, flashFluidState, flashGlobalMolarities);
187  Valgrind::CheckDefined(defect);
188 
189  // create field matrices and vectors out of the evaluation vector to solve
190  // the linear system of equations.
191  for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
192  for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
193  J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
194 
195  b[eqIdx] = defect[eqIdx].value();
196  }
197  Valgrind::CheckDefined(J);
198  Valgrind::CheckDefined(b);
199 
200  // Solve J*x = b
201  deltaX = 0;
202 
203  try { J.solve(deltaX, b); }
204  catch (Dune::FMatrixError e) {
205  throw Opm::NumericalProblem(e.what());
206  }
207  Valgrind::CheckDefined(deltaX);
208 
209  // update the fluid quantities.
210  Scalar relError = update_<MaterialLaw>(flashFluidState, flashParamCache, matParams, deltaX);
211 
212  if (relError < tolerance) {
213  assignOutputFluidState_(flashFluidState, fluidState);
214  return;
215  }
216  }
217 
218  OPM_THROW(Opm::NumericalProblem,
219  "ImmiscibleFlash solver failed: "
220  "{c_alpha^kappa} = {" << globalMolarities << "}, "
221  << "T = " << fluidState.temperature(/*phaseIdx=*/0));
222  }
223 
224 
225 protected:
226  template <class FluidState>
227  static void printFluidState_(const FluidState& fs)
228  {
229  std::cout << "saturations: ";
230  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
231  std::cout << fs.saturation(phaseIdx) << " ";
232  std::cout << "\n";
233 
234  std::cout << "pressures: ";
235  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
236  std::cout << fs.pressure(phaseIdx) << " ";
237  std::cout << "\n";
238 
239  std::cout << "densities: ";
240  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
241  std::cout << fs.density(phaseIdx) << " ";
242  std::cout << "\n";
243 
244  std::cout << "global component molarities: ";
245  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
246  Scalar sum = 0;
247  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
248  sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
249  }
250  std::cout << sum << " ";
251  }
252  std::cout << "\n";
253  }
254 
255  template <class MaterialLaw, class InputFluidState, class FlashFluidState>
256  static void assignFlashFluidState_(const InputFluidState& inputFluidState,
257  FlashFluidState& flashFluidState,
258  const typename MaterialLaw::Params& matParams,
259  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
260  {
261  typedef typename FlashFluidState::Scalar FlashEval;
262 
263  // copy the temperature: even though the model which uses the flash solver might
264  // be non-isothermal, the flash solver does not consider energy. (it could be
265  // modified to do so relatively easily, but it would come at increased
266  // computational cost and normally temperature instead of "total internal energy
267  // of the fluids" is specified.)
268  flashFluidState.setTemperature(inputFluidState.temperature(/*phaseIdx=*/0));
269 
270  // copy the saturations: the first N-1 phases are primary variables, the last one
271  // is one minus the sum of the former.
272  FlashEval Slast = 1.0;
273  for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
274  FlashEval S = inputFluidState.saturation(phaseIdx);
275  S.setDerivative(S0PvIdx + phaseIdx, 1.0);
276 
277  Slast -= S;
278 
279  flashFluidState.setSaturation(phaseIdx, S);
280  }
281  flashFluidState.setSaturation(numPhases - 1, Slast);
282 
283  // copy the pressures: the first pressure is the first primary variable, the
284  // remaining ones are given as p_beta = p_alpha + p_calpha,beta
285  FlashEval p0 = inputFluidState.pressure(0);
286  p0.setDerivative(p0PvIdx, 1.0);
287 
288  std::array<FlashEval, numPhases> pc;
289  MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
290  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
291  flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
292 
293  flashParamCache.updateAll(flashFluidState);
294 
295  // compute the density of each phase.
296  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
297  const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
298  flashFluidState.setDensity(phaseIdx, rho);
299  }
300  }
301 
302  template <class FlashFluidState, class OutputFluidState>
303  static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
304  OutputFluidState& outputFluidState)
305  {
306  outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value());
307 
308  // copy the saturations, pressures and densities
309  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
310  const auto& S = flashFluidState.saturation(phaseIdx).value();
311  outputFluidState.setSaturation(phaseIdx, S);
312 
313  const auto& p = flashFluidState.pressure(phaseIdx).value();
314  outputFluidState.setPressure(phaseIdx, p);
315 
316  const auto& rho = flashFluidState.density(phaseIdx).value();
317  outputFluidState.setDensity(phaseIdx, rho);
318  }
319  }
320 
321  template <class FluidState>
322  static void solveAllIncompressible_(FluidState& fluidState,
323  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
324  const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities)
325  {
326  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
327  Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
328  fluidState.setDensity(phaseIdx, rho);
329 
330  Scalar saturation =
331  globalMolarities[/*compIdx=*/phaseIdx]
332  / fluidState.molarDensity(phaseIdx);
333  fluidState.setSaturation(phaseIdx, saturation);
334  }
335  }
336 
337  template <class FluidState, class FlashDefectVector, class FlashComponentVector>
338  static void evalDefect_(FlashDefectVector& b,
339  const FluidState& fluidState,
340  const FlashComponentVector& globalMolarities)
341  {
342  // global molarities are given
343  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
344  b[phaseIdx] =
345  fluidState.saturation(phaseIdx)
346  * fluidState.molarity(phaseIdx, /*compIdx=*/phaseIdx);
347  b[phaseIdx] -= globalMolarities[/*compIdx=*/phaseIdx];
348  }
349  }
350 
351  template <class MaterialLaw, class FlashFluidState, class EvalVector>
352  static Scalar update_(FlashFluidState& fluidState,
353  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
354  const typename MaterialLaw::Params& matParams,
355  const EvalVector& deltaX)
356  {
357  // note that it is possible that FlashEval::Scalar is an Evaluation itself
358  typedef typename FlashFluidState::Scalar FlashEval;
359  typedef Opm::MathToolbox<FlashEval> FlashEvalToolbox;
360 
361  typedef typename FlashEvalToolbox::ValueType InnerEval;
362 
363 #ifndef NDEBUG
364  // make sure we don't swallow non-finite update vectors
365  assert(deltaX.dimension == numEq);
366  for (unsigned i = 0; i < numEq; ++i)
367  assert(std::isfinite(Opm::scalarValue(deltaX[i])));
368 #endif
369 
370  Scalar relError = 0;
371  for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
372  FlashEval tmp = getQuantity_(fluidState, pvIdx);
373  InnerEval delta = deltaX[pvIdx];
374 
375  relError = std::max(relError,
376  std::abs(Opm::scalarValue(delta))
377  * quantityWeight_(fluidState, pvIdx));
378 
379  if (isSaturationIdx_(pvIdx)) {
380  // dampen to at most 20% change in saturation per
381  // iteration
382  delta = Opm::min(0.25, Opm::max(-0.25, delta));
383  }
384  else if (isPressureIdx_(pvIdx)) {
385  // dampen to at most 30% change in pressure per
386  // iteration
387  delta = Opm::min(0.5*fluidState.pressure(0).value(),
388  Opm::max(-0.50*fluidState.pressure(0).value(), delta));
389  }
390 
391  tmp -= delta;
392  setQuantity_(fluidState, pvIdx, tmp);
393  }
394 
395  completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
396 
397  return relError;
398  }
399 
400  template <class MaterialLaw, class FlashFluidState>
401  static void completeFluidState_(FlashFluidState& flashFluidState,
402  typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
403  const typename MaterialLaw::Params& matParams)
404  {
405  typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
406 
407  typedef typename FlashFluidState::Scalar FlashEval;
408 
409  // calculate the saturation of the last phase as a function of
410  // the other saturations
411  FlashEval sumSat = 0.0;
412  for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
413  sumSat += flashFluidState.saturation(phaseIdx);
414  flashFluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
415 
416  // update the pressures using the material law (saturations
417  // and first pressure are already set because it is implicitly
418  // solved for.)
419  Dune::FieldVector<FlashEval, numPhases> pC;
420  MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
421  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
422  flashFluidState.setPressure(phaseIdx,
423  flashFluidState.pressure(0)
424  + (pC[phaseIdx] - pC[0]));
425 
426  // update the parameter cache
427  paramCache.updateAll(flashFluidState, /*except=*/ParamCache::Temperature|ParamCache::Composition);
428 
429  // update all densities
430  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
431  const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
432  flashFluidState.setDensity(phaseIdx, rho);
433  }
434  }
435 
436  static bool isPressureIdx_(unsigned pvIdx)
437  { return pvIdx == 0; }
438 
439  static bool isSaturationIdx_(unsigned pvIdx)
440  { return 1 <= pvIdx && pvIdx < numPhases; }
441 
442  // retrieves a quantity from the fluid state
443  template <class FluidState>
444  static const typename FluidState::Scalar& getQuantity_(const FluidState& fluidState, unsigned pvIdx)
445  {
446  assert(pvIdx < numEq);
447 
448  // first pressure
449  if (pvIdx < 1) {
450  unsigned phaseIdx = 0;
451  return fluidState.pressure(phaseIdx);
452  }
453  // saturations
454  else {
455  assert(pvIdx < numPhases);
456  unsigned phaseIdx = pvIdx - 1;
457  return fluidState.saturation(phaseIdx);
458  }
459  }
460 
461  // set a quantity in the fluid state
462  template <class FluidState>
463  static void setQuantity_(FluidState& fluidState,
464  unsigned pvIdx,
465  const typename FluidState::Scalar& value)
466  {
467  assert(pvIdx < numEq);
468 
469  // first pressure
470  if (pvIdx < 1) {
471  unsigned phaseIdx = 0;
472  fluidState.setPressure(phaseIdx, value);
473  }
474  // saturations
475  else {
476  assert(pvIdx < numPhases);
477  unsigned phaseIdx = pvIdx - 1;
478  fluidState.setSaturation(phaseIdx, value);
479  }
480  }
481 
482  template <class FluidState>
483  static Scalar quantityWeight_(const FluidState& /*fluidState*/, unsigned pvIdx)
484  {
485  // first pressure
486  if (pvIdx < 1)
487  return 1e-8;
488  // first M - 1 saturations
489  else {
490  assert(pvIdx < numPhases);
491  return 1.0;
492  }
493  }
494 };
495 
496 } // namespace Opm
497 
498 #endif
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: ImmiscibleFluidState.hpp:49
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &)
Guess initial values for all quantities.
Definition: ImmiscibleFlash.hpp:94
Definition: MathToolbox.hpp:48
Determines the pressures and saturations of all fluid phases given the total mass of all components...
Definition: ImmiscibleFlash.hpp:74
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:57
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)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: ImmiscibleFlash.hpp:113
Representation of an evaluation of a function and its derivatives w.r.t.