[diffcoeff] Fix occurences of diffCoeff(p,p,c) using FS::getMainComponent

FluidSystem::getMainComponent returns the main componet index of the phase.
This doesn't actually make sense for the tracer models which why this commit also
* removes getMainComponent for tracer fluid systems (Changelog amended)
* guards all occurences of getMainComponent in code that works for tracer too
with constexpr if

CHANGELOG.md

@@ -68,6 +68,8 @@ ffSol = partial(sol, ffFaceIdx, ffCellCenterIdx);
 Regarding changes made to the effective laws and diffusionCoefficient containters, Backwards-compatibility is maintined for a large extent, barring any volumevariable classes defined externally that inherit from the non-isothermal volumevariables.
 If you use a self defined volumevariables class that inherits from the non-isothermal volumevariables, please adapt the your volumevariables class to fit to the non-isothermal volumevariables, and include the new methods for accessing the diffusion and effective diffusion coefficients.
 
+- Tracer model: tracer fluid systems do no longer provide a `getMainComponent` function since this simply doesn't make sense -- the main bulk component is not modeled.
+
 ### Deprecated properties, to be removed after 3.2:
 
 ### Deprecated classes/files, to be removed after 3.2:

dumux/flux/box/fickslaw.hh

@@ -88,8 +88,6 @@ public:
                                     const int phaseIdx,
                                     const ElementFluxVariablesCache& elemFluxVarsCache)
     {
-        ComponentFluxVector componentFlux(0.0);
-
         // get inside and outside diffusion tensors and calculate the harmonic mean
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
         const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
@@ -103,32 +101,23 @@ public:
         for (auto&& scv : scvs(fvGeometry))
             rho += massOrMolarDensity(elemVolVars[scv], referenceSystem, phaseIdx)*shapeValues[scv.indexInElement()][0];
 
+        ComponentFluxVector componentFlux(0.0);
         for (int compIdx = 0; compIdx < numComponents; compIdx++)
         {
-            if (compIdx == FluidSystem::getMainComponent(phaseIdx))
-                continue;
-
-            // effective diffusion tensors
-            using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
-            auto insideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(insideVolVars, phaseIdx, phaseIdx, compIdx);
-            auto outsideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(outsideVolVars, phaseIdx, phaseIdx, compIdx);
-
-            // scale by extrusion factor
-            insideD *= insideVolVars.extrusionFactor();
-            outsideD *= outsideVolVars.extrusionFactor();
+            if constexpr (!FluidSystem::isTracerFluidSystem())
+                if (compIdx == FluidSystem::getMainComponent(phaseIdx))
+                    continue;
 
-            // the resulting averaged diffusion tensor
-            const auto D = problem.spatialParams().harmonicMean(insideD, outsideD, scvf.unitOuterNormal());
-
-            // the mole/mass fraction gradient
-            Dune::FieldVector<Scalar, dimWorld> gradX(0.0);
-            for (auto&& scv : scvs(fvGeometry))
-                gradX.axpy(massOrMoleFraction(elemVolVars[scv], referenceSystem, phaseIdx, compIdx), fluxVarsCache.gradN(scv.indexInElement()));
+            const auto D = averageDiffusionCoefficient_(phaseIdx, compIdx, insideVolVars, outsideVolVars, problem, scvf);
 
             // compute the diffusive flux
-            componentFlux[compIdx] = -1.0*rho*vtmv(scvf.unitOuterNormal(), D, gradX)*scvf.area();
-            if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
-                componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
+            const auto massOrMoleFrac = [&](const SubControlVolume& scv){ return massOrMoleFraction(elemVolVars[scv], referenceSystem, phaseIdx, compIdx); };
+            componentFlux[compIdx] = discreteFlux_(fvGeometry, scvf, fluxVarsCache, massOrMoleFrac, D, rho);
+
+            // if the main component is balanced subtract the same flux from there (conservation)
+            if constexpr (!FluidSystem::isTracerFluidSystem())
+                if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx))
+                    componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
         }
         return componentFlux;
     }
