#14565 Grid model results: Make SGAS for Two-phase Gas/Water models when missing in output

Simulators for two-phase gas/water models may report only SWAT. In this case no SGAS result entry was created, so SGAS was missing from the result list and TERNARY fell back to computing SOIL = 1 - SWAT, presenting the gas saturation as oil saturation.

Add RigSgasResultCalculator that creates a placeholder SGAS entry when the model has gas and water but no oil phase, SWAT is present and SGAS is missing. SGAS is then computed as SGAS = 1 - SWAT. The computation is dispatched from findOrLoadKnownScalarResult() and findOrLoadKnownScalarResultForTimeStep() the same way as SOIL.

The equivalent logic in RigCaseCellResultsData::testAndComputeSgasForTimeStep() was only reachable from RigSoilResultCalculator, which never runs for models without an oil phase. This function is replaced by computeSgasForTimeStep() delegating to the new calculator, and the unreachable call in RigSoilResultCalculator is removed.
This commit is contained in:
Magne Sjaastad committed 2026-08-21 14:36:22 +02:00
1 parent 51a62edec9
commit f13fbbbd87
8 files changed
+367 -75

No files matched your search

@@ -2,6 +2,7 @@ set(SOURCE_GROUP_SOURCE_FILES
${CMAKE_CURRENT_LIST_DIR}/RigEclipseResultCalculator.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSoilResultCalculator.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSwatResultCalculator.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSgasResultCalculator.cpp
${CMAKE_CURRENT_LIST_DIR}/RigFaultDistanceCalculator.cpp
${CMAKE_CURRENT_LIST_DIR}/RigFaultDistanceResultCalculator.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSelectedFaultDistanceResultCalculator.cpp
@@ -0,0 +1,129 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2026- Equinor 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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "RigSgasResultCalculator.h"
#include "RigCaseCellResultsData.h"
#include "RigEclipseCaseData.h"
#include "RigEclipseResultInfo.h"
#include "RiaPhaseTools.h"
#include "RiaResultNames.h"
//==================================================================================================
///
//==================================================================================================
RigSgasResultCalculator::RigSgasResultCalculator( RigCaseCellResultsData& resultsData )
: RigEclipseResultCalculator( resultsData )
{
}
//==================================================================================================
///
//==================================================================================================
RigSgasResultCalculator::~RigSgasResultCalculator()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigSgasResultCalculator::checkAndCreatePlaceholderEntry( const RigEclipseResultAddress& resVarAddr )
{
if ( !isMatching( resVarAddr ) ) return;
// Only compute SGAS for two-phase gas/water models, see issue #14565
if ( !isTwoPhaseGasWater() ) return;
if ( m_resultsData->hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ) ) )
{
return;
}
if ( !m_resultsData->hasResultEntry( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::swat() ) ) )
{
return;
}
bool needsToBeStored = false;
size_t sgasIndex = m_resultsData->findOrCreateScalarResultIndex( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE,
RiaResultNames::sgas() ),
needsToBeStored );
m_resultsData->setMustBeCalculated( sgasIndex );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigSgasResultCalculator::isMatching( const RigEclipseResultAddress& resVarAddr ) const
{
return resVarAddr.resultName() == RiaResultNames::sgas() && resVarAddr.resultCatType() == RiaDefines::ResultCatType::DYNAMIC_NATIVE;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigSgasResultCalculator::calculate( const RigEclipseResultAddress& resVarAddr, size_t timeStepIndex )
{
if ( !isTwoPhaseGasWater() ) return;
RigEclipseResultAddress swatAddr( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::swat() );
size_t scalarIndexSWAT = m_resultsData->findOrLoadKnownScalarResultForTimeStep( swatAddr, timeStepIndex );
if ( scalarIndexSWAT == cvf::UNDEFINED_SIZE_T ) return;
size_t swatTimeStepCount = m_resultsData->infoForEachResultIndex()[scalarIndexSWAT].timeStepInfos().size();
if ( timeStepIndex >= swatTimeStepCount ) return;
size_t swatResultValueCount = m_resultsData->cellScalarResults( swatAddr, timeStepIndex ).size();
if ( swatResultValueCount == 0 ) return;
// Make sure memory is allocated for the new SGAS results
size_t sgasResultScalarIndex = m_resultsData->findScalarResultIndexFromAddress( resVarAddr );
if ( sgasResultScalarIndex == cvf::UNDEFINED_SIZE_T ) return;
m_resultsData->m_cellScalarResults[sgasResultScalarIndex].resize( swatTimeStepCount );
if ( !m_resultsData->cellScalarResults( resVarAddr, timeStepIndex ).empty() )
{
// Data is computed and allocated, nothing more to do
return;
}
m_resultsData->m_cellScalarResults[sgasResultScalarIndex][timeStepIndex].resize( swatResultValueCount );
const std::vector<double>& swatForTimeStep = m_resultsData->cellScalarResults( swatAddr, timeStepIndex );
std::vector<double>* sgasForTimeStep = m_resultsData->modifiableCellScalarResult( resVarAddr, timeStepIndex );
if ( !sgasForTimeStep ) return;
#pragma omp parallel for
for ( int idx = 0; idx < static_cast<int>( swatResultValueCount ); idx++ )
{
sgasForTimeStep->at( idx ) = 1.0 - swatForTimeStep[idx];
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigSgasResultCalculator::isTwoPhaseGasWater() const
{
const RigEclipseCaseData* eclipseCaseData = m_resultsData->m_ownerCaseData;
if ( !eclipseCaseData ) return false;
return RiaPhaseTools::isTwoPhaseGasWater( eclipseCaseData->availablePhases() );
}
@@ -0,0 +1,44 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2026- Equinor 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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include <cstddef>
#include "RigEclipseResultCalculator.h"
class RigCaseCellResultsData;
class RigEclipseResultAddress;
//==================================================================================================
/// Compute SGAS for two-phase gas/water models where SGAS is not present in the simulator output.
/// SGAS is then computed as SGAS = 1 - SWAT
//==================================================================================================
class RigSgasResultCalculator : public RigEclipseResultCalculator
{
public:
RigSgasResultCalculator( RigCaseCellResultsData& resultsData );
~RigSgasResultCalculator() override;
void checkAndCreatePlaceholderEntry( const RigEclipseResultAddress& resVarAddr ) override;
bool isMatching( const RigEclipseResultAddress& resVarAddr ) const override;
void calculate( const RigEclipseResultAddress& resVarAddr, size_t timeStepIndex ) override;
private:
bool isTwoPhaseGasWater() const;
};
@@ -77,10 +77,6 @@ bool RigSoilResultCalculator::isMatching( const RigEclipseResultAddress& resVarA
//--------------------------------------------------------------------------------------------------
void RigSoilResultCalculator::calculate( const RigEclipseResultAddress& resVarAddr, size_t timeStepIndex )
{
// See similar function in RifReaderOpmRft::values, but the current implementation is not suitable for merging
// Compute SGAS based on SWAT if the simulation contains no oil
m_resultsData->testAndComputeSgasForTimeStep( timeStepIndex );
RigEclipseResultAddress SWATAddr( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::swat() );
RigEclipseResultAddress SGASAddr( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() );
RigEclipseResultAddress SSOLAddr( RiaDefines::ResultCatType::DYNAMIC_NATIVE, "SSOL" );
@@ -49,6 +49,7 @@
#include "RigNestedHybridGridResultTools.h"
#include "RigOilVolumeResultCalculator.h"
#include "RigPorvSoilSgasResultCalculator.h"
#include "RigSgasResultCalculator.h"
#include "RigSoilResultCalculator.h"
#include "RigStatisticsDataCache.h"
#include "RigStatisticsMath.h"
@@ -1041,6 +1042,13 @@ void RigCaseCellResultsData::createPlaceholderResultEntries()
RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::swat() ) );
}
// SGAS
{
RigSgasResultCalculator sgasCalculator( *this );
sgasCalculator.checkAndCreatePlaceholderEntry(
RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ) );
}
// Oil Volume
if ( RiaPreferencesSystem::current()->isFeatureEnabled( "oil-volume-result" ) )
{
@@ -1547,6 +1555,25 @@ size_t RigCaseCellResultsData::findOrLoadKnownScalarResult( const RigEclipseResu
return scalarResultIndex;
}
}
else if ( resultName == RiaResultNames::sgas() )
{
if ( mustBeCalculated( scalarResultIndex ) )
{
// Trigger loading of SWAT to establish time step count if no data has been loaded from file at this point
findOrLoadKnownScalarResult( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::swat() ) );
m_cellScalarResults[scalarResultIndex].resize( maxTimeStepCount() );
for ( size_t timeStepIdx = 0; timeStepIdx < maxTimeStepCount(); timeStepIdx++ )
{
if ( m_cellScalarResults[scalarResultIndex][timeStepIdx].empty() )
{
computeSgasForTimeStep( timeStepIdx );
}
}
return scalarResultIndex;
}
}
else if ( resultName == RiaResultNames::completionTypeResultName() )
{
if ( type == RiaDefines::ResultCatType::STATIC_NATIVE )
@@ -1749,6 +1776,22 @@ size_t RigCaseCellResultsData::findOrLoadKnownScalarResultForTimeStep( const Rig
return soilScalarResultIndex;
}
}
else if ( type == RiaDefines::ResultCatType::DYNAMIC_NATIVE && resultName.toUpper() == RiaResultNames::sgas() )
{
size_t sgasScalarResultIndex = findScalarResultIndexFromAddress( resVarAddr );
if ( mustBeCalculated( sgasScalarResultIndex ) )
{
m_cellScalarResults[sgasScalarResultIndex].resize( maxTimeStepCount() );
if ( m_cellScalarResults[sgasScalarResultIndex][timeStepIndex].empty() )
{
computeSgasForTimeStep( timeStepIndex );
}
return sgasScalarResultIndex;
}
}
else if ( type == RiaDefines::ResultCatType::DYNAMIC_NATIVE && resultName == RiaResultNames::completionTypeResultName() )
{
size_t completionTypeScalarResultIndex = findScalarResultIndexFromAddress( resVarAddr );
@@ -1851,77 +1894,11 @@ void RigCaseCellResultsData::computeSOILForTimeStep( size_t timeStepIndex )
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseCellResultsData::testAndComputeSgasForTimeStep( size_t timeStepIndex )
void RigCaseCellResultsData::computeSgasForTimeStep( size_t timeStepIndex )
{
size_t scalarIndexSWAT =
findOrLoadKnownScalarResultForTimeStep( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::swat() ),
timeStepIndex );
if ( scalarIndexSWAT == cvf::UNDEFINED_SIZE_T )
{
return;
}
if ( m_readerInterface.isNull() ) return;
std::set<RiaDefines::PhaseType> phases = m_readerInterface->availablePhases();
if ( phases.count( RiaDefines::PhaseType::GAS_PHASE ) == 0 ) return;
if ( phases.count( RiaDefines::PhaseType::OIL_PHASE ) > 0 ) return;
// Simulation type is gas and water. No SGAS is present, compute SGAS based on SWAT
size_t scalarIndexSGAS =
findOrCreateScalarResultIndex( RigEclipseResultAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() ), false );
if ( m_cellScalarResults[scalarIndexSGAS].size() > timeStepIndex )
{
std::vector<double>& values = m_cellScalarResults[scalarIndexSGAS][timeStepIndex];
if ( !values.empty() ) return;
}
size_t swatResultValueCount = 0;
size_t swatTimeStepCount = 0;
{
std::vector<double>& swatForTimeStep = m_cellScalarResults[scalarIndexSWAT][timeStepIndex];
if ( !swatForTimeStep.empty() )
{
swatResultValueCount = swatForTimeStep.size();
swatTimeStepCount = infoForEachResultIndex()[scalarIndexSWAT].timeStepInfos().size();
}
}
m_cellScalarResults[scalarIndexSGAS].resize( swatTimeStepCount );
if ( !m_cellScalarResults[scalarIndexSGAS][timeStepIndex].empty() )
{
return;
}
m_cellScalarResults[scalarIndexSGAS][timeStepIndex].resize( swatResultValueCount );
std::vector<double>* swatForTimeStep = nullptr;
{
swatForTimeStep = &( m_cellScalarResults[scalarIndexSWAT][timeStepIndex] );
if ( swatForTimeStep->empty() )
{
swatForTimeStep = nullptr;
}
}
std::vector<double>& sgasForTimeStep = m_cellScalarResults[scalarIndexSGAS][timeStepIndex];
#pragma omp parallel for
for ( int idx = 0; idx < static_cast<int>( swatResultValueCount ); idx++ )
{
double sgasValue = 1.0;
if ( swatForTimeStep )
{
sgasValue -= swatForTimeStep->at( idx );
}
sgasForTimeStep[idx] = sgasValue;
}
RigEclipseResultAddress SGASAddr( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() );
RigSgasResultCalculator calculator( *this );
calculator.calculate( SGASAddr, timeStepIndex );
}
//--------------------------------------------------------------------------------------------------
@@ -199,6 +199,7 @@ private:
friend class RimEclipseStatisticsCaseEvaluator;
friend class RigSoilResultCalculator;
friend class RigSwatResultCalculator;
friend class RigSgasResultCalculator;
friend class RigFaultDistanceResultCalculator;
friend class RigMobilePoreVolumeResultCalculator;
friend class RigIndexIjkResultCalculator;
@@ -218,7 +219,7 @@ private:
void setMustBeCalculated( size_t scalarResultIndex );
void computeSOILForTimeStep( size_t timeStepIndex );
void testAndComputeSgasForTimeStep( size_t timeStepIndex );
void computeSgasForTimeStep( size_t timeStepIndex );
bool hasCompleteTransmissibilityResults() const;
@@ -25,6 +25,7 @@ set(SOURCE_UNITTEST_FILES
${CMAKE_CURRENT_LIST_DIR}/RigEclipseResultTools-Test.cpp
${CMAKE_CURRENT_LIST_DIR}/RigEnsembleFractureStatisticsCalculator-Test.cpp
${CMAKE_CURRENT_LIST_DIR}/RigMobilePoreVolumeResultCalculator-Test.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSgasResultCalculator-Test.cpp
${CMAKE_CURRENT_LIST_DIR}/RigReservoir-Test.cpp
${CMAKE_CURRENT_LIST_DIR}/RigResdataGridConverter-Test.cpp
${CMAKE_CURRENT_LIST_DIR}/RigGridExportAdapter-Test.cpp
@@ -0,0 +1,143 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2026- Equinor 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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "gtest/gtest.h"
#include "RiaDefines.h"
#include "RiaResultNames.h"
#include "RifReaderMockModel.h"
#include "RigActiveCellInfo.h"
#include "RigCaseCellResultsData.h"
#include "RigEclipseCaseData.h"
#include "RigEclipseResultAddress.h"
#include "RigEclipseResultInfo.h"
#include "RigMainGrid.h"
#include "RimEclipseResultCase.h"
#include <QDateTime>
#include <memory>
#include <set>
namespace
{
struct MockCase
{
std::unique_ptr<RimEclipseResultCase> resultCase;
cvf::ref<RigEclipseCaseData> caseData;
};
MockCase createMockCase( const std::set<RiaDefines::PhaseType>& phases )
{
MockCase mockCase;
mockCase.resultCase.reset( new RimEclipseResultCase );
mockCase.caseData = new RigEclipseCaseData( mockCase.resultCase.get() );
cvf::ref<RifReaderMockModel> mockReader = new RifReaderMockModel;
mockReader->setWorldCoordinates( cvf::Vec3d( 0, 0, 0 ), cvf::Vec3d( 100, 100, 100 ) );
mockReader->setCellCounts( cvf::Vec3st( 2, 2, 2 ) );
mockReader->setResultInfo( 0, 0 );
mockReader->enableWellData( false );
mockReader->open( "", mockCase.caseData.p() );
mockCase.caseData->mainGrid()->computeCachedData();
mockCase.caseData->setAvailablePhases( phases );
mockCase.resultCase->setReservoirData( mockCase.caseData.p() );
return mockCase;
}
void addSwatResult( RigCaseCellResultsData* cellResults, const std::vector<std::vector<double>>& valuesForEachTimeStep )
{
RigEclipseResultAddress swatAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::swat() );
cellResults->createResultEntry( swatAddress, false );
std::vector<RigEclipseTimeStepInfo> timeStepInfos;
for ( int i = 0; i < static_cast<int>( valuesForEachTimeStep.size() ); i++ )
{
timeStepInfos.push_back( RigEclipseTimeStepInfo( QDateTime(), i, i ) );
}
cellResults->setTimeStepInfos( swatAddress, timeStepInfos );
auto* timeStepValues = cellResults->modifiableCellScalarResultTimesteps( swatAddress );
timeStepValues->resize( valuesForEachTimeStep.size() );
for ( size_t i = 0; i < valuesForEachTimeStep.size(); i++ )
{
( *timeStepValues )[i] = valuesForEachTimeStep[i];
}
}
} // namespace
//--------------------------------------------------------------------------------------------------
/// SGAS is not reported by the simulator for two-phase gas/water models, compute SGAS = 1 - SWAT
//--------------------------------------------------------------------------------------------------
TEST( RigSgasResultCalculatorTest, ComputeSgasForTwoPhaseGasWater )
{
MockCase mockCase = createMockCase( { RiaDefines::PhaseType::GAS_PHASE, RiaDefines::PhaseType::WATER_PHASE } );
RigCaseCellResultsData* cellResults = mockCase.caseData->results( RiaDefines::PorosityModelType::MATRIX_MODEL );
ASSERT_NE( cellResults, nullptr );
const size_t cellCount = cellResults->activeCellInfo()->reservoirActiveCellCount();
ASSERT_GT( cellCount, 0u );
const std::vector<std::vector<double>> swatValues = { std::vector<double>( cellCount, 0.2 ), std::vector<double>( cellCount, 0.7 ) };
addSwatResult( cellResults, swatValues );
cellResults->createPlaceholderResultEntries();
const RigEclipseResultAddress sgasAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() );
ASSERT_TRUE( cellResults->hasResultEntry( sgasAddress ) );
ASSERT_TRUE( cellResults->ensureKnownResultLoaded( sgasAddress ) );
for ( size_t timeStepIndex = 0; timeStepIndex < swatValues.size(); timeStepIndex++ )
{
const auto& sgasValues = cellResults->cellScalarResults( sgasAddress, timeStepIndex );
ASSERT_EQ( sgasValues.size(), cellCount );
for ( size_t i = 0; i < cellCount; i++ )
{
EXPECT_DOUBLE_EQ( sgasValues[i], 1.0 - swatValues[timeStepIndex][i] );
}
}
}
//--------------------------------------------------------------------------------------------------
/// SGAS must not be derived from SWAT when the oil phase is present
//--------------------------------------------------------------------------------------------------
TEST( RigSgasResultCalculatorTest, NoComputedSgasWhenOilPhaseIsPresent )
{
MockCase mockCase =
createMockCase( { RiaDefines::PhaseType::OIL_PHASE, RiaDefines::PhaseType::GAS_PHASE, RiaDefines::PhaseType::WATER_PHASE } );
RigCaseCellResultsData* cellResults = mockCase.caseData->results( RiaDefines::PorosityModelType::MATRIX_MODEL );
ASSERT_NE( cellResults, nullptr );
const size_t cellCount = cellResults->activeCellInfo()->reservoirActiveCellCount();
ASSERT_GT( cellCount, 0u );
addSwatResult( cellResults, { std::vector<double>( cellCount, 0.2 ) } );
cellResults->createPlaceholderResultEntries();
const RigEclipseResultAddress sgasAddress( RiaDefines::ResultCatType::DYNAMIC_NATIVE, RiaResultNames::sgas() );
EXPECT_FALSE( cellResults->hasResultEntry( sgasAddress ) );
}