diff --git a/opm/input/eclipse/EclipseState/Grid/FieldProps.cpp b/opm/input/eclipse/EclipseState/Grid/FieldProps.cpp index 1651643d4a2..19b02b1cb65 100644 --- a/opm/input/eclipse/EclipseState/Grid/FieldProps.cpp +++ b/opm/input/eclipse/EclipseState/Grid/FieldProps.cpp @@ -673,43 +673,6 @@ bool FieldProps::supported(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_tsecond.data[0])>>; - 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()); - } - this->double_data.erase(mult_iter); - } - multiplier_kw_infos_.clear(); -} - template <> Fieldprops::FieldData& FieldProps::init_get(const std::string& keyword_name, const Fieldprops::keywords::keyword_info& kw_info, bool multiplier_in_edit) { @@ -814,6 +777,63 @@ bool FieldProps::has(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(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& porv) { if (this->has("MULTPV")) { const auto& multpv = this->get("MULTPV"); - std::transform(porv_data.begin(), porv_data.end(), multpv.begin(), porv_data.begin(), std::multiplies()); + 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 index 26b187aba7e..6db998ae4d4 100644 --- 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> double_keywords = {{"DISPERC",keyword_info{}.unit_string("Length")}, {"MINPVV", keyword_info{}.unit_string("ReservoirVolume").global_kw(true)}, - {"MULTPV", keyword_info{}.init(1.0)}, + {"MULTPV", keyword_info{}.init(1.0).mult(true)}, {"NTG", keyword_info{}.init(1.0)}, {"PORO", keyword_info{}.distribute_top(true)}, {"PERMX", keyword_info{}.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> double_keywords = {{"MULTPV", keyword_info{}.init(1.0)}, +static const std::unordered_map> double_keywords = {{"MULTPV", keyword_info{}.init(1.0).mult(true)}, {"PORV", keyword_info{}.unit_string("ReservoirVolume")}, {"MULTX", keyword_info{}.init(1.0).mult(true)}, {"MULTX-", keyword_info{}.init(1.0).mult(true)}, diff --git a/opm/output/eclipse/WriteInit.cpp b/opm/output/eclipse/WriteInit.cpp index b1d0af312ef..f69ed59479c 100644 --- 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 index 4da64c52f84..9c550c69603 100644 --- 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> +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 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>& 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()); + 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("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> +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 expected(10, 1.); + expected.insert(expected.end(), 90, 1.5); + return { deckString, expected }; +} + +BOOST_AUTO_TEST_CASE(MULTPVBOXInit) +{ + checkMULTPV(createMULTPVBOXDECK()); +}