OilfieldToolsNet/ResInsight
4
0 Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2026-10-25 06:38:13 -05:00
Code Issues Packages Projects Releases Wiki Activity

Split NNC, aquifer and connection code out of RigSimulationInputTool (#15001)

Browse Source 
Move EDITNNC, AQU* and COMPDAT/COMPSEGS processing into RigSimulationInput{Nnc,Aquifer,Connection}Tools, and move their tests into matching test files.
This commit is contained in:
Kristian Bendiksen committed 2026-10-23 12:32:17 +02:00
1 parent 7dfa813c06
commit 9feace91ee
17 files changed
+2407 -1984
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

ApplicationLibCode/ReservoirDataModel/SimulationFile/CMakeLists_files.cmake
+3
View File
@@ -1,6 +1,9 @@
set(SOURCE_GROUP_SOURCE_FILES
${CMAKE_CURRENT_LIST_DIR}/RigModelPaddingSettings.cpp
${CMAKE_CURRENT_LIST_DIR}/RigPadModel.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputAquiferTools.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputConnectionTools.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputNncTools.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputSettings.cpp
${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputTool.cpp
)
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.cpp
+279
View File
@@ -0,0 +1,279 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2025- 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 "RigSimulationInputAquiferTools.h"
#include "RifOpmDeckTools.h"
#include "RifOpmFlowDeckFile.h"
#include "RimEclipseCase.h"
#include "RigBoundingBoxIjk.h"
#include "RigSimulationInputSettings.h"
#include "RigSimulationInputTool.h"
#include "opm/input/eclipse/Deck/DeckItem.hpp"
#include "opm/input/eclipse/Deck/DeckRecord.hpp"
#include "opm/input/eclipse/Parser/ParserKeywords/A.hpp"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<void, QString> RigSimulationInputAquiferTools::replaceAquconKeywordIndices( RimEclipseCase* eclipseCase,
const RigSimulationInputSettings& settings,
RifOpmFlowDeckFile& deckFile )
{
return RigSimulationInputTool::replaceKeywordWithBoxIndices( "AQUCON", eclipseCase, settings, deckFile, processAquconRecord );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<void, QString> RigSimulationInputAquiferTools::replaceAquanconKeywordIndices( RimEclipseCase* eclipseCase,
const RigSimulationInputSettings& settings,
RifOpmFlowDeckFile& deckFile )
{
return RigSimulationInputTool::replaceKeywordWithBoxIndices( "AQUANCON", eclipseCase, settings, deckFile, processAquanconRecord );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<void, QString> RigSimulationInputAquiferTools::replaceAqunumKeywordIndices( RimEclipseCase* eclipseCase,
const RigSimulationInputSettings& settings,
RifOpmFlowDeckFile& deckFile )
{
return RigSimulationInputTool::replaceKeywordWithBoxIndices( "AQUNUM", eclipseCase, settings, deckFile, processAqunumRecord );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<Opm::DeckRecord, QString> RigSimulationInputAquiferTools::processAquconRecord( const Opm::DeckRecord& record,
const caf::VecIjk0& min,
const caf::VecIjk0& max,
const RigRefinement& refinement )
{
// AQUCON format: ID I1 I2 J1 J2 K1 K2 CONNECT_FACE ...
// Items: 0=ID, 1=I1, 2=I2, 3=J1, 4=J2, 5=K1, 6=K2, 7+=other parameters
if ( record.size() < 1 )
{
return std::unexpected( QString( "AQUCON record has insufficient items (expected at least 1, got %1)" ).arg( record.size() ) );
}
using A = Opm::ParserKeywords::AQUCON;
std::vector<Opm::DeckItem> items;
// Copy aquifer ID (first item)
items.push_back( record.getItem( 0 ) );
if ( record.size() >= 7 && record.getItem( 1 ).hasValue( 0 ) && record.getItem( 2 ).hasValue( 0 ) && record.getItem( 3 ).hasValue( 0 ) &&
record.getItem( 4 ).hasValue( 0 ) && record.getItem( 5 ).hasValue( 0 ) && record.getItem( 6 ).hasValue( 0 ) )
{
// Transform IJK box coordinates (items 1-6: I1, I2, J1, J2, K1, K2)
// Note: AQUCON uses 1-based Eclipse coordinates
caf::VecIjk1 origMin = RigSimulationInputTool::extractIjk( record, 1, 3, 5 ); // I1, J1, K1
caf::VecIjk1 origMax = RigSimulationInputTool::extractIjk( record, 2, 4, 6 ); // I2, J2, K2
// Create input bounding box (0-based, inclusive)
RigBoundingBoxIjk<caf::VecIjk0> inputBox( origMin.toZeroBased(), origMax.toZeroBased() );
// Get aquifer ID for logging
int aquiferId = record.getItem( 0 ).get<int>( 0 );
QString recordIdent = QString( "for aquifer %1" ).arg( aquiferId );
// Transform box to sector coordinates
auto transformResult = RigSimulationInputTool::transformBoxToSectorCoordinates( inputBox,
min,
max,
refinement,
QString::fromStdString( A::keywordName ),
recordIdent );
if ( !transformResult )
{
return std::unexpected( transformResult.error() );
}
items.push_back( RifOpmDeckTools::item( A::I1::itemName, static_cast<int>( transformResult->min().x() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::I2::itemName, static_cast<int>( transformResult->max().x() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::J1::itemName, static_cast<int>( transformResult->min().y() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::J2::itemName, static_cast<int>( transformResult->max().y() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::K1::itemName, static_cast<int>( transformResult->min().z() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::K2::itemName, static_cast<int>( transformResult->max().z() + 1 ) ) );
// Copy remaining items (CONNECT_FACE and other parameters)
for ( size_t i = 7; i < record.size(); ++i )
{
items.push_back( record.getItem( i ) );
}
}
else
{
// No box definition or incomplete box definition - copy remaining items as-is
for ( size_t i = 1; i < record.size(); ++i )
{
items.push_back( record.getItem( i ) );
}
}
return Opm::DeckRecord{ std::move( items ) };
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<Opm::DeckRecord, QString> RigSimulationInputAquiferTools::processAquanconRecord( const Opm::DeckRecord& record,
const caf::VecIjk0& min,
const caf::VecIjk0& max,
const RigRefinement& refinement )
{
// AQUANCON format: ID I1 I2 J1 J2 K1 K2 AQU_FACE ...
// Items: 0=ID, 1=I1, 2=I2, 3=J1, 4=J2, 5=K1, 6=K2, 7+=other parameters
if ( record.size() < 1 )
{
return std::unexpected( QString( "AQUANCON record has insufficient items (expected at least 1, got %1)" ).arg( record.size() ) );
}
using A = Opm::ParserKeywords::AQUANCON;
std::vector<Opm::DeckItem> items;
// Copy aquifer ID (first item)
items.push_back( record.getItem( 0 ) );
if ( record.size() >= 7 && record.getItem( 1 ).hasValue( 0 ) && record.getItem( 2 ).hasValue( 0 ) && record.getItem( 3 ).hasValue( 0 ) &&
record.getItem( 4 ).hasValue( 0 ) && record.getItem( 5 ).hasValue( 0 ) && record.getItem( 6 ).hasValue( 0 ) )
{
// Transform IJK box coordinates (items 1-6: I1, I2, J1, J2, K1, K2)
// Note: AQUANCON uses 1-based Eclipse coordinates
caf::VecIjk1 origMin = RigSimulationInputTool::extractIjk( record, 1, 3, 5 ); // I1, J1, K1
caf::VecIjk1 origMax = RigSimulationInputTool::extractIjk( record, 2, 4, 6 ); // I2, J2, K2
// Create input bounding box (0-based, inclusive)
RigBoundingBoxIjk<caf::VecIjk0> inputBox( origMin.toZeroBased(), origMax.toZeroBased() );
// Get aquifer ID for logging
int aquiferId = record.getItem( 0 ).get<int>( 0 );
QString recordIdent = QString( "for aquifer %1" ).arg( aquiferId );
// Transform box to sector coordinates
auto transformResult = RigSimulationInputTool::transformBoxToSectorCoordinates( inputBox,
min,
max,
refinement,
QString::fromStdString( A::keywordName ),
recordIdent );
if ( !transformResult )
{
return std::unexpected( transformResult.error() );
}
items.push_back( RifOpmDeckTools::item( A::I1::itemName, static_cast<int>( transformResult->min().x() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::I2::itemName, static_cast<int>( transformResult->max().x() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::J1::itemName, static_cast<int>( transformResult->min().y() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::J2::itemName, static_cast<int>( transformResult->max().y() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::K1::itemName, static_cast<int>( transformResult->min().z() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::K2::itemName, static_cast<int>( transformResult->max().z() + 1 ) ) );
// Copy remaining items (AQUFACE and other parameters)
for ( size_t i = 7; i < record.size(); ++i )
{
items.push_back( record.getItem( i ) );
}
}
else
{
// No box definition or incomplete box definition - copy remaining items as-is
for ( size_t i = 1; i < record.size(); ++i )
{
items.push_back( record.getItem( i ) );
}
}
return Opm::DeckRecord{ std::move( items ) };
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<Opm::DeckRecord, QString> RigSimulationInputAquiferTools::processAqunumRecord( const Opm::DeckRecord& record,
const caf::VecIjk0& min,
const caf::VecIjk0& max,
const RigRefinement& refinement )
{
// AQUNUM format: ID I J K AREA ...
// Items: 0=ID, 1=I, 2=J, 3=K, 4+=other parameters
if ( record.size() < 1 )
{
return std::unexpected( QString( "AQUNUM record has insufficient items (expected at least 1, got %1)" ).arg( record.size() ) );
}
using A = Opm::ParserKeywords::AQUNUM;
std::vector<Opm::DeckItem> items;
// Copy aquifer ID (first item)
items.push_back( record.getItem( 0 ) );
if ( record.size() >= 4 && record.getItem( 1 ).hasValue( 0 ) && record.getItem( 2 ).hasValue( 0 ) && record.getItem( 3 ).hasValue( 0 ) )
{
// Transform IJK coordinates (items 1-6: I1, I2, J1, J2, K1, K2)
// Note: AQUNUM uses 1-based Eclipse coordinates
caf::VecIjk1 origCell = RigSimulationInputTool::extractIjk( record, 1, 2, 3 ); // I, J, K
// Create input bounding box (0-based, inclusive)
RigBoundingBoxIjk<caf::VecIjk0> inputBox( origCell.toZeroBased(), origCell.toZeroBased() );
// Get aquifer ID for logging
int aquiferId = record.getItem( 0 ).get<int>( 0 );
QString recordIdent = QString( "for aquifer %1" ).arg( aquiferId );
// Transform box to sector coordinates
auto transformResult = RigSimulationInputTool::transformBoxToSectorCoordinates( inputBox,
min,
max,
refinement,
QString::fromStdString( A::keywordName ),
recordIdent );
if ( !transformResult )
{
return std::unexpected( transformResult.error() );
}
items.push_back( RifOpmDeckTools::item( A::I::itemName, static_cast<int>( transformResult->min().x() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::J::itemName, static_cast<int>( transformResult->min().y() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( A::K::itemName, static_cast<int>( transformResult->min().z() + 1 ) ) );
// Copy remaining items (AREA and other parameters)
for ( size_t i = 4; i < record.size(); ++i )
{
items.push_back( record.getItem( i ) );
}
}
else
{
// No box definition or incomplete box definition - copy remaining items as-is
for ( size_t i = 1; i < record.size(); ++i )
{
items.push_back( record.getItem( i ) );
}
}
return Opm::DeckRecord{ std::move( items ) };
}
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.h
+63
View File
@@ -0,0 +1,63 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2025- 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 "RigRefinement.h"
#include "cafVecIjk.h"
#include <QString>
#include <expected>
class RimEclipseCase;
class RigSimulationInputSettings;
class RifOpmFlowDeckFile;
namespace Opm
{
class DeckRecord;
} // namespace Opm
//==================================================================================================
///
/// Aquifer keyword (AQUCON, AQUANCON, AQUNUM) handling when exporting simulation input files
///
//==================================================================================================
class RigSimulationInputAquiferTools
{
public:
static std::expected<void, QString>
replaceAquconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::expected<void, QString>
replaceAquanconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::expected<void, QString>
replaceAqunumKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::expected<Opm::DeckRecord, QString>
processAquconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement );
static std::expected<Opm::DeckRecord, QString>
processAquanconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement );
static std::expected<Opm::DeckRecord, QString>
processAqunumRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement );
};
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.cpp
+744
View File
@@ -0,0 +1,744 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2025- 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 "RigSimulationInputConnectionTools.h"
#include "RiaLogging.h"
#include "RifOpmDeckTools.h"
#include "RimEclipseCase.h"
#include "RigActiveCellInfo.h"
#include "RigCaseCellResultsData.h"
#include "RigEclipseCaseData.h"
#include "RigEclipseResultAddress.h"
#include "RigGridExportAdapter.h"
#include "RigMainGrid.h"
#include "RigSimulationInputSettings.h"
#include "RigSimulationInputTool.h"
#include "RigTransmissibilityEquations.h"
#include "opm/input/eclipse/Deck/DeckItem.hpp"
#include "opm/input/eclipse/Deck/DeckKeyword.hpp"
#include "opm/input/eclipse/Deck/DeckRecord.hpp"
#include "opm/input/eclipse/Parser/ParserKeywords/C.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <tuple>
//--------------------------------------------------------------------------------------------------
/// Refined sub-cell index and its length fraction of the original cell
//--------------------------------------------------------------------------------------------------
namespace
{
struct RefinedSubcell
{
size_t index = 0;
double fraction = 1.0;
};
} // namespace
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
static size_t centerSubcellIndex( const RigRefinement& refinement, RigRefinement::Dimension dim, size_t sectorIndex )
{
return refinement.cumulativeOffset( dim, sectorIndex ) + refinement.subcellCount( dim, sectorIndex ) / 2;
}
//--------------------------------------------------------------------------------------------------
/// All refined sub-cells of an original (sector-relative) cell index, with length fractions
//--------------------------------------------------------------------------------------------------
static std::vector<RefinedSubcell> allSubcells( const RigRefinement& refinement, RigRefinement::Dimension dim, size_t sectorIndex )
{
const size_t offset = refinement.cumulativeOffset( dim, sectorIndex );
const auto& cumulative = refinement.cumulativeFractions( dim, sectorIndex );
std::vector<RefinedSubcell> subcells;
double previous = 0.0;
for ( size_t s = 0; s < cumulative.size(); ++s )
{
subcells.push_back( { offset + s, cumulative[s] - previous } );
previous = cumulative[s];
}
return subcells;
}
//--------------------------------------------------------------------------------------------------
/// Penetration direction from COMPDAT DIR item (X -> I, Y -> J, Z or default -> K)
//--------------------------------------------------------------------------------------------------
static RigRefinement::Dimension compdatPenetrationDimension( const Opm::DeckRecord& record )
{
using C = Opm::ParserKeywords::COMPDAT;
if ( record.hasItem( C::DIR::itemName ) )
{
const auto& item = record.getItem( C::DIR::itemName );
if ( item.hasValue( 0 ) )
{
const auto dir = item.get<std::string>( 0 );
if ( !dir.empty() && ( dir[0] == 'X' || dir[0] == 'x' ) ) return RigRefinement::DimI;
if ( !dir.empty() && ( dir[0] == 'Y' || dir[0] == 'y' ) ) return RigRefinement::DimJ;
}
}
return RigRefinement::DimK;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
static bool hasExplicitDoubleValue( const Opm::DeckRecord& record, const std::string& itemName )
{
if ( !record.hasItem( itemName ) ) return false;
const auto& item = record.getItem( itemName );
return item.hasValue( 0 ) && !item.defaultApplied( 0 );
}
//--------------------------------------------------------------------------------------------------
/// Length fraction of the center sub-cell of an original (sector-relative) cell index
//--------------------------------------------------------------------------------------------------
static double centerSubcellFraction( const RigRefinement& refinement, RigRefinement::Dimension dim, size_t sectorIndex )
{
const auto& cumulative = refinement.cumulativeFractions( dim, sectorIndex );
const size_t center = refinement.subcellCount( dim, sectorIndex ) / 2;
if ( center >= cumulative.size() ) return 1.0;
return cumulative[center] - ( center > 0 ? cumulative[center - 1] : 0.0 );
}
namespace
{
struct LateralRefinementCorrection
{
double connectionFactorMultiplier = 1.0;
double pressureEquivalentRadius = 0.0;
};
} // namespace
//--------------------------------------------------------------------------------------------------
/// Correction of an explicit connection factor when the cell is refined perpendicular to the well.
/// Peaceman: CF = c * theta * K * h / ( ln( r0 / rw ) + S ). Lateral refinement keeps K * h, but reduces
/// the pressure equivalent radius r0, so CF_new = CF_old * ( ln( r0_old / rw ) + S ) / ( ln( r0_new / rw ) + S ).
/// r0_old is taken from the PR item if given, otherwise computed from the original cell.
//--------------------------------------------------------------------------------------------------
static std::optional<LateralRefinementCorrection>
lateralRefinementCorrection( const Opm::DeckRecord& record,
const std::string& wellName,
RigRefinement::Dimension penetrationDim,
double fraction1,
double fraction2,
const std::optional<RigSimulationInputConnectionTools::ConnectionCellProperties>& cellProperties )
{
using C = Opm::ParserKeywords::COMPDAT;
if ( fraction1 >= 1.0 && fraction2 >= 1.0 ) return std::nullopt;
if ( !hasExplicitDoubleValue( record, C::CONNECTION_TRANSMISSIBILITY_FACTOR::itemName ) ) return std::nullopt;
auto warning = [&]( const std::string& reason )
{
RiaLogging::warning(
std::format( "COMPDAT for well {}: connection factor not corrected for refinement perpendicular to the well ({})", wellName, reason ) );
};
if ( !hasExplicitDoubleValue( record, C::DIAMETER::itemName ) )
{
warning( "no well diameter" );
return std::nullopt;
}
const double wellRadius = 0.5 * record.getItem( C::DIAMETER::itemName ).get<double>( 0 );
double skin = 0.0;
if ( record.hasItem( C::SKIN::itemName ) && record.getItem( C::SKIN::itemName ).hasValue( 0 ) )
{
skin = record.getItem( C::SKIN::itemName ).get<double>( 0 );
}
// Peaceman radius of the original and refined cell, using the cell sizes and permeabilities perpendicular to the well
std::optional<std::pair<double, double>> computedRadii;
if ( cellProperties )
{
const int d1 = ( penetrationDim == RigRefinement::DimI ) ? 1 : 0;
const int d2 = ( penetrationDim == RigRefinement::DimK ) ? 1 : 2;
const double size1 = cellProperties->cellSize[d1];
const double size2 = cellProperties->cellSize[d2];
double perm1 = cellProperties->permeability[d1];
double perm2 = cellProperties->permeability[d2];
if ( perm1 <= 0.0 || perm2 <= 0.0 )
{
perm1 = 1.0;
perm2 = 1.0;
}
if ( size1 > 0.0 && size2 > 0.0 )
{
const double coarse = RigTransmissibilityEquations::peacemanRadius( perm1, perm2, size1, size2 );
const double fine = RigTransmissibilityEquations::peacemanRadius( perm1, perm2, fraction1 * size1, fraction2 * size2 );
computedRadii = std::make_pair( coarse, fine );
}
}
double originalRadius = 0.0;
if ( hasExplicitDoubleValue( record, C::PR::itemName ) )
{
originalRadius = record.getItem( C::PR::itemName ).get<double>( 0 );
}
else if ( computedRadii )
{
originalRadius = computedRadii->first;
}
else
{
warning( "no pressure equivalent radius or cell geometry" );
return std::nullopt;
}
// Peaceman radius is homogeneous in the cell sizes, so the ratio is exact when both directions are refined equally
const double radiusRatio = computedRadii ? computedRadii->second / computedRadii->first
: std::sqrt( 0.5 * ( fraction1 * fraction1 + fraction2 * fraction2 ) );
const double refinedRadius = originalRadius * radiusRatio;
const double originalDenominator = std::log( originalRadius / wellRadius ) + skin;
const double refinedDenominator = std::log( refinedRadius / wellRadius ) + skin;
if ( wellRadius <= 0.0 || originalDenominator <= 0.0 || refinedDenominator <= 0.0 )
{
warning(
std::format( "refined pressure equivalent radius {:.4g} is too small compared to well radius {:.4g}", refinedRadius, wellRadius ) );
return std::nullopt;
}
return LateralRefinementCorrection{ originalDenominator / refinedDenominator, refinedRadius };
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigSimulationInputConnectionTools::ConnectionCellPropertiesProvider
RigSimulationInputConnectionTools::connectionCellPropertiesProvider( RimEclipseCase* eclipseCase )
{
if ( !eclipseCase || !eclipseCase->eclipseCaseData() || !eclipseCase->eclipseCaseData()->mainGrid() ) return {};
auto* caseData = eclipseCase->eclipseCaseData();
auto* results = caseData->results( RiaDefines::PorosityModelType::MATRIX_MODEL );
auto* activeCellInfo = caseData->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
std::array<const std::vector<double>*, 3> permeabilities = { nullptr, nullptr, nullptr };
if ( results && activeCellInfo )
{
const std::array<QString, 3> permNames = { "PERMX", "PERMY", "PERMZ" };
for ( size_t i = 0; i < 3; ++i )
{
RigEclipseResultAddress address( RiaDefines::ResultCatType::STATIC_NATIVE, permNames[i] );
if ( results->ensureKnownResultLoaded( address ) && !results->cellScalarResults( address ).empty() )
{
permeabilities[i] = &results->cellScalarResults( address, 0 );
}
}
}
const RigMainGrid* mainGrid = caseData->mainGrid();
return [mainGrid, activeCellInfo, permeabilities]( const caf::VecIjk0& ijk ) -> std::optional<ConnectionCellProperties>
{
if ( !mainGrid->isCellValid( ijk.i(), ijk.j(), ijk.k() ) ) return std::nullopt;
const size_t cellIndex = mainGrid->cellIndexFromIJK( ijk.i(), ijk.j(), ijk.k() );
const RigCell& cell = mainGrid->cell( cellIndex );
if ( cell.isInvalid() ) return std::nullopt;
using F = cvf::StructGridInterface;
ConnectionCellProperties properties;
properties.cellSize = cvf::Vec3d( ( cell.faceCenter( F::NEG_I ) - cell.faceCenter( F::POS_I ) ).length(),
( cell.faceCenter( F::NEG_J ) - cell.faceCenter( F::POS_J ) ).length(),
( cell.faceCenter( F::NEG_K ) - cell.faceCenter( F::POS_K ) ).length() );
properties.permeability = cvf::Vec3d( 1.0, 1.0, 1.0 );
if ( activeCellInfo && permeabilities[0] && permeabilities[1] && permeabilities[2] )
{
const size_t resultIndex = activeCellInfo->cellResultIndex( ReservoirCellIndex( cellIndex ) ).value();
if ( resultIndex < permeabilities[0]->size() && resultIndex < permeabilities[1]->size() && resultIndex < permeabilities[2]->size() )
{
properties.permeability =
cvf::Vec3d( ( *permeabilities[0] )[resultIndex], ( *permeabilities[1] )[resultIndex], ( *permeabilities[2] )[resultIndex] );
}
}
return properties;
};
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigSimulationInputConnectionTools::ConnectionKey::operator<( const ConnectionKey& other ) const
{
return std::make_tuple( wellName, ijk.x(), ijk.y(), ijk.z() ) <
std::make_tuple( other.wellName, other.ijk.x(), other.ijk.y(), other.ijk.z() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigSimulationInputConnectionTools::ConnectionDirectionMap
RigSimulationInputConnectionTools::connectionDirections( const std::vector<Opm::DeckKeyword>& compdatKeywords )
{
ConnectionDirectionMap directions;
for ( const auto& keyword : compdatKeywords )
{
for ( size_t r = 0; r < keyword.size(); ++r )
{
const auto& record = keyword.getRecord( r );
if ( record.size() < 5 ) continue;
const auto& wellItem = record.getItem( 0 );
if ( !wellItem.hasValue( 0 ) || wellItem.getType() != Opm::type_tag::string ) continue;
const auto wellName = wellItem.get<std::string>( 0 );
const int i = record.getItem( 1 ).get<int>( 0 ) - 1;
const int j = record.getItem( 2 ).get<int>( 0 ) - 1;
const int k1 = record.getItem( 3 ).get<int>( 0 ) - 1;
const int k2 = record.getItem( 4 ).get<int>( 0 ) - 1;
if ( i < 0 || j < 0 || k1 < 0 || k2 < k1 ) continue;
const auto dim = compdatPenetrationDimension( record );
for ( int k = k1; k <= k2; ++k )
{
directions[{ wellName, caf::VecIjk0( i, j, k ) }] = dim;
}
}
}
return directions;
}
//--------------------------------------------------------------------------------------------------
/// Cells where the well travels towards decreasing index along the connection penetration direction,
/// determined from the neighbouring COMPSEGS cells on the same branch sorted by measured depth
//--------------------------------------------------------------------------------------------------
RigSimulationInputConnectionTools::ConnectionKeySet
RigSimulationInputConnectionTools::reversedSegmentCells( const std::vector<Opm::DeckKeyword>& compsegsKeywords,
const ConnectionDirectionMap& directions )
{
struct SegmentCell
{
double startMd = 0.0;
caf::VecIjk0 ijk;
};
ConnectionKeySet reversed;
for ( const auto& keyword : compsegsKeywords )
{
if ( keyword.size() < 2 ) continue;
const auto& wellItem = keyword.getRecord( 0 ).getItem( 0 );
if ( !wellItem.hasValue( 0 ) || wellItem.getType() != Opm::type_tag::string ) continue;
const auto wellName = wellItem.get<std::string>( 0 );
std::map<int, std::vector<SegmentCell>> cellsPerBranch;
for ( size_t r = 1; r < keyword.size(); ++r )
{
const auto& record = keyword.getRecord( r );
if ( record.size() < 6 ) continue;
if ( !record.getItem( 0 ).hasValue( 0 ) || !record.getItem( 1 ).hasValue( 0 ) || !record.getItem( 2 ).hasValue( 0 ) ) continue;
const int i = record.getItem( 0 ).get<int>( 0 ) - 1;
const int j = record.getItem( 1 ).get<int>( 0 ) - 1;
const int k = record.getItem( 2 ).get<int>( 0 ) - 1;
if ( i < 0 || j < 0 || k < 0 ) continue;
const int branch = record.getItem( 3 ).hasValue( 0 ) ? record.getItem( 3 ).get<int>( 0 ) : 0;
const double startMd = record.getItem( 4 ).hasValue( 0 ) ? record.getItem( 4 ).get<double>( 0 ) : 0.0;
cellsPerBranch[branch].push_back( { startMd, caf::VecIjk0( i, j, k ) } );
}
for ( auto& [branch, cells] : cellsPerBranch )
{
std::stable_sort( cells.begin(), cells.end(), []( const SegmentCell& a, const SegmentCell& b ) { return a.startMd < b.startMd; } );
for ( size_t c = 0; c < cells.size(); ++c )
{
const auto& ijk = cells[c].ijk;
const auto key = ConnectionKey{ wellName, ijk };
auto dim = RigRefinement::DimK;
if ( auto it = directions.find( key ); it != directions.end() ) dim = it->second;
// Prefer the next cell along the well, fall back to the previous cell
std::optional<bool> isReversed;
if ( c + 1 < cells.size() && cells[c + 1].ijk[dim] != ijk[dim] )
{
isReversed = cells[c + 1].ijk[dim] < ijk[dim];
}
else if ( c > 0 && cells[c - 1].ijk[dim] != ijk[dim] )
{
isReversed = ijk[dim] < cells[c - 1].ijk[dim];
}
if ( isReversed.value_or( false ) ) reversed.insert( key );
}
}
}
return reversed;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<std::vector<Opm::DeckRecord>, QString>
RigSimulationInputConnectionTools::processCompdatRecord( const Opm::DeckRecord& record,
const std::string& wellName,
const RigSimulationInputSettings& settings,
const ConnectionCellPropertiesProvider& cellProperties )
{
// COMPDAT format: WELL I J K1 K2 STATE ...
// Items: 0=WELL, 1=I, 2=J, 3=K1, 4=K2
if ( record.size() < 5 )
{
return std::unexpected( QString( "COMPDAT record for well %1 has insufficient items (expected at least 5, got %2)" )
.arg( wellName.c_str() )
.arg( record.size() ) );
}
int origI = record.getItem( 1 ).get<int>( 0 ) - 1; // Convert to 0-based
int origJ = record.getItem( 2 ).get<int>( 0 ) - 1;
int origK1 = record.getItem( 3 ).get<int>( 0 ) - 1;
int origK2 = record.getItem( 4 ).get<int>( 0 ) - 1;
const auto& refinement = settings.refinement();
auto transformResultK1 = RigGridExportAdapter::transformIjkToSectorCoordinates( caf::VecIjk0( origI, origJ, origK1 ),
settings.min(),
settings.max(),
refinement,
true );
if ( !transformResultK1 )
{
return std::unexpected(
QString( "COMPDAT K1 coordinate for well %1 is out of sector bounds: %2" ).arg( wellName.c_str() ).arg( transformResultK1.error() ) );
}
auto transformResultK2 = RigGridExportAdapter::transformIjkToSectorCoordinates( caf::VecIjk0( origI, origJ, origK2 ),
settings.min(),
settings.max(),
refinement,
true );
if ( !transformResultK2 )
{
return std::unexpected(
QString( "COMPDAT K2 coordinate for well %1 is out of sector bounds: %2" ).arg( wellName.c_str() ).arg( transformResultK2.error() ) );
}
if ( origK2 < origK1 )
{
return std::unexpected(
QString( "COMPDAT record for well %1 has K2 (%2) < K1 (%3)" ).arg( wellName.c_str() ).arg( origK2 + 1 ).arg( origK1 + 1 ) );
}
using C = Opm::ParserKeywords::COMPDAT;
const auto dim = compdatPenetrationDimension( record );
const size_t sectorI = origI - settings.min().x();
const size_t sectorJ = origJ - settings.min().y();
const size_t sectorK1 = origK1 - settings.min().z();
const size_t sectorK2 = origK2 - settings.min().z();
const bool hasExplicitFactors = hasExplicitDoubleValue( record, C::CONNECTION_TRANSMISSIBILITY_FACTOR::itemName ) ||
hasExplicitDoubleValue( record, C::Kh::itemName );
bool isKRefinedInRange = false;
for ( size_t k = sectorK1; k <= sectorK2; ++k )
{
if ( refinement.subcellCount( RigRefinement::DimK, k ) > 1 ) isKRefinedInRange = true;
}
struct KSpan
{
size_t k1 = 0;
size_t k2 = 0;
double fraction = 1.0;
size_t sectorK = 0; // Original (sector-relative) K, used for the lateral refinement correction
};
std::vector<KSpan> kSpans;
if ( !isKRefinedInRange )
{
kSpans.push_back( { refinement.cumulativeOffset( RigRefinement::DimK, sectorK1 ),
refinement.cumulativeOffset( RigRefinement::DimK, sectorK2 ),
1.0,
sectorK1 } );
}
else if ( dim == RigRefinement::DimK && !hasExplicitFactors )
{
// The simulator computes the connection factors, so the full refined range can be given in one record
kSpans.push_back( { refinement.cumulativeOffset( RigRefinement::DimK, sectorK1 ),
refinement.cumulativeOffset( RigRefinement::DimK, sectorK2 + 1 ) - 1,
1.0,
sectorK1 } );
}
else if ( dim == RigRefinement::DimK )
{
for ( size_t k = sectorK1; k <= sectorK2; ++k )
{
for ( const auto& subcell : allSubcells( refinement, RigRefinement::DimK, k ) )
{
kSpans.push_back( { subcell.index, subcell.index, subcell.fraction, k } );
}
}
}
else
{
for ( size_t k = sectorK1; k <= sectorK2; ++k )
{
const size_t centerK = centerSubcellIndex( refinement, RigRefinement::DimK, k );
kSpans.push_back( { centerK, centerK, 1.0, k } );
}
}
const std::vector<RefinedSubcell> iSubcells =
( dim == RigRefinement::DimI )
? allSubcells( refinement, RigRefinement::DimI, sectorI )
: std::vector<RefinedSubcell>{ { centerSubcellIndex( refinement, RigRefinement::DimI, sectorI ), 1.0 } };
const std::vector<RefinedSubcell> jSubcells =
( dim == RigRefinement::DimJ )
? allSubcells( refinement, RigRefinement::DimJ, sectorJ )
: std::vector<RefinedSubcell>{ { centerSubcellIndex( refinement, RigRefinement::DimJ, sectorJ ), 1.0 } };
// Refinement perpendicular to the well, applies to the center sub-cell
auto lateralCorrection = [&]( size_t sectorK ) -> std::optional<LateralRefinementCorrection>
{
double fraction1 = 1.0;
double fraction2 = 1.0;
if ( dim == RigRefinement::DimK )
{
fraction1 = centerSubcellFraction( refinement, RigRefinement::DimI, sectorI );
fraction2 = centerSubcellFraction( refinement, RigRefinement::DimJ, sectorJ );
}
else
{
fraction1 = ( dim == RigRefinement::DimI ) ? centerSubcellFraction( refinement, RigRefinement::DimJ, sectorJ )
: centerSubcellFraction( refinement, RigRefinement::DimI, sectorI );
fraction2 = centerSubcellFraction( refinement, RigRefinement::DimK, sectorK );
}
if ( fraction1 >= 1.0 && fraction2 >= 1.0 ) return std::nullopt;
std::optional<ConnectionCellProperties> properties;
if ( cellProperties )
{
properties = cellProperties( caf::VecIjk0( origI, origJ, sectorK + settings.min().z() ) );
}
return lateralRefinementCorrection( record, wellName, dim, fraction1, fraction2, properties );
};
auto createRecord =
[&]( size_t i, size_t j, size_t k1, size_t k2, double fraction, const std::optional<LateralRefinementCorrection>& correction )
{
std::vector<Opm::DeckItem> items;
items.push_back( record.getItem( 0 ) );
items.push_back( RifOpmDeckTools::item( C::I::itemName, static_cast<int>( i + 1 ) ) );
items.push_back( RifOpmDeckTools::item( C::J::itemName, static_cast<int>( j + 1 ) ) );
items.push_back( RifOpmDeckTools::item( C::K1::itemName, static_cast<int>( k1 + 1 ) ) );
items.push_back( RifOpmDeckTools::item( C::K2::itemName, static_cast<int>( k2 + 1 ) ) );
for ( size_t idx = 5; idx < record.size(); ++idx )
{
const auto& item = record.getItem( idx );
const bool isCF = item.name() == C::CONNECTION_TRANSMISSIBILITY_FACTOR::itemName;
const bool isScalable = isCF || item.name() == C::Kh::itemName;
const bool isExplicit = hasExplicitDoubleValue( record, item.name() );
if ( isCF && isExplicit && correction )
{
items.push_back(
RifOpmDeckTools::item( item.name(), item.get<double>( 0 ) * fraction * correction->connectionFactorMultiplier ) );
}
else if ( item.name() == C::PR::itemName && isExplicit && correction )
{
items.push_back( RifOpmDeckTools::item( item.name(), correction->pressureEquivalentRadius ) );
}
else if ( fraction != 1.0 && isScalable && isExplicit )
{
items.push_back( RifOpmDeckTools::item( item.name(), item.get<double>( 0 ) * fraction ) );
}
else
{
items.push_back( item );
}
}
return Opm::DeckRecord{ std::move( items ) };
};
std::vector<Opm::DeckRecord> records;
for ( const auto& kSpan : kSpans )
{
const auto correction = lateralCorrection( kSpan.sectorK );
for ( const auto& iSubcell : iSubcells )
{
for ( const auto& jSubcell : jSubcells )
{
records.push_back( createRecord( iSubcell.index,
jSubcell.index,
kSpan.k1,
kSpan.k2,
kSpan.fraction * iSubcell.fraction * jSubcell.fraction,
correction ) );
}
}
}
return records;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<std::vector<Opm::DeckRecord>, QString>
RigSimulationInputConnectionTools::processCompsegsRecord( const Opm::DeckRecord& record,
const std::string& wellName,
bool isWellNameRecord,
const RigSimulationInputSettings& settings,
const ConnectionDirectionMap& directions,
const ConnectionKeySet& reversedCells )
{
// COMPSEGS format: first record is well name, subsequent records are segment data
// Well name record: just copy as-is
if ( isWellNameRecord )
{
return std::vector<Opm::DeckRecord>{ Opm::DeckRecord( record ) };
}
// Segment record format: I J K BRANCH START_MD END_MD DIRECTION END_IJK ...
// Items: 0=I, 1=J, 2=K, 4=START_MD, 5=END_MD, 6=DIRECTION
if ( record.size() < 3 )
{
return std::unexpected( QString( "COMPSEGS segment record for well %1 has insufficient items (expected at least 3, got %2)" )
.arg( wellName.c_str() )
.arg( record.size() ) );
}
caf::VecIjk0 origIjk = RigSimulationInputTool::extractIjk( record, 0, 1, 2 ).toZeroBased();
const auto& refinement = settings.refinement();
auto transformResult = RigGridExportAdapter::transformIjkToSectorCoordinates( origIjk, settings.min(), settings.max(), refinement, true );
if ( !transformResult )
{
return std::unexpected(
QString( "COMPSEGS segment coordinate for well %1 is out of sector bounds: %2" ).arg( wellName.c_str() ).arg( transformResult.error() ) );
}
auto dim = RigRefinement::DimK;
if ( auto it = directions.find( { wellName, origIjk } ); it != directions.end() )
{
dim = it->second;
}
const size_t sectorIndex = origIjk[dim] - settings.min()[dim];
std::vector<RefinedSubcell> subcells = { { ( *transformResult )[dim], 1.0 } };
using C = Opm::ParserKeywords::COMPSEGS;
auto hasItemValue = [&record]( const std::string& itemName )
{ return record.hasItem( itemName ) && record.getItem( itemName ).hasValue( 0 ); };
// A cell range is given by END_IJK. DIRECTION alone (with an MD interval) describes a single cell.
const bool isCellRange = hasItemValue( C::END_IJK::itemName );
const bool hasEndMd = hasItemValue( C::DISTANCE_END::itemName );
// A defaulted start MD is zero (only allowed for the first segment record)
const double startMd = hasItemValue( C::DISTANCE_START::itemName ) ? record.getItem( C::DISTANCE_START::itemName ).get<double>( 0 ) : 0.0;
const double endMd = hasEndMd ? record.getItem( C::DISTANCE_END::itemName ).get<double>( 0 ) : 0.0;
// Whether the MD interval is divided between the sub-cells, otherwise the record is copied unchanged to each sub-cell
bool splitMdInterval = false;
if ( refinement.subcellCount( dim, sectorIndex ) > 1 )
{
if ( isCellRange )
{
RiaLogging::warning(
std::format( "COMPSEGS record for well {} at ({}, {}, {}) uses a cell range (END_IJK), which is not supported "
"for refined sector export. Only the center refined sub-cell is assigned to the segment.",
wellName,
origIjk.x() + 1,
origIjk.y() + 1,
origIjk.z() + 1 ) );
}
else
{
// All refined sub-cells along the connection are completed by COMPDAT, so all must be assigned to a segment
subcells = allSubcells( refinement, dim, sectorIndex );
// Assign the MD interval in the order the well travels through the sub-cells
if ( reversedCells.contains( { wellName, origIjk } ) )
{
std::reverse( subcells.begin(), subcells.end() );
}
splitMdInterval = hasEndMd;
if ( !hasEndMd )
{
RiaLogging::warning( std::format( "COMPSEGS record for well {} at ({}, {}, {}) has no end MD. "
"The record is copied unchanged to each refined sub-cell.",
wellName,
origIjk.x() + 1,
origIjk.y() + 1,
origIjk.z() + 1 ) );
}
}
}
std::vector<Opm::DeckRecord> records;
double startFraction = 0.0;
for ( const auto& subcell : subcells )
{
caf::VecIjk0 refinedIjk = *transformResult;
refinedIjk[dim] = subcell.index;
std::vector<Opm::DeckItem> items;
items.push_back( RifOpmDeckTools::item( "I", static_cast<int>( refinedIjk.x() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( "J", static_cast<int>( refinedIjk.y() + 1 ) ) );
items.push_back( RifOpmDeckTools::item( "K", static_cast<int>( refinedIjk.z() + 1 ) ) );
const double endFraction = startFraction + subcell.fraction;
for ( size_t i = 3; i < record.size(); ++i )
{
const auto& itemName = record.getItem( i ).name();
if ( splitMdInterval && ( itemName == C::DISTANCE_START::itemName || itemName == C::DISTANCE_END::itemName ) )
{
const double fraction = ( itemName == C::DISTANCE_START::itemName ) ? startFraction : endFraction;
items.push_back( RifOpmDeckTools::item( record.getItem( i ).name(), startMd + ( endMd - startMd ) * fraction ) );
}
else
{
items.push_back( record.getItem( i ) );
}
}
startFraction = endFraction;
records.push_back( Opm::DeckRecord{ std::move( items ) } );
}
return records;
}
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.h
+99
View File
@@ -0,0 +1,99 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2025- 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 "RigRefinement.h"
#include "cafVecIjk.h"
#include "cvfVector3.h"
#include <QString>
#include <expected>
#include <functional>
#include <map>
#include <optional>
#include <set>
#include <string>
#include <vector>
class RimEclipseCase;
class RigSimulationInputSettings;
namespace Opm
{
class DeckRecord;
class DeckKeyword;
} // namespace Opm
//==================================================================================================
///
/// Well connection keyword (COMPDAT, COMPSEGS) handling when exporting simulation input files
///
//==================================================================================================
class RigSimulationInputConnectionTools
{
public:
// Penetration direction of well connections, keyed on (well name, original 0-based IJK)
struct ConnectionKey
{
std::string wellName;
caf::VecIjk0 ijk;
bool operator<( const ConnectionKey& other ) const;
};
using ConnectionKeySet = std::set<ConnectionKey>;
using ConnectionDirectionMap = std::map<ConnectionKey, RigRefinement::Dimension>;
static ConnectionDirectionMap connectionDirections( const std::vector<Opm::DeckKeyword>& compdatKeywords );
// Segment cells (well name, original 0-based IJK) where the well travels towards decreasing index along the
// connection penetration direction, derived from neighbouring COMPSEGS cells on the same branch ordered by MD
static ConnectionKeySet reversedSegmentCells( const std::vector<Opm::DeckKeyword>& compsegsKeywords,
const ConnectionDirectionMap& directions );
// Size and permeability of an original grid cell, ordered I/X, J/Y, K/Z
struct ConnectionCellProperties
{
cvf::Vec3d cellSize;
cvf::Vec3d permeability;
};
using ConnectionCellPropertiesProvider = std::function<std::optional<ConnectionCellProperties>( const caf::VecIjk0& originalIjk )>;
static ConnectionCellPropertiesProvider connectionCellPropertiesProvider( RimEclipseCase* eclipseCase );
// Refined connections along the penetration direction are split into one record per sub-cell with
// explicit CF and KH scaled by the sub-cell length fraction. Perpendicular directions use the center sub-cell,
// and explicit CF is corrected for the smaller Peaceman pressure equivalent radius of the refined cell.
static std::expected<std::vector<Opm::DeckRecord>, QString>
processCompdatRecord( const Opm::DeckRecord& record,
const std::string& wellName,
const RigSimulationInputSettings& settings,
const ConnectionCellPropertiesProvider& cellProperties = {} );
// Segment records are split into one record per sub-cell along the connection penetration direction,
// with the MD interval divided by the sub-cell length fractions in the order the well travels through them.
// Without an end MD the record is copied unchanged to each sub-cell. Cell ranges (END_IJK) use the center sub-cell.
static std::expected<std::vector<Opm::DeckRecord>, QString> processCompsegsRecord( const Opm::DeckRecord& record,
const std::string& wellName,
bool isWellNameRecord,
const RigSimulationInputSettings& settings,
const ConnectionDirectionMap& directions,
const ConnectionKeySet& reversedCells = {} );
};
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.cpp
+279
View File
@@ -0,0 +1,279 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2025- 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 "RigSimulationInputNncTools.h"
#include "RiaLogging.h"
#include "RifOpmFlowDeckFile.h"
#include "ProjectDataModel/Jobs/RimKeywordFactory.h"
#include "RimEclipseCase.h"
#include "RigBoundingBoxIjk.h"
#include "RigEclipseCaseData.h"
#include "RigMainGrid.h"
#include "RigSimulationInputSettings.h"
#include "RigSimulationInputTool.h"
#include "opm/input/eclipse/Deck/DeckItem.hpp"
#include "opm/input/eclipse/Deck/DeckKeyword.hpp"
#include "opm/input/eclipse/Deck/DeckRecord.hpp"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RigSimulationInputNncTools::NNCConnection>
RigSimulationInputNncTools::extractDeckEditNncConnections( RifOpmFlowDeckFile& deckFile, const RigMainGrid& mainGrid )
{
std::vector<NNCConnection> connections;
// Find all EDITNNC keywords in deck (there may be multiple)
auto keywords = deckFile.findAllKeywordsWithIndices( "EDITNNC" );
for ( const auto& keywordPair : keywords )
{
const auto& keyword = keywordPair.second;
// Each keyword can have multiple records
for ( size_t recordIdx = 0; recordIdx < keyword.size(); ++recordIdx )
{
const auto& record = keyword.getRecord( recordIdx );
// Read IJK coordinates directly from deck (1-based) and convert to 0-based
caf::VecIjk1 ijk1( record.getItem( 0 ).get<int>( 0 ), record.getItem( 1 ).get<int>( 0 ), record.getItem( 2 ).get<int>( 0 ) );
caf::VecIjk1 ijk2( record.getItem( 3 ).get<int>( 0 ), record.getItem( 4 ).get<int>( 0 ), record.getItem( 5 ).get<int>( 0 ) );
double transMult = record.getItem( 6 ).get<double>( 0 );
// Convert IJK to global cell indices
size_t c1Idx = mainGrid.cellIndexFromIJK( ijk1.toZeroBased() );
size_t c2Idx = mainGrid.cellIndexFromIJK( ijk2.toZeroBased() );
// Validate indices
if ( c1Idx == cvf::UNDEFINED_SIZE_T || c2Idx == cvf::UNDEFINED_SIZE_T )
{
RiaLogging::warning( std::format( "Invalid EDITNNC connection in deck: ({})-({})", ijk1.toString(), ijk2.toString() ) );
continue;
}
connections.push_back( { c1Idx, c2Idx, transMult } );
}
}
return connections;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RigSimulationInputNncTools::NNCConnection>
RigSimulationInputNncTools::filterInternalSectorConnections( const std::vector<NNCConnection>& allConnections,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const caf::VecIjk0& max )
{
std::vector<NNCConnection> filteredConnections;
RigBoundingBoxIjk sectorBounds( min, max );
for ( const auto& conn : allConnections )
{
// Get IJK for both cells
auto ijk1 = mainGrid.ijkFromCellIndex( conn.c1GlobIdx );
auto ijk2 = mainGrid.ijkFromCellIndex( conn.c2GlobIdx );
// Only keep if BOTH cells are inside sector bounds
if ( ijk1 && ijk2 && sectorBounds.contains( *ijk1 ) && sectorBounds.contains( *ijk2 ) )
{
filteredConnections.push_back( conn );
}
}
return filteredConnections;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<RigSimulationInputNncTools::TransformedNNCConnection, QString>
RigSimulationInputNncTools::transformNNCToSectorCoordinates( const NNCConnection& connection,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const RigRefinement& refinement )
{
// Get global IJK for both cells
auto gijk1 = mainGrid.ijkFromCellIndex( connection.c1GlobIdx );
if ( !gijk1 )
{
return std::unexpected( QString( "Invalid cell index: %1" ).arg( connection.c1GlobIdx ) );
}
auto gijk2 = mainGrid.ijkFromCellIndex( connection.c2GlobIdx );
if ( !gijk2 )
{
return std::unexpected( QString( "Invalid cell index: %1" ).arg( connection.c2GlobIdx ) );
}
// Transform to sector coordinates using utility method
TransformedNNCConnection transformed;
transformed.cell1 = RigSimulationInputTool::transformToSectorCoordinates( *gijk1, min, refinement );
transformed.cell2 = RigSimulationInputTool::transformToSectorCoordinates( *gijk2, min, refinement );
transformed.transmissibility = connection.transmissibility;
return transformed;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RigSimulationInputNncTools::TransformedNNCConnection>
RigSimulationInputNncTools::refineEditNncConnection( const NNCConnection& connection,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const RigRefinement& refinement )
{
// Get global IJK for both cells
auto gijk1 = mainGrid.ijkFromCellIndex( connection.c1GlobIdx );
auto gijk2 = mainGrid.ijkFromCellIndex( connection.c2GlobIdx );
if ( !gijk1 || !gijk2 )
{
return {}; // Return empty on error
}
// For EDITNNC: Connect subcells at same relative positions within coarse cells
// Example: refinement (2,2,1) creates 4 connections:
// (0,0,0) → (0,0,0)
// (0,1,0) → (0,1,0)
// (1,0,0) → (1,0,0)
// (1,1,0) → (1,1,0)
// Each connection gets the SAME TRANSMUL (no division)
// Get base sector coordinates for both cells (start of refined block)
caf::VecIjk0 baseCell1 = RigSimulationInputTool::transformToSectorCoordinates( *gijk1, min, refinement );
caf::VecIjk0 baseCell2 = RigSimulationInputTool::transformToSectorCoordinates( *gijk2, min, refinement );
// Get subcell counts for each cell
size_t rel1I = gijk1->i() - min.i();
size_t rel1J = gijk1->j() - min.j();
size_t rel1K = gijk1->k() - min.k();
size_t subcellsI = refinement.subcellCount( RigRefinement::DimI, rel1I );
size_t subcellsJ = refinement.subcellCount( RigRefinement::DimJ, rel1J );
size_t subcellsK = refinement.subcellCount( RigRefinement::DimK, rel1K );
std::vector<TransformedNNCConnection> refined;
for ( size_t subI = 0; subI < subcellsI; ++subI )
{
for ( size_t subJ = 0; subJ < subcellsJ; ++subJ )
{
for ( size_t subK = 0; subK < subcellsK; ++subK )
{
// Add subcell offset to base coordinates
TransformedNNCConnection conn;
conn.cell1 = caf::VecIjk0( baseCell1.i() + subI, baseCell1.j() + subJ, baseCell1.k() + subK );
// SAME relative position in cell2
conn.cell2 = caf::VecIjk0( baseCell2.i() + subI, baseCell2.j() + subJ, baseCell2.k() + subK );
conn.transmissibility = connection.transmissibility; // Same TRAN_MULT for all
refined.push_back( conn );
}
}
}
return refined;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::expected<void, QString> RigSimulationInputNncTools::exportEditNncKeyword( RimEclipseCase* eclipseCase,
const RigSimulationInputSettings& settings,
RifOpmFlowDeckFile& deckFile )
{
RigMainGrid* mainGrid = eclipseCase->eclipseCaseData()->mainGrid();
if ( !mainGrid ) return std::unexpected( "No main grid available" );
// Step 1: Extract EDITNNC from deck only
auto allConnections = extractDeckEditNncConnections( deckFile, *mainGrid );
if ( allConnections.empty() )
{
RiaLogging::info( "No EDITNNC keywords found in deck - skipping EDITNNC export" );
return {};
}
// Step 2: Filter for internal sector
auto sectorConnections = filterInternalSectorConnections( allConnections, *mainGrid, settings.min(), settings.max() );
if ( sectorConnections.empty() )
{
RiaLogging::info( "No internal EDITNNC connections in sector - skipping EDITNNC export" );
return {};
}
RiaLogging::info(
std::format( "Found {} internal EDITNNC connections in sector (out of {} total)", sectorConnections.size(), allConnections.size() ) );
// Step 3: Transform and refine
bool hasRef = settings.hasRefinement();
std::vector<TransformedNNCConnection> transformedConnections;
if ( hasRef )
{
// Refine to subcells with corresponding positions
for ( const auto& conn : sectorConnections )
{
auto refined = refineEditNncConnection( conn, *mainGrid, settings.min(), settings.refinement() );
transformedConnections.insert( transformedConnections.end(), refined.begin(), refined.end() );
}
}
else
{
// Simple coordinate transformation, no refinement
for ( const auto& conn : sectorConnections )
{
auto result = transformNNCToSectorCoordinates( conn, *mainGrid, settings.min(), settings.refinement() );
if ( result )
{
transformedConnections.push_back( *result );
}
else
{
RiaLogging::warning( result.error().toStdString() );
}
}
}
if ( transformedConnections.empty() )
{
return std::unexpected( "Failed to transform any EDITNNC connections" );
}
RiaLogging::info( std::format( "Exporting {} EDITNNC connections to sector model", transformedConnections.size() ) );
// Step 4: Create keyword
auto editnncKw = RimKeywordFactory::editnncKeyword( transformedConnections );
// Step 5: Replace/add to EDIT section (not GRID!)
deckFile.replaceKeyword( "EDIT", editnncKw );
return {};
}
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.h
+75
View File
@@ -0,0 +1,75 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2025- 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 "RigRefinement.h"
#include "cafVecIjk.h"
#include <QString>
#include <expected>
#include <vector>
class RimEclipseCase;
class RigSimulationInputSettings;
class RifOpmFlowDeckFile;
class RigMainGrid;
//==================================================================================================
///
/// EDITNNC handling when exporting simulation input files (sector models)
///
//==================================================================================================
class RigSimulationInputNncTools
{
public:
/// Raw NNC connection with global cell indices
struct NNCConnection
{
size_t c1GlobIdx; // Global cell index for cell 1
size_t c2GlobIdx; // Global cell index for cell 2
double transmissibility; // Transmissibility value
};
/// Transformed NNC connection with sector-relative 0-based IJK coordinates
struct TransformedNNCConnection
{
caf::VecIjk0 cell1 = caf::VecIjk0::ZERO; // 0-based sector IJK for cell 1
caf::VecIjk0 cell2 = caf::VecIjk0::ZERO; // 0-based sector IJK for cell 2
double transmissibility = 0.0; // Transmissibility value
};
static std::expected<void, QString>
exportEditNncKeyword( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::vector<NNCConnection> extractDeckEditNncConnections( RifOpmFlowDeckFile& deckFile, const RigMainGrid& mainGrid );
static std::vector<NNCConnection> filterInternalSectorConnections( const std::vector<NNCConnection>& allConnections,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const caf::VecIjk0& max );
static std::expected<TransformedNNCConnection, QString> transformNNCToSectorCoordinates( const NNCConnection& connection,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const RigRefinement& refinement );
static std::vector<TransformedNNCConnection> refineEditNncConnection( const NNCConnection& connection,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const RigRefinement& refinement );
};
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp
+19 -1185
View File
File diff suppressed because it is too large. Load diff
ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.h
+22 -126
View File
@@ -28,11 +28,8 @@
#include <expected>
#include <functional>
#include <map>
#include <optional>
#include <set>
#include <string>
#include <tuple>
#include <vector>
class RimEclipseView;
@@ -60,24 +57,6 @@ class DeckKeyword;
class RigSimulationInputTool
{
public:
// Data structures for NNC export
/// Raw NNC connection with global cell indices
struct NNCConnection
{
size_t c1GlobIdx; // Global cell index for cell 1
size_t c2GlobIdx; // Global cell index for cell 2
double transmissibility; // Transmissibility value
};
/// Transformed NNC connection with sector-relative 0-based IJK coordinates
struct TransformedNNCConnection
{
caf::VecIjk0 cell1 = caf::VecIjk0::ZERO; // 0-based sector IJK for cell 1
caf::VecIjk0 cell2 = caf::VecIjk0::ZERO; // 0-based sector IJK for cell 2
double transmissibility = 0.0; // Transmissibility value
};
static std::expected<void, QString>
exportSimulationInput( RimEclipseCase& eclipseCase, const RigSimulationInputSettings& settings, cvf::UByteArray* visibility );
@@ -96,86 +75,38 @@ public:
static std::expected<Opm::DeckRecord, QString>
processAddRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement );
static std::expected<Opm::DeckRecord, QString>
processAquconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement );
static std::expected<Opm::DeckRecord, QString>
processAquanconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement );
static std::expected<Opm::DeckRecord, QString>
processAqunumRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement );
// Penetration direction of well connections, keyed on (well name, original 0-based IJK)
struct ConnectionKey
{
std::string wellName;
caf::VecIjk0 ijk;
bool operator<( const ConnectionKey& other ) const;
};
using ConnectionKeySet = std::set<ConnectionKey>;
using ConnectionDirectionMap = std::map<ConnectionKey, RigRefinement::Dimension>;
static ConnectionDirectionMap connectionDirections( const std::vector<Opm::DeckKeyword>& compdatKeywords );
// Segment cells (well name, original 0-based IJK) where the well travels towards decreasing index along the
// connection penetration direction, derived from neighbouring COMPSEGS cells on the same branch ordered by MD
static ConnectionKeySet reversedSegmentCells( const std::vector<Opm::DeckKeyword>& compsegsKeywords,
const ConnectionDirectionMap& directions );
// Size and permeability of an original grid cell, ordered I/X, J/Y, K/Z
struct ConnectionCellProperties
{
cvf::Vec3d cellSize;
cvf::Vec3d permeability;
};
using ConnectionCellPropertiesProvider = std::function<std::optional<ConnectionCellProperties>( const caf::VecIjk0& originalIjk )>;
static ConnectionCellPropertiesProvider connectionCellPropertiesProvider( RimEclipseCase* eclipseCase );
// Refined connections along the penetration direction are split into one record per sub-cell with
// explicit CF and KH scaled by the sub-cell length fraction. Perpendicular directions use the center sub-cell,
// and explicit CF is corrected for the smaller Peaceman pressure equivalent radius of the refined cell.
static std::expected<std::vector<Opm::DeckRecord>, QString>
processCompdatRecord( const Opm::DeckRecord& record,
const std::string& wellName,
const RigSimulationInputSettings& settings,
const ConnectionCellPropertiesProvider& cellProperties = {} );
// Segment records are split into one record per sub-cell along the connection penetration direction,
// with the MD interval divided by the sub-cell length fractions in the order the well travels through them.
// Without an end MD the record is copied unchanged to each sub-cell. Cell ranges (END_IJK) use the center sub-cell.
static std::expected<std::vector<Opm::DeckRecord>, QString> processCompsegsRecord( const Opm::DeckRecord& record,
const std::string& wellName,
bool isWellNameRecord,
const RigSimulationInputSettings& settings,
const ConnectionDirectionMap& directions,
const ConnectionKeySet& reversedCells = {} );
static std::vector<int> createRefinedVisibility( const RigGridExportAdapter& gridAdapter );
static std::expected<void, QString>
transformKeywordsInDeckFile( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
// Public EDITNNC methods for testing
static std::vector<NNCConnection> extractDeckEditNncConnections( RifOpmFlowDeckFile& deckFile, const RigMainGrid& mainGrid );
static std::vector<NNCConnection> filterInternalSectorConnections( const std::vector<NNCConnection>& allConnections,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const caf::VecIjk0& max );
static std::expected<TransformedNNCConnection, QString> transformNNCToSectorCoordinates( const NNCConnection& connection,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const RigRefinement& refinement );
static std::vector<TransformedNNCConnection> refineEditNncConnection( const NNCConnection& connection,
const RigMainGrid& mainGrid,
const caf::VecIjk0& min,
const RigRefinement& refinement );
// Utility: Transform global IJK to sector-relative coordinates with refinement
static caf::VecIjk0 transformToSectorCoordinates( const caf::VecIjk0& globalIjk, const caf::VecIjk0& min, const RigRefinement& refinement );
// Generic helper for processing keywords with box indices
using RecordProcessorFunc =
std::function<std::expected<Opm::DeckRecord, QString>( const Opm::DeckRecord&, const caf::VecIjk0&, const caf::VecIjk0&, const RigRefinement& )>;
static std::expected<void, QString> replaceKeywordWithBoxIndices( const std::string& keywordName,
RimEclipseCase* eclipseCase,
const RigSimulationInputSettings& settings,
RifOpmFlowDeckFile& deckFile,
RecordProcessorFunc processorFunc );
// Helper function to transform bounding box from global to sector coordinates
// Returns bounding box with 0-based sector-relative coordinates
static std::expected<RigBoundingBoxIjk<caf::VecIjk0>, QString>
transformBoxToSectorCoordinates( const RigBoundingBoxIjk<caf::VecIjk0>& inputBox,
const caf::VecIjk0& sectorMin,
const caf::VecIjk0& sectorMax,
const RigRefinement& refinement,
const QString& keywordName,
const QString& recordIdentifier );
// Helper function to extract IJK coordinates from a deck record
// Returns a VecIjk1 (1-based) constructed from values at the specified item indices
static caf::VecIjk1 extractIjk( const Opm::DeckRecord& record, size_t indexI, size_t indexJ, size_t indexK );
private:
static std::expected<void, QString> updateCornerPointGridInDeckFile( RimEclipseCase* eclipseCase,
const RigSimulationInputSettings& settings,
@@ -202,15 +133,6 @@ private:
cvf::ref<cvf::UByteArray> visibility,
RifOpmFlowDeckFile& deckFile );
static std::expected<void, QString>
replaceAquconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::expected<void, QString>
replaceAquanconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::expected<void, QString>
replaceAqunumKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::expected<void, QString>
addFaultsToDeckFile( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
@@ -225,14 +147,8 @@ private:
static std::set<std::string> wellNamesToInclude( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings );
static std::expected<void, QString>
exportEditNncKeyword( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile );
static std::expected<void, QString> applyModelPadding( RifOpmFlowDeckFile& deckFile, const RigModelPaddingSettings& paddingSettings );
// Utility: Transform global IJK to sector-relative coordinates with refinement
static caf::VecIjk0 transformToSectorCoordinates( const caf::VecIjk0& globalIjk, const caf::VecIjk0& min, const RigRefinement& refinement );
static std::expected<void, QString> addOperNumRegionAndOperater( RifOpmFlowDeckFile& deckFile,
const RigGridExportAdapter& gridAdapter,
const std::vector<int>& operNumResult,
@@ -243,24 +159,4 @@ private:
static std::expected<Opm::DeckRecord, QString>
processWelspecsRecord( const Opm::DeckRecord& record, const std::string& wellName, const RigSimulationInputSettings& settings );
static std::expected<void, QString> replaceKeywordWithBoxIndices( const std::string& keywordName,
RimEclipseCase* eclipseCase,
const RigSimulationInputSettings& settings,
RifOpmFlowDeckFile& deckFile,
RecordProcessorFunc processorFunc );
// Helper function to transform bounding box from global to sector coordinates
// Returns bounding box with 0-based sector-relative coordinates
static std::expected<RigBoundingBoxIjk<caf::VecIjk0>, QString>
transformBoxToSectorCoordinates( const RigBoundingBoxIjk<caf::VecIjk0>& inputBox,
const caf::VecIjk0& sectorMin,
const caf::VecIjk0& sectorMax,
const RigRefinement& refinement,
const QString& keywordName,
const QString& recordIdentifier );
// Helper function to extract IJK coordinates from a deck record
// Returns a VecIjk1 (1-based) constructed from values at the specified item indices
static caf::VecIjk1 extractIjk( const Opm::DeckRecord& record, size_t indexI, size_t indexJ, size_t indexK );
};