@@ -154,20 +143,11 @@ public:
         std::array<std::vector<Scalar>, numComponents> ti;
         for (int compIdx = 0; compIdx < numComponents; compIdx++)
         {
-            if(compIdx == FluidSystem::getMainComponent(phaseIdx))
-                continue;
-
-            // effective diffusion tensors
-            using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
-            auto insideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(insideVolVars, phaseIdx, phaseIdx, compIdx);
-            auto outsideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(outsideVolVars, phaseIdx, phaseIdx, compIdx);
-
-            // scale by extrusion factor
-            insideD *= insideVolVars.extrusionFactor();
-            outsideD *= outsideVolVars.extrusionFactor();
+            if constexpr (!FluidSystem::isTracerFluidSystem())
+                if (compIdx == FluidSystem::getMainComponent(phaseIdx))
+                    continue;
 
-            // the resulting averaged diffusion tensor
-            const auto D = problem.spatialParams().harmonicMean(insideD, outsideD, scvf.unitOuterNormal());
+            const auto D = averageDiffusionCoefficient_(phaseIdx, compIdx, insideVolVars, outsideVolVars, problem, scvf);
 
             ti[compIdx].resize(fvGeometry.numScv());
             for (auto&& scv : scvs(fvGeometry))
@@ -176,6 +156,57 @@ public:
 
         return ti;
     }
+
+private:
+    static Scalar averageDiffusionCoefficient_(const int phaseIdx, const int compIdx,
+                                               const VolumeVariables& insideVV, const VolumeVariables& outsideVV,
+                                               const Problem& problem,
+                                               const SubControlVolumeFace& scvf)
+    {
+        // effective diffusion tensors
+        auto [insideD, outsideD] = diffusionCoefficientsAtInterface_(phaseIdx, compIdx, insideVV, outsideVV);
+
+        // scale by extrusion factor
+        insideD *= insideVV.extrusionFactor();
+        outsideD *= outsideVV.extrusionFactor();
+
+        // the resulting averaged diffusion tensor
+        return problem.spatialParams().harmonicMean(insideD, outsideD, scvf.unitOuterNormal());
+    }
+
+    static std::pair<Scalar, Scalar>
+    diffusionCoefficientsAtInterface_(const int phaseIdx, const int compIdx,
+                                      const VolumeVariables& insideVV, const VolumeVariables& outsideVV)
+    {
+        if constexpr (!FluidSystem::isTracerFluidSystem())
+        {
+            using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+            const auto mainCompIdx = FluidSystem::getMainComponent(phaseIdx);
+            const auto insideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(insideVV, phaseIdx, mainCompIdx, compIdx);
+            const auto outsideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(outsideVV, phaseIdx, mainCompIdx, compIdx);
+            return { std::move(insideD), std::move(outsideD) };
+        }
+        else
+        {
+            using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+            const auto insideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(insideVV, 0, 0, compIdx);
+            const auto outsideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(outsideVV, 0, 0, compIdx);
+            return { std::move(insideD), std::move(outsideD) };
+        }
+    }
+
+    template<class EvaluateVariable, class Tensor>
+    static Scalar discreteFlux_(const FVElementGeometry& fvGeometry,
+                                const SubControlVolumeFace& scvf,
+                                const FluxVarCache& fluxVarsCache,
+                                const EvaluateVariable& massOrMoleFraction,
+                                const Tensor& D, const Scalar preFactor)
+    {
+        Dune::FieldVector<Scalar, dimWorld> gradX(0.0);
+        for (auto&& scv : scvs(fvGeometry))
+            gradX.axpy(massOrMoleFraction(scv), fluxVarsCache.gradN(scv.indexInElement()));
+        return -1.0*preFactor*vtmv(scvf.unitOuterNormal(), D, gradX)*scvf.area();
+    }
 };
 
 } // end namespace Dumux

dumux/flux/ccmpfa/fickslaw.hh

@@ -190,8 +190,9 @@ public:
         ComponentFluxVector componentFlux(0.0);
         for (int compIdx = 0; compIdx < numComponents; compIdx++)
         {
-            if(compIdx == FluidSystem::getMainComponent(phaseIdx))
-                continue;
+            if constexpr (!FluidSystem::isTracerFluidSystem())
+                if (compIdx == FluidSystem::getMainComponent(phaseIdx))
+                    continue;
 
             // calculate the density at the interface
             const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
@@ -210,9 +211,10 @@ public:
         }
 
         // accumulate the phase component flux
