#1698 Use opm-flowdiagnostics calculation of fluxes if not present in Eclipse restart file

ApplicationCode/ProjectDataModel/RimEclipseResultCase.cpp

@@ -136,14 +136,11 @@ bool RimEclipseResultCase::openEclipseGridFile()
     m_gridAndWellDataIsReadFromFile = true;
     m_activeCellInfoIsReadFromFile = true;
 
-    if (eclipseCaseData()->results(RifReaderInterface::MATRIX_RESULTS)->hasFlowDiagUsableFluxes())
+    m_flowDagSolverInterface = new RigFlowDiagSolverInterface(this);
-    {
-        m_flowDagSolverInterface = new RigFlowDiagSolverInterface(this);
 
-        if (m_flowDiagSolutions.size() == 0)
+    if (m_flowDiagSolutions.size() == 0)
-        {
+    {
-            m_flowDiagSolutions.push_back(new RimFlowDiagSolution());
+        m_flowDiagSolutions.push_back(new RimFlowDiagSolution());
-        }
     }
     
     return true;

ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h

@@ -81,6 +81,24 @@ namespace RigFlowDiagInterfaceTools {
         return extractFluxField(G, getFlux);
     }
 
+    inline Opm::FlowDiagnostics::ConnectionValues
+        calculateFluxField(const Opm::ECLGraph& G,
+                           const Opm::ECLInitFileData& init,
+                           const Opm::ECLRestartData& rstrt)
+    {
+        auto satfunc = Opm::ECLSaturationFunc(G, init);
+
+        Opm::ECLFluxCalc calc(G, std::move(satfunc));
+
+        auto getFlux = [&calc, &rstrt]
+        (const Opm::ECLPhaseIndex p)
+        {
+            return calc.flux(rstrt, p);
+        };
+
+        return extractFluxField(G, getFlux);
+    }
+
     template <class WellFluxes>
     Opm::FlowDiagnostics::CellSetValues
         extractWellFlows(const Opm::ECLGraph& G,

ApplicationCode/ReservoirDataModel/RigFlowDiagResults.cpp

@@ -698,5 +698,5 @@ RigFlowDiagSolverInterface::FlowCharacteristicsResultFrame RigFlowDiagResults::f
     const std::vector<double>* injectorResults = resultValues(injectorAddress, frameIndex);
     const std::vector<double>* producerResults = resultValues(producerAddress, frameIndex);
 
-    return solverInterface()->calculateFlowCharacteristics(*injectorResults, *producerResults, max_pv_fraction);
+    return solverInterface()->calculateFlowCharacteristics(injectorResults, producerResults, max_pv_fraction);
 }

ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp

@@ -164,14 +164,11 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI
         // Get set of files
         QString gridFileName = m_eclipseCase->gridFileName();
 
-        QStringList m_filesWithSameBaseName;
+        std::string initFileName = getInitFileName();
 
-        if ( !RifEclipseOutputFileTools::findSiblingFilesWithSameBaseName(gridFileName, &m_filesWithSameBaseName) ) return result;
+        if (initFileName.empty()) return result;
 
-        QString initFileName = RifEclipseOutputFileTools::firstFileNameOfType(m_filesWithSameBaseName, ECL_INIT_FILE);
+        m_opmFlowDiagStaticData = new RigOpmFlowDiagStaticData(gridFileName.toStdString(), initFileName);
-
-        m_opmFlowDiagStaticData = new RigOpmFlowDiagStaticData(gridFileName.toStdString(),
-                                               initFileName.toStdString());
 
         progressInfo.incrementProgress();
         progressInfo.setProgressDescription("Calculating Connectivities");
@@ -189,6 +186,9 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI
         m_opmFlowDiagStaticData->m_fldToolbox->assignPoreVolume( m_opmFlowDiagStaticData->m_poreVolume);
 
         // Look for unified restart file
+        QStringList m_filesWithSameBaseName;
+
+        if ( !RifEclipseOutputFileTools::findSiblingFilesWithSameBaseName(gridFileName, &m_filesWithSameBaseName) ) return result;
 
         QString restartFileName = RifEclipseOutputFileTools::firstFileNameOfType(m_filesWithSameBaseName, ECL_UNIFIED_RESTART_FILE);
         if ( !restartFileName.isEmpty() )
@@ -251,10 +251,19 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI
     Opm::FlowDiagnostics::CellSetValues sumWellFluxPrCell;
 
     {
-        Opm::FlowDiagnostics::ConnectionValues connectionsVals = RigFlowDiagInterfaceTools::extractFluxFieldFromRestartFile(*(m_opmFlowDiagStaticData->m_eclGraph),
+        if (m_eclipseCase->eclipseCaseData()->results(RifReaderInterface::MATRIX_RESULTS)->hasFlowDiagUsableFluxes())
-                                                                                                                            *currentRestartData);
+        {
+            Opm::FlowDiagnostics::ConnectionValues connectionsVals = RigFlowDiagInterfaceTools::extractFluxFieldFromRestartFile(*(m_opmFlowDiagStaticData->m_eclGraph),
+                *currentRestartData);
+            m_opmFlowDiagStaticData->m_fldToolbox->assignConnectionFlux(connectionsVals);
+        }
+        else
+        {
+            Opm::ECLInitFileData init(getInitFileName());
+            Opm::FlowDiagnostics::ConnectionValues connectionVals = RigFlowDiagInterfaceTools::calculateFluxField((*m_opmFlowDiagStaticData->m_eclGraph), init, *currentRestartData);
+            m_opmFlowDiagStaticData->m_fldToolbox->assignConnectionFlux(connectionVals);
+        }
 
-        m_opmFlowDiagStaticData->m_fldToolbox->assignConnectionFlux(connectionsVals);
 
         progressInfo.incrementProgress();
 
@@ -402,17 +411,22 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-RigFlowDiagSolverInterface::FlowCharacteristicsResultFrame RigFlowDiagSolverInterface::calculateFlowCharacteristics(const std::vector<double>& injector_tof,
+RigFlowDiagSolverInterface::FlowCharacteristicsResultFrame RigFlowDiagSolverInterface::calculateFlowCharacteristics(const std::vector<double>* injector_tof,
-                                                                                                                    const std::vector<double>& producer_tof,
+                                                                                                                    const std::vector<double>* producer_tof,
                                                                                                                     double max_pv_fraction)
 {
     using namespace Opm::FlowDiagnostics;
     RigFlowDiagSolverInterface::FlowCharacteristicsResultFrame result;
 
+    if (injector_tof == nullptr || producer_tof == nullptr)
+    {
+        return result;
+    }
+
     try
     {
-        Graph flowCapStorCapCurve = flowCapacityStorageCapacityCurve(injector_tof,
+        Graph flowCapStorCapCurve = flowCapacityStorageCapacityCurve(*injector_tof,
-                                                                     producer_tof,
+                                                                     *producer_tof,
                                                                      m_opmFlowDiagStaticData->m_poreVolume,
                                                                      max_pv_fraction);
 
@@ -428,6 +442,22 @@ RigFlowDiagSolverInterface::FlowCharacteristicsResultFrame RigFlowDiagSolverInte
     return result;
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+std::string RigFlowDiagSolverInterface::getInitFileName() const
+{
+    QString gridFileName = m_eclipseCase->gridFileName();
+
+    QStringList m_filesWithSameBaseName;
+
+    if (!RifEclipseOutputFileTools::findSiblingFilesWithSameBaseName(gridFileName, &m_filesWithSameBaseName)) return std::string();
+
+    QString initFileName = RifEclipseOutputFileTools::firstFileNameOfType(m_filesWithSameBaseName, ECL_INIT_FILE);
+
+    return initFileName.toStdString();
+}
+
 
 //--------------------------------------------------------------------------------------------------
 /// 

ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.h

@@ -83,11 +83,13 @@ public:
                                         std::map<std::string, std::vector<int> > injectorTracers, 
                                         std::map<std::string, std::vector<int> > producerTracers);
 
-    FlowCharacteristicsResultFrame calculateFlowCharacteristics(const std::vector<double>& injector_tof,
+    FlowCharacteristicsResultFrame calculateFlowCharacteristics(const std::vector<double>* injector_tof,
-                                                                const std::vector<double>& producer_tof,
+                                                                const std::vector<double>* producer_tof,
                                                                 double max_pv_fraction);
 
 private:
+    std::string getInitFileName() const;
+
     RimEclipseResultCase * m_eclipseCase;
 
     cvf::ref<RigOpmFlowDiagStaticData> m_opmFlowDiagStaticData;