From 85fdc400eb7b7e5f56e0eb8e768aea053aa668ca Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Wed, 25 Jul 2018 11:09:42 +0200
Subject: [PATCH] [cleanup][freeflow] calculate molar flux to mass flux in case
 of useMoles false in fluxvariables instead of in ficks law. diffusion laws
 should always return the molar flux

---
 .../staggered/freeflow/fickslaw.hh               | 16 ++++------------
 .../compositional/staggered/fluxvariables.hh     |  7 +++++--
 2 files changed, 9 insertions(+), 14 deletions(-)

diff --git a/dumux/discretization/staggered/freeflow/fickslaw.hh b/dumux/discretization/staggered/freeflow/fickslaw.hh
index 102ce6d88a..f56fde0bd6 100644
--- a/dumux/discretization/staggered/freeflow/fickslaw.hh
+++ b/dumux/discretization/staggered/freeflow/fickslaw.hh
@@ -19,7 +19,7 @@
 /*!
  * \file
  * \brief This file contains the data which is required to calculate
- *        diffusive mass fluxes due to molecular diffusion with Fick's law.
+ *        diffusive molar fluxes due to molecular diffusion with Fick's law.
  */
 #ifndef DUMUX_DISCRETIZATION_STAGGERED_FICKS_LAW_HH
 #define DUMUX_DISCRETIZATION_STAGGERED_FICKS_LAW_HH
@@ -118,19 +118,11 @@ public:
             flux[compIdx] = avgDensity * tij * (insideMoleFraction - outsideMoleFraction);
         }
 
-        const Scalar cumulativeFlux = std::accumulate(flux.begin(), flux.end(), 0.0);
-        flux[FluidSystem::getMainComponent(0)] = - cumulativeFlux;
-
         // Fick's law (for binary systems) states that the net flux of moles within the bulk phase has to be zero:
         // If a given amount of molecules A travel into one direction, the same amount of molecules B have to
-        // go into the opposite direction. This, however, means that the net mass flux whithin the bulk phase
-        // (given different molecular weigths of A and B) does not have to be zero. We account for this:
-        if(!useMoles)
-        {
-            //convert everything to a mass flux
-            for(int compIdx = 0; compIdx < numComponents; ++compIdx)
-                flux[compIdx] *= FluidSystem::molarMass(compIdx);
-        }
+        // go into the opposite direction.
+        const Scalar cumulativeFlux = std::accumulate(flux.begin(), flux.end(), 0.0);
+        flux[FluidSystem::getMainComponent(0)] = - cumulativeFlux;
 
         return flux;
     }
diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh
index 4bfaea1854..9a0d891210 100644
--- a/dumux/freeflow/compositional/staggered/fluxvariables.hh
+++ b/dumux/freeflow/compositional/staggered/fluxvariables.hh
@@ -52,6 +52,7 @@ class FreeflowNCFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
     using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
 public:
@@ -73,6 +74,8 @@ public:
     {
         CellCenterPrimaryVariables flux(0.0);
 
+        const auto diffusiveFluxes = MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf);
+
         for (int compIdx = 0; compIdx < numComponents; ++compIdx)
         {
             auto upwindTerm = [compIdx](const auto& volVars)
@@ -83,9 +86,9 @@ public:
             };
 
             flux[compIdx] = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTerm);
-        }
 
-        flux += MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf);
+            flux[compIdx] += useMoles ? diffusiveFluxes[compIdx] : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+        }
 
         // in case one balance is substituted by the total mass balance
         if (ModelTraits::replaceCompEqIdx() < numComponents)
-- 
GitLab