Merge pull request #3225 from akva2/internalize_thpresft

Internalize THPRESFT keyword

opm/input/eclipse/EclipseState/SimulationConfig/BCConfig.hpp

@@ -22,6 +22,7 @@
 
 #include <vector>
 #include <cstddef>
+#include <optional>
 
 #include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/GridDims.hpp>

opm/input/eclipse/EclipseState/SimulationConfig/RockConfig.hpp

@@ -21,6 +21,7 @@
 #define OPM_ROCK_CONFIG_HPP
 
 #include <string>
+#include <vector>
 
 namespace Opm {
 

opm/input/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp

@@ -71,6 +71,7 @@ namespace Opm {
         }
 
     private:
+        friend class EclipseState;
         ThresholdPressure m_ThresholdPressure;
         BCConfig m_bcconfig;
         RockConfig m_rock_config;

opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp

@@ -20,6 +20,7 @@
 #ifndef OPM_TRESHOLD_PRESSURES_HPP
 #define OPM_TRESHOLD_PRESSURES_HPP
 
+#include <cstddef>
 #include <map>
 #include <vector>
 
@@ -27,6 +28,7 @@ namespace Opm {
 
 
     class Deck;
+    class FaultCollection;
     class FieldPropsManager;
 
     class ThresholdPressure {
@@ -45,6 +47,11 @@ namespace Opm {
             , m_irreversible(false)
         {}
 
+
+        //! \brief Reads the THPRESFT keyword if present.
+        void readFaults(const Deck& deck,
+                        const FaultCollection& faults);
+
         //! \brief Returns an instance for serialization tests.
         static ThresholdPressure serializationTestObject();
 
@@ -73,6 +80,11 @@ namespace Opm {
           hasThresholdPressure(r1,r2) first to be safe.
         */
         double getThresholdPressure(int r1 , int r2) const;
+
+        //! \brief Returns threshold pressure for a fault.
+        double getThresholdPressureFault(int idx) const;
+
+        size_t ftSize() const;
         size_t size() const;
         bool active() const;
         bool restart() const;
@@ -89,6 +101,7 @@ namespace Opm {
             serializer(m_irreversible);
             serializer(m_thresholdPressureTable);
             serializer(m_pressureTable);
+            serializer(m_thresholdFaultTable);
         }
 
     private:
@@ -102,6 +115,7 @@ namespace Opm {
 
         std::vector<std::pair<bool,double>> m_thresholdPressureTable;
         std::map<std::pair<int,int> , std::pair<bool , double> > m_pressureTable;
+        std::vector<double> m_thresholdFaultTable;
     };
 } //namespace Opm
 

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

@@ -143,6 +143,7 @@ namespace Opm {
 
         this->applyMULTXYZ();
         this->initFaults(deck);
+        m_simulationConfig.m_ThresholdPressure.readFaults(deck,m_faults);
 
         if (this->getInitConfig().restartRequested()) {
             verify_consistent_restart_information(deck.get<ParserKeywords::RESTART>().back(),

src/opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.cpp

@@ -18,7 +18,9 @@
  */
 
 #include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/utility/shmatch.hpp>
 #include <opm/input/eclipse/Deck/DeckSection.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaultCollection.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
 #include <opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
 #include <opm/input/eclipse/Parser/ParserKeywords/E.hpp>
@@ -127,6 +129,38 @@ namespace Opm {
         }
     }
 
+    void ThresholdPressure::readFaults(const Deck& deck,
+                                       const FaultCollection& faults)
+    {
+        if (!DeckSection::hasGRID(deck) || !DeckSection::hasSOLUTION(deck))
+            return;
+
+        GRIDSection gridSection( deck );
+
+        const bool thpresftKeyword  = gridSection.hasKeyword<ParserKeywords::THPRESFT>();
+        if (!thpresftKeyword)
+            return;
+
+        // extract the multipliers from the deck keyword
+        m_thresholdFaultTable.resize(faults.size(), -1.0);
+        for (const auto& thpresft : gridSection.get<ParserKeywords::THPRESFT>()) {
+            for (size_t recordIdx = 0; recordIdx < thpresft.size(); ++ recordIdx) {
+                const DeckRecord& record = thpresft.getRecord(recordIdx);
+
+                const std::string& faultName = record.getItem("FAULT_NAME").getTrimmedString(0);
+                double thpresValue = record.getItem("VALUE").getSIDouble(0);
+
+                for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) {
+                    auto& fault = faults.getFault(faultIdx);
+                    if (!shmatch(faultName, fault.getName()))
+                        continue;
+
+                    m_thresholdFaultTable[faultIdx] = thpresValue;
+                }
+            }
+        }
+    }
+
     ThresholdPressure ThresholdPressure::serializationTestObject()
     {
         ThresholdPressure result;
@@ -166,9 +200,11 @@ namespace Opm {
             }
         }
 
-
     }
 
+    double ThresholdPressure::getThresholdPressureFault(int idx) const {
+        return m_thresholdFaultTable[idx];
+    }
 
     std::pair<int,int> ThresholdPressure::makeIndex(int r1 , int r2) const {
         if (this->m_irreversible)
@@ -199,6 +235,10 @@ namespace Opm {
         return m_pressureTable.size();
     }
 
+    size_t ThresholdPressure::ftSize() const {
+        return m_thresholdFaultTable.size();
+    }
+
     bool ThresholdPressure::active() const {
         return m_active;
     }

tests/parser/ThresholdPressureTest.cpp

@@ -26,7 +26,9 @@
 #include <boost/test/unit_test.hpp>
 
 #include <opm/input/eclipse/Deck/Deck.hpp>
+#include <opm/input/eclipse/Deck/DeckSection.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaultCollection.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
 #include <opm/input/eclipse/EclipseState/InitConfig/InitConfig.hpp>
 #include <opm/input/eclipse/EclipseState/Runspec.hpp>
@@ -54,6 +56,39 @@ THPRES
 /
 )";
 
+const std::string& inputStrft = R"(
+RUNSPEC
+EQLOPTS
+THPRES /
+DIMENS
+10 10 10 /
+
+REGIONS
+EQLNUM
+   120*3 /
+SOLUTION
+THPRES
+1 2 12.0/
+1 3 5.0/
+2 3 33.0 /
+2 3 7.0/
+/
+
+GRID
+DX
+1000*0.25 /
+DY
+1000*0.25 /
+DZ
+1000*0.25 /
+FAULTS
+  'F1'  1  1  1  4   1  4  'X' /
+/
+THPRESFT
+  'F1' 1.0 /
+/
+)";
+
 const std::string& inputStr_RESTART = R"(
 RUNSPEC
 EQLOPTS
@@ -235,7 +270,7 @@ struct Setup
     InitConfig initConfig;
     ThresholdPressure threshPres;
 
-    explicit Setup(const std::string& input) :
+    explicit Setup(const std::string& input, bool ft = false) :
             deck(createDeck(ParseContext(), input)),
             tablemanager(deck),
             grid(10, 3, 4),
@@ -243,6 +278,12 @@ struct Setup
             initConfig(deck),
             threshPres(initConfig.restartRequested(), deck, fp)
     {
+        if (ft) {
+            GRIDSection gridSection(deck);
+            GridDims dims(deck);
+            auto faults = FaultCollection(gridSection, dims);
+            threshPres.readFaults(deck, faults);
+        }
     }
 
     explicit Setup(const std::string& input, const ParseContext& parseContextArg) :
@@ -341,3 +382,9 @@ BOOST_AUTO_TEST_CASE(Irreversible) {
     const auto& thp = s.threshPres;
     BOOST_CHECK(thp.getThresholdPressure(1,2) == -thp.getThresholdPressure(2,1));
 }
+
+BOOST_AUTO_TEST_CASE(Faults) {
+    Setup s(inputStrft, true);
+    const auto& thp = s.threshPres;
+    BOOST_CHECK(thp.getThresholdPressureFault(0) == 100000.0);
+}