diff --git a/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.cpp b/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.cpp index 0edd15ce86..07c08023aa 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.cpp @@ -38,9 +38,8 @@ #include #include #include -#include #include -#include +#include namespace { @@ -48,16 +47,6 @@ 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; } - -// A reconstructed cell, recorded so that a deeper (nested) refinement level can find it as a parent. -// Keyed by the cell's TMP: a nested cell carries its immediate parent's TMP, so a child looks up its -// own TMP among the cells one level shallower. -struct BuiltRef -{ - RigGridBase* grid = nullptr; - size_t localCell = 0; - int level = 0; -}; } // namespace //-------------------------------------------------------------------------------------------------- @@ -85,9 +74,6 @@ bool RigNestedHybridGridReconstructor::reconstruct( RigEclipseCaseData* caseData { return setError( "OLDI/OLDJ/OLDK size does not match main grid cell count." ); } - const bool haveTmp = input.tmpI.size() == cellCount && input.tmpJ.size() == cellCount && input.tmpK.size() == cellCount; - if ( !haveTmp ) return setError( "TMPI/TMPJ/TMPK size does not match main grid cell count." ); - if ( grid->gridCount() > 1 ) { return setError( "Grid already contains local grids; skipping reconstruction." ); @@ -130,154 +116,118 @@ bool RigNestedHybridGridReconstructor::reconstruct( RigEclipseCaseData* caseData const size_t origCellCount = grid->totalCellCount(); // before any LGR cells are appended std::map sourceCells; // LGR global cell index -> source flat cell index int nextGridId = (int)grid->gridCount(); // main grid is 0 - size_t deferred = 0; - // Index every reconstructed cell by its TMP, so a deeper level can find its immediate parent. - int tmpMax[3] = { 0, 0, 0 }; - for ( const auto& [level, cells] : cellsByLevel ) - for ( size_t f : cells ) - { - tmpMax[0] = std::max( tmpMax[0], input.tmpI[f] ); - tmpMax[1] = std::max( tmpMax[1], input.tmpJ[f] ); - tmpMax[2] = std::max( tmpMax[2], input.tmpK[f] ); - } - const size_t tmpDimX = (size_t)tmpMax[0] + 1; - const size_t tmpDimY = (size_t)tmpMax[1] + 1; - auto tmpKey = [&]( size_t f ) { return naturalIndex( input.tmpI[f], input.tmpJ[f], input.tmpK[f], tmpDimX, tmpDimY ); }; - - std::map> tmpToBuilt; - - // Record the cells appended since beforeTotal (the just-built level), keyed by their TMP. - auto recordBuilt = [&]( size_t beforeTotal, int level ) - { - for ( size_t gc = beforeTotal; gc < grid->totalCellCount(); gc++ ) - { - auto it = sourceCells.find( gc ); - if ( it == sourceCells.end() ) continue; // hole - size_t localCell = 0; - RigGridBase* cellGrid = grid->gridAndGridLocalIdxFromGlobalCellIdx( gc, &localCell ); - tmpToBuilt[tmpKey( it->second )].push_back( { cellGrid, localCell, level } ); - } - }; + using CoarseIjk = std::array; + using Factor = std::array; - // Build refinement levels from shallow to deep. For each level, first try to resolve every cell's - // immediate parent among the already-built cells one level shallower (true nesting, any depth); - // otherwise treat the level as a direct, uniform refinement of the coarse grid. Cells whose parent - // cannot be resolved unambiguously are deferred (left as flat cells) rather than mis-nested. + // Build every refinement level directly below the main grid. OLDIJK identifies the coarse parent, + // while the distinct flat-grid coordinates within each parent determine the local cell ordering. + // Adjacent parents with equal refinement dimensions are combined into one regular LGR. for ( const auto& [level, cells] : cellsByLevel ) { - // Try to resolve each cell's immediate parent (a unique built cell one level shallower whose - // TMP equals this cell's TMP). - std::map> cellToParent; // flat cell -> (parent grid, parent local cell) + std::map> cellsByParent; for ( size_t f : cells ) + cellsByParent[{ input.oldI[f], input.oldJ[f], input.oldK[f] }].push_back( f ); + + std::map, 3>> coordinatesByParent; + std::map> parentsByFactor; + for ( const auto& [parent, parentCells] : cellsByParent ) { - auto it = tmpToBuilt.find( tmpKey( f ) ); - if ( it == tmpToBuilt.end() ) continue; - const BuiltRef* unique = nullptr; - for ( const BuiltRef& b : it->second ) + auto& coordinates = coordinatesByParent[parent]; + for ( size_t f : parentCells ) { - if ( b.level != level - 1 ) continue; - if ( unique ) // more than one candidate -> ambiguous - { - unique = nullptr; - break; - } - unique = &b; + const std::array flatIjk = { f % nx, ( f / nx ) % ny, f / ( nx * ny ) }; + for ( int axis = 0; axis < 3; axis++ ) + coordinates[axis].push_back( flatIjk[axis] ); } - if ( unique ) cellToParent[f] = { unique->grid, unique->localCell }; - } - const size_t beforeTotal = grid->totalCellCount(); - - if ( cellToParent.size() >= cells.size() * 95 / 100 && !cellToParent.empty() ) - { - // Nested level: build LGR(s) inside the resolved parent grid(s). - deferred += cells.size() - cellToParent.size(); // unresolved cells are left as flat cells - buildNestedLevel( caseData, cells, input, cellToParent, nextGridId, sourceCells ); - recordBuilt( beforeTotal, level ); - continue; - } - - // Primary level: a direct, uniform refinement of the coarse grid. - 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++ ) + Factor factor; + for ( int axis = 0; axis < 3; axis++ ) { - 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] ); + auto& values = coordinates[axis]; + std::sort( values.begin(), values.end() ); + values.erase( std::unique( values.begin(), values.end() ), values.end() ); + factor[axis] = values.size(); } + parentsByFactor[factor].insert( parent ); } - const int coarseDim[3] = { c1[0] - c0[0] + 1, c1[1] - c0[1] + 1, c1[2] - c0[2] + 1 }; - const int tmpDim[3] = { t1[0] - t0[0] + 1, t1[1] - t0[1] + 1, t1[2] - t0[2] + 1 }; - - int factor[3]; - for ( int a = 0; a < 3; a++ ) - factor[a] = std::max( 1, (int)std::lround( (double)tmpDim[a] / (double)coarseDim[a] ) ); - - // Verify the level is a uniform refinement of the coarse grid: the coarse cell implied by each - // cell's TMP (via the factor) must equal its OLD index. If not, the level cannot be placed and - // is deferred (left as flat cells). - 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 ) + int componentIndex = 0; + for ( const auto& [factor, parents] : parentsByFactor ) { - deferred += cells.size(); - continue; - } - - const cvf::Vec3st dims( (size_t)coarseDim[0] * factor[0], (size_t)coarseDim[1] * factor[1], (size_t)coarseDim[2] * factor[2] ); + std::set remaining = parents; + while ( !remaining.empty() ) + { + std::vector component; + std::vector stack = { *remaining.begin() }; + remaining.erase( stack.front() ); - // Map each LGR local cell to its parent (collapsed) coarse cell in the flat grid, and to its - // source flat cell (UNDEFINED = hole). - std::vector boxToParent( dims.x() * dims.y() * dims.z(), cvf::UNDEFINED_SIZE_T ); - for ( size_t lk = 0; lk < dims.z(); lk++ ) - for ( size_t lj = 0; lj < dims.y(); lj++ ) - for ( size_t li = 0; li < dims.x(); li++ ) + while ( !stack.empty() ) { - size_t local = naturalIndex( li, lj, lk, dims.x(), dims.y() ); - size_t oi = (size_t)c0[0] + li / factor[0]; - size_t oj = (size_t)c0[1] + lj / factor[1]; - size_t ok = (size_t)c0[2] + lk / factor[2]; - boxToParent[local] = naturalIndex( oi - 1, oj - 1, ( ok - 1 ) * kFactor, nx, ny ); + const CoarseIjk parent = stack.back(); + stack.pop_back(); + component.push_back( parent ); + + for ( int axis = 0; axis < 3; axis++ ) + for ( int direction : { -1, 1 } ) + { + CoarseIjk neighbour = parent; + neighbour[axis] += direction; + auto it = remaining.find( neighbour ); + if ( it == remaining.end() ) continue; + stack.push_back( *it ); + remaining.erase( it ); + } } - std::vector boxToFlat( dims.x() * dims.y() * dims.z(), cvf::UNDEFINED_SIZE_T ); - for ( size_t f : cells ) - { - size_t li = input.tmpI[f] - t0[0]; - size_t lj = input.tmpJ[f] - t0[1]; - size_t lk = input.tmpK[f] - t0[2]; - size_t local = naturalIndex( li, lj, lk, dims.x(), dims.y() ); - boxToFlat[local] = f; - } - - const QString gridName = QString( "LGR_NHG_L%1" ).arg( level ); - RigLocalGrid* levelGrid = buildLocalGrid( caseData, grid, nextGridId++, gridName, dims, boxToFlat, boxToParent, sourceCells ); + CoarseIjk c0 = { INT_MAX, INT_MAX, INT_MAX }; + CoarseIjk c1 = { 0, 0, 0 }; + for ( const CoarseIjk& parent : component ) + for ( int axis = 0; axis < 3; axis++ ) + { + c0[axis] = std::min( c0[axis], parent[axis] ); + c1[axis] = std::max( c1[axis], parent[axis] ); + } - // Record EVERY cell of this primary level (including holes), keyed by its TMP slot, so a deeper - // level can resolve its parent here - the parent is often a hole cell (one that was further - // refined and so has no geometry of its own until a child supplies it). - for ( size_t lk = 0; lk < dims.z(); lk++ ) - for ( size_t lj = 0; lj < dims.y(); lj++ ) - for ( size_t li = 0; li < dims.x(); li++ ) + const cvf::Vec3st dims( (size_t)( c1[0] - c0[0] + 1 ) * factor[0], + (size_t)( c1[1] - c0[1] + 1 ) * factor[1], + (size_t)( c1[2] - c0[2] + 1 ) * factor[2] ); + std::vector boxToParent( dims.x() * dims.y() * dims.z(), cvf::UNDEFINED_SIZE_T ); + for ( size_t lk = 0; lk < dims.z(); lk++ ) + for ( size_t lj = 0; lj < dims.y(); lj++ ) + for ( size_t li = 0; li < dims.x(); li++ ) + { + size_t local = naturalIndex( li, lj, lk, dims.x(), dims.y() ); + size_t oi = (size_t)c0[0] + li / factor[0]; + size_t oj = (size_t)c0[1] + lj / factor[1]; + size_t ok = (size_t)c0[2] + lk / factor[2]; + boxToParent[local] = naturalIndex( oi - 1, oj - 1, ( ok - 1 ) * kFactor, nx, ny ); + } + + std::vector boxToFlat( dims.x() * dims.y() * dims.z(), cvf::UNDEFINED_SIZE_T ); + for ( const CoarseIjk& parent : component ) { - size_t local = naturalIndex( li, lj, lk, dims.x(), dims.y() ); - size_t ti = (size_t)t0[0] + li, tj = (size_t)t0[1] + lj, tk = (size_t)t0[2] + lk; - tmpToBuilt[naturalIndex( ti, tj, tk, tmpDimX, tmpDimY )].push_back( { levelGrid, local, level } ); + const auto& coordinates = coordinatesByParent[parent]; + for ( size_t f : cellsByParent[parent] ) + { + const std::array flatIjk = { f % nx, ( f / nx ) % ny, f / ( nx * ny ) }; + size_t localIjk[3]; + for ( int axis = 0; axis < 3; axis++ ) + { + const auto coordinateIt = std::lower_bound( coordinates[axis].begin(), coordinates[axis].end(), flatIjk[axis] ); + localIjk[axis] = (size_t)( parent[axis] - c0[axis] ) * factor[axis] + + (size_t)( coordinateIt - coordinates[axis].begin() ); + } + boxToFlat[naturalIndex( localIjk[0], localIjk[1], localIjk[2], dims.x(), dims.y() )] = f; + } } + + QString gridName = QString( "LGR_NHG_L%1" ).arg( level ); + if ( componentIndex > 0 ) gridName += QString( "_%1" ).arg( componentIndex ); + componentIndex++; + buildLocalGrid( caseData, grid, nextGridId++, gridName, dims, boxToFlat, boxToParent, sourceCells ); + } + } } if ( sourceCells.empty() ) return setError( "No refined regions could be reconstructed." ); @@ -298,11 +248,8 @@ bool RigNestedHybridGridReconstructor::reconstruct( RigEclipseCaseData* caseData // The grid count changed - invalidate per-grid caches that were sized for the flat grid. caseData->clearWellCellsInGridCache(); - RiaLogging::info( QString( "Nested hybrid grid: created %1 LGRs (%2 cells)%3" ) - .arg( grid->gridCount() - 1 ) - .arg( sourceCells.size() ) - .arg( deferred > 0 ? QString( ", %1 nested-level cells left un-nested" ).arg( deferred ) : QString() ) - .toStdString() ); + RiaLogging::info( + QString( "Nested hybrid grid: created %1 LGRs (%2 cells)" ).arg( grid->gridCount() - 1 ).arg( sourceCells.size() ).toStdString() ); return true; } @@ -321,8 +268,7 @@ RigLocalGrid* RigNestedHybridGridReconstructor::buildLocalGrid( RigEclipseCaseDa const cvf::Vec3st& dims, const std::vector& boxToFlat, const std::vector& boxToParent, - std::map& sourceCells, - bool synthesizeHoleParentGeometry ) + std::map& sourceCells ) { RigMainGrid* grid = caseData->mainGrid(); @@ -347,9 +293,6 @@ RigLocalGrid* RigNestedHybridGridReconstructor::buildLocalGrid( RigEclipseCaseDa grid->nodes().resize( nodeStart + lgrCellCount * 8, cvf::Vec3d( 0, 0, 0 ) ); } - // Track the child cells of each refined parent (global indices), for parent geometry synthesis. - std::map> childrenByParentGlobal; - for ( size_t glc = 0; glc < lgrCellCount; glc++ ) { RigCell& lgrCell = grid->cell( cellStart + glc ); @@ -388,7 +331,6 @@ RigLocalGrid* RigNestedHybridGridReconstructor::buildLocalGrid( RigEclipseCaseDa if ( parentGlobal != cvf::UNDEFINED_SIZE_T ) { grid->cell( parentGlobal ).setSubGrid( localGrid ); - childrenByParentGlobal[parentGlobal].push_back( cellStart + glc ); } sourceCells[cellStart + glc] = flat; @@ -396,214 +338,9 @@ RigLocalGrid* RigNestedHybridGridReconstructor::buildLocalGrid( RigEclipseCaseDa localGrid->setParentGrid( parentGrid ); - // For nested LGRs the parent cell may itself be a hole in its (parent) grid - i.e. a coarser cell - // that was subdivided and therefore has no geometry. Give it an axis-aligned bounding-box geometry - // enclosing its children so picking / bounding boxes behave. - if ( synthesizeHoleParentGeometry ) - { - for ( const auto& [parentGlobal, children] : childrenByParentGlobal ) - { - RigCell& parentCell = grid->cell( parentGlobal ); - if ( !parentCell.isInvalid() ) continue; - - cvf::Vec3d mn( HUGE_VAL, HUGE_VAL, HUGE_VAL ), mx( -HUGE_VAL, -HUGE_VAL, -HUGE_VAL ); - for ( size_t child : children ) - { - for ( size_t c = 0; c < 8; c++ ) - { - const cvf::Vec3d& v = grid->nodes()[grid->cell( child ).cornerIndices()[c]]; - mn.x() = std::min( mn.x(), v.x() ); - mn.y() = std::min( mn.y(), v.y() ); - mn.z() = std::min( mn.z(), v.z() ); - mx.x() = std::max( mx.x(), v.x() ); - mx.y() = std::max( mx.y(), v.y() ); - mx.z() = std::max( mx.z(), v.z() ); - } - } - - const std::array corners = { cvf::Vec3d( mn.x(), mn.y(), mn.z() ), - cvf::Vec3d( mx.x(), mn.y(), mn.z() ), - cvf::Vec3d( mx.x(), mx.y(), mn.z() ), - cvf::Vec3d( mn.x(), mx.y(), mn.z() ), - cvf::Vec3d( mn.x(), mn.y(), mx.z() ), - cvf::Vec3d( mx.x(), mn.y(), mx.z() ), - cvf::Vec3d( mx.x(), mx.y(), mx.z() ), - cvf::Vec3d( mn.x(), mx.y(), mx.z() ) }; - - const size_t base = grid->nodes().size(); - grid->nodes().resize( base + 8, cvf::Vec3d( 0, 0, 0 ) ); - for ( size_t c = 0; c < 8; c++ ) - { - grid->nodes()[base + c] = corners[c]; - parentCell.cornerIndices()[c] = base + c; - } - parentCell.setInvalid( false ); - } - } - return localGrid; } -//-------------------------------------------------------------------------------------------------- -/// Build the nested LGR(s) for a level that refines another (parent) level. Each cell's parent cell is -/// (TMP - parentTmpOrigin) and its position within that parent cell comes from the flat-IJK sub-block. -/// The cells are merged into connected regions in the refined parent coordinate space, so a stack of -/// refined parent cells (e.g. across K) becomes a single LGR rather than one LGR per parent cell. One -/// compact LGR is built per connected region, placed inside parentGrid. Returns the number of cells -/// that could not be nested. -//-------------------------------------------------------------------------------------------------- -void RigNestedHybridGridReconstructor::buildNestedLevel( RigEclipseCaseData* caseData, - const std::vector& cells, - const NestedHybridInput& input, - const std::map>& cellToParent, - int& nextGridId, - std::map& sourceCells ) -{ - if ( cells.empty() || cellToParent.empty() ) return; - - RigMainGrid* grid = caseData->mainGrid(); - const size_t nx = grid->cellCountI(); - const size_t ny = grid->cellCountJ(); - const int level = input.refine[cells.front()]; - - // Group the resolved cells by their parent grid (a nested level may refine cells in several parent - // LGRs); each LGR can only refine cells of a single parent grid. - std::map> byParentGrid; - for ( size_t f : cells ) - { - auto it = cellToParent.find( f ); - if ( it != cellToParent.end() ) byParentGrid[it->second.first].push_back( f ); - } - - int component = 0; - for ( const auto& [parentGrid, gcells] : byParentGrid ) - { - const size_t pDimX = parentGrid->cellCountI(); - const size_t pDimY = parentGrid->cellCountJ(); - - // Parent-grid-local IJK of a cell's parent cell. - auto parentIjk = [&]( size_t f, size_t ijk[3] ) - { - size_t pl = cellToParent.at( f ).second; - ijk[0] = pl % pDimX; - ijk[1] = ( pl / pDimX ) % pDimY; - ijk[2] = pl / ( pDimX * pDimY ); - }; - - // Per parent cell, the min flat-IJK of its children, and the per-axis sub-refinement factor. - std::map> flatOrigin; - std::map> flatExtent; - for ( size_t f : gcells ) - { - size_t ijk[3]; - parentIjk( f, ijk ); - size_t plocal = naturalIndex( ijk[0], ijk[1], ijk[2], pDimX, pDimY ); - std::array fxyz = { f % nx, ( f / nx ) % ny, f / ( nx * ny ) }; - auto oit = flatOrigin.find( plocal ); - if ( oit == flatOrigin.end() ) - { - flatOrigin[plocal] = fxyz; - flatExtent[plocal] = fxyz; - } - else - { - for ( int a = 0; a < 3; a++ ) - { - oit->second[a] = std::min( oit->second[a], fxyz[a] ); - flatExtent[plocal][a] = std::max( flatExtent[plocal][a], fxyz[a] ); - } - } - } - - size_t sub[3] = { 1, 1, 1 }; - for ( const auto& [plocal, mn] : flatOrigin ) - for ( int a = 0; a < 3; a++ ) - sub[a] = std::max( sub[a], flatExtent[plocal][a] - mn[a] + 1 ); - - // Position of each nested cell in the refined parent coordinate space: - // q = parentLocalIjk * sub + subPos (subPos = flat IJK - the parent cell's flat-IJK origin) - const size_t qDimX = pDimX * sub[0]; - const size_t qDimY = pDimY * sub[1]; - std::vector> qOfCell( gcells.size() ); - std::map qToCell; // q linear index -> index into gcells - for ( size_t idx = 0; idx < gcells.size(); idx++ ) - { - size_t f = gcells[idx]; - size_t ijk[3]; - parentIjk( f, ijk ); - size_t plocal = naturalIndex( ijk[0], ijk[1], ijk[2], pDimX, pDimY ); - const std::array& flatMin = flatOrigin[plocal]; - std::array q = { ijk[0] * sub[0] + ( f % nx ) - flatMin[0], - ijk[1] * sub[1] + ( ( f / nx ) % ny ) - flatMin[1], - ijk[2] * sub[2] + ( f / ( nx * ny ) ) - flatMin[2] }; - qOfCell[idx] = q; - qToCell[naturalIndex( q[0], q[1], q[2], qDimX, qDimY )] = idx; - } - - // Merge cells into connected regions (6-connectivity in q space) - one LGR per region. - std::vector visited( gcells.size(), false ); - for ( size_t seed = 0; seed < gcells.size(); seed++ ) - { - if ( visited[seed] ) continue; - - std::vector regionCells; - std::vector stack = { seed }; - visited[seed] = true; - while ( !stack.empty() ) - { - size_t cur = stack.back(); - stack.pop_back(); - regionCells.push_back( cur ); - - const std::array& q = qOfCell[cur]; - const int nbr[6][3] = { { 1, 0, 0 }, { -1, 0, 0 }, { 0, 1, 0 }, { 0, -1, 0 }, { 0, 0, 1 }, { 0, 0, -1 } }; - for ( const auto& d : nbr ) - { - long ni = (long)q[0] + d[0], nj = (long)q[1] + d[1], nk = (long)q[2] + d[2]; - if ( ni < 0 || nj < 0 || nk < 0 ) continue; - auto it = qToCell.find( naturalIndex( (size_t)ni, (size_t)nj, (size_t)nk, qDimX, qDimY ) ); - if ( it != qToCell.end() && !visited[it->second] ) - { - visited[it->second] = true; - stack.push_back( it->second ); - } - } - } - - // Bounding box of the region in q space. - size_t qmin[3] = { SIZE_MAX, SIZE_MAX, SIZE_MAX }, qmax[3] = { 0, 0, 0 }; - for ( size_t idx : regionCells ) - for ( int a = 0; a < 3; a++ ) - { - qmin[a] = std::min( qmin[a], qOfCell[idx][a] ); - qmax[a] = std::max( qmax[a], qOfCell[idx][a] ); - } - const cvf::Vec3st dims( qmax[0] - qmin[0] + 1, qmax[1] - qmin[1] + 1, qmax[2] - qmin[2] + 1 ); - - std::vector boxToParent( dims.x() * dims.y() * dims.z(), cvf::UNDEFINED_SIZE_T ); - for ( size_t lk = 0; lk < dims.z(); lk++ ) - for ( size_t lj = 0; lj < dims.y(); lj++ ) - for ( size_t li = 0; li < dims.x(); li++ ) - { - size_t local = naturalIndex( li, lj, lk, dims.x(), dims.y() ); - size_t qi = qmin[0] + li, qj = qmin[1] + lj, qk = qmin[2] + lk; - boxToParent[local] = naturalIndex( qi / sub[0], qj / sub[1], qk / sub[2], pDimX, pDimY ); - } - - std::vector boxToFlat( dims.x() * dims.y() * dims.z(), cvf::UNDEFINED_SIZE_T ); - for ( size_t idx : regionCells ) - { - const std::array& q = qOfCell[idx]; - size_t local = naturalIndex( q[0] - qmin[0], q[1] - qmin[1], q[2] - qmin[2], dims.x(), dims.y() ); - boxToFlat[local] = gcells[idx]; - } - - const QString gridName = QString( "LGR_NHG_L%1_%2" ).arg( level ).arg( component++ ); - buildLocalGrid( caseData, parentGrid, nextGridId++, gridName, dims, boxToFlat, boxToParent, sourceCells, true ); - } - } -} - //-------------------------------------------------------------------------------------------------- /// Register all appended LGR cells as additional active cells, for both porosity models. Each active /// LGR cell gets a NEW result index appended after the existing active cells, which preserves the diff --git a/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.h b/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.h index 332c14337f..dd85c6a9bb 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.h +++ b/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridReconstructor.h @@ -41,14 +41,15 @@ class RigLocalGrid; /// copying the real refined geometry, linking each region to its parent (coarse) cell, and hiding /// the original scattered cells. /// -/// The parent of each refined cell is provided explicitly by sidecar properties (no HOSTNUM): -/// - REFINE : per-cell nesting level (1 = unrefined base, 2/3/4 = refined levels) +/// The coarse parent and refinement level of each refined cell are provided explicitly by sidecar +/// properties (no HOSTNUM): +/// - REFINE : per-cell refinement level (1 = unrefined base, 2/3/4 = refined levels) /// - OLDI/OLDJ/OLDK : the cell's parent COARSE cell IJK (1-based) -/// - TMPI/TMPJ/TMPK : the cell's local position in refined coordinate space /// /// The coarse host cell is the collapsed flat cell at (OLDI-1, OLDJ-1, (OLDK-1)*KF), where the K -/// refinement factor KF = NZ / coarseNZ. Refinement is per-level and non-uniform, footprints may be -/// non-rectangular (holes), and level-(n+1) regions nest inside a single level-n cell. +/// refinement factor KF = NZ / coarseNZ. Each refinement level is reconstructed as one or more LGRs +/// directly below the main grid. The flat-grid IJK positions determine the ordering within each +/// coarse parent. //================================================================================================== class RigNestedHybridGridReconstructor { @@ -60,9 +61,6 @@ class RigNestedHybridGridReconstructor std::vector oldI; // 1-based coarse parent I std::vector oldJ; // 1-based coarse parent J std::vector oldK; // 1-based coarse parent K - std::vector tmpI; // local refined I - std::vector tmpJ; // local refined J - std::vector tmpK; // local refined K }; // Returns true if at least one refined region was reconstructed. @@ -81,20 +79,7 @@ class RigNestedHybridGridReconstructor const cvf::Vec3st& dims, const std::vector& boxToFlat, const std::vector& boxToParent, - std::map& sourceCells, - bool synthesizeHoleParentGeometry = false ); - - // Build the nested LGR(s) for a level that refines another (parent) level rather than the coarse - // grid (e.g. level 4 inside level 3). cellToParent maps each resolved cell (flat index) to its - // immediate parent cell (parent grid + parent-grid-local cell index). Cells are grouped by parent - // grid and merged into connected regions; one LGR is created per region, placed inside its parent - // grid (true LGR-in-LGR, any depth). - static void buildNestedLevel( RigEclipseCaseData* caseData, - const std::vector& cells, - const NestedHybridInput& input, - const std::map>& cellToParent, - int& nextGridId, - std::map& sourceCells ); + std::map& sourceCells ); static void updateActiveCellInfo( RigEclipseCaseData* caseData, const std::map& sourceCells ); diff --git a/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridResultTools.cpp b/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridResultTools.cpp index 26b587119e..7c3d7c8ae9 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridResultTools.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigNestedHybridGridResultTools.cpp @@ -75,7 +75,7 @@ QString RigNestedHybridGridResultTools::refineSidecarFilePath( const QString& gr //-------------------------------------------------------------------------------------------------- /// Nested hybrid grid: the parent mapping is provided in a sidecar GRDECL file named /// "_OLDIJK.grdecl" next to the grid file. It holds, per flat cell, the original -/// coarse cell IJK (OLDI/OLDJ/OLDK) and the local refined coordinates (TMPI/TMPJ/TMPK). +/// coarse cell IJK (OLDI/OLDJ/OLDK). /// Returns its path if it exists. //-------------------------------------------------------------------------------------------------- QString RigNestedHybridGridResultTools::oldIjkSidecarFilePath( const QString& gridFileName ) @@ -120,8 +120,8 @@ void RigNestedHybridGridResultTools::importRefineSidecarIfPresent( const QString } //-------------------------------------------------------------------------------------------------- -/// Load the OLDIJK sidecar (OLDI/OLDJ/OLDK/TMPI/TMPJ/TMPK) as input properties so the parent-cell -/// mapping is visible and scriptable, mirroring the REFINE property. +/// Load the OLDIJK sidecar as input properties so the parent-cell mapping is visible and scriptable, +/// mirroring the REFINE property. //-------------------------------------------------------------------------------------------------- void RigNestedHybridGridResultTools::importOldIjkSidecarIfPresent( const QString& gridFileName, RimEclipseInputPropertyCollection* inputPropertyCollection, @@ -179,9 +179,6 @@ void RigNestedHybridGridResultTools::reconstructNestedHybridGridIfPresent( const input.oldI = readIntKeyword( oldIjkContent, "OLDI" ); input.oldJ = readIntKeyword( oldIjkContent, "OLDJ" ); input.oldK = readIntKeyword( oldIjkContent, "OLDK" ); - input.tmpI = readIntKeyword( oldIjkContent, "TMPI" ); - input.tmpJ = readIntKeyword( oldIjkContent, "TMPJ" ); - input.tmpK = readIntKeyword( oldIjkContent, "TMPK" ); QString errorMessage; RigNestedHybridGridReconstructor::reconstruct( eclipseCaseData, input, &errorMessage ); @@ -354,7 +351,7 @@ RigEclipseResultAddress RigNestedHybridGridResultTools::computeCoarseAggregate( //-------------------------------------------------------------------------------------------------- /// Per refinement level, compute the aggregate (pore-volume-weighted average or sum) of a source -/// result over the cells of each immediate parent and broadcast it back onto that level's cells. All +/// result over the cells of each coarse parent and broadcast it back onto that level's cells. All /// other cells are left undefined (blank) so each level's result shows only that level. One result /// "_COARSE_L" is created per level present (stored on the active refined cells; /// the parent cells are inactive). @@ -409,7 +406,7 @@ std::vector RigNestedHybridGridResultTools::computePerL auto activeIndex = [&]( size_t reservoirCell ) { return activeCellInfo->cellResultIndex( ReservoirCellIndex( reservoirCell ) ).value(); }; // Collect every active reconstructed-LGR cell with its refinement level (from REFINE) and its - // immediate parent cell (from the LGR hierarchy). + // coarse parent cell (from the sibling LGR hierarchy). struct CellRef { size_t resultIndex; @@ -480,7 +477,7 @@ std::vector RigNestedHybridGridResultTools::computePerL { const std::vector& src = sourceTs[ts]; - // Accumulation keyed by (level, immediate parent) - cells of different levels are never + // Accumulation keyed by (level, coarse parent) - cells of different levels are never // accumulated together. The zero-bulk-volume filter excludes hidden duplicates in both modes. std::map>> acc; for ( const CellRef& cr : cellRefs ) diff --git a/ApplicationLibCode/UnitTests/RigNestedHybridGridReconstructor-Test.cpp b/ApplicationLibCode/UnitTests/RigNestedHybridGridReconstructor-Test.cpp index c18dd57232..36e2cf459a 100644 --- a/ApplicationLibCode/UnitTests/RigNestedHybridGridReconstructor-Test.cpp +++ b/ApplicationLibCode/UnitTests/RigNestedHybridGridReconstructor-Test.cpp @@ -56,7 +56,7 @@ // The test model TestModels/NestedHybridGrid/DROGON_NESTED.* is a single flat 150x84x96 EGRID where // the refined cells of each coarse (15x24x12) cell are appended in per-level I bands. The sidecars: // DROGON_NESTED_REFINE.grdecl : per-cell nesting level (1 base, 2/3/4 refined) -// DROGON_NESTED_OLDIJK.grdecl : OLDI/OLDJ/OLDK (coarse parent IJK) + TMPI/TMPJ/TMPK (local coords) +// DROGON_NESTED_OLDIJK.grdecl : OLDI/OLDJ/OLDK (coarse parent IJK) // // Opening the ~1.2M cell EGRID, parsing the sidecars and reconstructing the LGR hierarchy is // expensive (several seconds), so it is done once in SetUpTestSuite() and shared read-only by all @@ -98,9 +98,6 @@ RigNestedHybridGridReconstructor::NestedHybridInput buildInput( const QString& d 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; } @@ -241,19 +238,21 @@ TEST_F( RigNestedHybridGridReconstructorTest, ReconstructFromOldIjk ) } ASSERT_GT( level4Parents.size(), 0u ); - // Count the level-4 LGRs (named "LGR_NHG_L4_"). They are merged into connected regions, - // so there are far fewer than the number of coarse parents that own level-4 cells. + // Count the level-4 LGRs. Adjacent coarse parents with equal refinement dimensions are merged. size_t numLevel4Lgrs = 0; for ( size_t i = 1; i < mainGrid->gridCount(); i++ ) { RigLocalGrid* lgr = dynamic_cast( mainGrid->gridByIndex( i ) ); - if ( lgr && QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4_" ) ) numLevel4Lgrs++; + ASSERT_NE( lgr, nullptr ); + EXPECT_TRUE( lgr->isReconstructedGrid() ); + EXPECT_EQ( lgr->parentGrid(), mainGrid ); + if ( QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4" ) ) numLevel4Lgrs++; } EXPECT_GT( numLevel4Lgrs, 0u ); EXPECT_LT( numLevel4Lgrs, level4Parents.size() ); // merged, not one-per-parent - // Level 2 and level 3 each become one LGR refining the coarse grid; level 4 nests inside level 3 - // as merged regions. None is counted as an on-file grid. + // Every refinement level is reconstructed as sibling LGRs below the main grid. None is counted as + // an on-file grid. EXPECT_EQ( mainGrid->gridCount(), 1u + 2u + numLevel4Lgrs ); EXPECT_EQ( mainGrid->gridCountOnFile(), 1u ); @@ -262,7 +261,7 @@ TEST_F( RigNestedHybridGridReconstructorTest, ReconstructFromOldIjk ) ASSERT_TRUE( lgr2 != nullptr ); EXPECT_EQ( lgr2->cellCountI(), 26u ); EXPECT_EQ( lgr2->cellCountJ(), 44u ); - EXPECT_EQ( lgr2->cellCountK(), 48u ); + EXPECT_EQ( lgr2->cellCountK(), 24u ); EXPECT_TRUE( lgr2->isReconstructedGrid() ); EXPECT_EQ( lgr2->parentGrid(), mainGrid ); @@ -273,25 +272,23 @@ TEST_F( RigNestedHybridGridReconstructorTest, ReconstructFromOldIjk ) EXPECT_EQ( lgr3->cellCountJ(), 76u ); EXPECT_EQ( lgr3->cellCountK(), 48u ); - // Level 4 nests inside the level-3 LGR (true LGR-in-LGR). Find a level-4 LGR and verify its parent - // grid is the level-3 LGR, and the level-3 cell it subdivides points back to it. + // Level 4 is also reconstructed directly below the main grid. RigLocalGrid* lgr4 = nullptr; for ( size_t i = 1; i < mainGrid->gridCount(); i++ ) { RigLocalGrid* lgr = dynamic_cast( mainGrid->gridByIndex( i ) ); - if ( lgr && QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4_" ) ) + if ( lgr && QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4" ) ) { lgr4 = lgr; break; } } ASSERT_TRUE( lgr4 != nullptr ); - EXPECT_EQ( lgr4->parentGrid(), lgr3 ); + EXPECT_EQ( lgr4->parentGrid(), mainGrid ); EXPECT_TRUE( lgr4->isReconstructedGrid() ); { - // Pick a real (sourced) level-4 cell and verify its parent is a level-3 LGR cell whose subgrid - // is this level-4 LGR (parentCellIndex must be local to the level-3 grid). + // Pick a real (sourced) level-4 cell and verify that OLDIJK linked it to a coarse cell. const size_t l4Begin = lgr4->reservoirCellIndex( 0 ); const size_t l4End = l4Begin + lgr4->cellCount(); const std::map& srcAll = mainGrid->nestedHybridLgrSourceCells(); @@ -300,9 +297,9 @@ TEST_F( RigNestedHybridGridReconstructorTest, ReconstructFromOldIjk ) { if ( lgrGlobal < l4Begin || lgrGlobal >= l4End ) continue; size_t parentLocal = mainGrid->cell( lgrGlobal ).parentCellIndex(); - size_t parentGlobal = lgr3->reservoirCellIndex( parentLocal ); + size_t parentGlobal = mainGrid->reservoirCellIndex( parentLocal ); EXPECT_EQ( mainGrid->cell( parentGlobal ).subGrid(), lgr4 ); - EXPECT_EQ( mainGrid->cell( parentGlobal ).hostGrid(), lgr3 ); + EXPECT_EQ( mainGrid->cell( parentGlobal ).hostGrid(), mainGrid ); checkedNest = true; break; } @@ -510,7 +507,7 @@ TEST_F( RigNestedHybridGridReconstructorTest, CoarsePoreVolumeWeightedAggregate //-------------------------------------------------------------------------------------------------- /// QC per refinement level: _COARSE_L4 holds, on each level-4 cell, the pore-volume-weighted -/// average over the level-4 cells of its immediate (level-3) parent. +/// average over the level-4 cells of its coarse parent. //-------------------------------------------------------------------------------------------------- TEST_F( RigNestedHybridGridReconstructorTest, PerLevelPoreVolumeWeightedAggregate ) { @@ -547,7 +544,7 @@ TEST_F( RigNestedHybridGridReconstructorTest, PerLevelPoreVolumeWeightedAggregat for ( size_t i = 1; i < mainGrid->gridCount(); i++ ) { RigLocalGrid* lgr = dynamic_cast( mainGrid->gridByIndex( i ) ); - if ( !lgr || !QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4_" ) ) continue; + if ( !lgr || !QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4" ) ) continue; for ( size_t c = 0; c < lgr->cellCount(); c++ ) if ( ai->isActive( ReservoirCellIndex( lgr->reservoirCellIndex( c ) ) ) ) activeLevel4++; } @@ -557,12 +554,12 @@ TEST_F( RigNestedHybridGridReconstructorTest, PerLevelPoreVolumeWeightedAggregat EXPECT_GT( activeLevel4, 0u ); EXPECT_EQ( definedCount, activeLevel4 ); - // Find a level-4 LGR and group its cells by their (level-3) parent cell. + // Find a level-4 LGR and group its cells by their coarse parent cell. RigLocalGrid* lgr4 = nullptr; for ( size_t i = 1; i < mainGrid->gridCount(); i++ ) { RigLocalGrid* lgr = dynamic_cast( mainGrid->gridByIndex( i ) ); - if ( lgr && QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4_" ) ) + if ( lgr && QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4" ) ) { lgr4 = lgr; break; @@ -570,7 +567,7 @@ TEST_F( RigNestedHybridGridReconstructorTest, PerLevelPoreVolumeWeightedAggregat } ASSERT_TRUE( lgr4 != nullptr ); - std::map> byParent; // parent local cell -> level-4 global cells + std::map> byParent; // coarse parent cell -> level-4 global cells for ( size_t c = 0; c < lgr4->cellCount(); c++ ) { size_t global = lgr4->reservoirCellIndex( c ); @@ -668,7 +665,7 @@ TEST_F( RigNestedHybridGridReconstructorTest, CoarseSumAggregate ) //-------------------------------------------------------------------------------------------------- /// Per-level SUM aggregate: _COARSE_L4 holds, on each level-4 cell, the sum over the level-4 -/// cells of its immediate (level-3) parent, and is blank everywhere else. +/// cells of its coarse parent, and is blank everywhere else. //-------------------------------------------------------------------------------------------------- TEST_F( RigNestedHybridGridReconstructorTest, PerLevelSumAggregate ) { @@ -696,7 +693,7 @@ TEST_F( RigNestedHybridGridReconstructorTest, PerLevelSumAggregate ) for ( size_t i = 1; i < mainGrid->gridCount(); i++ ) { RigLocalGrid* lgr = dynamic_cast( mainGrid->gridByIndex( i ) ); - if ( !lgr || !QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4_" ) ) continue; + if ( !lgr || !QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4" ) ) continue; for ( size_t c = 0; c < lgr->cellCount(); c++ ) if ( ai->isActive( ReservoirCellIndex( lgr->reservoirCellIndex( c ) ) ) ) activeLevel4++; } @@ -706,12 +703,12 @@ TEST_F( RigNestedHybridGridReconstructorTest, PerLevelSumAggregate ) EXPECT_GT( activeLevel4, 0u ); EXPECT_EQ( definedCount, activeLevel4 ); - // Find a level-4 LGR and group its cells by their (level-3) parent cell. + // Find a level-4 LGR and group its cells by their coarse parent cell. RigLocalGrid* lgr4 = nullptr; for ( size_t i = 1; i < mainGrid->gridCount(); i++ ) { RigLocalGrid* lgr = dynamic_cast( mainGrid->gridByIndex( i ) ); - if ( lgr && QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4_" ) ) + if ( lgr && QString::fromStdString( lgr->gridName() ).startsWith( "LGR_NHG_L4" ) ) { lgr4 = lgr; break; @@ -719,7 +716,7 @@ TEST_F( RigNestedHybridGridReconstructorTest, PerLevelSumAggregate ) } ASSERT_TRUE( lgr4 != nullptr ); - std::map> byParent; // parent local cell -> level-4 global cells + std::map> byParent; // coarse parent cell -> level-4 global cells for ( size_t c = 0; c < lgr4->cellCount(); c++ ) { size_t global = lgr4->reservoirCellIndex( c );