Merge pull request #847 from daavid00/dispersivity

Support for mechanical dispersion

opm/models/blackoil/blackoildispersionmodule.hh

@@ -0,0 +1,520 @@
+// -*- 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 mechanical dispersion.
+ */
+#ifndef EWOMS_DISPERSION_MODULE_HH
+#define EWOMS_DISPERSION_MODULE_HH
+
+#include <opm/models/discretization/common/fvbaseproperties.hh>
+
+#include <opm/material/common/Valgrind.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <stdexcept>
+
+namespace Opm {
+
+/*!
+ * \ingroup Dispersion
+ * \class Opm::BlackOilDispersionModule
+ * \brief Provides the auxiliary methods required for consideration of the
+ * dispersion equation. 
+ */
+template <class TypeTag, bool enableDispersion>
+class BlackOilDispersionModule;
+
+template <class TypeTag, bool enableDispersion>
+class BlackOilDispersionExtensiveQuantities;
+
+/*!
+ * \copydoc Opm::BlackOilDispersionModule
+ */
+template <class TypeTag>
+class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/false>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+
+    enum { numPhases = FluidSystem::numPhases };
+
+public:
+    using ExtensiveQuantities = BlackOilDispersionExtensiveQuantities<TypeTag,false>;
+    /*!
+     * \brief Adds the dispersive flux to the flux vector over a flux
+     *        integration point.
+     */
+    template <class Context>
+    static void addDispersiveFlux(RateVector&,
+                                  const Context&,
+                                  unsigned,
+                                  unsigned)
+    {}
+
+    template<class FluidState, class Scalar>
+    static void addDispersiveFlux(RateVector&,
+                                  const FluidState&,
+                                  const FluidState&,
+                                  const Evaluation&,
+                                  const Scalar&)
+    {}
+};
+
+/*!
+ * \copydoc Opm::BlackOilDispersionModule
+ */
+template <class TypeTag>
+class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/true>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Model = GetPropType<TypeTag, Properties::Model>;
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
+
+    using Toolbox = MathToolbox<Evaluation>;
+
+public:
+    using ExtensiveQuantities = BlackOilDispersionExtensiveQuantities<TypeTag,true>;
+    /*!
+     * \brief Adds the mass flux due to dispersion to the flux vector over the
+     *        flux integration point.
+     */
+    template <class Context>
+    static void addDispersiveFlux(RateVector& flux, const Context& context,
+                                  unsigned spaceIdx, unsigned timeIdx)
+    {
+        // Only work if dispersion is enabled by DISPERC in the deck
+        if (!context.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
+            return;
+        }
+        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();
+        const auto& dispersivity = extQuants.dispersivity();
+        const auto& normVelocityAvg = extQuants.normVelocityAvg();
+        addDispersiveFlux(flux, fluidStateI, fluidStateJ, dispersivity, normVelocityAvg);
+    }
+
+    /*!
+     * \brief Adds the mass flux due to dispersion to the flux vector over the
+     *        integration point. Following the notation in blackoilmodel.hh,
+     *        the dispersive flux for component \f$\kappa\f$ in phase \f$\alpha\f$
+     *        is given by: \f$-b_\alpha E||\mathrm{v}_\alpha||\mathbf{grad}X_\alpha^\kappa\f$,
+     *        where \f$b_\alpha\f$ is the shrinkage/expansion factor [-], E is the isotropic
+     *        dispersivity coefficient [L], \f$\mathrm{v}_\alpha\f$ is the filter velocity
+     *        [L/T], and \f$X_\alpha^\kappa\f$ the component mass fraction [-]. Each component mass
+     *        fraction can be computed using \f$R_s,\;R_v,\;R_{sw},\;R_{vw}\f$. For example,
+     *        \f$X_w^G=\frac{R_{sw}}{R_{sw}+\rho_w/\rho_g}\f$, where \f$\rho_w\f$ and \f$\rho_g\f$
+     *        are the reference densities.
+     *        Following the implementation of the diffusive flux (blackoildiffusionmodule.hh) and considering
+     *        the case for the water phase and gas component as an example, for cells i and j, the discrete version
+     *        of the dispersive flux at the face's integration point is given by
+     *        \f$-b_{w,ij}v_{w,ij}(\frac{1}{R_{sw,ij}+\rho_w/\rho_g})D_{ij}(R_{sw,i}-R_{sw,j})\f$
+     *        where \f$b_{w,ij}\f$, \f$v_{w,ij}\f$, and \f$R_{sw,ij}\f$ are computed using the arithmetic mean, and
+     *        the ratio \f$\frac{1}{R_{sw,ij}+\rho_w/\rho_g}\f$ is denoted as conversion factor. The dispersivity
+     *        \f$D_{ij}\f$ is computed in ecltransmissibility_impl.hh, using the dispersion coefficients \f$E_i\f$
+     *        and \f$E_j\f$.
+     */
+    template<class FluidState, class Scalar>
+    static void addDispersiveFlux(RateVector& flux,
+                                  const FluidState& fluidStateI,
+                                  const FluidState& fluidStateJ,
+                                  const Evaluation& dispersivity,
+                                  const Scalar& normVelocityAvg)
+    {
+        unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+
+            // no dispersion in water for blackoil models unless water can contain dissolved gas
+            if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
+                continue;
+            }
+
+            // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
+            if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) {
+                continue;
+            }
+
+            // arithmetic mean of the phase's b factor
+            Evaluation bAvg = fluidStateI.invB(phaseIdx);
+            bAvg += Toolbox::value(fluidStateJ.invB(phaseIdx));
+            bAvg /= 2;
+
+            Evaluation convFactor = 1.0;
+            Evaluation diffR = 0.0;
+            if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) {
+                Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2;
+                convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg);
+                diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs());
+            }
+            if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) {
+                Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2;
+                convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex));
+                diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv());
+            }
+            if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
+                Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2;
+                convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg);
+                diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw());
+            }
+            if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
+                Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2;
+                convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex));
+                diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw());
+            }
+
+            // mass flux of solvent component
+            unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
+            unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
+            flux[conti0EqIdx + activeSolventCompIdx] +=
+                    - bAvg
+                    * normVelocityAvg[phaseIdx]
+                    * convFactor
+                    * dispersivity
+                    * diffR;
+
+            // mass flux of solute component
+            unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
+            unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx);
+            flux[conti0EqIdx + activeSoluteCompIdx] +=
+                    bAvg
+                    * normVelocityAvg[phaseIdx]
+                    * convFactor
+                    * dispersivity
+                    * diffR;
+        }
+    }
+
+private:
+    static Scalar toMassFractionGasOil (unsigned regionIdx) {
+        Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
+        Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
+        return rhoO / rhoG;
+    }
+    static Scalar toMassFractionGasWater (unsigned regionIdx) {
+        Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
+        Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
+        return rhoW / rhoG;
+    }
+};
+
+/*!
+ * \ingroup Dispersion
+ * \class Opm::BlackOilDispersionIntensiveQuantities
+ *
+ * \brief Provides the volumetric quantities required for the
+ *        calculation of dispersive fluxes.
+ */
+template <class TypeTag, bool enableDispersion>
+class BlackOilDispersionIntensiveQuantities;
+
+/*!
+ * \copydoc Opm::DispersionIntensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/false>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+public:
+    /*!
+     * \brief Returns the max. norm of the filter velocity of the cell.
+     */
+    Scalar normVelocityCell(unsigned, unsigned) const
+    {
+        throw std::logic_error("Method normVelocityCell() "
+                               "does not make sense if dispersion is disabled");
+    }
+
+protected:
+    /*!
+     * \brief Update the quantities required to calculate dispersive
+     *        fluxes.
+     */
+    template<class ElementContext>
+    void update_(ElementContext&,
+                 unsigned,
+                 unsigned)
+    { }
+};
+
+/*!
+ * \copydoc Opm::DispersionIntensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/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>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+    enum { gasCompIdx = FluidSystem::gasCompIdx };
+    enum { oilCompIdx = FluidSystem::oilCompIdx };
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
+
+public:    
+    /*!
+     * \brief Returns the max. norm of the filter velocity of the cell.
+     */
+    Scalar normVelocityCell(unsigned phaseIdx) const
+    {
+
+        return normVelocityCell_[phaseIdx];
+    }
+
+protected:
+    /*!
+     * \brief Update the quantities required to calculate dispersive
+     *        mass fluxes. This considers the linear disperison model
+     *        described in the SPE CSP11 benchmark document (eq. 2.3)
+     *        https://github.com/Simulation-Benchmarks/11thSPE-CSP/
+     *        blob/main/description/spe_csp11_description.pdf
+     *        The maximum norm is used to compute the cell
+     *        filter velocity value of the corresponding phase.
+     */
+    template<class ElementContext>
+    void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    {
+        // Only work if dispersion is enabled by DISPERC in the deck
+        if (!elemCtx.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
+            return;
+        }
+        const auto& problem = elemCtx.simulator().problem();
+        if (problem.model().linearizer().getVelocityInfo().empty()) {
+            return;
+        }
+        const std::array<int, 3> phaseIdxs = { gasPhaseIdx, oilPhaseIdx, waterPhaseIdx };
+        const std::array<int, 3> compIdxs = { gasCompIdx, oilCompIdx, waterCompIdx };
+        const auto& velocityInf = problem.model().linearizer().getVelocityInfo();
+        unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
+        auto velocityInfos = velocityInf[globalDofIdx];
+        for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
+            normVelocityCell_[i] = 0;
+        }
+        for (auto& velocityInfo : velocityInfos) {
+            for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
+                if (FluidSystem::phaseIsActive(phaseIdxs[i])) {
+                    normVelocityCell_[phaseIdxs[i]] = max( normVelocityCell_[phaseIdxs[i]], 
+                        std::abs( velocityInfo.velocity[conti0EqIdx
+                            + Indices::canonicalToActiveComponentIndex(compIdxs[i])] ));
+                }
+            }
+        }
+    }
+
+private:
+    Scalar normVelocityCell_[numPhases];
+};
+
+/*!
+ * \ingroup Dispersion
+ * \class Opm::BlackOilDispersionExtensiveQuantities
+ *
+ * \brief Provides the quantities required to calculate dispersive mass fluxes.
+ */
+template <class TypeTag, bool enableDispersion>
+class BlackOilDispersionExtensiveQuantities;
+
+/*!
+ * \copydoc Opm::DispersionExtensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/false>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+
+    enum { numPhases = FluidSystem::numPhases };
+
+protected:
+    /*!
+     * \brief Update the quantities required to calculate
+     *        the dispersive fluxes.
+     */
+    void update_(const ElementContext&,
+                 unsigned,
+                 unsigned)
+    {}
+
+    template <class Context, class FluidState>
+    void updateBoundary_(const Context&,
+                         unsigned,
+                         unsigned,
+                         const FluidState&)
+    {}
+
+public:
+    using ScalarArray = Scalar[numPhases];
+
+    static void update(ScalarArray&,
+                       const IntensiveQuantities&,
+                       const IntensiveQuantities&)
+    {}
+
+    /*!
+     * \brief The dispersivity the face.
+     *
+     */
+    const Scalar& dispersivity() const
+    {
+        throw std::logic_error("The method dispersivity() does not "
+                               "make sense if dispersion is disabled.");
+    }
+
+    /*!
+     * \brief The effective filter velocity coefficient in a
+     *        fluid phase at the face's integration point
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     * \copydoc Doxygen::compIdxParam
+     */
+    const Scalar& normVelocityAvg(unsigned) const
+    {
+        throw std::logic_error("The method normVelocityAvg() "
+                               "does not make sense if dispersion is disabled.");
+    }
+
+};
+
+/*!
+ * \copydoc Opm::BlackOilDispersionExtensiveQuantities
+ */
+template <class TypeTag>
+class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/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>;
+    using Toolbox = MathToolbox<Evaluation>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+
+    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>;
+public:
+    using ScalarArray = Scalar[numPhases];
+    static void update(ScalarArray& normVelocityAvg,
+                       const IntensiveQuantities& intQuantsInside,
+                       const IntensiveQuantities& intQuantsOutside)
+    {
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+            // no dispersion in water for blackoil models unless water can contain dissolved gas
+            if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
+                continue;
+            }
+            // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
+            if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) {
+                continue;
+            }
+            // use the arithmetic average for the effective
+            // velocity coefficients at the face's integration point.
+            normVelocityAvg[phaseIdx] = 0.5 *
+                ( intQuantsInside.normVelocityCell(phaseIdx) +
+                    intQuantsOutside.normVelocityCell(phaseIdx) );
+            Valgrind::CheckDefined(normVelocityAvg[phaseIdx]);
+        }
+    }
+protected:
+    template <class Context, class FluidState>
+    void updateBoundary_(const Context&,
+                         unsigned,
+                         unsigned,
+                         const FluidState&)
+    {
+        throw std::runtime_error("Not implemented: Dispersion across boundary not implemented for blackoil");
+    }
+
+public:
+    /*!
+     * \brief The dispersivity of the face.
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     * \copydoc Doxygen::compIdxParam
+     */
+    const Scalar& dispersivity() const
+    { return dispersivity_; }
+
+    /*!
+     * \brief The effective velocity coefficient in a
+     *        fluid phase at the face's integration point
+     *
+     * \copydoc Doxygen::phaseIdxParam
+     * \copydoc Doxygen::compIdxParam
+     */
+    const Scalar& normVelocityAvg(unsigned phaseIdx) const
+    { return normVelocityAvg_[phaseIdx]; }
+
+    const auto& normVelocityAvg() const{
+        return normVelocityAvg_;
+    }
+
+private:
+    Scalar dispersivity_;
+    ScalarArray normVelocityAvg_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilintensivequantities.hh

