Support for mechanical dispersion

ebos/eclalugridvanguard.hh

@@ -189,7 +189,9 @@ public:
                                                                  getPropValue<TypeTag,
                                                                  Properties::EnableEnergy>(),
                                                                  getPropValue<TypeTag,
-                                                                 Properties::EnableDiffusion>());
+                                                                 Properties::EnableDiffusion>(),
+                                                                 getPropValue<TypeTag,
+                                                                 Properties::EnableDispersion>());
             // Re-ordering  for ALUGrid
             globalTrans_->update(false, [&](unsigned int i) { return gridEquilIdxToGridIdx(i);});
         }

ebos/eclcpgridvanguard.hh

@@ -281,7 +281,8 @@ protected:
                                                     this->grid(),
                                                     this->cellCentroids(),
                                                     getPropValue<TypeTag, Properties::EnableEnergy>(),
-                                                    getPropValue<TypeTag, Properties::EnableDiffusion>()));
+                                                    getPropValue<TypeTag, Properties::EnableDiffusion>(),
+                                                    getPropValue<TypeTag, Properties::EnableDispersion>()));
         globalTrans_->update(false);
     }
 

ebos/eclproblem.hh

@@ -145,6 +145,7 @@ class EclProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     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 { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
     enum { enableApiTracking = getPropValue<TypeTag, Properties::EnableApiTracking>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
@@ -278,7 +279,8 @@ public:
                               simulator.vanguard().grid(),
                               simulator.vanguard().cellCentroids(),
                               enableEnergy,
-                              enableDiffusion)
+                              enableDiffusion,
+                              enableDispersion)
         , thresholdPressures_(simulator)
         , wellModel_(simulator)
         , aquiferModel_(simulator)
@@ -824,6 +826,13 @@ public:
         return transmissibilities_.diffusivity(globalCellIn, globalCellOut);
     }
 
+    /*!
+     * give the dispersivity for a face i.e. pair.
+     */
+    Scalar dispersivity(const unsigned globalCellIn, const unsigned globalCellOut) const{
+        return transmissibilities_.dispersivity(globalCellIn, globalCellOut);
+    }
+
     /*!
      * \brief Direct access to a boundary transmissibility.
      */
@@ -2448,6 +2457,7 @@ private:
         ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;
         ConditionalStorage<enableEnergy, Scalar> thermalHalfTransOut;
         ConditionalStorage<enableDiffusion, Scalar> diffusivity;
+        ConditionalStorage<enableDispersion, Scalar> dispersivity;
         Scalar transmissibility;
     };
 
@@ -2473,6 +2483,8 @@ private:
                 }
                 if constexpr (enableDiffusion)
                     *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
+                if (enableDispersion)
+                    dofData.dispersivity = transmissibilities_.dispersivity(globalCenterElemIdx, globalElemIdx);
             }
         };
 

ebos/eclproblem_properties.hh

@@ -269,6 +269,12 @@ struct EnableDiffusion<TypeTag, TTag::EclBaseProblem> {
     static constexpr bool value = true;
 };
 
+// Enable dispersion
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::EclBaseProblem> {
+    static constexpr bool value = false;
+};
+
 // only write the solutions for the report steps to disk
 template<class TypeTag>
 struct EnableWriteAllSolutions<TypeTag, TTag::EclBaseProblem> {

ebos/ecltransmissibility.hh

@@ -66,7 +66,8 @@ public:
                         const Grid& grid,
                         std::function<std::array<double,dimWorld>(int)> centroids,
                         bool enableEnergy,
-                        bool enableDiffusivity);
+                        bool enableDiffusivity,
+                        bool enableDispersivity);
 
     /*!
      * \brief Return the permeability for an element.
@@ -107,6 +108,11 @@ public:
      */
     Scalar diffusivity(unsigned elemIdx1, unsigned elemIdx2) const;
 
