diff --git a/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.cpp b/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.cpp index 349b75df4c..59c1d47296 100644 --- a/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.cpp +++ b/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.cpp @@ -30,7 +30,7 @@ #include "RigEclipseResultTools.h" #include "RigFault.h" #include "RigMainGrid.h" -#include "RigSimulationInputTool.h" +#include "RigSimulationInputNncTools.h" #include "RimEclipseCase.h" #include "RimWellEventWellSpec.h" @@ -752,7 +752,7 @@ Opm::DeckKeyword operaterKeyword( std::string targetProperty, //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -Opm::DeckKeyword editnncKeyword( const std::vector& connections ) +Opm::DeckKeyword editnncKeyword( const std::vector& connections ) { using E = Opm::ParserKeywords::EDITNNC; diff --git a/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.h b/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.h index 250afb9ae2..67e131a198 100644 --- a/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.h +++ b/ApplicationLibCode/ProjectDataModel/Jobs/RimKeywordFactory.h @@ -18,7 +18,7 @@ #pragma once -#include "ReservoirDataModel/SimulationFile/RigSimulationInputTool.h" +#include "ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.h" #include "cvfVector3.h" @@ -79,7 +79,7 @@ Opm::DeckKeyword operaterKeyword( std::string targetProperty, std::string inputProperty, std::optional alpha = std::nullopt, std::optional beta = std::nullopt ); -Opm::DeckKeyword editnncKeyword( const std::vector& connections ); +Opm::DeckKeyword editnncKeyword( const std::vector& connections ); Opm::DeckKeyword datesKeyword( const QDateTime& date ); QString deckKeywordToString( const Opm::DeckKeyword& keyword ); diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/CMakeLists_files.cmake b/ApplicationLibCode/ReservoirDataModel/SimulationFile/CMakeLists_files.cmake index c788d0c0a6..1302838a61 100644 --- a/ApplicationLibCode/ReservoirDataModel/SimulationFile/CMakeLists_files.cmake +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/CMakeLists_files.cmake @@ -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 ) diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.cpp b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.cpp new file mode 100644 index 0000000000..d7de923719 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.cpp @@ -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 +// 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 RigSimulationInputAquiferTools::replaceAquconKeywordIndices( RimEclipseCase* eclipseCase, + const RigSimulationInputSettings& settings, + RifOpmFlowDeckFile& deckFile ) +{ + return RigSimulationInputTool::replaceKeywordWithBoxIndices( "AQUCON", eclipseCase, settings, deckFile, processAquconRecord ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::expected RigSimulationInputAquiferTools::replaceAquanconKeywordIndices( RimEclipseCase* eclipseCase, + const RigSimulationInputSettings& settings, + RifOpmFlowDeckFile& deckFile ) +{ + return RigSimulationInputTool::replaceKeywordWithBoxIndices( "AQUANCON", eclipseCase, settings, deckFile, processAquanconRecord ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::expected RigSimulationInputAquiferTools::replaceAqunumKeywordIndices( RimEclipseCase* eclipseCase, + const RigSimulationInputSettings& settings, + RifOpmFlowDeckFile& deckFile ) +{ + return RigSimulationInputTool::replaceKeywordWithBoxIndices( "AQUNUM", eclipseCase, settings, deckFile, processAqunumRecord ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::expected 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 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 inputBox( origMin.toZeroBased(), origMax.toZeroBased() ); + + // Get aquifer ID for logging + int aquiferId = record.getItem( 0 ).get( 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( transformResult->min().x() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::I2::itemName, static_cast( transformResult->max().x() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::J1::itemName, static_cast( transformResult->min().y() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::J2::itemName, static_cast( transformResult->max().y() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::K1::itemName, static_cast( transformResult->min().z() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::K2::itemName, static_cast( 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 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 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 inputBox( origMin.toZeroBased(), origMax.toZeroBased() ); + + // Get aquifer ID for logging + int aquiferId = record.getItem( 0 ).get( 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( transformResult->min().x() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::I2::itemName, static_cast( transformResult->max().x() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::J1::itemName, static_cast( transformResult->min().y() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::J2::itemName, static_cast( transformResult->max().y() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::K1::itemName, static_cast( transformResult->min().z() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::K2::itemName, static_cast( 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 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 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 inputBox( origCell.toZeroBased(), origCell.toZeroBased() ); + + // Get aquifer ID for logging + int aquiferId = record.getItem( 0 ).get( 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( transformResult->min().x() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::J::itemName, static_cast( transformResult->min().y() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( A::K::itemName, static_cast( 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 ) }; +} diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.h b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.h new file mode 100644 index 0000000000..7873159718 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputAquiferTools.h @@ -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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "RigRefinement.h" + +#include "cafVecIjk.h" + +#include + +#include + +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 + replaceAquconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); + + static std::expected + replaceAquanconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); + + static std::expected + replaceAqunumKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); + + static std::expected + processAquconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement ); + + static std::expected + processAquanconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement ); + + static std::expected + processAqunumRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement ); +}; diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.cpp b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.cpp new file mode 100644 index 0000000000..db6697db9f --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.cpp @@ -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 +// 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 +#include +#include +#include + +//-------------------------------------------------------------------------------------------------- +/// 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 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 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( 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( const Opm::DeckRecord& record, + const std::string& wellName, + RigRefinement::Dimension penetrationDim, + double fraction1, + double fraction2, + const std::optional& 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( 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( 0 ); + } + + // Peaceman radius of the original and refined cell, using the cell sizes and permeabilities perpendicular to the well + std::optional> 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( 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*, 3> permeabilities = { nullptr, nullptr, nullptr }; + if ( results && activeCellInfo ) + { + const std::array 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 + { + 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& 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( 0 ); + const int i = record.getItem( 1 ).get( 0 ) - 1; + const int j = record.getItem( 2 ).get( 0 ) - 1; + const int k1 = record.getItem( 3 ).get( 0 ) - 1; + const int k2 = record.getItem( 4 ).get( 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& 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( 0 ); + + std::map> 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( 0 ) - 1; + const int j = record.getItem( 1 ).get( 0 ) - 1; + const int k = record.getItem( 2 ).get( 0 ) - 1; + if ( i < 0 || j < 0 || k < 0 ) continue; + + const int branch = record.getItem( 3 ).hasValue( 0 ) ? record.getItem( 3 ).get( 0 ) : 0; + const double startMd = record.getItem( 4 ).hasValue( 0 ) ? record.getItem( 4 ).get( 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 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, 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( 0 ) - 1; // Convert to 0-based + int origJ = record.getItem( 2 ).get( 0 ) - 1; + int origK1 = record.getItem( 3 ).get( 0 ) - 1; + int origK2 = record.getItem( 4 ).get( 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 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 iSubcells = + ( dim == RigRefinement::DimI ) + ? allSubcells( refinement, RigRefinement::DimI, sectorI ) + : std::vector{ { centerSubcellIndex( refinement, RigRefinement::DimI, sectorI ), 1.0 } }; + const std::vector jSubcells = + ( dim == RigRefinement::DimJ ) + ? allSubcells( refinement, RigRefinement::DimJ, sectorJ ) + : std::vector{ { 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 + { + 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 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& correction ) + { + std::vector items; + items.push_back( record.getItem( 0 ) ); + items.push_back( RifOpmDeckTools::item( C::I::itemName, static_cast( i + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( C::J::itemName, static_cast( j + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( C::K1::itemName, static_cast( k1 + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( C::K2::itemName, static_cast( 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( 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( 0 ) * fraction ) ); + } + else + { + items.push_back( item ); + } + } + + return Opm::DeckRecord{ std::move( items ) }; + }; + + std::vector 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, 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( 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 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( 0 ) : 0.0; + const double endMd = hasEndMd ? record.getItem( C::DISTANCE_END::itemName ).get( 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 records; + + double startFraction = 0.0; + for ( const auto& subcell : subcells ) + { + caf::VecIjk0 refinedIjk = *transformResult; + refinedIjk[dim] = subcell.index; + + std::vector items; + items.push_back( RifOpmDeckTools::item( "I", static_cast( refinedIjk.x() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( "J", static_cast( refinedIjk.y() + 1 ) ) ); + items.push_back( RifOpmDeckTools::item( "K", static_cast( 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; +} diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.h b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.h new file mode 100644 index 0000000000..856b21f98b --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputConnectionTools.h @@ -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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "RigRefinement.h" + +#include "cafVecIjk.h" +#include "cvfVector3.h" + +#include + +#include +#include +#include +#include +#include +#include +#include + +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; + using ConnectionDirectionMap = std::map; + + static ConnectionDirectionMap connectionDirections( const std::vector& 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& 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( 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, 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, QString> processCompsegsRecord( const Opm::DeckRecord& record, + const std::string& wellName, + bool isWellNameRecord, + const RigSimulationInputSettings& settings, + const ConnectionDirectionMap& directions, + const ConnectionKeySet& reversedCells = {} ); +}; diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.cpp b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.cpp new file mode 100644 index 0000000000..cbc665377b --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.cpp @@ -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 +// 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::extractDeckEditNncConnections( RifOpmFlowDeckFile& deckFile, const RigMainGrid& mainGrid ) +{ + std::vector 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( 0 ), record.getItem( 1 ).get( 0 ), record.getItem( 2 ).get( 0 ) ); + caf::VecIjk1 ijk2( record.getItem( 3 ).get( 0 ), record.getItem( 4 ).get( 0 ), record.getItem( 5 ).get( 0 ) ); + double transMult = record.getItem( 6 ).get( 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::filterInternalSectorConnections( const std::vector& allConnections, + const RigMainGrid& mainGrid, + const caf::VecIjk0& min, + const caf::VecIjk0& max ) +{ + std::vector 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::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::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 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 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 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 {}; +} diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.h b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.h new file mode 100644 index 0000000000..bb32e1e293 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputNncTools.h @@ -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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "RigRefinement.h" + +#include "cafVecIjk.h" + +#include + +#include +#include + +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 + exportEditNncKeyword( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); + + static std::vector extractDeckEditNncConnections( RifOpmFlowDeckFile& deckFile, const RigMainGrid& mainGrid ); + static std::vector filterInternalSectorConnections( const std::vector& allConnections, + const RigMainGrid& mainGrid, + const caf::VecIjk0& min, + const caf::VecIjk0& max ); + static std::expected transformNNCToSectorCoordinates( const NNCConnection& connection, + const RigMainGrid& mainGrid, + const caf::VecIjk0& min, + const RigRefinement& refinement ); + static std::vector refineEditNncConnection( const NNCConnection& connection, + const RigMainGrid& mainGrid, + const caf::VecIjk0& min, + const RigRefinement& refinement ); +}; diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp index ec25b79c03..8ba6782104 100644 --- a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.cpp @@ -29,9 +29,7 @@ #include "ProjectDataModel/Jobs/RimKeywordFactory.h" #include "RimEclipseCase.h" -#include "RigActiveCellInfo.h" #include "RigBoundingBoxIjk.h" -#include "RigCaseCellResultsData.h" #include "RigEclipseCaseData.h" #include "RigEclipseResultAddress.h" #include "RigEclipseResultTools.h" @@ -41,8 +39,10 @@ #include "RigPadModel.h" #include "RigRefinement.h" #include "RigResdataGridConverter.h" +#include "RigSimulationInputAquiferTools.h" +#include "RigSimulationInputConnectionTools.h" +#include "RigSimulationInputNncTools.h" #include "RigSimulationInputSettings.h" -#include "RigTransmissibilityEquations.h" #include "Well/RigSimWellData.h" #include "Well/RigWellResultFrame.h" #include "Well/RigWellResultPoint.h" @@ -66,7 +66,8 @@ #include #include - +#include +#include #include //-------------------------------------------------------------------------------------------------- @@ -129,17 +130,17 @@ std::expected RigSimulationInputTool::exportSimulationInput( RimE return result; } - if ( auto result = replaceAquconKeywordIndices( &eclipseCase, settings, deckFile ); !result ) + if ( auto result = RigSimulationInputAquiferTools::replaceAquconKeywordIndices( &eclipseCase, settings, deckFile ); !result ) { return result; } - if ( auto result = replaceAquanconKeywordIndices( &eclipseCase, settings, deckFile ); !result ) + if ( auto result = RigSimulationInputAquiferTools::replaceAquanconKeywordIndices( &eclipseCase, settings, deckFile ); !result ) { return result; } - if ( auto result = replaceAqunumKeywordIndices( &eclipseCase, settings, deckFile ); !result ) + if ( auto result = RigSimulationInputAquiferTools::replaceAqunumKeywordIndices( &eclipseCase, settings, deckFile ); !result ) { return result; } @@ -171,7 +172,7 @@ std::expected RigSimulationInputTool::exportSimulationInput( RimE return result; } - if ( auto result = exportEditNncKeyword( &eclipseCase, settings, deckFile ); !result ) + if ( auto result = RigSimulationInputNncTools::exportEditNncKeyword( &eclipseCase, settings, deckFile ); !result ) { return result; } @@ -1065,36 +1066,6 @@ std::expected RigSimulationInputTool::replaceKeywordWithBoxIndice return {}; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::expected RigSimulationInputTool::replaceAquconKeywordIndices( RimEclipseCase* eclipseCase, - const RigSimulationInputSettings& settings, - RifOpmFlowDeckFile& deckFile ) -{ - return replaceKeywordWithBoxIndices( "AQUCON", eclipseCase, settings, deckFile, processAquconRecord ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::expected RigSimulationInputTool::replaceAquanconKeywordIndices( RimEclipseCase* eclipseCase, - const RigSimulationInputSettings& settings, - RifOpmFlowDeckFile& deckFile ) -{ - return replaceKeywordWithBoxIndices( "AQUANCON", eclipseCase, settings, deckFile, processAquanconRecord ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::expected RigSimulationInputTool::replaceAqunumKeywordIndices( RimEclipseCase* eclipseCase, - const RigSimulationInputSettings& settings, - RifOpmFlowDeckFile& deckFile ) -{ - return replaceKeywordWithBoxIndices( "AQUNUM", eclipseCase, settings, deckFile, processAqunumRecord ); -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -1271,701 +1242,6 @@ std::expected RigSimulationInputTool::processWelspecsR return Opm::DeckRecord{ std::move( items ) }; } -//-------------------------------------------------------------------------------------------------- -/// 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 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 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( 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( const Opm::DeckRecord& record, - const std::string& wellName, - RigRefinement::Dimension penetrationDim, - double fraction1, - double fraction2, - const std::optional& 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( 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( 0 ); - } - - // Peaceman radius of the original and refined cell, using the cell sizes and permeabilities perpendicular to the well - std::optional> 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( 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 }; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -RigSimulationInputTool::ConnectionCellPropertiesProvider RigSimulationInputTool::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*, 3> permeabilities = { nullptr, nullptr, nullptr }; - if ( results && activeCellInfo ) - { - const std::array 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 - { - 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 RigSimulationInputTool::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() ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -RigSimulationInputTool::ConnectionDirectionMap RigSimulationInputTool::connectionDirections( const std::vector& 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( 0 ); - const int i = record.getItem( 1 ).get( 0 ) - 1; - const int j = record.getItem( 2 ).get( 0 ) - 1; - const int k1 = record.getItem( 3 ).get( 0 ) - 1; - const int k2 = record.getItem( 4 ).get( 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 -//-------------------------------------------------------------------------------------------------- -RigSimulationInputTool::ConnectionKeySet RigSimulationInputTool::reversedSegmentCells( const std::vector& 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( 0 ); - - std::map> 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( 0 ) - 1; - const int j = record.getItem( 1 ).get( 0 ) - 1; - const int k = record.getItem( 2 ).get( 0 ) - 1; - if ( i < 0 || j < 0 || k < 0 ) continue; - - const int branch = record.getItem( 3 ).hasValue( 0 ) ? record.getItem( 3 ).get( 0 ) : 0; - const double startMd = record.getItem( 4 ).hasValue( 0 ) ? record.getItem( 4 ).get( 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 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, QString> - RigSimulationInputTool::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( 0 ) - 1; // Convert to 0-based - int origJ = record.getItem( 2 ).get( 0 ) - 1; - int origK1 = record.getItem( 3 ).get( 0 ) - 1; - int origK2 = record.getItem( 4 ).get( 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 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 iSubcells = - ( dim == RigRefinement::DimI ) - ? allSubcells( refinement, RigRefinement::DimI, sectorI ) - : std::vector{ { centerSubcellIndex( refinement, RigRefinement::DimI, sectorI ), 1.0 } }; - const std::vector jSubcells = - ( dim == RigRefinement::DimJ ) - ? allSubcells( refinement, RigRefinement::DimJ, sectorJ ) - : std::vector{ { 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 - { - 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 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& correction ) - { - std::vector items; - items.push_back( record.getItem( 0 ) ); - items.push_back( RifOpmDeckTools::item( C::I::itemName, static_cast( i + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( C::J::itemName, static_cast( j + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( C::K1::itemName, static_cast( k1 + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( C::K2::itemName, static_cast( 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( 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( 0 ) * fraction ) ); - } - else - { - items.push_back( item ); - } - } - - return Opm::DeckRecord{ std::move( items ) }; - }; - - std::vector 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, QString> RigSimulationInputTool::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( 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 = 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 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( 0 ) : 0.0; - const double endMd = hasEndMd ? record.getItem( C::DISTANCE_END::itemName ).get( 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 records; - - double startFraction = 0.0; - for ( const auto& subcell : subcells ) - { - caf::VecIjk0 refinedIjk = *transformResult; - refinedIjk[dim] = subcell.index; - - std::vector items; - items.push_back( RifOpmDeckTools::item( "I", static_cast( refinedIjk.x() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( "J", static_cast( refinedIjk.y() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( "K", static_cast( 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; -} - //-------------------------------------------------------------------------------------------------- /// Helper function to extract IJK coordinates from a deck record /// Returns a VecIjk1 (1-based) constructed from values at the specified item indices @@ -2222,211 +1498,6 @@ std::expected RigSimulationInputTool::processAddRecord return Opm::DeckRecord{ std::move( items ) }; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::expected RigSimulationInputTool::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 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 = extractIjk( record, 1, 3, 5 ); // I1, J1, K1 - caf::VecIjk1 origMax = extractIjk( record, 2, 4, 6 ); // I2, J2, K2 - - // Create input bounding box (0-based, inclusive) - RigBoundingBoxIjk inputBox( origMin.toZeroBased(), origMax.toZeroBased() ); - - // Get aquifer ID for logging - int aquiferId = record.getItem( 0 ).get( 0 ); - QString recordIdent = QString( "for aquifer %1" ).arg( aquiferId ); - - // Transform box to sector coordinates - auto transformResult = - 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( transformResult->min().x() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::I2::itemName, static_cast( transformResult->max().x() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::J1::itemName, static_cast( transformResult->min().y() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::J2::itemName, static_cast( transformResult->max().y() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::K1::itemName, static_cast( transformResult->min().z() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::K2::itemName, static_cast( 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 RigSimulationInputTool::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 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 = extractIjk( record, 1, 3, 5 ); // I1, J1, K1 - caf::VecIjk1 origMax = extractIjk( record, 2, 4, 6 ); // I2, J2, K2 - - // Create input bounding box (0-based, inclusive) - RigBoundingBoxIjk inputBox( origMin.toZeroBased(), origMax.toZeroBased() ); - - // Get aquifer ID for logging - int aquiferId = record.getItem( 0 ).get( 0 ); - QString recordIdent = QString( "for aquifer %1" ).arg( aquiferId ); - - // Transform box to sector coordinates - auto transformResult = - 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( transformResult->min().x() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::I2::itemName, static_cast( transformResult->max().x() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::J1::itemName, static_cast( transformResult->min().y() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::J2::itemName, static_cast( transformResult->max().y() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::K1::itemName, static_cast( transformResult->min().z() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::K2::itemName, static_cast( 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 RigSimulationInputTool::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 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 = extractIjk( record, 1, 2, 3 ); // I, J, K - - // Create input bounding box (0-based, inclusive) - RigBoundingBoxIjk inputBox( origCell.toZeroBased(), origCell.toZeroBased() ); - - // Get aquifer ID for logging - int aquiferId = record.getItem( 0 ).get( 0 ); - QString recordIdent = QString( "for aquifer %1" ).arg( aquiferId ); - - // Transform box to sector coordinates - auto transformResult = - 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( transformResult->min().x() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::J::itemName, static_cast( transformResult->min().y() + 1 ) ) ); - items.push_back( RifOpmDeckTools::item( A::K::itemName, static_cast( 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 ) }; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -2702,10 +1773,10 @@ std::expected RigSimulationInputTool::filterAndUpdateWellKeywords "WTEMP", "WTEST", "WTRACER", "WSEGVALV", "WSEGAICD", "WSEGSICD" }; // Collect connection directions from the original COMPDAT before it is transformed, used to split COMPSEGS - const auto directions = connectionDirections( deckFile.findAllKeywords( "COMPDAT" ) ); - const auto reversedCells = reversedSegmentCells( deckFile.findAllKeywords( "COMPSEGS" ), directions ); + const auto directions = RigSimulationInputConnectionTools::connectionDirections( deckFile.findAllKeywords( "COMPDAT" ) ); + const auto reversedCells = RigSimulationInputConnectionTools::reversedSegmentCells( deckFile.findAllKeywords( "COMPSEGS" ), directions ); - const auto cellProperties = connectionCellPropertiesProvider( eclipseCase ); + const auto cellProperties = RigSimulationInputConnectionTools::connectionCellPropertiesProvider( eclipseCase ); // Process each type of well keyword // Use findAllKeywordsWithIndices to get all occurrences with their positions @@ -2784,7 +1855,7 @@ std::expected RigSimulationInputTool::filterAndUpdateWellKeywords } else if ( keywordName == "COMPDAT" ) { - auto result = processCompdatRecord( record, wellName, settings, cellProperties ); + auto result = RigSimulationInputConnectionTools::processCompdatRecord( record, wellName, settings, cellProperties ); if ( result ) { for ( auto& rec : result.value() ) @@ -2800,7 +1871,12 @@ std::expected RigSimulationInputTool::filterAndUpdateWellKeywords } else if ( keywordName == "COMPSEGS" ) { - auto result = processCompsegsRecord( record, wellName, isWellNameRecord, settings, directions, reversedCells ); + auto result = RigSimulationInputConnectionTools::processCompsegsRecord( record, + wellName, + isWellNameRecord, + settings, + directions, + reversedCells ); if ( result ) { for ( auto& rec : result.value() ) @@ -2939,78 +2015,6 @@ std::expected RigSimulationInputTool::addOperNumRegionAndOperater return {}; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RigSimulationInputTool::extractDeckEditNncConnections( RifOpmFlowDeckFile& deckFile, - const RigMainGrid& mainGrid ) -{ - std::vector 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( 0 ), record.getItem( 1 ).get( 0 ), record.getItem( 2 ).get( 0 ) ); - caf::VecIjk1 ijk2( record.getItem( 3 ).get( 0 ), record.getItem( 4 ).get( 0 ), record.getItem( 5 ).get( 0 ) ); - double transMult = record.getItem( 6 ).get( 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 - RigSimulationInputTool::filterInternalSectorConnections( const std::vector& allConnections, - const RigMainGrid& mainGrid, - const caf::VecIjk0& min, - const caf::VecIjk0& max ) -{ - std::vector 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; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -3029,176 +2033,6 @@ caf::VecIjk0 RigSimulationInputTool::transformToSectorCoordinates( const caf::Ve refinement.cumulativeOffset( RigRefinement::DimK, sectorK ) ); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::expected - RigSimulationInputTool::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 = transformToSectorCoordinates( *gijk1, min, refinement ); - transformed.cell2 = transformToSectorCoordinates( *gijk2, min, refinement ); - transformed.transmissibility = connection.transmissibility; - - return transformed; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RigSimulationInputTool::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 = transformToSectorCoordinates( *gijk1, min, refinement ); - caf::VecIjk0 baseCell2 = 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 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 RigSimulationInputTool::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 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 {}; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.h b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.h index 5bbcc1ded8..1ec7472b00 100644 --- a/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.h +++ b/ApplicationLibCode/ReservoirDataModel/SimulationFile/RigSimulationInputTool.h @@ -28,11 +28,8 @@ #include #include -#include -#include #include #include -#include #include 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 exportSimulationInput( RimEclipseCase& eclipseCase, const RigSimulationInputSettings& settings, cvf::UByteArray* visibility ); @@ -96,86 +75,38 @@ public: static std::expected processAddRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement ); - static std::expected - processAquconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement ); - - static std::expected - processAquanconRecord( const Opm::DeckRecord& record, const caf::VecIjk0& min, const caf::VecIjk0& max, const RigRefinement& refinement ); - - static std::expected - 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; - using ConnectionDirectionMap = std::map; - - static ConnectionDirectionMap connectionDirections( const std::vector& 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& 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( 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, 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, QString> processCompsegsRecord( const Opm::DeckRecord& record, - const std::string& wellName, - bool isWellNameRecord, - const RigSimulationInputSettings& settings, - const ConnectionDirectionMap& directions, - const ConnectionKeySet& reversedCells = {} ); - static std::vector createRefinedVisibility( const RigGridExportAdapter& gridAdapter ); static std::expected transformKeywordsInDeckFile( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); - // Public EDITNNC methods for testing - static std::vector extractDeckEditNncConnections( RifOpmFlowDeckFile& deckFile, const RigMainGrid& mainGrid ); - static std::vector filterInternalSectorConnections( const std::vector& allConnections, - const RigMainGrid& mainGrid, - const caf::VecIjk0& min, - const caf::VecIjk0& max ); - static std::expected transformNNCToSectorCoordinates( const NNCConnection& connection, - const RigMainGrid& mainGrid, - const caf::VecIjk0& min, - const RigRefinement& refinement ); - static std::vector 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( const Opm::DeckRecord&, const caf::VecIjk0&, const caf::VecIjk0&, const RigRefinement& )>; + static std::expected 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, QString> + transformBoxToSectorCoordinates( const RigBoundingBoxIjk& 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 updateCornerPointGridInDeckFile( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, @@ -202,15 +133,6 @@ private: cvf::ref visibility, RifOpmFlowDeckFile& deckFile ); - static std::expected - replaceAquconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); - - static std::expected - replaceAquanconKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); - - static std::expected - replaceAqunumKeywordIndices( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); - static std::expected addFaultsToDeckFile( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); @@ -225,14 +147,8 @@ private: static std::set wellNamesToInclude( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings ); - static std::expected - exportEditNncKeyword( RimEclipseCase* eclipseCase, const RigSimulationInputSettings& settings, RifOpmFlowDeckFile& deckFile ); - static std::expected 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 addOperNumRegionAndOperater( RifOpmFlowDeckFile& deckFile, const RigGridExportAdapter& gridAdapter, const std::vector& operNumResult, @@ -243,24 +159,4 @@ private: static std::expected processWelspecsRecord( const Opm::DeckRecord& record, const std::string& wellName, const RigSimulationInputSettings& settings ); - - static std::expected 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, QString> - transformBoxToSectorCoordinates( const RigBoundingBoxIjk& 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 ); }; diff --git a/ApplicationLibCode/UnitTests/CMakeLists.txt b/ApplicationLibCode/UnitTests/CMakeLists.txt index 89c27298df..c4613fa2f7 100644 --- a/ApplicationLibCode/UnitTests/CMakeLists.txt +++ b/ApplicationLibCode/UnitTests/CMakeLists.txt @@ -10,6 +10,9 @@ set(SOURCE_UNITTEST_FILES ${CMAKE_CURRENT_LIST_DIR}/Ert-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigModelPaddingSettings-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigPadModel-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputAquiferTools-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputConnectionTools-Test.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputNncTools-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigSimulationInputTool-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigBoundingBoxIjk-Test.cpp ${CMAKE_CURRENT_LIST_DIR}/RigContourMapGrid-Test.cpp diff --git a/ApplicationLibCode/UnitTests/RigSimulationInputAquiferTools-Test.cpp b/ApplicationLibCode/UnitTests/RigSimulationInputAquiferTools-Test.cpp new file mode 100644 index 0000000000..11564add0c --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigSimulationInputAquiferTools-Test.cpp @@ -0,0 +1,153 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RigSimulationInputTestHelpers.h" + +#include "RifOpmDeckTools.h" + +#include "RigSimulationInputAquiferTools.h" + +#include "opm/input/eclipse/Deck/DeckItem.hpp" +#include "opm/input/eclipse/Deck/DeckRecord.hpp" + +#include +#include + +using RigSimulationInputTestHelpers::makeRef; + +//-------------------------------------------------------------------------------------------------- +/// Helper to create an AQUCON DeckRecord +/// AQUCON format: ID I1 I2 J1 J2 K1 K2 CONNECT_FACE ... (1-based Eclipse coordinates) +//-------------------------------------------------------------------------------------------------- +static Opm::DeckRecord createAquconRecord( int aquiferId, int i1, int i2, int j1, int j2, int k1, int k2, const std::string& connectFace ) +{ + std::vector items; + items.push_back( RifOpmDeckTools::item( "ID", aquiferId ) ); + items.push_back( RifOpmDeckTools::item( "I1", i1 ) ); + items.push_back( RifOpmDeckTools::item( "I2", i2 ) ); + items.push_back( RifOpmDeckTools::item( "J1", j1 ) ); + items.push_back( RifOpmDeckTools::item( "J2", j2 ) ); + items.push_back( RifOpmDeckTools::item( "K1", k1 ) ); + items.push_back( RifOpmDeckTools::item( "K2", k2 ) ); + items.push_back( RifOpmDeckTools::item( "CONNECT_FACE", connectFace ) ); + + return Opm::DeckRecord{ std::move( items ) }; +} + +//-------------------------------------------------------------------------------------------------- +/// Test AQUCON record with partial overlap (based on user-provided data) +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputAquiferTools, ProcessAquconRecord_PartialOverlap ) +{ + // Sector: min=(40, 27, 21), max=(100, 72, 58) (inclusive, 0-based) + // This corresponds to a sector with cells [40-100, 27-72, 21-58] + caf::VecIjk0 sectorMin( 40, 27, 21 ); + caf::VecIjk0 sectorMax( 100, 72, 58 ); + auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); + + // AQUCON record from user data: 1 57 57 28 36 46 58 'I+' (1-based Eclipse) + // Converts to 0-based: I[56-56], J[27-35], K[45-57] + // Sector is I[40-100], J[27-72], K[21-58] + // Overlap: I[56-56], J[27-35], K[45-57] - all inside, no clamping needed + auto record = createAquconRecord( 1, 57, 57, 28, 36, 46, 58, "I+" ); + + auto result = RigSimulationInputAquiferTools::processAquconRecord( record, sectorMin, sectorMax, refinement ); + + ASSERT_TRUE( result.has_value() ); + + // Check transformed coordinates (sector-relative, 1-based) + // I: (56-40)*1+1 = 17 to (56-40)*1+1 = 17 + // J: (27-27)*1+1 = 1 to (35-27)*1+1 = 9 + // K: (45-21)*1+1 = 25 to (57-21)*1+1 = 37 + EXPECT_EQ( 17, result->getItem( 1 ).get( 0 ) ); // I1 + EXPECT_EQ( 17, result->getItem( 2 ).get( 0 ) ); // I2 + EXPECT_EQ( 1, result->getItem( 3 ).get( 0 ) ); // J1 + EXPECT_EQ( 9, result->getItem( 4 ).get( 0 ) ); // J2 + EXPECT_EQ( 25, result->getItem( 5 ).get( 0 ) ); // K1 + EXPECT_EQ( 37, result->getItem( 6 ).get( 0 ) ); // K2 + + // Check aquifer ID and connect face are preserved + EXPECT_EQ( 1, result->getItem( 0 ).get( 0 ) ); + EXPECT_EQ( "I+", result->getItem( 7 ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test AQUCON record with partial K overlap requiring clamping +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputAquiferTools, ProcessAquconRecord_PartialOverlapWithClamping ) +{ + // Sector: min=(40, 27, 21), max=(100, 72, 58) (inclusive, 0-based) + caf::VecIjk0 sectorMin( 40, 27, 21 ); + caf::VecIjk0 sectorMax( 100, 72, 58 ); + auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); + + // AQUCON record: 1 79 79 41 67 5 11 'I+' (1-based Eclipse) + // Converts to 0-based: I[78-78], J[40-66], K[4-10] + // Sector K is [21-58], so K[4-10] is outside and should be skipped + auto record = createAquconRecord( 1, 79, 79, 41, 67, 5, 11, "I+" ); + + auto result = RigSimulationInputAquiferTools::processAquconRecord( record, sectorMin, sectorMax, refinement ); + + // This should fail because K range [4-10] doesn't overlap with sector K range [21-58] + EXPECT_FALSE( result.has_value() ); + EXPECT_TRUE( result.error().contains( "does not overlap" ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test AQUCON record completely inside sector +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputAquiferTools, ProcessAquconRecord_CompletelyInside ) +{ + // Sector: min=(40, 27, 21), max=(100, 72, 58) (inclusive, 0-based) + caf::VecIjk0 sectorMin( 40, 27, 21 ); + caf::VecIjk0 sectorMax( 100, 72, 58 ); + auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); + + // AQUCON record: 1 61 61 48 72 12 17 'I+' (1-based Eclipse) + // Converts to 0-based: I[60-60], J[47-71], K[11-16] + // However, K[11-16] is outside sector K[21-58], so this should fail + auto record = createAquconRecord( 1, 61, 61, 48, 72, 12, 17, "I+" ); + + auto result = RigSimulationInputAquiferTools::processAquconRecord( record, sectorMin, sectorMax, refinement ); + + EXPECT_FALSE( result.has_value() ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test AQUCON record without box definition (only aquifer ID) +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputAquiferTools, ProcessAquconRecord_NoBoxDefinition ) +{ + caf::VecIjk0 sectorMin( 0, 0, 0 ); + caf::VecIjk0 sectorMax( 9, 9, 9 ); + auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); + + // Create a record with only aquifer ID (no box coordinates) + std::vector items; + items.push_back( RifOpmDeckTools::item( "ID", 1 ) ); + + Opm::DeckRecord record{ std::move( items ) }; + + auto result = RigSimulationInputAquiferTools::processAquconRecord( record, sectorMin, sectorMax, refinement ); + + EXPECT_TRUE( result.has_value() ); + EXPECT_EQ( result->size(), 1U ); + EXPECT_EQ( 1, result->getItem( 0 ).get( 0 ) ); +} diff --git a/ApplicationLibCode/UnitTests/RigSimulationInputConnectionTools-Test.cpp b/ApplicationLibCode/UnitTests/RigSimulationInputConnectionTools-Test.cpp new file mode 100644 index 0000000000..5ba93d009b --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigSimulationInputConnectionTools-Test.cpp @@ -0,0 +1,408 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RigSimulationInputTestHelpers.h" + +#include "RigNonUniformRefinement.h" +#include "RigSimulationInputConnectionTools.h" +#include "RigSimulationInputSettings.h" +#include "RigTransmissibilityEquations.h" + +#include "opm/input/eclipse/Deck/Deck.hpp" +#include "opm/input/eclipse/Deck/DeckItem.hpp" +#include "opm/input/eclipse/Deck/DeckRecord.hpp" +#include "opm/input/eclipse/Parser/Parser.hpp" + +#include +#include +#include + +using RigSimulationInputTestHelpers::setSettingsRefinement; + +//-------------------------------------------------------------------------------------------------- +/// Helper: parse a keyword from deck text and return all its records +//-------------------------------------------------------------------------------------------------- +static std::vector parseKeywordRecords( const std::string& deckText, const std::string& keywordName ) +{ + Opm::Parser parser; + auto deck = parser.parseString( deckText ); + const auto& keyword = deck[keywordName].back(); + + std::vector records; + for ( size_t i = 0; i < keyword.size(); ++i ) + { + records.push_back( keyword.getRecord( i ) ); + } + return records; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +static RigSimulationInputSettings createSectorSettings( const cvf::Vec3st& refCounts ) +{ + RigSimulationInputSettings settings; + settings.setMin( caf::VecIjk0( 0, 0, 0 ) ); + settings.setMax( caf::VecIjk0( 9, 9, 9 ) ); + setSettingsRefinement( settings, refCounts ); + return settings; +} + +//-------------------------------------------------------------------------------------------------- +/// A Z-direction connection in a K-refined layer must be completed in every sub-cell, with explicit +/// CF and KH split between the sub-cells (regression: only every other refined layer was completed) +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_ZDirectionSplitAcrossKSubcells ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); + auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 100.0 2* 'Z' /\n/\n", "COMPDAT" ); + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + + for ( size_t s = 0; s < 2; ++s ) + { + const auto& rec = result->at( s ); + EXPECT_EQ( 3, rec.getItem( "I" ).get( 0 ) ); + EXPECT_EQ( 4, rec.getItem( "J" ).get( 0 ) ); + EXPECT_EQ( 3 + static_cast( s ), rec.getItem( "K1" ).get( 0 ) ); + EXPECT_EQ( 3 + static_cast( s ), rec.getItem( "K2" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 5.0, rec.getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 50.0, rec.getItem( "Kh" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 0.216, rec.getItem( "DIAMETER" ).get( 0 ) ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// With defaulted CF and KH the simulator computes the factors, so the full refined K range is kept in one record +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_DefaultedFactorsKeepsSingleRecord ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); + auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 3 'OPEN' 2* 0.216 /\n/\n", "COMPDAT" ); + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 1u, result->size() ); + EXPECT_EQ( 3, result->at( 0 ).getItem( "K1" ).get( 0 ) ); + EXPECT_EQ( 6, result->at( 0 ).getItem( "K2" ).get( 0 ) ); + EXPECT_FALSE( result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).hasValue( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// An X-direction connection is split across the I sub-cells; J and K use the center sub-cell +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_XDirectionSplitAcrossISubcells ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 3, 2, 2 ) ); + auto records = parseKeywordRecords( "COMPDAT\n 'W' 2 1 1 1 'OPEN' 1* 9.0 0.216 90.0 2* 'X' /\n/\n", "COMPDAT" ); + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 3u, result->size() ); + + for ( size_t s = 0; s < 3; ++s ) + { + const auto& rec = result->at( s ); + EXPECT_EQ( 4 + static_cast( s ), rec.getItem( "I" ).get( 0 ) ); + EXPECT_EQ( 2, rec.getItem( "J" ).get( 0 ) ); + EXPECT_EQ( 2, rec.getItem( "K1" ).get( 0 ) ); + EXPECT_EQ( 2, rec.getItem( "K2" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 3.0, rec.getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 30.0, rec.getItem( "Kh" ).get( 0 ) ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// Non-uniform K refinement splits the CF by the sub-cell thickness fractions +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_NonUniformKRefinementScalesByFraction ) +{ + RigSimulationInputSettings settings; + settings.setMin( caf::VecIjk0( 0, 0, 0 ) ); + settings.setMax( caf::VecIjk0( 9, 9, 9 ) ); + auto refinement = std::make_unique( cvf::Vec3st( 10, 10, 10 ) ); + refinement->setCumulativeFractions( RigRefinement::DimK, 1, { 0.25, 1.0 } ); + settings.setRefinement( std::move( refinement ) ); + + auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 8.0 0.216 /\n/\n", "COMPDAT" ); + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + EXPECT_EQ( 2, result->at( 0 ).getItem( "K1" ).get( 0 ) ); + EXPECT_EQ( 3, result->at( 1 ).getItem( "K1" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 2.0, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 6.0, result->at( 1 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); + EXPECT_TRUE( result->at( 0 ).getItem( "Kh" ).defaultApplied( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// Without refinement the record is passed through unchanged +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_NoRefinementUnchanged ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 1 ) ); + auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 5 'OPEN' 1* 10.0 0.216 100.0 2* 'Z' /\n/\n", "COMPDAT" ); + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 1u, result->size() ); + EXPECT_EQ( 2, result->at( 0 ).getItem( "K1" ).get( 0 ) ); + EXPECT_EQ( 5, result->at( 0 ).getItem( "K2" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 10.0, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// COMPSEGS cells are split along the connection direction with the MD interval divided between sub-cells +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompsegsRecord_SplitAcrossKSubcells ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); + + Opm::Parser parser; + auto deck = parser.parseString( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 /\n/\n" + "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 120.0 /\n/\n" ); + + std::vector compdatKeywords = { deck["COMPDAT"].back() }; + auto directions = RigSimulationInputConnectionTools::connectionDirections( compdatKeywords ); + + const auto& compsegs = deck["COMPSEGS"].back(); + auto result = RigSimulationInputConnectionTools::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, directions ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + + EXPECT_EQ( 3, result->at( 0 ).getItem( 2 ).get( 0 ) ); + EXPECT_EQ( 4, result->at( 1 ).getItem( 2 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 100.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 110.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 110.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 120.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// Refinement perpendicular to a Z well reduces the Peaceman radius given by PR, so CF increases +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_LateralRefinementCorrectsCFUsingPR ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 10, 10, 1 ) ); + auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 100.0 1.0 1* 'Z' 20.0 /\n/\n", "COMPDAT" ); + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 1u, result->size() ); + + const auto& rec = result->at( 0 ); + EXPECT_EQ( 26, rec.getItem( "I" ).get( 0 ) ); + EXPECT_EQ( 36, rec.getItem( "J" ).get( 0 ) ); + + const double rw = 0.108; + const double expected = 10.0 * ( std::log( 20.0 / rw ) + 1.0 ) / ( std::log( 2.0 / rw ) + 1.0 ); + EXPECT_NEAR( expected, rec.getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ), 1e-9 ); + EXPECT_NEAR( 2.0, rec.getItem( "PR" ).get( 0 ), 1e-9 ); + EXPECT_DOUBLE_EQ( 100.0, rec.getItem( "Kh" ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// Without PR the Peaceman radii are computed from the cell sizes and permeabilities, also for +/// refinement in only one of the perpendicular directions +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_LateralRefinementCorrectsCFUsingCellProperties ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 4, 1, 1 ) ); + auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 /\n/\n", "COMPDAT" ); + + RigSimulationInputConnectionTools::ConnectionCellProperties properties{ cvf::Vec3d( 100.0, 50.0, 5.0 ), cvf::Vec3d( 200.0, 50.0, 10.0 ) }; + auto provider = [&]( const caf::VecIjk0& ) -> std::optional + { return properties; }; + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings, provider ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 1u, result->size() ); + + const double rw = 0.108; + const double r0Coarse = RigTransmissibilityEquations::peacemanRadius( 200.0, 50.0, 100.0, 50.0 ); + const double r0Fine = RigTransmissibilityEquations::peacemanRadius( 200.0, 50.0, 25.0, 50.0 ); + const double expected = 10.0 * std::log( r0Coarse / rw ) / std::log( r0Fine / rw ); + EXPECT_NEAR( expected, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ), 1e-9 ); +} + +//-------------------------------------------------------------------------------------------------- +/// The correction is skipped when the refined Peaceman radius gets smaller than the well radius +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompdatRecord_LateralRefinementTooSmallRadiusKeepsCF ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 10, 10, 1 ) ); + auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 4* 0.5 /\n/\n", "COMPDAT" ); + + auto result = RigSimulationInputConnectionTools::processCompdatRecord( records[0], "W", settings ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 1u, result->size() ); + EXPECT_DOUBLE_EQ( 10.0, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 0.5, result->at( 0 ).getItem( "PR" ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// A well travelling towards decreasing I assigns the MD interval from the highest to the lowest I sub-cell +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompsegsRecord_ReversedTraversalAlongI ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 2, 1, 1 ) ); + + Opm::Parser parser; + auto deck = parser.parseString( "COMPDAT\n 'W' 4 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n" + " 'W' 3 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n/\n" + "COMPSEGS\n 'W' /\n 4 2 1 1 100.0 120.0 /\n 3 2 1 1 120.0 140.0 /\n/\n" ); + + std::vector compdatKeywords = { deck["COMPDAT"].back() }; + std::vector compsegsKeywords = { deck["COMPSEGS"].back() }; + + auto directions = RigSimulationInputConnectionTools::connectionDirections( compdatKeywords ); + auto reversedCells = RigSimulationInputConnectionTools::reversedSegmentCells( compsegsKeywords, directions ); + + EXPECT_TRUE( reversedCells.contains( { "W", caf::VecIjk0( 3, 1, 0 ) } ) ); + EXPECT_TRUE( reversedCells.contains( { "W", caf::VecIjk0( 2, 1, 0 ) } ) ); + + const auto& compsegs = deck["COMPSEGS"].back(); + + // Original I=4 -> refined I sub-cells 7 and 8, the well enters at 8 + auto result = + RigSimulationInputConnectionTools::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, directions, reversedCells ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + EXPECT_EQ( 8, result->at( 0 ).getItem( 0 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 100.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 110.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); + EXPECT_EQ( 7, result->at( 1 ).getItem( 0 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 110.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 120.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); + + // Original I=3 is the last cell, direction taken from the previous cell + result = + RigSimulationInputConnectionTools::processCompsegsRecord( compsegs.getRecord( 2 ), "W", false, settings, directions, reversedCells ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + EXPECT_EQ( 6, result->at( 0 ).getItem( 0 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 120.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); + EXPECT_EQ( 5, result->at( 1 ).getItem( 0 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 140.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// A well travelling towards increasing I keeps the sub-cell order +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ReversedSegmentCells_ForwardTraversalNotReversed ) +{ + Opm::Parser parser; + auto deck = parser.parseString( "COMPDAT\n 'W' 3 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n" + " 'W' 4 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n/\n" + "COMPSEGS\n 'W' /\n 3 2 1 1 100.0 120.0 /\n 4 2 1 1 120.0 140.0 /\n/\n" ); + + std::vector compdatKeywords = { deck["COMPDAT"].back() }; + std::vector compsegsKeywords = { deck["COMPSEGS"].back() }; + + auto directions = RigSimulationInputConnectionTools::connectionDirections( compdatKeywords ); + auto reversedCells = RigSimulationInputConnectionTools::reversedSegmentCells( compsegsKeywords, directions ); + + EXPECT_TRUE( reversedCells.empty() ); +} + +//-------------------------------------------------------------------------------------------------- +/// DIRECTION without END_IJK is a single cell, so the MD interval is still split between the sub-cells +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompsegsRecord_DirectionWithoutRangeIsSplit ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); + + Opm::Parser parser; + auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 120.0 'K' /\n/\n" ); + + const auto& compsegs = deck["COMPSEGS"].back(); + auto result = RigSimulationInputConnectionTools::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + EXPECT_EQ( 3, result->at( 0 ).getItem( 2 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 110.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); + EXPECT_EQ( 4, result->at( 1 ).getItem( 2 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 110.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); + EXPECT_EQ( "K", result->at( 1 ).getItem( 6 ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// A defaulted start MD is zero +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompsegsRecord_DefaultedStartMdIsZero ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); + + Opm::Parser parser; + auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 1* 20.0 /\n/\n" ); + + const auto& compsegs = deck["COMPSEGS"].back(); + auto result = RigSimulationInputConnectionTools::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + EXPECT_DOUBLE_EQ( 0.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 10.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 10.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 20.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); +} + +//-------------------------------------------------------------------------------------------------- +/// Without an end MD every sub-cell still gets a segment record, with the MD items unchanged +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompsegsRecord_NoEndMdCopiedToAllSubcells ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); + + Opm::Parser parser; + auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 1* 'K' /\n/\n" ); + + const auto& compsegs = deck["COMPSEGS"].back(); + auto result = RigSimulationInputConnectionTools::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 2u, result->size() ); + for ( size_t s = 0; s < 2; ++s ) + { + EXPECT_EQ( 3 + static_cast( s ), result->at( s ).getItem( 2 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 100.0, result->at( s ).getItem( 4 ).get( 0 ) ); + EXPECT_FALSE( result->at( s ).getItem( 5 ).hasValue( 0 ) ); + } +} + +//-------------------------------------------------------------------------------------------------- +/// A cell range (END_IJK) is not split, only the center sub-cell is used +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputConnectionTools, ProcessCompsegsRecord_CellRangeUsesCenterSubcell ) +{ + auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); + + Opm::Parser parser; + auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 120.0 'K' 3 /\n/\n" ); + + const auto& compsegs = deck["COMPSEGS"].back(); + auto result = RigSimulationInputConnectionTools::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); + ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); + ASSERT_EQ( 1u, result->size() ); + EXPECT_EQ( 4, result->at( 0 ).getItem( 2 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 100.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); + EXPECT_DOUBLE_EQ( 120.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); +} diff --git a/ApplicationLibCode/UnitTests/RigSimulationInputNncTools-Test.cpp b/ApplicationLibCode/UnitTests/RigSimulationInputNncTools-Test.cpp new file mode 100644 index 0000000000..d56d0bd185 --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigSimulationInputNncTools-Test.cpp @@ -0,0 +1,193 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "gtest/gtest.h" + +#include "RiaTestDataDirectory.h" + +#include "RifReaderEclipseOutput.h" + +#include "RigEclipseCaseData.h" +#include "RigMainGrid.h" +#include "RigSimulationInputNncTools.h" +#include "RigUniformRefinement.h" + +#include "RimEclipseResultCase.h" + +#include +#include +#include + +#include +#include + +//-------------------------------------------------------------------------------------------------- +/// Helper to load model5 test case +//-------------------------------------------------------------------------------------------------- +static cvf::ref loadModel5TestCase() +{ + QDir baseFolder( TEST_DATA_DIR ); + bool subFolderExists = baseFolder.cd( "RigSimulationInputTool/model5" ); + if ( !subFolderExists ) return nullptr; + + QString egridFilename( "0_BASE_MODEL5.EGRID" ); + QString egridFilePath = baseFolder.absoluteFilePath( egridFilename ); + if ( !QFile::exists( egridFilePath ) ) return nullptr; + + std::unique_ptr resultCase( new RimEclipseResultCase ); + cvf::ref caseData = new RigEclipseCaseData( resultCase.get() ); + + cvf::ref reader = new RifReaderEclipseOutput; + bool loadResult = reader->open( egridFilePath, caseData.p() ); + if ( !loadResult ) return nullptr; + + return caseData; +} + +//-------------------------------------------------------------------------------------------------- +/// Test filterInternalSectorConnections - both cells inside +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputNncTools, FilterNNCConnections_BothInside ) +{ + // Load model5 test data (20x30x10 grid) + auto caseData = loadModel5TestCase(); + ASSERT_TRUE( caseData.notNull() ); + ASSERT_TRUE( caseData->mainGrid() != nullptr ); + EXPECT_EQ( 20u, caseData->mainGrid()->cellCountI() ); + EXPECT_EQ( 30u, caseData->mainGrid()->cellCountJ() ); + EXPECT_EQ( 10u, caseData->mainGrid()->cellCountK() ); + + RigMainGrid* mainGrid = caseData->mainGrid(); + + // Create test connections using the real grid + std::vector allConnections; + + // Connection 1: (2,3,4) to (3,3,4) - both inside sector [1,5] x [1,5] x [1,5] + size_t c1Idx = mainGrid->cellIndexFromIJK( 2, 3, 4 ); + size_t c2Idx = mainGrid->cellIndexFromIJK( 3, 3, 4 ); + allConnections.push_back( { c1Idx, c2Idx, 1.5 } ); + + // Connection 2: (10,20,8) to (11,20,8) - both outside sector + c1Idx = mainGrid->cellIndexFromIJK( 10, 20, 8 ); + c2Idx = mainGrid->cellIndexFromIJK( 11, 20, 8 ); + allConnections.push_back( { c1Idx, c2Idx, 2.5 } ); + + // Connection 3: (1,2,3) to (15,15,8) - one inside, one outside + c1Idx = mainGrid->cellIndexFromIJK( 1, 2, 3 ); + c2Idx = mainGrid->cellIndexFromIJK( 15, 15, 8 ); + allConnections.push_back( { c1Idx, c2Idx, 3.5 } ); + + // Filter for sector [1,5] x [1,5] x [1,5] + caf::VecIjk0 sectorMin( 1, 1, 1 ); + caf::VecIjk0 sectorMax( 5, 5, 5 ); + + auto filtered = RigSimulationInputNncTools::filterInternalSectorConnections( allConnections, *mainGrid, sectorMin, sectorMax ); + + // Only connection 1 should remain + ASSERT_EQ( 1u, filtered.size() ); + EXPECT_EQ( 1.5, filtered[0].transmissibility ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test transformNNCToSectorCoordinates - no refinement +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputNncTools, TransformNNCToSectorCoordinates_NoRefinement ) +{ + // Load model5 test data (20x30x10 grid) + auto caseData = loadModel5TestCase(); + ASSERT_TRUE( caseData.notNull() ); + RigMainGrid* mainGrid = caseData->mainGrid(); + + // Connection between (10,15,5) and (11,15,5) - I-face neighbors + size_t c1Idx = mainGrid->cellIndexFromIJK( 10, 15, 5 ); + size_t c2Idx = mainGrid->cellIndexFromIJK( 11, 15, 5 ); + + RigSimulationInputNncTools::NNCConnection connection{ c1Idx, c2Idx, 4.5 }; + + // Sector starts at (5,10,2), grid is 20x30x10 + caf::VecIjk0 sectorMin( 5, 10, 2 ); + auto refinement = RigUniformRefinement( cvf::Vec3st( 1, 1, 1 ), cvf::Vec3st( 20, 30, 10 ) ); + + auto result = RigSimulationInputNncTools::transformNNCToSectorCoordinates( connection, *mainGrid, sectorMin, refinement ); + + ASSERT_TRUE( result.has_value() ); + + // Cell 1: (10,15,5) - sector min (5,10,2) = (5,5,3) in 0-based sector coords + EXPECT_EQ( 5, result->cell1.i() ); + EXPECT_EQ( 5, result->cell1.j() ); + EXPECT_EQ( 3, result->cell1.k() ); + + // Cell 2: (11,15,5) - sector min (5,10,2) = (6,5,3) in 0-based sector coords + EXPECT_EQ( 6, result->cell2.i() ); + EXPECT_EQ( 5, result->cell2.j() ); + EXPECT_EQ( 3, result->cell2.k() ); + + EXPECT_EQ( 4.5, result->transmissibility ); +} + +//-------------------------------------------------------------------------------------------------- +/// Test refineEDITNNCConnection - corresponding subcells with refinement (2,2,1) +//-------------------------------------------------------------------------------------------------- +TEST( RigSimulationInputNncTools, RefineEDITNNCConnection_CorrespondingSubcells ) +{ + // Load model5 test data (20x30x10 grid) + auto caseData = loadModel5TestCase(); + ASSERT_TRUE( caseData.notNull() ); + RigMainGrid* mainGrid = caseData->mainGrid(); + + // Connection between cells (3,4,5) and (6,8,7) - non-neighbors + size_t c1Idx = mainGrid->cellIndexFromIJK( 3, 4, 5 ); + size_t c2Idx = mainGrid->cellIndexFromIJK( 6, 8, 7 ); + + RigSimulationInputNncTools::NNCConnection connection{ c1Idx, c2Idx, 2.5 }; // TRAN_MULT = 2.5 + + caf::VecIjk0 sectorMin( 0, 0, 0 ); + auto refinement = RigUniformRefinement( cvf::Vec3st( 2, 2, 1 ), cvf::Vec3st( 20, 30, 10 ) ); // Refine by 2x2x1, grid is 20x30x10 + + auto refined = RigSimulationInputNncTools::refineEditNncConnection( connection, *mainGrid, sectorMin, refinement ); + + // For EDITNNC: creates connections between corresponding subcells + // Should create 2*2*1 = 4 connections (one per subcell pair at same relative position) + ASSERT_EQ( 4u, refined.size() ); + + // All connections should have same TRAN_MULT (no distribution) + for ( const auto& conn : refined ) + { + EXPECT_DOUBLE_EQ( 2.5, conn.transmissibility ); + } + + // Verify corresponding subcell connections + // Cell1 starts at sector coordinates (3*2, 4*2, 5*1) = (6, 8, 5) + // Cell2 starts at sector coordinates (6*2, 8*2, 7*1) = (12, 16, 7) + + // Connection 0: (0,0,0) -> (0,0,0): (6,8,5) -> (12,16,7) + EXPECT_EQ( caf::VecIjk0( 6, 8, 5 ), refined[0].cell1 ); + EXPECT_EQ( caf::VecIjk0( 12, 16, 7 ), refined[0].cell2 ); + + // Connection 1: (0,1,0) -> (0,1,0): (6,9,5) -> (12,17,7) + EXPECT_EQ( caf::VecIjk0( 6, 9, 5 ), refined[1].cell1 ); + EXPECT_EQ( caf::VecIjk0( 12, 17, 7 ), refined[1].cell2 ); + + // Connection 2: (1,0,0) -> (1,0,0): (7,8,5) -> (13,16,7) + EXPECT_EQ( caf::VecIjk0( 7, 8, 5 ), refined[2].cell1 ); + EXPECT_EQ( caf::VecIjk0( 13, 16, 7 ), refined[2].cell2 ); + + // Connection 3: (1,1,0) -> (1,1,0): (7,9,5) -> (13,17,7) + EXPECT_EQ( caf::VecIjk0( 7, 9, 5 ), refined[3].cell1 ); + EXPECT_EQ( caf::VecIjk0( 13, 17, 7 ), refined[3].cell2 ); +} diff --git a/ApplicationLibCode/UnitTests/RigSimulationInputTestHelpers.h b/ApplicationLibCode/UnitTests/RigSimulationInputTestHelpers.h new file mode 100644 index 0000000000..7eb1a8d519 --- /dev/null +++ b/ApplicationLibCode/UnitTests/RigSimulationInputTestHelpers.h @@ -0,0 +1,60 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "RigNoRefinement.h" +#include "RigSimulationInputSettings.h" +#include "RigUniformRefinement.h" + +#include "cafVecIjk.h" +#include "cvfVector3.h" + +#include + +namespace RigSimulationInputTestHelpers +{ +//-------------------------------------------------------------------------------------------------- +/// Uniform refinement for the sector given by min/max +//-------------------------------------------------------------------------------------------------- +inline RigUniformRefinement makeRef( const caf::VecIjk0& sectorMin, const caf::VecIjk0& sectorMax, const cvf::Vec3st& refCounts ) +{ + cvf::Vec3st sectorSize( sectorMax.i() - sectorMin.i() + 1, sectorMax.j() - sectorMin.j() + 1, sectorMax.k() - sectorMin.k() + 1 ); + return RigUniformRefinement( refCounts, sectorSize ); +} + +//-------------------------------------------------------------------------------------------------- +/// Set refinement on settings, computing sector size from min/max. +/// If refCounts is (1,1,1), creates RigNoRefinement; otherwise RigUniformRefinement. +//-------------------------------------------------------------------------------------------------- +inline void setSettingsRefinement( RigSimulationInputSettings& settings, const cvf::Vec3st& refCounts ) +{ + cvf::Vec3st sectorSize( settings.max().i() - settings.min().i() + 1, + settings.max().j() - settings.min().j() + 1, + settings.max().k() - settings.min().k() + 1 ); + + if ( refCounts.x() == 1 && refCounts.y() == 1 && refCounts.z() == 1 ) + { + settings.setRefinement( std::make_unique( sectorSize ) ); + } + else + { + settings.setRefinement( std::make_unique( refCounts, sectorSize ) ); + } +} +} // namespace RigSimulationInputTestHelpers diff --git a/ApplicationLibCode/UnitTests/RigSimulationInputTool-Test.cpp b/ApplicationLibCode/UnitTests/RigSimulationInputTool-Test.cpp index 60255ed81d..b021311bf9 100644 --- a/ApplicationLibCode/UnitTests/RigSimulationInputTool-Test.cpp +++ b/ApplicationLibCode/UnitTests/RigSimulationInputTool-Test.cpp @@ -31,6 +31,7 @@ #include "RigNoRefinement.h" #include "RigNonUniformRefinement.h" #include "RigSimulationInputSettings.h" +#include "RigSimulationInputTestHelpers.h" #include "RigSimulationInputTool.h" #include "RigTransmissibilityEquations.h" #include "RigUniformRefinement.h" @@ -52,31 +53,8 @@ #include #include -static RigUniformRefinement makeRef( const caf::VecIjk0& sectorMin, const caf::VecIjk0& sectorMax, const cvf::Vec3st& refCounts ) -{ - cvf::Vec3st sectorSize( sectorMax.i() - sectorMin.i() + 1, sectorMax.j() - sectorMin.j() + 1, sectorMax.k() - sectorMin.k() + 1 ); - return RigUniformRefinement( refCounts, sectorSize ); -} - -//-------------------------------------------------------------------------------------------------- -/// Helper: set refinement on settings, computing sector size from min/max. -/// If refCounts is (1,1,1), creates RigNoRefinement; otherwise RigUniformRefinement. -//-------------------------------------------------------------------------------------------------- -static void setSettingsRefinement( RigSimulationInputSettings& settings, const cvf::Vec3st& refCounts ) -{ - cvf::Vec3st sectorSize( settings.max().i() - settings.min().i() + 1, - settings.max().j() - settings.min().j() + 1, - settings.max().k() - settings.min().k() + 1 ); - - if ( refCounts.x() == 1 && refCounts.y() == 1 && refCounts.z() == 1 ) - { - settings.setRefinement( std::make_unique( sectorSize ) ); - } - else - { - settings.setRefinement( std::make_unique( refCounts, sectorSize ) ); - } -} +using RigSimulationInputTestHelpers::makeRef; +using RigSimulationInputTestHelpers::setSettingsRefinement; //-------------------------------------------------------------------------------------------------- /// Helper to create an EQUALS DeckRecord @@ -281,126 +259,6 @@ TEST( RigSimulationInputTool, ProcessEqualsRecord_AtBoundary ) EXPECT_EQ( 10, result->getItem( 7 ).get( 0 ) ); // K2 } -//-------------------------------------------------------------------------------------------------- -/// Helper to create an AQUCON DeckRecord -/// AQUCON format: ID I1 I2 J1 J2 K1 K2 CONNECT_FACE ... (1-based Eclipse coordinates) -//-------------------------------------------------------------------------------------------------- -static Opm::DeckRecord createAquconRecord( int aquiferId, int i1, int i2, int j1, int j2, int k1, int k2, const std::string& connectFace ) -{ - std::vector items; - items.push_back( RifOpmDeckTools::item( "ID", aquiferId ) ); - items.push_back( RifOpmDeckTools::item( "I1", i1 ) ); - items.push_back( RifOpmDeckTools::item( "I2", i2 ) ); - items.push_back( RifOpmDeckTools::item( "J1", j1 ) ); - items.push_back( RifOpmDeckTools::item( "J2", j2 ) ); - items.push_back( RifOpmDeckTools::item( "K1", k1 ) ); - items.push_back( RifOpmDeckTools::item( "K2", k2 ) ); - items.push_back( RifOpmDeckTools::item( "CONNECT_FACE", connectFace ) ); - - return Opm::DeckRecord{ std::move( items ) }; -} - -//-------------------------------------------------------------------------------------------------- -/// Test AQUCON record with partial overlap (based on user-provided data) -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessAquconRecord_PartialOverlap ) -{ - // Sector: min=(40, 27, 21), max=(100, 72, 58) (inclusive, 0-based) - // This corresponds to a sector with cells [40-100, 27-72, 21-58] - caf::VecIjk0 sectorMin( 40, 27, 21 ); - caf::VecIjk0 sectorMax( 100, 72, 58 ); - auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); - - // AQUCON record from user data: 1 57 57 28 36 46 58 'I+' (1-based Eclipse) - // Converts to 0-based: I[56-56], J[27-35], K[45-57] - // Sector is I[40-100], J[27-72], K[21-58] - // Overlap: I[56-56], J[27-35], K[45-57] - all inside, no clamping needed - auto record = createAquconRecord( 1, 57, 57, 28, 36, 46, 58, "I+" ); - - auto result = RigSimulationInputTool::processAquconRecord( record, sectorMin, sectorMax, refinement ); - - ASSERT_TRUE( result.has_value() ); - - // Check transformed coordinates (sector-relative, 1-based) - // I: (56-40)*1+1 = 17 to (56-40)*1+1 = 17 - // J: (27-27)*1+1 = 1 to (35-27)*1+1 = 9 - // K: (45-21)*1+1 = 25 to (57-21)*1+1 = 37 - EXPECT_EQ( 17, result->getItem( 1 ).get( 0 ) ); // I1 - EXPECT_EQ( 17, result->getItem( 2 ).get( 0 ) ); // I2 - EXPECT_EQ( 1, result->getItem( 3 ).get( 0 ) ); // J1 - EXPECT_EQ( 9, result->getItem( 4 ).get( 0 ) ); // J2 - EXPECT_EQ( 25, result->getItem( 5 ).get( 0 ) ); // K1 - EXPECT_EQ( 37, result->getItem( 6 ).get( 0 ) ); // K2 - - // Check aquifer ID and connect face are preserved - EXPECT_EQ( 1, result->getItem( 0 ).get( 0 ) ); - EXPECT_EQ( "I+", result->getItem( 7 ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// Test AQUCON record with partial K overlap requiring clamping -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessAquconRecord_PartialOverlapWithClamping ) -{ - // Sector: min=(40, 27, 21), max=(100, 72, 58) (inclusive, 0-based) - caf::VecIjk0 sectorMin( 40, 27, 21 ); - caf::VecIjk0 sectorMax( 100, 72, 58 ); - auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); - - // AQUCON record: 1 79 79 41 67 5 11 'I+' (1-based Eclipse) - // Converts to 0-based: I[78-78], J[40-66], K[4-10] - // Sector K is [21-58], so K[4-10] is outside and should be skipped - auto record = createAquconRecord( 1, 79, 79, 41, 67, 5, 11, "I+" ); - - auto result = RigSimulationInputTool::processAquconRecord( record, sectorMin, sectorMax, refinement ); - - // This should fail because K range [4-10] doesn't overlap with sector K range [21-58] - EXPECT_FALSE( result.has_value() ); - EXPECT_TRUE( result.error().contains( "does not overlap" ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// Test AQUCON record completely inside sector -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessAquconRecord_CompletelyInside ) -{ - // Sector: min=(40, 27, 21), max=(100, 72, 58) (inclusive, 0-based) - caf::VecIjk0 sectorMin( 40, 27, 21 ); - caf::VecIjk0 sectorMax( 100, 72, 58 ); - auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); - - // AQUCON record: 1 61 61 48 72 12 17 'I+' (1-based Eclipse) - // Converts to 0-based: I[60-60], J[47-71], K[11-16] - // However, K[11-16] is outside sector K[21-58], so this should fail - auto record = createAquconRecord( 1, 61, 61, 48, 72, 12, 17, "I+" ); - - auto result = RigSimulationInputTool::processAquconRecord( record, sectorMin, sectorMax, refinement ); - - EXPECT_FALSE( result.has_value() ); -} - -//-------------------------------------------------------------------------------------------------- -/// Test AQUCON record without box definition (only aquifer ID) -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessAquconRecord_NoBoxDefinition ) -{ - caf::VecIjk0 sectorMin( 0, 0, 0 ); - caf::VecIjk0 sectorMax( 9, 9, 9 ); - auto refinement = makeRef( sectorMin, sectorMax, cvf::Vec3st( 1, 1, 1 ) ); - - // Create a record with only aquifer ID (no box coordinates) - std::vector items; - items.push_back( RifOpmDeckTools::item( "ID", 1 ) ); - - Opm::DeckRecord record{ std::move( items ) }; - - auto result = RigSimulationInputTool::processAquconRecord( record, sectorMin, sectorMax, refinement ); - - EXPECT_TRUE( result.has_value() ); - EXPECT_EQ( result->size(), 1U ); - EXPECT_EQ( 1, result->getItem( 0 ).get( 0 ) ); -} - //-------------------------------------------------------------------------------------------------- /// Test exportSimulationInput with model5 data //-------------------------------------------------------------------------------------------------- @@ -1517,162 +1375,6 @@ TEST( RigSimulationInputTool, ExportModel5WithRefinement_3_5_1 ) EXPECT_TRUE( checkCompdatCoordinates( "B-2H", 1, 6 ) ) << "B-2H COMPDAT K coordinates not centered correctly"; } -//-------------------------------------------------------------------------------------------------- -/// Helper to load model5 test case -//-------------------------------------------------------------------------------------------------- -static cvf::ref loadModel5TestCase() -{ - QDir baseFolder( TEST_DATA_DIR ); - bool subFolderExists = baseFolder.cd( "RigSimulationInputTool/model5" ); - if ( !subFolderExists ) return nullptr; - - QString egridFilename( "0_BASE_MODEL5.EGRID" ); - QString egridFilePath = baseFolder.absoluteFilePath( egridFilename ); - if ( !QFile::exists( egridFilePath ) ) return nullptr; - - std::unique_ptr resultCase( new RimEclipseResultCase ); - cvf::ref caseData = new RigEclipseCaseData( resultCase.get() ); - - cvf::ref reader = new RifReaderEclipseOutput; - bool loadResult = reader->open( egridFilePath, caseData.p() ); - if ( !loadResult ) return nullptr; - - return caseData; -} - -//-------------------------------------------------------------------------------------------------- -/// Test filterInternalSectorConnections - both cells inside -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, FilterNNCConnections_BothInside ) -{ - // Load model5 test data (20x30x10 grid) - auto caseData = loadModel5TestCase(); - ASSERT_TRUE( caseData.notNull() ); - ASSERT_TRUE( caseData->mainGrid() != nullptr ); - EXPECT_EQ( 20u, caseData->mainGrid()->cellCountI() ); - EXPECT_EQ( 30u, caseData->mainGrid()->cellCountJ() ); - EXPECT_EQ( 10u, caseData->mainGrid()->cellCountK() ); - - RigMainGrid* mainGrid = caseData->mainGrid(); - - // Create test connections using the real grid - std::vector allConnections; - - // Connection 1: (2,3,4) to (3,3,4) - both inside sector [1,5] x [1,5] x [1,5] - size_t c1Idx = mainGrid->cellIndexFromIJK( 2, 3, 4 ); - size_t c2Idx = mainGrid->cellIndexFromIJK( 3, 3, 4 ); - allConnections.push_back( { c1Idx, c2Idx, 1.5 } ); - - // Connection 2: (10,20,8) to (11,20,8) - both outside sector - c1Idx = mainGrid->cellIndexFromIJK( 10, 20, 8 ); - c2Idx = mainGrid->cellIndexFromIJK( 11, 20, 8 ); - allConnections.push_back( { c1Idx, c2Idx, 2.5 } ); - - // Connection 3: (1,2,3) to (15,15,8) - one inside, one outside - c1Idx = mainGrid->cellIndexFromIJK( 1, 2, 3 ); - c2Idx = mainGrid->cellIndexFromIJK( 15, 15, 8 ); - allConnections.push_back( { c1Idx, c2Idx, 3.5 } ); - - // Filter for sector [1,5] x [1,5] x [1,5] - caf::VecIjk0 sectorMin( 1, 1, 1 ); - caf::VecIjk0 sectorMax( 5, 5, 5 ); - - auto filtered = RigSimulationInputTool::filterInternalSectorConnections( allConnections, *mainGrid, sectorMin, sectorMax ); - - // Only connection 1 should remain - ASSERT_EQ( 1u, filtered.size() ); - EXPECT_EQ( 1.5, filtered[0].transmissibility ); -} - -//-------------------------------------------------------------------------------------------------- -/// Test transformNNCToSectorCoordinates - no refinement -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, TransformNNCToSectorCoordinates_NoRefinement ) -{ - // Load model5 test data (20x30x10 grid) - auto caseData = loadModel5TestCase(); - ASSERT_TRUE( caseData.notNull() ); - RigMainGrid* mainGrid = caseData->mainGrid(); - - // Connection between (10,15,5) and (11,15,5) - I-face neighbors - size_t c1Idx = mainGrid->cellIndexFromIJK( 10, 15, 5 ); - size_t c2Idx = mainGrid->cellIndexFromIJK( 11, 15, 5 ); - - RigSimulationInputTool::NNCConnection connection{ c1Idx, c2Idx, 4.5 }; - - // Sector starts at (5,10,2), grid is 20x30x10 - caf::VecIjk0 sectorMin( 5, 10, 2 ); - auto refinement = RigUniformRefinement( cvf::Vec3st( 1, 1, 1 ), cvf::Vec3st( 20, 30, 10 ) ); - - auto result = RigSimulationInputTool::transformNNCToSectorCoordinates( connection, *mainGrid, sectorMin, refinement ); - - ASSERT_TRUE( result.has_value() ); - - // Cell 1: (10,15,5) - sector min (5,10,2) = (5,5,3) in 0-based sector coords - EXPECT_EQ( 5, result->cell1.i() ); - EXPECT_EQ( 5, result->cell1.j() ); - EXPECT_EQ( 3, result->cell1.k() ); - - // Cell 2: (11,15,5) - sector min (5,10,2) = (6,5,3) in 0-based sector coords - EXPECT_EQ( 6, result->cell2.i() ); - EXPECT_EQ( 5, result->cell2.j() ); - EXPECT_EQ( 3, result->cell2.k() ); - - EXPECT_EQ( 4.5, result->transmissibility ); -} - -//-------------------------------------------------------------------------------------------------- -/// Test refineEDITNNCConnection - corresponding subcells with refinement (2,2,1) -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, RefineEDITNNCConnection_CorrespondingSubcells ) -{ - // Load model5 test data (20x30x10 grid) - auto caseData = loadModel5TestCase(); - ASSERT_TRUE( caseData.notNull() ); - RigMainGrid* mainGrid = caseData->mainGrid(); - - // Connection between cells (3,4,5) and (6,8,7) - non-neighbors - size_t c1Idx = mainGrid->cellIndexFromIJK( 3, 4, 5 ); - size_t c2Idx = mainGrid->cellIndexFromIJK( 6, 8, 7 ); - - RigSimulationInputTool::NNCConnection connection{ c1Idx, c2Idx, 2.5 }; // TRAN_MULT = 2.5 - - caf::VecIjk0 sectorMin( 0, 0, 0 ); - auto refinement = RigUniformRefinement( cvf::Vec3st( 2, 2, 1 ), cvf::Vec3st( 20, 30, 10 ) ); // Refine by 2x2x1, grid is 20x30x10 - - auto refined = RigSimulationInputTool::refineEditNncConnection( connection, *mainGrid, sectorMin, refinement ); - - // For EDITNNC: creates connections between corresponding subcells - // Should create 2*2*1 = 4 connections (one per subcell pair at same relative position) - ASSERT_EQ( 4u, refined.size() ); - - // All connections should have same TRAN_MULT (no distribution) - for ( const auto& conn : refined ) - { - EXPECT_DOUBLE_EQ( 2.5, conn.transmissibility ); - } - - // Verify corresponding subcell connections - // Cell1 starts at sector coordinates (3*2, 4*2, 5*1) = (6, 8, 5) - // Cell2 starts at sector coordinates (6*2, 8*2, 7*1) = (12, 16, 7) - - // Connection 0: (0,0,0) -> (0,0,0): (6,8,5) -> (12,16,7) - EXPECT_EQ( caf::VecIjk0( 6, 8, 5 ), refined[0].cell1 ); - EXPECT_EQ( caf::VecIjk0( 12, 16, 7 ), refined[0].cell2 ); - - // Connection 1: (0,1,0) -> (0,1,0): (6,9,5) -> (12,17,7) - EXPECT_EQ( caf::VecIjk0( 6, 9, 5 ), refined[1].cell1 ); - EXPECT_EQ( caf::VecIjk0( 12, 17, 7 ), refined[1].cell2 ); - - // Connection 2: (1,0,0) -> (1,0,0): (7,8,5) -> (13,16,7) - EXPECT_EQ( caf::VecIjk0( 7, 8, 5 ), refined[2].cell1 ); - EXPECT_EQ( caf::VecIjk0( 13, 16, 7 ), refined[2].cell2 ); - - // Connection 3: (1,1,0) -> (1,1,0): (7,9,5) -> (13,17,7) - EXPECT_EQ( caf::VecIjk0( 7, 9, 5 ), refined[3].cell1 ); - EXPECT_EQ( caf::VecIjk0( 13, 17, 7 ), refined[3].cell2 ); -} - //-------------------------------------------------------------------------------------------------- /// Test exportSimulationInput with model5 data and model padding //-------------------------------------------------------------------------------------------------- @@ -3060,371 +2762,3 @@ TEST( RigSimulationInputTool, ExportModel5WithBcconBcpropAndPadding ) EXPECT_EQ( 1, minKForNegZ ) << "\"Z-\" exterior face must remain at K=1 (top of padded grid)"; EXPECT_EQ( 1, maxKForNegZ ) << "\"Z-\" exterior face must remain at K=1 (top of padded grid)"; } - -//-------------------------------------------------------------------------------------------------- -/// Helper: parse a keyword from deck text and return all its records -//-------------------------------------------------------------------------------------------------- -static std::vector parseKeywordRecords( const std::string& deckText, const std::string& keywordName ) -{ - Opm::Parser parser; - auto deck = parser.parseString( deckText ); - const auto& keyword = deck[keywordName].back(); - - std::vector records; - for ( size_t i = 0; i < keyword.size(); ++i ) - { - records.push_back( keyword.getRecord( i ) ); - } - return records; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -static RigSimulationInputSettings createSectorSettings( const cvf::Vec3st& refCounts ) -{ - RigSimulationInputSettings settings; - settings.setMin( caf::VecIjk0( 0, 0, 0 ) ); - settings.setMax( caf::VecIjk0( 9, 9, 9 ) ); - setSettingsRefinement( settings, refCounts ); - return settings; -} - -//-------------------------------------------------------------------------------------------------- -/// A Z-direction connection in a K-refined layer must be completed in every sub-cell, with explicit -/// CF and KH split between the sub-cells (regression: only every other refined layer was completed) -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_ZDirectionSplitAcrossKSubcells ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); - auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 100.0 2* 'Z' /\n/\n", "COMPDAT" ); - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - - for ( size_t s = 0; s < 2; ++s ) - { - const auto& rec = result->at( s ); - EXPECT_EQ( 3, rec.getItem( "I" ).get( 0 ) ); - EXPECT_EQ( 4, rec.getItem( "J" ).get( 0 ) ); - EXPECT_EQ( 3 + static_cast( s ), rec.getItem( "K1" ).get( 0 ) ); - EXPECT_EQ( 3 + static_cast( s ), rec.getItem( "K2" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 5.0, rec.getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 50.0, rec.getItem( "Kh" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 0.216, rec.getItem( "DIAMETER" ).get( 0 ) ); - } -} - -//-------------------------------------------------------------------------------------------------- -/// With defaulted CF and KH the simulator computes the factors, so the full refined K range is kept in one record -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_DefaultedFactorsKeepsSingleRecord ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); - auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 3 'OPEN' 2* 0.216 /\n/\n", "COMPDAT" ); - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 1u, result->size() ); - EXPECT_EQ( 3, result->at( 0 ).getItem( "K1" ).get( 0 ) ); - EXPECT_EQ( 6, result->at( 0 ).getItem( "K2" ).get( 0 ) ); - EXPECT_FALSE( result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).hasValue( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// An X-direction connection is split across the I sub-cells; J and K use the center sub-cell -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_XDirectionSplitAcrossISubcells ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 3, 2, 2 ) ); - auto records = parseKeywordRecords( "COMPDAT\n 'W' 2 1 1 1 'OPEN' 1* 9.0 0.216 90.0 2* 'X' /\n/\n", "COMPDAT" ); - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 3u, result->size() ); - - for ( size_t s = 0; s < 3; ++s ) - { - const auto& rec = result->at( s ); - EXPECT_EQ( 4 + static_cast( s ), rec.getItem( "I" ).get( 0 ) ); - EXPECT_EQ( 2, rec.getItem( "J" ).get( 0 ) ); - EXPECT_EQ( 2, rec.getItem( "K1" ).get( 0 ) ); - EXPECT_EQ( 2, rec.getItem( "K2" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 3.0, rec.getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 30.0, rec.getItem( "Kh" ).get( 0 ) ); - } -} - -//-------------------------------------------------------------------------------------------------- -/// Non-uniform K refinement splits the CF by the sub-cell thickness fractions -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_NonUniformKRefinementScalesByFraction ) -{ - RigSimulationInputSettings settings; - settings.setMin( caf::VecIjk0( 0, 0, 0 ) ); - settings.setMax( caf::VecIjk0( 9, 9, 9 ) ); - auto refinement = std::make_unique( cvf::Vec3st( 10, 10, 10 ) ); - refinement->setCumulativeFractions( RigRefinement::DimK, 1, { 0.25, 1.0 } ); - settings.setRefinement( std::move( refinement ) ); - - auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 8.0 0.216 /\n/\n", "COMPDAT" ); - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - EXPECT_EQ( 2, result->at( 0 ).getItem( "K1" ).get( 0 ) ); - EXPECT_EQ( 3, result->at( 1 ).getItem( "K1" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 2.0, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 6.0, result->at( 1 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); - EXPECT_TRUE( result->at( 0 ).getItem( "Kh" ).defaultApplied( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// Without refinement the record is passed through unchanged -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_NoRefinementUnchanged ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 1 ) ); - auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 5 'OPEN' 1* 10.0 0.216 100.0 2* 'Z' /\n/\n", "COMPDAT" ); - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 1u, result->size() ); - EXPECT_EQ( 2, result->at( 0 ).getItem( "K1" ).get( 0 ) ); - EXPECT_EQ( 5, result->at( 0 ).getItem( "K2" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 10.0, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// COMPSEGS cells are split along the connection direction with the MD interval divided between sub-cells -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompsegsRecord_SplitAcrossKSubcells ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); - - Opm::Parser parser; - auto deck = parser.parseString( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 /\n/\n" - "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 120.0 /\n/\n" ); - - std::vector compdatKeywords = { deck["COMPDAT"].back() }; - auto directions = RigSimulationInputTool::connectionDirections( compdatKeywords ); - - const auto& compsegs = deck["COMPSEGS"].back(); - auto result = RigSimulationInputTool::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, directions ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - - EXPECT_EQ( 3, result->at( 0 ).getItem( 2 ).get( 0 ) ); - EXPECT_EQ( 4, result->at( 1 ).getItem( 2 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 100.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 110.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 110.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 120.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// Refinement perpendicular to a Z well reduces the Peaceman radius given by PR, so CF increases -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_LateralRefinementCorrectsCFUsingPR ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 10, 10, 1 ) ); - auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 100.0 1.0 1* 'Z' 20.0 /\n/\n", "COMPDAT" ); - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 1u, result->size() ); - - const auto& rec = result->at( 0 ); - EXPECT_EQ( 26, rec.getItem( "I" ).get( 0 ) ); - EXPECT_EQ( 36, rec.getItem( "J" ).get( 0 ) ); - - const double rw = 0.108; - const double expected = 10.0 * ( std::log( 20.0 / rw ) + 1.0 ) / ( std::log( 2.0 / rw ) + 1.0 ); - EXPECT_NEAR( expected, rec.getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ), 1e-9 ); - EXPECT_NEAR( 2.0, rec.getItem( "PR" ).get( 0 ), 1e-9 ); - EXPECT_DOUBLE_EQ( 100.0, rec.getItem( "Kh" ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// Without PR the Peaceman radii are computed from the cell sizes and permeabilities, also for -/// refinement in only one of the perpendicular directions -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_LateralRefinementCorrectsCFUsingCellProperties ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 4, 1, 1 ) ); - auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 /\n/\n", "COMPDAT" ); - - RigSimulationInputTool::ConnectionCellProperties properties{ cvf::Vec3d( 100.0, 50.0, 5.0 ), cvf::Vec3d( 200.0, 50.0, 10.0 ) }; - auto provider = [&]( const caf::VecIjk0& ) -> std::optional { return properties; }; - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings, provider ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 1u, result->size() ); - - const double rw = 0.108; - const double r0Coarse = RigTransmissibilityEquations::peacemanRadius( 200.0, 50.0, 100.0, 50.0 ); - const double r0Fine = RigTransmissibilityEquations::peacemanRadius( 200.0, 50.0, 25.0, 50.0 ); - const double expected = 10.0 * std::log( r0Coarse / rw ) / std::log( r0Fine / rw ); - EXPECT_NEAR( expected, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ), 1e-9 ); -} - -//-------------------------------------------------------------------------------------------------- -/// The correction is skipped when the refined Peaceman radius gets smaller than the well radius -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompdatRecord_LateralRefinementTooSmallRadiusKeepsCF ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 10, 10, 1 ) ); - auto records = parseKeywordRecords( "COMPDAT\n 'W' 3 4 2 2 'OPEN' 1* 10.0 0.216 4* 0.5 /\n/\n", "COMPDAT" ); - - auto result = RigSimulationInputTool::processCompdatRecord( records[0], "W", settings ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 1u, result->size() ); - EXPECT_DOUBLE_EQ( 10.0, result->at( 0 ).getItem( "CONNECTION_TRANSMISSIBILITY_FACTOR" ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 0.5, result->at( 0 ).getItem( "PR" ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// A well travelling towards decreasing I assigns the MD interval from the highest to the lowest I sub-cell -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompsegsRecord_ReversedTraversalAlongI ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 2, 1, 1 ) ); - - Opm::Parser parser; - auto deck = parser.parseString( "COMPDAT\n 'W' 4 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n" - " 'W' 3 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n/\n" - "COMPSEGS\n 'W' /\n 4 2 1 1 100.0 120.0 /\n 3 2 1 1 120.0 140.0 /\n/\n" ); - - std::vector compdatKeywords = { deck["COMPDAT"].back() }; - std::vector compsegsKeywords = { deck["COMPSEGS"].back() }; - - auto directions = RigSimulationInputTool::connectionDirections( compdatKeywords ); - auto reversedCells = RigSimulationInputTool::reversedSegmentCells( compsegsKeywords, directions ); - - EXPECT_TRUE( reversedCells.contains( { "W", caf::VecIjk0( 3, 1, 0 ) } ) ); - EXPECT_TRUE( reversedCells.contains( { "W", caf::VecIjk0( 2, 1, 0 ) } ) ); - - const auto& compsegs = deck["COMPSEGS"].back(); - - // Original I=4 -> refined I sub-cells 7 and 8, the well enters at 8 - auto result = RigSimulationInputTool::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, directions, reversedCells ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - EXPECT_EQ( 8, result->at( 0 ).getItem( 0 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 100.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 110.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); - EXPECT_EQ( 7, result->at( 1 ).getItem( 0 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 110.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 120.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); - - // Original I=3 is the last cell, direction taken from the previous cell - result = RigSimulationInputTool::processCompsegsRecord( compsegs.getRecord( 2 ), "W", false, settings, directions, reversedCells ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - EXPECT_EQ( 6, result->at( 0 ).getItem( 0 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 120.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); - EXPECT_EQ( 5, result->at( 1 ).getItem( 0 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 140.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// A well travelling towards increasing I keeps the sub-cell order -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ReversedSegmentCells_ForwardTraversalNotReversed ) -{ - Opm::Parser parser; - auto deck = parser.parseString( "COMPDAT\n 'W' 3 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n" - " 'W' 4 2 1 1 'OPEN' 1* 10.0 0.216 3* 'X' /\n/\n" - "COMPSEGS\n 'W' /\n 3 2 1 1 100.0 120.0 /\n 4 2 1 1 120.0 140.0 /\n/\n" ); - - std::vector compdatKeywords = { deck["COMPDAT"].back() }; - std::vector compsegsKeywords = { deck["COMPSEGS"].back() }; - - auto directions = RigSimulationInputTool::connectionDirections( compdatKeywords ); - auto reversedCells = RigSimulationInputTool::reversedSegmentCells( compsegsKeywords, directions ); - - EXPECT_TRUE( reversedCells.empty() ); -} - -//-------------------------------------------------------------------------------------------------- -/// DIRECTION without END_IJK is a single cell, so the MD interval is still split between the sub-cells -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompsegsRecord_DirectionWithoutRangeIsSplit ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); - - Opm::Parser parser; - auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 120.0 'K' /\n/\n" ); - - const auto& compsegs = deck["COMPSEGS"].back(); - auto result = RigSimulationInputTool::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - EXPECT_EQ( 3, result->at( 0 ).getItem( 2 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 110.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); - EXPECT_EQ( 4, result->at( 1 ).getItem( 2 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 110.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); - EXPECT_EQ( "K", result->at( 1 ).getItem( 6 ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// A defaulted start MD is zero -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompsegsRecord_DefaultedStartMdIsZero ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); - - Opm::Parser parser; - auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 1* 20.0 /\n/\n" ); - - const auto& compsegs = deck["COMPSEGS"].back(); - auto result = RigSimulationInputTool::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - EXPECT_DOUBLE_EQ( 0.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 10.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 10.0, result->at( 1 ).getItem( 4 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 20.0, result->at( 1 ).getItem( 5 ).get( 0 ) ); -} - -//-------------------------------------------------------------------------------------------------- -/// Without an end MD every sub-cell still gets a segment record, with the MD items unchanged -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompsegsRecord_NoEndMdCopiedToAllSubcells ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); - - Opm::Parser parser; - auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 1* 'K' /\n/\n" ); - - const auto& compsegs = deck["COMPSEGS"].back(); - auto result = RigSimulationInputTool::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 2u, result->size() ); - for ( size_t s = 0; s < 2; ++s ) - { - EXPECT_EQ( 3 + static_cast( s ), result->at( s ).getItem( 2 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 100.0, result->at( s ).getItem( 4 ).get( 0 ) ); - EXPECT_FALSE( result->at( s ).getItem( 5 ).hasValue( 0 ) ); - } -} - -//-------------------------------------------------------------------------------------------------- -/// A cell range (END_IJK) is not split, only the center sub-cell is used -//-------------------------------------------------------------------------------------------------- -TEST( RigSimulationInputTool, ProcessCompsegsRecord_CellRangeUsesCenterSubcell ) -{ - auto settings = createSectorSettings( cvf::Vec3st( 1, 1, 2 ) ); - - Opm::Parser parser; - auto deck = parser.parseString( "COMPSEGS\n 'W' /\n 3 4 2 1 100.0 120.0 'K' 3 /\n/\n" ); - - const auto& compsegs = deck["COMPSEGS"].back(); - auto result = RigSimulationInputTool::processCompsegsRecord( compsegs.getRecord( 1 ), "W", false, settings, {} ); - ASSERT_TRUE( result.has_value() ) << result.error().toStdString(); - ASSERT_EQ( 1u, result->size() ); - EXPECT_EQ( 4, result->at( 0 ).getItem( 2 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 100.0, result->at( 0 ).getItem( 4 ).get( 0 ) ); - EXPECT_DOUBLE_EQ( 120.0, result->at( 0 ).getItem( 5 ).get( 0 ) ); -}