[mpfa][fluxcaching] use constant zero threshold for diffcoeffs

dumux/porousmediumflow/fluxvariablescachefiller.hh

@@ -194,6 +194,7 @@ class PorousMediumFluxVariablesCacheFillerImplementation<TypeTag, Discretization
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
     using Element = typename GridView::template Codim<0>::Entity;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
     using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
     using FVElementGeometry = typename GridGeometry::LocalView;
@@ -624,22 +625,23 @@ private:
                             return volVars.effectiveDiffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
                     };
 
-                    // Effective diffusion coefficients might get zero if saturation = 0.
+                    // Effective diffusion coefficients might be zero if saturation = 0.
                     // Compute epsilon to detect obsolete rows in the iv-local matrices during assembly
-                    const auto& scv = *scvs(fvGeometry()).begin();
-                    const auto& scvf = *scvfs(fvGeometry()).begin();
+                    static const auto zeroD = getParamFromGroup<Scalar>(
+                        problem().paramGroup(),
+                        "Mpfa.ZeroEffectiveDiffusionCoefficientThreshold",
+                        1e-16
+                    );
+
+                    // compute a representative transmissibility for this interaction volume, using
+                    // the threshold for zero diffusion coefficients, and use this as epsilon
+                    const auto& scv = fvGeometry().scv(iv.localScv(0).gridScvIndex());
+                    const auto& scvf = fvGeometry().scvf(iv.localScvf(0).gridScvfIndex());
                     const auto& vv = elemVolVars()[scv];
-                    const auto D = [&] ()
-                    {
-                        // diffusion coefficients below 1e-20 are treated as zeroes!!
-                        using std::max;
-                        if constexpr (!FluidSystem::isTracerFluidSystem())
-                            return max(1e-20, vv.diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx));
-                        else
-                            return max(1e-20, vv.diffusionCoefficient(0, 0, compIdx));
-                    } ();
+                    const auto eps = Extrusion::area(scvf)*computeTpfaTransmissibility(
+                        scvf, scv, zeroD, vv.extrusionFactor()
+                    );
 
-                    auto eps = 1e-7*computeTpfaTransmissibility(scvf, scv, D, vv.extrusionFactor())*Extrusion::area(scvf);
                     localAssembler.assembleMatrices(handle.diffusionHandle(), iv, getD, eps);
                 }