diff --git a/opm/output/eclipse/RestartValue.hpp b/opm/output/eclipse/RestartValue.hpp index 08602ffdc..e766920ea 100644 --- a/opm/output/eclipse/RestartValue.hpp +++ b/opm/output/eclipse/RestartValue.hpp @@ -71,6 +71,9 @@ namespace Opm { void addExtra(const std::string& key, UnitSystem::measure dimension, std::vector data); void addExtra(const std::string& key, std::vector data); const std::vector& getExtra(const std::string& key) const; + + void convertFromSI(const UnitSystem& units); + void convertToSI(const UnitSystem& units); }; } diff --git a/src/opm/output/eclipse/AggregateGroupData.cpp b/src/opm/output/eclipse/AggregateGroupData.cpp index 9d19308e3..c97764707 100644 --- a/src/opm/output/eclipse/AggregateGroupData.cpp +++ b/src/opm/output/eclipse/AggregateGroupData.cpp @@ -16,9 +16,11 @@ You should have received a copy of the GNU General Public License along with OPM. If not, see . */ -#include + #include +#include + #include #include #include @@ -146,14 +148,14 @@ namespace { //const std::vector< const Opm::Group* > groups = sched.getGroups(simStep); const auto& groups = sched.getGroups(simStep); // make group index for current report step - std::map groupIndexMap; + std::map indexGroupMap; for (const auto* group : groups) { int ind = (group->name() == "FIELD") ? ngmaxz(inteHead)-1 : group->seqIndex()-1; const std::pair groupPair = std::make_pair(static_cast(ind), group); - groupIndexMap.insert(groupPair); + indexGroupMap.insert(groupPair); } - return groupIndexMap; + return indexGroupMap; } std::map currentGroupMapNameIndex(const Opm::Schedule& sched, const size_t simStep, const std::vector& inteHead) @@ -242,6 +244,7 @@ namespace { void staticContrib(const Opm::Schedule& sched, const Opm::Group& group, const int nwgmax, + const int ngmaxz, const std::size_t simStep, IGrpArray& igrp, const std::vector& inteHead) @@ -249,34 +252,42 @@ namespace { // find the number of wells or child groups belonging to a group and store in // location nwgmax +1 in the igrp array - const auto& childGroups = sched.getChildGroups(group.name(), simStep); - const auto& childWells = sched.getChildWells(group.name(), simStep); + const auto childGroups = sched.getChildGroups(group.name(), simStep); + const auto childWells = sched.getChildWells(group.name(), simStep); //std::cout << "IGrpArray - staticContrib - before currentGroupMapIndexGroup(sched, simStep, inteHead)" << std::endl; //const auto& groupMapIndexGroup = currentGroupMapIndexGroup(sched, simStep, inteHead); - std::cout << "IGrpArray - staticContrib - before currentGroupMapNameIndex(sched, simStep, inteHead)" << std::endl; - const auto& groupMapNameIndex = currentGroupMapNameIndex(sched, simStep, inteHead); + //std::cout << "IGrpArray - staticContrib - before currentGroupMapNameIndex(sched, simStep, inteHead)" << std::endl; + const auto groupMapNameIndex = currentGroupMapNameIndex(sched, simStep, inteHead); if ((childGroups.size() != 0) && (childWells.size()!=0)) throw std::invalid_argument("group has both wells and child groups" + group.name()); int igrpCount = 0; if (childWells.size() != 0) { //group has child wells - for (const auto* well : childWells ) { + //std::cout << "IGrpArray - staticContrib: childwells for group.name(): " << group.name() << "childWells - size: " << childWells.size() << std::endl; + for (const auto well : childWells ) { + //std::cout << "Child well name: " << well->name() << " Well seqIndex(): " << well->seqIndex() << std::endl; igrp[igrpCount] = well->seqIndex()+1; igrpCount+=1; + //std::cout << "childWells: igrpCount after increment: " << igrpCount << std::endl; } } else if (childGroups.size() != 0) { //group has child groups //The field group always has seqIndex = 0 because it is always defined first //Hence the all groups except the Field group uses the seqIndex assigned - for (const auto& group : childGroups ) { - igrp[igrpCount] = group->seqIndex()+1; + //std::cout << "IGrpArray - staticContrib: childGroups for group.name(): " << group.name() << "childGroups - size: " << childGroups.size() << std::endl; + for (const auto grp : childGroups ) { + //std::cout << "Child Group name: " << grp->name() << " Group seqIndex(): " << grp->seqIndex()-1 << std::endl; + igrp[igrpCount] = grp->seqIndex(); igrpCount+=1; + //std::cout << "childGroups: igrpCount after increment: " << igrpCount << std::endl; } } - //assign the number of child wells or child groups to location - // nwgmax + //assign the number of child wells or child groups to + // location nwgmax + //std::cout << "IGrpArray - staticContrib: childGroups.size() " << childGroups.size() << "childWells.size(): " << childWells.size() << std::endl; + //std::cout << "IGrpArray - staticContrib: nwgmax: " << nwgmax << std::endl; igrp[nwgmax] = (childGroups.size() == 0) ? childWells.size() : childGroups.size(); @@ -290,11 +301,34 @@ namespace { const auto grpLevel = currentGroupLevel(sched, group, simStep); igrp[nwgmax+27] = grpLevel; + // set values for group probably connected to GCONPROD settings + // + if (group.name() != "FIELD") + { + igrp[nwgmax+ 5] = -1; + igrp[nwgmax+12] = -1; + igrp[nwgmax+17] = -1; + igrp[nwgmax+22] = -1; + + //assign values to group number (according to group sequence) + igrp[nwgmax+88] = group.seqIndex(); + igrp[nwgmax+89] = group.seqIndex(); + igrp[nwgmax+95] = group.seqIndex(); + igrp[nwgmax+96] = group.seqIndex(); + } + else + { + //assign values to group number (according to group sequence) + igrp[nwgmax+88] = ngmaxz; + igrp[nwgmax+89] = ngmaxz; + igrp[nwgmax+95] = ngmaxz; + igrp[nwgmax+96] = ngmaxz; + } //find parent group and store index of parent group in //location nwgmax + 28 const auto& groupTree = sched.getGroupTree( simStep ); const std::string& parent = groupTree.parent(group.name()); - std::cout << "IGrpArray - staticContrib - before groupMapNameIndex.find(parent)" << std::endl; + //std::cout << "IGrpArray - staticContrib - before groupMapNameIndex.find(parent)" << std::endl; if (group.name() == "FIELD") igrp[nwgmax+28] = 0; else { @@ -323,7 +357,52 @@ namespace { WV::WindowSize{ entriesPerGroup(inteHead) } }; } - } + + template + void staticContrib(SGrpArray& sGrp) + { + const auto dflt = -1.0e+20f; + const auto dflt_2 = -2.0e+20f; + const auto infty = 1.0e+20f; + const auto zero = 0.0f; + const auto one = 1.0f; + + const auto init = std::vector { // 132 Items (0..131) + // 0 1 2 3 4 + infty, infty, dflt , infty, zero , // 0.. 4 ( 0) + zero , infty, infty, infty , infty, // 5.. 9 ( 1) + infty, infty, infty, infty , dflt , // 10.. 14 ( 2) + infty, infty, infty, infty , dflt , // 15.. 19 ( 3) + infty, infty, infty, infty , dflt , // 20.. 24 ( 4) + zero , zero , zero , dflt_2, zero , // 24.. 29 ( 5) + zero , zero , zero , zero , zero , // 30.. 34 ( 6) + infty ,zero , zero , zero , infty, // 35.. 39 ( 7) + zero , zero , zero , zero , zero , // 40.. 44 ( 8) + zero , zero , zero , zero , zero , // 45.. 49 ( 9) + zero , infty, infty, infty , infty, // 50.. 54 (10) + infty, infty, infty, infty , infty, // 55.. 59 (11) + infty, infty, infty, infty , infty, // 60.. 64 (12) + infty, infty, infty, infty , zero , // 65.. 69 (13) + zero , zero , zero , zero , zero , // 70.. 74 (14) + zero , zero , zero , zero , infty, // 75.. 79 (15) + infty, zero , infty, zero , zero , // 80.. 84 (16) + zero , zero , zero , zero , zero , // 85.. 89 (17) + zero , zero , one , zero , zero , // 90.. 94 (18) + zero , zero , zero , zero , zero , // 95.. 99 (19) + zero , zero , zero , zero , zero , // 100..104 (20) + zero , zero , zero , zero , zero , // 105..109 (21) + zero , zero // 110..111 (22) + }; + + const auto sz = static_cast< + decltype(init.size())>(sGrp.size()); + + auto b = std::begin(init); + auto e = b + std::min(init.size(), sz); + + std::copy(b, e, std::begin(sGrp)); + } + } // SGrp #if 0 template @@ -506,52 +585,49 @@ captureDeclaredGroupData(const Schedule& sched, const std::size_t simStep, const std::vector& inteHead) { - std::cout << "captureDeclaredGroupData before creation groups" << std::endl; - //const auto& groups = sched.getGroups(simStep); std::map indexGroupMap = currentGroupMapIndexGroup(sched, simStep, inteHead); - std::cout << "captureDeclaredGroupData after creation groups" << std::endl; + std::map nameIndexMap = currentGroupMapNameIndex(sched, simStep, inteHead); std::vector curGroups; - curGroups.reserve(ngmaxz(inteHead)+1); + curGroups.resize(ngmaxz(inteHead)); + // + //initialize curgroups + for (auto* cg : curGroups) cg = nullptr; std::map ::iterator it = indexGroupMap.begin(); - while (it != indexGroupMap.end()) { - std::cout << "captureDeclaredGroupData before assign curGroups, it-first: " << it->first << " it-> second->name(): " << it->second->name() << std::endl; - curGroups[static_cast(it->first)] = it->second; - it++; - } + while (it != indexGroupMap.end()) + { + curGroups[static_cast(it->first)] = it->second; + it++; + } - // - std::cout << "captureDeclaredGroupData curGroups.size(): " << curGroups.size() << std::endl; - - std::cout << "curGroups -index: " << 0 << " names: " << curGroups[0]->name() << std::endl; - - //std::cout << "curGroups -index: " << i << " names: " << curGroups[i]->name() << std::endl; for (const auto* grp : curGroups) { - std::cout << "curGroups grp " << grp->name() << std::endl; + if (grp != nullptr) + { + const std::string gname = grp->name(); + const auto itr = nameIndexMap.find(gname); + const int ind = itr->second; + } } { - std::cout << "captureDeclaredGroupData before loop groupLoop" << std::endl; - groupLoop(curGroups, [sched, simStep, inteHead, this] + groupLoop(curGroups, [sched, simStep, inteHead, this] (const Group& group, const std::size_t groupID) -> void - { - auto ig = this->iGroup_[groupID]; - std::cout << "captureDeclaredGroupData: groupLoop, groupID: " << groupID << std::endl; - IGrp::staticContrib(sched, group, this->nGMaxz_, simStep, ig, inteHead); - }); + { + auto ig = this->iGroup_[groupID]; + IGrp::staticContrib(sched, group, this->nWGMax_, this->nGMaxz_, simStep, ig, inteHead); + }); } -#if 0 // Define Static Contributions to SWell Array. - wellLoop(wells, - [this](const Well& /* well */, const std::size_t wellID) -> void + groupLoop(curGroups, + [this](const Group& group, const std::size_t groupID) -> void { - auto sw = this->sWell_[wellID]; - - SWell::staticContrib(sw); + auto sw = this->sGroup_[groupID]; + SGrp::staticContrib(sw); }); +#if 0 // Define Static Contributions to ZWell Array. wellLoop(wells, [this](const Well& well, const std::size_t wellID) -> void diff --git a/src/opm/output/eclipse/RestartIO.cpp b/src/opm/output/eclipse/RestartIO.cpp index 36fb8d33d..ee3f41299 100644 --- a/src/opm/output/eclipse/RestartIO.cpp +++ b/src/opm/output/eclipse/RestartIO.cpp @@ -250,8 +250,8 @@ data::Wells restore_wells( const ::Opm::RestartIO::ecl_kw_type * opm_xwel, } const auto well_size = [&]( size_t acc, const Well* w ) { - return acc - + 2 + phases.size() + return acc + + 2 + phases.size() + ( w->getConnections( sim_step ).size() * (phases.size() + data::Connection::restart_size) ); }; @@ -308,7 +308,7 @@ data::Wells restore_wells( const ::Opm::RestartIO::ecl_kw_type * opm_xwel, return wells; } -//} +} //* should take grid as argument because it may be modified from the simulator */ @@ -359,7 +359,6 @@ RestartValue load( const std::string& filename, rst_value.addExtra(key, extra.dim, {data_ptr, end_ptr}); } else if (required) throw std::runtime_error("No such key in file: " + key); - } // Convert solution fields and extra data from user units to SI @@ -374,8 +373,8 @@ RestartValue load( const std::string& filename, return rst_value; } -//namespace { +namespace { std::vector serialize_ICON( int sim_step, int ncwmax, @@ -455,10 +454,10 @@ std::vector< int > serialize_OPM_IWEL( const data::Wells& wells, } std::vector< double > serialize_OPM_XWEL( const data::Wells& wells, - int sim_step, - const std::vector< const Well* >& sched_wells, - const Phases& phase_spec, - const EclipseGrid& grid ) { + int sim_step, + const std::vector< const Well* >& sched_wells, + const Phases& phase_spec, + const EclipseGrid& grid ) { using rt = data::Rates::opt; @@ -539,7 +538,7 @@ void write_kw(::Opm::RestartIO::ecl_rst_file_type * rst_file , Opm::RestartIO::E } /* void writeHeader(::Opm::RestartIO::ecl_rst_file_type * rst_file, - int sim_step, + int sim_step, int report_step, time_t posix_time, double sim_days, @@ -621,14 +620,14 @@ void writeGroup(::Opm::RestartIO::ecl_rst_file_type * rst_file, // find inteHead //const auto ih = Helpers::createInteHead(es, grid, schedule, simTime, report_step); // write IGRP to restart file - std::cout << "writeGroup before initializing groupData" << std::endl; + //std::cout << "writeGroup before initializing groupData" << std::endl; const size_t simStep = static_cast (sim_step); auto groupData = Helpers::AggregateGroupData(ih); - std::cout << "writeGroup before captureDeclaredGroupData" << std::endl; + //std::cout << "writeGroup before captureDeclaredGroupData" << std::endl; groupData.captureDeclaredGroupData(schedule, simStep, ih); - std::cout << "writeGroup before write_kw IGRP" << std::endl; - write_kw(rst_file, EclKW("IGRP", groupData.getIGroup())); - //write_kw(rst_file, EclKW("SGRP", groupData.getSGroup())); + //std::cout << "writeGroup before write_kw IGRP" << std::endl; + write_kw(rst_file, EclKW ("IGRP", groupData.getIGroup())); + write_kw(rst_file, EclKW("SGRP", groupData.getSGroup())); } @@ -714,13 +713,12 @@ void writeExtraData(::Opm::RestartIO::ecl_rst_file_type* rst_file, const Restart void writeWell(::Opm::RestartIO::ecl_rst_file_type* rst_file, int sim_step, const EclipseState& es , const EclipseGrid& grid, const Schedule& schedule, const data::Wells& wells) { - const auto sched_wells = schedule.getWells(sim_step); const auto& phases = es.runspec().phases(); const size_t ncwmax = schedule.getMaxNumConnectionsForWells(sim_step); const auto iwel_data = serialize_IWEL(sim_step, sched_wells, grid); - const auto icon_data = serialize_ICON(sim_step, ncwmax, sched_wells, grid); + const auto icon_data = serialize_ICON(sim_step , ncwmax, sched_wells, grid); const auto zwel_data = serialize_ZWEL( sched_wells ); ::Opm::RestartIO::write_kw( rst_file, ::Opm::RestartIO::EclKW< int >( IWEL_KW, iwel_data) ); @@ -759,7 +757,7 @@ void checkSaveArguments(const EclipseState& es, if (thpres.size() != num_regions * num_regions) throw std::runtime_error("THPRES vector has invalid size - should have num_region * num_regions."); - /* const RestartValue& restart_value, + const RestartValue& restart_value, const EclipseGrid& grid) { for (const auto& elm: restart_value.solution) @@ -767,7 +765,7 @@ void checkSaveArguments(const EclipseState& es, throw std::runtime_error("Wrong size on solution vector: " + elm.first); - if (es.getSimulationConfig().getThresholdPressure().size() > 0) { + /*if (es.getSimulationConfig().getThresholdPressure().size() > 0) { // If the the THPRES option is active the restart_value should have a // THPRES field. This is not enforced here because not all the opm // simulators have been updated to include the THPRES values. @@ -779,12 +777,15 @@ void checkSaveArguments(const EclipseState& es, size_t num_regions = es.getTableManager().getEqldims().getNumEquilRegions(); const auto& thpres = restart_value.getExtra("THPRES"); if (thpres.size() != num_regions * num_regions) - throw std::runtime_error("THPRES vector has invalid size - should have num_region * num_regions.");*/ + throw std::runtime_error("THPRES vector has invalid size - should have num_region * num_regions."); + + size_t num_regions = es.getTableManager().getEqldims().getNumEquilRegions(); + const auto& thpres = restart_value.getExtra("THPRES"); + if (thpres.size() != num_regions * num_regions) + throw std::runtime_error("THPRES vector has invalid size - should have num_region * num_regions.");*/ } } - - -//} // Anonymous namespace +} // Anonymous namespace void save(const std::string& filename, diff --git a/src/opm/output/eclipse/RestartValue.cpp b/src/opm/output/eclipse/RestartValue.cpp index 70155d40b..bcecaec10 100644 --- a/src/opm/output/eclipse/RestartValue.cpp +++ b/src/opm/output/eclipse/RestartValue.cpp @@ -81,5 +81,25 @@ namespace Opm { this->addExtra(key, UnitSystem::measure::identity, std::move(data)); } + void RestartValue::convertFromSI(const UnitSystem& units) { + this->solution.convertFromSI(units); + for (auto & extra_value : this->extra) { + const auto& restart_key = extra_value.first; + auto & data = extra_value.second; + + units.from_si(restart_key.dim, data); + } + } + + void RestartValue::convertToSI(const UnitSystem& units) { + this->solution.convertToSI(units); + for (auto & extra_value : this->extra) { + const auto& restart_key = extra_value.first; + auto & data = extra_value.second; + + units.to_si(restart_key.dim, data); + } + } + } diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp index 9bb45f56e..c02187bd7 100644 --- a/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp +++ b/src/opm/parser/eclipse/EclipseState/Schedule/Schedule.cpp @@ -1714,7 +1714,7 @@ namespace Opm { } return groups; } - + void Schedule::addWellToGroup( Group& newGroup, Well& well , size_t timeStep) { const std::string currentGroupName = well.getGroupName(timeStep); if (currentGroupName != "") { diff --git a/tests/parser/ParseContextTests.cpp b/tests/parser/ParseContextTests.cpp index 5f6da3459..8256de6ab 100644 --- a/tests/parser/ParseContextTests.cpp +++ b/tests/parser/ParseContextTests.cpp @@ -712,4 +712,4 @@ BOOST_AUTO_TEST_CASE( test_invalid_wtemplate_config ) { BOOST_CHECK_THROW( Schedule( deckUnSupported , grid , eclipseProperties, Phases(true, true, true) , parseContext), std::invalid_argument ); } -} \ No newline at end of file +} diff --git a/tests/test_serialize_ICON.cpp b/tests/test_serialize_ICON.cpp index 5e0e0d2f8..1c8253380 100644 --- a/tests/test_serialize_ICON.cpp +++ b/tests/test_serialize_ICON.cpp @@ -19,6 +19,10 @@ #include +#include // @@ +#include // @@ +#include // @@ + #define BOOST_TEST_MODULE serialize_ICON_TEST #include @@ -84,6 +88,11 @@ BOOST_AUTO_TEST_CASE( serialize_icon_test ) } w_offset += (ICONZ * ncwmax); } - } -}; + std::copy(icondata.begin(), + icondata.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << std::endl; + BOOST_CHECK_EQUAL(1, 1);// @@@ + } +}