@@ -36,6 +36,7 @@
 #include "blackoilbrinemodules.hh"
 #include "blackoilenergymodules.hh"
 #include "blackoildiffusionmodule.hh"
+#include "blackoildispersionmodule.hh"
 #include "blackoilmicpmodules.hh"
 
 #include <opm/common/TimingMacros.hpp>
@@ -70,6 +71,7 @@ class BlackOilIntensiveQuantities
     : public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
     , public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
     , public BlackOilDiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
+    , public BlackOilDispersionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDispersion>() >
     , public BlackOilSolventIntensiveQuantities<TypeTag>
     , public BlackOilExtboIntensiveQuantities<TypeTag>
     , public BlackOilPolymerIntensiveQuantities<TypeTag>
@@ -103,6 +105,7 @@ class BlackOilIntensiveQuantities
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
+    enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
@@ -124,6 +127,7 @@ class BlackOilIntensiveQuantities
     using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
     using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
     using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>;
+    using DispersionIntensiveQuantities = BlackOilDispersionIntensiveQuantities<TypeTag, enableDispersion>;
 
     using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
 
@@ -473,6 +477,9 @@ public:
         // update the diffusion specific quantities of the intensive quantities
         DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
 
+        // update the dispersion specific quantities of the intensive quantities
+        DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
+
 #ifndef NDEBUG
         // some safety checks in debug mode
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -36,6 +36,7 @@
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
 #include "blackoildiffusionmodule.hh"
