00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_COMPOSITION_FROM_FUGACITIES_HPP
00028 #define OPM_COMPOSITION_FROM_FUGACITIES_HPP
00029
00030 #include <opm/material/common/MathToolbox.hpp>
00031
00032
00033 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00034
00035 #include <dune/common/fvector.hh>
00036 #include <dune/common/fmatrix.hh>
00037
00038 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00039
00040
00041 #include <opm/common/ErrorMacros.hpp>
00042 #include <opm/common/Exceptions.hpp>
00043 #include <opm/common/Valgrind.hpp>
00044
00045 #include <limits>
00046
00047 namespace Opm {
00048
00053 template <class Scalar, class FluidSystem, class Evaluation = Scalar>
00054 class CompositionFromFugacities
00055 {
00056 enum { numComponents = FluidSystem::numComponents };
00057
00058 public:
00059 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
00060
00064 template <class FluidState>
00065 static void guessInitial(FluidState& fluidState,
00066 unsigned phaseIdx,
00067 const ComponentVector& )
00068 {
00069 if (FluidSystem::isIdealMixture(phaseIdx))
00070 return;
00071
00072
00073 for (unsigned i = 0; i < numComponents; ++ i) {
00074
00075 fluidState.setMoleFraction(phaseIdx,
00076 i,
00077 1.0/numComponents);
00078 }
00079 }
00080
00087 template <class FluidState>
00088 static void solve(FluidState& fluidState,
00089 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
00090 unsigned phaseIdx,
00091 const ComponentVector& targetFug)
00092 {
00093
00094
00095 if (FluidSystem::isIdealMixture(phaseIdx)) {
00096 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
00097 return;
00098 }
00099
00100
00101
00102
00103 Dune::FieldVector<Evaluation, numComponents> xInit;
00104 for (unsigned i = 0; i < numComponents; ++i) {
00105 xInit[i] = fluidState.moleFraction(phaseIdx, i);
00106 }
00107
00109
00111
00112
00113 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
00114
00115 Dune::FieldVector<Evaluation, numComponents> x;
00116
00117 Dune::FieldVector<Evaluation, numComponents> b;
00118
00119 paramCache.updatePhase(fluidState, phaseIdx);
00120
00121
00122 const int nMax = 25;
00123 for (int nIdx = 0; nIdx < nMax; ++nIdx) {
00124
00125 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
00126 Valgrind::CheckDefined(J);
00127 Valgrind::CheckDefined(b);
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 x = 0.0;
00142 try { J.solve(x, b); }
00143 catch (Dune::FMatrixError e)
00144 { throw Opm::NumericalProblem(e.what()); }
00145
00146
00147
00148 Valgrind::CheckDefined(x);
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
00168
00169 if (relError < 1e-9) {
00170 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
00171 fluidState.setDensity(phaseIdx, rho);
00172
00173
00174 return;
00175 }
00176 }
00177
00178 OPM_THROW(Opm::NumericalProblem,
00179 "Calculating the " << FluidSystem::phaseName(phaseIdx)
00180 << "Phase composition failed. Initial {x} = {"
00181 << xInit
00182 << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx)
00183 << ", T = " << fluidState.temperature(phaseIdx));
00184 }
00185
00186
00187 protected:
00188
00189
00190
00191 template <class FluidState>
00192 static void solveIdealMix_(FluidState& fluidState,
00193 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
00194 unsigned phaseIdx,
00195 const ComponentVector& fugacities)
00196 {
00197 for (unsigned i = 0; i < numComponents; ++ i) {
00198 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
00199 paramCache,
00200 phaseIdx,
00201 i);
00202 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
00203 Valgrind::CheckDefined(phi);
00204 Valgrind::CheckDefined(gamma);
00205 Valgrind::CheckDefined(fugacities[i]);
00206 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
00207 fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
00208 };
00209
00210 paramCache.updatePhase(fluidState, phaseIdx);
00211
00212 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
00213 fluidState.setDensity(phaseIdx, rho);
00214 return;
00215 }
00216
00217 template <class FluidState>
00218 static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents>& J,
00219 Dune::FieldVector<Evaluation, numComponents>& defect,
00220 FluidState& fluidState,
00221 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
00222 unsigned phaseIdx,
00223 const ComponentVector& targetFug)
00224 {
00225
00226 J = 0;
00227
00228 Scalar absError = 0;
00229
00230
00231 for (unsigned i = 0; i < numComponents; ++ i) {
00232 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
00233 paramCache,
00234 phaseIdx,
00235 i);
00236 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
00237 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
00238
00239 defect[i] = targetFug[i] - f;
00240 absError = std::max(absError, std::abs(Opm::scalarValue(defect[i])));
00241 }
00242
00243
00244 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
00245 for (unsigned i = 0; i < numComponents; ++ i) {
00247
00248
00249
00250
00251
00252
00253 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
00254 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
00255 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
00256
00257
00258
00259 for (unsigned j = 0; j < numComponents; ++j) {
00260
00261 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
00262 paramCache,
00263 phaseIdx,
00264 j);
00265
00266 const Evaluation& f =
00267 phi *
00268 fluidState.pressure(phaseIdx) *
00269 fluidState.moleFraction(phaseIdx, j);
00270
00271 const Evaluation& defJPlusEps = targetFug[j] - f;
00272
00273
00274
00275 J[j][i] = (defJPlusEps - defect[j])/eps;
00276 }
00277
00278
00279 fluidState.setMoleFraction(phaseIdx, i, xI);
00280 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
00281
00282
00284 }
00285
00286 return absError;
00287 }
00288
00289 template <class FluidState>
00290 static Scalar update_(FluidState& fluidState,
00291 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
00292 Dune::FieldVector<Evaluation, numComponents>& x,
00293 Dune::FieldVector<Evaluation, numComponents>& ,
00294 unsigned phaseIdx,
00295 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
00296 {
00297
00298 Dune::FieldVector<Evaluation, numComponents> origComp;
00299 Scalar relError = 0;
00300 Evaluation sumDelta = 0.0;
00301 Evaluation sumx = 0.0;
00302 for (unsigned i = 0; i < numComponents; ++i) {
00303 origComp[i] = fluidState.moleFraction(phaseIdx, i);
00304 relError = std::max(relError, std::abs(Opm::scalarValue(x[i])));
00305
00306 sumx += Opm::abs(fluidState.moleFraction(phaseIdx, i));
00307 sumDelta += Opm::abs(x[i]);
00308 }
00309
00310
00311 const Scalar maxDelta = 0.2;
00312 if (sumDelta > maxDelta)
00313 x /= (sumDelta/maxDelta);
00314
00315
00316 for (unsigned i = 0; i < numComponents; ++i) {
00317 Evaluation newx = origComp[i] - x[i];
00318
00319 if (targetFug[i] > 0)
00320 newx = Opm::max(0.0, newx);
00321
00322 else if (targetFug[i] < 0)
00323 newx = Opm::min(0.0, newx);
00324
00325 else
00326 newx = 0;
00327
00328 fluidState.setMoleFraction(phaseIdx, i, newx);
00329 }
00330
00331 paramCache.updateComposition(fluidState, phaseIdx);
00332
00333 return relError;
00334 }
00335
00336 template <class FluidState>
00337 static Scalar calculateDefect_(const FluidState& params,
00338 unsigned phaseIdx,
00339 const ComponentVector& targetFug)
00340 {
00341 Scalar result = 0.0;
00342 for (unsigned i = 0; i < numComponents; ++i) {
00343
00344
00345 result += std::abs(
00346 (targetFug[i] - params.fugacity(phaseIdx, i))
00347 /
00348 params.fugacityCoefficient(phaseIdx, i) );
00349 };
00350 return result;
00351 }
00352 };
00353
00354 }
00355
00356 #endif