Merge pull request #3570 from svenn-t/ppcwmax

Implementation of PPCWMAX

CMakeLists_files.cmake

@@ -246,6 +246,7 @@ if(ENABLE_ECL_INPUT)
     src/opm/input/eclipse/EclipseState/Tables/TableContainer.cpp
     src/opm/input/eclipse/EclipseState/Tables/TableIndex.cpp
     src/opm/input/eclipse/EclipseState/Tables/TLMixpar.cpp
+    src/opm/input/eclipse/EclipseState/Tables/Ppcwmax.cpp
     src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp
     src/opm/input/eclipse/EclipseState/Tables/TableSchema.cpp
     src/opm/input/eclipse/EclipseState/Tables/Tables.cpp
@@ -1124,6 +1125,7 @@ if(ENABLE_ECL_INPUT)
        opm/input/eclipse/EclipseState/Tables/SgcwmisTable.hpp
        opm/input/eclipse/EclipseState/Tables/Sof2Table.hpp
        opm/input/eclipse/EclipseState/Tables/TLMixpar.hpp
+       opm/input/eclipse/EclipseState/Tables/Ppcwmax.hpp
        opm/input/eclipse/EclipseState/Tables/TableManager.hpp
        opm/input/eclipse/EclipseState/Tables/SwfnTable.hpp
        opm/input/eclipse/EclipseState/Tables/EnptvdTable.hpp

opm/input/eclipse/EclipseState/Tables/Ppcwmax.hpp

@@ -0,0 +1,79 @@
+/*
+  Copyright (C) 2023 Equinor
+
+  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/>.
+ */
+#ifndef PPCWMAX_HPP
+#define PPCWMAX_HPP
+
+#include <cstddef>
+#include <vector>
+
+namespace Opm {
+class Deck;
+
+struct PpcwmaxRecord {
+    double max_cap_pres;
+    bool option;
+
+    PpcwmaxRecord() = default;
+    PpcwmaxRecord(double pres, bool optn) :
+        max_cap_pres(pres),
+        option(optn)
+    {};
+
+    bool operator==(const PpcwmaxRecord& other) const {
+        return this->max_cap_pres == other.max_cap_pres &&
+               this->option == other.option;
+    }
+
+    template<class Serializer>
+    void serializeOp(Serializer& serializer)
+    {
+        serializer(max_cap_pres);
+        serializer(option);
+    }
+};
+
+class Ppcwmax {
+public:
+    Ppcwmax() = default;
+    explicit Ppcwmax(const Deck& deck);
+    explicit Ppcwmax(std::initializer_list<PpcwmaxRecord> records);
+    static Ppcwmax serializationTestObject();
+    std::size_t size() const;
+    bool empty() const;
+    std::vector< PpcwmaxRecord >::const_iterator begin() const;
+    std::vector< PpcwmaxRecord >::const_iterator end() const;
+    const PpcwmaxRecord& operator[](const std::size_t index) const;
+
+    bool operator==(const Ppcwmax& other) const {
+        return this->data == other.data;
+    }
+
+    template<class Serializer>
+    void serializeOp(Serializer& serializer)
+    {
+        serializer(data);
+    }
+
+private:
+    std::vector<PpcwmaxRecord> data;  
+
+};
+}
+
+#endif

opm/input/eclipse/EclipseState/Tables/TableManager.hpp

