OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

#11708 Summary: Improve robustness for statistics calculations

Browse Source 
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.
This commit is contained in:
Magne Sjaastad 2024-09-19 09:49:56 +02:00
parent 48070f6539
commit 2c98438528
3 changed files with 44 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleCrossPlotStatisticsCase.cpp
View File

@ -130,7 +130,7 @@ void RimEnsembleCrossPlotStatisticsCase::calculate( const std::vector<RimSummary
std::vector<SampleData> 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 )

121
ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.cpp
View File

@ -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 <limits>
#include <set>
#include <vector>
//--------------------------------------------------------------------------------------------------
///
@ -129,55 +125,62 @@ RifSummaryReaderInterface* RimEnsembleStatisticsCase::summaryReader()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimEnsembleStatisticsCase::calculate( const std::vector<RimSummaryCase*>& sumCases,
void RimEnsembleStatisticsCase::calculate( const std::vector<RimSummaryCase*>& 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<time_t> allTimeSteps;
std::vector<std::vector<double>> 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<time_t>& 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<std::vector<double>> 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<double> valuesAtTimeStep;
valuesAtTimeStep.reserve( sumCases.size() );
valuesAtTimeStep.reserve( curveMerger.curveCount() );
for ( const std::vector<double>& 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<RimSummaryCase*> RimEnsembleStatisticsCase::validSummaryCases( const std::vector<RimSummaryCase*>& allSumCases,
const RifEclipseSummaryAddress& inputAddress,
bool includeIncompleteCurves )
std::pair<time_t, time_t> RimEnsembleStatisticsCase::findMinMaxTime( const std::vector<RimSummaryCase*>& sumCases,
const RifEclipseSummaryAddress& inputAddress )
{
std::vector<RimSummaryCase*> validCases;
std::vector<std::pair<RimSummaryCase*, time_t>> times;
time_t maxTimeStep = 0;
for ( auto& sumCase : allSumCases )
{
const auto& reader = sumCase->summaryReader();
if ( reader )
{
const std::vector<time_t>& 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<time_t, time_t> RimEnsembleStatisticsCase::findMinMaxTimeStep( const std::vector<RimSummaryCase*>& sumCases,
const RifEclipseSummaryAddress& inputAddress )
{
time_t minTimeStep = std::numeric_limits<time_t>::max();
time_t maxTimeStep = 0;
time_t minTime = std::numeric_limits<time_t>::max();
time_t maxTime = 0;
for ( const auto& sumCase : sumCases )
{
@ -278,13 +235,13 @@ std::pair<time_t, time_t> RimEnsembleStatisticsCase::findMinMaxTimeStep( const s
const std::vector<time_t>& 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 );
}
//--------------------------------------------------------------------------------------------------

11
ApplicationLibCode/ProjectDataModel/Summary/RimEnsembleStatisticsCase.h
View File

@ -43,18 +43,15 @@ public:
RifSummaryReaderInterface* summaryReader() override;
RiaDefines::EclipseUnitSystem unitSystem() const override;
void calculate( const std::vector<RimSummaryCase*>& sumCases, const RifEclipseSummaryAddress& inputAddress, bool includeIncompleteCurves );
void calculate( const std::vector<RimSummaryCase*>& summaryCases, const RifEclipseSummaryAddress& inputAddress, bool includeIncompleteCurves );
std::vector<time_t> timeSteps( const RifEclipseSummaryAddress& resultAddress ) const override;
std::pair<bool, std::vector<double>> values( const RifEclipseSummaryAddress& resultAddress ) const override;
std::string unitName( const RifEclipseSummaryAddress& resultAddress ) const override;
static std::vector<RimSummaryCase*> validSummaryCases( const std::vector<RimSummaryCase*>& allSumCases,
const RifEclipseSummaryAddress& inputAddress,
bool includeIncompleteCurves );
static std::pair<time_t, time_t> findMinMaxTimeStep( const std::vector<RimSummaryCase*>& sumCases,
const RifEclipseSummaryAddress& inputAddress );
static RiaDefines::DateTimePeriod findBestResamplingPeriod( time_t minTimeStep, time_t maxTimeStep );
static std::pair<time_t, time_t> findMinMaxTime( const std::vector<RimSummaryCase*>& sumCases,
const RifEclipseSummaryAddress& inputAddress );
static RiaDefines::DateTimePeriod findBestResamplingPeriod( time_t minTimeStep, time_t maxTimeStep );
private:
void clearData();