26 #ifndef OPM_COMPOSITION_FROM_FUGACITIES_HPP
27 #define OPM_COMPOSITION_FROM_FUGACITIES_HPP
31 #include <dune/common/fvector.hh>
32 #include <dune/common/fmatrix.hh>
34 #include <opm/common/ErrorMacros.hpp>
35 #include <opm/common/Exceptions.hpp>
46 template <
class Scalar,
class Flu
idSystem,
class Evaluation = Scalar>
49 enum { numComponents = FluidSystem::numComponents };
51 typedef typename FluidSystem::ParameterCache ParameterCache;
59 template <
class Flu
idState>
63 const ComponentVector &)
65 if (FluidSystem::isIdealMixture(phaseIdx))
69 for (
unsigned i = 0; i < numComponents; ++ i) {
71 fluidState.setMoleFraction(phaseIdx,
83 template <
class Flu
idState>
84 static void solve(FluidState &fluidState,
85 ParameterCache ¶mCache,
87 const ComponentVector &targetFug)
93 if (FluidSystem::isIdealMixture(phaseIdx)) {
101 Dune::FieldVector<Evaluation, numComponents> xInit;
102 for (
unsigned i = 0; i < numComponents; ++i) {
103 xInit[i] = fluidState.moleFraction(phaseIdx, i);
111 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
113 Dune::FieldVector<Evaluation, numComponents> x;
115 Dune::FieldVector<Evaluation, numComponents> b;
117 paramCache.updatePhase(fluidState, phaseIdx);
121 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
123 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
139 x = Toolbox::createConstant(0.0);
140 try { J.solve(x, b); }
141 catch (Dune::FMatrixError e)
142 {
throw Opm::NumericalProblem(e.what()); }
165 Scalar relError =
update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
167 if (relError < 1e-9) {
168 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
169 fluidState.setDensity(phaseIdx, rho);
176 OPM_THROW(Opm::NumericalProblem,
177 "Calculating the " << FluidSystem::phaseName(phaseIdx)
178 <<
"Phase composition failed. Initial {x} = {"
180 <<
"}, {fug_t} = {" << targetFug <<
"}, p = " << fluidState.pressure(phaseIdx)
181 <<
", T = " << fluidState.temperature(phaseIdx));
189 template <
class Flu
idState>
191 ParameterCache ¶mCache,
193 const ComponentVector &fugacities)
195 for (
unsigned i = 0; i < numComponents; ++ i) {
196 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
200 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
204 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
205 fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
208 paramCache.updatePhase(fluidState, phaseIdx);
210 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
211 fluidState.setDensity(phaseIdx, rho);
215 template <
class Flu
idState>
216 static Scalar
linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents> &J,
217 Dune::FieldVector<Evaluation, numComponents> &defect,
218 FluidState &fluidState,
219 ParameterCache ¶mCache,
221 const ComponentVector &targetFug)
231 for (
unsigned i = 0; i < numComponents; ++ i) {
232 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
236 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
237 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
239 defect[i] = targetFug[i] - f;
244 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
245 for (
unsigned i = 0; i < numComponents; ++ i) {
253 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
254 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
255 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
259 for (
unsigned j = 0; j < numComponents; ++j) {
261 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
266 const Evaluation& f =
268 fluidState.pressure(phaseIdx) *
269 fluidState.moleFraction(phaseIdx, j);
271 const Evaluation& defJPlusEps = targetFug[j] - f;
275 J[j][i] = (defJPlusEps - defect[j])/eps;
279 fluidState.setMoleFraction(phaseIdx, i, xI);
280 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
289 template <
class Flu
idState>
291 ParameterCache ¶mCache,
292 Dune::FieldVector<Evaluation, numComponents> &x,
293 Dune::FieldVector<Evaluation, numComponents> &,
295 const Dune::FieldVector<Evaluation, numComponents> &targetFug)
300 Dune::FieldVector<Evaluation, numComponents> origComp;
302 Evaluation sumDelta = Toolbox::createConstant(0.0);
303 Evaluation sumx = Toolbox::createConstant(0.0);
304 for (
unsigned i = 0; i < numComponents; ++i) {
305 origComp[i] = fluidState.moleFraction(phaseIdx, i);
308 sumx +=
Toolbox::abs(fluidState.moleFraction(phaseIdx, i));
313 const Scalar maxDelta = 0.2;
314 if (sumDelta > maxDelta)
315 x /= (sumDelta/maxDelta);
318 for (
unsigned i = 0; i < numComponents; ++i) {
319 Evaluation newx = origComp[i] - x[i];
321 if (targetFug[i] > 0)
324 else if (targetFug[i] < 0)
330 fluidState.setMoleFraction(phaseIdx, i, newx);
333 paramCache.updateComposition(fluidState, phaseIdx);
338 template <
class Flu
idState>
341 const ComponentVector &targetFug)
344 for (
unsigned i = 0; i < numComponents; ++i) {
348 (targetFug[i] - params.fugacity(phaseIdx, i))
350 params.fugacityCoefficient(phaseIdx, i) );
static Scalar update_(FluidState &fluidState, ParameterCache ¶mCache, Dune::FieldVector< Evaluation, numComponents > &x, Dune::FieldVector< Evaluation, numComponents > &, unsigned phaseIdx, const Dune::FieldVector< Evaluation, numComponents > &targetFug)
Definition: CompositionFromFugacities.hpp:290
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: CompositionFromFugacities.hpp:47
bool CheckDefined(const T &value OPM_UNUSED)
Make valgrind complain if any of the memory occupied by an object is undefined.
Definition: Valgrind.hpp:74
Definition: Air_Mesitylene.hpp:31
static Scalar calculateDefect_(const FluidState ¶ms, unsigned phaseIdx, const ComponentVector &targetFug)
Definition: CompositionFromFugacities.hpp:339
Some templates to wrap the valgrind client request macros.
Evaluation< Scalar, VarSetTag, numVars > max(const Evaluation< Scalar, VarSetTag, numVars > &x1, const Evaluation< Scalar, VarSetTag, numVars > &x2)
Definition: Math.hpp:114
static void solve(FluidState &fluidState, ParameterCache ¶mCache, unsigned phaseIdx, const ComponentVector &targetFug)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: CompositionFromFugacities.hpp:84
Dune::FieldVector< Evaluation, numComponents > ComponentVector
Definition: CompositionFromFugacities.hpp:54
static Scalar linearize_(Dune::FieldMatrix< Evaluation, numComponents, numComponents > &J, Dune::FieldVector< Evaluation, numComponents > &defect, FluidState &fluidState, ParameterCache ¶mCache, unsigned phaseIdx, const ComponentVector &targetFug)
Definition: CompositionFromFugacities.hpp:216
Evaluation< Scalar, VarSetTag, numVars > abs(const Evaluation< Scalar, VarSetTag, numVars > &)
Definition: Math.hpp:41
Evaluation< Scalar, VarSetTag, numVars > min(const Evaluation< Scalar, VarSetTag, numVars > &x1, const Evaluation< Scalar, VarSetTag, numVars > &x2)
Definition: Math.hpp:61
static void solveIdealMix_(FluidState &fluidState, ParameterCache ¶mCache, unsigned phaseIdx, const ComponentVector &fugacities)
Definition: CompositionFromFugacities.hpp:190
static void guessInitial(FluidState &fluidState, ParameterCache &, unsigned phaseIdx, const ComponentVector &)
Guess an initial value for the composition of the phase.
Definition: CompositionFromFugacities.hpp:60