Merge pull request #4021 from blattms/feature/ite-multpv-to-init-v2

Write MULTPV to INIT-file
OPM · May 13, 2024 · 8900b5d · 8900b5d
2 parents faefd4c + c579404
commit 8900b5d
Show file tree

Hide file tree

Showing 4 changed files with 269 additions and 42 deletions.
diff --git a/opm/input/eclipse/EclipseState/Grid/FieldProps.cpp b/opm/input/eclipse/EclipseState/Grid/FieldProps.cpp
@@ -673,43 +673,6 @@ bool FieldProps::supported<int>(const std::string& keyword) {
     return Fieldprops::keywords::isFipxxx(keyword);
 }
 
-void FieldProps::apply_multipliers()
-{
-    static const auto prefix = getMultiplierPrefix();
-
-    for(const auto& [mult_keyword, kw_info]: multiplier_kw_infos_)
-    {
-        const std::string keyword = mult_keyword.substr(prefix.size());
-        auto mult_iter = this->double_data.find(mult_keyword);
-        assert(mult_iter != this->double_data.end());
-        auto iter = this->double_data.find(keyword);
-        if (iter == this->double_data.end()) {
-            iter = this->double_data
-                .emplace(std::piecewise_construct,
-                         std::forward_as_tuple(keyword),
-                         std::forward_as_tuple(kw_info, this->active_size, kw_info.global ? this->global_size : 0))
-                .first;
-        }
-        using Scalar = typename std::remove_cv_t<std::remove_reference_t<decltype(iter->second.data[0])>>;
-        std::transform(iter->second.data.begin(), iter->second.data.end(),
-                       mult_iter->second.data.begin(), iter->second.data.begin(),
-                       std::multiplies<Scalar>());
-
-        // If data is global, then we also need to set the global_data. I think they should be the same at this stage, though!
-        if (kw_info.global)
-        {
-            assert(mult_iter->second.global_data.has_value() && iter->second.global_data.has_value());
-            std::transform(iter->second.global_data->begin(),
-                           iter->second.global_data->end(),
-                           mult_iter->second.global_data->begin(),
-                           iter->second.global_data->begin(),
-                           std::multiplies<Scalar>());
-        }
-        this->double_data.erase(mult_iter);
-    }
-    multiplier_kw_infos_.clear();
-}
-
 template <>
 Fieldprops::FieldData<double>& FieldProps::init_get(const std::string& keyword_name, const Fieldprops::keywords::keyword_info<double>& kw_info,
                                                     bool multiplier_in_edit) {
@@ -814,6 +777,63 @@ bool FieldProps::has<int>(const std::string& keyword) const {
 }
 
 
+void FieldProps::apply_multipliers()
+{
+    // We need to manually search for PORV in the map here instead of using
+    // the get method. The latter one will compute PORV from the cell volume
+    // using MULTPV, NTG, and PORO. Our intend here in the EDIT section is
+    // is to multiply exsisting MULTPV with new new ones and consistently
+    // change PORV as well. If PORV has been created before this is as easy
+    // as multiplying it and MULTPV with the additional MULTPV. In the other
+    // case we do not create PORV at all but just change MULTPV as PORV will
+    // be correctly computed from it.
+    // Hence we need to know whether PORV is already there
+    const bool hasPorvBefore =
+        (this->double_data.find(ParserKeywords::PORV::keywordName) ==
+         this->double_data.end());
+    static const auto prefix = getMultiplierPrefix();
+
+    for(const auto& [mult_keyword, kw_info]: multiplier_kw_infos_)
+    {
+        const std::string keyword = mult_keyword.substr(prefix.size());
+        auto mult_iter = this->double_data.find(mult_keyword);
+        assert(mult_iter != this->double_data.end());
+        auto iter = this->double_data.find(keyword);
+        if (iter == this->double_data.end()) {
+            iter = this->double_data
+                .emplace(std::piecewise_construct,
+                         std::forward_as_tuple(keyword),
+                         std::forward_as_tuple(kw_info, this->active_size, kw_info.global ? this->global_size : 0))
+                .first;
+        }
+
+        std::transform(iter->second.data.begin(), iter->second.data.end(),
+                       mult_iter->second.data.begin(), iter->second.data.begin(),
+                       std::multiplies<>());
+
+        // If data is global, then we also need to set the global_data. I think they should be the same at this stage, though!
+        if (kw_info.global)
+        {
+            assert(mult_iter->second.global_data.has_value() && iter->second.global_data.has_value());
+            std::transform(iter->second.global_data->begin(),
+                           iter->second.global_data->end(),
+                           mult_iter->second.global_data->begin(),
+                           iter->second.global_data->begin(),
+                           std::multiplies<>());
+        }
+        // If this is MULTPV we also need to apply the additional multiplier to PORV if that was initialized already.
+        // Note that the check for PORV is essential as otherwise the value constructed durig init_get will already apply
+        // current MULTPV.
+        if (keyword == ParserKeywords::MULTPV::keywordName && !hasPorvBefore) {
+            auto& porv = this->init_get<double>(ParserKeywords::PORV::keywordName);
+            auto& porv_data = porv.data;
+            std::transform(porv_data.begin(), porv_data.end(), mult_iter->second.data.begin(), porv_data.begin(), std::multiplies<>());
+        }
+        this->double_data.erase(mult_iter);
+    }
+    multiplier_kw_infos_.clear();
+}
+
 /*
   The ACTNUM and PORV keywords are special cased with quite extensive
   postprocessing, and should be access through the special ::porv() and
@@ -1245,7 +1265,7 @@ void FieldProps::init_porv(Fieldprops::FieldData<double>& porv) {
 
     if (this->has<double>("MULTPV")) {
         const auto& multpv = this->get<double>("MULTPV");
-        std::transform(porv_data.begin(), porv_data.end(), multpv.begin(), porv_data.begin(), std::multiplies<double>());
+        std::transform(porv_data.begin(), porv_data.end(), multpv.begin(), porv_data.begin(), std::multiplies<>());
     }
 
     for (const auto& mregp: this->multregp) {

diff --git a/opm/input/eclipse/EclipseState/Grid/FieldProps.hpp b/opm/input/eclipse/EclipseState/Grid/FieldProps.hpp
@@ -133,7 +133,7 @@ namespace ALIAS {
 namespace GRID {
 static const std::unordered_map<std::string, keyword_info<double>> double_keywords = {{"DISPERC",keyword_info<double>{}.unit_string("Length")},
                                                                                       {"MINPVV",  keyword_info<double>{}.unit_string("ReservoirVolume").global_kw(true)},
-                                                                                      {"MULTPV",  keyword_info<double>{}.init(1.0)},
+                                                                                      {"MULTPV",  keyword_info<double>{}.init(1.0).mult(true)},
                                                                                       {"NTG",     keyword_info<double>{}.init(1.0)},
                                                                                       {"PORO",    keyword_info<double>{}.distribute_top(true)},
                                                                                       {"PERMX",   keyword_info<double>{}.unit_string("Permeability").distribute_top(true)},
@@ -181,7 +181,7 @@ namespace EDIT {
   due to the special treatment of the TRAN fields.
 */
 
-static const std::unordered_map<std::string, keyword_info<double>> double_keywords = {{"MULTPV",  keyword_info<double>{}.init(1.0)},
+static const std::unordered_map<std::string, keyword_info<double>> double_keywords = {{"MULTPV",  keyword_info<double>{}.init(1.0).mult(true)},
                                                                                       {"PORV",    keyword_info<double>{}.unit_string("ReservoirVolume")},
                                                                                       {"MULTX",   keyword_info<double>{}.init(1.0).mult(true)},
                                                                                       {"MULTX-",  keyword_info<double>{}.init(1.0).mult(true)},

diff --git a/opm/output/eclipse/WriteInit.cpp b/opm/output/eclipse/WriteInit.cpp
@@ -667,10 +667,15 @@ void Opm::InitIO::write(const ::Opm::EclipseState&              es,
     {
         const auto writeAll = es.cfg().io().writeAllTransMultipliers();
 
-        const auto tranMult = es.getTransMult()
+        auto multipliers = es.getTransMult()
             .convertToSimProps(grid.getNumActive(), writeAll);
+        if (es.fieldProps().has_double("MULTPV")) {
+            multipliers.insert("MULTPV", UnitSystem::measure::identity,
+                        es.fieldProps().get_double("MULTPV"),
+                        data::TargetType::INIT);
+        }
 
-        writeSimulatorProperties(grid, tranMult, initFile);
+        writeSimulatorProperties(grid, multipliers, initFile);
     }
 
     writeTableData(es, units, initFile);

diff --git a/tests/test_EclipseIO.cpp b/tests/test_EclipseIO.cpp
@@ -739,3 +739,205 @@ BOOST_AUTO_TEST_CASE(MULTXYZInit)
     testMultxyz({ '3', '3' , '2'}, true);
     testMultxyz({ '3', '3' , '3'});
 }
+
+
+std::pair<std::string,std::vector<float>>
+createMULTPVDECK(bool edit)
+{
+    auto deckString = std::string { R"(RUNSPEC
+
+TITLE
+   1D OIL WATER
+
+
+
+DIMENS
+   100 1 1 /
+
+EQLDIMS
+/
+TABDIMS
+  2 1 100 /
+
+OIL
+WATER
+
+ENDSCALE
+/
+
+METRIC
+
+START
+   1 'JAN' 2024 /
+
+WELLDIMS
+   3 3 2 2 /
+
+UNIFIN
+UNIFOUT
+
+GRID
+
+INIT
+
+DX 
+  100*1 /
+DY
+        100*10 /
+DZ
+  100*1 /
+
+TOPS
+  100*2000 /
+
+PORO
+  100*0.3 /
+
+MULTPV
+  100*5 /  -- overwritten by the next MULTPV keyword
+
+MULTPV
+  100*10 /  -- partly overwritten by the next BOX statements
+
+BOX
+  69 75 1 1 1 1 /
+MULTPV
+  2*1 5*1.5 /
+ENDBOX
+
+BOX
+  76 85 1 1 1 1 /
+MULTPV
+  5*1.5 5*1 /
+ENDBOX
+
+EDIT
+)"};
+    std::vector<float> multpv(68, 10.);
+    multpv.insert(multpv.end(), 2, 1);
+    multpv.insert(multpv.end(), 10, 1.5);
+    multpv.insert(multpv.end(), 5, 1.);
+    multpv.insert(multpv.end(), 15, 10.);
+
+    if (edit) {
+        deckString += std::string { R"(MULTPV
+  75*10 25*10 / -- overwritten by the next
+MULTPV
+  75*0.8 25*1 /
+)"};
+        std::transform(multpv.begin(), multpv.begin() + 75, multpv.begin(),
+                       [](float f) { return f*.8; });
+    }
+    deckString += std::string { R"(PROPS
+
+SOLUTION
+SCHEDULE
+)"};
+    return { deckString, multpv };
+}
+
+void checkMULTPV(const std::pair<std::string, std::vector<float>>& deckAndValues)
+{
+    const auto [deckString, expectedMultPV] = deckAndValues;
+    const auto deck = Parser().parseString(deckString);
+    auto es = EclipseState( deck );
+    const auto& eclGrid = es.getInputGrid();
+    const Schedule schedule(deck, es, std::make_shared<Python>());
+    const SummaryConfig summary_config( deck, schedule, es.fieldProps(), es.aquifer());
+    const SummaryState st(TimeService::now());
+    es.getIOConfig().setBaseName( "MULTPVFOO" );
+    EclipseIO eclWriter( es, eclGrid , schedule, summary_config);
+    eclWriter.writeInitial( );
+
+    EclIO::EclFile initFile { "MULTPVFOO.INIT" };
+    BOOST_CHECK_MESSAGE( initFile.hasKey("MULTPV"),
+                         R"(INIT file must have MULTPV array)" );
+    const auto& multPVValues = initFile.get<float>("MULTPV");
+    auto expect = expectedMultPV.begin();
+    BOOST_CHECK(multPVValues.size() == expectedMultPV.size());
+
+    for (auto multVal = multPVValues.begin(); multVal != multPVValues.end();
+         ++multVal, ++expect) {
+        BOOST_CHECK_CLOSE(*multVal, *expect, 1e-8);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(MULTPVInit)
+{
+    checkMULTPV(createMULTPVDECK(false));
+    checkMULTPV(createMULTPVDECK(true));
+
+}
+
+std::pair<std::string,std::vector<float>>
+createMULTPVBOXDECK()
+{
+    auto deckString = std::string { R"(RUNSPEC
+
+TITLE
+   1D OIL WATER
+
+
+
+DIMENS
+   100 1 1 /
+
+EQLDIMS
+/
+TABDIMS
+  2 1 100 /
+
+OIL
+WATER
+
+ENDSCALE
+/
+
+METRIC
+
+START
+   1 'JAN' 2024 /
+
+WELLDIMS
+   3 3 2 2 /
+
+UNIFIN
+UNIFOUT
+
+GRID
+
+INIT
+
+DX 
+  100*1 /
+DY
+        100*10 /
+DZ
+  100*1 /
+
+TOPS
+  100*2000 /
+
+PORO
+  100*0.3 /
+
+BOX
+  11 100 1 1 1 1 /
+MULTPV
+  90*1.5 /
+ENDBOX
+
+EDIT
+PROPS
+SOLUTION
+SCHEDULE
+)"};
+    std::vector<float> expected(10, 1.);
+    expected.insert(expected.end(), 90, 1.5);
+    return { deckString, expected };
+}
+
+BOOST_AUTO_TEST_CASE(MULTPVBOXInit)
+{
+    checkMULTPV(createMULTPVBOXDECK());
+}