runs first steps of norne

ebos/eclthresholdpressure.hh

@@ -39,6 +39,7 @@
 #include <algorithm>
 #include <vector>
 
+#include <opm/models/discretization/common/smallelementcontext.hh>
 namespace Opm {
 
 /*!
@@ -61,11 +62,15 @@ class EclThresholdPressure : public EclGenericThresholdPressure<GetPropType<Type
                                                  GetPropType<TypeTag, Properties::GridView>,
                                                  GetPropType<TypeTag, Properties::ElementMapper>,
                                                  GetPropType<TypeTag, Properties::Scalar>>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
     using Simulator = GetPropType<TypeTag, Properties::Simulator>;
     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 GridView = GetPropType<TypeTag, Properties::GridView>;
+    enum { dimWorld = GridView::dimensionworld };
 
     enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
     enum { numPhases = FluidSystem::numPhases };
@@ -95,6 +100,95 @@ public:
     }
 
 private:
+    template<class Face,class Stencil,class ElemCtx>
+    double calculateMaxDp(Face& face, Stencil& stencil,
+                          ElemCtx& elemCtx,const unsigned& scvfIdx,
+                          const unsigned& i,const unsigned& j,const unsigned& insideElemIdx,const unsigned& outsideElemIdx){
+        typedef MathToolbox<Evaluation> Toolbox;
+        elemCtx.updateIntensiveQuantities(/*timeIdx=*/0);
+        elemCtx.updateExtensiveQuantities(/*timeIdx=*/0);
+        // determine the maximum difference of the pressure of any phase over the
+        // intersection
+        Scalar pth = 0.0;
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0);
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
+            const auto& up = elemCtx.intensiveQuantities(upIdx, /*timeIdx=*/0);
+            
+            if (up.mobility(phaseIdx) > 0.0) {
+                Scalar phaseVal = Toolbox::value(extQuants.pressureDifference(phaseIdx));
+                pth = std::max(pth, std::abs(phaseVal));
+            }
+        }
+        return pth;
+    }
+    
+    template<class Face,class Stencil>
+    double calculateMaxDp(Face& face, Stencil& stencil,
+                          SmallElementContext<TypeTag>& elemCtx,const unsigned& scvfIdx,
+                          const unsigned& i,const unsigned& j,
+                          const unsigned& insideElemIdx,const unsigned& outsideElemIdx){        
+        typedef MathToolbox<Evaluation> Toolbox;
+        // determine the maximum difference of the pressure of any phase over the
+        // intersection
+        Scalar pth = 0.0;
+        //const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0);
+        
+        Scalar Vin = elemCtx.dofVolume(i, /*timeIdx=*/0);
+        Scalar Vex = elemCtx.dofVolume(j, /*timeIdx=*/0);
+        
+        Scalar thpres = 0.0;//NB ??problem.thresholdPressure(globalIndexIn, globalIndexEx);
+        
+        // estimate the gravity correction: for performance reasons we use a simplified
+        // approach for this flux module that assumes that gravity is constant and always
+        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
+        const auto& problem = elemCtx.problem();
+        Scalar g = problem.gravity()[dimWorld - 1];
+        
+        const auto& intQuantsIn = elemCtx.intensiveQuantities(i, /*timeIdx*/0);
+        const auto& intQuantsEx = elemCtx.intensiveQuantities(j, /*timeIdx*/0);
+        
+        // this is quite hacky because the dune grid interface does not provide a
+        // cellCenterDepth() method (so we ask the problem to provide it). The "good"
+        // solution would be to take the Z coordinate of the element centroids, but since
+        // ECL seems to like to be inconsistent on that front, it needs to be done like
+        // here...
+        Scalar zIn = problem.dofCenterDepth(elemCtx, i, /*timeIdx*/0);
+        Scalar zEx = problem.dofCenterDepth(elemCtx, j, /*timeIdx*/0);
+        
+        // the distances from the DOF's depths. (i.e., the additional depth of the
+        // exterior DOF)
+        Scalar distZ = zIn - zEx;
+        
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            short dnIdx;
+            //
+            short upIdx;
+            Evaluation pressureDifference;
+            ExtensiveQuantities::calculatePhasePressureDiff_(upIdx,
+                                                             dnIdx,
+                                                             pressureDifference,
+                                                             intQuantsIn,
+                                                             intQuantsEx,
+                                                             scvfIdx,//input
+                                                             /*timeIdx*/0,//input
+                                                             phaseIdx,//input
+                                                             i,//input
+                                                             j,//intput
+                                                             Vin,
+                                                             Vex,
+                                                             insideElemIdx,
+                                                             outsideElemIdx,
+                                                             distZ*g,
+                                                             thpres);
+            const IntensiveQuantities& up = (upIdx == i) ? intQuantsIn : intQuantsEx;
+            if (up.mobility(phaseIdx) > 0.0) {
+                Scalar phaseVal = Toolbox::value(pressureDifference);
+                pth = std::max(pth, std::abs(phaseVal));
+            }
+        }
+        return pth;
+    }
     // compute the defaults of the threshold pressures using the initial condition
     void computeDefaultThresholdPressures_()
     {
@@ -107,16 +201,21 @@ private:
         auto elemIt = gridView.template begin</*codim=*/ 0>();
         const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
         ElementContext elemCtx(simulator_);
+        simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
         for (; elemIt != elemEndIt; ++elemIt) {
 
             const auto& elem = *elemIt;
             if (elem.partitionType() != Dune::InteriorEntity)
                 continue;
 
-            elemCtx.updateAll(elem);
+            elemCtx.updateStencil(elem);
+            //
+            
             const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
             
+
             for (unsigned scvfIdx = 0; scvfIdx < stencil.numInteriorFaces(); ++ scvfIdx) {
+                
                 const auto& face = stencil.interiorFace(scvfIdx);
 
                 unsigned i = face.interiorIndex();
@@ -131,26 +230,17 @@ private:
                 if (equilRegionInside == equilRegionOutside)
                     // the current face is not at the boundary between EQUIL regions!
                     continue;
-
+                const auto& problem = elemCtx.problem();
                 // don't include connections with negligible flow
-                const Evaluation& trans = simulator_.problem().transmissibility(elemCtx, i, j);
+                const Evaluation& trans = problem.transmissibility(elemCtx, i, j);
                 Scalar faceArea = face.area();
                 if (std::abs(faceArea*getValue(trans)) < 1e-18)
                     continue;
                 
-                // determine the maximum difference of the pressure of any phase over the
+                double pth = calculateMaxDp(face, stencil, elemCtx, scvfIdx,
-                // intersection
+                                            i, j,
-                Scalar pth = 0.0;
+                                            insideElemIdx, outsideElemIdx);
-                const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0);
+                // don't include connections with negligible flow
-                for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                    unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
-                    const auto& up = elemCtx.intensiveQuantities(upIdx, /*timeIdx=*/0);
-
-                    if (up.mobility(phaseIdx) > 0.0) {
-                        Scalar phaseVal = Toolbox::value(extQuants.pressureDifference(phaseIdx));
-                        pth = std::max(pth, std::abs(phaseVal));
-                    }
-                }
 
                 int offset1 = equilRegionInside*this->numEquilRegions_ + equilRegionOutside;
                 int offset2 = equilRegionOutside*this->numEquilRegions_ + equilRegionInside;

ebos/ecltracermodel.hh

@@ -86,6 +86,9 @@ class EclTracerModel : public EclGenericTracerModel<GetPropType<TypeTag, Propert
     enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
 
+    using Eval =  DenseAd::Evaluation<double, numEq>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
 public:
     EclTracerModel(Simulator& simulator)
         : BaseType(simulator.vanguard().gridView(),
@@ -224,7 +227,37 @@ protected:
             freeVolume = phaseVolume * variable<LhsEval>(1.0, 0);
 
     }
+    //template<class TypeTag>
+    void getVolumeFlux(unsigned& upIdx,
+                       Scalar& v,
+                       const FvBaseElementContext<TypeTag>& elemCtx,
+                       const int tracerPhaseIdx,
+                       unsigned scvfIdx
+        ){                       
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx*/ 0);
+        upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
+        v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
+    }
 
