///////////////////////////////////////////////////////////////////////////////// // // Copyright (C) 2017- Statoil ASA // // ResInsight is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY // WARRANTY; without even the implied warranty of MERCHANTABILITY or // FITNESS FOR A PARTICULAR PURPOSE. // // See the GNU General Public License at // for more details. // ///////////////////////////////////////////////////////////////////////////////// #include "RigNumberOfFloodedPoreVolumesCalculator.h" #include "RiaPorosityModel.h" #include "RigActiveCellInfo.h" #include "RigCaseCellResultsData.h" #include "RigEclipseCaseData.h" #include "RigMainGrid.h" #include "RigReservoirBuilderMock.h" #include "RimEclipseCase.h" #include "RimReservoirCellResultsStorage.h" #include #include #include "cafProgressInfo.h" //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- RigNumberOfFloodedPoreVolumesCalculator::RigNumberOfFloodedPoreVolumesCalculator(RimEclipseCase* caseToApply, const std::vector tracerNames) { RigMainGrid* mainGrid = caseToApply->eclipseCaseData()->mainGrid(); RigEclipseCaseData* eclipseCaseData = caseToApply->eclipseCaseData(); RigCaseCellResultsData* gridCellResults = caseToApply->results(RiaDefines::MATRIX_MODEL); RigActiveCellInfo* actCellInfo = caseToApply->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL); size_t resultCellCount = actCellInfo->reservoirCellResultCount(); size_t timeStepCount = caseToApply->eclipseCaseData()->results(RiaDefines::MATRIX_MODEL)->maxTimeStepCount(); size_t totalProgress = tracerNames.size() + 8 + timeStepCount + 2* timeStepCount; caf::ProgressInfo progress(totalProgress, "Calculating number of flooded mobile pore volumes." ); progress.setProgressDescription("Loading required results"); // PORV const std::vector* porvResults = nullptr; std::vector porvActiveCellsResultStorage; porvResults = RigCaseCellResultsData::getResultIndexableStaticResult(actCellInfo, gridCellResults, "PORV", porvActiveCellsResultStorage); progress.incrementProgress(); // SWCR if defined const std::vector* swcrResults = nullptr; swcrResults = RigCaseCellResultsData::getResultIndexableStaticResult(actCellInfo, gridCellResults, "SWCR", porvActiveCellsResultStorage); progress.incrementProgress(); std::vector scalarResultIndexTracers; for (QString tracerName : tracerNames) { scalarResultIndexTracers.push_back(gridCellResults->findOrLoadScalarResult(RiaDefines::DYNAMIC_NATIVE, tracerName)); progress.incrementProgress(); } std::vector > summedTracersAtAllTimesteps; //TODO: Option for Oil and Gas instead of water size_t scalarResultIndexFlowrateI = gridCellResults->findOrLoadScalarResult(RiaDefines::DYNAMIC_NATIVE, "FLRWATI+"); progress.incrementProgress(); size_t scalarResultIndexFlowrateJ = gridCellResults->findOrLoadScalarResult(RiaDefines::DYNAMIC_NATIVE, "FLRWATJ+"); progress.incrementProgress(); size_t scalarResultIndexFlowrateK = gridCellResults->findOrLoadScalarResult(RiaDefines::DYNAMIC_NATIVE, "FLRWATK+"); progress.incrementProgress(); std::vector* > flowrateIatAllTimeSteps; std::vector* > flowrateJatAllTimeSteps; std::vector* > flowrateKatAllTimeSteps; RigNNCData* nncData = eclipseCaseData->mainGrid()->nncData(); const std::vector connections = nncData->connections(); progress.incrementProgress(); //TODO: oil or gas flowrate std::vector* > flowrateNNCatAllTimeSteps; QString nncConnectionProperty = mainGrid->nncData()->propertyNameFluxWat(); progress.incrementProgress(); std::vector daysSinceSimulationStart = caseToApply->eclipseCaseData()->results(RiaDefines::MATRIX_MODEL)->daysSinceSimulationStart(); progress.incrementProgress(); for (size_t timeStep = 0; timeStep < daysSinceSimulationStart.size(); timeStep++) { const std::vector* flowrateI = nullptr; if (scalarResultIndexFlowrateI != cvf::UNDEFINED_SIZE_T) { flowrateI = &(eclipseCaseData->results(RiaDefines::MATRIX_MODEL)->cellScalarResults(scalarResultIndexFlowrateI, timeStep)); } flowrateIatAllTimeSteps.push_back(flowrateI); const std::vector* flowrateJ = nullptr; if (scalarResultIndexFlowrateJ != cvf::UNDEFINED_SIZE_T) { flowrateJ = &(eclipseCaseData->results(RiaDefines::MATRIX_MODEL)->cellScalarResults(scalarResultIndexFlowrateJ, timeStep)); } flowrateJatAllTimeSteps.push_back(flowrateJ); const std::vector* flowrateK = nullptr; if (scalarResultIndexFlowrateK != cvf::UNDEFINED_SIZE_T) { flowrateK = &(eclipseCaseData->results(RiaDefines::MATRIX_MODEL)->cellScalarResults(scalarResultIndexFlowrateK, timeStep)); } flowrateKatAllTimeSteps.push_back(flowrateK); size_t nativeTimeStepIndex = caseToApply->uiToNativeTimeStepIndex(timeStep); const std::vector* connectionFlowrate = nncData->dynamicConnectionScalarResultByName(nncConnectionProperty, nativeTimeStepIndex); flowrateNNCatAllTimeSteps.push_back(connectionFlowrate); //sum all tracers at current timestep std::vector summedTracerValues(resultCellCount); for (size_t tracerIndex : scalarResultIndexTracers) { if (tracerIndex != cvf::UNDEFINED_SIZE_T) { const std::vector* tracerResult = &(eclipseCaseData->results(RiaDefines::MATRIX_MODEL)->cellScalarResults(tracerIndex, timeStep)); for (size_t i = 0; i < summedTracerValues.size(); i++) { summedTracerValues[i] += tracerResult->at(i); } } } summedTracersAtAllTimesteps.push_back(summedTracerValues); progress.incrementProgress(); } progress.setNextProgressIncrement(2*timeStepCount); progress.setProgressDescription("Calculating"); calculate(mainGrid, caseToApply, daysSinceSimulationStart, porvResults, swcrResults, flowrateIatAllTimeSteps, flowrateJatAllTimeSteps, flowrateKatAllTimeSteps, connections, flowrateNNCatAllTimeSteps, summedTracersAtAllTimesteps); } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- std::vector>& RigNumberOfFloodedPoreVolumesCalculator::numberOfFloodedPorevolumes() { return m_cumWinflowPVAllTimeSteps; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RigNumberOfFloodedPoreVolumesCalculator::calculate(RigMainGrid* mainGrid, RimEclipseCase* caseToApply, std::vector daysSinceSimulationStart, const std::vector* porvResultsActiveCellsOnly, const std::vector* swcrResults, std::vector* > flowrateIatAllTimeSteps, std::vector* > flowrateJatAllTimeSteps, std::vector* > flowrateKatAllTimeSteps, const std::vector connections, std::vector* > flowrateNNCatAllTimeSteps, std::vector > summedTracersAtAllTimesteps) { //size_t totalNumberOfCells = mainGrid->globalCellArray().size(); RigActiveCellInfo* actCellInfo = caseToApply->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL); size_t resultCellCount = actCellInfo->reservoirCellResultCount(); caf::ProgressInfo progress(2*daysSinceSimulationStart.size() , ""); std::vector> cellQwInAtAllTimeSteps; std::vector cellQwInTimeStep0(resultCellCount); cellQwInAtAllTimeSteps.push_back(cellQwInTimeStep0); for (size_t timeStep = 1; timeStep < daysSinceSimulationStart.size(); timeStep++) { std::vector totoalFlowrateIntoCell(resultCellCount); //brukt result celle index / active antall i stedet if ( flowrateIatAllTimeSteps[timeStep-1] != nullptr && flowrateJatAllTimeSteps[timeStep-1] != nullptr && flowrateKatAllTimeSteps[timeStep-1] != nullptr) { const std::vector* flowrateI = flowrateIatAllTimeSteps[timeStep-1]; const std::vector* flowrateJ = flowrateJatAllTimeSteps[timeStep-1]; const std::vector* flowrateK = flowrateKatAllTimeSteps[timeStep-1]; if (flowrateI->size() > 0 && flowrateJ->size() > 0 && flowrateK->size() > 0) { distributeNeighbourCellFlow(mainGrid, caseToApply, summedTracersAtAllTimesteps[timeStep-1], flowrateI, flowrateJ, flowrateK, totoalFlowrateIntoCell); } } const std::vector* flowrateNNC = flowrateNNCatAllTimeSteps[timeStep-1]; if (flowrateNNC->size() > 0) { distributeNNCflow(connections, caseToApply, summedTracersAtAllTimesteps[timeStep-1], flowrateNNC, totoalFlowrateIntoCell); } std::vector CellQwIn(resultCellCount); double daysSinceSimStartNow = daysSinceSimulationStart[timeStep]; double daysSinceSimStartLastTimeStep = daysSinceSimulationStart[timeStep - 1]; double deltaT = daysSinceSimStartNow - daysSinceSimStartLastTimeStep; for (size_t cellResultIndex = 0; cellResultIndex < resultCellCount; cellResultIndex++) { CellQwIn[cellResultIndex] = cellQwInAtAllTimeSteps[timeStep - 1][cellResultIndex] + (totoalFlowrateIntoCell[cellResultIndex]) * deltaT; } cellQwInAtAllTimeSteps.push_back(CellQwIn); progress.incrementProgress(); } //Calculate number-of-cell-PV flooded std::vector cumWinflowPVTimeStep0(resultCellCount); m_cumWinflowPVAllTimeSteps.clear(); m_cumWinflowPVAllTimeSteps.push_back(cumWinflowPVTimeStep0); for (size_t timeStep = 1; timeStep < daysSinceSimulationStart.size(); timeStep++) { std::vector cumWinflowPV(resultCellCount); for (size_t cellResultIndex = 0; cellResultIndex < resultCellCount; cellResultIndex++) { double scaledPoreVolume = porvResultsActiveCellsOnly->at(cellResultIndex); if (swcrResults != nullptr && swcrResults->size() == resultCellCount) { scaledPoreVolume = scaledPoreVolume * (1 - swcrResults->at(cellResultIndex)); } cumWinflowPV[cellResultIndex] = cellQwInAtAllTimeSteps[timeStep][cellResultIndex] / scaledPoreVolume; } m_cumWinflowPVAllTimeSteps.push_back(cumWinflowPV); progress.incrementProgress(); } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RigNumberOfFloodedPoreVolumesCalculator::distributeNNCflow(std::vector connections, RimEclipseCase* caseToApply, std::vector summedTracerValues, const std::vector* flowrateNNC, std::vector &flowrateIntoCell) { RigActiveCellInfo* actCellInfo = caseToApply->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL); for (size_t connectionIndex = 0; connectionIndex < connections.size(); connectionIndex++) { RigConnection connection = connections[connectionIndex]; double connectionValue = flowrateNNC->at(connectionIndex); size_t cell1Index = connection.m_c1GlobIdx; size_t cell1ResultIndex = actCellInfo->cellResultIndex(cell1Index); size_t cell2Index = connection.m_c2GlobIdx; size_t cell2ResultIndex = actCellInfo->cellResultIndex(cell2Index); if (connectionValue > 0) { //Flow out of cell with cell1index, into cell cell2index flowrateIntoCell[cell2ResultIndex] += connectionValue * summedTracerValues[cell1ResultIndex]; } else if (connectionValue < 0) { //flow out of cell with cell2index, into cell cell1index flowrateIntoCell[cell1ResultIndex] += -1.0*connectionValue * summedTracerValues[cell2ResultIndex]; } } } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- void RigNumberOfFloodedPoreVolumesCalculator::distributeNeighbourCellFlow(RigMainGrid* mainGrid, RimEclipseCase* caseToApply, std::vector summedTracerValues, const std::vector* flrWatResultI, const std::vector* flrWatResultJ, const std::vector* flrWatResultK, std::vector &totalFlowrateIntoCell) { RigActiveCellInfo* actCellInfo = caseToApply->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL); for (size_t globalCellIndex = 0; globalCellIndex < mainGrid->globalCellArray().size(); globalCellIndex++) { if (!actCellInfo->isActive(globalCellIndex)) continue; const RigCell& cell = mainGrid->globalCellArray()[globalCellIndex]; RigGridBase* hostGrid = cell.hostGrid(); size_t gridLocalCellIndex = cell.gridLocalCellIndex(); size_t cellResultIndex = actCellInfo->cellResultIndex(globalCellIndex); size_t i, j, k; hostGrid->ijkFromCellIndex(gridLocalCellIndex, &i, &j, &k); if (i < (hostGrid->cellCountI()-1)) { size_t gridLocalCellIndexPosINeighbour = hostGrid->cellIndexFromIJK(i + 1, j, k); size_t reservoirCellIndexPosINeighbour = hostGrid->reservoirCellIndex(gridLocalCellIndexPosINeighbour); size_t cellResultIndexPosINeighbour = actCellInfo->cellResultIndex(reservoirCellIndexPosINeighbour); if (!actCellInfo->isActive(reservoirCellIndexPosINeighbour)) continue; if (hostGrid->cell(gridLocalCellIndexPosINeighbour).subGrid() != NULL) { //subgrid exists in cell, will be handled though NNCs continue; } if (flrWatResultI->at(cellResultIndex) > 0) { //Flow out of cell globalCellIndex, into cell i+1 totalFlowrateIntoCell[cellResultIndexPosINeighbour] += flrWatResultI->at(cellResultIndex) * summedTracerValues[cellResultIndex]; } else if (flrWatResultI->at(cellResultIndex) < 0) { //Flow into cell globelCellIndex, from cell i+1 totalFlowrateIntoCell[cellResultIndex] += (-1.0) * flrWatResultI->at(cellResultIndex) * summedTracerValues[cellResultIndexPosINeighbour]; } } if (j < (hostGrid->cellCountJ()-1)) { size_t gridLocalCellIndexPosJNeighbour = hostGrid->cellIndexFromIJK(i, j + 1, k); size_t reservoirCellIndexPosJNeighbour = hostGrid->reservoirCellIndex(gridLocalCellIndexPosJNeighbour); size_t cellResultIndexPosJNeighbour = actCellInfo->cellResultIndex(reservoirCellIndexPosJNeighbour); if (!actCellInfo->isActive(reservoirCellIndexPosJNeighbour)) continue; if (hostGrid->cell(gridLocalCellIndexPosJNeighbour).subGrid() != NULL) { //subgrid exists in cell, will be handled though NNCs continue; } if (flrWatResultJ->at(cellResultIndex) > 0) { //Flow out of cell globalCellIndex, into cell i+1 totalFlowrateIntoCell[cellResultIndexPosJNeighbour] += flrWatResultJ->at(cellResultIndex) * summedTracerValues[cellResultIndex]; } else if (flrWatResultJ->at(cellResultIndex) < 0) { //Flow into cell globelCellIndex, from cell i+1 totalFlowrateIntoCell[cellResultIndex] += (-1.0) * flrWatResultJ->at(cellResultIndex) * summedTracerValues[cellResultIndexPosJNeighbour]; } } if (k < (hostGrid->cellCountK()-1)) { size_t gridLocalCellIndexPosKNeighbour = hostGrid->cellIndexFromIJK(i, j, k + 1); size_t reservoirCellIndexPosKNeighbour = hostGrid->reservoirCellIndex(gridLocalCellIndexPosKNeighbour); size_t cellResultIndexPosKNeighbour = actCellInfo->cellResultIndex(reservoirCellIndexPosKNeighbour); if (!actCellInfo->isActive(reservoirCellIndexPosKNeighbour)) continue; if (hostGrid->cell(gridLocalCellIndexPosKNeighbour).subGrid() != NULL) { //subgrid exists in cell, will be handled though NNCs continue; } if (flrWatResultK->at(cellResultIndex) > 0) { //Flow out of cell globalCellIndex, into cell i+1 totalFlowrateIntoCell[cellResultIndexPosKNeighbour] += flrWatResultK->at(cellResultIndex) * summedTracerValues[cellResultIndex]; } else if (flrWatResultK->at(cellResultIndex) < 0) { //Flow into cell globelCellIndex, from cell i+1 totalFlowrateIntoCell[cellResultIndex] += (-1.0) * flrWatResultK->at(cellResultIndex) * summedTracerValues[cellResultIndexPosKNeighbour]; } } } }