All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
CompositionFromFugacities.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_COMPOSITION_FROM_FUGACITIES_HPP
28 #define OPM_COMPOSITION_FROM_FUGACITIES_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/ErrorMacros.hpp>
42 #include <opm/common/Exceptions.hpp>
43 #include <opm/common/Valgrind.hpp>
44 
45 #include <limits>
46 
47 namespace Opm {
48 
53 template <class Scalar, class FluidSystem, class Evaluation = Scalar>
55 {
56  enum { numComponents = FluidSystem::numComponents };
57 
58 public:
59  typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
60 
64  template <class FluidState>
65  static void guessInitial(FluidState& fluidState,
66  unsigned phaseIdx,
67  const ComponentVector& /*fugVec*/)
68  {
69  if (FluidSystem::isIdealMixture(phaseIdx))
70  return;
71 
72  // Pure component fugacities
73  for (unsigned i = 0; i < numComponents; ++ i) {
74  //std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
75  fluidState.setMoleFraction(phaseIdx,
76  i,
77  1.0/numComponents);
78  }
79  }
80 
87  template <class FluidState>
88  static void solve(FluidState& fluidState,
89  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
90  unsigned phaseIdx,
91  const ComponentVector& targetFug)
92  {
93  // use a much more efficient method in case the phase is an
94  // ideal mixture
95  if (FluidSystem::isIdealMixture(phaseIdx)) {
96  solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
97  return;
98  }
99 
100  //Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-25);
101 
102  // save initial composition in case something goes wrong
103  Dune::FieldVector<Evaluation, numComponents> xInit;
104  for (unsigned i = 0; i < numComponents; ++i) {
105  xInit[i] = fluidState.moleFraction(phaseIdx, i);
106  }
107 
109  // Newton method
111 
112  // Jacobian matrix
113  Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
114  // solution, i.e. phase composition
115  Dune::FieldVector<Evaluation, numComponents> x;
116  // right hand side
117  Dune::FieldVector<Evaluation, numComponents> b;
118 
119  paramCache.updatePhase(fluidState, phaseIdx);
120 
121  // maximum number of iterations
122  const int nMax = 25;
123  for (int nIdx = 0; nIdx < nMax; ++nIdx) {
124  // calculate Jacobian matrix and right hand side
125  linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
126  Valgrind::CheckDefined(J);
127  Valgrind::CheckDefined(b);
128 
129  /*
130  std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
131  for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
132  std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
133  std::cout << "\n";
134  std::cout << FluidSystem::phaseName(phaseIdx) << "Phase phi: ";
135  for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
136  std::cout << fluidState.fugacityCoefficient(phaseIdx, i) << " ";
137  std::cout << "\n";
138  */
139 
140  // Solve J*x = b
141  x = 0.0;
142  try { J.solve(x, b); }
143  catch (Dune::FMatrixError e)
144  { throw Opm::NumericalProblem(e.what()); }
145 
146  //std::cout << "original delta: " << x << "\n";
147 
148  Valgrind::CheckDefined(x);
149 
150  /*
151  std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
152  for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
153  std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
154  std::cout << "\n";
155  std::cout << "J: " << J << "\n";
156  std::cout << "rho: " << fluidState.density(phaseIdx) << "\n";
157  std::cout << "delta: " << x << "\n";
158  std::cout << "defect: " << b << "\n";
159 
160  std::cout << "J: " << J << "\n";
161 
162  std::cout << "---------------------------\n";
163  */
164 
165  // update the fluid composition. b is also used to store
166  // the defect for the next iteration.
167  Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
168 
169  if (relError < 1e-9) {
170  const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
171  fluidState.setDensity(phaseIdx, rho);
172 
173  //std::cout << "num iterations: " << nIdx << "\n";
174  return;
175  }
176  }
177 
178  OPM_THROW(Opm::NumericalProblem,
179  "Calculating the " << FluidSystem::phaseName(phaseIdx)
180  << "Phase composition failed. Initial {x} = {"
181  << xInit
182  << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx)
183  << ", T = " << fluidState.temperature(phaseIdx));
184  }
185 
186 
187 protected:
188  // update the phase composition in case the phase is an ideal
189  // mixture, i.e. the component's fugacity coefficients are
190  // independent of the phase's composition.
191  template <class FluidState>
192  static void solveIdealMix_(FluidState& fluidState,
193  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
194  unsigned phaseIdx,
195  const ComponentVector& fugacities)
196  {
197  for (unsigned i = 0; i < numComponents; ++ i) {
198  const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
199  paramCache,
200  phaseIdx,
201  i);
202  const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
203  Valgrind::CheckDefined(phi);
204  Valgrind::CheckDefined(gamma);
205  Valgrind::CheckDefined(fugacities[i]);
206  fluidState.setFugacityCoefficient(phaseIdx, i, phi);
207  fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
208  };
209 
210  paramCache.updatePhase(fluidState, phaseIdx);
211 
212  const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
213  fluidState.setDensity(phaseIdx, rho);
214  return;
215  }
216 
217  template <class FluidState>
218  static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents>& J,
219  Dune::FieldVector<Evaluation, numComponents>& defect,
220  FluidState& fluidState,
221  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
222  unsigned phaseIdx,
223  const ComponentVector& targetFug)
224  {
225  // reset jacobian
226  J = 0;
227 
228  Scalar absError = 0;
229  // calculate the defect (deviation of the current fugacities
230  // from the target fugacities)
231  for (unsigned i = 0; i < numComponents; ++ i) {
232  const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
233  paramCache,
234  phaseIdx,
235  i);
236  const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
237  fluidState.setFugacityCoefficient(phaseIdx, i, phi);
238 
239  defect[i] = targetFug[i] - f;
240  absError = std::max(absError, std::abs(Opm::scalarValue(defect[i])));
241  }
242 
243  // assemble jacobian matrix of the constraints for the composition
244  static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
245  for (unsigned i = 0; i < numComponents; ++ i) {
247  // approximately calculate partial derivatives of the
248  // fugacity defect of all components in regard to the mole
249  // fraction of the i-th component. This is done via
250  // forward differences
251 
252  // deviate the mole fraction of the i-th component
253  Evaluation xI = fluidState.moleFraction(phaseIdx, i);
254  fluidState.setMoleFraction(phaseIdx, i, xI + eps);
255  paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
256 
257  // compute new defect and derivative for all component
258  // fugacities
259  for (unsigned j = 0; j < numComponents; ++j) {
260  // compute the j-th component's fugacity coefficient ...
261  const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
262  paramCache,
263  phaseIdx,
264  j);
265  // ... and its fugacity ...
266  const Evaluation& f =
267  phi *
268  fluidState.pressure(phaseIdx) *
269  fluidState.moleFraction(phaseIdx, j);
270  // as well as the defect for this fugacity
271  const Evaluation& defJPlusEps = targetFug[j] - f;
272 
273  // use forward differences to calculate the defect's
274  // derivative
275  J[j][i] = (defJPlusEps - defect[j])/eps;
276  }
277 
278  // reset composition to original value
279  fluidState.setMoleFraction(phaseIdx, i, xI);
280  paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
281 
282  // end forward differences
284  }
285 
286  return absError;
287  }
288 
289  template <class FluidState>
290  static Scalar update_(FluidState& fluidState,
291  typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
292  Dune::FieldVector<Evaluation, numComponents>& x,
293  Dune::FieldVector<Evaluation, numComponents>& /*b*/,
294  unsigned phaseIdx,
295  const Dune::FieldVector<Evaluation, numComponents>& targetFug)
296  {
297  // store original composition and calculate relative error
298  Dune::FieldVector<Evaluation, numComponents> origComp;
299  Scalar relError = 0;
300  Evaluation sumDelta = 0.0;
301  Evaluation sumx = 0.0;
302  for (unsigned i = 0; i < numComponents; ++i) {
303  origComp[i] = fluidState.moleFraction(phaseIdx, i);
304  relError = std::max(relError, std::abs(Opm::scalarValue(x[i])));
305 
306  sumx += Opm::abs(fluidState.moleFraction(phaseIdx, i));
307  sumDelta += Opm::abs(x[i]);
308  }
309 
310  // chop update to at most 20% change in composition
311  const Scalar maxDelta = 0.2;
312  if (sumDelta > maxDelta)
313  x /= (sumDelta/maxDelta);
314 
315  // change composition
316  for (unsigned i = 0; i < numComponents; ++i) {
317  Evaluation newx = origComp[i] - x[i];
318  // only allow negative mole fractions if the target fugacity is negative
319  if (targetFug[i] > 0)
320  newx = Opm::max(0.0, newx);
321  // only allow positive mole fractions if the target fugacity is positive
322  else if (targetFug[i] < 0)
323  newx = Opm::min(0.0, newx);
324  // if the target fugacity is zero, the mole fraction must also be zero
325  else
326  newx = 0;
327 
328  fluidState.setMoleFraction(phaseIdx, i, newx);
329  }
330 
331  paramCache.updateComposition(fluidState, phaseIdx);
332 
333  return relError;
334  }
335 
336  template <class FluidState>
337  static Scalar calculateDefect_(const FluidState& params,
338  unsigned phaseIdx,
339  const ComponentVector& targetFug)
340  {
341  Scalar result = 0.0;
342  for (unsigned i = 0; i < numComponents; ++i) {
343  // sum of the fugacity defect weighted by the inverse
344  // fugacity coefficient
345  result += std::abs(
346  (targetFug[i] - params.fugacity(phaseIdx, i))
347  /
348  params.fugacityCoefficient(phaseIdx, i) );
349  };
350  return result;
351  }
352 }; // namespace Opm
353 
354 } // end namespace Opm
355 
356 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
static void solve(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache, unsigned phaseIdx, const ComponentVector &targetFug)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: CompositionFromFugacities.hpp:88
static void guessInitial(FluidState &fluidState, unsigned phaseIdx, const ComponentVector &)
Guess an initial value for the composition of the phase.
Definition: CompositionFromFugacities.hpp:65
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: CompositionFromFugacities.hpp:54