From 8efae160871c5b5b5367ab8c8bc40ae1ed16bb88 Mon Sep 17 00:00:00 2001 From: goncalvesmachadoc Date: Fri, 1 Nov 2019 12:56:20 +0100 Subject: [PATCH] =?UTF-8?q?=C3=A7onverted=20tabs=20to=20spaces?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ebos/ecloutputblackoilmodule.hh | 904 ++++++++++++----------- ebos/eclwriter.hh | 1207 +++++++++++++++---------------- 2 files changed, 1047 insertions(+), 1064 deletions(-) diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 16604d50c..5bbe04bbe 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -114,64 +114,64 @@ class EclOutputBlackOilModule enum WPId { WellLocationi = 0, //WLi - WellLocationj = 1, //WLj + WellLocationj = 1, //WLj OilRate = 2, //OR WaterRate = 3, //WR GasRate = 4, //GR - FluidResVol = 5, //FRV - WaterCut = 6, //WC - GasOilRatio = 7, //GOR - WatGasRatio = 8, //WGR - BHP = 9, //BHP - THP = 10, //THP - SteadyStatePI = 11, //SteadyStatePI - WellName = 0, //WName - CTRLMode = 1, //CTRL + FluidResVol = 5, //FRV + WaterCut = 6, //WC + GasOilRatio = 7, //GOR + WatGasRatio = 8, //WGR + BHP = 9, //BHP + THP = 10, //THP + SteadyStatePI = 11, //SteadyStatePI + WellName = 0, //WName + CTRLMode = 1, //CTRL }; static const int numWPValues = 12; - static const int numWPNames = 2; - }; - struct WellInjDataType + static const int numWPNames = 2; + }; + struct WellInjDataType { enum WIId { WellLocationi = 0, //WLi - WellLocationj = 1, //WLj + WellLocationj = 1, //WLj OilRate = 2, //OR WaterRate = 3, //WR GasRate = 4, //GR - FluidResVol = 5, //FRV - BHP = 6, //BHP - THP = 7, //THP - SteadyStateII = 8, //SteadyStateII - WellName = 0, //WName - CTRLModeOil = 1, //CTRLo - CTRLModeWat = 2, //CTRLw - CTRLModeGas = 3, //CTRLg + FluidResVol = 5, //FRV + BHP = 6, //BHP + THP = 7, //THP + SteadyStateII = 8, //SteadyStateII + WellName = 0, //WName + CTRLModeOil = 1, //CTRLo + CTRLModeWat = 2, //CTRLw + CTRLModeGas = 3, //CTRLg }; static const int numWIValues = 9; - static const int numWINames = 4; + static const int numWINames = 4; }; - struct WellCumDataType + struct WellCumDataType { enum WCId { WellLocationi = 0, //WLi - WellLocationj = 1, //WLj + WellLocationj = 1, //WLj OilProd = 2, //OP WaterProd = 3, //WP GasProd = 4, //GP - FluidResVolProd = 5, //FRVP - OilInj = 6, //OI + FluidResVolProd = 5, //FRVP + OilInj = 6, //OI WaterInj = 7, //WI GasInj = 8, //GI - FluidResVolInj = 9, //FRVI - WellName = 0, //WName - WellType = 1, //WType - WellCTRL = 2, //WCTRL + FluidResVolInj = 9, //FRVI + WellName = 0, //WName + WellType = 1, //WType + WellCTRL = 2, //WCTRL }; static const int numWCValues = 10; - static const int numWCNames = 3; + static const int numWCNames = 3; }; public: @@ -1107,385 +1107,381 @@ public: // write production report to output void outputProdLog(size_t reportStepNum, const bool substep, bool forceDisableProdOutput) { - if (!substep) { - - ScalarBuffer tmp_values(WellProdDataType::numWPValues, 0.0); - StringBuffer tmp_names(WellProdDataType::numWPNames, ""); - outputProductionReport_(tmp_values, tmp_names, forceDisableProdOutput); - - const auto& st = simulator_.vanguard().summaryState(); - const auto& schedule = simulator_.vanguard().schedule(); + if (!substep) { + + ScalarBuffer tmp_values(WellProdDataType::numWPValues, 0.0); + StringBuffer tmp_names(WellProdDataType::numWPNames, ""); + outputProductionReport_(tmp_values, tmp_names, forceDisableProdOutput); + + const auto& st = simulator_.vanguard().summaryState(); + const auto& schedule = simulator_.vanguard().schedule(); - for (const auto& gname: schedule.groupNames()) { + for (const auto& gname: schedule.groupNames()) { - auto gName = static_cast(gname); - auto get = [&st, &gName](const std::string& vector) - { - const auto key = vector + ':' + gName; + auto gName = static_cast(gname); + auto get = [&st, &gName](const std::string& vector) + { + const auto key = vector + ':' + gName; - return st.has(key) ? st.get(key) : 0.0; - }; + return st.has(key) ? st.get(key) : 0.0; + }; - tmp_names[0] = gname; + tmp_names[0] = gname; - tmp_values[2] = get("GOPR"); //WellProdDataType::OilRate - tmp_values[3] = get("GWPR"); //WellProdDataType::WaterRate - tmp_values[4] = get("GGPR"); //WellProdDataType::GasRate - tmp_values[5] = get("GVPR"); //WellProdDataType::FluidResVol - tmp_values[6] = get("GWCT"); //WellProdDataType::WaterCut - tmp_values[7] = get("GGOR"); //WellProdDataType::GasOilRatio - tmp_values[8] = get("GWPR")/get("GGPR"); //WellProdDataType::WaterGasRatio + tmp_values[2] = get("GOPR"); //WellProdDataType::OilRate + tmp_values[3] = get("GWPR"); //WellProdDataType::WaterRate + tmp_values[4] = get("GGPR"); //WellProdDataType::GasRate + tmp_values[5] = get("GVPR"); //WellProdDataType::FluidResVol + tmp_values[6] = get("GWCT"); //WellProdDataType::WaterCut + tmp_values[7] = get("GGOR"); //WellProdDataType::GasOilRatio + tmp_values[8] = get("GWPR")/get("GGPR"); //WellProdDataType::WaterGasRatio - outputProductionReport_(tmp_values, tmp_names, forceDisableProdOutput); - } + outputProductionReport_(tmp_values, tmp_names, forceDisableProdOutput); + } - for (const auto& wname: schedule.wellNames(reportStepNum)) { - - // don't bother with wells not on this process - const auto& defunctWellNames = simulator_.vanguard().defunctWellNames(); - if (defunctWellNames.find(wname) != defunctWellNames.end()) { - continue; - } + for (const auto& wname: schedule.wellNames(reportStepNum)) { + + // don't bother with wells not on this process + const auto& defunctWellNames = simulator_.vanguard().defunctWellNames(); + if (defunctWellNames.find(wname) != defunctWellNames.end()) { + continue; + } - const auto& well = schedule.getWell2(wname, reportStepNum); - - // Ignore injector wells - if (well.isInjector()){ - continue; - } - - tmp_names[0] = wname;//WellProdDataType::WellName - - - auto wName = static_cast(wname); - auto get = [&st, &wName](const std::string& vector) - { - const auto key = vector + ':' + wName; + const auto& well = schedule.getWell2(wname, reportStepNum); + + // Ignore injector wells + if (well.isInjector()){ + continue; + } + + tmp_names[0] = wname;//WellProdDataType::WellName + + + auto wName = static_cast(wname); + auto get = [&st, &wName](const std::string& vector) + { + const auto key = vector + ':' + wName; - return st.has(key) ? st.get(key) : 0.0; - }; - - const auto& controls = well.productionControls(st); - using CMode = ::Opm::Well2::ProducerCMode; - - auto fctl = [](const auto wmctl) -> std::string - { - switch (wmctl) { - case CMode::ORAT: return "ORAT"; - case CMode::WRAT: return "WRAT"; - case CMode::GRAT: return "GRAT"; - case CMode::LRAT: return "LRAT"; - case CMode::RESV: return "RESV"; - case CMode::THP: return "THP"; - case CMode::BHP: return "BHP"; - case CMode::CRAT: return "CRate"; - case CMode::GRUP: return "GRUP"; - default: - { - return "none"; - } - } - }; + return st.has(key) ? st.get(key) : 0.0; + }; + + const auto& controls = well.productionControls(st); + using CMode = ::Opm::Well2::ProducerCMode; + + auto fctl = [](const auto wmctl) -> std::string + { + switch (wmctl) { + case CMode::ORAT: return "ORAT"; + case CMode::WRAT: return "WRAT"; + case CMode::GRAT: return "GRAT"; + case CMode::LRAT: return "LRAT"; + case CMode::RESV: return "RESV"; + case CMode::THP: return "THP"; + case CMode::BHP: return "BHP"; + case CMode::CRAT: return "CRate"; + case CMode::GRUP: return "GRUP"; + default: + { + return "none"; + } + } + }; - tmp_names[1] = fctl(controls.cmode); //WellProdDataType::CTRLMode - - tmp_values[0] = well.getHeadI() + 1;//WellProdDataType::WellLocationi - tmp_values[1] = well.getHeadJ() + 1;//WellProdDataType::WellLocationj - tmp_values[2] = get("WOPR"); //WellProdDataType::OilRate - tmp_values[3] = get("WWPR"); //WellProdDataType::WaterRate - tmp_values[4] = get("WGPR"); //WellProdDataType::GasRate - tmp_values[5] = get("WVPR"); //WellProdDataType::FluidResVol - tmp_values[6] = get("WWCT"); //WellProdDataType::WaterCut - tmp_values[7] = get("WGOR"); //WellProdDataType::GasOilRatio - tmp_values[8] = get("WWPR")/get("WGPR"); //WellProdDataType::WaterGasRatio - tmp_values[9] = get("WBHP"); //WellProdDataType::BHP - tmp_values[10] = get("WTHP"); //WellProdDataType::THP - //tmp_values[11] = 0; //WellProdDataType::SteadyStatePI // - - outputProductionReport_(tmp_values, tmp_names, forceDisableProdOutput); - - } - } + tmp_names[1] = fctl(controls.cmode); //WellProdDataType::CTRLMode + + tmp_values[0] = well.getHeadI() + 1;//WellProdDataType::WellLocationi + tmp_values[1] = well.getHeadJ() + 1;//WellProdDataType::WellLocationj + tmp_values[2] = get("WOPR"); //WellProdDataType::OilRate + tmp_values[3] = get("WWPR"); //WellProdDataType::WaterRate + tmp_values[4] = get("WGPR"); //WellProdDataType::GasRate + tmp_values[5] = get("WVPR"); //WellProdDataType::FluidResVol + tmp_values[6] = get("WWCT"); //WellProdDataType::WaterCut + tmp_values[7] = get("WGOR"); //WellProdDataType::GasOilRatio + tmp_values[8] = get("WWPR")/get("WGPR"); //WellProdDataType::WaterGasRatio + tmp_values[9] = get("WBHP"); //WellProdDataType::BHP + tmp_values[10] = get("WTHP"); //WellProdDataType::THP + //tmp_values[11] = 0; //WellProdDataType::SteadyStatePI // + + outputProductionReport_(tmp_values, tmp_names, forceDisableProdOutput); + + } + } } // write injection report to output void outputInjLog(size_t reportStepNum, const bool substep, bool forceDisableInjOutput) { - if (!substep) { - ScalarBuffer tmp_values(WellInjDataType::numWIValues, 0.0); - StringBuffer tmp_names(WellInjDataType::numWINames, ""); - outputInjectionReport_(tmp_values, tmp_names, forceDisableInjOutput); - - const auto& st = simulator_.vanguard().summaryState(); - const auto& schedule = simulator_.vanguard().schedule(); - for (const auto& gname: schedule.groupNames()) { + if (!substep) { + ScalarBuffer tmp_values(WellInjDataType::numWIValues, 0.0); + StringBuffer tmp_names(WellInjDataType::numWINames, ""); + outputInjectionReport_(tmp_values, tmp_names, forceDisableInjOutput); + + const auto& st = simulator_.vanguard().summaryState(); + const auto& schedule = simulator_.vanguard().schedule(); + for (const auto& gname: schedule.groupNames()) { - auto gName = static_cast(gname); - auto get = [&st, &gName](const std::string& vector) - { - const auto key = vector + ':' + gName; + auto gName = static_cast(gname); + auto get = [&st, &gName](const std::string& vector) + { + const auto key = vector + ':' + gName; - return st.has(key) ? st.get(key) : 0.0; - }; + return st.has(key) ? st.get(key) : 0.0; + }; - tmp_names[0] = gname; + tmp_names[0] = gname; - tmp_values[2] = get("GOIR");//WellInjDataType::OilRate - tmp_values[3] = get("GWIR"); //WellInjDataType::WaterRate - tmp_values[4] = get("GGIR"); //WellInjDataType::GasRate - tmp_values[5] = get("GVIR");//WellInjDataType::FluidResVol + tmp_values[2] = get("GOIR");//WellInjDataType::OilRate + tmp_values[3] = get("GWIR"); //WellInjDataType::WaterRate + tmp_values[4] = get("GGIR"); //WellInjDataType::GasRate + tmp_values[5] = get("GVIR");//WellInjDataType::FluidResVol - outputInjectionReport_(tmp_values, tmp_names, forceDisableInjOutput); - } + outputInjectionReport_(tmp_values, tmp_names, forceDisableInjOutput); + } - for (const auto& wname: schedule.wellNames(reportStepNum)) { - - // don't bother with wells not on this process - const auto& defunctWellNames = simulator_.vanguard().defunctWellNames(); - if (defunctWellNames.find(wname) != defunctWellNames.end()) { - continue; - } + for (const auto& wname: schedule.wellNames(reportStepNum)) { + + // don't bother with wells not on this process + const auto& defunctWellNames = simulator_.vanguard().defunctWellNames(); + if (defunctWellNames.find(wname) != defunctWellNames.end()) { + continue; + } - const auto& well = schedule.getWell2(wname, reportStepNum); - - // Ignore Producer wells - if (well.isProducer()){ - continue; - } + const auto& well = schedule.getWell2(wname, reportStepNum); + + // Ignore Producer wells + if (well.isProducer()){ + continue; + } - tmp_names[0] = wname; //WellInjDataType::WellName + tmp_names[0] = wname; //WellInjDataType::WellName - auto wName = static_cast(wname); - auto get = [&st, &wName](const std::string& vector) - { - const auto key = vector + ':' + wName; + auto wName = static_cast(wname); + auto get = [&st, &wName](const std::string& vector) + { + const auto key = vector + ':' + wName; - return st.has(key) ? st.get(key) : 0.0; - }; - - const auto& controls = well.injectionControls(st); - const auto ctlMode = controls.cmode; - const auto injType = controls.injector_type; - using CMode = ::Opm::Well2::InjectorCMode; - using WType = ::Opm::Well2::InjectorType; - - auto ftype = [](const auto wtype) -> std::string - { - switch (wtype) { - case WType::OIL: return "Oil"; - case WType::WATER: return "Wat"; - case WType::GAS: return "Gas"; - case WType::MULTI: return "Multi"; - - default: - { - return ""; - } - } - }; - - auto fctl = [](const auto wmctl) -> std::string - { - switch (wmctl) { - case CMode::RATE: return "RATE"; - case CMode::RESV: return "RESV"; - case CMode::THP: return "THP"; - case CMode::BHP: return "BHP"; - case CMode::GRUP: return "GRUP"; - - default: - { - return ""; - } - } - }; + return st.has(key) ? st.get(key) : 0.0; + }; + + const auto& controls = well.injectionControls(st); + const auto ctlMode = controls.cmode; + const auto injType = controls.injector_type; + using CMode = ::Opm::Well2::InjectorCMode; + using WType = ::Opm::Well2::InjectorType; + + auto ftype = [](const auto wtype) -> std::string + { + switch (wtype) { + case WType::OIL: return "Oil"; + case WType::WATER: return "Wat"; + case WType::GAS: return "Gas"; + case WType::MULTI: return "Multi"; + default: + { + return ""; + } + } + }; + + auto fctl = [](const auto wmctl) -> std::string + { + switch (wmctl) { + case CMode::RATE: return "RATE"; + case CMode::RESV: return "RESV"; + case CMode::THP: return "THP"; + case CMode::BHP: return "BHP"; + case CMode::GRUP: return "GRUP"; + default: + { + return ""; + } + } + }; - const auto flowtype = ftype(injType); - const auto flowctl = fctl(ctlMode); - if(flowtype == "Oil") //WellInjDataType::CTRLModeOil - { - if (flowctl == "RATE"){ tmp_names[1] = "ORAT"; } - else { tmp_names[1] = flowctl; } - } - else if (flowtype == "Wat") //WellInjDataType::CTRLModeWat - { - if (flowctl == "RATE"){ tmp_names[3] = "WRAT"; } - else { tmp_names[2] = flowctl; } - } - else if (flowtype == "Gas") //WellInjDataType::CTRLModeGas - { - if (flowctl == "RATE"){ tmp_names[3] = "GRAT"; } - else { tmp_names[3] = flowctl; } - } + const auto flowtype = ftype(injType); + const auto flowctl = fctl(ctlMode); + if(flowtype == "Oil") //WellInjDataType::CTRLModeOil + { + if (flowctl == "RATE"){ tmp_names[1] = "ORAT"; } + else { tmp_names[1] = flowctl; } + } + else if (flowtype == "Wat") //WellInjDataType::CTRLModeWat + { + if (flowctl == "RATE"){ tmp_names[3] = "WRAT"; } + else { tmp_names[2] = flowctl; } + } + else if (flowtype == "Gas") //WellInjDataType::CTRLModeGas + { + if (flowctl == "RATE"){ tmp_names[3] = "GRAT"; } + else { tmp_names[3] = flowctl; } + } - tmp_values[0] = well.getHeadI() + 1; //WellInjDataType::wellLocationi - tmp_values[1] = well.getHeadJ() + 1; //WellInjDataType::wellLocationj - tmp_values[2] = get("WOIR"); //WellInjDataType::OilRate - tmp_values[3] = get("WWIR"); //WellInjDataType::WaterRate - tmp_values[4] = get("WGIR"); //WellInjDataType::GasRate - tmp_values[5] = get("WVIR");//WellInjDataType::FluidResVol - tmp_values[6] = get("WBHP"); //WellInjDataType::BHP - tmp_values[7] = get("WTHP"); //WellInjDataType::THP - //tmp_values[8] = 0; //WellInjDataType::SteadyStateII - - outputInjectionReport_(tmp_values, tmp_names, forceDisableInjOutput); - - } - } + tmp_values[0] = well.getHeadI() + 1; //WellInjDataType::wellLocationi + tmp_values[1] = well.getHeadJ() + 1; //WellInjDataType::wellLocationj + tmp_values[2] = get("WOIR"); //WellInjDataType::OilRate + tmp_values[3] = get("WWIR"); //WellInjDataType::WaterRate + tmp_values[4] = get("WGIR"); //WellInjDataType::GasRate + tmp_values[5] = get("WVIR");//WellInjDataType::FluidResVol + tmp_values[6] = get("WBHP"); //WellInjDataType::BHP + tmp_values[7] = get("WTHP"); //WellInjDataType::THP + //tmp_values[8] = 0; //WellInjDataType::SteadyStateII + + outputInjectionReport_(tmp_values, tmp_names, forceDisableInjOutput); + + } + } } // write cumulative production and injection reports to output - void outputCumLog(size_t reportStepNum, const bool substep, bool forceDisableCumOutput) + void outputCumLog(size_t reportStepNum, const bool substep, bool forceDisableCumOutput) { - if (!substep) { - ScalarBuffer tmp_values(WellCumDataType::numWCValues, 0.0); - StringBuffer tmp_names(WellCumDataType::numWCNames, ""); - outputCumulativeReport_(tmp_values, tmp_names, forceDisableCumOutput); - - const auto& st = simulator_.vanguard().summaryState(); - const auto& schedule = simulator_.vanguard().schedule(); + if (!substep) { + ScalarBuffer tmp_values(WellCumDataType::numWCValues, 0.0); + StringBuffer tmp_names(WellCumDataType::numWCNames, ""); + outputCumulativeReport_(tmp_values, tmp_names, forceDisableCumOutput); + + const auto& st = simulator_.vanguard().summaryState(); + const auto& schedule = simulator_.vanguard().schedule(); - for (const auto& gname: schedule.groupNames()) { + for (const auto& gname: schedule.groupNames()) { - auto gName = static_cast(gname); - auto get = [&st, &gName](const std::string& vector) - { - const auto key = vector + ':' + gName; + auto gName = static_cast(gname); + auto get = [&st, &gName](const std::string& vector) + { + const auto key = vector + ':' + gName; - return st.has(key) ? st.get(key) : 0.0; - }; + return st.has(key) ? st.get(key) : 0.0; + }; + + tmp_names[0] = gname; - tmp_names[0] = gname; + tmp_values[2] = get("GOPT"); //WellCumDataType::OilProd + tmp_values[3] = get("GWPT"); //WellCumDataType::WaterProd + tmp_values[4] = get("GGPT"); //WellCumDataType::GasProd + tmp_values[5] = get("GVPT");//WellCumDataType::FluidResVolProd + tmp_values[2] = get("GOIT"); //WellCumDataType::OilInj + tmp_values[3] = get("GWIT"); //WellCumDataType::WaterInj + tmp_values[4] = get("GGIT"); //WellCumDataType::GasInj + tmp_values[5] = get("GVIT");//WellCumDataType::FluidResVolInj - tmp_values[2] = get("GOPT"); //WellCumDataType::OilProd - tmp_values[3] = get("GWPT"); //WellCumDataType::WaterProd - tmp_values[4] = get("GGPT"); //WellCumDataType::GasProd - tmp_values[5] = get("GVPT");//WellCumDataType::FluidResVolProd - tmp_values[2] = get("GOIT"); //WellCumDataType::OilInj - tmp_values[3] = get("GWIT"); //WellCumDataType::WaterInj - tmp_values[4] = get("GGIT"); //WellCumDataType::GasInj - tmp_values[5] = get("GVIT");//WellCumDataType::FluidResVolInj - - outputCumulativeReport_(tmp_values, tmp_names, forceDisableCumOutput); - } + outputCumulativeReport_(tmp_values, tmp_names, forceDisableCumOutput); + } - for (const auto& wname : schedule.wellNames(reportStepNum)) { - - // don't bother with wells not on this process - const auto& defunctWellNames = simulator_.vanguard().defunctWellNames(); - if (defunctWellNames.find(wname) != defunctWellNames.end()) { - continue; - } + for (const auto& wname : schedule.wellNames(reportStepNum)) { + + // don't bother with wells not on this process + const auto& defunctWellNames = simulator_.vanguard().defunctWellNames(); + if (defunctWellNames.find(wname) != defunctWellNames.end()) { + continue; + } - const auto& well = schedule.getWell2(wname, reportStepNum); + const auto& well = schedule.getWell2(wname, reportStepNum); - tmp_names[0] = wname; //WellCumDataType::WellName + tmp_names[0] = wname; //WellCumDataType::WellName - auto wName = static_cast(wname); - auto get = [&st, &wName](const std::string& vector) - { - const auto key = vector + ':' + wName; + auto wName = static_cast(wname); + auto get = [&st, &wName](const std::string& vector) + { + const auto key = vector + ':' + wName; - return st.has(key) ? st.get(key) : 0.0; - }; - - if (well.isInjector()) { - - const auto& controls = well.injectionControls(st); - const auto ctlMode = controls.cmode; - const auto injType = controls.injector_type; - using CMode = ::Opm::Well2::InjectorCMode; - using WType = ::Opm::Well2::InjectorType; - - auto ftype = [](const auto wtype) -> std::string - { - switch (wtype) { - case WType::OIL: return "Oil"; - case WType::WATER: return "Wat"; - case WType::GAS: return "Gas"; - case WType::MULTI: return "Multi"; - - default: - { - return ""; - } - } - }; - - auto fctl = [](const auto wmctl) -> std::string - { - switch (wmctl) { - case CMode::RATE: return "RATE"; - case CMode::RESV: return "RESV"; - case CMode::THP: return "THP"; - case CMode::BHP: return "BHP"; - case CMode::GRUP: return "GRUP"; + return st.has(key) ? st.get(key) : 0.0; + }; + + if (well.isInjector()) { + + const auto& controls = well.injectionControls(st); + const auto ctlMode = controls.cmode; + const auto injType = controls.injector_type; + using CMode = ::Opm::Well2::InjectorCMode; + using WType = ::Opm::Well2::InjectorType; + + auto ftype = [](const auto wtype) -> std::string + { + switch (wtype) { + case WType::OIL: return "Oil"; + case WType::WATER: return "Wat"; + case WType::GAS: return "Gas"; + case WType::MULTI: return "Multi"; + default: + { + return ""; + } + } + }; + + auto fctl = [](const auto wmctl) -> std::string + { + switch (wmctl) { + case CMode::RATE: return "RATE"; + case CMode::RESV: return "RESV"; + case CMode::THP: return "THP"; + case CMode::BHP: return "BHP"; + case CMode::GRUP: return "GRUP"; + default: + { + return ""; + } + } + }; - default: - { - return ""; - } - } - }; - - tmp_names[1] = "INJ"; //WellCumDataType::WellType - const auto flowctl = fctl(ctlMode); - if (flowctl == "RATE") //WellCumDataType::WellCTRL - { - const auto flowtype = ftype(injType); - if(flowtype == "Oil"){ tmp_names[2] = "ORAT"; } - else if(flowtype == "Wat"){ tmp_names[2] = "WRAT"; } - else if(flowtype == "Gas"){ tmp_names[2] = "GRAT"; } - } - else - { - tmp_names[2] = flowctl; - } + tmp_names[1] = "INJ"; //WellCumDataType::WellType + const auto flowctl = fctl(ctlMode); + if (flowctl == "RATE") //WellCumDataType::WellCTRL + { + const auto flowtype = ftype(injType); + if(flowtype == "Oil"){ tmp_names[2] = "ORAT"; } + else if(flowtype == "Wat"){ tmp_names[2] = "WRAT"; } + else if(flowtype == "Gas"){ tmp_names[2] = "GRAT"; } + } + else + { + tmp_names[2] = flowctl; + } - } - else if (well.isProducer()) { - - const auto& controls = well.productionControls(st); - using CMode = ::Opm::Well2::ProducerCMode; - - auto fctl = [](const auto wmctl) -> std::string - { - switch (wmctl) { - case CMode::ORAT: return "ORAT"; - case CMode::WRAT: return "WRAT"; - case CMode::GRAT: return "GRAT"; - case CMode::LRAT: return "LRAT"; - case CMode::RESV: return "RESV"; - case CMode::THP: return "THP"; - case CMode::BHP: return "BHP"; - case CMode::CRAT: return "CRAT"; - case CMode::GRUP: return "GRUP"; - default: - { - return "none"; - } - } - }; - tmp_names[1] = "PROD"; //WellProdDataType::CTRLMode - tmp_names[2] = fctl(controls.cmode); //WellProdDataType::CTRLMode + } + else if (well.isProducer()) { + + const auto& controls = well.productionControls(st); + using CMode = ::Opm::Well2::ProducerCMode; + + auto fctl = [](const auto wmctl) -> std::string + { + switch (wmctl) { + case CMode::ORAT: return "ORAT"; + case CMode::WRAT: return "WRAT"; + case CMode::GRAT: return "GRAT"; + case CMode::LRAT: return "LRAT"; + case CMode::RESV: return "RESV"; + case CMode::THP: return "THP"; + case CMode::BHP: return "BHP"; + case CMode::CRAT: return "CRAT"; + case CMode::GRUP: return "GRUP"; + default: + { + return "none"; + } + } + }; + tmp_names[1] = "PROD"; //WellProdDataType::CTRLMode + tmp_names[2] = fctl(controls.cmode); //WellProdDataType::CTRLMode + } + + tmp_values[0] = well.getHeadI() + 1; //WellCumDataType::wellLocationi + tmp_values[1] = well.getHeadJ() + 1; //WellCumDataType::wellLocationj + tmp_values[2] = get("WOPT"); //WellCumDataType::OilProd + tmp_values[3] = get("WWPT"); //WellCumDataType::WaterProd + tmp_values[4] = get("WGPT"); //WellCumDataType::GasProd + tmp_values[5] = get("WVPT");//WellCumDataType::FluidResVolProd + tmp_values[2] = get("WOIT"); //WellCumDataType::OilInj + tmp_values[3] = get("WWIT"); //WellCumDataType::WaterInj + tmp_values[4] = get("WGIT"); //WellCumDataType::GasInj + tmp_values[5] = get("WVIT");//WellCumDataType::FluidResVolInj + + outputCumulativeReport_(tmp_values, tmp_names, forceDisableCumOutput); + + } } - - tmp_values[0] = well.getHeadI() + 1; //WellCumDataType::wellLocationi - tmp_values[1] = well.getHeadJ() + 1; //WellCumDataType::wellLocationj - tmp_values[2] = get("WOPT"); //WellCumDataType::OilProd - tmp_values[3] = get("WWPT"); //WellCumDataType::WaterProd - tmp_values[4] = get("WGPT"); //WellCumDataType::GasProd - tmp_values[5] = get("WVPT");//WellCumDataType::FluidResVolProd - tmp_values[2] = get("WOIT"); //WellCumDataType::OilInj - tmp_values[3] = get("WWIT"); //WellCumDataType::WaterInj - tmp_values[4] = get("WGIT"); //WellCumDataType::GasInj - tmp_values[5] = get("WVIT");//WellCumDataType::FluidResVolInj - - outputCumulativeReport_(tmp_values, tmp_names, forceDisableCumOutput); - - } - } } void setRestart(const Opm::data::Solution& sol, unsigned elemIdx, unsigned globalDofIndex) @@ -1902,9 +1898,9 @@ private: void outputProductionReport_(const ScalarBuffer& wellProd, const StringBuffer& wellProdNames, const bool forceDisableProdOutput) { - if(forceDisableProdOutput) - return; - + if(forceDisableProdOutput) + return; + const Opm::UnitSystem& units = simulator_.vanguard().eclState().getUnits(); std::ostringstream ss; if (wellProdNames[WellProdDataType::WellName].empty()) { @@ -1912,17 +1908,17 @@ private: << ": WELL : LOCATION :CTRL: OIL : WATER : GAS : FLUID : WATER : GAS/OIL : WAT/GAS : BHP OR : THP OR :\n"// STEADY-ST PI :\n" << ": NAME : (I,J,K) :MODE: RATE : RATE : RATE : RES.VOL. : CUT : RATIO : RATIO : CON.PR.: BLK.PR.:\n";// OR POTN OF PREF. PH:\n"; if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC) { - ss << ": : : : SCM/DAY : SCM/DAY : SCM/DAY : RCM/DAY : SCM/SCM : SCM/SCM : SCM/SCM : BARSA : BARSA :\n";// :\n"; - } - if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { - ss << ": : : : STB/DAY : STB/DAY : MSCF/DAY : RB/DAY : : MSCF/STB : STB/MSCF : PSIA : PSIA :\n";// :\n"; - } - if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_LAB) { - ss << ": : : : SCC/HR : SCC/HR : SCC/HR : RCC : SCC/SCC : SCC/SCC : SCC/SCC : ATMA : ATMA :\n";// :\n"; - } - ss << "=================================================================================================================================\n";//=================== \n"; + ss << ": : : : SCM/DAY : SCM/DAY : SCM/DAY : RCM/DAY : SCM/SCM : SCM/SCM : SCM/SCM : BARSA : BARSA :\n";// :\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { + ss << ": : : : STB/DAY : STB/DAY : MSCF/DAY : RB/DAY : : MSCF/STB : STB/MSCF : PSIA : PSIA :\n";// :\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_LAB) { + ss << ": : : : SCC/HR : SCC/HR : SCC/HR : RCC : SCC/SCC : SCC/SCC : SCC/SCC : ATMA : ATMA :\n";// :\n"; + } + ss << "=================================================================================================================================\n";//=================== \n"; } - else { + else { if (wellProd[WellProdDataType::WellLocationi] < 1) { ss << std::right << std::fixed << ":" << std::setw (8) << wellProdNames[WellProdDataType::WellName] << ":" << std::setprecision(0) << std::setw(11) << "" << ":" << std::setw(4) << wellProdNames[WellProdDataType::CTRLMode] << ":" << std::setprecision(1) << std::setw(11) << wellProd[WellProdDataType::OilRate] << ":" << std::setw(11) << wellProd[WellProdDataType::WaterRate] << ":" << std::setw(11)<< wellProd[WellProdDataType::GasRate] << ":" << std::setw(11) << wellProd[WellProdDataType::FluidResVol] << std::setprecision(3) << ":" << std::setw(11) << wellProd[WellProdDataType::WaterCut] << std::setprecision(2) << ":" << std::setw(10) << wellProd[WellProdDataType::GasOilRatio] << std::setprecision(4) << ":" << std::setw(12) << wellProd[WellProdDataType::WatGasRatio] << std::setprecision(1) << ":" << std::setw(8) << "" << ":" << std::setw(8) << "" << ": \n";//wellProd[WellProdDataType::SteadyStatePI] << std::setw(10) << "\n" } @@ -1931,14 +1927,14 @@ private: } ss << ":"<< std::setfill ('-') << std::setw (9) << ":" << std::setfill ('-') << std::setw (12) << ":" << std::setfill ('-') << std::setw (5) << ":" << std::setfill ('-') << std::setw (12) << ":" << std::setfill ('-') << std::setw (12) << ":" << std::setfill ('-') << std::setw (12) << ":" << std::setfill ('-') << std::setw (12) << ":" << std::setfill ('-') << std::setw (12) << ":" << std::setfill ('-') << std::setw (11) << ":" << std::setfill ('-') << std::setw (13) << ":" << std::setfill ('-') << std::setw (9) << ":" << std::setfill ('-') << std::setw (9) << ":" << "\n"; } - Opm::OpmLog::note(ss.str()); - } + Opm::OpmLog::note(ss.str()); + } void outputInjectionReport_(const ScalarBuffer& wellInj, const StringBuffer& wellInjNames, const bool forceDisableInjOutput) { - if(forceDisableInjOutput) - return; - + if(forceDisableInjOutput) + return; + const Opm::UnitSystem& units = simulator_.vanguard().eclState().getUnits(); std::ostringstream ss; if (wellInjNames[WellInjDataType::WellName].empty()) { @@ -1946,17 +1942,17 @@ private: << ": WELL : LOCATION : CTRL : CTRL : CTRL : OIL : WATER : GAS : FLUID : BHP OR : THP OR :\n"// STEADY-ST II :\n" << ": NAME : (I,J,K) : MODE : MODE : MODE : RATE : RATE : RATE : RES.VOL. : CON.PR.: BLK.PR.:\n";// OR POTENTIAL :\n"; if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC) { - ss << ": : : OIL : WAT : GAS : SCM/DAY : SCM/DAY : SCM/DAY : RCM/DAY : BARSA : BARSA :\n";// :\n"; - } - if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { - ss << ": : : OIL : WAT : GAS : STB/DAY : STB/DAY : MSCF/DAY : RB/DAY : PSIA : PSIA :\n";// :\n"; - } - if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_LAB) { - ss << ": : : OIL : WAT : GAS : SCC/HR : SCC/HR : SCC/HR : RCC/HR : ATMA : ATMA :\n";// :\n"; - } - ss << "==============================================================================================================\n";//===================== \n"; + ss << ": : : OIL : WAT : GAS : SCM/DAY : SCM/DAY : SCM/DAY : RCM/DAY : BARSA : BARSA :\n";// :\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { + ss << ": : : OIL : WAT : GAS : STB/DAY : STB/DAY : MSCF/DAY : RB/DAY : PSIA : PSIA :\n";// :\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_LAB) { + ss << ": : : OIL : WAT : GAS : SCC/HR : SCC/HR : SCC/HR : RCC/HR : ATMA : ATMA :\n";// :\n"; + } + ss << "==============================================================================================================\n";//===================== \n"; } - else { + else { if (wellInj[WellInjDataType::WellLocationi] < 1) { ss << std::right << std::fixed << std::setprecision(0) << ":" << std::setw (8) << wellInjNames[WellInjDataType::WellName] << ":" << std::setw(11) << "" << ":" << std::setw(6) << wellInjNames[WellInjDataType::CTRLModeOil] << ":" << std::setw(6) << wellInjNames[WellInjDataType::CTRLModeWat] << ":" << std::setw(6) << wellInjNames[WellInjDataType::CTRLModeGas] << ":" << std::setprecision(1) << std::setw(11) << wellInj[WellInjDataType::OilRate] << ":" << std::setw(11) << wellInj[WellInjDataType::WaterRate] << ":" << std::setw(11)<< wellInj[WellInjDataType::GasRate] << ":" << std::setw(11) << wellInj[WellInjDataType::FluidResVol] << ":" << std::setw(8)<< "" << ":" << std::setw(8)<< "" << ": \n";//wellInj[WellInjDataType::SteadyStateII] << std::setw(10) << "\n" } @@ -1965,14 +1961,14 @@ private: } ss << ":--------:-----------:------:------:------:------------:----------:-----------:-----------:--------:--------: \n";//--------------------:\n"; } - Opm::OpmLog::note(ss.str()); + Opm::OpmLog::note(ss.str()); } - - void outputCumulativeReport_(const ScalarBuffer& wellCum, const StringBuffer& wellCumNames, const bool forceDisableCumOutput) + + void outputCumulativeReport_(const ScalarBuffer& wellCum, const StringBuffer& wellCumNames, const bool forceDisableCumOutput) { - if(forceDisableCumOutput) - return; - + if(forceDisableCumOutput) + return; + const Opm::UnitSystem& units = simulator_.vanguard().eclState().getUnits(); std::ostringstream ss; if (wellCumNames[WellCumDataType::WellName].empty()) { @@ -1980,17 +1976,17 @@ private: << ": WELL : LOCATION : WELL :CTRL: OIL : WATER : GAS : Prod : OIL : WATER : GAS : INJ :\n" << ": NAME : (I,J,K) : TYPE :MODE: PROD : PROD : PROD : RES.VOL. : INJ : INJ : INJ : RES.VOL. :\n"; if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_METRIC) { - ss << ": : : : : MSCM : MSCM : MMSCM : MRCM : MSCM : MSCM : MMSCM : MRCM :\n"; - } - if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { - ss << ": : : : : MSTB : MSTB : MMSCF : MRB : MSTB : MSTB : MMSCF : MRB :\n"; - } - if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_LAB) { - ss << ": : : : : MSCC : MSCC : MMSCC : MRCC : MSCC : MSCC : MMSCC : MRCC :\n"; - } - ss << "====================================================================================================================================\n"; + ss << ": : : : : MSCM : MSCM : MMSCM : MRCM : MSCM : MSCM : MMSCM : MRCM :\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_FIELD) { + ss << ": : : : : MSTB : MSTB : MMSCF : MRB : MSTB : MSTB : MMSCF : MRB :\n"; + } + if (units.getType() == Opm::UnitSystem::UnitType::UNIT_TYPE_LAB) { + ss << ": : : : : MSCC : MSCC : MMSCC : MRCC : MSCC : MSCC : MMSCC : MRCC :\n"; + } + ss << "====================================================================================================================================\n"; } - else { + else { if (wellCum[WellCumDataType::WellLocationi] < 1) { ss << std::right << std::fixed << std::setprecision(0) << ":" << std::setw (8) << wellCumNames[WellCumDataType::WellName] << ":" << std::setw(11) << "" << ":" << std::setw(8) << wellCumNames[WellCumDataType::WellType] << ":" << std::setw(4) << wellCumNames[WellCumDataType::WellCTRL] << ":" << std::setprecision(1) << std::setw(11) << wellCum[WellCumDataType::OilProd] << ":" << std::setw(11) << wellCum[WellCumDataType::WaterProd] << ":" << std::setw(11)<< wellCum[WellCumDataType::GasProd]/1000 << ":" << std::setw(11) << wellCum[WellCumDataType::FluidResVolProd] << ":" << std::setw(11) << wellCum[WellCumDataType::OilProd] << ":" << std::setw(11) << wellCum[WellCumDataType::WaterInj] << ":" << std::setw(11) << wellCum[WellCumDataType::GasInj]/1000 << ":" << std::setw(11) << wellCum[WellCumDataType::FluidResVolInj] << ": \n"; } @@ -1999,7 +1995,7 @@ private: } ss << ":--------:-----------:--------:----:------------:----------:-----------:-----------:------------:----------:-----------:-----------: \n"; } - Opm::OpmLog::note(ss.str()); + Opm::OpmLog::note(ss.str()); } std::string fipEnumToString_(int i) @@ -2023,62 +2019,62 @@ private: typedef typename WellProdDataType::WPId WPId; switch(static_cast(i)) { - case WellProdDataType::WellName: return "WName"; - case WellProdDataType::WellLocationi: return "WLi"; - case WellProdDataType::WellLocationj: return "WLj"; - case WellProdDataType::CTRLMode: return "CTRL"; - case WellProdDataType::OilRate: return "OR"; - case WellProdDataType::WaterRate: return "WR"; - case WellProdDataType::GasRate: return "GR"; - case WellProdDataType::FluidResVol: return "FRV"; - case WellProdDataType::WaterCut: return "WC"; - case WellProdDataType::GasOilRatio: return "GOR"; - case WellProdDataType::WatGasRatio: return "WGR"; - case WellProdDataType::BHP: return "BHP"; - case WellProdDataType::THP: return "THP"; - case WellProdDataType::SteadyStatePI: return "SteadyStatePI"; + case WellProdDataType::WellName: return "WName"; + case WellProdDataType::WellLocationi: return "WLi"; + case WellProdDataType::WellLocationj: return "WLj"; + case WellProdDataType::CTRLMode: return "CTRL"; + case WellProdDataType::OilRate: return "OR"; + case WellProdDataType::WaterRate: return "WR"; + case WellProdDataType::GasRate: return "GR"; + case WellProdDataType::FluidResVol: return "FRV"; + case WellProdDataType::WaterCut: return "WC"; + case WellProdDataType::GasOilRatio: return "GOR"; + case WellProdDataType::WatGasRatio: return "WGR"; + case WellProdDataType::BHP: return "BHP"; + case WellProdDataType::THP: return "THP"; + case WellProdDataType::SteadyStatePI: return "SteadyStatePI"; } return "ERROR"; - } - std::string WIEnumToString_(int i) + } + std::string WIEnumToString_(int i) { typedef typename WellInjDataType::WIId WIId; switch(static_cast(i)) { - case WellInjDataType::WellName: return "WName"; - case WellInjDataType::WellLocationi: return "WLi"; - case WellInjDataType::WellLocationj: return "WLj"; - case WellInjDataType::CTRLModeOil: return "CTRLo"; - case WellInjDataType::CTRLModeWat: return "CTRLw"; - case WellInjDataType::CTRLModeGas: return "CTRLg"; - case WellInjDataType::OilRate: return "OR"; - case WellInjDataType::WaterRate: return "WR"; - case WellInjDataType::GasRate: return "GR"; - case WellInjDataType::FluidResVol: return "FRV"; - case WellInjDataType::BHP: return "BHP"; - case WellInjDataType::THP: return "THP"; - case WellInjDataType::SteadyStateII: return "SteadyStateII"; + case WellInjDataType::WellName: return "WName"; + case WellInjDataType::WellLocationi: return "WLi"; + case WellInjDataType::WellLocationj: return "WLj"; + case WellInjDataType::CTRLModeOil: return "CTRLo"; + case WellInjDataType::CTRLModeWat: return "CTRLw"; + case WellInjDataType::CTRLModeGas: return "CTRLg"; + case WellInjDataType::OilRate: return "OR"; + case WellInjDataType::WaterRate: return "WR"; + case WellInjDataType::GasRate: return "GR"; + case WellInjDataType::FluidResVol: return "FRV"; + case WellInjDataType::BHP: return "BHP"; + case WellInjDataType::THP: return "THP"; + case WellInjDataType::SteadyStateII: return "SteadyStateII"; } return "ERROR"; } - std::string WCEnumToString_(int i) + std::string WCEnumToString_(int i) { typedef typename WellCumDataType::WCId WCId; switch(static_cast(i)) { - case WellCumDataType::WellName: return "WName"; - case WellCumDataType::WellLocationi: return "WLi"; - case WellCumDataType::WellLocationj: return "WLj"; - case WellCumDataType::WellType: return "WType"; - case WellCumDataType::WellCTRL: return "WCTRL"; - case WellCumDataType::OilProd: return "OP"; - case WellCumDataType::WaterProd: return "WP"; - case WellCumDataType::GasProd: return "GP"; - case WellCumDataType::FluidResVolProd: return "FRVP"; - case WellCumDataType::OilInj: return "OI"; - case WellCumDataType::WaterInj: return "WI"; - case WellCumDataType::GasInj: return "GI"; - case WellCumDataType::FluidResVolInj: return "FRVI"; + case WellCumDataType::WellName: return "WName"; + case WellCumDataType::WellLocationi: return "WLi"; + case WellCumDataType::WellLocationj: return "WLj"; + case WellCumDataType::WellType: return "WType"; + case WellCumDataType::WellCTRL: return "WCTRL"; + case WellCumDataType::OilProd: return "OP"; + case WellCumDataType::WaterProd: return "WP"; + case WellCumDataType::GasProd: return "GP"; + case WellCumDataType::FluidResVolProd: return "FRVP"; + case WellCumDataType::OilInj: return "OI"; + case WellCumDataType::WaterInj: return "WI"; + case WellCumDataType::GasInj: return "GI"; + case WellCumDataType::FluidResVolInj: return "FRVI"; } return "ERROR"; } diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index a9c289b94..ecc05d994 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -43,19 +43,19 @@ #include #include +#include #include #include -#include -#include #include +#include #include -#include #include -#include #include +#include +#include #ifdef HAVE_MPI #include @@ -71,9 +71,11 @@ END_PROPERTIES namespace Opm { -template class EclWriter; +template +class EclWriter; -template class EclOutputBlackOilModule; +template +class EclOutputBlackOilModule; /*! * \brief Detect whether two cells are direct vertical neighbours. @@ -84,45 +86,44 @@ template class EclOutputBlackOilModule; * \tparam CM The type of the cartesian index mapper. * \param cartMapper The mapper onto cartesian indices. * \param cartesianToActive The mapping of cartesian indices to active indices. - * \param smallGlobalIndex The cartesian cell index of the cell with smaller - * index - * \param largeGlobalIndex The cartesian cell index of the cell with larger - * index - * \return True if the cells have the same i and j indices and all cartesian - * cells + * \param smallGlobalIndex The cartesian cell index of the cell with smaller index + * \param largeGlobalIndex The cartesian cell index of the cell with larger index + * \return True if the cells have the same i and j indices and all cartesian cells * between them are inactive. */ -inline bool -directVerticalNeighbors(const std::array &cartDims, - const std::unordered_map &cartesianToActive, - int smallGlobalIndex, int largeGlobalIndex) { - assert(smallGlobalIndex <= largeGlobalIndex); - std::array ijk1, ijk2; - auto globalToIjk = [cartDims](int gc) { - std::array ijk; - ijk[0] = gc % cartDims[0]; - gc /= cartDims[0]; - ijk[1] = gc % cartDims[1]; - ijk[2] = gc / cartDims[1]; - return ijk; - }; - ijk1 = globalToIjk(smallGlobalIndex); - ijk2 = globalToIjk(largeGlobalIndex); - assert(ijk2[2] >= ijk1[2]); +inline +bool directVerticalNeighbors(const std::array& cartDims, + const std::unordered_map& cartesianToActive, + int smallGlobalIndex, int largeGlobalIndex) +{ + assert(smallGlobalIndex <= largeGlobalIndex); + std::array ijk1, ijk2; + auto globalToIjk = [cartDims](int gc) { + std::array ijk; + ijk[0] = gc % cartDims[0]; + gc /= cartDims[0]; + ijk[1] = gc % cartDims[1]; + ijk[2] = gc / cartDims[1]; + return ijk; + }; + ijk1 = globalToIjk(smallGlobalIndex); + ijk2 = globalToIjk(largeGlobalIndex); + assert(ijk2[2]>=ijk1[2]); - if (ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) { - assert((largeGlobalIndex - smallGlobalIndex) % - (cartDims[0] * cartDims[1]) == - 0); - for (int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; - gi < largeGlobalIndex; gi += cartDims[0] * cartDims[1]) { - if (cartesianToActive.find(gi) != cartesianToActive.end()) { + if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1) + { + assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0); + for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex; + gi += cartDims[0] * cartDims[1] ) + { + if ( cartesianToActive.find( gi ) != cartesianToActive.end() ) + { + return false; + } + } + return true; + } else return false; - } - } - return true; - } else - return false; } /*! @@ -140,614 +141,600 @@ directVerticalNeighbors(const std::array &cartDims, * - This class requires to use the black oil model with the element * centered finite volume discretization. */ -template class EclWriter { - typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; - typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::template Codim<0>::Iterator ElementIterator; +template +class EclWriter +{ + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GridView::template Codim<0>::Iterator ElementIterator; - typedef CollectDataToIORank CollectDataToIORankType; + typedef CollectDataToIORank CollectDataToIORankType; - typedef std::vector ScalarBuffer; + typedef std::vector ScalarBuffer; + + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; - enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; - enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; public: - static void registerParameters() { - EclOutputBlackOilModule::registerParameters(); - - EWOMS_REGISTER_PARAM(TypeTag, bool, EnableAsyncEclOutput, - "Write the ECL-formated results in a non-blocking way " - "(i.e., using a separate thread)."); - } - - // The Simulator object should preferably have been const - the - // only reason that is not the case is due to the SummaryState - // object owned deep down by the vanguard. - EclWriter(Simulator &simulator) - : simulator_(simulator), collectToIORank_(simulator_.vanguard()), - eclOutputModule_(simulator, collectToIORank_) { - if (collectToIORank_.isIORank()) { - globalGrid_ = simulator_.vanguard().grid(); - globalGrid_.switchToGlobalView(); - eclIO_.reset(new Opm::EclipseIO( - simulator_.vanguard().eclState(), - Opm::UgGridHelpers::createEclipseGrid( - globalGrid_, simulator_.vanguard().eclState().getInputGrid()), - simulator_.vanguard().schedule(), - simulator_.vanguard().summaryConfig())); - } - - // create output thread if enabled and rank is I/O rank - // async output is enabled by default if pthread are enabled - bool enableAsyncOutput = - EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput); - int numWorkerThreads = 0; - if (enableAsyncOutput && collectToIORank_.isIORank()) - numWorkerThreads = 1; - taskletRunner_.reset(new TaskletRunner(numWorkerThreads)); - } - - ~EclWriter() {} - - const Opm::EclipseIO &eclIO() const { - assert(eclIO_); - return *eclIO_; - } - - void writeInit() { - if (collectToIORank_.isIORank()) { - std::map> integerVectors; - if (collectToIORank_.isParallel()) - integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); - auto cartMap = Opm::cartesianToCompressed( - globalGrid_.size(0), Opm::UgGridHelpers::globalCell(globalGrid_)); - eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, - exportNncStructure_(cartMap)); - } - } - - /*! - * \brief collect and pass data and pass it to eclIO writer - */ - - void evalSummaryState(bool isSubStep) { - int reportStepNum = simulator_.episodeIndex() + 1; - /* - The summary data is not evaluated for timestep 0, that is - implemented with a: - - if (time_step == 0) - return; - - check somewhere in the summary code. When the summary code was - split in separate methods Summary::eval() and - Summary::add_timestep() it was necessary to pull this test out - here to ensure that the well and group related keywords in the - restart file, like XWEL and XGRP were "correct" also in the - initial report step. - - "Correct" in this context means unchanged behavior, might very - well be more correct to actually remove this if test. - */ - if (reportStepNum == 0) - return; - - Scalar curTime = simulator_.time() + simulator_.timeStepSize(); - Scalar totalCpuTime = simulator_.executionTimer().realTimeElapsed() + - simulator_.setupTimer().realTimeElapsed() + - simulator_.vanguard().externalSetupTime(); - - Opm::data::Wells localWellData = - simulator_.problem().wellModel().wellData(); - - const auto &gridView = simulator_.vanguard().gridView(); - int numElements = gridView.size(/*codim=*/0); - bool log = collectToIORank_.isIORank(); - eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, - /*isRestart*/ false); - - ElementContext elemCtx(simulator_); - ElementIterator elemIt = gridView.template begin(); - const ElementIterator &elemEndIt = gridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const Element &elem = *elemIt; - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - eclOutputModule_.processElement(elemCtx); - } - - if (collectToIORank_.isParallel()) - collectToIORank_.collect({}, eclOutputModule_.getBlockData(), - localWellData); - - std::map miscSummaryData; - std::map> regionData; - eclOutputModule_.outputFipLog(miscSummaryData, regionData, isSubStep); - - bool forceDisableProdOutput = false; - bool forceDisableInjOutput = false; - bool forceDisableCumOutput = false; - eclOutputModule_.outputProdLog(reportStepNum, isSubStep, forceDisableProdOutput); - eclOutputModule_.outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput); - eclOutputModule_.outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput); - - std::vector buffer; - if (collectToIORank_.isIORank()) { - const auto &summary = eclIO_->summary(); - const auto &eclState = simulator_.vanguard().eclState(); - - // Add TCPU - if (totalCpuTime != 0.0) - miscSummaryData["TCPU"] = totalCpuTime; - - const Opm::data::Wells &wellData = collectToIORank_.isParallel() - ? collectToIORank_.globalWellData() - : localWellData; - - const std::map, double> &blockData = - collectToIORank_.isParallel() ? collectToIORank_.globalBlockData() - : eclOutputModule_.getBlockData(); - - summary.eval(summaryState(), reportStepNum, curTime, eclState, schedule(), - wellData, miscSummaryData, regionData, blockData); - buffer = summaryState().serialize(); - } - - if (collectToIORank_.isParallel()) { -#ifdef HAVE_MPI - unsigned long buffer_size = buffer.size(); - MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, - MPI_COMM_WORLD); - if (!collectToIORank_.isIORank()) - buffer.resize(buffer_size); - - MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, - MPI_COMM_WORLD); - if (!collectToIORank_.isIORank()) { - Opm::SummaryState &st = summaryState(); - st.deserialize(buffer); - } -#endif - } - } - - void writeOutput(bool isSubStep) { - Scalar curTime = simulator_.time() + simulator_.timeStepSize(); - Scalar nextStepSize = simulator_.problem().nextTimeStepSize(); - - // output using eclWriter if enabled - Opm::data::Wells localWellData = - simulator_.problem().wellModel().wellData(); - - int reportStepNum = simulator_.episodeIndex() + 1; - const auto &gridView = simulator_.vanguard().gridView(); - int numElements = gridView.size(/*codim=*/0); - bool log = collectToIORank_.isIORank(); - eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, - /*isRestart*/ false); - - ElementContext elemCtx(simulator_); - ElementIterator elemIt = gridView.template begin(); - const ElementIterator &elemEndIt = gridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const Element &elem = *elemIt; - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - eclOutputModule_.processElement(elemCtx); - } - eclOutputModule_.outputErrorLog(); - - // collect all data to I/O rank and assign to sol - Opm::data::Solution localCellData = {}; - if (!isSubStep) - eclOutputModule_.assignToSolution(localCellData); - - // add cell data to perforations for Rft output - if (!isSubStep) - eclOutputModule_.addRftDataToWells(localWellData, reportStepNum); - - if (collectToIORank_.isParallel()) - collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(), - localWellData); - - if (collectToIORank_.isIORank()) { - const auto &eclState = simulator_.vanguard().eclState(); - const auto &simConfig = eclState.getSimulationConfig(); - - bool enableDoublePrecisionOutput = - EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision); - const Opm::data::Solution &cellData = - collectToIORank_.isParallel() ? collectToIORank_.globalCellData() - : localCellData; - const Opm::data::Wells &wellData = collectToIORank_.isParallel() - ? collectToIORank_.globalWellData() - : localWellData; - Opm::RestartValue restartValue(cellData, wellData); - - if (simConfig.useThresholdPressure()) - restartValue.addExtra("THRESHPR", Opm::UnitSystem::measure::pressure, - simulator_.problem().thresholdPressure().data()); - - // Add suggested next timestep to extra data. - if (!isSubStep) - restartValue.addExtra("OPMEXTRA", std::vector(1, nextStepSize)); - - // first, create a tasklet to write the data for the current time step to - // disk - auto eclWriteTasklet = std::make_shared( - summaryState(), *eclIO_, reportStepNum, isSubStep, curTime, - restartValue, enableDoublePrecisionOutput); - - // then, make sure that the previous I/O request has been completed and - // the - // number of incomplete tasklets does not increase between time steps - taskletRunner_->barrier(); - - // finally, start a new output writing job - taskletRunner_->dispatch(eclWriteTasklet); - } - } - - void beginRestart() { - bool enableHysteresis = - simulator_.problem().materialLawManager()->enableHysteresis(); - bool enableSwatinit = simulator_.vanguard() - .eclState() - .get3DProperties() - .hasDeckDoubleGridProperty("SWATINIT"); - std::vector solutionKeys{ - {"PRESSURE", Opm::UnitSystem::measure::pressure}, - {"SWAT", Opm::UnitSystem::measure::identity, - static_cast( - FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))}, - {"SGAS", Opm::UnitSystem::measure::identity, - static_cast( - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))}, - {"TEMP", Opm::UnitSystem::measure::temperature, enableEnergy}, - {"SSOLVENT", Opm::UnitSystem::measure::identity, enableSolvent}, - {"RS", Opm::UnitSystem::measure::gas_oil_ratio, - FluidSystem::enableDissolvedGas()}, - {"RV", Opm::UnitSystem::measure::oil_gas_ratio, - FluidSystem::enableVaporizedOil()}, - {"SOMAX", Opm::UnitSystem::measure::identity, - simulator_.problem().vapparsActive()}, - {"PCSWM_OW", Opm::UnitSystem::measure::identity, enableHysteresis}, - {"KRNSW_OW", Opm::UnitSystem::measure::identity, enableHysteresis}, - {"PCSWM_GO", Opm::UnitSystem::measure::identity, enableHysteresis}, - {"KRNSW_GO", Opm::UnitSystem::measure::identity, enableHysteresis}, - {"PPCW", Opm::UnitSystem::measure::pressure, enableSwatinit}}; - - const auto &inputThpres = - eclState().getSimulationConfig().getThresholdPressure(); - std::vector extraKeys = { - {"OPMEXTRA", Opm::UnitSystem::measure::identity, false}, - {"THRESHPR", Opm::UnitSystem::measure::pressure, inputThpres.active()}}; - - // The episodeIndex is rewined one back before beginRestart is called - // and can not be used here. - // We just ask the initconfig directly to be sure that we use the correct - // index. - const auto &initconfig = simulator_.vanguard().eclState().getInitConfig(); - int restartStepIdx = initconfig.getRestartStep(); - - const auto &gridView = simulator_.vanguard().gridView(); - unsigned numElements = gridView.size(/*codim=*/0); - eclOutputModule_.allocBuffers(numElements, restartStepIdx, - /*isSubStep=*/false, /*log=*/false, - /*isRestart*/ true); - + static void registerParameters() { - Opm::SummaryState &summaryState = simulator_.vanguard().summaryState(); - auto restartValues = - loadParallelRestart(eclIO_.get(), summaryState, solutionKeys, - extraKeys, gridView.grid().comm()); + EclOutputBlackOilModule::registerParameters(); - for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { - unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx); - eclOutputModule_.setRestart(restartValues.solution, elemIdx, globalIdx); - } - - if (inputThpres.active()) { - Simulator &mutableSimulator = const_cast(simulator_); - auto &thpres = mutableSimulator.problem().thresholdPressure(); - const auto &thpresValues = restartValues.getExtra("THRESHPR"); - thpres.setFromRestart(thpresValues); - } - restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0]; - - // initialize the well model from restart values - simulator_.problem().wellModel().initFromRestartFile(restartValues); + EWOMS_REGISTER_PARAM(TypeTag, bool, EnableAsyncEclOutput, + "Write the ECL-formated results in a non-blocking way (i.e., using a separate thread)."); } - } - void endRestart() {} + // The Simulator object should preferably have been const - the + // only reason that is not the case is due to the SummaryState + // object owned deep down by the vanguard. + EclWriter(Simulator& simulator) + : simulator_(simulator) + , collectToIORank_(simulator_.vanguard()) + , eclOutputModule_(simulator, collectToIORank_) + { + if (collectToIORank_.isIORank()) { + globalGrid_ = simulator_.vanguard().grid(); + globalGrid_.switchToGlobalView(); + eclIO_.reset(new Opm::EclipseIO(simulator_.vanguard().eclState(), + Opm::UgGridHelpers::createEclipseGrid(globalGrid_, simulator_.vanguard().eclState().getInputGrid()), + simulator_.vanguard().schedule(), + simulator_.vanguard().summaryConfig())); + } - const EclOutputBlackOilModule &eclOutputModule() const { - return eclOutputModule_; - } + // create output thread if enabled and rank is I/O rank + // async output is enabled by default if pthread are enabled + bool enableAsyncOutput = EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncEclOutput); + int numWorkerThreads = 0; + if (enableAsyncOutput && collectToIORank_.isIORank()) + numWorkerThreads = 1; + taskletRunner_.reset(new TaskletRunner(numWorkerThreads)); + } + + ~EclWriter() + { } + + const Opm::EclipseIO& eclIO() const + { + assert(eclIO_); + return *eclIO_; + } + + void writeInit() + { + if (collectToIORank_.isIORank()) { + std::map > integerVectors; + if (collectToIORank_.isParallel()) + integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); + auto cartMap = Opm::cartesianToCompressed(globalGrid_.size(0), + Opm::UgGridHelpers::globalCell(globalGrid_)); + eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap)); + } + } + + /*! + * \brief collect and pass data and pass it to eclIO writer + */ + + void evalSummaryState(bool isSubStep) + { + int reportStepNum = simulator_.episodeIndex() + 1; + /* + The summary data is not evaluated for timestep 0, that is + implemented with a: + + if (time_step == 0) + return; + + check somewhere in the summary code. When the summary code was + split in separate methods Summary::eval() and + Summary::add_timestep() it was necessary to pull this test out + here to ensure that the well and group related keywords in the + restart file, like XWEL and XGRP were "correct" also in the + initial report step. + + "Correct" in this context means unchanged behavior, might very + well be more correct to actually remove this if test. + */ + if (reportStepNum == 0) + return; + + Scalar curTime = simulator_.time() + simulator_.timeStepSize(); + Scalar totalCpuTime = + simulator_.executionTimer().realTimeElapsed() + + simulator_.setupTimer().realTimeElapsed() + + simulator_.vanguard().externalSetupTime(); + + Opm::data::Wells localWellData = simulator_.problem().wellModel().wellData(); + + const auto& gridView = simulator_.vanguard().gridView(); + int numElements = gridView.size(/*codim=*/0); + bool log = collectToIORank_.isIORank(); + eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, /*isRestart*/ false); + + ElementContext elemCtx(simulator_); + ElementIterator elemIt = gridView.template begin(); + const ElementIterator& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + eclOutputModule_.processElement(elemCtx); + } + + if (collectToIORank_.isParallel()) + collectToIORank_.collect({}, eclOutputModule_.getBlockData(), localWellData); + + std::map miscSummaryData; + std::map> regionData; + eclOutputModule_.outputFipLog(miscSummaryData, regionData, isSubStep); + + bool forceDisableProdOutput = false; + bool forceDisableInjOutput = false; + bool forceDisableCumOutput = false; + eclOutputModule_.outputProdLog(reportStepNum, isSubStep, forceDisableProdOutput); + eclOutputModule_.outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput); + eclOutputModule_.outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput); + + std::vector buffer; + if (collectToIORank_.isIORank()) { + const auto& summary = eclIO_->summary(); + const auto& eclState = simulator_.vanguard().eclState(); + + // Add TCPU + if (totalCpuTime != 0.0) + miscSummaryData["TCPU"] = totalCpuTime; + + const Opm::data::Wells& wellData = collectToIORank_.isParallel() ? collectToIORank_.globalWellData() : localWellData; + + const std::map, double>& blockData + = collectToIORank_.isParallel() + ? collectToIORank_.globalBlockData() + : eclOutputModule_.getBlockData(); + + summary.eval(summaryState(), + reportStepNum, + curTime, + eclState, + schedule(), + wellData, + miscSummaryData, + regionData, + blockData); + buffer = summaryState().serialize(); + } + + if (collectToIORank_.isParallel()) { +#ifdef HAVE_MPI + unsigned long buffer_size = buffer.size(); + MPI_Bcast(&buffer_size, 1, MPI_UNSIGNED_LONG, collectToIORank_.ioRank, MPI_COMM_WORLD); + if (!collectToIORank_.isIORank()) + buffer.resize( buffer_size ); + + MPI_Bcast(buffer.data(), buffer_size, MPI_CHAR, collectToIORank_.ioRank, MPI_COMM_WORLD); + if (!collectToIORank_.isIORank()) { + Opm::SummaryState& st = summaryState(); + st.deserialize(buffer); + } +#endif + } + } + + + void writeOutput(bool isSubStep) + { + Scalar curTime = simulator_.time() + simulator_.timeStepSize(); + Scalar nextStepSize = simulator_.problem().nextTimeStepSize(); + + // output using eclWriter if enabled + Opm::data::Wells localWellData = simulator_.problem().wellModel().wellData(); + + int reportStepNum = simulator_.episodeIndex() + 1; + const auto& gridView = simulator_.vanguard().gridView(); + int numElements = gridView.size(/*codim=*/0); + bool log = collectToIORank_.isIORank(); + eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log, /*isRestart*/ false); + + ElementContext elemCtx(simulator_); + ElementIterator elemIt = gridView.template begin(); + const ElementIterator& elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + eclOutputModule_.processElement(elemCtx); + } + eclOutputModule_.outputErrorLog(); + + // collect all data to I/O rank and assign to sol + Opm::data::Solution localCellData = {}; + if (!isSubStep) + eclOutputModule_.assignToSolution(localCellData); + + // add cell data to perforations for Rft output + if (!isSubStep) + eclOutputModule_.addRftDataToWells(localWellData, reportStepNum); + + if (collectToIORank_.isParallel()) + collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(), localWellData); + + + if (collectToIORank_.isIORank()) { + const auto& eclState = simulator_.vanguard().eclState(); + const auto& simConfig = eclState.getSimulationConfig(); + + bool enableDoublePrecisionOutput = EWOMS_GET_PARAM(TypeTag, bool, EclOutputDoublePrecision); + const Opm::data::Solution& cellData = collectToIORank_.isParallel() ? collectToIORank_.globalCellData() : localCellData; + const Opm::data::Wells& wellData = collectToIORank_.isParallel() ? collectToIORank_.globalWellData() : localWellData; + Opm::RestartValue restartValue(cellData, wellData); + + if (simConfig.useThresholdPressure()) + restartValue.addExtra("THRESHPR", Opm::UnitSystem::measure::pressure, simulator_.problem().thresholdPressure().data()); + + // Add suggested next timestep to extra data. + if (!isSubStep) + restartValue.addExtra("OPMEXTRA", std::vector(1, nextStepSize)); + + // first, create a tasklet to write the data for the current time step to disk + auto eclWriteTasklet = std::make_shared(summaryState(), + *eclIO_, + reportStepNum, + isSubStep, + curTime, + restartValue, + enableDoublePrecisionOutput); + + // then, make sure that the previous I/O request has been completed and the + // number of incomplete tasklets does not increase between time steps + taskletRunner_->barrier(); + + // finally, start a new output writing job + taskletRunner_->dispatch(eclWriteTasklet); + } + } + + void beginRestart() + { + bool enableHysteresis = simulator_.problem().materialLawManager()->enableHysteresis(); + bool enableSwatinit = simulator_.vanguard().eclState().get3DProperties().hasDeckDoubleGridProperty("SWATINIT"); + std::vector solutionKeys{ + {"PRESSURE", Opm::UnitSystem::measure::pressure}, + {"SWAT", Opm::UnitSystem::measure::identity, static_cast(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))}, + {"SGAS", Opm::UnitSystem::measure::identity, static_cast(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))}, + {"TEMP" , Opm::UnitSystem::measure::temperature, enableEnergy}, + {"SSOLVENT" , Opm::UnitSystem::measure::identity, enableSolvent}, + {"RS", Opm::UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()}, + {"RV", Opm::UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()}, + {"SOMAX", Opm::UnitSystem::measure::identity, simulator_.problem().vapparsActive()}, + {"PCSWM_OW", Opm::UnitSystem::measure::identity, enableHysteresis}, + {"KRNSW_OW", Opm::UnitSystem::measure::identity, enableHysteresis}, + {"PCSWM_GO", Opm::UnitSystem::measure::identity, enableHysteresis}, + {"KRNSW_GO", Opm::UnitSystem::measure::identity, enableHysteresis}, + {"PPCW", Opm::UnitSystem::measure::pressure, enableSwatinit} + }; + + const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure(); + std::vector extraKeys = {{"OPMEXTRA", Opm::UnitSystem::measure::identity, false}, + {"THRESHPR", Opm::UnitSystem::measure::pressure, inputThpres.active()}}; + + // The episodeIndex is rewined one back before beginRestart is called + // and can not be used here. + // We just ask the initconfig directly to be sure that we use the correct + // index. + const auto& initconfig = simulator_.vanguard().eclState().getInitConfig(); + int restartStepIdx = initconfig.getRestartStep(); + + const auto& gridView = simulator_.vanguard().gridView(); + unsigned numElements = gridView.size(/*codim=*/0); + eclOutputModule_.allocBuffers(numElements, restartStepIdx, /*isSubStep=*/false, /*log=*/false, /*isRestart*/ true); + + { + Opm::SummaryState& summaryState = simulator_.vanguard().summaryState(); + auto restartValues = loadParallelRestart(eclIO_.get(), summaryState, solutionKeys, extraKeys, + gridView.grid().comm()); + + for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) { + unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx); + eclOutputModule_.setRestart(restartValues.solution, elemIdx, globalIdx); + } + + if (inputThpres.active()) { + Simulator& mutableSimulator = const_cast(simulator_); + auto& thpres = mutableSimulator.problem().thresholdPressure(); + const auto& thpresValues = restartValues.getExtra("THRESHPR"); + thpres.setFromRestart(thpresValues); + } + restartTimeStepSize_ = restartValues.getExtra("OPMEXTRA")[0]; + + // initialize the well model from restart values + simulator_.problem().wellModel().initFromRestartFile(restartValues); + } + } + + void endRestart() + {} + + const EclOutputBlackOilModule& eclOutputModule() const + { return eclOutputModule_; } + + Scalar restartTimeStepSize() const + { return restartTimeStepSize_; } - Scalar restartTimeStepSize() const { return restartTimeStepSize_; } private: - static bool enableEclOutput_() { - return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); - } + static bool enableEclOutput_() + { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); } - Opm::data::Solution - computeTrans_(const std::unordered_map &cartesianToActive) const { - const auto &cartMapper = simulator_.vanguard().cartesianIndexMapper(); - const auto &cartDims = cartMapper.cartesianDimensions(); - const int globalSize = cartDims[0] * cartDims[1] * cartDims[2]; + Opm::data::Solution computeTrans_(const std::unordered_map& cartesianToActive) const + { + const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper(); + const auto& cartDims = cartMapper.cartesianDimensions(); + const int globalSize = cartDims[0]*cartDims[1]*cartDims[2]; - Opm::data::CellData tranx = {Opm::UnitSystem::measure::transmissibility, - std::vector(globalSize), - Opm::data::TargetType::INIT}; - Opm::data::CellData trany = {Opm::UnitSystem::measure::transmissibility, - std::vector(globalSize), - Opm::data::TargetType::INIT}; - Opm::data::CellData tranz = {Opm::UnitSystem::measure::transmissibility, - std::vector(globalSize), - Opm::data::TargetType::INIT}; + Opm::data::CellData tranx = {Opm::UnitSystem::measure::transmissibility, std::vector(globalSize), Opm::data::TargetType::INIT}; + Opm::data::CellData trany = {Opm::UnitSystem::measure::transmissibility, std::vector(globalSize), Opm::data::TargetType::INIT}; + Opm::data::CellData tranz = {Opm::UnitSystem::measure::transmissibility, std::vector(globalSize), Opm::data::TargetType::INIT}; - for (size_t i = 0; i < tranx.data.size(); ++i) { - tranx.data[0] = 0.0; - trany.data[0] = 0.0; - tranz.data[0] = 0.0; - } + for (size_t i = 0; i < tranx.data.size(); ++i) { + tranx.data[0] = 0.0; + trany.data[0] = 0.0; + tranz.data[0] = 0.0; + } - typedef typename Grid::LeafGridView GlobalGridView; - const GlobalGridView &globalGridView = globalGrid_.leafGridView(); -#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) - typedef Dune::MultipleCodimMultipleGeomTypeMapper - ElementMapper; - ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + typedef typename Grid :: LeafGridView GlobalGridView; + const GlobalGridView& globalGridView = globalGrid_.leafGridView(); +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); #else - typedef Dune::MultipleCodimMultipleGeomTypeMapper - ElementMapper; - ElementMapper globalElemMapper(globalGridView); + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + ElementMapper globalElemMapper(globalGridView); #endif - const auto &cartesianCellIdx = globalGrid_.globalCell(); - const EclTransmissibility *globalTrans; + const auto& cartesianCellIdx = globalGrid_.globalCell(); + const EclTransmissibility* globalTrans; - if (!collectToIORank_.isParallel()) { - // in the sequential case we must use the transmissibilites defined by - // the problem. (because in the sequential case, the grid manager does - // not compute "global" transmissibilities for performance reasons. in - // the parallel case, the problem's transmissibilities can't be used - // because this object refers to the distributed grid and we need the - // sequential version here.) - globalTrans = &simulator_.problem().eclTransmissibilities(); - } else { - globalTrans = &(simulator_.vanguard().globalTransmissibility()); - } - - auto elemIt = globalGridView.template begin(); - const auto &elemEndIt = globalGridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const auto &elem = *elemIt; - - auto isIt = globalGridView.ibegin(elem); - const auto &isEndIt = globalGridView.iend(elem); - for (; isIt != isEndIt; ++isIt) { - const auto &is = *isIt; - - if (!is.neighbor()) - continue; // intersection is on the domain boundary - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - continue; // we only need to handle each connection once, thank you. - - // Ordering of compressed and uncompressed index should be the same - assert(cartesianCellIdx[c1] <= cartesianCellIdx[c2]); - int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]); - int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]); - - if (gc2 - gc1 == 1) { - tranx.data[gc1] = globalTrans->transmissibility(c1, c2); - continue; // skip other if clauses as they are false, last one needs - // some computation + if (!collectToIORank_.isParallel()) + { + // in the sequential case we must use the transmissibilites defined by + // the problem. (because in the sequential case, the grid manager does + // not compute "global" transmissibilities for performance reasons. in + // the parallel case, the problem's transmissibilities can't be used + // because this object refers to the distributed grid and we need the + // sequential version here.) + globalTrans = &simulator_.problem().eclTransmissibilities(); + } + else + { + globalTrans = &(simulator_.vanguard().globalTransmissibility()); } - if (gc2 - gc1 == cartDims[0]) { - trany.data[gc1] = globalTrans->transmissibility(c1, c2); - continue; // skipt next if clause as it needs some computation + auto elemIt = globalGridView.template begin(); + const auto& elemEndIt = globalGridView.template end(); + for (; elemIt != elemEndIt; ++ elemIt) { + const auto& elem = *elemIt; + + auto isIt = globalGridView.ibegin(elem); + const auto& isEndIt = globalGridView.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + const auto& is = *isIt; + + if (!is.neighbor()) + continue; // intersection is on the domain boundary + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + continue; // we only need to handle each connection once, thank you. + + // Ordering of compressed and uncompressed index should be the same + assert(cartesianCellIdx[c1] <= cartesianCellIdx[c2]); + int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]); + int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]); + + if (gc2 - gc1 == 1) { + tranx.data[gc1] = globalTrans->transmissibility(c1, c2); + continue; // skip other if clauses as they are false, last one needs some computation + } + + if (gc2 - gc1 == cartDims[0]) { + trany.data[gc1] = globalTrans->transmissibility(c1, c2); + continue; // skipt next if clause as it needs some computation + } + + if ( gc2 - gc1 == cartDims[0]*cartDims[1] || + directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) + tranz.data[gc1] = globalTrans->transmissibility(c1, c2); + } } - if (gc2 - gc1 == cartDims[0] * cartDims[1] || - directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2)) - tranz.data[gc1] = globalTrans->transmissibility(c1, c2); - } + return {{"TRANX", tranx}, + {"TRANY", trany}, + {"TRANZ", tranz}}; } - return {{"TRANX", tranx}, {"TRANY", trany}, {"TRANZ", tranz}}; - } + Opm::NNC exportNncStructure_(const std::unordered_map& cartesianToActive) const + { + std::size_t nx = eclState().getInputGrid().getNX(); + std::size_t ny = eclState().getInputGrid().getNY(); + auto nncData = sortNncAndApplyEditnnc(eclState().getInputNNC().nncdata(), + eclState().getInputEDITNNC().data()); + const auto& unitSystem = simulator_.vanguard().deck().getActiveUnitSystem(); + std::vector outputNnc; + std::size_t index = 0; - Opm::NNC exportNncStructure_( - const std::unordered_map &cartesianToActive) const { - std::size_t nx = eclState().getInputGrid().getNX(); - std::size_t ny = eclState().getInputGrid().getNY(); - auto nncData = sortNncAndApplyEditnnc(eclState().getInputNNC().nncdata(), - eclState().getInputEDITNNC().data()); - const auto &unitSystem = simulator_.vanguard().deck().getActiveUnitSystem(); - std::vector outputNnc; - std::size_t index = 0; + for( const auto& entry : nncData ) { + // test whether NNC is not a neighboring connection + // cell2>=cell1 holds due to sortNncAndApplyEditnnc + assert( entry.cell2 >= entry.cell1 ); + auto cellDiff = entry.cell2 - entry.cell1; - for (const auto &entry : nncData) { - // test whether NNC is not a neighboring connection - // cell2>=cell1 holds due to sortNncAndApplyEditnnc - assert(entry.cell2 >= entry.cell1); - auto cellDiff = entry.cell2 - entry.cell1; + if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) { + auto tt = unitSystem.from_si(Opm::UnitSystem::measure::transmissibility, entry.trans); + // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6 + if ( tt >= 1.0e-6 ) + outputNnc.emplace_back(entry.cell1, entry.cell2, entry.trans); + } + ++index; + } - if (cellDiff != 1 && cellDiff != nx && cellDiff != nx * ny) { - auto tt = unitSystem.from_si(Opm::UnitSystem::measure::transmissibility, - entry.trans); - // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems - // like the threshold is 1.0e-6 - if (tt >= 1.0e-6) - outputNnc.emplace_back(entry.cell1, entry.cell2, entry.trans); - } - ++index; + auto nncCompare = []( const Opm::NNCdata& nnc1, const Opm::NNCdata& nnc2){ + return nnc1.cell1 < nnc2.cell1 || + ( nnc1.cell1 == nnc2.cell1 && nnc1.cell2 < nnc2.cell2);}; + // Sort the nncData values from the deck as they need to be + // Checked when writing NNC transmissibilities from the simulation. + std::sort(nncData.begin(), nncData.end(), nncCompare); + + typedef typename Grid :: LeafGridView GlobalGridView; + const GlobalGridView& globalGridView = globalGrid_.leafGridView(); +#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + +#else + typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; + ElementMapper globalElemMapper(globalGridView); +#endif + + const EclTransmissibility* globalTrans; + if (!collectToIORank_.isParallel()) { + // in the sequential case we must use the transmissibilites defined by + // the problem. (because in the sequential case, the grid manager does + // not compute "global" transmissibilities for performance reasons. in + // the parallel case, the problem's transmissibilities can't be used + // because this object refers to the distributed grid and we need the + // sequential version here.) + globalTrans = &simulator_.problem().eclTransmissibilities(); + } + else + { + globalTrans = &(simulator_.vanguard().globalTransmissibility()); + } + + auto cartDims = simulator_.vanguard().cartesianIndexMapper().cartesianDimensions(); + auto elemIt = globalGridView.template begin(); + const auto& elemEndIt = globalGridView.template end(); + for (; elemIt != elemEndIt; ++ elemIt) { + const auto& elem = *elemIt; + + auto isIt = globalGridView.ibegin(elem); + const auto& isEndIt = globalGridView.iend(elem); + for (; isIt != isEndIt; ++ isIt) { + const auto& is = *isIt; + + if (!is.neighbor()) + continue; // intersection is on the domain boundary + + unsigned c1 = globalElemMapper.index(is.inside()); + unsigned c2 = globalElemMapper.index(is.outside()); + + if (c1 > c2) + continue; // we only need to handle each connection once, thank you. + + // TODO (?): use the cartesian index mapper to make this code work + // with grids other than Dune::CpGrid. The problem is that we need + // the a mapper for the sequential grid, not for the distributed one. + std::size_t cc1 = globalGrid_.globalCell()[c1]; + std::size_t cc2 = globalGrid_.globalCell()[c2]; + + if ( cc2 < cc1 ) + std::swap(cc1, cc2); + + auto cellDiff = cc2 - cc1; + + if (cellDiff != 1 && + cellDiff != nx && + cellDiff != nx*ny && + ! directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { + // We need to check whether an NNC for this face was also specified + // via the NNC keyword in the deck (i.e. in the first origNncSize entries. + auto t = globalTrans->transmissibility(c1, c2); + auto candidate = std::lower_bound(nncData.begin(), nncData.end(), Opm::NNCdata(cc1, cc2, 0.0), nncCompare); + + while ( candidate != nncData.end() && candidate->cell1 == cc1 + && candidate->cell2 == cc2) { + t -= candidate->trans; + ++candidate; + } + // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC + // with corresponding EDITNNC above). In addition we do set small transmissibilties + // to zero when setting up the simulator. These will be ignored here, too. + auto tt = unitSystem.from_si(Opm::UnitSystem::measure::transmissibility, std::abs(t)); + if ( tt > 1e-12 ) + outputNnc.push_back({cc1, cc2, t}); + } + } + } + Opm::NNC ret; + for(const auto& nncItem: outputNnc) + ret.addNNC(nncItem.cell1, nncItem.cell2, nncItem.trans); + return ret; } - auto nncCompare = [](const Opm::NNCdata &nnc1, const Opm::NNCdata &nnc2) { - return nnc1.cell1 < nnc2.cell1 || - (nnc1.cell1 == nnc2.cell1 && nnc1.cell2 < nnc2.cell2); + struct EclWriteTasklet + : public TaskletInterface + { + Opm::SummaryState summaryState_; + Opm::EclipseIO& eclIO_; + int reportStepNum_; + bool isSubStep_; + double secondsElapsed_; + Opm::RestartValue restartValue_; + bool writeDoublePrecision_; + + explicit EclWriteTasklet(const Opm::SummaryState& summaryState, + Opm::EclipseIO& eclIO, + int reportStepNum, + bool isSubStep, + double secondsElapsed, + Opm::RestartValue restartValue, + bool writeDoublePrecision) + : summaryState_(summaryState) + , eclIO_(eclIO) + , reportStepNum_(reportStepNum) + , isSubStep_(isSubStep) + , secondsElapsed_(secondsElapsed) + , restartValue_(restartValue) + , writeDoublePrecision_(writeDoublePrecision) + { } + + // callback to eclIO serial writeTimeStep method + void run() + { + eclIO_.writeTimeStep(summaryState_, + reportStepNum_, + isSubStep_, + secondsElapsed_, + restartValue_, + writeDoublePrecision_); + } }; - // Sort the nncData values from the deck as they need to be - // Checked when writing NNC transmissibilities from the simulation. - std::sort(nncData.begin(), nncData.end(), nncCompare); - typedef typename Grid::LeafGridView GlobalGridView; - const GlobalGridView &globalGridView = globalGrid_.leafGridView(); -#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 6) - typedef Dune::MultipleCodimMultipleGeomTypeMapper - ElementMapper; - ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); + const Opm::EclipseState& eclState() const + { return simulator_.vanguard().eclState(); } -#else - typedef Dune::MultipleCodimMultipleGeomTypeMapper - ElementMapper; - ElementMapper globalElemMapper(globalGridView); -#endif + Opm::SummaryState& summaryState() + { return simulator_.vanguard().summaryState(); } - const EclTransmissibility *globalTrans; - if (!collectToIORank_.isParallel()) { - // in the sequential case we must use the transmissibilites defined by - // the problem. (because in the sequential case, the grid manager does - // not compute "global" transmissibilities for performance reasons. in - // the parallel case, the problem's transmissibilities can't be used - // because this object refers to the distributed grid and we need the - // sequential version here.) - globalTrans = &simulator_.problem().eclTransmissibilities(); - } else { - globalTrans = &(simulator_.vanguard().globalTransmissibility()); - } + const Opm::Schedule& schedule() const + { return simulator_.vanguard().schedule(); } - auto cartDims = - simulator_.vanguard().cartesianIndexMapper().cartesianDimensions(); - auto elemIt = globalGridView.template begin(); - const auto &elemEndIt = globalGridView.template end(); - for (; elemIt != elemEndIt; ++elemIt) { - const auto &elem = *elemIt; + Simulator& simulator_; + CollectDataToIORankType collectToIORank_; + EclOutputBlackOilModule eclOutputModule_; + std::unique_ptr eclIO_; + Grid globalGrid_; + std::unique_ptr taskletRunner_; + Scalar restartTimeStepSize_; - auto isIt = globalGridView.ibegin(elem); - const auto &isEndIt = globalGridView.iend(elem); - for (; isIt != isEndIt; ++isIt) { - const auto &is = *isIt; - if (!is.neighbor()) - continue; // intersection is on the domain boundary - - unsigned c1 = globalElemMapper.index(is.inside()); - unsigned c2 = globalElemMapper.index(is.outside()); - - if (c1 > c2) - continue; // we only need to handle each connection once, thank you. - - // TODO (?): use the cartesian index mapper to make this code work - // with grids other than Dune::CpGrid. The problem is that we need - // the a mapper for the sequential grid, not for the distributed one. - std::size_t cc1 = globalGrid_.globalCell()[c1]; - std::size_t cc2 = globalGrid_.globalCell()[c2]; - - if (cc2 < cc1) - std::swap(cc1, cc2); - - auto cellDiff = cc2 - cc1; - - if (cellDiff != 1 && cellDiff != nx && cellDiff != nx * ny && - !directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) { - // We need to check whether an NNC for this face was also specified - // via the NNC keyword in the deck (i.e. in the first origNncSize - // entries. - auto t = globalTrans->transmissibility(c1, c2); - auto candidate = - std::lower_bound(nncData.begin(), nncData.end(), - Opm::NNCdata(cc1, cc2, 0.0), nncCompare); - - while (candidate != nncData.end() && candidate->cell1 == cc1 && - candidate->cell2 == cc2) { - t -= candidate->trans; - ++candidate; - } - // eclipse ignores NNCs with zero transmissibility (different - // threshold than for NNC - // with corresponding EDITNNC above). In addition we do set small - // transmissibilties - // to zero when setting up the simulator. These will be ignored here, - // too. - auto tt = unitSystem.from_si( - Opm::UnitSystem::measure::transmissibility, std::abs(t)); - if (tt > 1e-12) - outputNnc.push_back({cc1, cc2, t}); - } - } - } - Opm::NNC ret; - for (const auto &nncItem : outputNnc) - ret.addNNC(nncItem.cell1, nncItem.cell2, nncItem.trans); - return ret; - } - - struct EclWriteTasklet : public TaskletInterface { - Opm::SummaryState summaryState_; - Opm::EclipseIO &eclIO_; - int reportStepNum_; - bool isSubStep_; - double secondsElapsed_; - Opm::RestartValue restartValue_; - bool writeDoublePrecision_; - - explicit EclWriteTasklet(const Opm::SummaryState &summaryState, - Opm::EclipseIO &eclIO, int reportStepNum, - bool isSubStep, double secondsElapsed, - Opm::RestartValue restartValue, - bool writeDoublePrecision) - : summaryState_(summaryState), eclIO_(eclIO), - reportStepNum_(reportStepNum), isSubStep_(isSubStep), - secondsElapsed_(secondsElapsed), restartValue_(restartValue), - writeDoublePrecision_(writeDoublePrecision) {} - - // callback to eclIO serial writeTimeStep method - void run() { - eclIO_.writeTimeStep(summaryState_, reportStepNum_, isSubStep_, - secondsElapsed_, restartValue_, - writeDoublePrecision_); - } - }; - - const Opm::EclipseState &eclState() const { - return simulator_.vanguard().eclState(); - } - - Opm::SummaryState &summaryState() { - return simulator_.vanguard().summaryState(); - } - - const Opm::Schedule &schedule() const { - return simulator_.vanguard().schedule(); - } - - Simulator &simulator_; - CollectDataToIORankType collectToIORank_; - EclOutputBlackOilModule eclOutputModule_; - std::unique_ptr eclIO_; - Grid globalGrid_; - std::unique_ptr taskletRunner_; - Scalar restartTimeStepSize_; }; } // namespace Opm