27 #ifndef OPM_COMPOSITION_FROM_FUGACITIES_HPP
28 #define OPM_COMPOSITION_FROM_FUGACITIES_HPP
34 #include <dune/common/fvector.hh>
35 #include <dune/common/fmatrix.hh>
45 template <
class Scalar,
class Flu
idSystem,
class Evaluation = Scalar>
48 enum { numComponents = FluidSystem::numComponents };
51 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
56 template <
class Flu
idState>
59 const ComponentVector& )
61 if (FluidSystem::isIdealMixture(phaseIdx))
65 for (
unsigned i = 0; i < numComponents; ++ i) {
67 fluidState.setMoleFraction(phaseIdx,
79 template <
class Flu
idState>
80 static void solve(FluidState& fluidState,
81 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
83 const ComponentVector& targetFug)
87 if (FluidSystem::isIdealMixture(phaseIdx)) {
88 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
93 Dune::FieldVector<Evaluation, numComponents> xInit;
94 for (
unsigned i = 0; i < numComponents; ++i) {
95 xInit[i] = fluidState.moleFraction(phaseIdx, i);
103 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
105 Dune::FieldVector<Evaluation, numComponents> x;
107 Dune::FieldVector<Evaluation, numComponents> b;
109 paramCache.updatePhase(fluidState, phaseIdx);
113 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
115 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
116 Valgrind::CheckDefined(J);
117 Valgrind::CheckDefined(b);
132 try { J.solve(x, b); }
133 catch (
const Dune::FMatrixError& e)
138 Valgrind::CheckDefined(x);
157 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
159 if (relError < 1e-9) {
160 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
161 fluidState.setDensity(phaseIdx, rho);
168 std::ostringstream oss;
169 oss <<
"Calculating the " << FluidSystem::phaseName(phaseIdx)
170 <<
"Phase composition failed. Initial {x} = {"
172 <<
"}, {fug_t} = {" << targetFug <<
"}, p = " << fluidState.pressure(phaseIdx)
173 <<
", T = " << fluidState.temperature(phaseIdx);
182 template <
class Flu
idState>
183 static void solveIdealMix_(FluidState& fluidState,
184 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
186 const ComponentVector& fugacities)
188 for (
unsigned i = 0; i < numComponents; ++ i) {
189 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
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);
201 paramCache.updatePhase(fluidState, phaseIdx);
203 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
204 fluidState.setDensity(phaseIdx, rho);
208 template <
class Flu
idState>
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,
214 const ComponentVector& targetFug)
222 for (
unsigned i = 0; i < numComponents; ++ i) {
223 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
227 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
228 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
230 defect[i] = targetFug[i] - f;
231 absError = std::max(absError, std::abs(scalarValue(defect[i])));
235 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
236 for (
unsigned i = 0; i < numComponents; ++ i) {
244 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
245 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
246 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
250 for (
unsigned j = 0; j < numComponents; ++j) {
252 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
257 const Evaluation& f =
259 fluidState.pressure(phaseIdx) *
260 fluidState.moleFraction(phaseIdx, j);
262 const Evaluation& defJPlusEps = targetFug[j] - f;
266 J[j][i] = (defJPlusEps - defect[j])/eps;
270 fluidState.setMoleFraction(phaseIdx, i, xI);
271 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
280 template <
class Flu
idState>
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>& ,
286 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
289 Dune::FieldVector<Evaluation, numComponents> origComp;
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])));
297 sumx += abs(fluidState.moleFraction(phaseIdx, i));
298 sumDelta += abs(x[i]);
302 const Scalar maxDelta = 0.2;
303 if (sumDelta > maxDelta)
304 x /= (sumDelta/maxDelta);
307 for (
unsigned i = 0; i < numComponents; ++i) {
308 Evaluation newx = origComp[i] - x[i];
310 if (targetFug[i] > 0)
311 newx = max(0.0, newx);
313 else if (targetFug[i] < 0)
314 newx = min(0.0, newx);
319 fluidState.setMoleFraction(phaseIdx, i, newx);
322 paramCache.updateComposition(fluidState, phaseIdx);
327 template <
class Flu
idState>
328 static Scalar calculateDefect_(
const FluidState& params,
330 const ComponentVector& targetFug)
333 for (
unsigned i = 0; i < numComponents; ++i) {
337 (targetFug[i] - params.fugacity(phaseIdx, i))
339 params.fugacityCoefficient(phaseIdx, i) );
Provides the opm-material specific exception classes.
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 > ¶mCache, 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