-        for(int compIdx = 0; compIdx < numComponents; compIdx++)
-            if(compIdx != FluidSystem::getMainComponent(phaseIdx) && BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
-                componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
+        for (int compIdx = 0; compIdx < numComponents; compIdx++)
+            if constexpr (!FluidSystem::isTracerFluidSystem())
+                if (compIdx != FluidSystem::getMainComponent(phaseIdx) && BalanceEqOpts::mainComponentIsBalanced(phaseIdx))
+                    componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
 
         return componentFlux;
     }

dumux/flux/cctpfa/fickslaw.hh

@@ -141,8 +141,9 @@ public:
         ComponentFluxVector componentFlux(0.0);
         for (int compIdx = 0; compIdx < numComponents; compIdx++)
         {
-            if(compIdx == FluidSystem::getMainComponent(phaseIdx))
-                continue;
+            if constexpr (!FluidSystem::isTracerFluidSystem())
+                if (compIdx == FluidSystem::getMainComponent(phaseIdx))
+                    continue;
 
             // diffusion tensors are always solution dependent
             Scalar tij = elemFluxVarsCache[scvf].diffusionTij(phaseIdx, compIdx);
@@ -165,8 +166,9 @@ public:
                                                         : branchingFacetDensity(elemVolVars, scvf, phaseIdx, rhoInside);
 
             componentFlux[compIdx] = rho*tij*(xInside - xOutside);
-            if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
-                componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
+            if constexpr (!FluidSystem::isTracerFluidSystem())
+                if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx))
+                    componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
         }
 
         return componentFlux;
