Merge pull request #3023 from totto82/addFlowCo2

Add diffusion in Flow

CMakeLists.txt

@@ -339,6 +339,26 @@ opm_add_test(flow_onephase
   DEPENDS opmsimulators
   LIBRARIES opmsimulators)
 
+opm_add_test(flow_co2_diffusion
+  ONLY_COMPILE
+  DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}
+  SOURCES
+  flow/flow_co2_diffusion.cpp
+  $<TARGET_OBJECTS:moduleVersion>
+  EXE_NAME flow_co2_diffusion
+  DEPENDS opmsimulators
+  LIBRARIES opmsimulators)
+
+opm_add_test(flow_bo_diffusion
+  ONLY_COMPILE
+  DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}
+  SOURCES
+  flow/flow_bo_diffusion.cpp
+  $<TARGET_OBJECTS:moduleVersion>
+  EXE_NAME flow_bo_diffusion
+  DEPENDS opmsimulators
+  LIBRARIES opmsimulators)
+
 opm_add_test(flow_onephase_energy
   ONLY_COMPILE
   DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}

ebos/eclproblem.hh

@@ -622,6 +622,7 @@ class EclProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
     enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
     enum { enableApiTracking = getPropValue<TypeTag, Properties::EnableApiTracking>() };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
@@ -1372,6 +1373,18 @@ public:
         return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
     }
 
+    /*!
+     * \copydoc EclTransmissiblity::diffusivity
+     */
+    template <class Context>
+    Scalar diffusivity(const Context& context,
+                       unsigned OPM_OPTIM_UNUSED fromDofLocalIdx,
+                       unsigned toDofLocalIdx) const
+    {
+        assert(fromDofLocalIdx == 0);
+        return *pffDofData_.get(context.element(), toDofLocalIdx).diffusivity;
+    }
+
     /*!
      * \copydoc EclTransmissiblity::transmissibilityBoundary
      */
@@ -3072,6 +3085,7 @@ private:
     struct PffDofData_
     {
         Opm::ConditionalStorage<enableEnergy, Scalar> thermalHalfTrans;
+        Opm::ConditionalStorage<enableDiffusion, Scalar> diffusivity;
         Scalar transmissibility;
     };
 
@@ -3093,6 +3107,8 @@ private:
 
                 if (enableEnergy)
                     *dofData.thermalHalfTrans = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
+                if (enableDiffusion)
+                    *dofData.diffusivity = transmissibilities_.diffusivity(globalCenterElemIdx, globalElemIdx);
             }
         };
 

ebos/ecltransmissibility.hh

@@ -72,6 +72,7 @@ class EclTransmissibility
     using Intersection = typename GridView::Intersection;
 
     static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
+    static const bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
 
     // Grid and world dimension
     enum { dimWorld = GridView::dimensionworld };
@@ -181,6 +182,13 @@ public:
             thermalHalfTransBoundary_.clear();
         }
 
+        // if diffusion is enabled, let's do the same for the "diffusivity"
+        if (enableDiffusion) {
+            diffusivity_->clear();
+            diffusivity_->reserve(numElements*3*1.05);
+            extractPorosity_();
+        }
+
         // 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
@@ -386,6 +394,42 @@ public:
                                                        faceDir);
 
                 trans_[isId_(elemIdx, outsideElemIdx)] = trans;
+
+                // update the "thermal half transmissibility" for the intersection
+                if (enableDiffusion) {
+
+                    Scalar halfDiffusivity1;
+                    Scalar halfDiffusivity2;
+
+                    computeHalfDiffusivity_(halfDiffusivity1,
+                                            faceAreaNormal,
+                                            distanceVector_(faceCenterInside,
+                                                            intersection.indexInInside(),
+                                                            elemIdx,
+                                                            axisCentroids),
+                                            porosity_[elemIdx]);
+                    computeHalfDiffusivity_(halfDiffusivity2,
+                                            faceAreaNormal,
+                                            distanceVector_(faceCenterOutside,
+                                                            intersection.indexInOutside(),
+                                                            outsideElemIdx,
+                                                            axisCentroids),
+                                            porosity_[outsideElemIdx]);
+
+                    applyNtg_(halfDiffusivity1, insideFaceIdx, elemIdx, ntg);
+                    applyNtg_(halfDiffusivity2, outsideFaceIdx, outsideElemIdx, ntg);
+
+                    //TODO Add support for multipliers
+                    Scalar diffusivity;
+                    if (std::abs(halfDiffusivity1) < 1e-30 || std::abs(halfDiffusivity2) < 1e-30)
+                        // avoid division by zero
+                        diffusivity = 0.0;
+                    else
+                        diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2);
+
+
+                    (*diffusivity_)[isId_(elemIdx, outsideElemIdx)] = diffusivity;
+               }
             }
         }
 
@@ -450,6 +494,12 @@ public:
     Scalar thermalHalfTransBoundary(unsigned insideElemIdx, unsigned boundaryFaceIdx) const
     { return thermalHalfTransBoundary_.at(std::make_pair(insideElemIdx, boundaryFaceIdx)); }
 
+    /*!
+     * \brief Return the diffusivity for the intersection between two elements.
+     */
+    Scalar diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
+    { return diffusivity_->at(isId_(elemIdx1, elemIdx2)); }
+
 private:
 
     void removeSmallNonCartesianTransmissibilities_()