@@ -57,6 +57,7 @@
 #include <opm/input/eclipse/EclipseState/Tables/Eqldims.hpp>
 #include <opm/input/eclipse/EclipseState/Tables/Regdims.hpp>
 #include <opm/input/eclipse/EclipseState/Tables/TLMixpar.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/Ppcwmax.hpp>
 
 namespace Opm {
 
@@ -78,6 +79,7 @@ namespace Opm {
         const Aqudims& getAqudims() const;
         const Regdims& getRegdims() const;
         const TLMixpar& getTLMixpar() const;
+        const Ppcwmax& getPpcwmax() const;
         /*
           WIll return max{ Tabdims::NTFIP , Regdims::NTFIP }.
         */
@@ -260,6 +262,7 @@ namespace Opm {
             serializer(m_rtemp);
             serializer(m_salinity);
             serializer(m_tlmixpar);
+            serializer(m_ppcwmax);
             if (!serializer.isSerializing()) {
                 m_simpleTables = simpleTables;
                 if (split.plyshMax > 0) {
@@ -391,6 +394,7 @@ namespace Opm {
         Eqldims m_eqldims;
         Aqudims m_aqudims;
         TLMixpar m_tlmixpar;
+        Ppcwmax m_ppcwmax;
 
         bool hasImptvd = false;// if deck has keyword IMPTVD
         bool hasEnptvd = false;// if deck has keyword ENPTVD

opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp

@@ -252,13 +252,17 @@ public:
      * that the capillary pressure given depends on the particuars of how the simulator
      * calculates its initial condition.
      */
-    Scalar applySwatinit(unsigned elemIdx,
+    std::pair<Scalar, bool>
+    applySwatinit(unsigned elemIdx,
                          Scalar pcow,
                          Scalar Sw);
 
     bool enableEndPointScaling() const
     { return enableEndPointScaling_; }
 
+    bool enablePpcwmax() const
+    { return enablePpcwmax_; }
+
     bool enableHysteresis() const
     { return hysteresisConfig_->enableHysteresis(); }
 
@@ -410,6 +414,10 @@ private:
     std::vector<int> imbnumRegionArray_;
     std::vector<Scalar> stoneEtas_;
 
+    bool enablePpcwmax_;
+    std::vector<Scalar> maxAllowPc_;
+    std::vector<bool> modifySwl_;
+
     bool hasGas;
     bool hasOil;
     bool hasWater;

src/opm/input/eclipse/EclipseState/Tables/Ppcwmax.cpp

@@ -0,0 +1,81 @@
+/*
+  Copyright 2023 Equinor
+
+  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 <opm/input/eclipse/Deck/Deck.hpp>
+#include <opm/input/eclipse/Parser/ParserKeywords/P.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/Ppcwmax.hpp>
+
+#include <iostream>
+
+namespace Opm {
+
+Ppcwmax::Ppcwmax(const Deck& deck) {
+    using PPCWMAX = ParserKeywords::PPCWMAX;
+    if (!deck.hasKeyword<PPCWMAX>())
+        return;
+    
+    const auto& keyword = deck.get<PPCWMAX>().back();
+    for (const auto& record : keyword) {
+        // Get first column value, which is max. allowable capillary pressure
+        const double max_cap_pres = record.getItem<PPCWMAX::MAXIMUM_CAPILLARY_PRESSURE>().getSIDouble(0);
+
+        // Store second column value as bool (yes/no = true/false)
+        bool option;
+        const auto& optn = record.getItem<PPCWMAX::MODIFY_CONNATE_SATURATION>().get< std::string >(0);
+        if (optn == "YES")
+            option = true;
+        else if (optn == "NO")
+            option = false;
+        else
+            throw std::runtime_error("Second column input in PPCWMAX must be YES or NO!");
+        
+        // Store input data
+        this->data.emplace_back(max_cap_pres, option);
+    }
+}
+
+Ppcwmax::Ppcwmax(std::initializer_list<PpcwmaxRecord> records)
+    : data{ records }
+{}
+
+bool Ppcwmax::empty() const {
+    return this->data.empty();
+}
+
+std::size_t Ppcwmax::size() const {
+    return this->data.size();
+}
+
+std::vector< PpcwmaxRecord >::const_iterator Ppcwmax::begin() const {
+    return this->data.begin();
+}
+
+std::vector< PpcwmaxRecord >::const_iterator Ppcwmax::end() const {
+    return this->data.end();
+}
+
+Ppcwmax Ppcwmax::serializationTestObject() {
+    return Ppcwmax({{45.0, "YES"}});
+}
+
+
+const PpcwmaxRecord& Ppcwmax::operator[](const std::size_t index) const {
+    return this->data.at(index);
+}
+
+}

src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp

@@ -135,6 +135,7 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
         m_tabdims( Tabdims(deck)),
         m_aqudims( Aqudims(deck)),
         m_tlmixpar( deck ),
+        m_ppcwmax( deck ),
         hasImptvd (deck.hasKeyword("IMPTVD")),
         hasEnptvd (deck.hasKeyword("ENPTVD")),
         hasEqlnum (deck.hasKeyword("EQLNUM")),
@@ -322,6 +323,7 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
         result.m_rtemp = 1.0;
         result.m_salinity = 1.0;
         result.m_tlmixpar = TLMixpar::serializationTestObject();
+        result.m_ppcwmax = Ppcwmax::serializationTestObject();
         return result;
     }
 
@@ -1174,6 +1176,10 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
         return m_tlmixpar;
     }
 
+    const Ppcwmax& TableManager::getPpcwmax() const {
+        return m_ppcwmax;
+    }
+
     const JFunc& TableManager::getJFunc() const {
         if (!this->jfunc.has_value())
             throw std::invalid_argument("Cannot get JFUNC table when JFUNC not in deck");
@@ -1275,6 +1281,7 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
                m_skprwatTables == data.m_skprwatTables &&
                m_skprpolyTables == data.m_skprpolyTables &&
                m_tlmixpar == data.m_tlmixpar &&
+               m_ppcwmax == data.m_ppcwmax &&
                m_tabdims == data.m_tabdims &&
                m_regdims == data.m_regdims &&
                m_eqldims == data.m_eqldims &&

src/opm/input/eclipse/share/keywords/000_Eclipse100/P/PPCWMAX

@@ -11,6 +11,7 @@
     {
       "name": "MAXIMUM_CAPILLARY_PRESSURE",
       "value_type": "DOUBLE",
+      "dimension": "Pressure",
       "default": 1e+20
     },
     {

src/opm/material/fluidmatrixinteractions/EclMaterialLawManager.cpp

@@ -81,6 +81,22 @@ initFromState(const EclipseState& eclState)
                 stoneEtas_.push_back(table.eta);
             }
         }
+        
+        const auto& ppcwmaxTables = tables.getPpcwmax();
+        this->enablePpcwmax_ = !ppcwmaxTables.empty();
+
+        if (this->enablePpcwmax_) {
+            maxAllowPc_.clear();
+            modifySwl_.clear();
+
+            maxAllowPc_.reserve(numSatRegions);
+            modifySwl_.reserve(numSatRegions);
+
+            for (const auto& table : ppcwmaxTables) {
+                maxAllowPc_.push_back(table.max_cap_pres);
+                modifySwl_.push_back(table.option);
+            }
+        }
     }
 
     this->unscaledEpsInfo_.resize(numSatRegions);
@@ -123,12 +139,13 @@ initParamsForElements(const EclipseState& eclState, size_t numCompressedElems)
 }
 
 template<class TraitsT>
-typename TraitsT::Scalar EclMaterialLawManager<TraitsT>::
+std::pair<typename TraitsT::Scalar, bool> EclMaterialLawManager<TraitsT>::
 applySwatinit(unsigned elemIdx,
               Scalar pcow,
               Scalar Sw)
 {
     auto& elemScaledEpsInfo = oilWaterScaledEpsInfoDrainage_[elemIdx];
+    bool newSwatInit = false;
 
     // TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74
 
@@ -163,7 +180,27 @@ applySwatinit(unsigned elemIdx,
         constexpr const Scalar pcowAtSwThreshold = 1.0; //Pascal
         // avoid divison by very small number
         if (std::abs(pcowAtSw) > pcowAtSwThreshold) {
-            elemScaledEpsInfo.maxPcow *= pcow/pcowAtSw;
+            // Scale max. capillary pressure to honor SWATINIT value
+            Scalar newMaxPcow = elemScaledEpsInfo.maxPcow * (pcow/pcowAtSw);
+
+            // Limit max. capillary pressure with PPCWMAX
+            int satRegionIdx = satnumRegionIdx(elemIdx);
+            if (enablePpcwmax() && (newMaxPcow > maxAllowPc_[satRegionIdx])) {
+                // Two options in PPCWMAX to modify connate Sw or not.  In both cases, init. Sw needs to be
+                // re-calculated (done in opm-simulators)
+                newSwatInit = true;
+                if (modifySwl_[satRegionIdx] == false) {
+                    // Max. cap. pressure set to PCWO in PPCWMAX
+                    elemScaledEpsInfo.maxPcow = maxAllowPc_[satRegionIdx];
+                }
+                else {
+                    // Max. cap. pressure remains unscaled and connate Sw is set to SWATINIT value
+                    elemScaledEpsInfo.Swl = Sw;
+                }
+            }
+            // Max. cap. pressure adjusted from SWATINIT data
+            else
+                elemScaledEpsInfo.maxPcow = newMaxPcow;
             auto& elemEclEpsScalingPoints = oilWaterScaledEpsPointsDrainage(elemIdx);
             elemEclEpsScalingPoints.init(elemScaledEpsInfo,
                                          *oilWaterEclEpsConfig_,
@@ -171,7 +208,7 @@ applySwatinit(unsigned elemIdx,
         }
     }
 
-    return Sw;
+    return {Sw, newSwatInit};
 }
 
 template<class TraitsT>

tests/parser/TableManagerTests.cpp

@@ -2374,6 +2374,29 @@ BOOST_AUTO_TEST_CASE( TestParseDIFFCWATGAS ) {
 }
 
 
+BOOST_AUTO_TEST_CASE(TestParsePPCWMAX) {
+    const std::string data = R"(
+        TABDIMS
+        2 /
+
+        PPCWMAX
+        10.0 /
+        1* YES/
+    )";
+
+    Opm::Parser parser;
+    auto deck = parser.parseString(data);
+    Opm::TableManager tables(deck);
+
+    const auto& ppcwmax = tables.getPpcwmax();
+    BOOST_CHECK_CLOSE(10.0e5, ppcwmax[0].max_cap_pres, epsilon());
+    BOOST_CHECK_EQUAL(false, ppcwmax[0].option);
+
+    BOOST_CHECK_CLOSE(1e+25, ppcwmax[1].max_cap_pres, epsilon());
+    BOOST_CHECK_EQUAL(true, ppcwmax[1].option);
+}
+
+
 BOOST_AUTO_TEST_CASE( TestParseROCK ) {
     const std::string data = R"(
       TABDIMS

tests/test_Serialization.cpp

@@ -227,6 +227,7 @@ TEST_FOR_TYPE(EndpointScaling)
 TEST_FOR_TYPE(Eqldims)
 TEST_FOR_TYPE(Equil)
 TEST_FOR_TYPE(TLMixpar)
+TEST_FOR_TYPE(Ppcwmax)
 TEST_FOR_TYPE(Events)
 TEST_FOR_TYPE(Fault)
 TEST_FOR_TYPE(FaultCollection)