@@ -180,12 +182,21 @@ public:
                                             const SubControlVolumeFace& scvf,
                                             const int phaseIdx, const int compIdx)
     {
-        using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+
 
         const auto insideScvIdx = scvf.insideScvIdx();
         const auto& insideScv = fvGeometry.scv(insideScvIdx);
         const auto& insideVolVars = elemVolVars[insideScvIdx];
-        const auto insideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(insideVolVars, phaseIdx, phaseIdx, compIdx);
+        const auto getDiffCoeff = [&](const auto& vv)
+        {
+            using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+            if constexpr (FluidSystem::isTracerFluidSystem())
+                return Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(vv, 0, 0, compIdx);
+            else
+                return Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(vv, phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
+        };
+
+        const auto insideD = getDiffCoeff(insideVolVars);
 
         const Scalar ti = computeTpfaTransmissibility(scvf, insideScv, insideD, insideVolVars.extrusionFactor());
 
@@ -200,7 +211,7 @@ public:
             const auto outsideScvIdx = scvf.outsideScvIdx();
             const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
             const auto& outsideVolVars = elemVolVars[outsideScvIdx];
-            const auto outsideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(outsideVolVars, phaseIdx, phaseIdx, compIdx);
+            const auto outsideD = getDiffCoeff(outsideVolVars);
 
             Scalar tj;
             if (dim == dimWorld)

dumux/flux/staggered/freeflow/fickslaw.hh

@@ -149,7 +149,8 @@ private:
     static Scalar getEffectiveDiffusionCoefficient_(const VolumeVariables& volVars, const int phaseIdx, const int compIdx)
     {
         if constexpr (Dumux::Deprecated::hasEffDiffCoeff<VolumeVariables>)
-            return volVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx);
+            return volVars.effectiveDiffusionCoefficient(phaseIdx,
+                    VolumeVariables::FluidSystem::getMainComponent(phaseIdx), compIdx);
         else
         {
             // TODO: remove this else clause after release 3.2!

dumux/flux/staggered/freeflow/maxwellstefanslaw.hh

@@ -246,7 +246,8 @@ private:
     static Scalar getEffectiveDiffusionCoefficient_(const VolumeVariables& volVars, const int phaseIdx, const int compIdx)
     {
         if constexpr (Dumux::Deprecated::hasEffDiffCoeff<VolumeVariables>)
-            return volVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx);
+            return volVars.effectiveDiffusionCoefficient(phaseIdx,
+                    VolumeVariables::FluidSystem::getMainComponent(phaseIdx), compIdx);
         else
         {
             // TODO: remove this else clause after release 3.2!

dumux/freeflow/compositional/iofields.hh

@@ -79,7 +79,8 @@ struct FreeflowNCIOFields
     static double getEffectiveDiffusionCoefficient_(const VolumeVariables& volVars, const int phaseIdx, const int compIdx)
     {
         if constexpr (Dumux::Deprecated::hasEffDiffCoeff<VolumeVariables>)
-            return volVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx);
+            return volVars.effectiveDiffusionCoefficient(phaseIdx,
+                    VolumeVariables::FluidSystem::getMainComponent(phaseIdx), compIdx);
         else
         {
             // TODO: remove this else clause after release 3.2!

dumux/material/fluidsystems/base.hh

@@ -63,12 +63,12 @@ public:
      * \brief Get the main component of a given phase if possible
      *
      * \param phaseIdx The index of the fluid phase to consider
-     * \note This method has to can assert at compile time if the fluid system doesn't assume a
-     *       main phase. Then using e.g. Fick's law will fail compiling.
      * \todo Unfortunately we currently still have the assumption in some volume variables (e.g. 1pnc, 2pnc)
      *       that the main component index of a phase is equal to the phase index of that phase. This means
      *       changing this only works if the volume variables are written accordingly.
+     * \note This only makes sense if this is not a tracer fluid system (then the bulk component is not balanced)
      */
+    template<class I = Implementation, std::enable_if_t<!I::isTracerFluidSystem(), int> = 0>
     static constexpr int getMainComponent(int phaseIdx)
     { return phaseIdx; }
 

dumux/porousmediumflow/fluxvariablescachefiller.hh

@@ -156,10 +156,15 @@ private:
         static constexpr int numComponents = ModelTraits::numFluidComponents();
 
         // forward to the filler of the diffusive quantities
-        for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-            for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
-                if (compIdx != FluidSystem::getMainComponent(phaseIdx))
+        if constexpr (FluidSystem::isTracerFluidSystem())
+            for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
                     DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this);
+        else
+            for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
+                    if (compIdx != FluidSystem::getMainComponent(phaseIdx))
+                        DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element, fvGeometry, elemVolVars, scvf, *this);
     }
 
     //! method to fill the quantities related to heat conduction
@@ -447,8 +452,9 @@ private:
             for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
             {
                 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-                if (compIdx == FluidSystem::getMainComponent(phaseIdx))
-                    continue;
+                if constexpr (!FluidSystem::isTracerFluidSystem())
+                    if (compIdx == FluidSystem::getMainComponent(phaseIdx))
+                        continue;
 
                 // fill diffusion caches
                 for (unsigned int i = 0; i < iv.localFaceData().size(); ++i)
@@ -592,8 +598,9 @@ private:
             {
                 // skip main component
                 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-                if (compIdx == FluidSystem::getMainComponent(phaseIdx))
-                    continue;
+                if constexpr (!FluidSystem::isTracerFluidSystem())
+                    if (compIdx == FluidSystem::getMainComponent(phaseIdx))
+                        continue;
 
                 // fill data in the handle
                 handle.diffusionHandle().setPhaseIndex(phaseIdx);
@@ -607,13 +614,18 @@ private:
                 using IvLocalAssembler = typename Traits::template LocalAssembler<Problem, FVElementGeometry, ElementVolumeVariables>;
                 IvLocalAssembler localAssembler(problem(), fvGeometry(), elemVolVars());
 
-                // lambda to obtain diffusion coefficient
-                auto getD = [phaseIdx, compIdx] (const auto& volVars)
-                { return Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(volVars, phaseIdx, phaseIdx, compIdx); };
-
                 // maybe (re-)assemble matrices
                 if (forceUpdateAll || diffusionIsSolDependent)
                 {
+                    // lambda to obtain diffusion coefficient
+                    const auto getD = [phaseIdx, compIdx] (const auto& volVars)
+                    {
+                        if constexpr (FluidSystem::isTracerFluidSystem())
+                            return Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(volVars, 0, 0, compIdx);
+                        else
+                            return Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(volVars, phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
+                    };
+
                     // Effective diffusion coefficients might get zero if saturation = 0.
                     // Compute epsilon to detect obsolete rows in the iv-local matrices during assembly
                     if constexpr (ModelTraits::numFluidPhases() > 1)
@@ -621,9 +633,8 @@ private:
                         const auto& scv = *scvs(fvGeometry()).begin();
                         const auto& scvf = *scvfs(fvGeometry()).begin();
                         const auto& vv = elemVolVars()[scv];
-                        const auto& D = vv.diffusionCoefficient(phaseIdx, phaseIdx, compIdx);
+                        const auto& D = vv.diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
                         const auto tij = computeTpfaTransmissibility(scvf, scv, D, vv.extrusionFactor())*scvf.area();
-
                         // use transmissibility with molecular coefficient for epsilon estimate
                         localAssembler.assembleMatrices(handle.diffusionHandle(), iv, getD, tij*1e-7);
                     }

dumux/porousmediumflow/nonequilibrium/localresidual.hh

@@ -337,7 +337,7 @@ public:
                     //additionally get equilibrium values from volume variables
                     const Scalar xEquil = volVars.xEquil(phaseIdx, compIdx);
                     //get the diffusion coefficient
-                    const Scalar diffCoeff = volVars.diffusionCoefficient(phaseIdx, phaseIdx, compIdx);
+                    const Scalar diffCoeff = volVars.diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
 
                     //now compute the flux
                     const Scalar compFluxIntoOtherPhase = factorMassTransfer * (xEquil-xNonEquil)/characteristicLength * awn * volVars.molarDensity(phaseIdx) * diffCoeff * sherwoodNumber;

examples/1ptracer/problem_tracer.hh

@@ -108,10 +108,6 @@ public:
     static constexpr bool isTracerFluidSystem()
     { return true; }
 
-    // and that no component is the main component
-    static constexpr int getMainComponent(int phaseIdx)
-    { return -1; }
-
     // We define the number of components of this fluid system (one single tracer component)
     static constexpr int numComponents = 1;
 
@@ -133,8 +129,7 @@ public:
     // We set the value for the binary diffusion coefficient. This
     // might depend on spatial parameters like pressure / temperature.
     // But, in this case we neglect diffusion and return 0.0:
-    static Scalar binaryDiffusionCoefficient(unsigned int compIdxI,
-                                             unsigned int compIdxJ,
+    static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
                                              const Problem& problem,
                                              const Element& element,
                                              const SubControlVolume& scv)

test/multidomain/facet/tracer_tracer/tracerfluidsystem.hh

@@ -47,10 +47,6 @@ public:
     static constexpr bool isTracerFluidSystem()
     { return true; }
 
-    //! No component is the main component
-    static constexpr int getMainComponent(int phaseIdx)
-    { return -1; }
-
     //! The number of components
     static constexpr int numComponents = 1;
 
@@ -68,8 +64,7 @@ public:
 
     //! Binary diffusion coefficient
     //! (might depend on spatial parameters like pressure / temperature)
-    static Scalar binaryDiffusionCoefficient(unsigned int compIIdx,
-                                             unsigned int compJIdx,
+    static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
                                              const Problem& problem,
                                              const Element& element,
                                              const SubControlVolume& scv)

test/porousmediumflow/tracer/2ptracer/problem_tracer.hh

@@ -88,10 +88,6 @@ public:
     static constexpr bool isTracerFluidSystem()
     { return true; }
 
-    //! No component is the main component
-    static constexpr int getMainComponent(int phaseIdx)
-    { return -1; }
-
     //! The number of components
     static constexpr int numComponents = 1;
 

test/porousmediumflow/tracer/constvel/problem.hh

@@ -104,10 +104,6 @@ public:
     static constexpr bool isTracerFluidSystem()
     { return true; }
 
-    //! None of the components are the main component of the phase
-    static constexpr int getMainComponent(int phaseIdx)
-    { return -1; }
-
     //! The number of components
     static constexpr int numComponents = 2;
     static constexpr int numPhases = 1;