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