+#include "blackoildispersionmodule.hh"
 #include "blackoilmicpmodules.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
@@ -92,6 +93,7 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
     static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
     static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
+    static constexpr bool enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>();
     static constexpr bool enableMICP = getPropValue<TypeTag, Properties::EnableMICP>();
 
     using SolventModule = BlackOilSolventModule<TypeTag>;
@@ -101,6 +103,7 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     using FoamModule = BlackOilFoamModule<TypeTag>;
     using BrineModule = BlackOilBrineModule<TypeTag>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+    using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
     using MICPModule = BlackOilMICPModule<TypeTag>;
 
     using Toolbox = MathToolbox<Evaluation>;
@@ -118,7 +121,8 @@ public:
         double Vex;
         double inAlpha;
         double outAlpha;
-	double diffusivity;
+        double diffusivity;
+        double dispersivity;
     };
     /*!
      * \copydoc FvBaseLocalResidual::computeStorage
@@ -295,9 +299,10 @@ public:
         // for thermal harmonic mean of half trans
         const Scalar inAlpha = problem.thermalHalfTransmissibility(globalIndexIn, globalIndexEx);
         const Scalar outAlpha = problem.thermalHalfTransmissibility(globalIndexEx, globalIndexIn);
-	const Scalar diffusivity = problem.diffusivity(globalIndexEx, globalIndexIn);
+        const Scalar diffusivity = problem.diffusivity(globalIndexEx, globalIndexIn);
+        const Scalar dispersivity = problem.dispersivity(globalIndexEx, globalIndexIn);
 
-        const ResidualNBInfo res_nbinfo {trans, faceArea, thpres, distZ * g, facedir, Vin, Vex, inAlpha, outAlpha, diffusivity};
+        const ResidualNBInfo res_nbinfo {trans, faceArea, thpres, distZ * g, facedir, Vin, Vex, inAlpha, outAlpha, diffusivity, dispersivity};
 
         calculateFluxes_(flux,
                          darcy,
@@ -441,7 +446,7 @@ public:
         static_assert(!enableBrine, "Relevant computeFlux() method must be implemented for this module before enabling.");
         // BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
-        // deal with diffusion (if present)
+        // deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity).
         if constexpr(enableDiffusion){
             typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
             DiffusionModule::ExtensiveQuantities::update(effectiveDiffusionCoefficient, intQuantsIn, intQuantsEx);
@@ -453,6 +458,19 @@ public:
                                               tmpdiffusivity,
                                               effectiveDiffusionCoefficient);
 
+        }
+        // deal with dispersion (if present). opm-models expects per area flux (added in the tmpdispersivity).
+        if constexpr(enableDispersion){
+            typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg;
+            DispersionModule::ExtensiveQuantities::update(normVelocityAvg, intQuantsIn, intQuantsEx);
+            const Scalar dispersivity = nbInfo.dispersivity;
+            const Scalar tmpdispersivity = dispersivity / faceArea;
+            DispersionModule::addDispersiveFlux(flux,
+                                                intQuantsIn.fluidState(),
+                                                intQuantsEx.fluidState(),
+                                                tmpdispersivity,
+                                                normVelocityAvg);
+
         }
         // deal with micp (if present)
         static_assert(!enableMICP, "Relevant computeFlux() method must be implemented for this module before enabling.");

opm/models/blackoil/blackoilmodel.hh

@@ -53,6 +53,7 @@
 #include <opm/models/io/vtkcompositionmodule.hh>
 #include <opm/models/io/vtkblackoilmodule.hh>
 #include "blackoildiffusionmodule.hh"
+#include "blackoildispersionmodule.hh"
 #include <opm/models/io/vtkdiffusionmodule.hh>
 
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
@@ -177,6 +178,10 @@ struct EnableEnergy<TypeTag, TTag::BlackOilModel> { static constexpr bool value
 template<class TypeTag>
 struct EnableDiffusion<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
 
+//! disable disperison by default
+template<class TypeTag>
+struct EnableDispersion<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
@@ -227,7 +232,7 @@ namespace Opm {
  * \f]
  *
  * Since the gas and water phases are assumed to be immiscible, this
- * is sufficint to calculate their density. For the formation volume
+ * is sufficient to calculate their density. For the formation volume
  * factor of the the oil phase \f$B_o\f$ determines the density of
  * *saturated* oil, i.e. the density of the oil phase if some gas
  * phase is present.
@@ -288,6 +293,7 @@ private:
     enum { numComponents = FluidSystem::numComponents };
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
+    enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
 
     static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
     static constexpr bool waterEnabled = Indices::waterEnabled;
@@ -297,6 +303,7 @@ private:
     using PolymerModule = BlackOilPolymerModule<TypeTag>;
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+    using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
     using MICPModule = BlackOilMICPModule<TypeTag>;
 
 public:

opm/models/common/multiphasebaseproperties.hh

@@ -80,6 +80,9 @@ struct EnableGravity { using type = UndefinedProperty; };
 //! Enable diffusive fluxes?
 template<class TypeTag, class MyTypeTag>
 struct EnableDiffusion { using type = UndefinedProperty; };
+//! Enable dispersive fluxes?
+template<class TypeTag, class MyTypeTag>
+struct EnableDispersion { using type = UndefinedProperty; };
 
 } // namespace Opm::Properties
 

opm/models/discretization/common/tpfalinearizer.hh

@@ -329,6 +329,16 @@ public:
         return floresInfo_;
     }
 
+    /*!
+     * \brief Return constant reference to the velocityInfo.
+     *
+     * (This object is only non-empty if the DISPERC keyword is true.)
+     */
+    const auto& getVelocityInfo() const{
+
+        return velocityInfo_;
+    }
+
     void updateDiscretizationParameters()
     {
         updateStoredTransmissibilities();
@@ -451,6 +461,7 @@ private:
                         Scalar outAlpha {0.};
                         FaceDirection dirId = FaceDirection::Unknown;
                         Scalar diffusivity {0.};
+                        Scalar dispersivity {0.};
                         if constexpr(enableEnergy){
                             inAlpha = problem_().thermalHalfTransmissibility(myIdx, neighborIdx);
                             outAlpha = problem_().thermalHalfTransmissibility(neighborIdx, myIdx);
@@ -458,10 +469,13 @@ private:
                         if constexpr(enableDiffusion){
                             diffusivity = problem_().diffusivity(myIdx, neighborIdx);
                         }
+                        if (simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
+                            dispersivity = problem_().dispersivity(myIdx, neighborIdx);
+                        }
                         if (materialLawManager->hasDirectionalRelperms()) {
                             dirId = scvf.faceDirFromDirId();
                         }
-                        loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, {trans, area, thpres, dZg, dirId, Vin, Vex, inAlpha, outAlpha, diffusivity}, nullptr};
+                        loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, {trans, area, thpres, dZg, dirId, Vin, Vex, inAlpha, outAlpha, diffusivity, dispersivity}, nullptr};
 
                     }
                 }
