From fc69d1530803e019ec1edb1c02e143ad108cd82b Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 7 Apr 2020 16:44:19 +0200
Subject: [PATCH] [diffusion] Do not store free diffusion coefficient in the
 volvars

---
 .../porousmediumflow/1pnc/volumevariables.hh  | 28 ++-------
 .../porousmediumflow/2p2c/volumevariables.hh  | 26 ++-------
 .../porousmediumflow/2pnc/volumevariables.hh  | 24 ++------
 .../porousmediumflow/3p3c/volumevariables.hh  | 14 ++---
 .../3pwateroil/volumevariables.hh             | 20 +++----
 dumux/porousmediumflow/co2/volumevariables.hh | 22 ++-----
 .../porousmediumflow/mpnc/volumevariables.hh  | 58 ++++++-------------
 .../richards/volumevariables.hh               | 25 ++------
 .../richardsnc/volumevariables.hh             | 21 ++-----
 9 files changed, 66 insertions(+), 172 deletions(-)

diff --git a/dumux/porousmediumflow/1pnc/volumevariables.hh b/dumux/porousmediumflow/1pnc/volumevariables.hh
index 2013dc8914..d23e1508c5 100644
--- a/dumux/porousmediumflow/1pnc/volumevariables.hh
+++ b/dumux/porousmediumflow/1pnc/volumevariables.hh
@@ -103,28 +103,11 @@ public:
         permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
         EnergyVolVars::updateEffectiveThermalConductivity();
 
-        // Second instance of a parameter cache.
-        // Could be avoided if diffusion coefficients also
-        // became part of the fluid state.
-        typename FluidSystem::ParameterCache paramCache;
-
-        paramCache.updatePhase(fluidState_, 0);
-
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            0,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
             return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
         };
 
-        diffCoeff_.update(getDiffusionCoefficient);
         effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
     }
 