@@ -881,6 +931,21 @@ private:
                                    "(The PERM{X,Y,Z} keywords are missing)");
     }
 
+    void extractPorosity_()
+    {
+        // read the intrinsic porosity from the eclState. Note that all arrays
+        // provided by eclState are one-per-cell of "uncompressed" grid, whereas the
+        // simulation grid might remove a few elements. (e.g. because it is distributed
+        // over several processes.)
+        const auto& fp = vanguard_.eclState().fieldProps();
+        if (fp.has_double("PORO")) {
+            porosity_ = fp.get_double("PORO");
+        }
+        else
+            throw std::logic_error("Can't read the porosityfrom the ecl state. "
+                                   "(The PORO keywords are missing)");
+    }
+
     std::uint64_t isId_(std::uint32_t elemIdx1, std::uint32_t elemIdx2) const
     {
         std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
@@ -924,6 +989,20 @@ private:
         halfTrans /= distance.two_norm2();
     }
 
+    void computeHalfDiffusivity_(Scalar& halfDiff,
+                                 const DimVector& areaNormal,
+                                 const DimVector& distance,
+                                 const Scalar& poro) const
+    {
+        halfDiff = poro;
+        Scalar val = 0;
+        for (unsigned i = 0; i < areaNormal.size(); ++i)
+            val += areaNormal[i]*distance[i];
+
+        halfDiff *= std::abs(val);
+        halfDiff /= distance.two_norm2();
+    }
+
     DimVector distanceVector_(const DimVector& center,
                               int faceIdx, // in the reference element that contains the intersection
                               unsigned elemIdx,
@@ -1000,11 +1079,14 @@ private:
     const Vanguard& vanguard_;
     Scalar transmissibilityThreshold_;
     std::vector<DimMatrix> permeability_;
+    std::vector<Scalar> porosity_;
     std::unordered_map<std::uint64_t, Scalar> trans_;
     std::map<std::pair<unsigned, unsigned>, Scalar> transBoundary_;
     std::map<std::pair<unsigned, unsigned>, Scalar> thermalHalfTransBoundary_;
     Opm::ConditionalStorage<enableEnergy,
                             std::unordered_map<std::uint64_t, Scalar> > thermalHalfTrans_;
+    Opm::ConditionalStorage<enableDiffusion,
+                            std::unordered_map<std::uint64_t, Scalar> > diffusivity_;
 };
 
 } // namespace Opm

flow/flow_bo_diffusion.cpp

@@ -0,0 +1,43 @@
+/*
+  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 3 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/>.
+*/
+#include "config.h"
+#include <opm/simulators/flow/Main.hpp>
+
+namespace Opm {
+namespace Properties {
+namespace TTag {
+struct EclFlowDiffusionProblem {
+    using InheritsFrom = std::tuple<EclFlowProblem>;
+};
+}
+
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::EclFlowDiffusionProblem> {
+    static constexpr bool value = true;
+};
+
+
+}}
+
+
+// ----------------- Main program -----------------
+int main(int argc, char** argv)
+{
+    using TypeTag = Opm::Properties::TTag::EclFlowDiffusionProblem;
+    auto mainObject = Opm::Main(argc, argv);
+    return mainObject.runStatic<TypeTag>();
+}

flow/flow_co2_diffusion.cpp

@@ -0,0 +1,64 @@
+/*
+  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 3 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/>.
+*/
+#include "config.h"
+#include <opm/simulators/flow/Main.hpp>
+#include <opm/models/blackoil/blackoiltwophaseindices.hh>
+
+namespace Opm {
+namespace Properties {
+namespace TTag {
+struct EclFlowCo2Problem {
+    using InheritsFrom = std::tuple<EclFlowProblem>;
+};
+}
+
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::EclFlowCo2Problem> {
+    static constexpr bool value = true;
+};
+
+//! The indices required by the model
+template<class TypeTag>
+struct Indices<TypeTag, TTag::EclFlowCo2Problem>
+{
+private:
+    // it is unfortunately not possible to simply use 'TypeTag' here because this leads
+    // to cyclic definitions of some properties. if this happens the compiler error
+    // messages unfortunately are *really* confusing and not really helpful.
+    using BaseTypeTag = TTag::EclFlowProblem;
+    using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
+
+public:
+    typedef Opm::BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                         getPropValue<TypeTag, Properties::EnableExtbo>(),
+                                         getPropValue<TypeTag, Properties::EnablePolymer>(),
+                                         getPropValue<TypeTag, Properties::EnableEnergy>(),
+                                         getPropValue<TypeTag, Properties::EnableFoam>(),
+                                         getPropValue<TypeTag, Properties::EnableBrine>(),
+                                         /*PVOffset=*/0,
+                                         /*disabledCompIdx=*/FluidSystem::waterCompIdx> type;
+};
+}}
+
+
+// ----------------- Main program -----------------
+int main(int argc, char** argv)
+{
+    using TypeTag = Opm::Properties::TTag::EclFlowCo2Problem;
+    auto mainObject = Opm::Main(argc, argv);
+    return mainObject.runStatic<TypeTag>();
+}