50 enum { numComponents = FluidSystem::numComponents };
53 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
58 template <
class Flu
idState>
61 const ComponentVector& )
63 if (FluidSystem::isIdealMixture(phaseIdx))
67 for (
unsigned i = 0; i < numComponents; ++ i) {
69 fluidState.setMoleFraction(phaseIdx,
81 template <
class Flu
idState>
82 static void solve(FluidState& fluidState,
83 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
85 const ComponentVector& targetFug)
89 if (FluidSystem::isIdealMixture(phaseIdx)) {
90 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
95 Dune::FieldVector<Evaluation, numComponents> xInit;
96 for (
unsigned i = 0; i < numComponents; ++i) {
97 xInit[i] = fluidState.moleFraction(phaseIdx, i);
105 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
107 Dune::FieldVector<Evaluation, numComponents> x;
109 Dune::FieldVector<Evaluation, numComponents> b;
111 paramCache.updatePhase(fluidState, phaseIdx);
115 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
117 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
118 Valgrind::CheckDefined(J);
119 Valgrind::CheckDefined(b);
134 try { J.solve(x, b); }
135 catch (
const Dune::FMatrixError& e)
140 Valgrind::CheckDefined(x);
159 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
161 if (relError < 1e-9) {
162 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
163 fluidState.setDensity(phaseIdx, rho);
170 auto cast = [](
const auto d)
173 if constexpr (std::is_same_v<
decltype(d),
const quad>)
174 return static_cast<double>(d);
181 std::string(
"Calculating the ")
182 + FluidSystem::phaseName(phaseIdx).data()
183 +
"Phase composition failed. Initial {x} = {";
184 for (
const auto& v : xInit)
185 msg +=
" " + std::to_string(cast(getValue(v)));
186 msg +=
" }, {fug_t} = {";
187 for (
const auto& v : targetFug)
188 msg +=
" " + std::to_string(cast(getValue(v)));
189 msg +=
" }, p = " + std::to_string(cast(getValue(fluidState.pressure(phaseIdx))))
190 +
", T = " + std::to_string(cast(getValue(fluidState.temperature(phaseIdx))));
199 template <
class Flu
idState>
200 static void solveIdealMix_(FluidState& fluidState,
201 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
203 const ComponentVector& fugacities)
205 for (
unsigned i = 0; i < numComponents; ++ i) {
206 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
210 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
211 Valgrind::CheckDefined(phi);
212 Valgrind::CheckDefined(gamma);
213 Valgrind::CheckDefined(fugacities[i]);
214 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
215 fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
218 paramCache.updatePhase(fluidState, phaseIdx);
220 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
221 fluidState.setDensity(phaseIdx, rho);
225 template <
class Flu
idState>
226 static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents>& J,
227 Dune::FieldVector<Evaluation, numComponents>& defect,
228 FluidState& fluidState,
229 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
231 const ComponentVector& targetFug)
239 for (
unsigned i = 0; i < numComponents; ++ i) {
240 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
244 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
245 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
247 defect[i] = targetFug[i] - f;
248 absError = std::max(absError, std::abs(scalarValue(defect[i])));
252 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
253 for (
unsigned i = 0; i < numComponents; ++ i) {
261 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
262 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
263 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
267 for (
unsigned j = 0; j < numComponents; ++j) {
269 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
274 const Evaluation& f =
276 fluidState.pressure(phaseIdx) *
277 fluidState.moleFraction(phaseIdx, j);
279 const Evaluation& defJPlusEps = targetFug[j] - f;
283 J[j][i] = (defJPlusEps - defect[j])/eps;
287 fluidState.setMoleFraction(phaseIdx, i, xI);
288 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
297 template <
class Flu
idState>
298 static Scalar update_(FluidState& fluidState,
299 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
300 Dune::FieldVector<Evaluation, numComponents>& x,
301 Dune::FieldVector<Evaluation, numComponents>& ,
303 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
306 Dune::FieldVector<Evaluation, numComponents> origComp;
308 Evaluation sumDelta = 0.0;
309 Evaluation sumx = 0.0;
310 for (
unsigned i = 0; i < numComponents; ++i) {
311 origComp[i] = fluidState.moleFraction(phaseIdx, i);
312 relError = std::max(relError, std::abs(scalarValue(x[i])));
314 sumx += abs(fluidState.moleFraction(phaseIdx, i));
315 sumDelta += abs(x[i]);
319 const Scalar maxDelta = 0.2;
320 if (sumDelta > maxDelta)
321 x /= (sumDelta/maxDelta);
324 for (
unsigned i = 0; i < numComponents; ++i) {
325 Evaluation newx = origComp[i] - x[i];
327 if (targetFug[i] > 0)
328 newx = max(0.0, newx);
330 else if (targetFug[i] < 0)
331 newx = min(0.0, newx);
336 fluidState.setMoleFraction(phaseIdx, i, newx);
339 paramCache.updateComposition(fluidState, phaseIdx);
344 template <
class Flu
idState>
345 static Scalar calculateDefect_(
const FluidState& params,
347 const ComponentVector& targetFug)
350 for (
unsigned i = 0; i < numComponents; ++i) {
354 (targetFug[i] - params.fugacity(phaseIdx, i))
356 params.fugacityCoefficient(phaseIdx, i) );