49 enum { numComponents = FluidSystem::numComponents };
52 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
57 template <
class Flu
idState>
60 const ComponentVector& )
62 if (FluidSystem::isIdealMixture(phaseIdx))
66 for (
unsigned i = 0; i < numComponents; ++ i) {
68 fluidState.setMoleFraction(phaseIdx,
80 template <
class Flu
idState>
81 static void solve(FluidState& fluidState,
82 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
84 const ComponentVector& targetFug)
88 if (FluidSystem::isIdealMixture(phaseIdx)) {
89 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
94 Dune::FieldVector<Evaluation, numComponents> xInit;
95 for (
unsigned i = 0; i < numComponents; ++i) {
96 xInit[i] = fluidState.moleFraction(phaseIdx, i);
104 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
106 Dune::FieldVector<Evaluation, numComponents> x;
108 Dune::FieldVector<Evaluation, numComponents> b;
110 paramCache.updatePhase(fluidState, phaseIdx);
114 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
116 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
117 Valgrind::CheckDefined(J);
118 Valgrind::CheckDefined(b);
133 try { J.solve(x, b); }
134 catch (
const Dune::FMatrixError& e)
139 Valgrind::CheckDefined(x);
158 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
160 if (relError < 1e-9) {
161 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
162 fluidState.setDensity(phaseIdx, rho);
169 auto cast = [](
const auto d)
172 if constexpr (std::is_same_v<
decltype(d),
const quad>)
173 return static_cast<double>(d);
180 std::string(
"Calculating the ")
181 + FluidSystem::phaseName(phaseIdx)
182 +
"Phase composition failed. Initial {x} = {";
183 for (
const auto& v : xInit)
184 msg +=
" " + std::to_string(cast(getValue(v)));
185 msg +=
" }, {fug_t} = {";
186 for (
const auto& v : targetFug)
187 msg +=
" " + std::to_string(cast(getValue(v)));
188 msg +=
" }, p = " + std::to_string(cast(getValue(fluidState.pressure(phaseIdx))))
189 +
", T = " + std::to_string(cast(getValue(fluidState.temperature(phaseIdx))));
198 template <
class Flu
idState>
199 static void solveIdealMix_(FluidState& fluidState,
200 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
202 const ComponentVector& fugacities)
204 for (
unsigned i = 0; i < numComponents; ++ i) {
205 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
209 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
210 Valgrind::CheckDefined(phi);
211 Valgrind::CheckDefined(gamma);
212 Valgrind::CheckDefined(fugacities[i]);
213 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
214 fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
217 paramCache.updatePhase(fluidState, phaseIdx);
219 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
220 fluidState.setDensity(phaseIdx, rho);
224 template <
class Flu
idState>
225 static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents>& J,
226 Dune::FieldVector<Evaluation, numComponents>& defect,
227 FluidState& fluidState,
228 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
230 const ComponentVector& targetFug)
238 for (
unsigned i = 0; i < numComponents; ++ i) {
239 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
243 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
244 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
246 defect[i] = targetFug[i] - f;
247 absError = std::max(absError, std::abs(scalarValue(defect[i])));
251 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
252 for (
unsigned i = 0; i < numComponents; ++ i) {
260 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
261 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
262 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
266 for (
unsigned j = 0; j < numComponents; ++j) {
268 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
273 const Evaluation& f =
275 fluidState.pressure(phaseIdx) *
276 fluidState.moleFraction(phaseIdx, j);
278 const Evaluation& defJPlusEps = targetFug[j] - f;
282 J[j][i] = (defJPlusEps - defect[j])/eps;
286 fluidState.setMoleFraction(phaseIdx, i, xI);
287 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
296 template <
class Flu
idState>
297 static Scalar update_(FluidState& fluidState,
298 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
299 Dune::FieldVector<Evaluation, numComponents>& x,
300 Dune::FieldVector<Evaluation, numComponents>& ,
302 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
305 Dune::FieldVector<Evaluation, numComponents> origComp;
307 Evaluation sumDelta = 0.0;
308 Evaluation sumx = 0.0;
309 for (
unsigned i = 0; i < numComponents; ++i) {
310 origComp[i] = fluidState.moleFraction(phaseIdx, i);
311 relError = std::max(relError, std::abs(scalarValue(x[i])));
313 sumx += abs(fluidState.moleFraction(phaseIdx, i));
314 sumDelta += abs(x[i]);
318 const Scalar maxDelta = 0.2;
319 if (sumDelta > maxDelta)
320 x /= (sumDelta/maxDelta);
323 for (
unsigned i = 0; i < numComponents; ++i) {
324 Evaluation newx = origComp[i] - x[i];
326 if (targetFug[i] > 0)
327 newx = max(0.0, newx);
329 else if (targetFug[i] < 0)
330 newx = min(0.0, newx);
335 fluidState.setMoleFraction(phaseIdx, i, newx);
338 paramCache.updateComposition(fluidState, phaseIdx);
343 template <
class Flu
idState>
344 static Scalar calculateDefect_(
const FluidState& params,
346 const ComponentVector& targetFug)
349 for (
unsigned i = 0; i < numComponents; ++i) {
353 (targetFug[i] - params.fugacity(phaseIdx, i))
355 params.fugacityCoefficient(phaseIdx, i) );