Merge branch 'fix/old-staggered-fluxes' into 'master'

Fix/old staggered fluxes

See merge request !3572

CHANGELOG.md

@@ -20,6 +20,8 @@ be evaluated with the function values provided in the same order as the names wh
 - __Hyperelastic__: Added a hyperelastic model (large deformations) and a test (in geomechanics)
 - __Property system__: Added a function `inheritsFrom` that allows to query whether a type tag is present in the inheritance hierarchy of another type tag
 - __PDESolver__: The LinearPDESolver/NewtonSolver now has an `apply` method that returns a bool (true if converged) instead of throwing an exception
+- __Freeflow Nonisothermal__: An issue with the compositional heat flux's upwinding calculations has been fixed
+- __Compositional Staggered__: The staggered discretization method has been updated to fix inconsistencies in handling density when evaluating diffusive boundary fluxes. In addition, a total mass balance is now always considered, rather than a total mole balance.
 
 ### Immediate interface changes not allowing/requiring a deprecation period:
 

dumux/flux/staggered/freeflow/fickslaw.hh

@@ -115,7 +115,10 @@ public:
 
             if (scvf.boundary())
             {
-                flux[compIdx] = insideDensity * insideDiffCoeff
+                // When the face lies on a boundary, use the average density
+                const Scalar outsideDensity = massOrMolarDensity(outsideVolVars, referenceSystem, phaseIdx);
+                const Scalar avgDensity = 0.5*(insideDensity + outsideDensity);
+                flux[compIdx] = avgDensity * insideDiffCoeff
                                 * (massOrMoleFractionInside - massOrMoleFractionOutside) / insideDistance;
             }
             else

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;

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);
 

dumux/freeflow/nonisothermal/localresidual.hh

@@ -145,14 +145,15 @@ public:
         auto diffusiveFlux = FluxVariables::MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf);
         for (int compIdx = 0; compIdx < FluxVariables::numComponents; ++compIdx)
         {
-            const bool insideIsUpstream = scvf.directionSign() == sign(diffusiveFlux[compIdx]);
+            // define the upstream direction according to the sign of the diffusive flux
+            using std::signbit;
+            const bool insideIsUpstream = !signbit(diffusiveFlux[compIdx]);
             const auto& upstreamVolVars = insideIsUpstream ? elemVolVars[scvf.insideScvIdx()] : elemVolVars[scvf.outsideScvIdx()];
 
             if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
                 flux[localEnergyBalanceIdx] += diffusiveFlux[compIdx] * upstreamVolVars.componentEnthalpy(compIdx);
             else
                 flux[localEnergyBalanceIdx] += diffusiveFlux[compIdx] * upstreamVolVars.componentEnthalpy(compIdx)* elemVolVars[scvf.insideScvIdx()].molarMass(compIdx);
-
         }
     }
 };

test/freeflow/ransnc/CMakeLists.txt

@@ -15,7 +15,7 @@ dumux_add_test(NAME test_ff_rans2c_kepsilon_twolayer
                             COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_rans2c_kepsilon_twolayer-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_kepsilon_twolayer-00032.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_kepsilon_twolayer-00031.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_kepsilon params.input
                              -Problem.Name test_ff_rans2c_kepsilon_twolayer")
 
@@ -38,7 +38,7 @@ dumux_add_test(NAME test_ff_rans2c_sst
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS  --script fuzzy
                         --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_rans2c_sst-reference.vtu
-                                ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_sst-00035.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_sst-00034.vtu
                         --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_sst params.input
                         -Problem.Name test_ff_rans2c_sst
                         -RANS.WriteFlatWallBoundedFields False
@@ -52,7 +52,7 @@ dumux_add_test(NAME test_ff_rans2c_komega
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_rans2c_komega-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_komega-00036.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_komega-00035.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_komega params.input
                              -Problem.Name test_ff_rans2c_komega")
 
@@ -89,7 +89,7 @@ dumux_add_test(NAME test_ff_rans2c_zeroeq
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_rans2c_zeroeq-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_zeroeq-00030.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_zeroeq-00029.vtu
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2c_zeroeq params.input
                              -Problem.Name test_ff_rans2c_zeroeq")
 
@@ -133,7 +133,7 @@ dumux_add_test(NAME test_ff_rans2cni_sst
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS  --script fuzzy
                         --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_rans2cni_sst-reference.vtu
-                                ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2cni_sst-00031.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2cni_sst-00030.vtu
                         --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_rans2cni_sst params_nonisothermal.input
                         -Problem.Name test_ff_rans2cni_sst
                         -RANS.WriteFlatWallBoundedFields False

test/multidomain/boundary/stokesdarcy/1p2c_2p2c/CMakeLists.txt

@@ -25,9 +25,9 @@ dumux_add_test(NAME test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal_stokes-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal_stokes-00039.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal_stokes-00040.vtu
                                      ${CMAKE_SOURCE_DIR}/test/references/test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal_darcy-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal_darcy-00039.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal_darcy-00040.vtu
 
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal params_nonisothermal.input
                              -Vtk.OutputName test_md_boundary_darcy2p2cni_stokes1p2cni_horizontal")