Merge pull request #375 from dr-robertk/PR/add-keyword-SOF2-SWFN

Added keywords SOF2 and SWFN to allow for different saturation functions.

opm/parser/eclipse/CMakeLists.txt

@@ -173,6 +173,8 @@ EclipseState/Tables/PlymaxTable.hpp
 EclipseState/Tables/PvtgTable.hpp
 EclipseState/Tables/PlyrockTable.hpp
 EclipseState/Tables/SwofTable.hpp
+EclipseState/Tables/SwfnTable.hpp
+EclipseState/Tables/Sof2Table.hpp
 EclipseState/Tables/EnptvdTable.hpp
 EclipseState/Tables/FullTable.hpp
 EclipseState/Tables/PlyviscTable.hpp

opm/parser/eclipse/EclipseState/EclipseState.cpp

@@ -39,23 +39,23 @@ namespace Opm {
         class DistributeTopLayer : public GridPropertyBasePostProcessor<double>
         {
         public:
-            DistributeTopLayer(const EclipseState& eclipseState) : 
+            DistributeTopLayer(const EclipseState& eclipseState) :
                 m_eclipseState( eclipseState )
             { }
-            
+
-            
+
             void apply(std::vector<double>& values) const {
                 EclipseGridConstPtr grid = m_eclipseState.getEclipseGrid();
                 size_t layerSize = grid->getNX() * grid->getNY();
                 size_t gridSize  = grid->getCartesianSize();
-                
+
                 for (size_t globalIndex = layerSize; globalIndex < gridSize; globalIndex++) {
                     if (std::isnan( values[ globalIndex ] ))
                         values[globalIndex] = values[globalIndex - layerSize];
                 }
             }
-            
+
-            
+
         private:
             const EclipseState& m_eclipseState;
         };
@@ -65,11 +65,11 @@ namespace Opm {
         class InitPORV : public GridPropertyBasePostProcessor<double>
         {
         public:
-            InitPORV(const EclipseState& eclipseState) : 
+            InitPORV(const EclipseState& eclipseState) :
                 m_eclipseState( eclipseState )
             { }
-            
+
-            
+
             void apply(std::vector<double>& ) const {
                 EclipseGridConstPtr grid = m_eclipseState.getEclipseGrid();
                 /*
@@ -77,10 +77,10 @@ namespace Opm {
                   values input vector, instead it fetches the PORV
                   property one more time, and then manipulates that.
                 */
-                auto porv = m_eclipseState.getDoubleGridProperty("PORV");   
+                auto porv = m_eclipseState.getDoubleGridProperty("PORV");
                 if (porv->containsNaN()) {
                     auto poro = m_eclipseState.getDoubleGridProperty("PORO");
-                    auto ntg = m_eclipseState.getDoubleGridProperty("NTG");   
+                    auto ntg = m_eclipseState.getDoubleGridProperty("NTG");
                     if (poro->containsNaN())
                         throw std::logic_error("Do not have information for the PORV keyword - some defaulted values in PORO");
                     {
@@ -92,25 +92,25 @@ namespace Opm {
                                 porv->iset( globalIndex , cell_poro * cell_volume * cell_ntg);
                             }
                         }
-                    } 
+                    }
                 }
-                
+
                 if (m_eclipseState.hasDoubleGridProperty("MULTPV")) {
-                    auto multpv = m_eclipseState.getDoubleGridProperty("MULTPV");   
+                    auto multpv = m_eclipseState.getDoubleGridProperty("MULTPV");
                     porv->multiplyWith( *multpv );
                 }
             }
-            
+
-            
+
         private:
             const EclipseState& m_eclipseState;
         };
 
-        
+
     }
 
 
