- faster updateProperty

- refactored for making local updating in inhereted classes

ebos/eclproblem.hh

@@ -1052,22 +1052,7 @@ public:
 
         // update maximum water saturation and minimum pressure
         // used when ROCKCOMP is activated
-        const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
+        asImp_().updateExplicitQuantitites_();
-        const bool invalidateFromMinPressure = updateMinPressure_();
-
-        // update hysteresis and max oil saturation used in vappars
-        const bool invalidateFromHyst = updateHysteresis_();
-        const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
-
-        // the derivatives may have change
-        bool invalidateIntensiveQuantities = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
-        if (invalidateIntensiveQuantities){
-            OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
-            this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
-        }
-        
-        if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
-            updateMaxPolymerAdsorption_();
 
         wellModel_.beginTimeStep();
         if (enableAquifers_)
@@ -1126,9 +1111,10 @@ public:
         tracerModel_.endTimeStep();
 
         // deal with DRSDT and DRVDT
-        updateCompositionChangeLimits_();
+        asImp_().updateCompositionChangeLimits_();
         {
         OPM_TIMEBLOCK(driftCompansation);
+        asImp_().updateCompositionChangeLimits_();
         if (enableDriftCompensation_) {
             const auto& residual = this->model().linearizer().residual();
             for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
@@ -2081,24 +2067,50 @@ public:
         serializer(tracerModel_);
         serializer(*materialLawManager_);
         serializer(*eclWriter_);
+    }    
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation *>(this); }
+    
+    void updateExplicitQuantitites_()
+    {
+        OPM_TIMEBLOCK(updateExplicitQuantities);
+        const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
+        const bool invalidateFromMinPressure = updateMinPressure_();
+
+        // update hysteresis and max oil saturation used in vappars
+        const bool invalidateFromHyst = updateHysteresis_();
+        const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
+
+        // the derivatives may have change
+        bool invalidateIntensiveQuantities
+            = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
+        if (invalidateIntensiveQuantities) {
+            OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
+            this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
+        }
+
+        if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
+            updateMaxPolymerAdsorption_();
     }
 
-private:
     template<class UpdateFunc>
     void updateProperty_(const std::string& failureMsg,
                          UpdateFunc func)
     {
         OPM_TIMEBLOCK(updateProperty);
-        ElementContext elemCtx(this->simulator());
+        const auto& problem = this->simulator().problem();
+        const auto& model = this->simulator().model();
+        const auto& primaryVars = model.solution(/*timeIdx*/0);
         const auto& vanguard = this->simulator().vanguard();
+        size_t numGridDof = primaryVars.size();
         OPM_BEGIN_PARALLEL_TRY_CATCH();
-        for (const auto& elem : elements(vanguard.gridView())) {
+#ifdef _OPENMP
-            elemCtx.updatePrimaryStencil(elem);
+#pragma omp parallel for
-            elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+#endif
-
+        for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
-            unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+                const auto& iq = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/ 0);
-            const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
+                func(dofIdx, iq);
-            func(compressedDofIdx, iq);
         }
         OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
     }
@@ -2119,60 +2131,71 @@ private:
         this->updateProperty_("EclProblem::updateCompositionChangeLimits_()) failed:",
                               [this,episodeIdx,active](unsigned compressedDofIdx, const IntensiveQuantities& iq)
                               {
-                                  auto& simulator = this->simulator();
+                                  this->updateCompositionChangeLimits_(compressedDofIdx,
-                                  auto& vanguard = simulator.vanguard();
+                                                                       iq,
-                                  if (active[0]) {
+                                                                       episodeIdx,
-                                      // This implements the convective DRSDT as described in
+                                                                       active);
-                                      // Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow simulator"
+                              }
-                                      // Submitted to TCCS 11, 2021
+            );
-                                      const Scalar g = this->gravity_[dim - 1];
+    }
-                                      const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
-                                      const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
-                                      const Scalar distZ = vanguard.cellThickness(compressedDofIdx);
-                                      const auto& fs = iq.fluidState();
-                                      const Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx));
-                                      const Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
-                                      const Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx));
-                                      const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(),t,p);
-                                      const Scalar saturatedInvB = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p);
-                                      const Scalar rsZero = 0.0;
-                                      const Scalar pureDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p,rsZero) * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
-                                      const Scalar saturatedDensity = saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex()) + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
-                                      const Scalar deltaDensity = saturatedDensity - pureDensity;
-                                      const Scalar rs = getValue(fs.Rs());
-                                      const Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(),t,p,rs);
-                                      const Scalar poro =  getValue(iq.porosity());
-                                      // Note that for so = 0 this gives no limits (inf) for the dissolution rate
-                                      // Also we restrict the effect of convective mixing to positive density differences
-                                      // i.e. we only allow for fingers moving downward
-                                      this->convectiveDrs_[compressedDofIdx] = permz * rssat * max(0.0, deltaDensity) * g / ( so * visc * distZ * poro);
-                                  }
 
