From c45c329620fe9daf880169d07de783528bd30928 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Wed, 4 Jul 2018 11:29:54 +0200
Subject: [PATCH] [flash] Introduce simple preconditioning to make ncp flash
more robust
---
dumux/material/constraintsolvers/ncpflash.hh | 26 +++++++++++++++++---
1 file changed, 23 insertions(+), 3 deletions(-)
diff --git a/dumux/material/constraintsolvers/ncpflash.hh b/dumux/material/constraintsolvers/ncpflash.hh
index f1308db95e..359350c7f4 100644
--- a/dumux/material/constraintsolvers/ncpflash.hh
+++ b/dumux/material/constraintsolvers/ncpflash.hh
@@ -182,10 +182,9 @@ public:
const int nMax = 50; // <- maximum number of newton iterations
for (int nIdx = 0; nIdx < nMax; ++nIdx) {
// regularize fluid state
- for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
if (fluidState.averageMolarMass(phaseIdx) < 1e-10)
fluidState.setMoleFraction(phaseIdx, /*compIdx=*/0, 1e-10);
- }
// calculate Jacobian matrix and right hand side
linearize_<MaterialLaw>(J, b, fluidState, paramCache, matParams, globalMolarities);
@@ -195,7 +194,23 @@ public:
// Solve J*x = b
deltaX = 0;
- try { J.solve(deltaX, b); }
+ try {
+
+ // simple preconditioning of the linear system to reduce condition number
+ int eqIdx = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
+ {
+ J[eqIdx] *= 10e-7; // roughly the magnitude of the fugacities
+ b[eqIdx] *= 10e-7;
+ ++eqIdx;
+ }
+ }
+
+ J.solve(deltaX, b);
+
+ }
catch (Dune::FMatrixError &e)
{
/*
@@ -271,6 +286,11 @@ protected:
std::cout << fs.density(phaseIdx) << " ";
std::cout << "\n";
+ std::cout << "molar densities: ";
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ std::cout << fs.molarDensity(phaseIdx) << " ";
+ std::cout << "\n";
+
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
--
GitLab