Merge pull request #5081 from totto82/fix_gaswater_rockcomp

Use gas pressure for gas-water cases in rockcomp

ebos/eclgenericproblem.hh

@@ -282,7 +282,7 @@ public:
         serializer(maxOilSaturation_);
         serializer(polymer_);
         serializer(maxWaterSaturation_);
-        serializer(minOilPressure_);
+        serializer(minRefPressure_);
         serializer(overburdenPressure_);
         serializer(solventSaturation_);
         serializer(solventRsw_);
@@ -363,7 +363,7 @@ protected:
     PolymerSolutionContainer<Scalar> polymer_;
     std::vector<Scalar> maxOilSaturation_;
     std::vector<Scalar> maxWaterSaturation_;
-    std::vector<Scalar> minOilPressure_;
+    std::vector<Scalar> minRefPressure_;
     std::vector<Scalar> overburdenPressure_;
     std::vector<Scalar> solventSaturation_;
     std::vector<Scalar> solventRsw_;

ebos/eclgenericproblem_impl.hh

@@ -103,7 +103,7 @@ serializationTestObject(const EclipseState& eclState,
     EclGenericProblem result(eclState, schedule, gridView);
     result.maxOilSaturation_ = {1.0, 2.0};
     result.maxWaterSaturation_ = {6.0};
-    result.minOilPressure_ = {7.0, 8.0, 9.0, 10.0};
+    result.minRefPressure_ = {7.0, 8.0, 9.0, 10.0};
     result.overburdenPressure_ = {11.0};
     result.solventSaturation_ = {15.0};
     result.solventRsw_ = {18.0};
@@ -238,7 +238,7 @@ readRockCompactionParameters_()
     case RockConfig::Hysteresis::IRREVERS:
         // interpolate the porv volume multiplier using the minimum pressure in the cell
         // i.e. don't allow re-inflation.
-        minOilPressure_.resize(numElem, 1e99);
+        minRefPressure_.resize(numElem, 1e99);
         break;
     default:
         throw std::runtime_error("Not support ROCKOMP hysteresis option ");
@@ -597,10 +597,10 @@ template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 minOilPressure(unsigned globalDofIdx) const
 {
-    if (minOilPressure_.empty())
+    if (minRefPressure_.empty())
         return 0.0;
 
-    return minOilPressure_[globalDofIdx];
+    return minRefPressure_[globalDofIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
@@ -775,7 +775,7 @@ bool EclGenericProblem<GridView,FluidSystem,Scalar>::
 operator==(const EclGenericProblem& rhs) const
 {
     return this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
-           this->minOilPressure_ == rhs.minOilPressure_ &&
+           this->minRefPressure_ == rhs.minRefPressure_ &&
            this->overburdenPressure_ == rhs.overburdenPressure_ &&
            this->solventSaturation_ == rhs.solventSaturation_ &&
            this->solventRsw_ == rhs.solventRsw_ &&

ebos/eclproblem.hh

@@ -1544,17 +1544,7 @@ public:
                         throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
                         break;
                 }
-                int phaseIndex;
+                double pressure = initialFluidStates_[globalDofIdx].pressure(refPressurePhaseIdx_());
-                if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
-                    phaseIndex = oilPhaseIdx;
-                }
-                else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
-                    phaseIndex = gasPhaseIdx;
-                }
-                else {
-                    phaseIndex = waterPhaseIdx;
-                }
-                double pressure = initialFluidStates_[globalDofIdx].pressure(phaseIndex);
                 const auto pressure_input = bc.pressure;
                 if (pressure_input) {
                     pressure = *pressure_input;
@@ -1578,7 +1568,7 @@ public:
                         fluidState.setPressure(phaseIdx, pressure);
                 }
                 
-                double temperature = initialFluidStates_[globalDofIdx].temperature(phaseIndex);
+                double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature
                 const auto temperature_input = bc.temperature;
                 if(temperature_input)
                     temperature = *temperature_input;
@@ -1656,19 +1646,19 @@ public:
             tableIdx = this->rockTableIdx_[elementIdx];
 
         const auto& fs = intQuants.fluidState();
-        LhsEval effectiveOilPressure = decay<LhsEval>(fs.pressure(oilPhaseIdx));
+        LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
-        if (!this->minOilPressure_.empty())
+        if (!this->minRefPressure_.empty())
             // The pore space change is irreversible
-            effectiveOilPressure =
+            effectivePressure =
-                min(decay<LhsEval>(fs.pressure(oilPhaseIdx)),
+                min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
-                                   this->minOilPressure_[elementIdx]);
+                                   this->minRefPressure_[elementIdx]);
 
         if (!this->overburdenPressure_.empty())
