mirror of
https://github.com/OPM/ResInsight.git
synced 2026-10-24 22:27:55 -05:00
#14510 Nested hybrid grid: prototype FIPNEST parent-child arrays through the INIT file
Adds a compact FIPNEST/FIPSLOT/REFINE encoding of the nested hybrid grid parent-child relations that can travel inside the INIT file instead of the OLDIJK sidecar: FIPNEST is the 1-based flat natural index of each cell's immediate parent (refined-away host slots carry their own link so every chain terminates at a coarse cell), FIPSLOT packs the within-parent i/j/k offsets, and REFINE is the existing per-cell level. FIPNEST alone is provably insufficient: the refinement level is not the parent-chain depth (L2 and L3 both refine the coarse grid directly in the DROGON dataset), and within-parent placement is ambiguous from flat coordinates when a level has hole layers. RigNestedHybridGridFipnestCodec derives the arrays from the sidecar input by mirroring the reconstructor's level classification, and synthesizes an equivalent reconstruction input back from them (identical up to per-level constant TMP shifts, which reconstruct() is invariant to), so reusing the reconstructor unchanged guarantees that a grid imported through FIPNEST equals the sidecar-based import. A successful sidecar reconstruction auto-exports <gridbase>_FIPNEST.grdecl next to the grid file, scripts/nested_hybrid_grid/ insert_fipnest_into_init.py inserts the arrays into an INIT file using resfo, and RimEclipseResultCase prefers the INIT-embedded arrays over the sidecars when both exist. Verified by unit tests covering the codec round-trip, array invariants, cellwise grid equality of the two import paths, the GRDECL writer, and a hermetic end-to-end pass that stages the model, appends the INTE arrays to the INIT copy and imports without sidecars.
This commit is contained in:
Kristian Bendiksen committed
1 parent
f742fdcac4
commit
d6722f4ba8
10 files changed
+1774
-4
No files matched your search
@@ -294,7 +294,8 @@ bool RimEclipseResultCase::importGridAndResultMetaData( bool showTimeStepFilter
|
||||
|
||||
const bool hasNestedHybridSidecars = !RigNestedHybridGridResultTools::refineSidecarFilePath( gridFileName() ).isEmpty() &&
|
||||
!RigNestedHybridGridResultTools::oldIjkSidecarFilePath( gridFileName() ).isEmpty();
|
||||
if ( hasNestedHybridSidecars )
|
||||
const bool hasFipnestInInit = RigNestedHybridGridResultTools::initFileHasFipnest( gridFileName() );
|
||||
if ( hasNestedHybridSidecars || hasFipnestInInit )
|
||||
{
|
||||
// The flat nested hybrid grid piles the refined and collapsed coarse cells into the same
|
||||
// physical space, which makes the geometric fault/NNC computation in computeCachedData()
|
||||
@@ -304,8 +305,19 @@ bool RimEclipseResultCase::importGridAndResultMetaData( bool showTimeStepFilter
|
||||
// extended to cover the new LGR cells.
|
||||
loadAndSynchronizeInputProperties( false );
|
||||
RigNestedHybridGridResultTools::importRefineSidecarIfPresent( gridFileName(), inputPropertyCollection(), eclipseCaseData() );
|
||||
RigNestedHybridGridResultTools::importOldIjkSidecarIfPresent( gridFileName(), inputPropertyCollection(), eclipseCaseData() );
|
||||
RigNestedHybridGridResultTools::reconstructNestedHybridGridIfPresent( gridFileName(), eclipseCaseData() );
|
||||
|
||||
// The FIPNEST/FIPSLOT/REFINE arrays embedded in the INIT file take precedence over the
|
||||
// OLDIJK sidecar files; the sidecars remain as the fallback (#14510).
|
||||
bool reconstructed = false;
|
||||
if ( hasFipnestInInit )
|
||||
{
|
||||
reconstructed = RigNestedHybridGridResultTools::reconstructNestedHybridGridFromInitFile( gridFileName(), eclipseCaseData() );
|
||||
}
|
||||
if ( !reconstructed && hasNestedHybridSidecars )
|
||||
{
|
||||
RigNestedHybridGridResultTools::importOldIjkSidecarIfPresent( gridFileName(), inputPropertyCollection(), eclipseCaseData() );
|
||||
RigNestedHybridGridResultTools::reconstructNestedHybridGridIfPresent( gridFileName(), eclipseCaseData() );
|
||||
}
|
||||
computeCachedData();
|
||||
}
|
||||
else
|
||||
|
||||
@@ -89,6 +89,7 @@ set(SOURCE_GROUP_SOURCE_FILES
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigNonUniformRefinement.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigNestedHybridGridReconstructor.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigNestedHybridGridResultTools.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigNestedHybridGridFipnestCodec.cpp
|
||||
)
|
||||
|
||||
list(APPEND CODE_SOURCE_FILES ${SOURCE_GROUP_SOURCE_FILES})
|
||||
@@ -0,0 +1,672 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2026- Equinor ASA
|
||||
//
|
||||
// ResInsight is free software: you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation, either version 3 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
// FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
||||
// for more details.
|
||||
//
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#include "RigNestedHybridGridFipnestCodec.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <climits>
|
||||
#include <cmath>
|
||||
#include <map>
|
||||
#include <numeric>
|
||||
#include <set>
|
||||
#include <vector>
|
||||
|
||||
namespace
|
||||
{
|
||||
size_t naturalIndex( size_t i, size_t j, size_t k, size_t nx, size_t ny )
|
||||
{
|
||||
return i + j * nx + k * nx * ny;
|
||||
}
|
||||
|
||||
std::array<int, 3> flatIjk( size_t f, size_t nx, size_t ny )
|
||||
{
|
||||
return { (int)( f % nx ), (int)( ( f / nx ) % ny ), (int)( f / ( nx * ny ) ) };
|
||||
}
|
||||
|
||||
// A primary (coarse-refining) level, recorded so the next level can be encoded/decoded against it.
|
||||
struct PrimaryLevelRecord
|
||||
{
|
||||
int level = 0;
|
||||
int t0[3] = { 0, 0, 0 }; // min TMP of the level's cells
|
||||
int c0[3] = { 0, 0, 0 }; // min OLD of the level's cells
|
||||
int factor[3] = { 1, 1, 1 }; // per-axis refinement factor
|
||||
int dims[3] = { 0, 0, 0 }; // TMP-box dimensions (coarseDim * factor)
|
||||
|
||||
// Per-axis linear map from the level's TMP space to the flat grid: flat = s*tmp + c. Lets a
|
||||
// refined-away (hole) host slot be located as a flat cell even though no real cell occupies it.
|
||||
bool bandValid = false;
|
||||
long long bandS[3] = { 0, 0, 0 };
|
||||
long long bandC[3] = { 0, 0, 0 };
|
||||
|
||||
std::map<std::array<int, 3>, size_t> realCellByTmp; // TMP triple -> flat cell of the level's real cells
|
||||
};
|
||||
} // namespace
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
int RigNestedHybridGridFipnestCodec::packSlot( int offI, int offJ, int offK )
|
||||
{
|
||||
return 1 + offI + 100 * offJ + 10000 * offK;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RigNestedHybridGridFipnestCodec::unpackSlot( int slot, int& offI, int& offJ, int& offK )
|
||||
{
|
||||
const int v = slot - 1;
|
||||
offI = v % 100;
|
||||
offJ = ( v / 100 ) % 100;
|
||||
offK = v / 10000;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Mirrors RigNestedHybridGridReconstructor::reconstruct()'s level classification: levels are
|
||||
/// processed shallow to deep; a level whose cells' TMPs fall inside the previous primary level's
|
||||
/// TMP box (>= 95%) is nested under that level, otherwise it is a direct refinement of the coarse
|
||||
/// grid. The immediate parent written to FIPNEST is the coarse host slot for primary-level cells
|
||||
/// and the previous level's (possibly refined-away) host slot for nested cells.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RigNestedHybridGridFipnestCodec::ParentChildArrays
|
||||
RigNestedHybridGridFipnestCodec::computeParentChildArrays( const RigNestedHybridGridReconstructor::NestedHybridInput& input,
|
||||
size_t nx,
|
||||
size_t ny,
|
||||
size_t nz )
|
||||
{
|
||||
ParentChildArrays result;
|
||||
|
||||
const size_t cellCount = nx * ny * nz;
|
||||
if ( input.refine.size() != cellCount || input.oldI.size() != cellCount || input.oldJ.size() != cellCount ||
|
||||
input.oldK.size() != cellCount || input.tmpI.size() != cellCount || input.tmpJ.size() != cellCount || input.tmpK.size() != cellCount )
|
||||
{
|
||||
return result; // empty arrays signal invalid input
|
||||
}
|
||||
|
||||
result.fipnest.assign( cellCount, 0 );
|
||||
result.fipslot.assign( cellCount, 0 );
|
||||
|
||||
// Same coarse-dimension and K-factor derivation as the reconstructor.
|
||||
int coarseNz = 0;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
coarseNz = std::max( coarseNz, input.oldK[f] );
|
||||
const size_t kFactor = ( coarseNz > 0 && nz % (size_t)coarseNz == 0 ) ? nz / (size_t)coarseNz : 1;
|
||||
|
||||
std::map<int, std::vector<size_t>> cellsByLevel;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
const int level = input.refine[f];
|
||||
if ( level <= 1 ) continue;
|
||||
if ( input.oldI[f] < 1 || input.oldJ[f] < 1 || input.oldK[f] < 1 )
|
||||
{
|
||||
result.unresolvedRefinedCells++; // padding / unmapped; skipped by the reconstructor too
|
||||
continue;
|
||||
}
|
||||
cellsByLevel[level].push_back( f );
|
||||
}
|
||||
|
||||
auto coarseHostFlatIndex = [&]( int oi, int oj, int ok )
|
||||
{ return naturalIndex( (size_t)( oi - 1 ), (size_t)( oj - 1 ), (size_t)( ok - 1 ) * kFactor, nx, ny ); };
|
||||
|
||||
std::vector<PrimaryLevelRecord> primaries;
|
||||
std::set<int> nestedLevels;
|
||||
|
||||
for ( const auto& [level, cells] : cellsByLevel )
|
||||
{
|
||||
if ( nestedLevels.count( level - 1 ) )
|
||||
{
|
||||
// The reconstructor would nest this level under the (nested) level below it, but the
|
||||
// encoding only supports nesting under a primary level (see class comment). Leaving the
|
||||
// cells unresolved is safer than mis-encoding them as a primary refinement.
|
||||
result.unresolvedRefinedCells += cells.size();
|
||||
continue;
|
||||
}
|
||||
|
||||
const PrimaryLevelRecord* prev = nullptr;
|
||||
for ( const PrimaryLevelRecord& p : primaries )
|
||||
if ( p.level == level - 1 ) prev = &p;
|
||||
|
||||
auto tmpOf = [&]( size_t f ) { return std::array<int, 3>{ input.tmpI[f], input.tmpJ[f], input.tmpK[f] }; };
|
||||
auto inPrevBox = [&]( const std::array<int, 3>& t )
|
||||
{
|
||||
if ( !prev ) return false;
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
if ( t[a] < prev->t0[a] || t[a] >= prev->t0[a] + prev->dims[a] ) return false;
|
||||
return true;
|
||||
};
|
||||
|
||||
size_t inBox = 0;
|
||||
for ( size_t f : cells )
|
||||
if ( inPrevBox( tmpOf( f ) ) ) inBox++;
|
||||
|
||||
if ( prev && inBox >= cells.size() * 95 / 100 && inBox > 0 )
|
||||
{
|
||||
// Nested level: each cell's immediate parent is the previous level's cell (or hole slot)
|
||||
// at the cell's own TMP - the same lookup the reconstructor performs via tmpToBuilt.
|
||||
std::map<std::array<int, 3>, std::array<int, 3>> groupFlatMin; // parent TMP -> min flat IJK of its children
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const auto t = tmpOf( f );
|
||||
if ( !inPrevBox( t ) ) continue;
|
||||
const auto fxyz = flatIjk( f, nx, ny );
|
||||
auto it = groupFlatMin.find( t );
|
||||
if ( it == groupFlatMin.end() )
|
||||
groupFlatMin[t] = fxyz;
|
||||
else
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
it->second[a] = std::min( it->second[a], fxyz[a] );
|
||||
}
|
||||
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const auto t = tmpOf( f );
|
||||
if ( !inPrevBox( t ) )
|
||||
{
|
||||
result.unresolvedRefinedCells++;
|
||||
continue;
|
||||
}
|
||||
|
||||
size_t parentFlat = 0;
|
||||
bool parentOk = false;
|
||||
if ( auto it = prev->realCellByTmp.find( t ); it != prev->realCellByTmp.end() )
|
||||
{
|
||||
parentFlat = it->second;
|
||||
parentOk = true;
|
||||
}
|
||||
else if ( prev->bandValid )
|
||||
{
|
||||
// The host was refined away; locate its collapsed flat slot through the band map
|
||||
// and give that slot its own chain link back to its coarse host.
|
||||
long long p[3];
|
||||
parentOk = true;
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
p[a] = prev->bandS[a] * t[a] + prev->bandC[a];
|
||||
const long long dim = ( a == 0 ) ? (long long)nx : ( a == 1 ) ? (long long)ny : (long long)nz;
|
||||
if ( p[a] < 0 || p[a] >= dim ) parentOk = false;
|
||||
}
|
||||
if ( parentOk )
|
||||
{
|
||||
parentFlat = naturalIndex( (size_t)p[0], (size_t)p[1], (size_t)p[2], nx, ny );
|
||||
|
||||
int hostOld[3], hostOff[3];
|
||||
bool hostOffOk = true;
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
hostOld[a] = prev->c0[a] + ( t[a] - prev->t0[a] ) / prev->factor[a];
|
||||
hostOff[a] = ( t[a] - prev->t0[a] ) % prev->factor[a];
|
||||
if ( hostOff[a] >= 100 ) hostOffOk = false; // beyond FIPSLOT digit capacity
|
||||
}
|
||||
const size_t hostFlat = coarseHostFlatIndex( hostOld[0], hostOld[1], hostOld[2] );
|
||||
|
||||
if ( input.refine[parentFlat] > 1 || !hostOffOk || hostFlat >= cellCount || input.refine[hostFlat] > 1 )
|
||||
{
|
||||
parentOk = false; // slot occupied by a refined cell, or host not encodable
|
||||
}
|
||||
else
|
||||
{
|
||||
const int expectedLink = (int)hostFlat + 1;
|
||||
const int expectedSlot = packSlot( hostOff[0], hostOff[1], hostOff[2] );
|
||||
if ( result.fipnest[parentFlat] == 0 )
|
||||
{
|
||||
result.fipnest[parentFlat] = expectedLink;
|
||||
result.fipslot[parentFlat] = expectedSlot;
|
||||
}
|
||||
else if ( result.fipnest[parentFlat] != expectedLink || result.fipslot[parentFlat] != expectedSlot )
|
||||
{
|
||||
parentOk = false; // the slot already carries a different chain link
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const auto fxyz = flatIjk( f, nx, ny );
|
||||
const auto& mn = groupFlatMin[t];
|
||||
const int off[3] = { fxyz[0] - mn[0], fxyz[1] - mn[1], fxyz[2] - mn[2] };
|
||||
if ( !parentOk || off[0] >= 100 || off[1] >= 100 || off[2] >= 100 )
|
||||
{
|
||||
result.unresolvedRefinedCells++;
|
||||
continue;
|
||||
}
|
||||
|
||||
result.fipnest[f] = (int)parentFlat + 1;
|
||||
result.fipslot[f] = packSlot( off[0], off[1], off[2] );
|
||||
}
|
||||
nestedLevels.insert( level );
|
||||
continue;
|
||||
}
|
||||
|
||||
// Primary level: a direct, uniform refinement of the coarse grid (same derivation and
|
||||
// uniformity gate as the reconstructor).
|
||||
int c0[3] = { INT_MAX, INT_MAX, INT_MAX }, c1[3] = { 0, 0, 0 };
|
||||
int t0[3] = { INT_MAX, INT_MAX, INT_MAX }, t1[3] = { 0, 0, 0 };
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const int t[3] = { input.tmpI[f], input.tmpJ[f], input.tmpK[f] };
|
||||
const int c[3] = { input.oldI[f], input.oldJ[f], input.oldK[f] };
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
t0[a] = std::min( t0[a], t[a] );
|
||||
t1[a] = std::max( t1[a], t[a] );
|
||||
c0[a] = std::min( c0[a], c[a] );
|
||||
c1[a] = std::max( c1[a], c[a] );
|
||||
}
|
||||
}
|
||||
|
||||
int factor[3];
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
const int coarseDim = c1[a] - c0[a] + 1;
|
||||
const int tmpDim = t1[a] - t0[a] + 1;
|
||||
factor[a] = std::max( 1, (int)std::lround( (double)tmpDim / (double)coarseDim ) );
|
||||
}
|
||||
|
||||
size_t matched = 0;
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const int pi = c0[0] + ( input.tmpI[f] - t0[0] ) / factor[0];
|
||||
const int pj = c0[1] + ( input.tmpJ[f] - t0[1] ) / factor[1];
|
||||
const int pk = c0[2] + ( input.tmpK[f] - t0[2] ) / factor[2];
|
||||
if ( pi == input.oldI[f] && pj == input.oldJ[f] && pk == input.oldK[f] ) matched++;
|
||||
}
|
||||
if ( matched < cells.size() * 95 / 100 )
|
||||
{
|
||||
result.unresolvedRefinedCells += cells.size(); // the reconstructor defers this level
|
||||
continue;
|
||||
}
|
||||
|
||||
PrimaryLevelRecord record;
|
||||
record.level = level;
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
record.t0[a] = t0[a];
|
||||
record.c0[a] = c0[a];
|
||||
record.factor[a] = factor[a];
|
||||
record.dims[a] = ( c1[a] - c0[a] + 1 ) * factor[a];
|
||||
}
|
||||
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const int t[3] = { input.tmpI[f], input.tmpJ[f], input.tmpK[f] };
|
||||
const int old[3] = { input.oldI[f], input.oldJ[f], input.oldK[f] };
|
||||
int off[3];
|
||||
bool offOk = true;
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
off[a] = ( t[a] - t0[a] ) - ( old[a] - c0[a] ) * factor[a];
|
||||
if ( off[a] < 0 || off[a] >= factor[a] || off[a] >= 100 ) offOk = false;
|
||||
}
|
||||
// The chain must terminate at the coarse host, so the flat cell occupying that slot has
|
||||
// to be unrefined (a refined cell there would carry its own, unrelated link).
|
||||
const size_t hostFlat = coarseHostFlatIndex( old[0], old[1], old[2] );
|
||||
if ( !offOk || hostFlat >= cellCount || input.refine[hostFlat] > 1 )
|
||||
{
|
||||
result.unresolvedRefinedCells++;
|
||||
continue;
|
||||
}
|
||||
|
||||
result.fipnest[f] = (int)hostFlat + 1;
|
||||
result.fipslot[f] = packSlot( off[0], off[1], off[2] );
|
||||
record.realCellByTmp[{ t[0], t[1], t[2] }] = f;
|
||||
}
|
||||
|
||||
// Fit the per-axis linear TMP -> flat band map from the level's real cells.
|
||||
record.bandValid = true;
|
||||
for ( int a = 0; a < 3 && record.bandValid; a++ )
|
||||
{
|
||||
long long flatAtT0 = -1, flatAtT1 = -1;
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const int t[3] = { input.tmpI[f], input.tmpJ[f], input.tmpK[f] };
|
||||
const int fa = flatIjk( f, nx, ny )[a];
|
||||
if ( t[a] == t0[a] ) flatAtT0 = fa;
|
||||
if ( t[a] == t1[a] ) flatAtT1 = fa;
|
||||
}
|
||||
if ( flatAtT0 < 0 || flatAtT1 < 0 )
|
||||
{
|
||||
record.bandValid = false;
|
||||
break;
|
||||
}
|
||||
if ( t1[a] == t0[a] )
|
||||
{
|
||||
record.bandS[a] = 0;
|
||||
record.bandC[a] = flatAtT0;
|
||||
}
|
||||
else
|
||||
{
|
||||
const long long dt = t1[a] - t0[a];
|
||||
const long long df = flatAtT1 - flatAtT0;
|
||||
if ( df % dt != 0 )
|
||||
{
|
||||
record.bandValid = false;
|
||||
break;
|
||||
}
|
||||
record.bandS[a] = df / dt;
|
||||
record.bandC[a] = flatAtT0 - record.bandS[a] * t0[a];
|
||||
}
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const int t[3] = { input.tmpI[f], input.tmpJ[f], input.tmpK[f] };
|
||||
if ( record.bandS[a] * t[a] + record.bandC[a] != flatIjk( f, nx, ny )[a] )
|
||||
{
|
||||
record.bandValid = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
primaries.push_back( record );
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Synthesizes OLDI/J/K from each cell's parent chain and TMPI/J/K from the recovered per-level
|
||||
/// refinement factors and the FIPSLOT offsets. The synthesized TMPs equal the original sidecar TMPs
|
||||
/// up to a constant per-level per-axis shift, which reconstruct() is invariant to; disjoint shifts
|
||||
/// are chosen so no TMP triple of one level collides with another level's slot registration.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RigNestedHybridGridReconstructor::NestedHybridInput
|
||||
RigNestedHybridGridFipnestCodec::buildInputFromParentChildArrays( const std::vector<int>& fipnest,
|
||||
const std::vector<int>& fipslot,
|
||||
const std::vector<int>& refine,
|
||||
size_t nx,
|
||||
size_t ny,
|
||||
size_t nz,
|
||||
QString* warnings )
|
||||
{
|
||||
RigNestedHybridGridReconstructor::NestedHybridInput input;
|
||||
|
||||
auto warn = [&]( const QString& msg )
|
||||
{
|
||||
if ( warnings ) *warnings += ( warnings->isEmpty() ? "" : "\n" ) + msg;
|
||||
};
|
||||
|
||||
const size_t cellCount = nx * ny * nz;
|
||||
if ( fipnest.size() != cellCount || fipslot.size() != cellCount || refine.size() != cellCount )
|
||||
{
|
||||
warn( "FIPNEST/FIPSLOT/REFINE size does not match the grid cell count." );
|
||||
return input;
|
||||
}
|
||||
|
||||
input.refine = refine;
|
||||
input.oldI.assign( cellCount, 0 );
|
||||
input.oldJ.assign( cellCount, 0 );
|
||||
input.oldK.assign( cellCount, 0 );
|
||||
input.tmpI.assign( cellCount, 0 );
|
||||
input.tmpJ.assign( cellCount, 0 );
|
||||
input.tmpK.assign( cellCount, 0 );
|
||||
|
||||
// Walk each chain to its FIPNEST == 0 root (the coarse host). Memoized; cycles yield no root.
|
||||
std::vector<long long> rootOfCell( cellCount, -2 ); // -2 unvisited, -1 invalid, else root flat index
|
||||
auto rootOf = [&]( size_t f ) -> long long
|
||||
{
|
||||
std::vector<size_t> path;
|
||||
size_t cur = f;
|
||||
while ( rootOfCell[cur] == -2 )
|
||||
{
|
||||
const int link = fipnest[cur];
|
||||
if ( link == 0 )
|
||||
{
|
||||
rootOfCell[cur] = (long long)cur;
|
||||
break;
|
||||
}
|
||||
if ( link < 1 || (size_t)link > cellCount || path.size() > 16 )
|
||||
{
|
||||
rootOfCell[cur] = -1; // out of range or unreasonably deep (cycle)
|
||||
break;
|
||||
}
|
||||
path.push_back( cur );
|
||||
cur = (size_t)link - 1;
|
||||
}
|
||||
const long long root = rootOfCell[cur];
|
||||
for ( size_t c : path )
|
||||
rootOfCell[c] = root;
|
||||
return rootOfCell[f];
|
||||
};
|
||||
|
||||
// K refinement factor: coarse hosts occupy every kFactor'th flat K layer.
|
||||
long long kFactor = 0;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
if ( refine[f] <= 1 || fipnest[f] == 0 ) continue;
|
||||
const long long root = rootOf( f );
|
||||
if ( root < 0 ) continue;
|
||||
kFactor = std::gcd( kFactor, (long long)( (size_t)root / ( nx * ny ) ) );
|
||||
}
|
||||
// gcd with NZ guarantees a divisor of NZ; with no refined K layers above 0 this yields NZ itself,
|
||||
// which maps every host to coarse K = 1 - consistent with how the reconstructor re-derives it.
|
||||
kFactor = std::gcd( kFactor, (long long)nz );
|
||||
|
||||
std::map<int, std::vector<size_t>> cellsByLevel;
|
||||
size_t unresolved = 0;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
if ( refine[f] <= 1 ) continue;
|
||||
const long long root = fipnest[f] != 0 ? rootOf( f ) : -1;
|
||||
if ( root < 0 || refine[(size_t)root] > 1 )
|
||||
{
|
||||
unresolved++; // no parent chain; the cell is left out (stays a flat cell)
|
||||
continue;
|
||||
}
|
||||
const auto r = flatIjk( (size_t)root, nx, ny );
|
||||
input.oldI[f] = r[0] + 1;
|
||||
input.oldJ[f] = r[1] + 1;
|
||||
input.oldK[f] = r[2] / (int)kFactor + 1;
|
||||
cellsByLevel[refine[f]].push_back( f );
|
||||
}
|
||||
if ( unresolved > 0 )
|
||||
{
|
||||
warn( QString( "%1 refined cells have no usable FIPNEST chain and are left un-nested." ).arg( unresolved ) );
|
||||
}
|
||||
|
||||
// Give every unrefined cell its own coarse position as OLD. The reconstructor re-derives the
|
||||
// K refinement factor from the maximum OLDK, so covering all NZ layers here pins that factor to
|
||||
// the one used for the chain roots above (the refined hosts alone may stop below the top layer,
|
||||
// which would silently shift every parent's K slot).
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
if ( refine[f] > 1 ) continue;
|
||||
const auto r = flatIjk( f, nx, ny );
|
||||
input.oldI[f] = r[0] + 1;
|
||||
input.oldJ[f] = r[1] + 1;
|
||||
input.oldK[f] = r[2] / (int)kFactor + 1;
|
||||
}
|
||||
|
||||
// Per-axis refinement factor and TMP shift of each synthesized primary level.
|
||||
struct LevelParams
|
||||
{
|
||||
int factor[3] = { 1, 1, 1 };
|
||||
int shift[3] = { 0, 0, 0 };
|
||||
};
|
||||
std::map<int, LevelParams> primaryParams;
|
||||
|
||||
int nextShift[3] = { 0, 0, 0 };
|
||||
|
||||
// Offsets above the FIPSLOT digit capacity cannot come from the encoder; treat them as corrupt
|
||||
// input rather than letting huge decoded offsets poison the factor recovery (or overflow).
|
||||
const int maxValidSlot = packSlot( 99, 99, 99 );
|
||||
auto slotValid = [&]( size_t f ) { return fipslot[f] > 0 && fipslot[f] <= maxValidSlot; };
|
||||
|
||||
auto offsetsOf = [&]( size_t f )
|
||||
{
|
||||
std::array<int, 3> off = { 0, 0, 0 };
|
||||
if ( slotValid( f ) ) unpackSlot( fipslot[f], off[0], off[1], off[2] );
|
||||
return off;
|
||||
};
|
||||
|
||||
size_t droppedCells = 0;
|
||||
|
||||
for ( const auto& [level, cells] : cellsByLevel )
|
||||
{
|
||||
// A level is nested when its cells' immediate parents are themselves refined-slot cells
|
||||
// (they carry a FIPNEST link of their own); primary-level parents are coarse hosts.
|
||||
size_t nestedVotes = 0;
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
const size_t parent = (size_t)fipnest[f] - 1;
|
||||
if ( fipnest[parent] != 0 ) nestedVotes++;
|
||||
}
|
||||
const bool isNested = nestedVotes >= cells.size() * 95 / 100 && nestedVotes > 0;
|
||||
|
||||
if ( isNested )
|
||||
{
|
||||
const auto paramsIt = primaryParams.find( level - 1 );
|
||||
if ( paramsIt == primaryParams.end() )
|
||||
{
|
||||
warn( QString( "Level %1 nests under a level that was not built as a primary refinement; its cells are left un-nested." )
|
||||
.arg( level ) );
|
||||
for ( size_t f : cells )
|
||||
input.oldI[f] = input.oldJ[f] = input.oldK[f] = 0;
|
||||
continue;
|
||||
}
|
||||
const LevelParams& pp = paramsIt->second;
|
||||
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
// The cell carries its immediate parent's TMP: recompute that slot's TMP from the
|
||||
// parent's own coarse host and FIPSLOT offsets.
|
||||
const size_t parent = (size_t)fipnest[f] - 1;
|
||||
if ( fipnest[parent] == 0 || !slotValid( parent ) || rootOf( parent ) < 0 )
|
||||
{
|
||||
input.oldI[f] = input.oldJ[f] = input.oldK[f] = 0;
|
||||
droppedCells++;
|
||||
continue;
|
||||
}
|
||||
const auto pr = flatIjk( (size_t)rootOf( parent ), nx, ny );
|
||||
const int pOld[3] = { pr[0] + 1, pr[1] + 1, pr[2] / (int)kFactor + 1 };
|
||||
const auto pOff = offsetsOf( parent );
|
||||
|
||||
input.tmpI[f] = ( pOld[0] - 1 ) * pp.factor[0] + pOff[0] + pp.shift[0];
|
||||
input.tmpJ[f] = ( pOld[1] - 1 ) * pp.factor[1] + pOff[1] + pp.shift[1];
|
||||
input.tmpK[f] = ( pOld[2] - 1 ) * pp.factor[2] + pOff[2] + pp.shift[2];
|
||||
}
|
||||
continue;
|
||||
}
|
||||
|
||||
// Primary level: recover the per-axis refinement factor. Within one coarse column the flat
|
||||
// coordinate advances by the band stride s per offset step, and between adjacent coarse
|
||||
// columns the column base advances by s*factor - so factor = (base stride) / s. This is what
|
||||
// distinguishes e.g. offsets {0,2} of factor 4 (hole layers) from {0,1} of factor 2.
|
||||
LevelParams params;
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
std::map<int, std::map<int, long long>> columns; // old -> off -> flat coordinate
|
||||
const int* oldArr = ( a == 0 ) ? input.oldI.data() : ( a == 1 ) ? input.oldJ.data() : input.oldK.data();
|
||||
int maxOff = 0;
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
if ( oldArr[f] < 1 || !slotValid( f ) ) continue;
|
||||
const auto off = offsetsOf( f );
|
||||
columns[oldArr[f]][off[a]] = flatIjk( f, nx, ny )[a];
|
||||
maxOff = std::max( maxOff, off[a] );
|
||||
}
|
||||
|
||||
long long s = 0;
|
||||
for ( const auto& [old, byOff] : columns )
|
||||
{
|
||||
if ( byOff.size() < 2 ) continue;
|
||||
const auto first = byOff.begin();
|
||||
const auto last = std::prev( byOff.end() );
|
||||
const long long dOff = last->first - first->first;
|
||||
const long long dFlat = last->second - first->second;
|
||||
if ( dOff > 0 && dFlat % dOff == 0 )
|
||||
{
|
||||
s = dFlat / dOff;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
long long baseStride = 0;
|
||||
if ( s > 0 )
|
||||
{
|
||||
std::map<int, long long> baseByOld;
|
||||
for ( const auto& [old, byOff] : columns )
|
||||
baseByOld[old] = byOff.begin()->second - s * byOff.begin()->first;
|
||||
for ( auto it = std::next( baseByOld.begin() ); it != baseByOld.end(); ++it )
|
||||
{
|
||||
const auto prevIt = std::prev( it );
|
||||
const long long dOld = it->first - prevIt->first;
|
||||
const long long dBase = it->second - prevIt->second;
|
||||
if ( dOld > 0 && dBase % dOld == 0 )
|
||||
{
|
||||
const long long candidate = dBase / dOld;
|
||||
if ( baseStride == 0 )
|
||||
baseStride = candidate;
|
||||
else if ( baseStride != candidate )
|
||||
{
|
||||
warn( QString( "Level %1 axis %2: inconsistent column stride; using the first one." ).arg( level ).arg( a ) );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ( s > 0 && baseStride > 0 && baseStride % s == 0 && baseStride / s > maxOff )
|
||||
{
|
||||
params.factor[a] = (int)( baseStride / s );
|
||||
}
|
||||
else
|
||||
{
|
||||
params.factor[a] = maxOff + 1;
|
||||
if ( columns.size() > 1 )
|
||||
{
|
||||
warn( QString( "Level %1 axis %2: refinement factor not recoverable from column strides; assuming %3." )
|
||||
.arg( level )
|
||||
.arg( a )
|
||||
.arg( params.factor[a] ) );
|
||||
}
|
||||
}
|
||||
params.shift[a] = nextShift[a];
|
||||
}
|
||||
|
||||
int maxTmp[3] = { 0, 0, 0 };
|
||||
for ( size_t f : cells )
|
||||
{
|
||||
if ( !slotValid( f ) )
|
||||
{
|
||||
input.oldI[f] = input.oldJ[f] = input.oldK[f] = 0;
|
||||
droppedCells++;
|
||||
continue;
|
||||
}
|
||||
const auto off = offsetsOf( f );
|
||||
const int old[3] = { input.oldI[f], input.oldJ[f], input.oldK[f] };
|
||||
const int tmp[3] = { ( old[0] - 1 ) * params.factor[0] + off[0] + params.shift[0],
|
||||
( old[1] - 1 ) * params.factor[1] + off[1] + params.shift[1],
|
||||
( old[2] - 1 ) * params.factor[2] + off[2] + params.shift[2] };
|
||||
input.tmpI[f] = tmp[0];
|
||||
input.tmpJ[f] = tmp[1];
|
||||
input.tmpK[f] = tmp[2];
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
maxTmp[a] = std::max( maxTmp[a], tmp[a] );
|
||||
}
|
||||
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
nextShift[a] = maxTmp[a] + 1;
|
||||
|
||||
primaryParams[level] = params;
|
||||
}
|
||||
|
||||
if ( droppedCells > 0 )
|
||||
{
|
||||
warn( QString( "%1 refined cells have an unusable parent link or slot and are left un-nested." ).arg( droppedCells ) );
|
||||
}
|
||||
|
||||
return input;
|
||||
}
|
||||
@@ -0,0 +1,90 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2026- Equinor ASA
|
||||
//
|
||||
// ResInsight is free software: you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation, either version 3 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
// FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
||||
// for more details.
|
||||
//
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#pragma once
|
||||
|
||||
#include "RigNestedHybridGridReconstructor.h"
|
||||
|
||||
#include <QString>
|
||||
|
||||
#include <vector>
|
||||
|
||||
//==================================================================================================
|
||||
/// Converts between the OLDIJK/TMP sidecar description of a nested hybrid grid and the compact
|
||||
/// FIPNEST/FIPSLOT/REFINE parent-child encoding intended to travel inside the INIT file (#14510).
|
||||
///
|
||||
/// Encoding (all arrays full-length, one value per flat cell in natural order):
|
||||
/// - FIPNEST : 0 for unrefined (level-1) cells; otherwise the 1-based flat natural index of the
|
||||
/// cell's IMMEDIATE parent - the coarse host cell for a primary-level cell, or the
|
||||
/// (collapsed) flat slot of the host one level shallower for a nested cell. A
|
||||
/// referenced host slot that is itself refined away carries its own FIPNEST value
|
||||
/// (its coarse host), so every chain terminates at a FIPNEST == 0 cell.
|
||||
/// - FIPSLOT : 1 + offI + 100*offJ + 10000*offK, the cell's 0-based position within its
|
||||
/// immediate parent; 0 where FIPNEST is 0.
|
||||
/// - REFINE : the per-cell nesting level, unchanged from the existing sidecar/result.
|
||||
///
|
||||
/// FIPNEST alone is not sufficient to rebuild the grid: the refinement level is not the parent
|
||||
/// chain depth (several levels can refine the coarse grid directly), and within-parent placement
|
||||
/// is ambiguous from flat coordinates when a level has hole layers. REFINE and FIPSLOT carry
|
||||
/// exactly those two missing pieces.
|
||||
///
|
||||
/// The decoder synthesizes a NestedHybridInput that reproduces the sidecar input up to per-level
|
||||
/// constant TMP shifts, which RigNestedHybridGridReconstructor::reconstruct() is invariant to -
|
||||
/// reusing the reconstructor unchanged is what guarantees that a grid imported through FIPNEST
|
||||
/// equals the grid imported through the sidecars.
|
||||
///
|
||||
/// Known prototype limitations (beyond the L5+ one above):
|
||||
/// - The decoder infers the K refinement factor as the gcd of the chain roots' flat K positions.
|
||||
/// If every refined host happens to sit on a common coarse-K stride (e.g. only every other
|
||||
/// layer hosts a refinement), the gcd overestimates the factor and the coarse K axis is
|
||||
/// compressed self-consistently; the reconstruction still places every parent at its correct
|
||||
/// flat cell, but LGR K padding can differ from the sidecar import. Storing the factor
|
||||
/// explicitly would remove the ambiguity.
|
||||
/// - Cells the encoder cannot express (unresolvedRefinedCells) stay flat on the FIPNEST path,
|
||||
/// while the sidecar path may still place them - the two imports only match exactly when the
|
||||
/// encoder reports zero unresolved cells.
|
||||
//==================================================================================================
|
||||
class RigNestedHybridGridFipnestCodec
|
||||
{
|
||||
public:
|
||||
struct ParentChildArrays
|
||||
{
|
||||
std::vector<int> fipnest;
|
||||
std::vector<int> fipslot;
|
||||
size_t unresolvedRefinedCells = 0;
|
||||
};
|
||||
|
||||
// Compute the FIPNEST/FIPSLOT arrays from a sidecar input, mirroring the reconstructor's own
|
||||
// level classification and parent resolution. Refined cells whose parent cannot be encoded are
|
||||
// left at 0 and counted in unresolvedRefinedCells.
|
||||
static ParentChildArrays
|
||||
computeParentChildArrays( const RigNestedHybridGridReconstructor::NestedHybridInput& input, size_t nx, size_t ny, size_t nz );
|
||||
|
||||
// Synthesize a reconstruction input from the FIPNEST/FIPSLOT/REFINE arrays. Appends any
|
||||
// non-fatal observations to warnings (one per line) if given.
|
||||
static RigNestedHybridGridReconstructor::NestedHybridInput buildInputFromParentChildArrays( const std::vector<int>& fipnest,
|
||||
const std::vector<int>& fipslot,
|
||||
const std::vector<int>& refine,
|
||||
size_t nx,
|
||||
size_t ny,
|
||||
size_t nz,
|
||||
QString* warnings = nullptr );
|
||||
|
||||
static int packSlot( int offI, int offJ, int offK );
|
||||
static void unpackSlot( int slot, int& offI, int& offJ, int& offK );
|
||||
};
|
||||
@@ -21,8 +21,10 @@
|
||||
#include "RiaDefines.h"
|
||||
#include "RiaLogging.h"
|
||||
#include "RiaResultNames.h"
|
||||
#include "RiaStringEncodingTools.h"
|
||||
|
||||
#include "RifEclipseKeywordContent.h"
|
||||
#include "RifEclipseOutputFileTools.h"
|
||||
#include "RifEclipseTextFileReader.h"
|
||||
#include "RifInputPropertyLoader.h"
|
||||
|
||||
@@ -32,15 +34,19 @@
|
||||
#include "RigEclipseCaseData.h"
|
||||
#include "RigLocalGrid.h"
|
||||
#include "RigMainGrid.h"
|
||||
#include "RigNestedHybridGridFipnestCodec.h"
|
||||
#include "RigNestedHybridGridReconstructor.h"
|
||||
#include "RigTypeSafeIndex.h"
|
||||
|
||||
#include "RimEclipseInputProperty.h"
|
||||
#include "RimEclipseInputPropertyCollection.h"
|
||||
|
||||
#include "ert/ecl/ecl_file.hpp"
|
||||
|
||||
#include <QDir>
|
||||
#include <QFile>
|
||||
#include <QFileInfo>
|
||||
#include <QTextStream>
|
||||
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
@@ -96,6 +102,89 @@ QString RigNestedHybridGridResultTools::oldIjkSidecarFilePath( const QString& gr
|
||||
return {};
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Nested hybrid grid: the compact parent-child encoding (FIPNEST/FIPSLOT/REFINE, see
|
||||
/// RigNestedHybridGridFipnestCodec) in a sidecar GRDECL file named "<grid-basename>_FIPNEST.grdecl".
|
||||
/// Returns its path if it exists.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
QString RigNestedHybridGridResultTools::fipnestSidecarFilePath( const QString& gridFileName )
|
||||
{
|
||||
QFileInfo gridFileInfo( gridFileName );
|
||||
if ( !gridFileInfo.exists() ) return {};
|
||||
|
||||
QDir dir = gridFileInfo.absoluteDir();
|
||||
const QString baseName = gridFileInfo.completeBaseName();
|
||||
|
||||
const QStringList candidates = { baseName + "_FIPNEST.grdecl", baseName + "_FIPNEST.GRDECL" };
|
||||
for ( const QString& candidate : candidates )
|
||||
{
|
||||
QString path = dir.absoluteFilePath( candidate );
|
||||
if ( QFile::exists( path ) ) return path;
|
||||
}
|
||||
|
||||
return {};
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Write integer keywords to a GRDECL text file, run-length encoding repeated values as "N*V".
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool RigNestedHybridGridResultTools::writeIntKeywordsToGrdeclFile( const QString& filePath,
|
||||
const std::vector<std::pair<QString, const std::vector<int>*>>& keywords )
|
||||
{
|
||||
QFile file( filePath );
|
||||
if ( !file.open( QIODevice::WriteOnly | QIODevice::Text ) ) return false;
|
||||
|
||||
QTextStream out( &file );
|
||||
out << "-- Nested hybrid grid parent-child arrays exported by ResInsight (#14510)\n";
|
||||
|
||||
const int maxTokensPerLine = 12;
|
||||
|
||||
for ( const auto& [keyword, values] : keywords )
|
||||
{
|
||||
if ( !values ) continue;
|
||||
|
||||
out << keyword << "\n";
|
||||
|
||||
int tokensOnLine = 0;
|
||||
auto emitToken = [&]( const QString& token )
|
||||
{
|
||||
out << ' ' << token;
|
||||
if ( ++tokensOnLine == maxTokensPerLine )
|
||||
{
|
||||
out << "\n";
|
||||
tokensOnLine = 0;
|
||||
}
|
||||
};
|
||||
|
||||
size_t i = 0;
|
||||
while ( i < values->size() )
|
||||
{
|
||||
const int v = ( *values )[i];
|
||||
size_t run = 1;
|
||||
while ( i + run < values->size() && ( *values )[i + run] == v )
|
||||
run++;
|
||||
|
||||
// Only run-length encode non-negative values; "N*-V" is not universally parsed.
|
||||
if ( run >= 4 && v >= 0 )
|
||||
{
|
||||
emitToken( QString( "%1*%2" ).arg( run ).arg( v ) );
|
||||
}
|
||||
else
|
||||
{
|
||||
for ( size_t r = 0; r < run; r++ )
|
||||
emitToken( QString::number( v ) );
|
||||
}
|
||||
i += run;
|
||||
}
|
||||
|
||||
if ( tokensOnLine != 0 ) out << "\n";
|
||||
out << " /\n";
|
||||
}
|
||||
|
||||
out.flush();
|
||||
return out.status() == QTextStream::Ok && file.error() == QFile::NoError;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
@@ -143,6 +232,60 @@ void RigNestedHybridGridResultTools::importOldIjkSidecarIfPresent( const QString
|
||||
RifInputPropertyLoader::loadAndSynchronizeInputProperties( inputPropertyCollection, eclipseCaseData, std::vector<QString>{ sidecarPath }, false );
|
||||
}
|
||||
|
||||
namespace
|
||||
{
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Auto-export the compact FIPNEST/FIPSLOT/REFINE parent-child encoding next to the grid file after
|
||||
/// a successful sidecar-based reconstruction (#14510), unless the file already exists. The arrays
|
||||
/// are computed from the same sidecar input the reconstruction used.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void exportFipnestSidecarIfAbsent( const QString& gridFileName,
|
||||
const RigNestedHybridGridReconstructor::NestedHybridInput& input,
|
||||
const RigEclipseCaseData* eclipseCaseData )
|
||||
{
|
||||
const QString existing = RigNestedHybridGridResultTools::fipnestSidecarFilePath( gridFileName );
|
||||
if ( !existing.isEmpty() )
|
||||
{
|
||||
RiaLogging::info( QString( "Nested hybrid grid: FIPNEST sidecar already present: %1" ).arg( existing ).toStdString() );
|
||||
return;
|
||||
}
|
||||
|
||||
const RigMainGrid* mainGrid = eclipseCaseData->mainGrid();
|
||||
const size_t nx = mainGrid->cellCountI();
|
||||
const size_t ny = mainGrid->cellCountJ();
|
||||
const size_t nz = mainGrid->cellCountK();
|
||||
|
||||
const auto arrays = RigNestedHybridGridFipnestCodec::computeParentChildArrays( input, nx, ny, nz );
|
||||
if ( arrays.fipnest.empty() )
|
||||
{
|
||||
RiaLogging::warning( "Nested hybrid grid: could not compute FIPNEST arrays; sidecar not exported." );
|
||||
return;
|
||||
}
|
||||
if ( arrays.unresolvedRefinedCells > 0 )
|
||||
{
|
||||
RiaLogging::warning(
|
||||
QString( "Nested hybrid grid: %1 refined cells could not be encoded in FIPNEST." ).arg( arrays.unresolvedRefinedCells ).toStdString() );
|
||||
}
|
||||
|
||||
// Note: RifEclipseTextFileReader parses GRDECL values as float, so FIPNEST indices written here
|
||||
// stay exact only up to 2^24 cells (~16.7M); the INIT-embedded INTE path has no such limit.
|
||||
QFileInfo gridFileInfo( gridFileName );
|
||||
const QString path = gridFileInfo.absoluteDir().absoluteFilePath( gridFileInfo.completeBaseName() + "_FIPNEST.grdecl" );
|
||||
|
||||
const std::vector<std::pair<QString, const std::vector<int>*>> keywords = { { "FIPNEST", &arrays.fipnest },
|
||||
{ "FIPSLOT", &arrays.fipslot },
|
||||
{ RiaResultNames::refine(), &input.refine } };
|
||||
if ( RigNestedHybridGridResultTools::writeIntKeywordsToGrdeclFile( path, keywords ) )
|
||||
{
|
||||
RiaLogging::info( QString( "Nested hybrid grid: exported FIPNEST sidecar to %1" ).arg( path ).toStdString() );
|
||||
}
|
||||
else
|
||||
{
|
||||
RiaLogging::warning( QString( "Nested hybrid grid: failed to write FIPNEST sidecar %1" ).arg( path ).toStdString() );
|
||||
}
|
||||
}
|
||||
} // namespace
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
@@ -184,12 +327,127 @@ void RigNestedHybridGridResultTools::reconstructNestedHybridGridIfPresent( const
|
||||
input.tmpK = readIntKeyword( oldIjkContent, "TMPK" );
|
||||
|
||||
QString errorMessage;
|
||||
RigNestedHybridGridReconstructor::reconstruct( eclipseCaseData, input, &errorMessage );
|
||||
if ( RigNestedHybridGridReconstructor::reconstruct( eclipseCaseData, input, &errorMessage ) )
|
||||
{
|
||||
exportFipnestSidecarIfAbsent( gridFileName, input, eclipseCaseData );
|
||||
}
|
||||
|
||||
// The caller computes grid caches (search tree, faults, NNCs) once, after this reconstruction, so
|
||||
// that the expensive geometric passes run on the clean grid rather than the flat overlapping one.
|
||||
}
|
||||
|
||||
namespace
|
||||
{
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Open the INIT file next to the grid file, or nullptr if there is none. The caller owns the handle.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
ecl_file_type* openInitFileNextToGrid( const QString& gridFileName, QString* initFileNameOut = nullptr )
|
||||
{
|
||||
QStringList fileSet;
|
||||
if ( !RifEclipseOutputFileTools::findSiblingFilesWithSameBaseName( gridFileName, &fileSet ) ) return nullptr;
|
||||
|
||||
const QString initFileName = RifEclipseOutputFileTools::firstFileNameOfType( fileSet, ECL_INIT_FILE );
|
||||
if ( initFileName.isEmpty() ) return nullptr;
|
||||
|
||||
if ( initFileNameOut ) *initFileNameOut = initFileName;
|
||||
return ecl_file_open( RiaStringEncodingTools::toNativeEncoded( initFileName ).data(), ECL_FILE_CLOSE_STREAM );
|
||||
}
|
||||
} // namespace
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool RigNestedHybridGridResultTools::initFileHasFipnest( const QString& gridFileName )
|
||||
{
|
||||
ecl_file_type* initFile = openInitFileNextToGrid( gridFileName );
|
||||
if ( !initFile ) return false;
|
||||
|
||||
const bool hasFipnest = ecl_file_has_kw( initFile, "FIPNEST" );
|
||||
ecl_file_close( initFile );
|
||||
|
||||
return hasFipnest;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Reconstruct the nested hybrid grid LGR hierarchy from the FIPNEST/FIPSLOT/REFINE arrays embedded
|
||||
/// in the INIT file - the sidecar-free import path prototyped in #14510. The REFINE levels are
|
||||
/// stored as an input-property result before the reconstruction so the per-level result helpers
|
||||
/// (and the LGR-cell extension) see them, mirroring the sidecar path.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool RigNestedHybridGridResultTools::reconstructNestedHybridGridFromInitFile( const QString& gridFileName, RigEclipseCaseData* eclipseCaseData )
|
||||
{
|
||||
if ( !eclipseCaseData || !eclipseCaseData->mainGrid() ) return false;
|
||||
|
||||
QString initFileName;
|
||||
ecl_file_type* initFile = openInitFileNextToGrid( gridFileName, &initFileName );
|
||||
if ( !initFile ) return false;
|
||||
|
||||
auto readIntKeyword = [&]( const char* keyword, std::vector<int>* values )
|
||||
{ return ecl_file_has_kw( initFile, keyword ) && RifEclipseOutputFileTools::keywordData( initFile, keyword, 0, values ); };
|
||||
|
||||
std::vector<int> fipnest, fipslot, refine;
|
||||
const bool haveFipnest = readIntKeyword( "FIPNEST", &fipnest );
|
||||
const bool haveAll = haveFipnest && readIntKeyword( "FIPSLOT", &fipslot ) &&
|
||||
readIntKeyword( RiaResultNames::refine().toLatin1().data(), &refine );
|
||||
ecl_file_close( initFile );
|
||||
|
||||
if ( !haveFipnest ) return false;
|
||||
if ( !haveAll )
|
||||
{
|
||||
RiaLogging::warning( QString( "Nested hybrid grid: %1 contains FIPNEST but not its FIPSLOT/REFINE companions; "
|
||||
"falling back to sidecar import." )
|
||||
.arg( initFileName )
|
||||
.toStdString() );
|
||||
return false;
|
||||
}
|
||||
|
||||
const RigMainGrid* mainGrid = eclipseCaseData->mainGrid();
|
||||
const size_t nx = mainGrid->cellCountI();
|
||||
const size_t ny = mainGrid->cellCountJ();
|
||||
const size_t nz = mainGrid->cellCountK();
|
||||
const size_t cellCount = nx * ny * nz;
|
||||
if ( fipnest.size() != cellCount || fipslot.size() != cellCount || refine.size() != cellCount )
|
||||
{
|
||||
RiaLogging::warning( QString( "Nested hybrid grid: FIPNEST/FIPSLOT/REFINE in %1 do not cover all %2 cells; "
|
||||
"falling back to sidecar import." )
|
||||
.arg( initFileName )
|
||||
.arg( cellCount )
|
||||
.toStdString() );
|
||||
return false;
|
||||
}
|
||||
|
||||
RiaLogging::info(
|
||||
QString( "Nested hybrid grid: reconstructing from the FIPNEST/FIPSLOT/REFINE arrays in %1" ).arg( initFileName ).toStdString() );
|
||||
|
||||
QString warnings;
|
||||
const auto input = RigNestedHybridGridFipnestCodec::buildInputFromParentChildArrays( fipnest, fipslot, refine, nx, ny, nz, &warnings );
|
||||
for ( const QString& line : warnings.split( '\n', Qt::SkipEmptyParts ) )
|
||||
{
|
||||
RiaLogging::warning( ( "Nested hybrid grid: " + line ).toStdString() );
|
||||
}
|
||||
|
||||
// Make the refinement levels available as the REFINE result (used by the per-level aggregation
|
||||
// helpers) unless the sidecar import already provided it. Created before the reconstruction so
|
||||
// the full-length array is extended to cover the new LGR cells.
|
||||
if ( RigCaseCellResultsData* results = eclipseCaseData->results( RiaDefines::PorosityModelType::MATRIX_MODEL ) )
|
||||
{
|
||||
RigEclipseResultAddress refineAddr( RiaDefines::ResultCatType::INPUT_PROPERTY,
|
||||
RiaDefines::ResultDataType::INTEGER,
|
||||
RiaResultNames::refine() );
|
||||
if ( !results->hasResultEntry( refineAddr ) )
|
||||
{
|
||||
results->createResultEntry( refineAddr, false );
|
||||
if ( auto* timesteps = results->modifiableCellScalarResultTimesteps( refineAddr ) )
|
||||
{
|
||||
timesteps->push_back( std::vector<double>( refine.begin(), refine.end() ) );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
QString errorMessage;
|
||||
return RigNestedHybridGridReconstructor::reconstruct( eclipseCaseData, input, &errorMessage );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
|
||||
@@ -22,6 +22,7 @@
|
||||
|
||||
#include <QString>
|
||||
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
class RigCaseCellResultsData;
|
||||
@@ -46,6 +47,22 @@ public:
|
||||
// Path to the "<grid-basename>_OLDIJK.grdecl" sidecar next to the grid file, or empty if absent.
|
||||
static QString oldIjkSidecarFilePath( const QString& gridFileName );
|
||||
|
||||
// Path to the "<grid-basename>_FIPNEST.grdecl" parent-child sidecar, or empty if absent.
|
||||
static QString fipnestSidecarFilePath( const QString& gridFileName );
|
||||
|
||||
// Write integer keywords to a GRDECL text file (run-length encoded). Used for the auto-exported
|
||||
// FIPNEST/FIPSLOT/REFINE parent-child sidecar (#14510).
|
||||
static bool writeIntKeywordsToGrdeclFile( const QString& filePath,
|
||||
const std::vector<std::pair<QString, const std::vector<int>*>>& keywords );
|
||||
|
||||
// True if the INIT file next to the grid file contains the FIPNEST keyword.
|
||||
static bool initFileHasFipnest( const QString& gridFileName );
|
||||
|
||||
// Reconstruct the nested hybrid grid from the FIPNEST/FIPSLOT/REFINE arrays embedded in the INIT
|
||||
// file (#14510). Returns false if the arrays are absent or unusable, so the caller can fall back
|
||||
// to the sidecar-based path.
|
||||
static bool reconstructNestedHybridGridFromInitFile( const QString& gridFileName, RigEclipseCaseData* eclipseCaseData );
|
||||
|
||||
static void importRefineSidecarIfPresent( const QString& gridFileName,
|
||||
RimEclipseInputPropertyCollection* inputPropertyCollection,
|
||||
RigEclipseCaseData* eclipseCaseData );
|
||||
|
||||
@@ -29,6 +29,7 @@ set(SOURCE_UNITTEST_FILES
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigGridExportAdapter-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigNonUniformRefinement-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigNestedHybridGridReconstructor-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigNestedHybridGridFipnest-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigStatisticsMath-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigStimPlanModelTools-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigWellPathIntersectionTools-Test.cpp
|
||||
|
||||
@@ -0,0 +1,539 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2026- Equinor ASA
|
||||
//
|
||||
// ResInsight is free software: you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation, either version 3 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
// FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
||||
// for more details.
|
||||
//
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#include "gtest/gtest.h"
|
||||
|
||||
#include "RiaStringEncodingTools.h"
|
||||
#include "RiaTestDataDirectory.h"
|
||||
|
||||
#include "RifEclipseKeywordContent.h"
|
||||
#include "RifEclipseTextFileReader.h"
|
||||
#include "RifReaderEclipseOutput.h"
|
||||
|
||||
#include "RigEclipseCaseData.h"
|
||||
#include "RigLocalGrid.h"
|
||||
#include "RigMainGrid.h"
|
||||
#include "RigNestedHybridGridFipnestCodec.h"
|
||||
#include "RigNestedHybridGridReconstructor.h"
|
||||
#include "RigNestedHybridGridResultTools.h"
|
||||
|
||||
#include "RimEclipseResultCase.h"
|
||||
|
||||
#include <ert/ecl/ecl_endian_flip.h>
|
||||
#include <ert/ecl/ecl_kw.h>
|
||||
#include <ert/ecl/fortio.h>
|
||||
|
||||
#include <QDir>
|
||||
#include <QFile>
|
||||
#include <QTemporaryDir>
|
||||
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <map>
|
||||
#include <memory>
|
||||
#include <set>
|
||||
|
||||
//==================================================================================================
|
||||
// Tests for the FIPNEST/FIPSLOT/REFINE parent-child encoding of a nested hybrid grid (#14510):
|
||||
// the codec that derives the arrays from the OLDIJK/TMP sidecar description and back, the GRDECL
|
||||
// sidecar auto-export, and the sidecar-free import path that reads the arrays from the INIT file.
|
||||
// The requirement backing most of these tests: a grid reconstructed through FIPNEST must equal the
|
||||
// grid reconstructed through the sidecars.
|
||||
//==================================================================================================
|
||||
|
||||
namespace
|
||||
{
|
||||
const size_t NX = 150, NY = 84, NZ = 96;
|
||||
|
||||
QString nestedHybridModelDir()
|
||||
{
|
||||
QDir baseFolder( TEST_MODEL_DIR );
|
||||
baseFolder.cd( "NestedHybridGrid" );
|
||||
return baseFolder.absolutePath();
|
||||
}
|
||||
|
||||
// Read all integer keywords of a GRDECL file in a single pass.
|
||||
std::map<QString, std::vector<int>> readIntKeywords( const QString& filePath )
|
||||
{
|
||||
std::map<QString, std::vector<int>> result;
|
||||
auto content = RifEclipseTextFileReader::readKeywordAndValues( filePath.toStdString() );
|
||||
for ( const auto& kw : content )
|
||||
{
|
||||
std::vector<int> values;
|
||||
values.reserve( kw.values.size() );
|
||||
for ( float v : kw.values )
|
||||
values.push_back( (int)std::lround( v ) );
|
||||
result[QString::fromStdString( kw.keyword )] = std::move( values );
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
RigNestedHybridGridReconstructor::NestedHybridInput buildSidecarInput( const QString& dir )
|
||||
{
|
||||
RigNestedHybridGridReconstructor::NestedHybridInput input;
|
||||
auto refine = readIntKeywords( dir + "/DROGON_NESTED_REFINE.grdecl" );
|
||||
auto oldIjk = readIntKeywords( dir + "/DROGON_NESTED_OLDIJK.grdecl" );
|
||||
input.refine = refine["REFINE"];
|
||||
input.oldI = oldIjk["OLDI"];
|
||||
input.oldJ = oldIjk["OLDJ"];
|
||||
input.oldK = oldIjk["OLDK"];
|
||||
input.tmpI = oldIjk["TMPI"];
|
||||
input.tmpJ = oldIjk["TMPJ"];
|
||||
input.tmpK = oldIjk["TMPK"];
|
||||
return input;
|
||||
}
|
||||
|
||||
// A freshly opened (not yet reconstructed) DROGON_NESTED case.
|
||||
struct OpenedCase
|
||||
{
|
||||
std::unique_ptr<RimEclipseResultCase> resultCase;
|
||||
cvf::ref<RigEclipseCaseData> caseData;
|
||||
cvf::ref<RifReaderEclipseOutput> reader;
|
||||
};
|
||||
|
||||
bool openCase( const QString& gridFile, OpenedCase& openedCase )
|
||||
{
|
||||
openedCase.resultCase = std::make_unique<RimEclipseResultCase>();
|
||||
openedCase.caseData = new RigEclipseCaseData( openedCase.resultCase.get() );
|
||||
openedCase.reader = new RifReaderEclipseOutput;
|
||||
if ( !openedCase.reader->open( gridFile, openedCase.caseData.p() ) ) return false;
|
||||
openedCase.caseData->mainGrid()->computeCachedData();
|
||||
return true;
|
||||
}
|
||||
|
||||
// Compare the reconstructed LGR structure of two cases: same grids (name, dimensions, cell range,
|
||||
// parent), the same LGR-cell -> source-flat-cell mapping, and coarse-parent maps that partition the
|
||||
// refined cells identically (the numeric coarse indices may legitimately differ, see the codec).
|
||||
void expectEqualReconstruction( RigMainGrid* a, RigMainGrid* b )
|
||||
{
|
||||
ASSERT_EQ( a->gridCount(), b->gridCount() );
|
||||
ASSERT_EQ( a->totalCellCount(), b->totalCellCount() );
|
||||
|
||||
for ( size_t i = 1; i < a->gridCount(); i++ )
|
||||
{
|
||||
RigLocalGrid* lgrA = dynamic_cast<RigLocalGrid*>( a->gridByIndex( i ) );
|
||||
RigLocalGrid* lgrB = dynamic_cast<RigLocalGrid*>( b->gridByIndex( i ) );
|
||||
ASSERT_TRUE( lgrA != nullptr && lgrB != nullptr );
|
||||
EXPECT_EQ( lgrA->gridName(), lgrB->gridName() );
|
||||
EXPECT_EQ( lgrA->cellCountI(), lgrB->cellCountI() );
|
||||
EXPECT_EQ( lgrA->cellCountJ(), lgrB->cellCountJ() );
|
||||
EXPECT_EQ( lgrA->cellCountK(), lgrB->cellCountK() );
|
||||
EXPECT_EQ( lgrA->reservoirCellIndex( 0 ), lgrB->reservoirCellIndex( 0 ) );
|
||||
EXPECT_EQ( lgrA->parentGrid()->gridName(), lgrB->parentGrid()->gridName() );
|
||||
}
|
||||
|
||||
EXPECT_EQ( a->nestedHybridLgrSourceCells(), b->nestedHybridLgrSourceCells() );
|
||||
|
||||
// Coarse-parent maps: same key set, and the values partition the keys identically.
|
||||
const std::map<size_t, size_t>& cpA = a->nestedHybridCoarseParents();
|
||||
const std::map<size_t, size_t>& cpB = b->nestedHybridCoarseParents();
|
||||
ASSERT_EQ( cpA.size(), cpB.size() );
|
||||
std::map<size_t, size_t> valueMap; // A-parent-value -> B-parent-value
|
||||
for ( const auto& [flatCell, parentA] : cpA )
|
||||
{
|
||||
auto itB = cpB.find( flatCell );
|
||||
ASSERT_TRUE( itB != cpB.end() ) << "coarse-parent key sets differ at flat cell " << flatCell;
|
||||
auto [mapped, inserted] = valueMap.try_emplace( parentA, itB->second );
|
||||
EXPECT_EQ( mapped->second, itB->second ) << "coarse-parent partitions differ at flat cell " << flatCell;
|
||||
}
|
||||
std::set<size_t> mappedBValues;
|
||||
for ( const auto& [va, vb] : valueMap )
|
||||
mappedBValues.insert( vb );
|
||||
EXPECT_EQ( mappedBValues.size(), valueMap.size() ) << "coarse-parent value mapping is not a bijection";
|
||||
}
|
||||
|
||||
// Append one INTE keyword to an existing unformatted Eclipse output file (e.g. an INIT file).
|
||||
bool appendIntKeywordToEclipseFile( const QString& filePath, const char* keyword, const std::vector<int>& values )
|
||||
{
|
||||
fortio_type* fortio = fortio_open_append( RiaStringEncodingTools::toNativeEncoded( filePath ).data(), false, ECL_ENDIAN_FLIP );
|
||||
if ( !fortio ) return false;
|
||||
|
||||
ecl_kw_type* kw = ecl_kw_alloc_new( keyword, (int)values.size(), ECL_INT, values.data() );
|
||||
const bool ok = ecl_kw_fwrite( kw, fortio );
|
||||
ecl_kw_free( kw );
|
||||
fortio_fclose( fortio );
|
||||
return ok;
|
||||
}
|
||||
} // namespace
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Shared fixture: the sidecar input arrays and the codec-derived FIPNEST/FIPSLOT arrays are
|
||||
/// computed once (array-only work, no grid).
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
class RigNestedHybridGridFipnestTest : public ::testing::Test
|
||||
{
|
||||
protected:
|
||||
static void SetUpTestSuite()
|
||||
{
|
||||
const QString dir = nestedHybridModelDir();
|
||||
if ( !QFile::exists( dir + "/DROGON_NESTED.EGRID" ) )
|
||||
{
|
||||
s_setupError = "Missing test model: " + dir + "/DROGON_NESTED.EGRID";
|
||||
return;
|
||||
}
|
||||
|
||||
s_input = buildSidecarInput( dir );
|
||||
if ( s_input.refine.size() != NX * NY * NZ )
|
||||
{
|
||||
s_setupError = "Unexpected sidecar size";
|
||||
return;
|
||||
}
|
||||
|
||||
s_arrays = RigNestedHybridGridFipnestCodec::computeParentChildArrays( s_input, NX, NY, NZ );
|
||||
s_setupOk = true;
|
||||
}
|
||||
|
||||
static void TearDownTestSuite()
|
||||
{
|
||||
s_input = {};
|
||||
s_arrays = {};
|
||||
s_setupOk = false;
|
||||
}
|
||||
|
||||
void SetUp() override { ASSERT_TRUE( s_setupOk ) << s_setupError.toStdString(); }
|
||||
|
||||
static RigNestedHybridGridReconstructor::NestedHybridInput s_input;
|
||||
static RigNestedHybridGridFipnestCodec::ParentChildArrays s_arrays;
|
||||
static bool s_setupOk;
|
||||
static QString s_setupError;
|
||||
};
|
||||
|
||||
RigNestedHybridGridReconstructor::NestedHybridInput RigNestedHybridGridFipnestTest::s_input;
|
||||
RigNestedHybridGridFipnestCodec::ParentChildArrays RigNestedHybridGridFipnestTest::s_arrays;
|
||||
bool RigNestedHybridGridFipnestTest::s_setupOk = false;
|
||||
QString RigNestedHybridGridFipnestTest::s_setupError;
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// Synthetic round-trip with a K refinement factor of 2 where the refined hosts stop below the top
|
||||
/// coarse-K layer. The reconstructor re-derives the K factor from the maximum OLDK of its input, so
|
||||
/// the decoder must emit OLD for the unrefined cells too - otherwise the factor comes out as 4 here
|
||||
/// and every parent lands on the wrong K layer.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigNestedHybridGridFipnestCodecTest, SyntheticKFactorRoundTrip )
|
||||
{
|
||||
// Flat 8x1x8 grid; coarse 4x1x4 embedded at (i, 0, k*2); one level-2 region refining coarse
|
||||
// cells (1..2, 1, 1..2) by (2,1,2), its cells appended as a band at i 4..7, k 0..3.
|
||||
const size_t nx = 8, ny = 1, nz = 8;
|
||||
|
||||
RigNestedHybridGridReconstructor::NestedHybridInput input;
|
||||
input.refine.assign( nx * ny * nz, 1 );
|
||||
input.oldI.assign( nx * ny * nz, 1 );
|
||||
input.oldJ.assign( nx * ny * nz, 1 );
|
||||
input.oldK.assign( nx * ny * nz, 1 );
|
||||
input.tmpI.assign( nx * ny * nz, 0 );
|
||||
input.tmpJ.assign( nx * ny * nz, 0 );
|
||||
input.tmpK.assign( nx * ny * nz, 0 );
|
||||
|
||||
for ( size_t k = 0; k < nz; k++ )
|
||||
for ( size_t i = 0; i < 4; i++ )
|
||||
{
|
||||
const size_t f = i + k * nx * ny;
|
||||
input.oldI[f] = (int)i + 1;
|
||||
input.oldK[f] = (int)k / 2 + 1;
|
||||
}
|
||||
|
||||
for ( int tk = 0; tk < 4; tk++ )
|
||||
for ( int ti = 0; ti < 4; ti++ )
|
||||
{
|
||||
const size_t f = (size_t)( ti + 4 ) + (size_t)tk * nx * ny;
|
||||
input.refine[f] = 2;
|
||||
input.oldI[f] = ti / 2 + 1;
|
||||
input.oldK[f] = tk / 2 + 1;
|
||||
input.tmpI[f] = ti;
|
||||
input.tmpK[f] = tk;
|
||||
}
|
||||
|
||||
const auto arrays = RigNestedHybridGridFipnestCodec::computeParentChildArrays( input, nx, ny, nz );
|
||||
ASSERT_EQ( arrays.fipnest.size(), nx * ny * nz );
|
||||
EXPECT_EQ( arrays.unresolvedRefinedCells, 0u );
|
||||
|
||||
QString warnings;
|
||||
const auto translated =
|
||||
RigNestedHybridGridFipnestCodec::buildInputFromParentChildArrays( arrays.fipnest, arrays.fipslot, input.refine, nx, ny, nz, &warnings );
|
||||
EXPECT_TRUE( warnings.isEmpty() ) << warnings.toStdString();
|
||||
ASSERT_EQ( translated.refine.size(), nx * ny * nz );
|
||||
|
||||
// The coarse host (OLD) of every refined cell must be recovered exactly, and the K factor the
|
||||
// reconstructor re-derives (NZ / max OLDK) must match the true factor of 2.
|
||||
int maxOldK = 0;
|
||||
for ( size_t f = 0; f < nx * ny * nz; f++ )
|
||||
{
|
||||
maxOldK = std::max( maxOldK, translated.oldK[f] );
|
||||
if ( input.refine[f] <= 1 ) continue;
|
||||
EXPECT_EQ( translated.oldI[f], input.oldI[f] );
|
||||
EXPECT_EQ( translated.oldJ[f], input.oldJ[f] );
|
||||
EXPECT_EQ( translated.oldK[f], input.oldK[f] );
|
||||
}
|
||||
ASSERT_GT( maxOldK, 0 );
|
||||
ASSERT_EQ( nz % (size_t)maxOldK, 0u );
|
||||
EXPECT_EQ( nz / (size_t)maxOldK, 2u );
|
||||
|
||||
// Relative TMP must be reproduced cellwise (the reconstructor is shift-invariant).
|
||||
int minOrig[3] = { INT_MAX, INT_MAX, INT_MAX }, minTrans[3] = { INT_MAX, INT_MAX, INT_MAX };
|
||||
for ( size_t f = 0; f < nx * ny * nz; f++ )
|
||||
{
|
||||
if ( input.refine[f] <= 1 ) continue;
|
||||
const int orig[3] = { input.tmpI[f], input.tmpJ[f], input.tmpK[f] };
|
||||
const int trans[3] = { translated.tmpI[f], translated.tmpJ[f], translated.tmpK[f] };
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
minOrig[a] = std::min( minOrig[a], orig[a] );
|
||||
minTrans[a] = std::min( minTrans[a], trans[a] );
|
||||
}
|
||||
}
|
||||
for ( size_t f = 0; f < nx * ny * nz; f++ )
|
||||
{
|
||||
if ( input.refine[f] <= 1 ) continue;
|
||||
EXPECT_EQ( input.tmpI[f] - minOrig[0], translated.tmpI[f] - minTrans[0] );
|
||||
EXPECT_EQ( input.tmpJ[f] - minOrig[1], translated.tmpJ[f] - minTrans[1] );
|
||||
EXPECT_EQ( input.tmpK[f] - minOrig[2], translated.tmpK[f] - minTrans[2] );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// T2: structural invariants of the encoded arrays on the DROGON_NESTED dataset.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST_F( RigNestedHybridGridFipnestTest, ArrayInvariants )
|
||||
{
|
||||
const size_t cellCount = NX * NY * NZ;
|
||||
ASSERT_EQ( s_arrays.fipnest.size(), cellCount );
|
||||
ASSERT_EQ( s_arrays.fipslot.size(), cellCount );
|
||||
EXPECT_EQ( s_arrays.unresolvedRefinedCells, 0u );
|
||||
|
||||
std::map<int, size_t> refinedByLevel;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
if ( s_input.refine[f] > 1 ) refinedByLevel[s_input.refine[f]]++;
|
||||
EXPECT_EQ( refinedByLevel[2], 4768u );
|
||||
EXPECT_EQ( refinedByLevel[3], 74708u );
|
||||
EXPECT_EQ( refinedByLevel[4], 5984u );
|
||||
|
||||
size_t virtualSlots = 0;
|
||||
std::map<int, size_t> level4ChildrenByParent;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
const int link = s_arrays.fipnest[f];
|
||||
if ( s_input.refine[f] > 1 )
|
||||
{
|
||||
// Every refined cell has a parent link and a slot, and its chain ends at an unrefined
|
||||
// cell within a few steps.
|
||||
ASSERT_GT( link, 0 ) << "refined cell " << f << " has no parent link";
|
||||
ASSERT_GT( s_arrays.fipslot[f], 0 );
|
||||
|
||||
size_t cur = f;
|
||||
int steps = 0;
|
||||
while ( s_arrays.fipnest[cur] != 0 )
|
||||
{
|
||||
ASSERT_LE( ++steps, 4 );
|
||||
cur = (size_t)s_arrays.fipnest[cur] - 1;
|
||||
ASSERT_LT( cur, cellCount );
|
||||
}
|
||||
EXPECT_EQ( s_input.refine[cur], 1 );
|
||||
|
||||
if ( s_input.refine[f] == 4 ) level4ChildrenByParent[link]++;
|
||||
}
|
||||
else if ( link != 0 )
|
||||
{
|
||||
// An unrefined cell with a link is a refined-away host slot: it must point straight at
|
||||
// an unrefined coarse cell with no further link.
|
||||
virtualSlots++;
|
||||
EXPECT_EQ( s_input.refine[(size_t)link - 1], 1 );
|
||||
EXPECT_EQ( s_arrays.fipnest[(size_t)link - 1], 0 );
|
||||
EXPECT_GT( s_arrays.fipslot[f], 0 );
|
||||
}
|
||||
}
|
||||
|
||||
// The 5984 level-4 cells nest in 748 level-3 host slots, 8 children (2x2x2) each; every one of
|
||||
// those hosts was refined away and is represented by a virtual slot.
|
||||
EXPECT_EQ( level4ChildrenByParent.size(), 748u );
|
||||
for ( const auto& [parent, count] : level4ChildrenByParent )
|
||||
EXPECT_EQ( count, 8u );
|
||||
EXPECT_EQ( virtualSlots, 748u );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// T1: the translator reproduces the sidecar input from the encoded arrays - REFINE and OLD exactly,
|
||||
/// TMP up to the per-level per-axis constant shift the reconstructor is invariant to.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST_F( RigNestedHybridGridFipnestTest, CodecRoundTripReproducesSidecarInput )
|
||||
{
|
||||
QString warnings;
|
||||
const auto translated =
|
||||
RigNestedHybridGridFipnestCodec::buildInputFromParentChildArrays( s_arrays.fipnest, s_arrays.fipslot, s_input.refine, NX, NY, NZ, &warnings );
|
||||
EXPECT_TRUE( warnings.isEmpty() ) << warnings.toStdString();
|
||||
|
||||
const size_t cellCount = NX * NY * NZ;
|
||||
ASSERT_EQ( translated.refine.size(), cellCount );
|
||||
EXPECT_EQ( translated.refine, s_input.refine );
|
||||
|
||||
// OLD (the coarse host) must match exactly on every refined cell.
|
||||
size_t oldMismatches = 0;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
if ( s_input.refine[f] <= 1 ) continue;
|
||||
if ( translated.oldI[f] != s_input.oldI[f] || translated.oldJ[f] != s_input.oldJ[f] || translated.oldK[f] != s_input.oldK[f] )
|
||||
oldMismatches++;
|
||||
}
|
||||
EXPECT_EQ( oldMismatches, 0u );
|
||||
|
||||
// TMP: per level and axis, the level-relative coordinates must match cellwise.
|
||||
std::map<int, std::array<int, 6>> tmpMins; // level -> min original tmp[3], min translated tmp[3]
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
if ( s_input.refine[f] <= 1 ) continue;
|
||||
const int level = s_input.refine[f];
|
||||
auto it = tmpMins.find( level );
|
||||
const int orig[3] = { s_input.tmpI[f], s_input.tmpJ[f], s_input.tmpK[f] };
|
||||
const int trans[3] = { translated.tmpI[f], translated.tmpJ[f], translated.tmpK[f] };
|
||||
if ( it == tmpMins.end() )
|
||||
{
|
||||
tmpMins[level] = { orig[0], orig[1], orig[2], trans[0], trans[1], trans[2] };
|
||||
}
|
||||
else
|
||||
{
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
{
|
||||
it->second[a] = std::min( it->second[a], orig[a] );
|
||||
it->second[3 + a] = std::min( it->second[3 + a], trans[a] );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
size_t tmpMismatches = 0;
|
||||
for ( size_t f = 0; f < cellCount; f++ )
|
||||
{
|
||||
if ( s_input.refine[f] <= 1 ) continue;
|
||||
const auto& mins = tmpMins[s_input.refine[f]];
|
||||
const int orig[3] = { s_input.tmpI[f], s_input.tmpJ[f], s_input.tmpK[f] };
|
||||
const int trans[3] = { translated.tmpI[f], translated.tmpJ[f], translated.tmpK[f] };
|
||||
for ( int a = 0; a < 3; a++ )
|
||||
if ( orig[a] - mins[a] != trans[a] - mins[3 + a] ) tmpMismatches++;
|
||||
}
|
||||
EXPECT_EQ( tmpMismatches, 0u );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// T4: the GRDECL writer round-trips through the standard text file reader, including run-length
|
||||
/// encoded stretches and large values.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST_F( RigNestedHybridGridFipnestTest, GrdeclWriterRoundTrip )
|
||||
{
|
||||
QTemporaryDir tempDir;
|
||||
ASSERT_TRUE( tempDir.isValid() );
|
||||
const QString path = tempDir.filePath( "roundtrip_FIPNEST.grdecl" );
|
||||
|
||||
std::vector<int> a = { 0, 0, 0, 0, 0, 7, 7, 7, 7, 1209600, 3, 3, 1, 2 };
|
||||
std::vector<int> b( 1000, 42 );
|
||||
b[999] = -5;
|
||||
|
||||
const std::vector<std::pair<QString, const std::vector<int>*>> keywords = { { "FIPNEST", &a }, { "FIPSLOT", &b } };
|
||||
ASSERT_TRUE( RigNestedHybridGridResultTools::writeIntKeywordsToGrdeclFile( path, keywords ) );
|
||||
|
||||
auto readBack = readIntKeywords( path );
|
||||
EXPECT_EQ( readBack["FIPNEST"], a );
|
||||
EXPECT_EQ( readBack["FIPSLOT"], b );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// T5: the stock DROGON_NESTED INIT file has no FIPNEST, so the INIT-based path must not trigger.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST_F( RigNestedHybridGridFipnestTest, StockInitFileHasNoFipnest )
|
||||
{
|
||||
const QString gridFile = nestedHybridModelDir() + "/DROGON_NESTED.EGRID";
|
||||
EXPECT_FALSE( RigNestedHybridGridResultTools::initFileHasFipnest( gridFile ) );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// T3: a grid reconstructed from the translated FIPNEST arrays equals the grid reconstructed from
|
||||
/// the sidecar input.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST_F( RigNestedHybridGridFipnestTest, ReconstructionFromFipnestEqualsSidecarReconstruction )
|
||||
{
|
||||
const QString gridFile = nestedHybridModelDir() + "/DROGON_NESTED.EGRID";
|
||||
|
||||
OpenedCase sidecarCase;
|
||||
ASSERT_TRUE( openCase( gridFile, sidecarCase ) );
|
||||
QString err;
|
||||
ASSERT_TRUE( RigNestedHybridGridReconstructor::reconstruct( sidecarCase.caseData.p(), s_input, &err ) ) << err.toStdString();
|
||||
|
||||
QString warnings;
|
||||
const auto translated =
|
||||
RigNestedHybridGridFipnestCodec::buildInputFromParentChildArrays( s_arrays.fipnest, s_arrays.fipslot, s_input.refine, NX, NY, NZ, &warnings );
|
||||
|
||||
OpenedCase fipnestCase;
|
||||
ASSERT_TRUE( openCase( gridFile, fipnestCase ) );
|
||||
ASSERT_TRUE( RigNestedHybridGridReconstructor::reconstruct( fipnestCase.caseData.p(), translated, &err ) ) << err.toStdString();
|
||||
|
||||
expectEqualReconstruction( sidecarCase.caseData->mainGrid(), fipnestCase.caseData->mainGrid() );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// T6: hermetic end-to-end pass through the real files. Phase 1 stages the model with sidecars in a
|
||||
/// temp dir and verifies the sidecar-based reconstruction auto-exports the FIPNEST GRDECL. Phase 2
|
||||
/// stages EGRID + INIT only (no sidecars) in a second dir, inserts the FIPNEST/FIPSLOT/REFINE INTE
|
||||
/// arrays into the INIT copy, imports through the INIT-based path and compares against the
|
||||
/// sidecar-based reconstruction.
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST_F( RigNestedHybridGridFipnestTest, EndToEndInitFileImport )
|
||||
{
|
||||
const QString dir = nestedHybridModelDir();
|
||||
|
||||
// Phase 1: sidecar import in a staged dir -> auto-exported FIPNEST sidecar.
|
||||
QTemporaryDir stagedSidecarDir;
|
||||
ASSERT_TRUE( stagedSidecarDir.isValid() );
|
||||
for ( const QString& f : { QString( "DROGON_NESTED.EGRID" ),
|
||||
QString( "DROGON_NESTED.INIT" ),
|
||||
QString( "DROGON_NESTED_REFINE.grdecl" ),
|
||||
QString( "DROGON_NESTED_OLDIJK.grdecl" ) } )
|
||||
{
|
||||
ASSERT_TRUE( QFile::copy( dir + "/" + f, stagedSidecarDir.filePath( f ) ) );
|
||||
}
|
||||
const QString stagedGridFile = stagedSidecarDir.filePath( "DROGON_NESTED.EGRID" );
|
||||
|
||||
OpenedCase sidecarCase;
|
||||
ASSERT_TRUE( openCase( stagedGridFile, sidecarCase ) );
|
||||
RigNestedHybridGridResultTools::reconstructNestedHybridGridIfPresent( stagedGridFile, sidecarCase.caseData.p() );
|
||||
ASSERT_GT( sidecarCase.caseData->mainGrid()->gridCount(), 1u );
|
||||
|
||||
const QString exportedFipnest = RigNestedHybridGridResultTools::fipnestSidecarFilePath( stagedGridFile );
|
||||
ASSERT_FALSE( exportedFipnest.isEmpty() ) << "sidecar import did not auto-export the FIPNEST GRDECL";
|
||||
|
||||
auto exportedArrays = readIntKeywords( exportedFipnest );
|
||||
ASSERT_EQ( exportedArrays["FIPNEST"].size(), NX * NY * NZ );
|
||||
ASSERT_EQ( exportedArrays["FIPSLOT"].size(), NX * NY * NZ );
|
||||
ASSERT_EQ( exportedArrays["REFINE"].size(), NX * NY * NZ );
|
||||
|
||||
// Phase 2: EGRID + augmented INIT only - no sidecars anywhere near the grid file.
|
||||
QTemporaryDir stagedInitDir;
|
||||
ASSERT_TRUE( stagedInitDir.isValid() );
|
||||
ASSERT_TRUE( QFile::copy( dir + "/DROGON_NESTED.EGRID", stagedInitDir.filePath( "DROGON_NESTED.EGRID" ) ) );
|
||||
ASSERT_TRUE( QFile::copy( dir + "/DROGON_NESTED.INIT", stagedInitDir.filePath( "DROGON_NESTED.INIT" ) ) );
|
||||
|
||||
const QString augmentedInit = stagedInitDir.filePath( "DROGON_NESTED.INIT" );
|
||||
const QString fipnestGridFile = stagedInitDir.filePath( "DROGON_NESTED.EGRID" );
|
||||
ASSERT_TRUE( appendIntKeywordToEclipseFile( augmentedInit, "FIPNEST", exportedArrays["FIPNEST"] ) );
|
||||
ASSERT_TRUE( appendIntKeywordToEclipseFile( augmentedInit, "FIPSLOT", exportedArrays["FIPSLOT"] ) );
|
||||
ASSERT_TRUE( appendIntKeywordToEclipseFile( augmentedInit, "REFINE", exportedArrays["REFINE"] ) );
|
||||
|
||||
EXPECT_TRUE( RigNestedHybridGridResultTools::initFileHasFipnest( fipnestGridFile ) );
|
||||
|
||||
OpenedCase fipnestCase;
|
||||
ASSERT_TRUE( openCase( fipnestGridFile, fipnestCase ) );
|
||||
ASSERT_TRUE( RigNestedHybridGridResultTools::reconstructNestedHybridGridFromInitFile( fipnestGridFile, fipnestCase.caseData.p() ) );
|
||||
|
||||
expectEqualReconstruction( sidecarCase.caseData->mainGrid(), fipnestCase.caseData->mainGrid() );
|
||||
}
|
||||
@@ -0,0 +1,40 @@
|
||||
# Nested hybrid grid: FIPNEST INIT-file tooling (#14510)
|
||||
|
||||
ResInsight can import a "nested hybrid grid" (a single flat EGRID with the refined cells
|
||||
appended in per-level I bands) and rebuild it as a true LGR hierarchy. The parent-child
|
||||
description normally comes from the `<gridbase>_REFINE.grdecl` + `<gridbase>_OLDIJK.grdecl`
|
||||
sidecar files. Issue #14510 prototypes a compact alternative that travels inside the INIT
|
||||
file instead, as three full-length (NX*NY*NZ) INTE arrays:
|
||||
|
||||
- `FIPNEST` — per cell, the 1-based flat natural index of the immediate parent cell
|
||||
(0 for unrefined cells). Refined-away host slots carry their own link to their coarse
|
||||
host, so every chain terminates in an unrefined cell.
|
||||
- `FIPSLOT` — `1 + offI + 100*offJ + 10000*offK`, the cell's position within its parent.
|
||||
- `REFINE` — the per-cell nesting level (unchanged from the existing sidecar).
|
||||
|
||||
FIPNEST alone is not sufficient: the refinement level is not the parent-chain depth
|
||||
(several levels can refine the coarse grid directly), and within-parent placement is
|
||||
ambiguous from flat coordinates when a level has hole layers. See
|
||||
`RigNestedHybridGridFipnestCodec` for the full reasoning.
|
||||
|
||||
## Workflow
|
||||
|
||||
1. Import the grid into ResInsight once with the REFINE/OLDIJK sidecars present. On a
|
||||
successful reconstruction ResInsight auto-exports `<gridbase>_FIPNEST.grdecl` next to
|
||||
the grid file (skipped if the file already exists).
|
||||
2. Insert the arrays into a copy of the INIT file:
|
||||
|
||||
```bash
|
||||
python insert_fipnest_into_init.py DROGON_NESTED_FIPNEST.grdecl \
|
||||
DROGON_NESTED.INIT staged/DROGON_NESTED.INIT
|
||||
```
|
||||
|
||||
Requires `resfo` and `numpy` (`pip install resfo`).
|
||||
3. Place the EGRID and the augmented INIT in a directory *without* the sidecar files and
|
||||
open it in ResInsight. The log shows
|
||||
`Nested hybrid grid: reconstructing from the FIPNEST/FIPSLOT/REFINE arrays in ...` and
|
||||
the reconstructed LGR hierarchy is identical to the sidecar-based import (covered by
|
||||
the `RigNestedHybridGridFipnestTest` unit tests).
|
||||
|
||||
When both the sidecars and a FIPNEST-augmented INIT are present, the INIT arrays take
|
||||
precedence and the sidecars remain as the fallback.
|
||||
@@ -0,0 +1,140 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Insert the nested-hybrid-grid parent-child arrays into an Eclipse INIT file.
|
||||
|
||||
Reads the FIPNEST/FIPSLOT/REFINE keywords from a GRDECL file (typically the
|
||||
"<gridbase>_FIPNEST.grdecl" sidecar that ResInsight auto-exports when importing a
|
||||
nested hybrid grid through the OLDIJK sidecars, see issue #14510) and appends them
|
||||
as INTE arrays to a copy of an existing INIT file. ResInsight detects FIPNEST in
|
||||
the INIT file and reconstructs the nested hybrid grid without any sidecar files.
|
||||
|
||||
Example:
|
||||
python insert_fipnest_into_init.py DROGON_NESTED_FIPNEST.grdecl \\
|
||||
DROGON_NESTED.INIT /tmp/staged/DROGON_NESTED.INIT
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
import resfo
|
||||
|
||||
DEFAULT_KEYWORDS = ["FIPNEST", "FIPSLOT", "REFINE"]
|
||||
|
||||
|
||||
def read_grdecl_int_keywords(path, wanted):
|
||||
"""Parse integer keywords from a GRDECL text file.
|
||||
|
||||
Handles '--' comments and 'N*V' run-length encoded values. Returns a dict
|
||||
keyword -> list of ints for the keywords in `wanted` that are present.
|
||||
"""
|
||||
wanted = {kw.upper() for kw in wanted}
|
||||
result = {}
|
||||
current = None
|
||||
values = []
|
||||
|
||||
with open(path, "r", encoding="ascii", errors="replace") as f:
|
||||
for line in f:
|
||||
comment = line.find("--")
|
||||
if comment >= 0:
|
||||
line = line[:comment]
|
||||
|
||||
for token in line.split():
|
||||
if current is None:
|
||||
if token.upper() in wanted:
|
||||
current = token.upper()
|
||||
values = []
|
||||
continue
|
||||
|
||||
if token == "/" or token.endswith("/"):
|
||||
token = token[:-1]
|
||||
if token:
|
||||
values.extend(_expand_token(token))
|
||||
result[current] = values
|
||||
current = None
|
||||
continue
|
||||
|
||||
values.extend(_expand_token(token))
|
||||
|
||||
if current is not None:
|
||||
raise ValueError(f"keyword {current} in {path} is not terminated with '/'")
|
||||
return result
|
||||
|
||||
|
||||
def _expand_token(token):
|
||||
if "*" in token:
|
||||
count, value = token.split("*", 1)
|
||||
return [int(value)] * int(count)
|
||||
return [int(token)]
|
||||
|
||||
|
||||
def grid_cell_count_from_init(contents):
|
||||
"""NX*NY*NZ from the INTEHEAD array (items 9-11, 1-based) of an INIT file."""
|
||||
for kw, array in contents:
|
||||
if kw.strip() == "INTEHEAD":
|
||||
return int(array[8]) * int(array[9]) * int(array[10])
|
||||
return None
|
||||
|
||||
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(
|
||||
description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
|
||||
)
|
||||
parser.add_argument(
|
||||
"grdecl", help="GRDECL file holding the FIPNEST/FIPSLOT/REFINE keywords"
|
||||
)
|
||||
parser.add_argument("init_in", help="existing INIT file to augment")
|
||||
parser.add_argument(
|
||||
"init_out", help="augmented INIT file to write (must differ from init_in)"
|
||||
)
|
||||
parser.add_argument(
|
||||
"--keywords",
|
||||
nargs="+",
|
||||
default=DEFAULT_KEYWORDS,
|
||||
help=f"keywords to transfer (default: {' '.join(DEFAULT_KEYWORDS)})",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--formatted",
|
||||
action="store_true",
|
||||
help="write a formatted (ASCII) INIT file instead of an unformatted one",
|
||||
)
|
||||
args = parser.parse_args()
|
||||
|
||||
if args.init_in == args.init_out:
|
||||
parser.error("init_out must differ from init_in")
|
||||
|
||||
arrays = read_grdecl_int_keywords(args.grdecl, args.keywords)
|
||||
missing = [kw for kw in args.keywords if kw.upper() not in arrays]
|
||||
if missing:
|
||||
sys.exit(f"error: keyword(s) {', '.join(missing)} not found in {args.grdecl}")
|
||||
|
||||
contents = [(kw, array) for kw, array in resfo.read(args.init_in)]
|
||||
|
||||
cell_count = grid_cell_count_from_init(contents)
|
||||
if cell_count is not None:
|
||||
for kw in args.keywords:
|
||||
n = len(arrays[kw.upper()])
|
||||
if n != cell_count:
|
||||
sys.exit(
|
||||
f"error: {kw} has {n} values, the INIT grid has {cell_count} cells (NX*NY*NZ)"
|
||||
)
|
||||
|
||||
# Replace any existing occurrences, then append in the requested order.
|
||||
upper_keywords = {kw.upper() for kw in args.keywords}
|
||||
contents = [
|
||||
(kw, array) for kw, array in contents if kw.strip() not in upper_keywords
|
||||
]
|
||||
for kw in args.keywords:
|
||||
contents.append(
|
||||
(kw.upper().ljust(8), np.asarray(arrays[kw.upper()], dtype=np.int32))
|
||||
)
|
||||
|
||||
file_format = resfo.Format.FORMATTED if args.formatted else resfo.Format.UNFORMATTED
|
||||
resfo.write(args.init_out, contents, fileformat=file_format)
|
||||
|
||||
print(
|
||||
f"wrote {args.init_out} with {', '.join(kw.upper() for kw in args.keywords)} ({len(contents)} arrays total)"
|
||||
)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in new issue
Block a user