From 2c98438528d5098bc412fcd300488345be15f927 Mon Sep 17 00:00:00 2001 From: Magne Sjaastad Date: Thu, 19 Sep 2024 09:49:56 +0200 Subject: [PATCH] #11708 Summary: Improve robustness for statistics calculations There can be incomplete realizations in an ensemble. Make sure that the calculation is robust for all variants of incomplete curves. Use RiaTimeHistoryCurveMerger to ensure one definition of time steps and extraction of curve values to be used to compute ensemble statistics. Remove unused function. --- .../RimEnsembleCrossPlotStatisticsCase.cpp | 2 +- .../Summary/RimEnsembleStatisticsCase.cpp | 121 ++++++------------ .../Summary/RimEnsembleStatisticsCase.h | 11 +- 3 files changed, 44 insertions(+), 90 deletions(-) diff --git a/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleCrossPlotStatisticsCase.cpp b/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleCrossPlotStatisticsCase.cpp index 6ffb4259f2..fdeb533540 100644 --- a/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleCrossPlotStatisticsCase.cpp +++ b/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleCrossPlotStatisticsCase.cpp @@ -130,7 +130,7 @@ void RimEnsembleCrossPlotStatisticsCase::calculate( const std::vector sampleData; - auto [minTimeStep, maxTimeStep] = RimEnsembleStatisticsCase::findMinMaxTimeStep( sumCases, inputAddressX ); + auto [minTimeStep, maxTimeStep] = RimEnsembleStatisticsCase::findMinMaxTime( sumCases, inputAddressX ); RiaDefines::DateTimePeriod period = RimEnsembleStatisticsCase::findBestResamplingPeriod( minTimeStep, maxTimeStep ); for ( const auto& sumCase : sumCases ) diff --git a/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.cpp b/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.cpp index 82e3be9ebb..ee8416675a 100644 --- a/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.cpp +++ b/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.cpp @@ -18,17 +18,13 @@ #include "RimEnsembleStatisticsCase.h" +#include "RiaCurveMerger.h" #include "RiaSummaryTools.h" #include "RiaTimeHistoryCurveResampler.h" #include "RigStatisticsMath.h" -#include "RimEnsembleCurveSet.h" -#include "RimSummaryEnsemble.h" - #include -#include -#include //-------------------------------------------------------------------------------------------------- /// @@ -129,55 +125,62 @@ RifSummaryReaderInterface* RimEnsembleStatisticsCase::summaryReader() //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RimEnsembleStatisticsCase::calculate( const std::vector& sumCases, +void RimEnsembleStatisticsCase::calculate( const std::vector& summaryCases, const RifEclipseSummaryAddress& inputAddress, bool includeIncompleteCurves ) { clearData(); + if ( !inputAddress.isValid() ) return; - if ( sumCases.empty() ) return; + if ( summaryCases.empty() ) return; // Use first summary case to get unit system and other meta data - m_firstSummaryCase = sumCases.front(); + m_firstSummaryCase = summaryCases.front(); - auto [minTimeStep, maxTimeStep] = findMinMaxTimeStep( sumCases, inputAddress ); - RiaDefines::DateTimePeriod period = findBestResamplingPeriod( minTimeStep, maxTimeStep ); + const auto [minTime, maxTime] = findMinMaxTime( summaryCases, inputAddress ); + RiaDefines::DateTimePeriod period = findBestResamplingPeriod( minTime, maxTime ); - std::vector allTimeSteps; - std::vector> caseAndTimeStepValues; - caseAndTimeStepValues.reserve( sumCases.size() ); - for ( const auto& sumCase : sumCases ) + // The last time step for the individual realizations in an ensemble is usually identical. Add a small threshold to improve robustness. + const auto timeThreshold = maxTime - ( maxTime - minTime ) * 0.01; + + RiaTimeHistoryCurveMerger curveMerger; + + for ( const auto& sumCase : summaryCases ) { const auto& reader = sumCase->summaryReader(); if ( reader ) { const std::vector& timeSteps = reader->timeSteps( inputAddress ); - auto [isOk, values] = reader->values( inputAddress ); + const auto [isOk, values] = reader->values( inputAddress ); - if ( values.empty() ) continue; + if ( values.empty() || timeSteps.empty() ) continue; - if ( !includeIncompleteCurves && timeSteps.size() != values.size() ) continue; + if ( !includeIncompleteCurves && ( timeSteps.back() < timeThreshold ) ) continue; - auto [resampledTimeSteps, resampledValues] = RiaSummaryTools::resampledValuesForPeriod( inputAddress, timeSteps, values, period ); - - if ( allTimeSteps.empty() ) allTimeSteps = resampledTimeSteps; - caseAndTimeStepValues.push_back( resampledValues ); + const auto [resampledTimeSteps, resampledValues] = + RiaSummaryTools::resampledValuesForPeriod( inputAddress, timeSteps, values, period ); + curveMerger.addCurveData( resampledTimeSteps, resampledValues ); } } - m_timeSteps = allTimeSteps; + curveMerger.computeInterpolatedValues(); - for ( size_t timeStepIndex = 0; timeStepIndex < allTimeSteps.size(); timeStepIndex++ ) + std::vector> curveValues; + for ( size_t i = 0; i < curveMerger.curveCount(); i++ ) + { + curveValues.push_back( curveMerger.interpolatedYValuesForAllXValues( i ) ); + } + + m_timeSteps = curveMerger.allXValues(); + + for ( size_t timeStepIndex = 0; timeStepIndex < m_timeSteps.size(); timeStepIndex++ ) { std::vector valuesAtTimeStep; - valuesAtTimeStep.reserve( sumCases.size() ); + valuesAtTimeStep.reserve( curveMerger.curveCount() ); - for ( const std::vector& caseValues : caseAndTimeStepValues ) + for ( size_t curveIdx = 0; curveIdx < curveMerger.curveCount(); curveIdx++ ) { - if ( timeStepIndex < caseValues.size() ) - { - valuesAtTimeStep.push_back( caseValues[timeStepIndex] ); - } + valuesAtTimeStep.push_back( curveValues[curveIdx][timeStepIndex] ); } double p10, p50, p90, mean; @@ -218,57 +221,11 @@ void RimEnsembleStatisticsCase::clearData() //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::vector RimEnsembleStatisticsCase::validSummaryCases( const std::vector& allSumCases, - const RifEclipseSummaryAddress& inputAddress, - bool includeIncompleteCurves ) +std::pair RimEnsembleStatisticsCase::findMinMaxTime( const std::vector& sumCases, + const RifEclipseSummaryAddress& inputAddress ) { - std::vector validCases; - std::vector> times; - - time_t maxTimeStep = 0; - - for ( auto& sumCase : allSumCases ) - { - const auto& reader = sumCase->summaryReader(); - if ( reader ) - { - const std::vector& timeSteps = reader->timeSteps( inputAddress ); - if ( !timeSteps.empty() ) - { - time_t lastTimeStep = timeSteps.back(); - - maxTimeStep = std::max( lastTimeStep, maxTimeStep ); - times.push_back( std::make_pair( sumCase, lastTimeStep ) ); - } - } - } - - for ( const auto& [sumCase, lastTimeStep] : times ) - { - // Previous versions tested on identical first time step, this test is now removed. For large simulations with - // numerical issues the first time step can be slightly different - // - // https://github.com/OPM/ResInsight/issues/7774 - // - if ( !includeIncompleteCurves && lastTimeStep != maxTimeStep ) - { - continue; - } - - validCases.push_back( sumCase ); - } - - return validCases; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::pair RimEnsembleStatisticsCase::findMinMaxTimeStep( const std::vector& sumCases, - const RifEclipseSummaryAddress& inputAddress ) -{ - time_t minTimeStep = std::numeric_limits::max(); - time_t maxTimeStep = 0; + time_t minTime = std::numeric_limits::max(); + time_t maxTime = 0; for ( const auto& sumCase : sumCases ) { @@ -278,13 +235,13 @@ std::pair RimEnsembleStatisticsCase::findMinMaxTimeStep( const s const std::vector& timeSteps = reader->timeSteps( inputAddress ); if ( !timeSteps.empty() ) { - minTimeStep = std::min( timeSteps.front(), minTimeStep ); - maxTimeStep = std::max( timeSteps.back(), maxTimeStep ); + minTime = std::min( timeSteps.front(), minTime ); + maxTime = std::max( timeSteps.back(), maxTime ); } } } - return std::make_pair( minTimeStep, maxTimeStep ); + return std::make_pair( minTime, maxTime ); } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.h b/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.h index 7322f30e67..9b3fcc518f 100644 --- a/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.h +++ b/ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.h @@ -43,18 +43,15 @@ public: RifSummaryReaderInterface* summaryReader() override; RiaDefines::EclipseUnitSystem unitSystem() const override; - void calculate( const std::vector& sumCases, const RifEclipseSummaryAddress& inputAddress, bool includeIncompleteCurves ); + void calculate( const std::vector& summaryCases, const RifEclipseSummaryAddress& inputAddress, bool includeIncompleteCurves ); std::vector timeSteps( const RifEclipseSummaryAddress& resultAddress ) const override; std::pair> values( const RifEclipseSummaryAddress& resultAddress ) const override; std::string unitName( const RifEclipseSummaryAddress& resultAddress ) const override; - static std::vector validSummaryCases( const std::vector& allSumCases, - const RifEclipseSummaryAddress& inputAddress, - bool includeIncompleteCurves ); - static std::pair findMinMaxTimeStep( const std::vector& sumCases, - const RifEclipseSummaryAddress& inputAddress ); - static RiaDefines::DateTimePeriod findBestResamplingPeriod( time_t minTimeStep, time_t maxTimeStep ); + static std::pair findMinMaxTime( const std::vector& sumCases, + const RifEclipseSummaryAddress& inputAddress ); + static RiaDefines::DateTimePeriod findBestResamplingPeriod( time_t minTimeStep, time_t maxTimeStep ); private: void clearData();