@@ -522,15 +536,17 @@ private:
         jacobian_->clear();
     }
 
-    // Initialize the flows and flores sparse tables
+    // Initialize the flows, flores, and velocity sparse tables
     void createFlows_()
     {
         OPM_TIMEBLOCK(createFlows);
         // If FLOWS/FLORES is set in any RPTRST in the schedule, then we initializate the sparse tables
         // For now, do the same also if any block flows are requested (TODO: only save requested cells...)
+        // If DISPERC is in the deck, we initialize the sparse table here as well.
         const bool anyFlows = simulator_().problem().eclWriter()->eclOutputModule().anyFlows();
-        const bool anyFlores = simulator_().problem().eclWriter()->eclOutputModule().anyFlores();
+        const bool anyFlores = simulator_().problem().eclWriter()->eclOutputModule().anyFlores();                         
-        if ((!anyFlows || !flowsInfo_.empty())  && (!anyFlores || !floresInfo_.empty())) {
+        const bool enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
+        if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && !enableDispersion) {
             return;
         }
         const auto& model = model_();
@@ -539,6 +555,7 @@ private:
         unsigned numCells = model.numTotalDof();
         std::unordered_multimap<int, std::pair<int, int>> nncIndices;
         std::vector<FlowInfo> loc_flinfo;
+        std::vector<VelocityInfo> loc_vlinfo;
         unsigned int nncId = 0;
         VectorBlock flow(0.0);
 
@@ -555,12 +572,16 @@ private:
         if (anyFlores) {
             floresInfo_.reserve(numCells, 6 * numCells);
         }
+        if (enableDispersion) {
+            velocityInfo_.reserve(numCells, 6 * numCells);
+        }
 
         for (const auto& elem : elements(gridView_())) {
             stencil.update(elem);
             for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
                 unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
                 loc_flinfo.resize(stencil.numDof() - 1);
+                loc_vlinfo.resize(stencil.numDof() - 1);
                 for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
                     unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
                     if (dofIdx > 0) {
@@ -579,6 +600,7 @@ private:
                             }
                         }
                         loc_flinfo[dofIdx - 1] = FlowInfo{faceId, flow, nncId};
+                        loc_vlinfo[dofIdx - 1] = VelocityInfo{flow};
                     }
                 }
                 if (anyFlows) {
@@ -587,6 +609,9 @@ private:
                 if (anyFlores) {
                     floresInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end());
                 }
