27 #ifndef OPM_COMPOSITION_FROM_FUGACITIES_HPP 28 #define OPM_COMPOSITION_FROM_FUGACITIES_HPP 33 #include <opm/common/utility/platform_dependent/disable_warnings.h> 35 #include <dune/common/fvector.hh> 36 #include <dune/common/fmatrix.hh> 38 #include <opm/common/utility/platform_dependent/reenable_warnings.h> 41 #include <opm/common/ErrorMacros.hpp> 42 #include <opm/common/Exceptions.hpp> 43 #include <opm/common/Valgrind.hpp> 53 template <
class Scalar,
class Flu
idSystem,
class Evaluation = Scalar>
56 enum { numComponents = FluidSystem::numComponents };
59 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
64 template <
class Flu
idState>
67 const ComponentVector& )
69 if (FluidSystem::isIdealMixture(phaseIdx))
73 for (
unsigned i = 0; i < numComponents; ++ i) {
75 fluidState.setMoleFraction(phaseIdx,
87 template <
class Flu
idState>
88 static void solve(FluidState& fluidState,
89 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
91 const ComponentVector& targetFug)
95 if (FluidSystem::isIdealMixture(phaseIdx)) {
96 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
103 Dune::FieldVector<Evaluation, numComponents> xInit;
104 for (
unsigned i = 0; i < numComponents; ++i) {
105 xInit[i] = fluidState.moleFraction(phaseIdx, i);
113 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
115 Dune::FieldVector<Evaluation, numComponents> x;
117 Dune::FieldVector<Evaluation, numComponents> b;
119 paramCache.updatePhase(fluidState, phaseIdx);
123 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
125 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
126 Valgrind::CheckDefined(J);
127 Valgrind::CheckDefined(b);
142 try { J.solve(x, b); }
143 catch (Dune::FMatrixError e)
144 {
throw Opm::NumericalProblem(e.what()); }
148 Valgrind::CheckDefined(x);
167 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
169 if (relError < 1e-9) {
170 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
171 fluidState.setDensity(phaseIdx, rho);
178 OPM_THROW(Opm::NumericalProblem,
179 "Calculating the " << FluidSystem::phaseName(phaseIdx)
180 <<
"Phase composition failed. Initial {x} = {" 182 <<
"}, {fug_t} = {" << targetFug <<
"}, p = " << fluidState.pressure(phaseIdx)
183 <<
", T = " << fluidState.temperature(phaseIdx));
191 template <
class Flu
idState>
192 static void solveIdealMix_(FluidState& fluidState,
193 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
195 const ComponentVector& fugacities)
197 for (
unsigned i = 0; i < numComponents; ++ i) {
198 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
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);
210 paramCache.updatePhase(fluidState, phaseIdx);
212 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
213 fluidState.setDensity(phaseIdx, rho);
217 template <
class Flu
idState>
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,
223 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;
240 absError = std::max(absError, std::abs(Opm::scalarValue(defect[i])));
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>
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>& ,
295 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
298 Dune::FieldVector<Evaluation, numComponents> origComp;
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])));
306 sumx += Opm::abs(fluidState.moleFraction(phaseIdx, i));
307 sumDelta += Opm::abs(x[i]);
311 const Scalar maxDelta = 0.2;
312 if (sumDelta > maxDelta)
313 x /= (sumDelta/maxDelta);
316 for (
unsigned i = 0; i < numComponents; ++i) {
317 Evaluation newx = origComp[i] - x[i];
319 if (targetFug[i] > 0)
320 newx = Opm::max(0.0, newx);
322 else if (targetFug[i] < 0)
323 newx = Opm::min(0.0, newx);
328 fluidState.setMoleFraction(phaseIdx, i, newx);
331 paramCache.updateComposition(fluidState, phaseIdx);
336 template <
class Flu
idState>
337 static Scalar calculateDefect_(
const FluidState& params,
339 const ComponentVector& targetFug)
342 for (
unsigned i = 0; i < numComponents; ++i) {
346 (targetFug[i] - params.fugacity(phaseIdx, i))
348 params.fugacityCoefficient(phaseIdx, i) );
static void solve(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache, 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
Definition: Air_Mesitylene.hpp:33
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: CompositionFromFugacities.hpp:54