From fbb5b23e0605f37f269a3121b13e60d8020e36d5 Mon Sep 17 00:00:00 2001
From: Mathis Kelm <mathis.kelm@iws.uni-stuttgart.de>
Date: Wed, 7 Jun 2023 16:46:11 +0200
Subject: [PATCH] [freeflow] use total mass balance in old staggered

---
 .../freeflow/compositional/staggered/fluxvariables.hh | 11 +++++++++--
 .../freeflow/compositional/staggered/localresidual.hh |  4 ++--
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh
index 2767cf4cf7..2b47c53038 100644
--- a/dumux/freeflow/compositional/staggered/fluxvariables.hh
+++ b/dumux/freeflow/compositional/staggered/fluxvariables.hh
@@ -43,6 +43,7 @@ class FreeflowNCFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered>
     using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
 public:
     static constexpr auto numComponents = ModelTraits::numFluidComponents();
@@ -78,7 +79,7 @@ public:
 
             flux[compIdx] = ParentType::advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTerm);
 
-            //check for the reference system and adapt units of the diffusive flux accordingly.
+            // check for the reference system and adapt units of the diffusive flux accordingly.
             if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
             {
                 flux[compIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
@@ -89,10 +90,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 mass balance
         if (ModelTraits::replaceCompEqIdx() < numComponents)
         {
-            flux[ModelTraits::replaceCompEqIdx()] = std::accumulate(flux.begin(), flux.end(), 0.0);
+            // accumulate fluxes to a total mass based flux
+            Scalar totalMassFlux = 0.0;
+            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                totalMassFlux += useMoles ? flux[compIdx]*FluidSystem::molarMass(compIdx) : flux[compIdx];
+
+            flux[ModelTraits::replaceCompEqIdx()] = totalMassFlux;
         }
 
         return flux;
diff --git a/dumux/freeflow/compositional/staggered/localresidual.hh b/dumux/freeflow/compositional/staggered/localresidual.hh
index 93eea2edde..571d9b87a0 100644
--- a/dumux/freeflow/compositional/staggered/localresidual.hh
+++ b/dumux/freeflow/compositional/staggered/localresidual.hh
@@ -73,9 +73,9 @@ public:
                 storage[eqIdx] += s;
         }
 
-        // in case one balance is substituted by the total mass balance
+        // in case one balance is substituted by the total mass balance (use the mass density)
         if(ModelTraits::replaceCompEqIdx() < numComponents)
-            storage[ModelTraits::replaceCompEqIdx()] = density;
+            storage[ModelTraits::replaceCompEqIdx()] = volVars.density();
 
         EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
 
-- 
GitLab