-                                  if (active[1]) {
+    void updateCompositionChangeLimits_(unsigned compressedDofIdx, const IntensiveQuantities& iq,int episodeIdx, const std::array<bool,3>& active)
-                                      const auto& fs = iq.fluidState();
+    {
+        auto& simulator = this->simulator();
+        auto& vanguard = simulator.vanguard();
+        if (active[0]) {
+            // This implements the convective DRSDT as described in
+            // Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow
+            // simulator" Submitted to TCCS 11, 2021
+            const Scalar g = this->gravity_[dim - 1];
+            const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
+            const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
+            const Scalar distZ = vanguard.cellThickness(compressedDofIdx);
+            const auto& fs = iq.fluidState();
+            const Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx));
+            const Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
+            const Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx));
+            const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), t, p);
+            const Scalar saturatedInvB
+                = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p);
+            const Scalar rsZero = 0.0;
+            const Scalar pureDensity
+                = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p, rsZero)
+                * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
+            const Scalar saturatedDensity = saturatedInvB
+                * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex())
+                   + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
+            const Scalar deltaDensity = saturatedDensity - pureDensity;
+            const Scalar rs = getValue(fs.Rs());
+            const Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(), t, p, rs);
+            const Scalar poro = getValue(iq.porosity());
+            // Note that for so = 0 this gives no limits (inf) for the dissolution rate
+            // Also we restrict the effect of convective mixing to positive density differences
+            // i.e. we only allow for fingers moving downward
+            this->convectiveDrs_[compressedDofIdx]
+                = permz * rssat * max(0.0, deltaDensity) * g / (so * visc * distZ * poro);
+        }
 
-                                      using FluidState = typename std::decay<decltype(fs)>::type;
+        if (active[1]) {
+            const auto& fs = iq.fluidState();
 
-                                      int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
+            using FluidState = typename std::decay<decltype(fs)>::type;
-                                      const auto& oilVaporizationControl = vanguard.schedule()[episodeIdx].oilvap();
-                                      if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
-                                          this->lastRs_[compressedDofIdx] =
-                                            BlackOil::template getRs_<FluidSystem,
-                                                                      FluidState,
-                                                                      Scalar>(fs, iq.pvtRegionIndex());
-                                      else
-                                        this->lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
-                                  }
 
-                                  if (active[2]) {
+            int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
-                                      const auto& fs = iq.fluidState();
+            const auto& oilVaporizationControl = vanguard.schedule()[episodeIdx].oilvap();
-                                      using FluidState = typename std::decay<decltype(fs)>::type;
+            if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
-                                      this->lastRv_[compressedDofIdx] =
+                this->lastRs_[compressedDofIdx]
-                                          BlackOil::template getRv_<FluidSystem,
+                    = BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
-                                                                    FluidState,
+            else
-                                                                    Scalar>(fs, iq.pvtRegionIndex());
+                this->lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
-                                  }
+        }
-                              });
+
+        if (active[2]) {
+            const auto& fs = iq.fluidState();
+            using FluidState = typename std::decay<decltype(fs)>::type;
+            this->lastRv_[compressedDofIdx]
+                = BlackOil::template getRv_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
+        }
     }
 
     bool updateMaxOilSaturation_()
