From 1d871b47ac235112cbec4ed8a3aa4627434f9cb4 Mon Sep 17 00:00:00 2001
From: Mathis Kelm <mathis.kelm@iws.uni-stuttgart.de>
Date: Fri, 9 Jun 2023 17:57:40 +0200
Subject: [PATCH] [freeflow] always use total mass balance, not mole

---
 .../freeflow/navierstokes/mass/1pnc/advectiveflux.hh  |  9 +++------
 .../freeflow/navierstokes/mass/1pnc/fluxvariables.hh  | 11 +++++------
 .../freeflow/navierstokes/mass/1pnc/localresidual.hh  |  6 +++---
 3 files changed, 11 insertions(+), 15 deletions(-)

diff --git a/dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh b/dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh
index e28440bb9c..46c490d453 100644
--- a/dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh
+++ b/dumux/freeflow/navierstokes/mass/1pnc/advectiveflux.hh
@@ -56,7 +56,7 @@ struct AdvectiveFlux<NavierStokesMassOnePNCModelTraits<nComp, useM, repCompEqIdx
         static constexpr bool useMoles = ModelTraits::useMoles();
         static constexpr auto numComponents = ModelTraits::numFluidComponents();
         static constexpr auto replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
-        static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < numComponents;
+        static constexpr bool useTotalMassBalance = replaceCompEqIdx < numComponents;
 
         for (int compIdx = 0; compIdx < numComponents; ++compIdx)
         {
@@ -78,14 +78,11 @@ struct AdvectiveFlux<NavierStokesMassOnePNCModelTraits<nComp, useM, repCompEqIdx
         }
 
         // in case one balance is substituted by the total mole balance
-        if constexpr(useTotalMoleOrMassBalance)
+        if constexpr(useTotalMassBalance)
         {
             auto upwindTerm = [&]()
             {
-                if constexpr (useMoles)
-                    return [](const auto& volVars) { return volVars.molarDensity(); };
-                else
-                    return [](const auto& volVars) { return volVars.density(); };
+                return [](const auto& volVars) { return volVars.density(); };
             }();
 
             flux[replaceCompEqIdx] += upwind(upwindTerm);
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh b/dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh
index 550a7f87d6..109ed2639b 100644
--- a/dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/mass/1pnc/fluxvariables.hh
@@ -59,7 +59,7 @@ class NavierStokesMassOnePNCFluxVariables
 
     static constexpr bool enableMolecularDiffusion = ModelTraits::enableMolecularDiffusion();
     static constexpr auto replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
-    static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < ModelTraits::numFluidComponents();
+    static constexpr bool useTotalMassBalance = replaceCompEqIdx < ModelTraits::numFluidComponents();
 
     using FluidSystem = typename VolumeVariables::FluidSystem;
 
@@ -112,17 +112,16 @@ public:
                     DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
             }
 
-            // in case one balance is substituted by the total mole balance
-            if constexpr(useTotalMoleOrMassBalance)
+            // in case one balance is substituted by the total mass balance
+            if constexpr(useTotalMassBalance)
             {
                 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
                 {
                     //check for the reference system and adapt units of the diffusive flux accordingly.
                     if constexpr (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
-                        result[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
+                        result[replaceCompEqIdx] += diffusiveFluxes[compIdx];
                     else if constexpr (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
-                        result[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]
-                                                : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+                        result[replaceCompEqIdx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
                     else
                         DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
                 }
diff --git a/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh b/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
index eea71d29d8..1cfd8a4803 100644
--- a/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
+++ b/dumux/freeflow/navierstokes/mass/1pnc/localresidual.hh
@@ -64,7 +64,7 @@ class NavierStokesMassOnePNCLocalResidual : public GetPropType<TypeTag, Properti
     static constexpr auto numComponents = ModelTraits::numFluidComponents();
 
     static constexpr int replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
-    static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < numComponents;
+    static constexpr bool useTotalMassBalance = replaceCompEqIdx < numComponents;
 
 public:
     //! Use the parent type's constructor
@@ -94,8 +94,8 @@ public:
         }
 
         // in case one balance is substituted by the total mass balance
-        if constexpr (useTotalMoleOrMassBalance)
-            storage[ModelTraits::replaceCompEqIdx()] = density;
+        if constexpr (useTotalMassBalance)
+            storage[ModelTraits::replaceCompEqIdx()] = volVars.density();
 
         // consider energy storage for non-isothermal models
         if constexpr (ModelTraits::enableEnergyBalance())
-- 
GitLab