Skip to content

Commit

Permalink
Lock lgr radial grid to outer pillars
Browse files Browse the repository at this point in the history
  • Loading branch information
@jonjenssen
jonjenssen committed May 24, 2024
1 parent ebd92c7 commit 9d30a75
Show file tree
Hide file tree
Showing 3 changed files with 110 additions and 61 deletions.
134 changes: 78 additions & 56 deletions ApplicationLibCode/FileInterface/RifOpmRadialGridTools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,72 +188,94 @@ void RifOpmRadialGridTools::transferCoordinatesRadial( Opm::EclIO::EGrid& opmMai
// First grid dimension is radius, check if cell has are at the outer-most slice
if ( !hostCellGlobalIndices.empty() && ( gridDimension[0] - 1 == ijkCell[0] ) )
{
std::array<double, 8> cellRadius{};
std::array<double, 8> cellTheta{};
std::array<double, 8> cellZ{};
opmGrid.getRadialCellCorners( ijkCell, cellRadius, cellTheta, cellZ );
auto hostCellIndex = hostCellGlobalIndices[opmCellIndex];

double maxRadius = *std::max_element( cellRadius.begin(), cellRadius.end() );
lockToHostPillars( riNode, opmMainGrid, opmGrid, ijkCell, hostCellIndex, opmCellIndex, opmNodeIndex, xCenterCoordOpm, yCenterCoordOpm );
}
}
}
}