+    //template <class TypeTag>
+    void getVolumeFlux(unsigned& upIdx,
+                       Scalar& v,
+                       const SmallElementContext<TypeTag>& elemCtx,
+                       const int tracerPhaseIdx,
+                       unsigned scvfIdx
+        ){                       
+        short upIdxV[numPhases];
+        Eval volumFlux[numPhases];
+        Eval pressureDifferences[numPhases];
+        ExtensiveQuantities::volumeAndPhasePressureDifferences(upIdxV ,
+                                                               volumFlux,
+                                                               pressureDifferences,
+                                                               elemCtx,
+                                                               scvfIdx,
+                                                               /*timeIdx*/ 0);
+        v = decay<Scalar>(volumFlux[tracerPhaseIdx]);
+        upIdx = upIdxV[tracerPhaseIdx] ;
+    }
     // evaluate the flux(es) over one face
     void computeFlux_(TracerEvaluation & freeFlux,
                       bool & isUpFree,
@@ -236,17 +269,18 @@ protected:
     {
         const auto& stencil = elemCtx.stencil(timeIdx);
         const auto& scvf = stencil.interiorFace(scvfIdx);
-
+        unsigned inIdx  = scvf.interiorIndex();
-        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        unsigned upIdx;
-        unsigned inIdx = extQuants.interiorIndex();
+        Scalar v;
-
+        getVolumeFlux(upIdx,
-        unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
+                      v,
-
+                      elemCtx,
+                      tracerPhaseIdx,
+                      scvfIdx);        
         const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
         const auto& fs = intQuants.fluidState();
-
         Scalar A = scvf.area();
-        Scalar v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
+        
         Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx));
 
         if (inIdx == upIdx) {
@@ -397,7 +431,9 @@ protected:
         auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
         auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
         for (; elemIt != elemEndIt; ++ elemIt) {
-            elemCtx.updateAll(*elemIt);
+            //elemCtx.updateAll(*elemIt);
+            elemCtx.updatePrimaryStencil(*elemIt);
+            elemCtx.updatePrimaryIntensiveQuantities(/*timIdx*/ 0.0);
             int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timIdx=*/0);
             Scalar fVolume;
             computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/0);