+                if (enableDispersion) {
+                    velocityInfo_.appendRow(loc_vlinfo.begin(), loc_vlinfo.end());
+                }
             }
         }
     }
@@ -626,7 +651,7 @@ private:
         // We do not call resetSystem_() here, since that will set
         // the full system to zero, not just our part.
         // Instead, that must be called before starting the linearization.
-
+        const bool& enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
         const bool& enableFlows = simulator_().problem().eclWriter()->eclOutputModule().hasFlows() ||
                                     simulator_().problem().eclWriter()->eclOutputModule().hasBlockFlows();
         const bool& enableFlores = simulator_().problem().eclWriter()->eclOutputModule().hasFlores();
@@ -661,6 +686,11 @@ private:
                 const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
                 LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo);
                 adres *= nbInfo.res_nbinfo.faceArea;
+                if (enableDispersion) {
+                    for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
+                        velocityInfo_[globI][loc].velocity[phaseIdx] = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
+                    }
+                }
                 if (enableFlows) {
                     for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
                         flowsInfo_[globI][loc].flow[phaseIdx] = adres[phaseIdx].value();
@@ -823,6 +853,12 @@ private:
     SparseTable<FlowInfo> flowsInfo_;
     SparseTable<FlowInfo> floresInfo_;
 
+    struct VelocityInfo
+    {
+        VectorBlock velocity;
+    };
+    SparseTable<VelocityInfo> velocityInfo_;
+
     using ScalarFluidState = typename IntensiveQuantities::ScalarFluidState;
     struct BoundaryConditionData
     {