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>
45template <
class Scalar,
class Flu
idSystem,
class Evaluation = Scalar>
48 enum { numComponents = FluidSystem::numComponents };
56 template <
class Flu
idState>
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,
87 if (FluidSystem::isIdealMixture(phaseIdx)) {
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);
132 try { J.solve(x, b); }
133 catch (
const Dune::FMatrixError& e)
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>
184 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
188 for (
unsigned i = 0; i < numComponents; ++ i) {
189 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
193 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
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,
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;
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>
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);
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>
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 Scalar update_(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache, Dune::FieldVector< Evaluation, numComponents > &x, Dune::FieldVector< Evaluation, numComponents > &, unsigned phaseIdx, const Dune::FieldVector< Evaluation, numComponents > &targetFug)
Definition: CompositionFromFugacities.hpp:281
static Scalar calculateDefect_(const FluidState ¶ms, unsigned phaseIdx, const ComponentVector &targetFug)
Definition: CompositionFromFugacities.hpp:328
static Scalar linearize_(Dune::FieldMatrix< Evaluation, numComponents, numComponents > &J, Dune::FieldVector< Evaluation, numComponents > &defect, FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache, unsigned phaseIdx, const ComponentVector &targetFug)
Definition: CompositionFromFugacities.hpp:209
Dune::FieldVector< Evaluation, numComponents > ComponentVector
Definition: CompositionFromFugacities.hpp:51
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 solveIdealMix_(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache, unsigned phaseIdx, const ComponentVector &fugacities)
Definition: CompositionFromFugacities.hpp:183
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
bool CheckDefined(const T &value)
Make valgrind complain if any of the memory occupied by an object is undefined.
Definition: Valgrind.hpp:74
Definition: Air_Mesitylene.hpp:34
ReturnEval_< Evaluation1, Evaluation2 >::type min(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:346
ReturnEval_< Evaluation1, Evaluation2 >::type max(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:341
auto scalarValue(const Evaluation &val) -> decltype(MathToolbox< Evaluation >::scalarValue(val))
Definition: MathToolbox.hpp:335
Evaluation abs(const Evaluation &value)
Definition: MathToolbox.hpp:350