@@ -333,13 +316,17 @@ public:
      */
     [[deprecated("Will be removed after release 3.2. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
-    { return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx); }
+    { return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx); }
 
     /*!
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -382,9 +369,6 @@ protected:
 private:
     PermeabilityType permeability_;
 
-    // Binary diffusion coefficient
-    DiffusionCoefficients diffCoeff_;
-
     // Effective diffusion coefficients for the phases
     DiffusionCoefficients effectiveDiffCoeff_;
 };
diff --git a/dumux/porousmediumflow/2p2c/volumevariables.hh b/dumux/porousmediumflow/2p2c/volumevariables.hh
index d0aa354de2..d37485441d 100644
--- a/dumux/porousmediumflow/2p2c/volumevariables.hh
+++ b/dumux/porousmediumflow/2p2c/volumevariables.hh
@@ -139,11 +139,6 @@ public:
         ParentType::update(elemSol, problem, element, scv);
         asImp_().completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
 
-        // Second instance of a parameter cache. Could be avoided if
-        // diffusion coefficients also became part of the fluid state.
-        typename FluidSystem::ParameterCache paramCache;
-        paramCache.updateAll(fluidState_);
-
         using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
 
@@ -159,21 +154,11 @@ public:
         EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
         permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
 
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            phaseIdx,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
             return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
         };
 
-        diffCoeff_.update(getDiffusionCoefficient);
         effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
 
         EnergyVolVars::updateEffectiveThermalConductivity();
@@ -390,13 +375,17 @@ public:
      */
     [[deprecated("Will be removed after release 3.2. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
-    { return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx); }
+    { return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx); }
 
     /*!
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -420,9 +409,6 @@ private:
     // Relative permeability within the control volume
     std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_;
 
-    // Binary diffusion coefficient
-    DiffusionCoefficients diffCoeff_;
-
     // Effective diffusion coefficients for the phases
     DiffusionCoefficients effectiveDiffCoeff_;
 
diff --git a/dumux/porousmediumflow/2pnc/volumevariables.hh b/dumux/porousmediumflow/2pnc/volumevariables.hh
index fc19b96e24..564571423c 100644
--- a/dumux/porousmediumflow/2pnc/volumevariables.hh
+++ b/dumux/porousmediumflow/2pnc/volumevariables.hh
@@ -137,11 +137,6 @@ public:
         /////////////
         // calculate the remaining quantities
         /////////////
-        // Second instance of a parameter cache.
-        // Could be avoided if diffusion coefficients also
-        // became part of the fluid state.
-        typename FluidSystem::ParameterCache paramCache;
-        paramCache.updateAll(fluidState_);
 
         using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         const auto& matParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
@@ -155,21 +150,11 @@ public:
         //update porosity before calculating the effective properties depending on it
         updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
 
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            phaseIdx,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
             return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
         };
 
-        diffCoeff_.update(getDiffusionCoefficient);
         effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
 
         // calculate the remaining quantities
@@ -458,7 +443,7 @@ public:
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
     {
         if (compIdx != phaseIdx)
-            return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
+            return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
         else
             DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
     }
@@ -467,7 +452,11 @@ public:
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -509,7 +498,6 @@ private:
     Scalar porosity_;               // Effective porosity within the control volume
     PermeabilityType permeability_; // Effective permeability within the control volume
     Scalar mobility_[ModelTraits::numFluidPhases()]; // Effective mobility within the control volume
-    DiffusionCoefficients diffCoeff_;
     DiffusionCoefficients effectiveDiffCoeff_;
 
 };
diff --git a/dumux/porousmediumflow/3p3c/volumevariables.hh b/dumux/porousmediumflow/3p3c/volumevariables.hh
index 1bf86058b2..3347dbd3ed 100644
--- a/dumux/porousmediumflow/3p3c/volumevariables.hh
+++ b/dumux/porousmediumflow/3p3c/volumevariables.hh
@@ -561,10 +561,6 @@ public:
          *       We can then add a normalized tensorial component
          *       e.g. obtained from DTI from the spatial params (currently not implemented)
          */
-         auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return FluidSystem::diffusionCoefficient(fluidState_, paramCache, phaseIdx, compJIdx);
-        };
 
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
@@ -574,7 +570,6 @@ public:
         // porosity & permeabilty
         updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
 
-        diffCoeff_.update(getDiffusionCoefficient);
         effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
 
         EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
@@ -717,7 +712,7 @@ public:
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
     {
         if (compIdx != phaseIdx)
-            return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
+            return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
         else
             DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
     }
@@ -726,7 +721,11 @@ public:
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::diffusionCoefficient(fluidState_, paramCache, phaseIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -749,7 +748,6 @@ private:
     Scalar mobility_[ModelTraits::numFluidPhases()];  //!< Effective mobility within the control volume
     Scalar bulkDensTimesAdsorpCoeff_; //!< the basis for calculating adsorbed NAPL
 
-    DiffusionCoefficients diffCoeff_;
     DiffusionCoefficients effectiveDiffCoeff_;
 };
 
diff --git a/dumux/porousmediumflow/3pwateroil/volumevariables.hh b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
index 1c7e703bec..7af61a9a95 100644
--- a/dumux/porousmediumflow/3pwateroil/volumevariables.hh
+++ b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
@@ -756,20 +756,11 @@ public:
         updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
         EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
 
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            if (phaseIdx != nPhaseIdx)
-                return FluidSystem::diffusionCoefficient(fluidState_, phaseIdx);
-            else
-                return 1.e-10;
-        };
-
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
             return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
         };
 
-        diffCoeff_.update(getDiffusionCoefficient);
         effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
 
         // permeability
@@ -908,13 +899,18 @@ public:
      */
     [[deprecated("Will be removed after release 3.2. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
-    { return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx); }
+    { return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx); }
 
     /*
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        if (phaseIdx != nPhaseIdx)
+            return FluidSystem::diffusionCoefficient(fluidState_, phaseIdx);
+        else
+            return 1.e-10;
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -960,9 +956,7 @@ private:
     Scalar mobility_[numPs];  //!< Effective mobility within the control volume
     Scalar bulkDensTimesAdsorpCoeff_; //!< the basis for calculating adsorbed NAPL
 
-    /* We need a tensor here !! */
     //!< Binary diffusion coefficients of the 3 components in the phases
-    DiffusionCoefficients diffCoeff_;
     DiffusionCoefficients effectiveDiffCoeff_;
 
 };