opm/simulators/aquifers/AquiferNumerical.hpp

@@ -49,12 +49,14 @@ public:
     using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
 
     enum { dimWorld = GridView::dimensionworld };
-
+    enum { numPhases = FluidSystem::numPhases };
     static const int numEq = BlackoilIndices::numEq;
 
     using Eval =  DenseAd::Evaluation<double, numEq>;
     using Toolbox = MathToolbox<Eval>;
 
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
     // Constructor
     AquiferNumerical(const SingleNumericalAquifer& aquifer,
                      const std::unordered_map<int, int>& cartesian_to_compressed,
@@ -250,6 +252,30 @@ private:
         return sum_pressure_watervolume / sum_watervolume;
     }
 
+    template<class ElemCtx>
+    const double getWaterFlux(ElemCtx& elem_ctx,unsigned face_idx) const{
+        const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
+        const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
+        return water_flux;
+    }
+    
+    const double getWaterFlux(SmallElementContext<TypeTag>& elem_ctx,unsigned face_idx) const{
+        short upIdx[numPhases];
+        
+        Eval volumFlux[numPhases];
+        Eval pressureDifferences[numPhases];
+        ExtensiveQuantities::volumeAndPhasePressureDifferences(upIdx ,
+                                                               volumFlux,
+                                                               pressureDifferences,
+                                                               elem_ctx,
+                                                               face_idx,
+                                                               /*timeIdx*/ 0);
+        return Toolbox::value(volumFlux[this->phaseIdx_()]);
+    }
+                
+                
+
+    
     double calculateAquiferFluxRate() const
     {
         double aquifer_flux = 0.0;
@@ -276,8 +302,6 @@ private:
             if (idx != 0) {
                 continue;
             }
-            elem_ctx.updateAllIntensiveQuantities();
-            elem_ctx.updateAllExtensiveQuantities();
     
             const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
             // const auto &problem = elem_ctx.problem();
@@ -300,8 +324,10 @@ private:
                 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
                     continue;
                 }
-                const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
+                elem_ctx.updateAllIntensiveQuantities();
-                const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
+                elem_ctx.updateAllExtensiveQuantities();
+
+                const double water_flux = getWaterFlux(elem_ctx,face_idx);
                 
                 const std::size_t up_id = water_flux >= 0.0 ? i : j;
                 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
@@ -316,6 +342,7 @@ private:
 
         return aquifer_flux;
     }
+    
 };
 } // namespace Opm
 #endif