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