MiscibleMultiPhaseComposition.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_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
28 #define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
29 
31 
32 
33 #include <opm/common/utility/platform_dependent/disable_warnings.h>
34 
35 #include <dune/common/fvector.hh>
36 #include <dune/common/fmatrix.hh>
37 
38 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
39 
40 
41 #include <opm/common/Exceptions.hpp>
42 #include <opm/common/ErrorMacros.hpp>
43 #include <opm/common/Valgrind.hpp>
44 
45 namespace Opm {
46 
54 template <class Scalar>
56 {
57 public:
59  {}
60 
61  MMPCAuxConstraint(unsigned phaseIndex, unsigned compIndex, Scalar val)
62  : phaseIdx_(phaseIndex)
63  , compIdx_(compIndex)
64  , value_(val)
65  {}
66 
70  void set(unsigned phaseIndex, unsigned compIndex, Scalar val)
71  {
72  phaseIdx_ = phaseIndex;
73  compIdx_ = compIndex;
74  value_ = val;
75  }
76 
81  unsigned phaseIdx() const
82  { return phaseIdx_; }
83 
88  unsigned compIdx() const
89  { return compIdx_; }
90 
95  Scalar value() const
96  { return value_; }
97 
98 private:
99  unsigned phaseIdx_;
100  unsigned compIdx_;
101  Scalar value_;
102 };
103 
127 template <class Scalar, class FluidSystem, class Evaluation = Scalar>
129 {
130  static const int numPhases = FluidSystem::numPhases;
131  static const int numComponents = FluidSystem::numComponents;
132 
134 
135  static_assert(numPhases <= numComponents,
136  "This solver requires that the number fluid phases is smaller or equal "
137  "to the number of components");
138 
139 
140 public:
164  template <class FluidState, class ParameterCache>
165  static void solve(FluidState& fluidState,
166  ParameterCache& paramCache,
167  int phasePresence,
168  const MMPCAuxConstraint<Evaluation>* auxConstraints,
169  unsigned numAuxConstraints,
170  bool setViscosity,
171  bool setInternalEnergy)
172  {
173  static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
174  "The scalar type of the fluid state must be 'Evaluation'");
175 
176 #ifndef NDEBUG
177  // currently this solver can only handle fluid systems which
178  // assume ideal mixtures of all fluids. TODO: relax this
179  // (requires solving a non-linear system of equations, i.e. using
180  // newton method.)
181  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
182  assert(FluidSystem::isIdealMixture(phaseIdx));
183  }
184 #endif
185 
186  // compute all fugacity coefficients
187  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
188  paramCache.updatePhase(fluidState, phaseIdx);
189 
190  // since we assume ideal mixtures, the fugacity
191  // coefficients of the components cannot depend on
192  // composition, i.e. the parameters in the cache are valid
193  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
194  Evaluation fugCoeff = Opm::decay<Evaluation>(
195  FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
196  fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
197  }
198  }
199 
200  // create the linear system of equations which defines the
201  // mole fractions
202  static const int numEq = numComponents*numPhases;
203  Dune::FieldMatrix<Evaluation, numEq, numEq> M(0.0);
204  Dune::FieldVector<Evaluation, numEq> x(0.0);
205  Dune::FieldVector<Evaluation, numEq> b(0.0);
206 
207  // assemble the equations expressing the fact that the
208  // fugacities of each component are equal in all phases
209  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
210  const Evaluation& entryCol1 =
211  fluidState.fugacityCoefficient(/*phaseIdx=*/0, compIdx)
212  *fluidState.pressure(/*phaseIdx=*/0);
213  unsigned col1Idx = compIdx;
214 
215  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
216  unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
217  unsigned col2Idx = phaseIdx*numComponents + compIdx;
218 
219  const Evaluation& entryCol2 =
220  fluidState.fugacityCoefficient(phaseIdx, compIdx)
221  *fluidState.pressure(phaseIdx);
222 
223  M[rowIdx][col1Idx] = entryCol1;
224  M[rowIdx][col2Idx] = -entryCol2;
225  }
226  }
227 
228  // assemble the equations expressing the assumption that the
229  // sum of all mole fractions in each phase must be 1 for the
230  // phases present.
231  unsigned presentPhases = 0;
232  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
233  if (!(phasePresence& (1 << phaseIdx)))
234  continue;
235 
236  unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
237  presentPhases += 1;
238 
239  b[rowIdx] = 1.0;
240  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
241  unsigned colIdx = phaseIdx*numComponents + compIdx;
242 
243  M[rowIdx][colIdx] = 1.0;
244  }
245  }
246 
247  assert(presentPhases + numAuxConstraints == numComponents);
248 
249  // incorperate the auxiliary equations, i.e., the explicitly given mole fractions
250  for (unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
251  unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
252  b[rowIdx] = auxConstraints[auxEqIdx].value();
253 
254  unsigned colIdx = auxConstraints[auxEqIdx].phaseIdx()*numComponents + auxConstraints[auxEqIdx].compIdx();
255  M[rowIdx][colIdx] = 1.0;
256  }
257 
258  // solve for all mole fractions
259  try {
260  Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-50);
261  M.solve(x, b);
262  }
263  catch (const Dune::FMatrixError& e) {
264  OPM_THROW(NumericalProblem,
265  "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() << "; M="<<M);
266  }
267  catch (...) {
268  throw;
269  }
270 
271 
272  // set all mole fractions and the additional quantities in
273  // the fluid state
274  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
275  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
276  unsigned rowIdx = phaseIdx*numComponents + compIdx;
277  fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
278  }
279  paramCache.updateComposition(fluidState, phaseIdx);
280 
281  const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
282  fluidState.setDensity(phaseIdx, rho);
283 
284  if (setViscosity) {
285  const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
286  fluidState.setViscosity(phaseIdx, mu);
287  }
288 
289  if (setInternalEnergy) {
290  const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
291  fluidState.setEnthalpy(phaseIdx, h);
292  }
293  }
294  }
295 
303  template <class FluidState, class ParameterCache>
304  static void solve(FluidState& fluidState,
305  ParameterCache& paramCache,
306  bool setViscosity,
307  bool setInternalEnergy)
308  {
309  solve(fluidState,
310  paramCache,
311  /*phasePresence=*/0xffffff,
312  /*numAuxConstraints=*/0,
313  /*auxConstraints=*/0,
314  setViscosity,
315  setInternalEnergy);
316  }
317 };
318 
319 } // namespace Opm
320 
321 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:128
Scalar value() const
Returns value of the mole fraction of the auxiliary constraint.
Definition: MiscibleMultiPhaseComposition.hpp:95
unsigned phaseIdx() const
Returns the index of the fluid phase for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:81
Definition: Air_Mesitylene.hpp:33
Definition: MathToolbox.hpp:48
static void solve(FluidState &fluidState, ParameterCache &paramCache, int phasePresence, const MMPCAuxConstraint< Evaluation > *auxConstraints, unsigned numAuxConstraints, bool setViscosity, bool setInternalEnergy)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:165
unsigned compIdx() const
Returns the index of the component for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:88
static void solve(FluidState &fluidState, ParameterCache &paramCache, bool setViscosity, bool setInternalEnergy)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:304
Specifies an auxiliary constraint for the MiscibleMultiPhaseComposition constraint solver...
Definition: MiscibleMultiPhaseComposition.hpp:55