Contour Statistics Map: Compute results

ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp

@@ -142,6 +142,7 @@ double RimEclipseContourMapProjection::sampleSpacing() const
 
     return 0.0;
 }
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ProjectDataModel/RimStatisticsContourMap.cpp

@@ -19,17 +19,27 @@
 #include "RimStatisticsContourMap.h"
 
 #include "RiaLogging.h"
+#include "RiaStatisticsTools.h"
 
 #include "RigCaseCellResultsData.h"
+#include "RigContourMapGrid.h"
 #include "RigEclipseCaseData.h"
-#include "RimEclipseCaseEnsemble.h"
-#include "RimProject.h"
+#include "RigEclipseContourMapProjection.h"
+#include "RigEclipseResultAddress.h"
+#include "RigMainGrid.h"
+#include "RigStatisticsMath.h"
 
+#include "RimEclipseCase.h"
+#include "RimEclipseCaseEnsemble.h"
+#include "RimEclipseContourMapProjection.h"
 #include "RimEclipseResultDefinition.h"
+#include "RimProject.h"
 
 #include "cafPdmUiPushButtonEditor.h"
 #include "cafProgressInfo.h"
 
+#include <limits>
+
 CAF_PDM_SOURCE_INIT( RimStatisticsContourMap, "RimStatisticalContourMap" );
 
 //--------------------------------------------------------------------------------------------------
@@ -115,4 +125,85 @@ void RimStatisticsContourMap::initAfterRead()
 void RimStatisticsContourMap::computeStatistics()
 {
     RiaLogging::info( "Computing statistics" );
+    auto ensemble = firstAncestorOrThisOfType<RimEclipseCaseEnsemble>();
+    if ( !ensemble ) return;
+
+    if ( ensemble->cases().empty() ) return;
+
+    RimEclipseCase* firstEclipseCase = ensemble->cases().front();
+    firstEclipseCase->ensureReservoirCaseIsOpen();
+
+    double                                         relativeSampleSpacing = 0.9;
+    int                                            timeStep              = 0;
+    RigContourMapCalculator::ResultAggregationEnum resultAggregation = RigContourMapCalculator::ResultAggregationEnum::RESULTS_MEAN_VALUE;
+
+    cvf::BoundingBox gridBoundingBox = firstEclipseCase->activeCellsBoundingBox();
+
+    auto computeSampleSpacing = []( auto ec, double relativeSampleSpacing )
+    {
+        if ( ec )
+        {
+            if ( auto mainGrid = ec->mainGrid() )
+            {
+                return relativeSampleSpacing * mainGrid->characteristicIJCellSize();
+            }
+        }
+
+        return 0.0;
+    };
+
+    double sampleSpacing = computeSampleSpacing( firstEclipseCase, relativeSampleSpacing );
+
+    auto contourMapGrid = std::make_unique<RigContourMapGrid>( gridBoundingBox, sampleSpacing );
+
+    std::vector<std::vector<double>> results;
+    for ( RimEclipseCase* eclipseCase : ensemble->cases() )
+    {
+        if ( eclipseCase->ensureReservoirCaseIsOpen() )
+        {
+            RiaLogging::info( QString( "Grid: %1" ).arg( eclipseCase->caseUserDescription() ) );
+
+            auto eclipseCaseData = eclipseCase->eclipseCaseData();
+            auto resultData      = eclipseCaseData->results( RiaDefines::PorosityModelType::MATRIX_MODEL );
+
+            RigEclipseContourMapProjection contourMapProjection( *contourMapGrid, *eclipseCaseData, *resultData );
+            contourMapProjection.generateGridMapping( resultAggregation, {} );
+
+            std::vector<double> result =
+                contourMapProjection.generateResults( m_resultDefinition()->eclipseResultAddress(), resultAggregation, timeStep );
+            results.push_back( result );
+        }
+    }
+
+    if ( !results.empty() )
+    {
+        int                 nCells = static_cast<int>( results[0].size() );
+        std::vector<double> p10Results( nCells, std::numeric_limits<double>::infinity() );
+        std::vector<double> p50Results( nCells, std::numeric_limits<double>::infinity() );
+        std::vector<double> p90Results( nCells, std::numeric_limits<double>::infinity() );
+        std::vector<double> meanResults( nCells, std::numeric_limits<double>::infinity() );
+        std::vector<double> minResults( nCells, std::numeric_limits<double>::infinity() );
+        std::vector<double> maxResults( nCells, std::numeric_limits<double>::infinity() );
+
+#pragma omp parallel for
+        for ( int i = 0; i < nCells; i++ )
+        {
+            size_t              numSamples = results.size();
+            std::vector<double> samples( numSamples, 0.0 );
+            for ( size_t s = 0; s < numSamples; s++ )
+                samples[s] = results[s][i];
+
+            double p10, p50, p90, mean;
+
+            RigStatisticsMath::calculateStatisticsCurves( samples, &p10, &p50, &p90, &mean, RigStatisticsMath::PercentileStyle::SWITCHED );
+
+            p10Results[i]  = p10;
+            p50Results[i]  = p50;
+            p90Results[i]  = p90;
+            meanResults[i] = mean;
+
+            minResults[i] = RiaStatisticsTools::minimumValue( samples );
+            maxResults[i] = RiaStatisticsTools::maximumValue( samples );
+        }
+    }
 }

