51 static constexpr int numComponents = FluidSystem::numComponents;
54 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
59 template <
class Flu
idState>
62 const ComponentVector& )
64 if (FluidSystem::isIdealMixture(phaseIdx))
68 for (
unsigned i = 0; i < numComponents; ++ i) {
70 fluidState.setMoleFraction(phaseIdx,
82 template <
class Flu
idState>
83 static void solve(FluidState& fluidState,
84 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
86 const ComponentVector& targetFug)
88 assert (phaseIdx <
static_cast<unsigned int>(FluidSystem::numPhases));
92 if (FluidSystem::isIdealMixture(phaseIdx)) {
93 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
98 Dune::FieldVector<Evaluation, numComponents> xInit;
99 for (
unsigned i = 0; i < numComponents; ++i) {
100 xInit[i] = fluidState.moleFraction(phaseIdx, i);
108 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
110 Dune::FieldVector<Evaluation, numComponents> x;
112 Dune::FieldVector<Evaluation, numComponents> b;
114 paramCache.updatePhase(fluidState, phaseIdx);
118 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
120 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
121 Valgrind::CheckDefined(J);
122 Valgrind::CheckDefined(b);
137 try { J.solve(x, b); }
138 catch (
const Dune::FMatrixError& e)
143 Valgrind::CheckDefined(x);
162 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
164 if (relError < 1e-9) {
165 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
166 fluidState.setDensity(phaseIdx, rho);
173 auto cast = [](
const auto d)
176 if constexpr (std::is_same_v<
decltype(d),
const quad>)
177 return static_cast<double>(d);
184 std::string(
"Calculating the ")
185 + FluidSystem::phaseName(phaseIdx).data()
186 +
"Phase composition failed. Initial {x} = {";
187 for (
const auto& v : xInit)
188 msg +=
" " + std::to_string(cast(getValue(v)));
189 msg +=
" }, {fug_t} = {";
190 for (
const auto& v : targetFug)
191 msg +=
" " + std::to_string(cast(getValue(v)));
192 msg +=
" }, p = " + std::to_string(cast(getValue(fluidState.pressure(phaseIdx))))
193 +
", T = " + std::to_string(cast(getValue(fluidState.temperature(phaseIdx))));
202 template <
class Flu
idState>
203 static void solveIdealMix_(FluidState& fluidState,
204 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
206 const ComponentVector& fugacities)
208 for (
unsigned i = 0; i < numComponents; ++ i) {
209 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
213 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
214 Valgrind::CheckDefined(phi);
215 Valgrind::CheckDefined(gamma);
216 Valgrind::CheckDefined(fugacities[i]);
217 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
218 fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
221 paramCache.updatePhase(fluidState, phaseIdx);
223 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
224 fluidState.setDensity(phaseIdx, rho);
228 template <
class Flu
idState>
229 static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents>& J,
230 Dune::FieldVector<Evaluation, numComponents>& defect,
231 FluidState& fluidState,
232 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
234 const ComponentVector& targetFug)
242 for (
unsigned i = 0; i < numComponents; ++ i) {
243 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
247 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
248 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
250 defect[i] = targetFug[i] - f;
251 absError = std::max(absError, std::abs(scalarValue(defect[i])));
255 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
256 for (
unsigned i = 0; i < numComponents; ++ i) {
264 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
265 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
266 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
270 for (
unsigned j = 0; j < numComponents; ++j) {
272 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
277 const Evaluation& f =
279 fluidState.pressure(phaseIdx) *
280 fluidState.moleFraction(phaseIdx, j);
282 const Evaluation& defJPlusEps = targetFug[j] - f;
286 J[j][i] = (defJPlusEps - defect[j])/eps;
290 fluidState.setMoleFraction(phaseIdx, i, xI);
291 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
300 template <
class Flu
idState>
301 static Scalar update_(FluidState& fluidState,
302 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
303 Dune::FieldVector<Evaluation, numComponents>& x,
304 Dune::FieldVector<Evaluation, numComponents>& ,
306 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
309 Dune::FieldVector<Evaluation, numComponents> origComp;
311 Evaluation sumDelta = 0.0;
312 Evaluation sumx = 0.0;
313 for (
unsigned i = 0; i < numComponents; ++i) {
314 origComp[i] = fluidState.moleFraction(phaseIdx, i);
315 relError = std::max(relError, std::abs(scalarValue(x[i])));
317 sumx += abs(fluidState.moleFraction(phaseIdx, i));
318 sumDelta += abs(x[i]);
322 const Scalar maxDelta = 0.2;
323 if (sumDelta > maxDelta)
324 x /= (sumDelta/maxDelta);
327 for (
unsigned i = 0; i < numComponents; ++i) {
328 Evaluation newx = origComp[i] - x[i];
330 if (targetFug[i] > 0)
331 newx = max(0.0, newx);
333 else if (targetFug[i] < 0)
334 newx = min(0.0, newx);
339 fluidState.setMoleFraction(phaseIdx, i, newx);
342 paramCache.updateComposition(fluidState, phaseIdx);
347 template <
class Flu
idState>
348 static Scalar calculateDefect_(
const FluidState& params,
350 const ComponentVector& targetFug)
353 for (
unsigned i = 0; i < numComponents; ++i) {
357 (targetFug[i] - params.fugacity(phaseIdx, i))
359 params.fugacityCoefficient(phaseIdx, i) );