Add blackoil diffusion model

opm/models/blackoil/blackoildiffusionmodule.hh

@@ -0,0 +1,524 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief Classes required for molecular diffusion.
+ */
+#ifndef EWOMS_DIFFUSION_MODULE_HH
+#define EWOMS_DIFFUSION_MODULE_HH
+
+#include <opm/models/discretization/common/fvbaseproperties.hh>
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+
+#include <dune/common/fvector.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup Diffusion
+ * \class Opm::BlackOilDiffusionModule
+ * \brief Provides the auxiliary methods required for consideration of the
+ * diffusion equation.
+ */
+template <class TypeTag, bool enableDiffusion>
+class BlackOilDiffusionModule;
+
+/*!
+ * \copydoc Opm::BlackOilDiffusionModule
+ */
+template <class TypeTag>
+class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/false>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+
+public:
+    /*!
+     * \brief Register all run-time parameters for the diffusion module.
+     */
+    static void registerParameters()
+    {}
+
+    /*!
+     * \brief Adds the diffusive mass flux flux to the flux vector over a flux
+     *        integration point.
+      */
+    template <class Context>
+    static void addDiffusiveFlux(RateVector& flux OPM_UNUSED,
+                                 const Context& context OPM_UNUSED,
+                                 unsigned spaceIdx OPM_UNUSED,
+                                 unsigned timeIdx OPM_UNUSED)
+    {}
+};
+
+/*!
+ * \copydoc Opm::BlackOilDiffusionModule
+ */
+template <class TypeTag>
+class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/true>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+
+    using Toolbox = Opm::MathToolbox<Evaluation>;
+
+public:
+    /*!
+     * \brief Register all run-time parameters for the diffusion module.
+     */
+    static void registerParameters()
+    {}
+
+    /*!
+     * \brief Adds the mass flux due to molecular diffusion to the flux vector over the
+     *        flux integration point.
+     */
+    template <class Context>
+    static void addDiffusiveFlux(RateVector& flux, const Context& context,
+                                 unsigned spaceIdx, unsigned timeIdx)
+    {
+        const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
+
+        const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
+        const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+
+            // no diffusion in water for blackoil models
+            if (FluidSystem::waterPhaseIdx == phaseIdx) {
+                continue;
+            }
+
+            // arithmetic mean of the phase's b factor
+            Evaluation bAvg = fluidStateI.invB(phaseIdx);
+            bAvg += Toolbox::value(fluidStateJ.invB(phaseIdx));
+            bAvg /= 2;
+
+            // mass flux of solvent component (oil in oil or gas in gas)
+            unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
+            unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
+            flux[conti0EqIdx + activeSolventCompIdx] +=
+                    -bAvg
+                    * extQuants.moleFractionGradientNormal(phaseIdx, solventCompIdx)
+                    * extQuants.effectiveDiffusionCoefficient(phaseIdx, solventCompIdx);
+
+            // mass flux of solute component (gas in oil or oil in gas)
+            unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
+            unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx);
+            flux[conti0EqIdx + activeSoluteCompIdx] +=
+                    -bAvg
+                    * extQuants.moleFractionGradientNormal(phaseIdx, soluteCompIdx)
+                    * extQuants.effectiveDiffusionCoefficient(phaseIdx, soluteCompIdx);
+        }
+    }
+};
+
+/*!
+ * \ingroup Diffusion
+ * \class Opm::BlackOilDiffusionIntensiveQuantities
+ *
+ * \brief Provides the volumetric quantities required for the
+ *        calculation of molecular diffusive fluxes.
+ */
+template <class TypeTag, bool enableDiffusion>
+class BlackOilDiffusionIntensiveQuantities;
+
+/*!
+ * \copydoc Opm::DiffusionIntensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+public:
+    /*!
+     * \brief Returns the tortuousity of the sub-domain of a fluid
+     *        phase in the porous medium.
+     */
+    Scalar tortuosity(unsigned phaseIdx OPM_UNUSED) const
+    {
+        throw std::logic_error("Method tortuosity() does not make sense "
+                               "if diffusion is disabled");
+    }
+
+    /*!
+     * \brief Returns the molecular diffusion coefficient for a
+     *        component in a phase.
+     */
+    Scalar diffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
+    {
+        throw std::logic_error("Method diffusionCoefficient() does not "
+                               "make sense if diffusion is disabled");
+    }
+
+    /*!
+     * \brief Returns the effective molecular diffusion coefficient of
+     *        the porous medium for a component in a phase.
+     */
+    Scalar effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
+    {
+        throw std::logic_error("Method effectiveDiffusionCoefficient() "
+                               "does not make sense if diffusion is disabled");
+    }
+
+protected:
+    /*!
+     * \brief Update the quantities required to calculate diffusive
+     *        mass fluxes.
+     */
+    template <class FluidState>
+    void update_(FluidState& fs OPM_UNUSED,
+                 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache OPM_UNUSED,
+                 const ElementContext& elemCtx OPM_UNUSED,
+                 unsigned dofIdx OPM_UNUSED,
+                 unsigned timeIdx OPM_UNUSED)
+    { }
+};
+
+/*!
+ * \copydoc Opm::DiffusionIntensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
+
+public:
+    /*!
+     * \brief Returns the molecular diffusion coefficient for a
+     *        component in a phase.
+     */
+    Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
+    { return diffusionCoefficient_[phaseIdx][compIdx]; }
+
+    /*!
+     * \brief Returns the tortuousity of the sub-domain of a fluid
+     *        phase in the porous medium.
+     */
+    Evaluation tortuosity(unsigned phaseIdx) const
+    { return tortuosity_[phaseIdx]; }
+
+    /*!
+     * \brief Returns the effective molecular diffusion coefficient of
+     *        the porous medium for a component in a phase.
+     */
+    Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
+    { return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; }
+
+protected:
+    /*!
+     * \brief Update the quantities required to calculate diffusive
+     *        mass fluxes.
+     */
+    template <class FluidState>
+    void update_(FluidState& fluidState,
+                 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
+                 const ElementContext& elemCtx,
+                 unsigned dofIdx,
+                 unsigned timeIdx)
+    {
+        using Toolbox = Opm::MathToolbox<Evaluation>;
+
+        const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+
+            // no diffusion in water for blackoil models
+            if (FluidSystem::waterPhaseIdx == phaseIdx) {
+                continue;
+            }
+            // TODO: let the problem do this (this is a constitutive
+            // relation of which the model should be free of from the
+            // abstraction POV!)
+            const Evaluation& base =
+                Toolbox::max(0.0001,
+                             intQuants.porosity()
+                             * intQuants.fluidState().saturation(phaseIdx));
+            tortuosity_[phaseIdx] =
+                1.0 / (intQuants.porosity() * intQuants.porosity())
+                * Toolbox::pow(base, 7.0/3.0);
+
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                diffusionCoefficient_[phaseIdx][compIdx] =
+                    FluidSystem::diffusionCoefficient(fluidState,
+                                                      paramCache,
+                                                      phaseIdx,
+                                                      compIdx);
+            }
+        }
+    }
+
+private:
+    Evaluation tortuosity_[numPhases];
+    Evaluation diffusionCoefficient_[numPhases][numComponents];
+};
+
+/*!
+ * \ingroup Diffusion
+ * \class Opm::BlackOilDiffusionExtensiveQuantities
+ *
+ * \brief Provides the quantities required to calculate diffusive mass fluxes.
+ */
+template <class TypeTag, bool enableDiffusion>
+class BlackOilDiffusionExtensiveQuantities;
+
+/*!
+ * \copydoc Opm::DiffusionExtensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/false>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+
+protected:
+    /*!
+     * \brief Update the quantities required to calculate
+     *        the diffusive mass fluxes.
+     */
+    void update_(const ElementContext& elemCtx OPM_UNUSED,
+                 unsigned faceIdx OPM_UNUSED,
+                 unsigned timeIdx OPM_UNUSED)
+    {}
+
+    template <class Context, class FluidState>
+    void updateBoundary_(const Context& context OPM_UNUSED,
+                         unsigned bfIdx OPM_UNUSED,
+                         unsigned timeIdx OPM_UNUSED,
+                         const FluidState& fluidState OPM_UNUSED)
+    {}
+
+public:
+    /*!
+     * \brief The the gradient of the mole fraction times the face normal.
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     * \copydoc Doxygen::compIdxParam
+     */
+    const Evaluation& moleFractionGradientNormal(unsigned phaseIdx OPM_UNUSED,
+                                                 unsigned compIdx OPM_UNUSED) const
+    {
+        throw std::logic_error("The method moleFractionGradient() does not "
+                               "make sense if diffusion is disabled.");
+    }
+
+    /*!
+     * \brief The effective diffusion coeffcient of a component in a
+     *        fluid phase at the face's integration point
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     * \copydoc Doxygen::compIdxParam
+     */
+    const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED,
+                                                    unsigned compIdx OPM_UNUSED) const
+    {
+        throw std::logic_error("The method effectiveDiffusionCoefficient() "
+                               "does not make sense if diffusion is disabled.");
+    }
+};
+
+/*!
+ * \copydoc Opm::BlackOilDiffusionExtensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+    enum { dimWorld = GridView::dimensionworld };
+    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
+    enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
+
+    using DimVector = Dune::FieldVector<Scalar, dimWorld>;
+    using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
+
+protected:
+    /*!
+     * \brief Update the quantities required to calculate
+     *        the diffusive mass fluxes.
+     */
+    void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
+    {
+        // This uses a simple finite difference approach and is only correct for a Cartesian grids.
+        // TODO compute diffusive "transmissibilities" by re-using the logic in EclTransmissibility
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& face = stencil.interiorFace(faceIdx);
+        const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
+
+        const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
+        const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
+
+        // distance between the centers
+        const auto& insideScv = stencil.subControlVolume(face.interiorIndex());
+        const auto& outsideScv = stencil.subControlVolume(face.exteriorIndex());
+
+        DimVector distVec = outsideScv.globalPos();
+        distVec -= insideScv.globalPos();
+        Scalar dist = distVec.two_norm();
+
+        // if the following assertation triggers, the center of the
+        // center of the interior SCV was not inside the element!
+        assert(dist > 0);
+
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+            // no diffusion in water for blackoil models
+            if (FluidSystem::waterPhaseIdx == phaseIdx) {
+                continue;
+            }
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                moleFractionGradientNormal_[phaseIdx][compIdx] =
+                    (intQuantsOutside.fluidState().moleFraction(phaseIdx, compIdx)
+                     -
+                     intQuantsInside.fluidState().moleFraction(phaseIdx, compIdx))
+                    / dist;
+                Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
+
+                // use the arithmetic average for the effective
+                // diffusion coefficients.
+                effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
+                    (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
+                     +
+                     intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
+                    / 2;
+
+                Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
+            }
+        }
+    }
+
+    template <class Context, class FluidState>
+    void updateBoundary_(const Context& context,
+                         unsigned bfIdx,
+                         unsigned timeIdx,
+                         const FluidState& fluidState)
+    {
+        // This uses a simple finite difference approach and is only correct for a Cartesian grids.
+        // TODO compute diffusive "transmissibilities" by re-using the logic in EclTransmissibility
+        const auto& stencil = context.stencil(timeIdx);
+        const auto& face = stencil.boundaryFace(bfIdx);
+
+        const auto& elemCtx = context.elementContext();
+        unsigned insideScvIdx = face.interiorIndex();
+        const auto& insideScv = stencil.subControlVolume(insideScvIdx);
+
+        const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
+        const auto& fluidStateInside = intQuantsInside.fluidState();
+
+        // distance between the center of the SCV and center of the boundary face
+        DimVector distVec = face.integrationPos();
+        distVec -= context.element().geometry().global(insideScv.localGeometry().center());
+
+        Scalar dist = distVec.two_norm();
+
+        // if the following assertation triggers, the center of the
+        // center of the interior SCV was not inside the element!
+        assert(dist > 0);
+
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+            // no diffusion in water for blackoil models
+            if (FluidSystem::waterPhaseIdx == phaseIdx) {
+                continue;
+            }
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                // calculate mole fraction gradient using two-point
+                // gradients
+                moleFractionGradientNormal_[phaseIdx][compIdx] =
+                    (fluidState.moleFraction(phaseIdx, compIdx)
+                     -
+                     fluidStateInside.moleFraction(phaseIdx, compIdx))
+                    / dist;
+                Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
+
+                // use effective diffusion coefficients of the interior finite
+                // volume.
+                effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
+                    intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx);
+
+                Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
+            }
+        }
+    }
+
+public:
+    /*!
+     * \brief The the gradient of the mole fraction times the face normal.
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     * \copydoc Doxygen::compIdxParam
+     */
+    const Evaluation& moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
+    { return moleFractionGradientNormal_[phaseIdx][compIdx]; }
+
+    /*!
+     * \brief The effective diffusion coeffcient of a component in a
+     *        fluid phase at the face's integration point
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     * \copydoc Doxygen::compIdxParam
+     */
+    const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
+    { return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
+
+private:
+    Evaluation moleFractionGradientNormal_[numPhases][numComponents];
+    Evaluation effectiveDiffusionCoefficient_[numPhases][numComponents];
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilextensivequantities.hh

@@ -32,7 +32,7 @@
 #include "blackoilsolventmodules.hh"
 #include "blackoilpolymermodules.hh"
 #include "blackoilenergymodules.hh"
-
+#include "blackoildiffusionmodule.hh"
 #include <opm/models/common/multiphasebaseextensivequantities.hh>
 
 namespace Opm {
@@ -54,6 +54,7 @@ class BlackOilExtensiveQuantities
     , public BlackOilSolventExtensiveQuantities<TypeTag>
     , public BlackOilPolymerExtensiveQuantities<TypeTag>
     , public BlackOilEnergyExtensiveQuantities<TypeTag>
+    , public BlackOilDiffusionExtensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>()>
 {
     using MultiPhaseParent = MultiPhaseBaseExtensiveQuantities<TypeTag>;
 
@@ -61,6 +62,10 @@ class BlackOilExtensiveQuantities
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
+    enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
+    using DiffusionExtensiveQuantities = Opm::BlackOilDiffusionExtensiveQuantities<TypeTag, enableDiffusion>;
+
+
 public:
     /*!
      * \brief Update the extensive quantities for a given sub-control-volume-face.
@@ -77,6 +82,7 @@ public:
         asImp_().updateSolvent(elemCtx, scvfIdx, timeIdx);
         asImp_().updatePolymer(elemCtx, scvfIdx, timeIdx);
         asImp_().updateEnergy(elemCtx, scvfIdx, timeIdx);
+        DiffusionExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx);
     }
 
     template <class Context, class FluidState>
@@ -88,6 +94,8 @@ public:
         MultiPhaseParent::updateBoundary(ctx, bfIdx, timeIdx, fluidState);
 
         asImp_().updateEnergyBoundary(ctx, bfIdx, timeIdx, fluidState);
+
+        DiffusionExtensiveQuantities::updateBoundary_(ctx, bfIdx, timeIdx, fluidState);
     }
 
 protected:

opm/models/blackoil/blackoilintensivequantities.hh

@@ -35,9 +35,9 @@
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
 #include "blackoilenergymodules.hh"
+#include "blackoildiffusionmodule.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 #include <opm/material/common/Valgrind.hpp>
-
 #include <dune/common/fmatrix.hh>
 
 #include <cstring>
@@ -55,6 +55,7 @@ template <class TypeTag>
 class BlackOilIntensiveQuantities
     : public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
     , public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
+    , public BlackOilDiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
     , public BlackOilSolventIntensiveQuantities<TypeTag>
     , public BlackOilExtboIntensiveQuantities<TypeTag>
     , public BlackOilPolymerIntensiveQuantities<TypeTag>
@@ -83,6 +84,7 @@ class BlackOilIntensiveQuantities
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
     enum { waterCompIdx = FluidSystem::waterCompIdx };
@@ -101,6 +103,7 @@ class BlackOilIntensiveQuantities
     using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
     using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
     using FluidState = Opm::BlackOilFluidState<Evaluation, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled,  enableBrine, Indices::numPhases >;
+    using DiffusionIntensiveQuantities = Opm::BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>;
 
 public:
     BlackOilIntensiveQuantities()
@@ -386,6 +389,9 @@ public:
         // velocity model
         FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
 
+        // update the diffusion specific quantities of the intensive quantities
+        DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
+
 #ifndef NDEBUG
         // some safety checks in debug mode
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {

opm/models/blackoil/blackoillocalresidual.hh

@@ -35,7 +35,7 @@
 #include "blackoilenergymodules.hh"
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
-
+#include "blackoildiffusionmodule.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 
 namespace Opm {
@@ -78,6 +78,7 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
 
     static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
     static constexpr bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
+    static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
 
     using Toolbox = Opm::MathToolbox<Evaluation>;
     using SolventModule = BlackOilSolventModule<TypeTag>;
@@ -86,6 +87,7 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
     using FoamModule = BlackOilFoamModule<TypeTag>;
     using BrineModule = BlackOilBrineModule<TypeTag>;
+    using DiffusionModule = Opm::BlackOilDiffusionModule<TypeTag, enableDiffusion>;
 
 public:
     /*!
@@ -205,6 +207,8 @@ public:
 
         // deal with salt (if present)
         BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
     }
 
     /*!

opm/models/blackoil/blackoilmodel.hh

@@ -51,6 +51,8 @@
 #include <opm/models/common/multiphasebasemodel.hh>
 #include <opm/models/io/vtkcompositionmodule.hh>
 #include <opm/models/io/vtkblackoilmodule.hh>
+#include "blackoildiffusionmodule.hh"
+#include <opm/models/io/vtkdiffusionmodule.hh>
 
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 #include <opm/material/common/Unused.hpp>
@@ -72,6 +74,7 @@ namespace TTag {
 //! The type tag for the black-oil problems
 struct BlackOilModel { using InheritsFrom = std::tuple<VtkComposition,
                                                        VtkBlackOilEnergy,
+                                                       VtkDiffusion,
                                                        VtkBlackOilPolymer,
                                                        VtkBlackOilSolvent,
                                                        VtkBlackOil,
@@ -162,6 +165,10 @@ struct EnableTemperature<TypeTag, TTag::BlackOilModel> { static constexpr bool v
 template<class TypeTag>
 struct EnableEnergy<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
 
+//! disable diffusion by default
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
+
 //! by default, scale the energy equation by the inverse of the energy required to heat
 //! up one kg of water by 30 Kelvin. If we conserve surface volumes, this must be divided
 //! by the weight of one cubic meter of water. This is required to make the "dumb" linear
@@ -270,6 +277,7 @@ class BlackOilModel
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = FluidSystem::numComponents };
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+    enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
 
     static const bool compositionSwitchEnabled = Indices::gasEnabled;
     static const bool waterEnabled = Indices::waterEnabled;
@@ -278,6 +286,8 @@ class BlackOilModel
     using ExtboModule = BlackOilExtboModule<TypeTag>;
     using PolymerModule = BlackOilPolymerModule<TypeTag>;
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
+    using DiffusionModule = Opm::BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+
 public:
     BlackOilModel(Simulator& simulator)
         : ParentType(simulator)
@@ -294,10 +304,14 @@ public:
         ExtboModule::registerParameters();
         PolymerModule::registerParameters();
         EnergyModule::registerParameters();
+        DiffusionModule::registerParameters();
 
         // register runtime parameters of the VTK output modules
         Opm::VtkBlackOilModule<TypeTag>::registerParameters();
         Opm::VtkCompositionModule<TypeTag>::registerParameters();
+
+        if (enableDiffusion)
+            Opm::VtkDiffusionModule<TypeTag>::registerParameters();
     }
 
     /*!
@@ -584,6 +598,9 @@ protected:
 
         this->addOutputModule(new Opm::VtkBlackOilModule<TypeTag>(this->simulator_));
         this->addOutputModule(new Opm::VtkCompositionModule<TypeTag>(this->simulator_));
+
+        if (enableDiffusion)
+            this->addOutputModule(new Opm::VtkDiffusionModule<TypeTag>(this->simulator_));
     }
 
 private: