Update release version to final 2019.10

improvements related to Aquifer modeling (Backport of pr #1068)

dune.module

@@ -5,8 +5,8 @@
 
 Module: opm-common
 Description: Open Porous Media Initiative shared infrastructure
-Version: 2019.10-rc2
-Label: 2019.10-rc2
+Version: 2019.10
+Label: 2019.10
 Maintainer: opm@opm-project.org
 MaintainerName: OPM community
 Url: http://opm-project.org

opm/parser/eclipse/EclipseState/Aquancon.hpp

@@ -56,14 +56,20 @@ namespace Opm {
 
         private:
 
-            std::vector<Aquancon::AquanconOutput> logic_application(const std::vector<Aquancon::AquanconOutput>& original_vector);
+            static std::vector<Aquancon::AquanconOutput> logic_application(const std::vector<Aquancon::AquanconOutput>& original_vector);
 
-            void collate_function(std::vector<Aquancon::AquanconOutput>& output_vector,
+            static void collate_function(std::vector<Aquancon::AquanconOutput>& output_vector,
                                   std::vector<Opm::AquanconRecord>& m_aqurecord,
-                                  std::vector<int> m_aquiferID_per_record, int m_maxAquID);
+                                  const std::vector<int>& m_aquiferID_per_record, int m_maxAquID);
 
-            void convert_record_id_to_aquifer_id(std::vector<int>& record_indices_matching_id, int i,
-                                                 std::vector<int> m_aquiferID_per_record);
+            static void convert_record_id_to_aquifer_id(std::vector<int>& record_indices_matching_id, int i,
+                                                        const std::vector<int>& m_aquiferID_per_record);
+
+            // for a cell to be inside reservoir, its indices need to be within the reservoir grid dimension range,
+            // and it needs to be active
+            static bool cellInsideReservoirAndActive(const EclipseGrid& grid, int i, int j, int k);
+
+            static bool neighborCellInsideReservoirAndActive(const EclipseGrid& grid, int i, int j, int k, FaceDir::DirEnum faceDir);
 
             std::vector<Aquancon::AquanconOutput> m_aquoutput;
     };

redhat/opm-common.spec

@@ -2,7 +2,7 @@
 # spec file for package opm-common
 #
 
-%define tag rc2
+%define tag final
 
 Name:           opm-common
 Version:        2019.10

src/opm/parser/eclipse/EclipseState/Aquancon.cpp

@@ -26,9 +26,6 @@
 namespace Opm {
     namespace{
         struct AquanconRecord{
-            // Grid cell box definition to connect aquifer
-            int i1, i2, j1, j2, k1, k2;
-
             std::vector<size_t> global_index_per_record;
 
             // Variables constants
@@ -55,54 +52,63 @@ namespace Opm {
         aqurecords.resize(aquanconKeyword.size());
         aquiferID_per_record.resize(aquanconKeyword.size());
 
-        // We now do a loop over each record entry in aquancon
         for (size_t aquanconRecordIdx = 0; aquanconRecordIdx < aquanconKeyword.size(); ++aquanconRecordIdx)
         {
             const auto& aquanconRecord = aquanconKeyword.getRecord(aquanconRecordIdx);
-            aquiferID_per_record.at(aquanconRecordIdx) = aquanconRecord.getItem("AQUIFER_ID").template get<int>(0);
+            aquiferID_per_record[aquanconRecordIdx] = aquanconRecord.getItem("AQUIFER_ID").template get<int>(0);
 
-            aqurecords.at(aquanconRecordIdx).i1 = aquanconRecord.getItem("I1").template get<int>(0);
-            aqurecords.at(aquanconRecordIdx).i2 = aquanconRecord.getItem("I2").template get<int>(0);
-            aqurecords.at(aquanconRecordIdx).j1 = aquanconRecord.getItem("J1").template get<int>(0);
-            aqurecords.at(aquanconRecordIdx).j2 = aquanconRecord.getItem("J2").template get<int>(0);
-            aqurecords.at(aquanconRecordIdx).k1 = aquanconRecord.getItem("K1").template get<int>(0);
-            aqurecords.at(aquanconRecordIdx).k2 = aquanconRecord.getItem("K2").template get<int>(0);
+            // offset the indices
+            const int i1 = aquanconRecord.getItem("I1").template get<int>(0) - 1;
+            const int i2 = aquanconRecord.getItem("I2").template get<int>(0) - 1;
+            const int j1 = aquanconRecord.getItem("J1").template get<int>(0) - 1;
+            const int j2 = aquanconRecord.getItem("J2").template get<int>(0) - 1;
+            const int k1 = aquanconRecord.getItem("K1").template get<int>(0) - 1;
+            const int k2 = aquanconRecord.getItem("K2").template get<int>(0) - 1;
 
-            aquiferID_per_record.at(aquanconRecordIdx) = aquanconRecord.getItem("AQUIFER_ID").template get<int>(0);
-            m_maxAquID = (m_maxAquID < aquiferID_per_record.at(aquanconRecordIdx) )?
-                            aquiferID_per_record.at(aquanconRecordIdx) : m_maxAquID;
+            m_maxAquID = (m_maxAquID < aquiferID_per_record[aquanconRecordIdx] )?
+                            aquiferID_per_record[aquanconRecordIdx] : m_maxAquID;
 
-            double influx_mult = aquanconRecord.getItem("INFLUX_MULT").getSIDouble(0);
+            const FaceDir::DirEnum faceDir = FaceDir::FromString(aquanconRecord.getItem("FACE").getTrimmedString(0));
 
-            FaceDir::DirEnum faceDir = FaceDir::FromString(aquanconRecord.getItem("FACE").getTrimmedString(0));
+            // whether allow aquifer connections to locate inside reservoir
+            const std::string& str_inside_reservoir  = aquanconRecord.getItem("CONNECT_ADJOINING_ACTIVE_CELL").getTrimmedString(0);
+            // not sure whether we should give a warning when input other than "YES" or "NO"
+            const bool allow_aquifer_inside_reservoir = str_inside_reservoir == "YES" ? true : false;
 
+            auto& aqurecord = aqurecords[aquanconRecordIdx];
             // Loop over the cartesian indices to convert to the global grid index
-            for (int k=aqurecords.at(aquanconRecordIdx).k1; k <= aqurecords.at(aquanconRecordIdx).k2; k++) {
-                for (int j=aqurecords.at(aquanconRecordIdx).j1; j <= aqurecords.at(aquanconRecordIdx).j2; j++)
-                    for (int i=aqurecords.at(aquanconRecordIdx).i1; i <= aqurecords.at(aquanconRecordIdx).i2; i++)
-                        aqurecords.at(aquanconRecordIdx).global_index_per_record.push_back
-                                                            (
-                                                                grid.getGlobalIndex(i-1, j-1, k-1)
-                                                            );
+            for (int k = k1; k <= k2; k++) {
+                for (int j = j1; j <= j2; j++) {
+                    for (int i = i1; i <= i2; i++) {
+                        if ( grid.cellActive(i, j, k) ) { // the cell itself needs to be active
+                            if (  allow_aquifer_inside_reservoir
+                               || !neighborCellInsideReservoirAndActive(grid, i, j, k, faceDir) ) {
+                                     aqurecord.global_index_per_record.push_back(grid.getGlobalIndex(i, j, k));
+                             }
+                        }
+                    }
+                }
             }
-            size_t global_index_per_record_size = aqurecords.at(aquanconRecordIdx).global_index_per_record.size();
+            const size_t global_index_per_record_size = aqurecord.global_index_per_record.size();
 
-            aqurecords.at(aquanconRecordIdx).influx_coeff_per_record.resize(global_index_per_record_size, nullptr);
+            aqurecord.influx_coeff_per_record.resize(global_index_per_record_size, nullptr);
 
             if (aquanconRecord.getItem("INFLUX_COEFF").hasValue(0))
             {
                 const double influx_coeff = aquanconRecord.getItem("INFLUX_COEFF").getSIDouble(0);
 
-                for (auto& influx: aqurecords.at(aquanconRecordIdx).influx_coeff_per_record)
+                for (auto& influx: aqurecord.influx_coeff_per_record)
                 {
                     influx.reset(new double(influx_coeff));
                 }
             }
 
 
-            aqurecords.at(aquanconRecordIdx).influx_mult_per_record.resize(global_index_per_record_size,influx_mult);
-            aqurecords.at(aquanconRecordIdx).face_per_record.resize(global_index_per_record_size,faceDir);
-            aqurecords.at(aquanconRecordIdx).record_index_per_record.resize(global_index_per_record_size,aquanconRecordIdx);
+            const double influx_mult = aquanconRecord.getItem("INFLUX_MULT").getSIDouble(0);
+
+            aqurecord.influx_mult_per_record.resize(global_index_per_record_size, influx_mult);
+            aqurecord.face_per_record.resize(global_index_per_record_size, faceDir);
+            aqurecord.record_index_per_record.resize(global_index_per_record_size, aquanconRecordIdx);
         }
 
         // Collate_function
@@ -113,10 +119,45 @@ namespace Opm {
     }
 
 
+    bool Aquancon::cellInsideReservoirAndActive(const Opm::EclipseGrid& grid, const int i, const int j, const int k)
+    {
+        if ( i < 0 || j < 0 || k < 0
+            || size_t(i) > grid.getNX() - 1
+            || size_t(j) > grid.getNY() - 1
+            || size_t(k) > grid.getNZ() - 1 )
+        {
+            return false;
+        }
+
+        return grid.cellActive(i, j, k );
+    }
+
+    bool Aquancon::neighborCellInsideReservoirAndActive(const Opm::EclipseGrid& grid,
+           const int i, const int j, const int k, const Opm::FaceDir::DirEnum faceDir)
+    {
+        switch(faceDir) {
+        case FaceDir::XMinus:
+            return cellInsideReservoirAndActive(grid, i - 1, j, k);
+        case FaceDir::XPlus:
+            return cellInsideReservoirAndActive(grid, i + 1, j, k);
+        case FaceDir::YMinus:
+            return cellInsideReservoirAndActive(grid, i, j - 1, k);
+        case FaceDir::YPlus:
+            return cellInsideReservoirAndActive(grid, i, j + 1, k);
+        case FaceDir::ZMinus:
+            return cellInsideReservoirAndActive(grid, i, j, k - 1);
+        case FaceDir::ZPlus:
+            return cellInsideReservoirAndActive(grid, i, j, k + 1);
+        default:
+            throw std::runtime_error("Unknown FaceDir enum " + std::to_string(faceDir));
+        }
+    }
+
     // This function is used to convert from a per record vector to a per aquifer ID vector.
     void Aquancon::collate_function(std::vector<Aquancon::AquanconOutput>& output_vector,
                                     std::vector<Opm::AquanconRecord>& aqurecords,
-                                    std::vector<int> aquiferID_per_record, int m_maxAquID)
+                                    const std::vector<int>& aquiferID_per_record,
+                                    const int m_maxAquID)
     {
         output_vector.resize(m_maxAquID);
         // Find record indices at which the aquifer ids are located in
@@ -257,7 +298,7 @@ namespace Opm {
     }
 
     void Aquancon::convert_record_id_to_aquifer_id(std::vector<int>& record_indices_matching_id,
-                                                   int i, std::vector<int> aquiferID_per_record)
+                                                   const int i, const std::vector<int>& aquiferID_per_record)
     {
         auto it = std::find_if( aquiferID_per_record.begin(), aquiferID_per_record.end(),
                                     [&](int id) {

src/opm/parser/eclipse/share/keywords/000_Eclipse100/A/AQUFETP

@@ -6,6 +6,6 @@
         {"name" : "V0" ,              		"value_type" : "DOUBLE" ,                        "dimension" : "Length*Length*Length"},
         {"name" : "C_T" , 			"value_type" : "DOUBLE" ,                        "dimension" : "1/Pressure"},
         {"name" : "PI" ,              		"value_type" : "DOUBLE" ,                        "dimension" : "ReservoirVolume/Pressure*Time"},
-        {"name" : "TABLE_NUM_WATER_PRESS" ,     "value_type" : "INT"    , "default_value" : 1},
-        {"name" : "SALINITY" ,        		"value_type" : "DOUBLE" , "default_value" : 0 , "dimension" : "Salinity"},
+        {"name" : "TABLE_NUM_WATER_PRESS" ,     "value_type" : "INT"    , "default" : 1},
+        {"name" : "SALINITY" ,        		"value_type" : "DOUBLE" , "default" : 0 , "dimension" : "Salinity"},
         {"name" : "TEMP" ,            		"value_type" : "DOUBLE" , "                      dimension" : "Temperature"}]}

src/opm/parser/eclipse/share/keywords/000_Eclipse100/M/MULT_XYZ

@@ -1,5 +1,5 @@
 {"name" : "MULT_XYZ" ,
- "sections" : ["GRID"],
+ "sections" : ["GRID", "EDIT"],
 
  "deck_names" : [
   "MULTTHT",

tests/parser/AquanconTests.cpp

@@ -49,8 +49,8 @@ inline Deck createAQUANCONDeck_DEFAULT_INFLUX2() {
         "SOLUTION\n"
         "\n"
         "AQUANCON\n"
-        "   1      1  1  1    1   1  1  J-  1.0 /\n"
-        "   1      1  1  1    1   1  1  J-   /\n"
+        "   1      2  2  1    1   1  1  J-  1.0 /\n"
+        "   1      2  2  1    1   1  1  J-   /\n"
         "/ \n";
 
     Parser parser;
@@ -170,3 +170,53 @@ BOOST_AUTO_TEST_CASE(AquanconTest_DEFAULT_INFLUX) {
     EclipseState eclState2( deck2 );
     BOOST_CHECK_THROW(Aquancon( eclState2.getInputGrid(), deck2), std::invalid_argument);
 }
+
+// allowing aquifer exists inside the reservoir
+inline Deck createAQUANCONDeck_ALLOW_INSIDE_AQUAN_OR_NOT() {
+    const char *deckData =
+        "DIMENS\n"
+        "3 3 3 /\n"
+        "\n"
+        "GRID\n"
+        "\n"
+        "ACTNUM\n"
+        " 0 8*1 0 8*1 0 8*1 /\n"
+        "DXV\n"
+        "1 1 1 /\n"
+        "\n"
+        "DYV\n"
+        "1 1 1 /\n"
+        "\n"
+        "DZV\n"
+        "1 1 1 /\n"
+        "\n"
+        "TOPS\n"
+        "9*100 /\n"
+        "\n"
+        "SOLUTION\n"
+        "\n"
+        "AQUFETP\n"
+        "  1 20.0 1000.0 2000. 0.000001 200.0 /\n"
+        "  2 20.0 1000.0 2000. 0.000001 200.0 /\n"
+        "/\n"
+        "AQUANCON\n"
+        "   1      1  1   1   1   1  1  J-  2* YES /\n"
+        "   1      2  2   1   1   1  1  J-  2* YES /\n"
+        "   1      2  2   2   2   1  1  J-  2* YES /\n"
+        "   2      1  1   1   1   3  3  J-  2* NO /\n"
+        "   2      2  2   1   1   3  3  J-  2* NO /\n"
+        "   2      2  2   2   2   3  3  J-  2* NO /\n"
+        "/ \n";
+
+    Parser parser;
+    return parser.parseString(deckData);
+}
+
+BOOST_AUTO_TEST_CASE(AquanconTest_ALLOW_AQUIFER_INSIDE_OR_NOT) {
+    auto deck = createAQUANCONDeck_ALLOW_INSIDE_AQUAN_OR_NOT();
+    const EclipseState eclState( deck );
+    const Aquancon aqucon( eclState.getInputGrid(), deck);
+    const std::vector<Aquancon::AquanconOutput>& aquifer_cons = aqucon.getAquOutput();
+    BOOST_CHECK_EQUAL(aquifer_cons[0].global_index.size(), 2);
+    BOOST_CHECK_EQUAL(aquifer_cons[1].global_index.size(), 1);
+}