From 8eb55f7e3b5201585cc3695a4fd31a2f12d3885c Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Mon, 2 Dec 2019 09:44:38 +0100
Subject: [PATCH] [3p3c][localresidual] fix calculation of diffusive fluxes.
 Sum of fluxes is flux of the remaining phase only in the right reference
 system. That is why this has to be calculated and then the units for all
 fluxes need to be adapted

---
 dumux/porousmediumflow/3p3c/localresidual.hh | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/dumux/porousmediumflow/3p3c/localresidual.hh b/dumux/porousmediumflow/3p3c/localresidual.hh
index 335843787c..2dbf5cd49b 100644
--- a/dumux/porousmediumflow/3p3c/localresidual.hh
+++ b/dumux/porousmediumflow/3p3c/localresidual.hh
@@ -165,26 +165,28 @@ public:
         const auto diffusionFluxesWPhase = fluxVars.molecularDiffusionFlux(wPhaseIdx);
         Scalar jGW = diffusionFluxesWPhase[gCompIdx];
         Scalar jNW = diffusionFluxesWPhase[nCompIdx];
+        Scalar jWW = -(jGW+jNW);
 
         //check for the reference system and adapt units of the diffusive flux accordingly.
         if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
         {
             jGW /= FluidSystem::molarMass(gCompIdx);
             jNW /= FluidSystem::molarMass(nCompIdx);
+            jWW /= FluidSystem::molarMass(wCompIdx);
         }
-        const Scalar jWW = -(jGW+jNW);
 
         const auto diffusionFluxesGPhase = fluxVars.molecularDiffusionFlux(gPhaseIdx);
         Scalar jWG = diffusionFluxesGPhase[wCompIdx];
         Scalar jNG = diffusionFluxesGPhase[nCompIdx];
+        Scalar jGG = -(jWG+jNG);
 
         //check for the reference system and adapt units of the diffusive flux accordingly.
         if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
         {
             jWG /= FluidSystem::molarMass(wCompIdx);
             jNG /= FluidSystem::molarMass(nCompIdx);
+            jGG /= FluidSystem::molarMass(gCompIdx);
         }
-        const Scalar jGG = -(jWG+jNG);
 
         // At the moment we do not consider diffusion in the NAPL phase
         const Scalar jWN = 0.0;
-- 
GitLab