-    EclipseState::EclipseState(DeckConstPtr deck, ParserLogPtr parserLog)    
+    EclipseState::EclipseState(DeckConstPtr deck, ParserLogPtr parserLog)
     {
         m_deckUnitSystem = deck->getActiveUnitSystem();
 
@@ -208,10 +208,18 @@ namespace Opm {
         return m_sgofTables;
     }
 
+    const std::vector<Sof2Table>& EclipseState::getSof2Tables() const {
+        return m_sof2Tables;
+    }
+
     const std::vector<SwofTable>& EclipseState::getSwofTables() const {
         return m_swofTables;
     }
 
+    const std::vector<SwfnTable>& EclipseState::getSwfnTables() const {
+        return m_swfnTables;
+    }
+
     ScheduleConstPtr EclipseState::getSchedule() const {
         return schedule;
     }
@@ -242,7 +250,9 @@ namespace Opm {
         initSimpleTables(deck, parserLog, "RSVD", m_rsvdTables);
         initSimpleTables(deck, parserLog, "RVVD", m_rvvdTables);
         initSimpleTables(deck, parserLog, "SGOF", m_sgofTables);
+        initSimpleTables(deck, parserLog, "SOF2", m_sof2Tables);
         initSimpleTables(deck, parserLog, "SWOF", m_swofTables);
+        initSimpleTables(deck, parserLog, "SWFN", m_swfnTables);
 
         // the ROCKTAB table comes with additional fun because the number of columns
         //depends on the presence of the RKTRMDIR keyword...
@@ -305,8 +315,8 @@ namespace Opm {
                 int K2 = faultRecord->getItem(6)->getInt(0) - 1;
                 FaceDir::DirEnum faceDir = FaceDir::FromString( faultRecord->getItem(7)->getString(0) );
                 std::shared_ptr<const FaultFace> face = std::make_shared<const FaultFace>(grid->getNX() , grid->getNY() , grid->getNZ(),
-                                                                                          static_cast<size_t>(I1) , static_cast<size_t>(I2) , 
+                                                                                          static_cast<size_t>(I1) , static_cast<size_t>(I2) ,
-                                                                                          static_cast<size_t>(J1) , static_cast<size_t>(J2) , 
+                                                                                          static_cast<size_t>(J1) , static_cast<size_t>(J2) ,
                                                                                           static_cast<size_t>(K1) , static_cast<size_t>(K2) ,
                                                                                           faceDir);
                 if (!m_faults->hasFault(faultName)) {
@@ -320,7 +330,7 @@ namespace Opm {
                 }
             }
         }
-        
+
         setMULTFLT(gridSection, parserLog);
 
         if (Section::hasEDIT(deck)) {
@@ -340,14 +350,14 @@ namespace Opm {
                 DeckRecordConstPtr faultRecord = *iter;
                 const std::string& faultName = faultRecord->getItem(0)->getString(0);
                 double multFlt = faultRecord->getItem(1)->getRawDouble(0);
-                
+
                 m_faults->setTransMult( faultName , multFlt );
             }
         }
     }
 
 
-    
+
     void EclipseState::initMULTREGT(DeckConstPtr deck, ParserLogPtr /*parserLog*/) {
         EclipseGridConstPtr grid = getEclipseGrid();
         std::shared_ptr<MULTREGTScanner> scanner = std::make_shared<MULTREGTScanner>();
@@ -358,7 +368,7 @@ namespace Opm {
                 DeckKeywordConstPtr multregtKeyword = gridSection->getKeyword("MULTREGT" , index);
                 scanner->addKeyword( multregtKeyword );
             }
-        }            
+        }
 
 
         if (Section::hasEDIT(deck)) {
@@ -371,7 +381,7 @@ namespace Opm {
 
         m_transMult->applyMULTREGT( scanner , m_intGridProperties);
     }
-    
+
 
 
     void EclipseState::initEclipseGrid(DeckConstPtr deck, ParserLogPtr parserLog) {
@@ -456,17 +466,17 @@ namespace Opm {
 
     bool EclipseState::hasIntGridProperty(const std::string& keyword) const {
          return m_intGridProperties->hasKeyword( keyword );
-    }   
+    }
 
     bool EclipseState::hasDoubleGridProperty(const std::string& keyword) const {
          return m_doubleGridProperties->hasKeyword( keyword );
-    }   
+    }
 
 
     /*
       1. The public methods getIntGridProperty & getDoubleGridProperty
          will invoke and run the property post processor (if any is
-         registered); the post processor will only run one time. 
+         registered); the post processor will only run one time.
 
          It is important that post processor is not run prematurely,
          internal functions in EclipseState should therefor ask for
@@ -520,8 +530,8 @@ namespace Opm {
                                 "Tried to load unsupported grid property from keyword: " + deckKeyword->name());
         }
     }
-        
+
-    
+
     void EclipseState::initProperties(DeckConstPtr deck, ParserLogPtr parserLog) {
         typedef GridProperties<int>::SupportedKeywordInfo SupportedIntKeywordInfo;
         std::shared_ptr<std::vector<SupportedIntKeywordInfo> > supportedIntKeywords(new std::vector<SupportedIntKeywordInfo>{
@@ -734,7 +744,7 @@ namespace Opm {
     }
 
     void EclipseState::processGridProperties(Opm::DeckConstPtr deck, ParserLogPtr parserLog, int enabledTypes) {
-        
+
         if (Section::hasGRID(deck)) {
             std::shared_ptr<Opm::GRIDSection> gridSection(new Opm::GRIDSection(deck) );
             scanSection(gridSection, parserLog, enabledTypes);
@@ -774,22 +784,22 @@ namespace Opm {
             else {
                 if (deckKeyword->name() == "ADD")
                     handleADDKeyword(deckKeyword, parserLog, boxManager, enabledTypes);
-            
+
                 if (deckKeyword->name() == "BOX")
                     handleBOXKeyword(deckKeyword, parserLog, boxManager);
-            
+
                 if (deckKeyword->name() == "COPY")
                     handleCOPYKeyword(deckKeyword, parserLog, boxManager, enabledTypes);
-            
+
                 if (deckKeyword->name() == "EQUALS")
                     handleEQUALSKeyword(deckKeyword, parserLog, boxManager, enabledTypes);
-            
+
                 if (deckKeyword->name() == "ENDBOX")
                     handleENDBOXKeyword(boxManager);
-            
+
                 if (deckKeyword->name() == "MULTIPLY")
                     handleMULTIPLYKeyword(deckKeyword, parserLog, boxManager, enabledTypes);
-            
+
                 boxManager.endKeyword();
             }
         }
@@ -807,7 +817,7 @@ namespace Opm {
         int J2 = record->getItem("J2")->getInt(0) - 1;
         int K1 = record->getItem("K1")->getInt(0) - 1;
         int K2 = record->getItem("K2")->getInt(0) - 1;
-        
+
         boxManager.setInputBox( I1 , I2 , J1 , J2 , K1 , K2 );
     }
 
@@ -822,9 +832,9 @@ namespace Opm {
             DeckRecordConstPtr record = deckKeyword->getRecord(recordIdx);
             const std::string& field = record->getItem("field")->getString(0);
             double      scaleFactor  = record->getItem("factor")->getRawDouble(0);
-            
+
             setKeywordBox(deckKeyword, recordIdx, parserLog, boxManager);
-            
+
             if (m_intGridProperties->hasKeyword( field )) {
                 if (enabledTypes & IntProperties) {
                     int intFactor = static_cast<int>(scaleFactor);
@@ -847,16 +857,16 @@ namespace Opm {
     /*
       The fine print of the manual says the ADD keyword should support
       some state dependent semantics regarding endpoint scaling arrays
-      in the PROPS section. That is not supported. 
+      in the PROPS section. That is not supported.
     */
     void EclipseState::handleADDKeyword(DeckKeywordConstPtr deckKeyword, ParserLogPtr parserLog, BoxManager& boxManager, int enabledTypes) {
         for (size_t recordIdx = 0; recordIdx < deckKeyword->size(); ++recordIdx) {
             DeckRecordConstPtr record = deckKeyword->getRecord(recordIdx);
             const std::string& field = record->getItem("field")->getString(0);
             double      shiftValue  = record->getItem("shift")->getRawDouble(0);
-            
+
             setKeywordBox(deckKeyword, recordIdx, parserLog, boxManager);
-            
+
             if (m_intGridProperties->hasKeyword( field )) {
                 if (enabledTypes & IntProperties) {
                     int intShift = static_cast<int>(shiftValue);
@@ -884,9 +894,9 @@ namespace Opm {
             DeckRecordConstPtr record = deckKeyword->getRecord(recordIdx);
             const std::string& field = record->getItem("field")->getString(0);
             double      value  = record->getItem("value")->getRawDouble(0);
-            
+
             setKeywordBox(deckKeyword, recordIdx, parserLog, boxManager);
-            
+
             if (m_intGridProperties->supportsKeyword( field )) {
                 if (enabledTypes & IntProperties) {
                     int intValue = static_cast<int>(value);
@@ -915,9 +925,9 @@ namespace Opm {
             DeckRecordConstPtr record = deckKeyword->getRecord(recordIdx);
             const std::string& srcField = record->getItem("src")->getString(0);
             const std::string& targetField = record->getItem("target")->getString(0);
-            
+
             setKeywordBox(deckKeyword, recordIdx, parserLog, boxManager);
-            
+
             if (m_intGridProperties->hasKeyword( srcField )) {
                 if (enabledTypes & IntProperties)
                     copyIntKeyword( srcField , targetField , boxManager.getActiveBox());
@@ -932,7 +942,7 @@ namespace Opm {
         }
     }
 
-    
+
     void EclipseState::copyIntKeyword(const std::string& srcField , const std::string& targetField , std::shared_ptr<const Box> inputBox) {
         std::shared_ptr<const GridProperty<int> > src = m_intGridProperties->getKeyword( srcField );
         std::shared_ptr<GridProperty<int> > target    = m_intGridProperties->getKeyword( targetField );
@@ -961,7 +971,7 @@ namespace Opm {
         DeckItemConstPtr K2Item = deckRecord->getItem("K2");
 
         size_t setCount = 0;
-        
+
         if (!I1Item->defaultApplied(0))
             setCount++;
 
@@ -979,7 +989,7 @@ namespace Opm {
 
         if (!K2Item->defaultApplied(0))
             setCount++;
-        
+
         if (setCount == 6) {
             boxManager.setKeywordBox( I1Item->getInt(0) - 1,
                                       I2Item->getInt(0) - 1,

opm/parser/eclipse/EclipseState/EclipseState.hpp

@@ -49,7 +49,9 @@
 #include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SgofTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/SwfnTable.hpp>
 
 #include <set>
 #include <memory>
@@ -107,7 +109,9 @@ namespace Opm {
         const std::vector<RvvdTable>& getRvvdTables() const;
         const std::vector<RtempvdTable>& getRtempvdTables() const;
         const std::vector<SgofTable>& getSgofTables() const;
+        const std::vector<Sof2Table>& getSof2Tables() const;
         const std::vector<SwofTable>& getSwofTables() const;
+        const std::vector<SwfnTable>& getSwfnTables() const;
 
         // the unit system used by the deck. note that it is rarely needed to convert
         // units because internally to opm-parser everything is represented by SI
@@ -195,7 +199,7 @@ namespace Opm {
         void handleENDBOXKeyword(BoxManager& boxManager);
         void handleEQUALSKeyword(DeckKeywordConstPtr deckKeyword , ParserLogPtr parserLog, BoxManager& boxManager, int enabledTypes);
         void handleMULTIPLYKeyword(DeckKeywordConstPtr deckKeyword , ParserLogPtr parserLog, BoxManager& boxManager, int enabledTypes);
-        
+
         void setKeywordBox(DeckKeywordConstPtr deckKeyword, size_t recordIdx, ParserLogPtr parserLog, BoxManager& boxManager);
 
         void copyIntKeyword(const std::string& srcField , const std::string& targetField , std::shared_ptr<const Box> inputBox);
@@ -223,7 +227,9 @@ namespace Opm {
         std::vector<RvvdTable> m_rvvdTables;
         std::vector<RtempvdTable> m_rtempvdTables;
         std::vector<SgofTable> m_sgofTables;
+        std::vector<Sof2Table> m_sof2Tables;
         std::vector<SwofTable> m_swofTables;
+        std::vector<SwfnTable> m_swfnTables;
 
         std::set<enum Phase::PhaseEnum> phases;
         std::string m_title;

opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp

@@ -0,0 +1,74 @@
+/*
+  Copyright (C) 2012 IRIS AS
+
+  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 OPM_PARSER_SOF2_TABLE_HPP
+#define OPM_PARSER_SOF2_TABLE_HPP
+
+#include "SingleRecordTable.hpp"
+
+namespace Opm {
+    // forward declaration
+    class EclipseState;
+
+    class Sof2Table : protected SingleRecordTable {
+        typedef SingleRecordTable ParentType;
+
+        friend class EclipseState;
+
+        /*!
+         * \brief Read the SOF2 keyword and provide some convenience
+         *        methods for it.
+         */
+        void init(Opm::DeckKeywordConstPtr keyword,
+                  int recordIdx)
+        {
+            ParentType::init(keyword,
+                             std::vector<std::string>{"SO", "KRO" },
+                             recordIdx,
+                             /*firstEntityOffset=*/0);
+
+            ParentType::checkNonDefaultable("SO");
+            ParentType::checkNonDefaultable("KRO");
+            ParentType::checkMonotonic("SO", /*isAscending=*/true);
+            ParentType::checkMonotonic("KRO", /*isAscending=*/true, /*strict*/false);
+        }
+
+    public:
+        Sof2Table() = default;
+
+#ifdef BOOST_TEST_MODULE
+        // DO NOT TRY TO CALL THIS METHOD! it is only for the unit tests!
+        void initFORUNITTESTONLY(Opm::DeckKeywordConstPtr keyword, size_t tableIdx)
+        { init(keyword, tableIdx); }
+#endif
+
+        using ParentType::numTables;
+        using ParentType::numRows;
+        using ParentType::numColumns;
+        using ParentType::evaluate;
+
+        const std::vector<double> &getSoColumn() const
+        { return ParentType::getColumn(0); }
+
+        const std::vector<double> &getKroColumn() const
+        { return ParentType::getColumn(1); }
+    };
+}
+
+#endif
+

opm/parser/eclipse/EclipseState/Tables/SwfnTable.hpp

@@ -0,0 +1,81 @@
+/*
+  Copyright (C) 2014 IRIS AS
+
+  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 OPM_PARSER_SWFN_TABLE_HPP
+#define OPM_PARSER_SWFN_TABLE_HPP
+
+#include "SingleRecordTable.hpp"
+
+namespace Opm {
+    // forward declaration
+    class EclipseState;
+
+    class SwfnTable : protected SingleRecordTable {
+        typedef SingleRecordTable ParentType;
+
+        friend class EclipseState;
+
+        /*!
+         * \brief Read the SWFN keyword and provide some convenience
+         *        methods for it.
+         */
+        void init(Opm::DeckKeywordConstPtr keyword,
+                  int recordIdx)
+        {
+            ParentType::init(keyword,
+                             std::vector<std::string>{"SW", "KRW", "PCOW"},
+                             recordIdx,
+                             /*firstEntityOffset=*/0);
+
+            ParentType::checkNonDefaultable("SW");
+            ParentType::checkMonotonic("SW",   /*isAscending=*/true);
+            ParentType::applyDefaultsLinear("KRW");
+            ParentType::applyDefaultsLinear("PCOW");
+            ParentType::checkMonotonic("KRW",  /*isAscending=*/true,  /*strict=*/false);
+            ParentType::checkMonotonic("PCOW", /*isAscending=*/false, /*strict=*/false);
+        }
+
+    public:
+        SwfnTable() = default;
+
+#ifdef BOOST_TEST_MODULE
+        // DO NOT TRY TO CALL THIS METHOD! it is only for the unit tests!
+        void initFORUNITTESTONLY(Opm::DeckKeywordConstPtr keyword, size_t tableIdx)
+        { init(keyword, tableIdx); }
+#endif
+
+        using ParentType::numTables;
+        using ParentType::numRows;
+        using ParentType::numColumns;
+        using ParentType::evaluate;
+
+        const std::vector<double> &getSwColumn() const
+        { return ParentType::getColumn(0); }
+
+        const std::vector<double> &getKrwColumn() const
+        { return ParentType::getColumn(1); }
+
+        // this column is p_o - p_w (non-wetting phase pressure minus
+        // wetting phase pressure for a given water saturation)
+        const std::vector<double> &getPcowColumn() const
+        { return ParentType::getColumn(2); }
+    };
+}
+
+#endif
+

opm/parser/share/keywords/000_Eclipse100/S/SOF2

@@ -0,0 +1,5 @@
+{"name" : "SOF2" , "sections" : ["PROPS"], "size" : {"keyword" : "TABDIMS" , "item" : "NTSFUN"},
+        "items" : [
+            {"name":"DATA", "value_type":"DOUBLE", "size_type" : "ALL" , "dimension" : ["1","1"]}
+         ]
+}

opm/parser/share/keywords/000_Eclipse100/S/SWFN

@@ -0,0 +1,5 @@
+{"name" : "SWFN" , "sections" : ["PROPS"], "size" : {"keyword" : "TABDIMS" , "item" : "NTSFUN"},
+        "items" : [
+            {"name":"DATA", "value_type":"DOUBLE", "size_type" : "ALL" , "dimension" : ["1","1","Pressure"]}
+         ]
+}