diff --git a/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h b/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h index 8942f7beae..b5292e26eb 100644 --- a/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h +++ b/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h @@ -18,56 +18,83 @@ namespace RigFlowDiagInterfaceTools { + template inline Opm::FlowDiagnostics::ConnectionValues - extractFluxField(const Opm::ECLGraph& G, const bool compute_fluxes) + extractFluxField(const Opm::ECLGraph& G, + FluxCalc&& getFlux) { using ConnVals = Opm::FlowDiagnostics::ConnectionValues; + + const auto actPh = G.activePhases(); + auto flux = ConnVals(ConnVals::NumConnections{ G.numConnections() }, - ConnVals::NumPhases{ 3 }); - + ConnVals::NumPhases{ actPh.size() }); + auto phas = ConnVals::PhaseID{ 0 }; - - Opm::ECLFluxCalc calc(G); - - const auto phases = { Opm::ECLGraph::PhaseIndex::Aqua , - Opm::ECLGraph::PhaseIndex::Liquid , - Opm::ECLGraph::PhaseIndex::Vapour }; - - for ( const auto& p : phases ) - { - const auto pflux = compute_fluxes ? calc.flux(p) : G.flux(p); - if ( ! pflux.empty() ) - { - assert (pflux.size() == flux.numConnections()); + + for (const auto& p : actPh) { + const auto pflux = getFlux(p); + + if (!pflux.empty()) { + assert(pflux.size() == flux.numConnections()); + auto conn = ConnVals::ConnID{ 0 }; - for ( const auto& v : pflux ) - { + for (const auto& v : pflux) { flux(conn, phas) = v; conn.id += 1; } } + phas.id += 1; } - + return flux; } + inline Opm::FlowDiagnostics::ConnectionValues + extractFluxField(const Opm::ECLGraph& G, + const Opm::ECLRestartData& rstrt, + const bool compute_fluxes) + { + if (compute_fluxes) { + Opm::ECLFluxCalc calc(G); + + auto getFlux = [&calc, &rstrt] + (const Opm::ECLGraph::PhaseIndex p) + { + return calc.flux(rstrt, p); + }; + + return extractFluxField(G, getFlux); + } + + auto getFlux = [&G, &rstrt] + (const Opm::ECLGraph::PhaseIndex p) + { + return G.flux(rstrt, p); + }; + + return extractFluxField(G, getFlux); + } + template Opm::FlowDiagnostics::CellSetValues - extractWellFlows(const Opm::ECLGraph& G, - const WellFluxes& well_fluxes) + extractWellFlows(const Opm::ECLGraph& G, + const WellFluxes& well_fluxes) { Opm::FlowDiagnostics::CellSetValues inflow; for (const auto& well : well_fluxes) { for (const auto& completion : well.completions) { - const int grid_index = completion.grid_index; + const auto& gridName = completion.gridName; const auto& ijk = completion.ijk; - const int cell_index = G.activeCell(ijk, grid_index); + const int cell_index = G.activeCell(ijk, gridName); if (cell_index >= 0) { - auto iterIsInsertedPair = inflow.emplace(cell_index, completion.reservoir_inflow_rate); - if (!iterIsInsertedPair.second){ // Not inserted, We had something there already - iterIsInsertedPair.first->second += completion.reservoir_inflow_rate; // Accumulate the flow at this connection - } + // Since inflow is a std::map, if the key was not + // already present operator[] will insert a + // value-initialized value (as in T() for a type + // T), which is zero for built-in numerical types, + // including double. + inflow[cell_index] += completion.reservoir_inflow_rate; } } } @@ -75,7 +102,6 @@ namespace RigFlowDiagInterfaceTools { return inflow; } - } diff --git a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp index 9fd5b6dfa1..b9e34ea5c5 100644 --- a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp +++ b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp @@ -104,19 +104,25 @@ void RigFlowDiagTimeStepResult::addResult(const RigFlowDiagResultAddress& resAdd } -class RigOpmFldStaticData : public cvf::Object +class RigOpmFlowDiagStaticData : public cvf::Object { public: - RigOpmFldStaticData(const std::string& grid, const std::string& init) : m_eclGraph(Opm::ECLGraph::load(grid, init)), m_hasUnifiedRestartFile(false) + RigOpmFlowDiagStaticData(const std::string& grid, const std::string& init) { - m_poreVolume = m_eclGraph.poreVolume(); + Opm::ECLInitFileData initData(init); + + m_eclGraph.reset(new Opm::ECLGraph(Opm::ECLGraph::load(grid, initData))); + + m_hasUnifiedRestartFile = false; + m_poreVolume = m_eclGraph->poreVolume(); } - Opm::ECLGraph m_eclGraph; + std::unique_ptr m_eclGraph; std::vector m_poreVolume; std::unique_ptr m_fldToolbox; bool m_hasUnifiedRestartFile; - QStringList m_restartFileNames; + std::vector m_singleRestartDataTimeSteps; + std::unique_ptr m_unifiedRestartData; }; @@ -150,9 +156,8 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI RigFlowDiagTimeStepResult result(m_eclipseCase->eclipseCaseData()->activeCellInfo(RifReaderInterface::MATRIX_RESULTS)->reservoirActiveCellCount()); caf::ProgressInfo progressInfo(8, "Calculating Flow Diagnostics"); - - if ( m_opmFldData.isNull() ) + if ( m_opmFlowDiagStaticData.isNull() ) { progressInfo.setProgressDescription("Grid access"); @@ -165,66 +170,77 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI QString initFileName = RifEclipseOutputFileTools::firstFileNameOfType(m_filesWithSameBaseName, ECL_INIT_FILE); - m_opmFldData = new RigOpmFldStaticData(gridFileName.toStdString(), + m_opmFlowDiagStaticData = new RigOpmFlowDiagStaticData(gridFileName.toStdString(), initFileName.toStdString()); progressInfo.incrementProgress(); progressInfo.setProgressDescription("Calculating Connectivities"); const Opm::FlowDiagnostics::ConnectivityGraph connGraph = - Opm::FlowDiagnostics::ConnectivityGraph{ static_cast(m_opmFldData->m_eclGraph.numCells()), - m_opmFldData->m_eclGraph.neighbours() }; + Opm::FlowDiagnostics::ConnectivityGraph{ static_cast(m_opmFlowDiagStaticData->m_eclGraph->numCells()), + m_opmFlowDiagStaticData->m_eclGraph->neighbours() }; progressInfo.incrementProgress(); progressInfo.setProgressDescription("Initialize Solver"); // Create the Toolbox. - m_opmFldData->m_fldToolbox.reset(new Opm::FlowDiagnostics::Toolbox{ connGraph }); - m_opmFldData->m_fldToolbox->assignPoreVolume( m_opmFldData->m_poreVolume); + m_opmFlowDiagStaticData->m_fldToolbox.reset(new Opm::FlowDiagnostics::Toolbox{ connGraph }); + m_opmFlowDiagStaticData->m_fldToolbox->assignPoreVolume( m_opmFlowDiagStaticData->m_poreVolume); // Look for unified restart file QString restartFileName = RifEclipseOutputFileTools::firstFileNameOfType(m_filesWithSameBaseName, ECL_UNIFIED_RESTART_FILE); if ( !restartFileName.isEmpty() ) { - m_opmFldData->m_eclGraph.assignFluxDataSource(restartFileName.toStdString()); - m_opmFldData->m_hasUnifiedRestartFile = true; + m_opmFlowDiagStaticData->m_unifiedRestartData.reset(new Opm::ECLRestartData(Opm::ECLRestartData(restartFileName.toStdString()))); + m_opmFlowDiagStaticData->m_hasUnifiedRestartFile = true; } else { + QStringList restartFileNames = RifEclipseOutputFileTools::filterFileNamesOfType(m_filesWithSameBaseName, ECL_RESTART_FILE); - m_opmFldData->m_restartFileNames = RifEclipseOutputFileTools::filterFileNamesOfType(m_filesWithSameBaseName, ECL_RESTART_FILE); - - size_t restartFileCount = static_cast(m_opmFldData->m_restartFileNames.size()); + size_t restartFileCount = static_cast(restartFileNames.size()); size_t maxTimeStepCount = m_eclipseCase->eclipseCaseData()->results(RifReaderInterface::MATRIX_RESULTS)->maxTimeStepCount(); - if (restartFileCount <= timeStepIndex && restartFileCount != maxTimeStepCount ) + if (restartFileCount <= timeStepIndex && restartFileCount != maxTimeStepCount ) { QMessageBox::critical(nullptr, "ResInsight", "Flow Diagnostics: Could not find all the restart files. Results will not be loaded."); return result; } - m_opmFldData->m_restartFileNames.sort(); // To make sure they are sorted in increasing *.X000N order. Hack. Should probably be actual time stored on file. - m_opmFldData->m_hasUnifiedRestartFile = false; + restartFileNames.sort(); // To make sure they are sorted in increasing *.X000N order. Hack. Should probably be actual time stored on file. + m_opmFlowDiagStaticData->m_hasUnifiedRestartFile = false; + + for (auto restartFileName : restartFileNames) + { + m_opmFlowDiagStaticData->m_singleRestartDataTimeSteps.push_back(Opm::ECLRestartData(restartFileName.toStdString())); + } } } progressInfo.setProgress(3); progressInfo.setProgressDescription("Assigning Flux Field"); - if ( ! m_opmFldData->m_hasUnifiedRestartFile ) + Opm::ECLRestartData* currentRestartData = nullptr; + + if ( ! m_opmFlowDiagStaticData->m_hasUnifiedRestartFile ) { - QString restartFileName = m_opmFldData->m_restartFileNames[static_cast(timeStepIndex)]; - m_opmFldData->m_eclGraph.assignFluxDataSource(restartFileName.toStdString()); + currentRestartData = &(m_opmFlowDiagStaticData->m_singleRestartDataTimeSteps[timeStepIndex]); } + else + { + currentRestartData = m_opmFlowDiagStaticData->m_unifiedRestartData.get(); + } + + CVF_ASSERT(currentRestartData); size_t resultIndexWithMaxTimeSteps = cvf::UNDEFINED_SIZE_T; m_eclipseCase->eclipseCaseData()->results(RifReaderInterface::MATRIX_RESULTS)->maxTimeStepCount(&resultIndexWithMaxTimeSteps); int reportStepNumber = m_eclipseCase->eclipseCaseData()->results(RifReaderInterface::MATRIX_RESULTS)->reportStepNumber(resultIndexWithMaxTimeSteps, timeStepIndex); - if ( ! m_opmFldData->m_eclGraph.selectReportStep(reportStepNumber) ) + if ( !currentRestartData->selectReportStep(reportStepNumber) ) { QMessageBox::critical(nullptr, "ResInsight", "Flow Diagnostics: Could not find the requested timestep in the result file. Results will not be loaded."); return result; @@ -235,22 +251,25 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI Opm::FlowDiagnostics::CellSetValues sumWellFluxPrCell; { - Opm::FlowDiagnostics::ConnectionValues connectionsVals = RigFlowDiagInterfaceTools::extractFluxField(m_opmFldData->m_eclGraph, false); + Opm::FlowDiagnostics::ConnectionValues connectionsVals = RigFlowDiagInterfaceTools::extractFluxField(*(m_opmFlowDiagStaticData->m_eclGraph), + *currentRestartData, + false); - m_opmFldData->m_fldToolbox->assignConnectionFlux(connectionsVals); + m_opmFlowDiagStaticData->m_fldToolbox->assignConnectionFlux(connectionsVals); progressInfo.incrementProgress(); Opm::ECLWellSolution wsol = Opm::ECLWellSolution{-1.0 , false}; - const std::vector well_fluxes = - wsol.solution(m_opmFldData->m_eclGraph.rawResultData(), m_opmFldData->m_eclGraph.numGrids()); + std::vector gridNames = m_opmFlowDiagStaticData->m_eclGraph->activeGrids(); - sumWellFluxPrCell = RigFlowDiagInterfaceTools::extractWellFlows(m_opmFldData->m_eclGraph, well_fluxes); + const std::vector well_fluxes = wsol.solution(*currentRestartData, gridNames); - m_opmFldData->m_fldToolbox->assignInflowFlux(sumWellFluxPrCell); + sumWellFluxPrCell = RigFlowDiagInterfaceTools::extractWellFlows(*(m_opmFlowDiagStaticData->m_eclGraph), well_fluxes); - // Filter connection cells with inconsistent well in flow direction (Hack, we should do something better) + m_opmFlowDiagStaticData->m_fldToolbox->assignInflowFlux(sumWellFluxPrCell); + + // Start Hack: Filter connection cells with inconsistent well in flow direction (Hack, we should do something better) for ( auto& tracerCellIdxsPair: injectorTracers ) { @@ -259,6 +278,8 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI for (int activeCellIdx : tracerCellIdxsPair.second) { auto activeCellIdxFluxPair = sumWellFluxPrCell.find(activeCellIdx); + CVF_TIGHT_ASSERT(activeCellIdxFluxPair != sumWellFluxPrCell.end()); + if (activeCellIdxFluxPair->second > 0 ) { filteredCellIndices.push_back(activeCellIdx); @@ -275,6 +296,8 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI for (int activeCellIdx : tracerCellIdxsPair.second) { auto activeCellIdxFluxPair = sumWellFluxPrCell.find(activeCellIdx); + CVF_TIGHT_ASSERT(activeCellIdxFluxPair != sumWellFluxPrCell.end()); + if (activeCellIdxFluxPair->second < 0 ) { filteredCellIndices.push_back(activeCellIdx); @@ -282,6 +305,8 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI } if (tracerCellIdxsPair.second.size() != filteredCellIndices.size()) tracerCellIdxsPair.second = filteredCellIndices; } + + // End Hack } progressInfo.incrementProgress(); @@ -299,7 +324,7 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI std::unique_ptr injectorSolution; try { - injectorSolution.reset(new Toolbox::Forward( m_opmFldData->m_fldToolbox->computeInjectionDiagnostics(injectorCellSets))); + injectorSolution.reset(new Toolbox::Forward( m_opmFlowDiagStaticData->m_fldToolbox->computeInjectionDiagnostics(injectorCellSets))); } catch (const std::exception& e) { @@ -329,7 +354,7 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI std::unique_ptr producerSolution; try { - producerSolution.reset(new Toolbox::Reverse(m_opmFldData->m_fldToolbox->computeProductionDiagnostics(prodjCellSets))); + producerSolution.reset(new Toolbox::Reverse(m_opmFlowDiagStaticData->m_fldToolbox->computeProductionDiagnostics(prodjCellSets))); } catch ( const std::exception& e ) { @@ -375,7 +400,7 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI { Graph flowCapStorCapCurve = flowCapacityStorageCapacityCurve(*(injectorSolution.get()), *(producerSolution.get()), - m_opmFldData->m_poreVolume); + m_opmFlowDiagStaticData->m_poreVolume); result.setFlowCapStorageCapCurve(flowCapStorCapCurve); result.setSweepEfficiencyCurve(sweepEfficiency(flowCapStorCapCurve)); diff --git a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.h b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.h index 653d9e54d2..81a154b715 100644 --- a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.h +++ b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.h @@ -70,7 +70,7 @@ private: class RigEclipseCaseData; -class RigOpmFldStaticData; +class RigOpmFlowDiagStaticData; class RigFlowDiagSolverInterface : public cvf::Object { @@ -85,7 +85,7 @@ public: private: RimEclipseResultCase * m_eclipseCase; - cvf::ref m_opmFldData; + cvf::ref m_opmFlowDiagStaticData; }; diff --git a/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp b/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp index 4810098233..5279eff299 100644 --- a/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp +++ b/ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp @@ -2,6 +2,7 @@ const std::string casePath = "\\\\csfiles\\Store\\ProjectData\\StatoilReservoir\\ReferenceCases\\simple_FlowDiag_Model\\"; +/* #include "exampleSetup.hpp" TEST(opm_flowdiagnostics_test, basic_construction) @@ -9,12 +10,14 @@ TEST(opm_flowdiagnostics_test, basic_construction) try { - + Opm::ECLRestartData rstrt(casePath + "SIMPLE.UNRST"); + Opm::ECLInitFileData initData(casePath + "SIMPLE.INIT"); + Opm::ECLGraph graph = Opm::ECLGraph::load(casePath + "SIMPLE.EGRID", - casePath + "SIMPLE.INIT"); - graph.assignFluxDataSource(casePath + "SIMPLE.UNRST"); + initData); + int step = 2; - if ( ! graph.selectReportStep(step) ) + if ( ! rstrt.selectReportStep(step) ) { std::ostringstream os; @@ -44,3 +47,4 @@ TEST(opm_flowdiagnostics_test, basic_construction) std::cerr << "Caught exception: " << e.what() << '\n'; } } +*/