mirror of
https://github.com/OPM/ResInsight.git
synced 2026-10-25 06:38:13 -05:00
Kristian Bendiksen
bcfa6e9fea
Well connections read from the restart file referenced the original refined cells in the main grid. These cells are hidden after reconstruction, but the well cell visualization does not check the invalid flag, so they were still drawn. Re-point well connections and well heads to the corresponding LGR cells and recompute the well cells per grid.
561 lines
26 KiB
C++
561 lines
26 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// 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 "RigNestedHybridGridReconstructor.h"
|
|
|
|
#include "RiaLogging.h"
|
|
|
|
#include "RigActiveCellInfo.h"
|
|
#include "RigCaseCellResultsData.h"
|
|
#include "RigCell.h"
|
|
#include "RigEclipseCaseData.h"
|
|
#include "RigEclipseResultAddress.h"
|
|
#include "RigGridBase.h"
|
|
#include "RigLocalGrid.h"
|
|
#include "RigMainGrid.h"
|
|
#include "RigNNCData.h"
|
|
#include "RigNestedHybridGridLayout.h"
|
|
#include "RigNestedHybridGridResultTools.h"
|
|
#include "RigNncConnection.h"
|
|
#include "Well/RigSimWellData.h"
|
|
#include "Well/RigWellResultBranch.h"
|
|
#include "Well/RigWellResultFrame.h"
|
|
#include "Well/RigWellResultPoint.h"
|
|
|
|
#include "RiaDefines.h"
|
|
|
|
#include <QRegularExpression>
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <climits>
|
|
#include <cmath>
|
|
#include <map>
|
|
#include <set>
|
|
|
|
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;
|
|
}
|
|
} // namespace
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigNestedHybridGridReconstructor::reconstruct( RigEclipseCaseData* caseData, const NestedHybridInput& input, QString* errorMessage )
|
|
{
|
|
auto setError = [&]( const QString& msg )
|
|
{
|
|
if ( errorMessage ) *errorMessage = msg;
|
|
RiaLogging::warning( ( "Nested hybrid grid: " + msg ).toStdString() );
|
|
return false;
|
|
};
|
|
|
|
if ( !caseData || !caseData->mainGrid() ) return setError( "No grid data." );
|
|
|
|
RigMainGrid* grid = caseData->mainGrid();
|
|
const size_t cellCount = grid->cellCount();
|
|
|
|
if ( input.refine.size() != cellCount )
|
|
{
|
|
return setError( QString( "REFINE size %1 does not match main grid cell count %2." ).arg( input.refine.size() ).arg( cellCount ) );
|
|
}
|
|
if ( input.oldI.size() != cellCount || input.oldJ.size() != cellCount || input.oldK.size() != cellCount )
|
|
{
|
|
return setError( "OLDI/OLDJ/OLDK size does not match main grid cell count." );
|
|
}
|
|
if ( grid->gridCount() > 1 )
|
|
{
|
|
return setError( "Grid already contains local grids; skipping reconstruction." );
|
|
}
|
|
|
|
const size_t nx = grid->cellCountI();
|
|
const size_t ny = grid->cellCountJ();
|
|
const size_t nz = grid->cellCountK();
|
|
|
|
// Coarse grid dimensions are the maximum 1-based OLD indices. The flat grid embeds each coarse
|
|
// cell at (OLDI-1, OLDJ-1, (OLDK-1)*KF), where KF is the global K refinement factor.
|
|
int coarseNx = 0, coarseNy = 0, coarseNz = 0;
|
|
for ( size_t f = 0; f < cellCount; f++ )
|
|
{
|
|
coarseNx = std::max( coarseNx, input.oldI[f] );
|
|
coarseNy = std::max( coarseNy, input.oldJ[f] );
|
|
coarseNz = std::max( coarseNz, input.oldK[f] );
|
|
}
|
|
const size_t kFactor = ( coarseNz > 0 && nz % (size_t)coarseNz == 0 ) ? nz / (size_t)coarseNz : 1;
|
|
|
|
const RigNestedHybridGridLayout layout( cvf::Vec3st( (size_t)coarseNx, (size_t)coarseNy, (size_t)coarseNz ), kFactor );
|
|
|
|
// Group all refined cells (REFINE > 1) by level, and record each refined cell's parent COARSE
|
|
// cell natural index (for the volume-weighted QC aggregate).
|
|
std::map<int, std::vector<size_t>> cellsByLevel;
|
|
std::map<size_t, size_t> coarseParents; // flat cell -> coarse natural index
|
|
for ( size_t f = 0; f < cellCount; f++ )
|
|
{
|
|
const int level = input.refine[f];
|
|
if ( level <= 1 ) continue;
|
|
const int oi = input.oldI[f];
|
|
const int oj = input.oldJ[f];
|
|
const int ok = input.oldK[f];
|
|
if ( oi < 1 || oj < 1 || ok < 1 ) continue; // padding / unmapped
|
|
|
|
cellsByLevel[level].push_back( f );
|
|
coarseParents[f] = naturalIndex( (size_t)( oi - 1 ), (size_t)( oj - 1 ), (size_t)( ok - 1 ), (size_t)coarseNx, (size_t)coarseNy );
|
|
}
|
|
|
|
if ( cellsByLevel.empty() ) return setError( "No refined cells (REFINE all <= 1)." );
|
|
|
|
const size_t origCellCount = grid->totalCellCount(); // before any LGR cells are appended
|
|
std::map<size_t, size_t> sourceCells; // LGR global cell index -> source flat cell index
|
|
int nextGridId = (int)grid->gridCount(); // main grid is 0
|
|
|
|
using CoarseIjk = std::array<int, 3>;
|
|
using Factor = std::array<size_t, 3>;
|
|
|
|
// 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 )
|
|
{
|
|
std::map<CoarseIjk, std::vector<size_t>> cellsByParent;
|
|
for ( size_t f : cells )
|
|
cellsByParent[{ input.oldI[f], input.oldJ[f], input.oldK[f] }].push_back( f );
|
|
|
|
std::map<CoarseIjk, std::array<std::vector<size_t>, 3>> coordinatesByParent;
|
|
std::map<Factor, std::set<CoarseIjk>> parentsByFactor;
|
|
for ( const auto& [parent, parentCells] : cellsByParent )
|
|
{
|
|
auto& coordinates = coordinatesByParent[parent];
|
|
for ( size_t f : parentCells )
|
|
{
|
|
const std::array<size_t, 3> flatIjk = { f % nx, ( f / nx ) % ny, f / ( nx * ny ) };
|
|
for ( int axis = 0; axis < 3; axis++ )
|
|
coordinates[axis].push_back( flatIjk[axis] );
|
|
}
|
|
|
|
Factor factor;
|
|
for ( int axis = 0; axis < 3; axis++ )
|
|
{
|
|
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 );
|
|
}
|
|
|
|
// Collect the connected components first so the component count is known before naming.
|
|
std::vector<std::pair<Factor, std::vector<CoarseIjk>>> components;
|
|
for ( const auto& [factor, parents] : parentsByFactor )
|
|
{
|
|
std::set<CoarseIjk> remaining = parents;
|
|
while ( !remaining.empty() )
|
|
{
|
|
std::vector<CoarseIjk> component;
|
|
std::vector<CoarseIjk> stack = { *remaining.begin() };
|
|
remaining.erase( stack.front() );
|
|
|
|
while ( !stack.empty() )
|
|
{
|
|
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 );
|
|
}
|
|
}
|
|
|
|
components.push_back( { factor, std::move( component ) } );
|
|
}
|
|
}
|
|
|
|
for ( size_t componentIndex = 0; componentIndex < components.size(); componentIndex++ )
|
|
{
|
|
const auto& [factor, component] = components[componentIndex];
|
|
|
|
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] );
|
|
}
|
|
|
|
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<size_t> 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];
|
|
caf::VecIjk0 host = layout.flatIjkFromOriginalIjk( caf::VecIjk0( oi - 1, oj - 1, ok - 1 ) );
|
|
boxToParent[local] = naturalIndex( host.x(), host.y(), host.z(), nx, ny );
|
|
}
|
|
|
|
std::vector<size_t> boxToFlat( dims.x() * dims.y() * dims.z(), cvf::UNDEFINED_SIZE_T );
|
|
for ( const CoarseIjk& parent : component )
|
|
{
|
|
const auto& coordinates = coordinatesByParent[parent];
|
|
for ( size_t f : cellsByParent[parent] )
|
|
{
|
|
const std::array<size_t, 3> 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;
|
|
}
|
|
}
|
|
|
|
// Single component on a level keeps the bare level name; several get 1-based suffixes.
|
|
QString gridName = QString( "LGR_NHG_L%1" ).arg( level );
|
|
if ( components.size() > 1 ) gridName += QString( "_%1" ).arg( componentIndex + 1 );
|
|
buildLocalGrid( caseData, grid, nextGridId++, gridName, dims, boxToFlat, boxToParent, sourceCells );
|
|
}
|
|
}
|
|
|
|
if ( sourceCells.empty() ) return setError( "No refined regions could be reconstructed." );
|
|
|
|
grid->setNestedHybridLgrSourceCells( sourceCells );
|
|
grid->setNestedHybridCoarseParents( coarseParents );
|
|
grid->setNestedHybridLayout( layout );
|
|
|
|
updateActiveCellInfo( caseData, sourceCells );
|
|
extendFullLengthResults( caseData, origCellCount, sourceCells );
|
|
repointNncConnections( caseData, origCellCount, sourceCells );
|
|
repointWellResultPoints( caseData, origCellCount, sourceCells );
|
|
|
|
// Fill LGR cells of any already-loaded active-cell-indexed results from their source flat cells.
|
|
for ( auto model : { RiaDefines::PorosityModelType::MATRIX_MODEL, RiaDefines::PorosityModelType::FRACTURE_MODEL } )
|
|
{
|
|
if ( RigCaseCellResultsData* results = caseData->results( model ) ) RigNestedHybridGridResultTools::extendLgrResults( results );
|
|
}
|
|
|
|
// 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)" ).arg( grid->gridCount() - 1 ).arg( sourceCells.size() ).toStdString() );
|
|
|
|
return true;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
QString RigNestedHybridGridReconstructor::levelGroupName( const QString& gridName )
|
|
{
|
|
static const QRegularExpression re( "^(LGR_NHG_L\\d+)(_\\d+)?$" );
|
|
const auto match = re.match( gridName );
|
|
return match.hasMatch() ? match.captured( 1 ) : gridName;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Build one RigLocalGrid by appending contiguous cells/nodes to the global arrays, copying the real
|
|
/// geometry from the scattered flat refined cells, hiding the originals, and linking every LGR cell
|
|
/// to the parent (coarse) cell it subdivides. boxToFlat maps each LGR local cell to its source flat
|
|
/// cell (UNDEFINED for holes); boxToParent maps it to its parent cell. The LGR may span many parent
|
|
/// cells; each parent owning a real cell gets its subGrid set. Mirrors RicCreateTemporaryLgrFeature.
|
|
//--------------------------------------------------------------------------------------------------
|
|
RigLocalGrid* RigNestedHybridGridReconstructor::buildLocalGrid( RigEclipseCaseData* caseData,
|
|
RigGridBase* parentGrid,
|
|
int gridId,
|
|
const QString& gridName,
|
|
const cvf::Vec3st& dims,
|
|
const std::vector<size_t>& boxToFlat,
|
|
const std::vector<size_t>& boxToParent,
|
|
std::map<size_t, size_t>& sourceCells )
|
|
{
|
|
RigMainGrid* grid = caseData->mainGrid();
|
|
|
|
const size_t lgrCellCount = dims.x() * dims.y() * dims.z();
|
|
if ( lgrCellCount == 0 ) return nullptr;
|
|
|
|
RigLocalGrid* localGrid = new RigLocalGrid( grid );
|
|
localGrid->setGridId( gridId );
|
|
localGrid->setGridName( gridName.toStdString() );
|
|
localGrid->setCellCounts( dims );
|
|
localGrid->setAsReconstructedGrid( true ); // synthesized in memory, not present in the result file
|
|
|
|
const size_t cellStart = grid->totalCellCount();
|
|
const size_t nodeStart = grid->nodes().size();
|
|
localGrid->setIndexToStartOfCells( cellStart );
|
|
grid->addLocalGrid( localGrid );
|
|
|
|
{
|
|
RigCell defaultCell;
|
|
defaultCell.setHostGrid( localGrid );
|
|
grid->reservoirCells().resize( cellStart + lgrCellCount, defaultCell );
|
|
grid->nodes().resize( nodeStart + lgrCellCount * 8, cvf::Vec3d( 0, 0, 0 ) );
|
|
}
|
|
|
|
for ( size_t glc = 0; glc < lgrCellCount; glc++ )
|
|
{
|
|
RigCell& lgrCell = grid->cell( cellStart + glc );
|
|
lgrCell.setHostGrid( localGrid );
|
|
lgrCell.setGridLocalCellIndex( glc );
|
|
|
|
// boxToParent holds the parent cell index LOCAL to parentGrid; setParentCellIndex expects the
|
|
// parent-grid-local index, while subGrid is set on the parent's global reservoir cell.
|
|
const size_t parentLocal = boxToParent[glc];
|
|
const size_t parentGlobal = parentLocal != cvf::UNDEFINED_SIZE_T ? parentGrid->reservoirCellIndex( parentLocal )
|
|
: cvf::UNDEFINED_SIZE_T;
|
|
lgrCell.setParentCellIndex( parentLocal );
|
|
|
|
const size_t flat = boxToFlat[glc];
|
|
if ( flat == cvf::UNDEFINED_SIZE_T )
|
|
{
|
|
lgrCell.setInvalid( true ); // hole - keep the regular IJK box but no geometry
|
|
continue;
|
|
}
|
|
|
|
RigCell& flatCell = grid->cell( flat );
|
|
|
|
// Copy the real corner geometry from the scattered flat refined cell.
|
|
for ( size_t c = 0; c < 8; c++ )
|
|
{
|
|
size_t newNodeIdx = nodeStart + glc * 8 + c;
|
|
grid->nodes()[newNodeIdx] = grid->nodes()[flatCell.cornerIndices()[c]];
|
|
lgrCell.cornerIndices()[c] = newNodeIdx;
|
|
}
|
|
lgrCell.setInvalid( flatCell.isInvalid() );
|
|
|
|
// Hide the original scattered cell - the LGR cell now represents it.
|
|
flatCell.setInvalid( true );
|
|
|
|
// This parent cell is really refined by this LGR.
|
|
if ( parentGlobal != cvf::UNDEFINED_SIZE_T )
|
|
{
|
|
grid->cell( parentGlobal ).setSubGrid( localGrid );
|
|
}
|
|
|
|
sourceCells[cellStart + glc] = flat;
|
|
}
|
|
|
|
localGrid->setParentGrid( parentGrid );
|
|
|
|
return localGrid;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// 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
|
|
/// invariant "result index == position in the active-cell list" that result calculators rely on.
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigNestedHybridGridReconstructor::updateActiveCellInfo( RigEclipseCaseData* caseData, const std::map<size_t, size_t>& sourceCells )
|
|
{
|
|
RigMainGrid* grid = caseData->mainGrid();
|
|
|
|
const std::array<RiaDefines::PorosityModelType, 2> models = { RiaDefines::PorosityModelType::MATRIX_MODEL,
|
|
RiaDefines::PorosityModelType::FRACTURE_MODEL };
|
|
|
|
for ( auto model : models )
|
|
{
|
|
RigActiveCellInfo* actInfo = caseData->activeCellInfo( model );
|
|
if ( !actInfo ) continue;
|
|
|
|
// Skip porosity models that are not in use (e.g. the fracture model in a single-porosity
|
|
// case). Adding cells to an unused model makes its active count non-zero, which makes the
|
|
// result reader treat the case as dual-porosity and mis-split every dynamic result.
|
|
const size_t oldActiveCellCount = actInfo->reservoirActiveCellCount();
|
|
if ( oldActiveCellCount == 0 ) continue;
|
|
|
|
actInfo->setReservoirCellCount( grid->totalCellCount() );
|
|
actInfo->setGridCount( grid->gridCount() );
|
|
|
|
size_t nextActiveIndex = oldActiveCellCount;
|
|
for ( size_t gridIdx = 1; gridIdx < grid->gridCount(); gridIdx++ )
|
|
{
|
|
RigGridBase* lgr = grid->gridByIndex( gridIdx );
|
|
size_t lgrActiveCount = 0;
|
|
for ( size_t local = 0; local < lgr->cellCount(); local++ )
|
|
{
|
|
size_t global = lgr->reservoirCellIndex( local );
|
|
auto it = sourceCells.find( global );
|
|
if ( it == sourceCells.end() ) continue;
|
|
if ( !actInfo->isActive( ReservoirCellIndex( it->second ) ) ) continue;
|
|
|
|
actInfo->setCellResultIndex( ReservoirCellIndex( global ), ActiveCellIndex( nextActiveIndex ) );
|
|
nextActiveIndex++;
|
|
lgrActiveCount++;
|
|
}
|
|
actInfo->setGridActiveCellCounts( gridIdx, lgrActiveCount );
|
|
}
|
|
|
|
actInfo->computeDerivedData();
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Full-length (all-cells) result arrays are indexed by global reservoir cell index, so the new LGR
|
|
/// cells (appended at the end) need their own entries. Copy each LGR cell's value from the source
|
|
/// flat cell it was built from. Active-cell-indexed results are addressed via the result index that
|
|
/// was already transferred in updateActiveCellInfo, so they are left untouched.
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigNestedHybridGridReconstructor::extendFullLengthResults( RigEclipseCaseData* caseData,
|
|
size_t origCellCount,
|
|
const std::map<size_t, size_t>& sourceCells )
|
|
{
|
|
const size_t newTotal = caseData->mainGrid()->totalCellCount();
|
|
|
|
const std::array<RiaDefines::PorosityModelType, 2> models = { RiaDefines::PorosityModelType::MATRIX_MODEL,
|
|
RiaDefines::PorosityModelType::FRACTURE_MODEL };
|
|
|
|
for ( auto model : models )
|
|
{
|
|
RigCaseCellResultsData* results = caseData->results( model );
|
|
if ( !results ) continue;
|
|
|
|
for ( const RigEclipseResultAddress& addr : results->existingResults() )
|
|
{
|
|
std::vector<std::vector<double>>* timesteps = results->modifiableCellScalarResultTimesteps( addr );
|
|
if ( !timesteps ) continue;
|
|
|
|
for ( std::vector<double>& values : *timesteps )
|
|
{
|
|
// Only full-length (all-cells) arrays need extending: their length equals the
|
|
// pre-reconstruction cell count. Active-cell-indexed arrays have the (smaller)
|
|
// active-cell length and are handled via the transferred result index.
|
|
if ( values.size() != origCellCount ) continue;
|
|
|
|
values.resize( newTotal, HUGE_VAL );
|
|
for ( const auto& [global, flat] : sourceCells )
|
|
{
|
|
if ( flat < origCellCount ) values[global] = values[flat];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// The file NNCs connect coarse cells to the original scattered refined cells. After those cells
|
|
/// were moved into the LGR and hidden, re-point each affected connection to the corresponding LGR
|
|
/// cell. Connection result data (transmissibilities) and polygons are indexed/located by position
|
|
/// and geometry, both unchanged, so only the cell indices need updating.
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigNestedHybridGridReconstructor::repointNncConnections( RigEclipseCaseData* caseData,
|
|
size_t origCellCount,
|
|
const std::map<size_t, size_t>& sourceCells )
|
|
{
|
|
RigMainGrid* grid = caseData->mainGrid();
|
|
if ( !grid ) return;
|
|
|
|
RigNNCData* nncData = grid->nncData();
|
|
if ( !nncData ) return;
|
|
|
|
// Inverse map: original flat cell index -> new LGR global cell index.
|
|
std::vector<size_t> flatToLgr( origCellCount, cvf::UNDEFINED_SIZE_T );
|
|
for ( const auto& [global, flat] : sourceCells )
|
|
{
|
|
if ( flat < origCellCount ) flatToLgr[flat] = global;
|
|
}
|
|
|
|
auto mapped = [&]( size_t idx ) { return ( idx < origCellCount && flatToLgr[idx] != cvf::UNDEFINED_SIZE_T ) ? flatToLgr[idx] : idx; };
|
|
|
|
std::vector<RigConnection>& connections = nncData->allConnections();
|
|
for ( size_t i = 0; i < connections.size(); i++ )
|
|
{
|
|
RigConnection& conn = connections[i];
|
|
size_t a = conn.c1GlobIdx();
|
|
size_t b = conn.c2GlobIdx();
|
|
size_t na = mapped( a );
|
|
size_t nb = mapped( b );
|
|
if ( na != a || nb != b )
|
|
{
|
|
connections[i] = RigConnection( na, nb, conn.face(), conn.polygon() );
|
|
}
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// The well connections read from the restart file reference the original scattered refined cells
|
|
/// in the main grid. Those cells are hidden after reconstruction, but the well cell geometry does not
|
|
/// check the invalid flag, so the stale connections would show the hidden flat cells. Re-point each
|
|
/// such connection and well head to the corresponding LGR cell.
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigNestedHybridGridReconstructor::repointWellResultPoints( RigEclipseCaseData* caseData,
|
|
size_t origCellCount,
|
|
const std::map<size_t, size_t>& sourceCells )
|
|
{
|
|
RigMainGrid* grid = caseData->mainGrid();
|
|
if ( !grid || caseData->wellResults().empty() ) return;
|
|
|
|
std::vector<size_t> flatToLgr( origCellCount, cvf::UNDEFINED_SIZE_T );
|
|
for ( const auto& [global, flat] : sourceCells )
|
|
{
|
|
if ( flat < origCellCount ) flatToLgr[flat] = global;
|
|
}
|
|
|
|
auto repoint = [&]( RigWellResultPoint& point )
|
|
{
|
|
// The flat grid has no LGRs of its own, so all original connections are in the main grid
|
|
if ( !point.isCell() || point.gridIndex() != 0 || point.cellIndex() >= origCellCount ) return;
|
|
|
|
const size_t global = flatToLgr[point.cellIndex()];
|
|
if ( global == cvf::UNDEFINED_SIZE_T ) return;
|
|
|
|
const RigCell& lgrCell = grid->cell( global );
|
|
point.setGridIndex( lgrCell.hostGrid()->gridIndex() );
|
|
point.setGridCellIndex( lgrCell.gridLocalCellIndex() );
|
|
};
|
|
|
|
cvf::Collection<RigSimWellData> wells = caseData->wellResults();
|
|
for ( auto& well : wells )
|
|
{
|
|
for ( RigWellResultFrame& frame : well->m_wellCellsTimeSteps )
|
|
{
|
|
RigWellResultPoint wellHead = frame.wellHead();
|
|
repoint( wellHead );
|
|
frame.setWellHead( wellHead );
|
|
|
|
std::vector<RigWellResultBranch> branches = frame.wellResultBranches();
|
|
for ( RigWellResultBranch& branch : branches )
|
|
{
|
|
std::vector<RigWellResultPoint> points = branch.branchResultPoints();
|
|
for ( RigWellResultPoint& point : points )
|
|
repoint( point );
|
|
branch.setBranchResultPoints( points );
|
|
}
|
|
frame.setWellResultBranches( branches );
|
|
}
|
|
}
|
|
|
|
// Recomputes the well cells per grid for the new grid set
|
|
caseData->setSimWellData( wells );
|
|
}
|