// Check if the radius is at the outer surface of the radial grid
// Adjust the outer nodes to match the corner pillars of the host cell
const double epsilon = 0.15;
if ( fabs( maxRadius - cellRadius[opmNodeIndex] ) < epsilon * cellRadius[opmNodeIndex] )
{
const auto hostCellIndex = hostCellGlobalIndices[opmCellIndex];
//--------------------------------------------------------------------------------------------------
//
//--------------------------------------------------------------------------------------------------
void RifOpmRadialGridTools::lockToHostPillars( cvf::Vec3d& riNode,
Opm::EclIO::EGrid& opmMainGrid,
Opm::EclIO::EGrid& opmGrid,
std::array<int, 3>& ijkCell,
int hostCellIndex,
int opmCellIndex,
size_t opmNodeIndex,
double xCenterCoordOpm,
double yCenterCoordOpm )
{
std::array<double, 8> cellRadius{};
std::array<double, 8> cellTheta{};
std::array<double, 8> cellZ{};
opmGrid.getRadialCellCorners( ijkCell, cellRadius, cellTheta, cellZ );

double closestPillarDistance = std::numeric_limits<double>::max();
int closestPillarIndex = -1;
double maxRadius = *std::max_element( cellRadius.begin(), cellRadius.end() );

const auto cylinderCoordX = opmX[opmNodeIndex] + xCenterCoordOpm;
const auto cylinderCoordY = opmY[opmNodeIndex] + yCenterCoordOpm;
const auto cylinderCoordZ = opmZ[opmNodeIndex];
// Check if the radius is at the outer surface of the radial grid
// Adjust the outer nodes to match the corner pillars of the host cell
const double epsilon = 0.15;
if ( fabs( maxRadius - cellRadius[opmNodeIndex] ) < epsilon * cellRadius[opmNodeIndex] )
{
std::array<double, 8> opmX{};
std::array<double, 8> opmY{};
std::array<double, 8> opmZ{};

const cvf::Vec3d coordinateOnCylinder = cvf::Vec3d( cylinderCoordX, cylinderCoordY, cylinderCoordZ );
opmGrid.getCellCorners( opmCellIndex, opmX, opmY, opmZ );

const auto candidates = computeSnapToCoordinates( opmMainGrid, opmGrid, hostCellIndex, opmCellIndex );
for ( int pillarIndex = 0; pillarIndex < static_cast<int>( candidates.size() ); pillarIndex++ )
{
for ( const auto& c : candidates[pillarIndex] )
{
double distance = coordinateOnCylinder.pointDistance( c );
if ( distance < closestPillarDistance )
{
closestPillarDistance = distance;
closestPillarIndex = pillarIndex;
}
}
}
double closestPillarDistance = std::numeric_limits<double>::max();
int closestPillarIndex = -1;

if ( closestPillarDistance < std::numeric_limits<double>::max() )
{
const auto& pillarCordinates = candidates[closestPillarIndex];

int layerCount = static_cast<int>( pillarCordinates.size() / 2 );
int layerIndexInMainGridCell = ijkCell[2] % layerCount;
int localNodeIndex = opmNodeIndex % 8;

cvf::Vec3d closestPillarCoord;
if ( localNodeIndex < 4 )
{
// Top of cell
int pillarCoordIndex = layerIndexInMainGridCell * 2;
closestPillarCoord = pillarCordinates[pillarCoordIndex];
}
else
{
// Bottom of cell
int pillarCoordIndex = layerIndexInMainGridCell * 2 + 1;
closestPillarCoord = pillarCordinates[pillarCoordIndex];
}

riNode.x() = closestPillarCoord.x();
riNode.y() = closestPillarCoord.y();
riNode.z() = -closestPillarCoord.z();
}
const auto cylinderCoordX = opmX[opmNodeIndex] + xCenterCoordOpm;
const auto cylinderCoordY = opmY[opmNodeIndex] + yCenterCoordOpm;
const auto cylinderCoordZ = opmZ[opmNodeIndex];

const cvf::Vec3d coordinateOnCylinder = cvf::Vec3d( cylinderCoordX, cylinderCoordY, cylinderCoordZ );

const auto candidates = computeSnapToCoordinates( opmMainGrid, opmGrid, hostCellIndex, opmCellIndex );
for ( int pillarIndex = 0; pillarIndex < static_cast<int>( candidates.size() ); pillarIndex++ )
{
for ( const auto& c : candidates[pillarIndex] )
{
double distance = coordinateOnCylinder.pointDistance( c );
if ( distance < closestPillarDistance )
{
closestPillarDistance = distance;
closestPillarIndex = pillarIndex;
}
}
}

if ( closestPillarDistance < std::numeric_limits<double>::max() )
{
const auto& pillarCordinates = candidates[closestPillarIndex];

int layerCount = static_cast<int>( pillarCordinates.size() / 2 );
int layerIndexInMainGridCell = ijkCell[2] % layerCount;
int localNodeIndex = opmNodeIndex % 8;

cvf::Vec3d closestPillarCoord;
if ( localNodeIndex < 4 )
{
// Top of cell
int pillarCoordIndex = layerIndexInMainGridCell * 2;
closestPillarCoord = pillarCordinates[pillarCoordIndex];
}
else
{
// Bottom of cell
int pillarCoordIndex = layerIndexInMainGridCell * 2 + 1;
closestPillarCoord = pillarCordinates[pillarCoordIndex];
}

riNode.x() = closestPillarCoord.x();
riNode.y() = closestPillarCoord.y();
riNode.z() = -closestPillarCoord.z();
}
}
}

Expand Down
10 changes: 10 additions & 0 deletions ApplicationLibCode/FileInterface/RifOpmRadialGridTools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,4 +53,14 @@ class RifOpmRadialGridTools

static std::vector<std::vector<cvf::Vec3d>>
computeSnapToCoordinates( Opm::EclIO::EGrid& opmMainGrid, Opm::EclIO::EGrid& opmGrid, int mainGridCellIndex, int lgrCellIndex );

static void lockToHostPillars( cvf::Vec3d& riNode,
Opm::EclIO::EGrid& opmMainGrid,
Opm::EclIO::EGrid& opmGrid,
std::array<int, 3>& ijkCell,
int hostCellIndex,
int opmCellIndex,
size_t opmNodeIndex,
double xCenterCoordOpm,
double yCenterCoordOpm );
};
27 changes: 22 additions & 5 deletions ApplicationLibCode/FileInterface/RifReaderOpmCommon.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,15 +262,16 @@ void RifReaderOpmCommon::transferGeometry( Opm::EclIO::EGrid& opmMainGrid,

const bool isRadialGrid = opmGrid.is_radial();

const auto& gridDimension = opmGrid.dimension();
const auto& hostCellGlobalIndices = opmGrid.hostCellsGlobalIndex();

// Compute the center of the LGR radial grid cells for each K layer
std::map<int, std::pair<double, double>> radialGridCenterTopLayerOpm =
RifOpmRadialGridTools::computeXyCenterForTopOfCells( opmMainGrid, opmGrid, localGrid );

// same mapping as libecl
const size_t cellMappingECLRi[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };

const auto host_parentcell = opmGrid.hostCellsGlobalIndex();

#pragma omp parallel for
for ( int opmCellIndex = 0; opmCellIndex < static_cast<int>( localGrid->cellCount() ); opmCellIndex++ )
{
Expand Down Expand Up @@ -304,9 +305,9 @@ void RifReaderOpmCommon::transferGeometry( Opm::EclIO::EGrid& opmMainGrid,
}

// parent cell index
if ( ( host_parentcell.size() > (size_t)opmCellIndex ) && host_parentcell[opmCellIndex] >= 0 )
if ( ( hostCellGlobalIndices.size() > (size_t)opmCellIndex ) && hostCellGlobalIndices[opmCellIndex] >= 0 )
{
cell.setParentCellIndex( host_parentcell[opmCellIndex] );
cell.setParentCellIndex( hostCellGlobalIndices[opmCellIndex] );
}
else
{
Expand All @@ -333,6 +334,22 @@ void RifReaderOpmCommon::transferGeometry( Opm::EclIO::EGrid& opmMainGrid,
riNode.z() = -opmZ[opmNodeIndex];

cell.cornerIndices()[riCornerIndex] = riNodeIndex;

// First grid dimension is radius, check if cell has are at the outer-most slice
if ( isRadialGrid && !hostCellGlobalIndices.empty() && ( gridDimension[0] - 1 == opmIJK[0] ) )
{
auto hostCellIndex = hostCellGlobalIndices[opmCellIndex];

RifOpmRadialGridTools::lockToHostPillars( riNode,
opmMainGrid,
opmGrid,
opmIJK,
hostCellIndex,
opmCellIndex,
opmNodeIndex,
xCenterCoordOpm,
yCenterCoordOpm );
}
}

cell.setInvalid( cell.isLongPyramidCell() );
Expand All @@ -344,7 +361,7 @@ void RifReaderOpmCommon::transferGeometry( Opm::EclIO::EGrid& opmMainGrid,

if ( parentGrid != nullptr )
{
for ( auto localCellInGlobalIdx : opmGrid.hostCellsGlobalIndex() )
for ( auto localCellInGlobalIdx : hostCellGlobalIndices )
{
parentGrid->cell( localCellInGlobalIdx ).setSubGrid( realLocalGrid );
}
Expand Down

0 comments on commit 9d30a75

Please sign in to comment.