@@ -2183,19 +2206,30 @@ private:
         // we use VAPPARS
         if (this->vapparsActive(episodeIdx)) {
             this->updateProperty_("EclProblem::updateMaxOilSaturation_() failed:",
-                                 [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
+                                  [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
-                                 {
+                                  {
-                                     const auto& fs = iq.fluidState();
+                                      this->updateMaxOilSaturation_(compressedDofIdx,iq);
-                                     const Scalar So = decay<Scalar>(fs.saturation(oilPhaseIdx));
+                                  });
-                                     auto& mos = this->maxOilSaturation_;
-                                     mos[compressedDofIdx] = std::max(mos[compressedDofIdx], So);
-                                 });
             return true;
         }
 
         return false;
     }
 
+    bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq) 
+    {
+        OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
+        const auto& fs = iq.fluidState();
+        const Scalar So = decay<Scalar>(fs.saturation(oilPhaseIdx));
+        auto& mos = this->maxOilSaturation_;
+        if(mos[compressedDofIdx], So){
+            mos[compressedDofIdx] = So;
+            return true;
+        }else{
+            return false;
+        }
+    }
+
     bool updateMaxWaterSaturation_()
     {
         OPM_TIMEBLOCK(updateMaxWaterSaturation);
@@ -2207,14 +2241,26 @@ private:
         this->updateProperty_("EclProblem::updateMaxWaterSaturation_() failed:",
                               [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
                               {
-                                  const auto& fs = iq.fluidState();
+                                  this->updateMaxWaterSaturation_(compressedDofIdx,iq);
-                                  const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
-                                  auto& mow = this->maxWaterSaturation_;
-                                  mow[compressedDofIdx] = std::max(mow[compressedDofIdx], Sw);
                                });
         return true;
     }
 
+
+    bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
+    {
+        OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation);
+        const auto& fs = iq.fluidState();
+        const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
+        auto& mow = this->maxWaterSaturation_;
+        if(mow[compressedDofIdx]< Sw){
+            mow[compressedDofIdx] = Sw;
+            return true;
+        }else{
+            return false;
+        }
+    }
+
     bool updateMinPressure_()
     {
         OPM_TIMEBLOCK(updateMinPressure);
@@ -2225,14 +2271,24 @@ private:
         this->updateProperty_("EclProblem::updateMinPressure_() failed:",
                               [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
                               {
-                                  const auto& fs = iq.fluidState();
+                                  this->updateMinPressure_(compressedDofIdx,iq);
-                                  const Scalar mo = getValue(fs.pressure(oilPhaseIdx));
-                                  auto& mos = this->minOilPressure_;
-                                  mos[compressedDofIdx] = std::min(mos[compressedDofIdx], mo);
                               });
         return true;
     }
 
+    bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
+        OPM_TIMEBLOCK_LOCAL(updateMinPressure);
+        const auto& fs = iq.fluidState();
+        const Scalar mo = getValue(fs.pressure(oilPhaseIdx));
+        auto& mos = this->minOilPressure_;
+        if(mos[compressedDofIdx]> mo){
+            mos[compressedDofIdx] = mo;
+            return true;
+        }else{
+            return false;
+        }
+    }
+
     void readMaterialParameters_()
     {
         OPM_TIMEBLOCK(readMaterialParameters);
@@ -2718,18 +2774,37 @@ private:
         return true;
     }
 
+
+    bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
+    {
+        OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
+        materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
+        //TODO change materials to give a bool
+        return true;
+    }
+
     void updateMaxPolymerAdsorption_()
     {
         // we need to update the max polymer adsoption data for all elements
         this->updateProperty_("EclProblem::updateMaxPolymerAdsorption_() failed:",
                               [this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
                               {
-                                  const Scalar pa = scalarValue(iq.polymerAdsorption());
+                                  this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
-                                  auto& mpa = this->maxPolymerAdsorption_;
-                                  mpa[compressedDofIdx] = std::max(mpa[compressedDofIdx], pa);
                               });
     }
 
+    bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
+    {
+        const Scalar pa = scalarValue(iq.polymerAdsorption());
+        auto& mpa = this->maxPolymerAdsorption_;
+        if(mpa[compressedDofIdx]<pa){
+            mpa[compressedDofIdx] = pa;
+            return true;
+        }else{
+            return false;
+        }
+    }
+
     struct PffDofData_
     {
         ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;