-            effectiveOilPressure -= this->overburdenPressure_[elementIdx];
+            effectivePressure -= this->overburdenPressure_[elementIdx];
 
 
         if (!this->rockCompPoroMult_.empty()) {
-            return this->rockCompPoroMult_[tableIdx].eval(effectiveOilPressure, /*extrapolation=*/true);
+            return this->rockCompPoroMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
         }
 
         // water compaction
@@ -1676,7 +1666,7 @@ public:
         LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
         LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
 
-        return this->rockCompPoroMultWc_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true);
+        return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
     }
 
     /*!
@@ -1696,26 +1686,26 @@ public:
             tableIdx = this->rockTableIdx_[elementIdx];
 
         const auto& fs = intQuants.fluidState();
-        LhsEval effectiveOilPressure = decay<LhsEval>(fs.pressure(oilPhaseIdx));
+        LhsEval effectivePressure = decay<LhsEval>(fs.pressure(refPressurePhaseIdx_()));
 
-        if (!this->minOilPressure_.empty())
+        if (!this->minRefPressure_.empty())
             // The pore space change is irreversible
-            effectiveOilPressure =
+            effectivePressure =
-                min(decay<LhsEval>(fs.pressure(oilPhaseIdx)),
+                min(decay<LhsEval>(fs.pressure(refPressurePhaseIdx_())),
-                    this->minOilPressure_[elementIdx]);
+                    this->minRefPressure_[elementIdx]);
 
         if (!this->overburdenPressure_.empty())
-            effectiveOilPressure -= this->overburdenPressure_[elementIdx];
+            effectivePressure -= this->overburdenPressure_[elementIdx];
 
         if (!this->rockCompTransMult_.empty())
-            return this->rockCompTransMult_[tableIdx].eval(effectiveOilPressure, /*extrapolation=*/true);
+            return this->rockCompTransMult_[tableIdx].eval(effectivePressure, /*extrapolation=*/true);
 
         // water compaction
         assert(!this->rockCompTransMultWc_.empty());
         LhsEval SwMax = max(decay<LhsEval>(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]);
         LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx);
 
-        return this->rockCompTransMultWc_[tableIdx].eval(effectiveOilPressure, SwDeltaMax, /*extrapolation=*/true);
+        return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true);
     }
 
     /*!
@@ -1919,7 +1909,7 @@ protected:
     {
         OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
         const auto& fs = iq.fluidState();
-        const Scalar So = decay<Scalar>(fs.saturation(oilPhaseIdx));
+        const Scalar So = decay<Scalar>(fs.saturation(refPressurePhaseIdx_()));
         auto& mos = this->maxOilSaturation_;
         if(mos[compressedDofIdx] < So){
             mos[compressedDofIdx] = So;
@@ -1964,7 +1954,7 @@ protected:
     {
         OPM_TIMEBLOCK(updateMinPressure);
         // IRREVERS option is used in ROCKCOMP
-        if (this->minOilPressure_.empty())
+        if (this->minRefPressure_.empty())
             return false;
 
         this->updateProperty_("EclProblem::updateMinPressure_() failed:",
@@ -1978,10 +1968,10 @@ protected:
     bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
         OPM_TIMEBLOCK_LOCAL(updateMinPressure);
         const auto& fs = iq.fluidState();
-        const Scalar mo = getValue(fs.pressure(oilPhaseIdx));
+        const Scalar min_pressure = getValue(fs.pressure(refPressurePhaseIdx_()));
-        auto& mos = this->minOilPressure_;
+        auto& min_pressures = this->minRefPressure_;
-        if(mos[compressedDofIdx]> mo){
+        if(min_pressures[compressedDofIdx]> min_pressure){
-            mos[compressedDofIdx] = mo;
+            min_pressures[compressedDofIdx] = min_pressure;
             return true;
         }else{
             return false;
@@ -2123,8 +2113,8 @@ protected:
                 this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx));
             if (!this->maxOilSaturation_.empty())
                 this->maxOilSaturation_[elemIdx] = std::max(this->maxOilSaturation_[elemIdx], fs.saturation(oilPhaseIdx));
-            if (!this->minOilPressure_.empty())
+            if (!this->minRefPressure_.empty())
-                this->minOilPressure_[elemIdx] = std::min(this->minOilPressure_[elemIdx], fs.pressure(oilPhaseIdx));
+                this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_()));
         }
 
 
@@ -2694,6 +2684,18 @@ private:
         }
     }
 
+    int refPressurePhaseIdx_() const {
+        if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+            return oilPhaseIdx;
+        }
+        else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
+            return gasPhaseIdx;
+        }
+        else {
+            return waterPhaseIdx;
+        }
+    }
+
     typename Vanguard::TransmissibilityType transmissibilities_;
 
     std::shared_ptr<EclMaterialLawManager> materialLawManager_;