ApplicationLibCode/ReservoirDataModel/RigContourMapCalculator.cpp

@@ -300,7 +300,7 @@ std::vector<RigContourMapCalculator::CellIndexAndResult>
     auto cellGridIdxVisibility = contourMapProjection.getCellVisibility();
     for ( size_t globalCellIdx : allCellIndices )
     {
-        if ( ( *cellGridIdxVisibility )[globalCellIdx] )
+        if ( cellGridIdxVisibility.isNull() || ( *cellGridIdxVisibility )[globalCellIdx] )
         {
             kLayerCellIndexVector[contourMapProjection.kLayer( globalCellIdx )].push_back( globalCellIdx );
         }

ApplicationLibCode/ReservoirDataModel/RigContourMapProjection.cpp

@@ -18,6 +18,7 @@
 
 #include "RigContourMapProjection.h"
 
+#include "RiaStatisticsTools.h"
 #include "RiaWeightedMeanCalculator.h"
 
 #include "RigContourMapCalculator.h"
@@ -28,6 +29,7 @@
 
 #include <algorithm>
 #include <array>
+#include <limits>
 
 //--------------------------------------------------------------------------------------------------
 ///
@@ -210,16 +212,8 @@ double RigContourMapProjection::calculateValueInMapCell( unsigned int
 //--------------------------------------------------------------------------------------------------
 double RigContourMapProjection::maxValue( const std::vector<double>& aggregatedResults )
 {
-    double maxV = -std::numeric_limits<double>::infinity();
-
-    for ( size_t index = 0; index < aggregatedResults.size(); ++index )
-    {
-        if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
-        {
-            maxV = std::max( maxV, aggregatedResults[index] );
-        }
-    }
-    return maxV;
+    double maxV = RiaStatisticsTools::maximumValue( aggregatedResults );
+    return maxV != -std::numeric_limits<double>::max() ? maxV : -std::numeric_limits<double>::infinity();
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -227,16 +221,8 @@ double RigContourMapProjection::maxValue( const std::vector<double>& aggregatedR
 //--------------------------------------------------------------------------------------------------
 double RigContourMapProjection::minValue( const std::vector<double>& aggregatedResults )
 {
-    double minV = std::numeric_limits<double>::infinity();
-
-    for ( size_t index = 0; index < aggregatedResults.size(); ++index )
-    {
-        if ( aggregatedResults[index] != std::numeric_limits<double>::infinity() )
-        {
-            minV = std::min( minV, aggregatedResults[index] );
-        }
-    }
-    return minV;
+    double minV = RiaStatisticsTools::minimumValue( aggregatedResults );
+    return minV != std::numeric_limits<double>::max() ? minV : std::numeric_limits<double>::infinity();
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ReservoirDataModel/RigEclipseContourMapProjection.cpp

@@ -129,6 +129,7 @@ std::pair<bool, std::vector<double>>
             // TODO: this was RimEclipseCellColors->hasStaticResult()
             if ( resultAddress.resultCatType() == RiaDefines::ResultCatType::STATIC_NATIVE && timeStep > 0 ) timeStep = 0;
 
+            resultData.ensureKnownResultLoaded( resultAddress );
             // When loading a project file, grid calculator results are not computed the first time this function is
             // called. Must check if result is loaded. See RimReloadCaseTools::updateAll3dViews()
             if ( resultAddress.isValid() && resultData.hasResultEntry( resultAddress ) && resultData.isResultLoaded( resultAddress ) )