diff --git a/dumux/porousmediumflow/co2/volumevariables.hh b/dumux/porousmediumflow/co2/volumevariables.hh
index 4c6613cba7..0e0b80d7c3 100644
--- a/dumux/porousmediumflow/co2/volumevariables.hh
+++ b/dumux/porousmediumflow/co2/volumevariables.hh
@@ -148,22 +148,11 @@ public:
         EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
         permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
 
-        // update the binary diffusion and effective diffusion coefficients
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            phaseIdx,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
             return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
         };
 
-        diffCoeff_.update(getDiffusionCoefficient);
         effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
 
         EnergyVolVars::updateEffectiveThermalConductivity();
@@ -473,14 +462,18 @@ public:
         if (phaseIdx == compIdx)
             DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
         else
-            return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
+            return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
     }
 
     /*!
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -504,9 +497,6 @@ private:
     // Relative permeability within the control volume
     std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_;
 
-    // Binary diffusion coefficient
-    DiffusionCoefficients diffCoeff_;
-
     // Effective diffusion coefficients for the phases
     DiffusionCoefficients effectiveDiffCoeff_;
 };
diff --git a/dumux/porousmediumflow/mpnc/volumevariables.hh b/dumux/porousmediumflow/mpnc/volumevariables.hh
index 30eac0424d..2e994c6b12 100644
--- a/dumux/porousmediumflow/mpnc/volumevariables.hh
+++ b/dumux/porousmediumflow/mpnc/volumevariables.hh
@@ -122,23 +122,11 @@ public:
         // porosity
         updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
 
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+        if constexpr (enableDiffusion)
         {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            phaseIdx,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
-        auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
-        };
+            auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+            { return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); };
 
-        if (enableDiffusion)
-        {
-            diffCoeff_.update(getDiffusionCoefficient);
             effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
         }
 
@@ -438,7 +426,7 @@ public:
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
     {
         if (compIdx != phaseIdx)
-            return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
+            return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
         else
             DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
     }
@@ -447,7 +435,11 @@ public:
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -506,9 +498,6 @@ protected:
     FluidState fluidState_;
     SolidState solidState_;
 
-    // Binary diffusion coefficient
-    DiffusionCoefficients diffCoeff_;
-
     // Effective diffusion coefficients for the phases
     DiffusionCoefficients effectiveDiffCoeff_;
 };
@@ -585,22 +574,11 @@ public:
 
         updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
 
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+        if constexpr (enableDiffusion)
         {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            phaseIdx,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
-        auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
-        };
-        if (enableDiffusion)
-        {
-            diffCoeff_.update(getDiffusionCoefficient);
+            auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+            { return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); };
+
             effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
         }
 
@@ -955,7 +933,7 @@ public:
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
     {
         if (compIdx != phaseIdx)
-            return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
+            return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
         else
             DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
     }
@@ -964,7 +942,11 @@ public:
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -1023,8 +1005,6 @@ protected:
     FluidState fluidState_;
     SolidState solidState_;
 
-    // Binary diffusion coefficient
-    DiffusionCoefficients diffCoeff_;
     // Effective diffusion coefficients for the phases
     DiffusionCoefficients effectiveDiffCoeff_;
 };
diff --git a/dumux/porousmediumflow/richards/volumevariables.hh b/dumux/porousmediumflow/richards/volumevariables.hh
index 3ed8fc01d6..7fe5a2f6bd 100644
--- a/dumux/porousmediumflow/richards/volumevariables.hh
+++ b/dumux/porousmediumflow/richards/volumevariables.hh
@@ -115,15 +115,6 @@ public:
         minPc_ = MaterialLaw::endPointPc(materialParams);
 
         typename FluidSystem::ParameterCache paramCache;
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            phaseIdx,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
             return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
@@ -172,9 +163,6 @@ public:
 
             //binary diffusion coefficients
             paramCache.updateAll(fluidState_);
-            diffCoeff_ = getDiffusionCoefficient(FluidSystem::gasPhaseIdx,
-                                                 FluidSystem::comp1Idx,
-                                                 FluidSystem::comp0Idx);
             effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(FluidSystem::gasPhaseIdx,
                                                                    FluidSystem::comp1Idx,
                                                                    FluidSystem::comp0Idx);
@@ -201,9 +189,6 @@ public:
 
                 // binary diffusion coefficients
                 paramCache.updateAll(fluidState_);
-                diffCoeff_ = getDiffusionCoefficient(FluidSystem::gasPhaseIdx,
-                                                     FluidSystem::comp1Idx,
-                                                     FluidSystem::comp0Idx);
                 effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(FluidSystem::gasPhaseIdx,
                                                                        FluidSystem::comp1Idx,
                                                                        FluidSystem::comp0Idx);
@@ -223,7 +208,6 @@ public:
                 massFraction_[FluidSystem::gasPhaseIdx] = 0.0;
 
                 // binary diffusion coefficients (none required for liquid phase only)
-                diffCoeff_ = 0.0;
                 effectiveDiffCoeff_ = 0.0;
             }
         }
@@ -500,7 +484,7 @@ public:
     Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
     {
         assert(enableWaterDiffusionInAir() && phaseIdx == FluidSystem::gasPhaseIdx && compIdx == FluidSystem::comp0Idx);
-        return diffCoeff_;
+        return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
     }
 
     /*!
@@ -511,7 +495,9 @@ public:
         assert(enableWaterDiffusionInAir());
         assert(phaseIdx == FluidSystem::gasPhaseIdx);
         assert(compIIdx != compJIdx);
-        return diffCoeff_;
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
     }
 
     /*!
@@ -535,9 +521,6 @@ protected:
     Scalar massFraction_[numPhases]; //!< The water mass fractions in water and air
     Scalar molarDensity_[numPhases]; //!< The molar density of water and air
 
-    // Binary diffusion coefficient
-    Scalar diffCoeff_;
-
     // Effective diffusion coefficients for the phases
     Scalar effectiveDiffCoeff_;
 
diff --git a/dumux/porousmediumflow/richardsnc/volumevariables.hh b/dumux/porousmediumflow/richardsnc/volumevariables.hh
index c991d6b613..cd4dbdead4 100644
--- a/dumux/porousmediumflow/richardsnc/volumevariables.hh
+++ b/dumux/porousmediumflow/richardsnc/volumevariables.hh
@@ -113,21 +113,11 @@ public:
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updatePhase(fluidState_, 0);
 
-        auto getDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
-        {
-            return FluidSystem::binaryDiffusionCoefficient(this->fluidState_,
-                                                            paramCache,
-                                                            phaseIdx,
-                                                            compIIdx,
-                                                            compJIdx);
-        };
-
         auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
         {
             return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
         };
 
-        diffCoeff_.update(getDiffusionCoefficient);
         effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
 
         // calculate the remaining quantities
@@ -412,13 +402,17 @@ public:
      */
     [[deprecated("Will be removed after release 3.2. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
     Scalar diffusionCoefficient(const int phaseIdx, const int compIdx) const
-    { return diffCoeff_(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);}
+    { return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);}
 
     /*!
      * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
      */
     Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
-    { return diffCoeff_(phaseIdx, compIIdx, compJIdx); }
+    {
+        typename FluidSystem::ParameterCache paramCache;
+        paramCache.updatePhase(fluidState_, phaseIdx);
+        return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+    }
 
     /*!
      * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
@@ -430,9 +424,6 @@ protected:
     FluidState fluidState_; //!< the fluid state
 
 private:
-    // Binary diffusion coefficient
-    DiffusionCoefficients diffCoeff_;
-
     // Effective diffusion coefficients for the phases
     DiffusionCoefficients effectiveDiffCoeff_;
 
-- 
GitLab