+    /*!
+     * \brief Return the dispersivity for the intersection between two elements.
+     */
+    Scalar dispersivity(unsigned elemIdx1, unsigned elemIdx2) const;
+
     /*!
      * \brief Actually compute the transmissibility over a face as a pre-compute step.
      *
@@ -230,6 +236,8 @@ protected:
 
     void extractPorosity_();
 
+    void extractDispersion_();
+
     void computeHalfTrans_(Scalar& halfTrans,
                            const DimVector& areaNormal,
                            int faceIdx, // in the reference element that contains the intersection
@@ -258,6 +266,7 @@ protected:
 
     std::vector<DimMatrix> permeability_;
     std::vector<Scalar> porosity_;
+    std::vector<Scalar> dispersion_;
     std::unordered_map<std::uint64_t, Scalar> trans_;
     const EclipseState& eclState_;
     const GridView& gridView_;
@@ -269,8 +278,10 @@ protected:
     std::map<std::pair<unsigned, unsigned>, Scalar> thermalHalfTransBoundary_;
     bool enableEnergy_;
     bool enableDiffusivity_;
+    bool enableDispersivity_;
     std::unordered_map<std::uint64_t, Scalar> thermalHalfTrans_; //NB this is based on direction map size is ca 2*trans_ (diffusivity_)
     std::unordered_map<std::uint64_t, Scalar> diffusivity_;
+    std::unordered_map<std::uint64_t, Scalar> dispersivity_;
 
     const LookUpData<Grid,GridView> lookUpData_;
     const LookUpCartesianData<Grid,GridView> lookUpCartesianData_;

ebos/ecltransmissibility_impl.hh

@@ -93,7 +93,8 @@ EclTransmissibility(const EclipseState& eclState,
                     const Grid& grid,
                     std::function<std::array<double,dimWorld>(int)> centroids,
                     bool enableEnergy,
-                    bool enableDiffusivity)
+                    bool enableDiffusivity,
+                    bool enableDispersivity)
       : eclState_(eclState)
       , gridView_(gridView)
       , cartMapper_(cartMapper)
@@ -101,6 +102,7 @@ EclTransmissibility(const EclipseState& eclState,
       , centroids_(centroids)
       , enableEnergy_(enableEnergy)
       , enableDiffusivity_(enableDiffusivity)
+      , enableDispersivity_(enableDispersivity)
       , lookUpData_(gridView)
       , lookUpCartesianData_(gridView, cartMapper)
 {
@@ -147,6 +149,17 @@ diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
 
 }
 
+template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
+Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
+dispersivity(unsigned elemIdx1, unsigned elemIdx2) const
+{
+    if (dispersivity_.empty())
+        return 0.0;
+
+    return dispersivity_.at(isId(elemIdx1, elemIdx2));
+
+}
+
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
 void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 update(bool global, const std::function<unsigned int(unsigned int)>& map, const bool applyNncMultregT)
@@ -159,7 +172,8 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
     unsigned numElements = elemMapper.size();
     // get the ntg values, the ntg values are modified for the cells merged with minpv
     const std::vector<double>& ntg = this->lookUpData_.assignFieldPropsDoubleOnLeaf(eclState_.fieldProps(), "NTG", numElements);
-    const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive() && enableDiffusivity_;
+    const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive();
+    const bool updateDispersivity = eclState_.getSimulationConfig().rock_config().dispersion();
 
     if (map)
         extractPermeability_(map);
@@ -209,6 +223,13 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
         extractPorosity_();
     }
 
+    // if dispersion is enabled, let's do the same for the "dispersivity"
+    if (updateDispersivity) {
+        dispersivity_.clear();
+        dispersivity_.reserve(numElements*3*1.05);
+        extractDispersion_();
+    }
+
     // The MULTZ needs special case if the option is ALL
     // Then the smallest multiplier is applied.
     // Default is to apply the top and bottom multiplier
@@ -314,6 +335,9 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
                 if (updateDiffusivity) {
                     diffusivity_[isId(elemIdx, outsideElemIdx)] = 0.0;
                 }
+                if (updateDispersivity) {
+                    dispersivity_[isId(elemIdx, outsideElemIdx)] = 0.0;
+                }
                 continue;
             }
 
@@ -485,6 +509,42 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
 
                 diffusivity_[isId(elemIdx, outsideElemIdx)] = diffusivity;
            }
+
+           // update the "dispersivity half transmissibility" for the intersection
+            if (updateDispersivity) {
+
+                Scalar halfDispersivity1;
+                Scalar halfDispersivity2;
+
+                computeHalfDiffusivity_(halfDispersivity1,
+                                        faceAreaNormal,
+                                        distanceVector_(faceCenterInside,
+                                                        intersection.indexInInside(),
+                                                        elemIdx,
+                                                        axisCentroids),
+                                        dispersion_[elemIdx]);
+                computeHalfDiffusivity_(halfDispersivity2,
+                                        faceAreaNormal,
+                                        distanceVector_(faceCenterOutside,
+                                                        intersection.indexInOutside(),
+                                                        outsideElemIdx,
+                                                        axisCentroids),
+                                        dispersion_[outsideElemIdx]);
+
+                applyNtg_(halfDispersivity1, insideFaceIdx, elemIdx, ntg);
+                applyNtg_(halfDispersivity2, outsideFaceIdx, outsideElemIdx, ntg);
+
+                //TODO Add support for multipliers
+                Scalar dispersivity;
+                if (std::abs(halfDispersivity1) < 1e-30 || std::abs(halfDispersivity2) < 1e-30)
+                    // avoid division by zero
+                    dispersivity = 0.0;
+                else
+                    dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2);
+
+
+                dispersivity_[isId(elemIdx, outsideElemIdx)] = dispersivity;
+           }
         }
     }
 
@@ -615,6 +675,18 @@ extractPorosity_()
                                "(The PORO keywords are missing)");
 }
 
+template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
+void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
+extractDispersion_()
+{
+    if (!enableDispersivity_) {
+        throw std::runtime_error("Dispersion disabled at compile time, but the deck "
+                                 "contains the DISPERC keyword.");
+    }
+    const auto& fp = eclState_.fieldProps();
+    dispersion_ = fp.get_double("DISPERC");
+}
+
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
 void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 removeSmallNonCartesianTransmissibilities_()

flow/flow_ebos_gasoil_energy.cpp

@@ -70,6 +70,9 @@ struct LocalResidual<TypeTag, TTag::EclFlowGasOilEnergyProblem> { using type = B
 template<class TypeTag>
 struct EnableDiffusion<TypeTag, TTag::EclFlowGasOilEnergyProblem> { static constexpr bool value = true; };
 
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::EclFlowGasOilEnergyProblem> { static constexpr bool value = true; };
+
 }}
 
 namespace Opm {

flow/flow_ebos_gasoildiffuse.cpp

@@ -41,10 +41,12 @@ struct Linearizer<TypeTag, TTag::EclFlowGasOilDiffuseProblem> { using type = Tpf
 template<class TypeTag>
 struct LocalResidual<TypeTag, TTag::EclFlowGasOilDiffuseProblem> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
 
-
 template<class TypeTag>
 struct EnableDiffusion<TypeTag, TTag::EclFlowGasOilDiffuseProblem> { static constexpr bool value = true; };
 
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::EclFlowGasOilDiffuseProblem> { static constexpr bool value = true; };
+
 //! The indices required by the model
 template<class TypeTag>
 struct Indices<TypeTag, TTag::EclFlowGasOilDiffuseProblem>

flow/flow_ebos_gaswater_dissolution_diffuse.cpp

@@ -48,6 +48,9 @@ struct LocalResidual<TypeTag, TTag::EclFlowGasWaterDissolutionDiffuseProblem> {
 template<class TypeTag>
 struct EnableDiffusion<TypeTag, TTag::EclFlowGasWaterDissolutionDiffuseProblem> { static constexpr bool value = true; };
 
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::EclFlowGasWaterDissolutionDiffuseProblem> { static constexpr bool value = true; };
+
 template<class TypeTag>
 struct EnableDisgasInWater<TypeTag, TTag::EclFlowGasWaterDissolutionDiffuseProblem> {
     static constexpr bool value = true;

flow/flow_ebos_gaswater_energy.cpp

@@ -46,7 +46,10 @@ template<class TypeTag>
 struct LocalResidual<TypeTag, TTag::EclFlowGasWaterEnergyProblem> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
 
 template<class TypeTag>
-struct EnableDiffusion<TypeTag, TTag::EclFlowGasWaterEnergyProblem> { static constexpr bool value = false; };
+struct EnableDiffusion<TypeTag, TTag::EclFlowGasWaterEnergyProblem> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::EclFlowGasWaterEnergyProblem> { static constexpr bool value = true; };
 
 template<class TypeTag>
 struct EnableEnergy<TypeTag, TTag::EclFlowGasWaterEnergyProblem> {

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -133,6 +133,10 @@ template<class TypeTag>
 struct EnableMICP<TypeTag, TTag::EclFlowProblem> {
     static constexpr bool value = false;
 };
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::EclFlowProblem> {
+    static constexpr bool value = false;
+};
 
 template<class TypeTag>
 struct EclWellModel<TypeTag, TTag::EclFlowProblem> {