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