mirror of
https://github.com/OPM/ResInsight.git
synced 2026-11-13 09:37:29 -06:00
Magne Sjaastad
11d8b926bb
Encapsulate the node table of RigMainGrid behind nodes() and modifiableNodes() returning spans, allocateNodes() and clearNodes(). Readers and LGR builders fill the node table in place through allocateNodes() instead of building a temporary vector, which reduces peak import memory by about 1.6 GiB for a 9.8M cell grid. The radial grid conversion and map axes transform update nodes in place.
420 lines
18 KiB
C++
420 lines
18 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (C) 2025- Equinor ASA
|
|
//
|
|
// ResInsight is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License, or
|
|
// (at your option) any later version.
|
|
//
|
|
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
// FITNESS FOR A PARTICULAR PURPOSE.
|
|
//
|
|
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
|
// for more details.
|
|
//
|
|
/////////////////////////////////////////////////////////////////////////////////
|
|
|
|
#include "RimCornerPointCase.h"
|
|
|
|
#include "RiaApplication.h"
|
|
#include "RiaDefines.h"
|
|
#include "RiaLogging.h"
|
|
#include "RiaPreferencesGrid.h"
|
|
#include "RiaQStringFormatter.h"
|
|
|
|
#include "RifInputPropertyLoader.h"
|
|
#include "RifRoffFileTools.h"
|
|
|
|
#include "RigActiveCellInfo.h"
|
|
#include "RigCaseCellResultsData.h"
|
|
#include "RigEclipseCaseData.h"
|
|
#include "RigEclipseResultAddress.h"
|
|
#include "RigMainGrid.h"
|
|
#include "RimReloadCaseTools.h"
|
|
|
|
#include "RimEclipseInputProperty.h"
|
|
#include "RimReservoirCellResultsStorage.h"
|
|
|
|
#include "cafPdmObjectScriptingCapability.h"
|
|
#include "cafProgressInfo.h"
|
|
|
|
#include <QDir>
|
|
#include <QFileInfo>
|
|
|
|
#include <algorithm>
|
|
#include <chrono>
|
|
|
|
#ifdef USE_OPENMP
|
|
#include <omp.h>
|
|
#endif
|
|
|
|
using namespace std::chrono;
|
|
|
|
CAF_PDM_SOURCE_INIT( RimCornerPointCase, "CornerPointCase", "RimCornerPointCase" );
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RimCornerPointCase::RimCornerPointCase()
|
|
: RimEclipseCase()
|
|
{
|
|
CAF_PDM_InitScriptableObject( "RimCornerPointCase", ":/EclipseInput48x48.png" );
|
|
setReservoirData( new RigEclipseCaseData( this ) );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RimCornerPointCase::~RimCornerPointCase()
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RimCornerPointCase::openEclipseGridFile()
|
|
{
|
|
return true;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::expected<RimCornerPointCase*, QString> RimCornerPointCase::createFromCoordinatesArray( const int nx,
|
|
const int ny,
|
|
const int nz,
|
|
const std::vector<float>& coord,
|
|
const std::vector<float>& zcorn,
|
|
const std::vector<float>& actnum )
|
|
{
|
|
CAF_ASSERT( nx > 0 );
|
|
CAF_ASSERT( ny > 0 );
|
|
CAF_ASSERT( nz > 0 );
|
|
|
|
size_t ncoord = ( nx + 1 ) * ( ny + 1 ) * 2 * 3;
|
|
size_t nzcorn = nx * ny * nz * 8;
|
|
size_t ntot = nx * ny * nz;
|
|
|
|
if ( coord.size() != ncoord )
|
|
return std::unexpected( QString( "Wrong size of coord array. Expected %1, but got %2" ).arg( ncoord ).arg( coord.size() ) );
|
|
|
|
if ( zcorn.size() != nzcorn )
|
|
return std::unexpected( QString( "Wrong size of zcorn array. Expected %1, but got %2" ).arg( nzcorn ).arg( zcorn.size() ) );
|
|
|
|
if ( actnum.size() != ntot )
|
|
return std::unexpected( QString( "Wrong size of actnum array. Expected %1, but got %2" ).arg( ntot ).arg( actnum.size() ) );
|
|
|
|
auto cornerPointCase = new RimCornerPointCase;
|
|
|
|
buildGrid( *cornerPointCase->eclipseCaseData(), nx, ny, nz, coord, zcorn, actnum );
|
|
cornerPointCase->computeCachedData();
|
|
computeDepthRelatedResults( *cornerPointCase );
|
|
|
|
return cornerPointCase;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::expected<void, QString> RimCornerPointCase::replaceGridFromCoordinatesArray( RimCornerPointCase& cornerPointCase,
|
|
const int nx,
|
|
const int ny,
|
|
const int nz,
|
|
const std::vector<float>& coord,
|
|
const std::vector<float>& zcorn,
|
|
const std::vector<float>& actnum )
|
|
{
|
|
RigActiveCellInfo* activeCellInfo = cornerPointCase.eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
|
|
CAF_ASSERT( activeCellInfo );
|
|
activeCellInfo->clear();
|
|
|
|
RigActiveCellInfo* fractureActiveCellInfo =
|
|
cornerPointCase.eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::FRACTURE_MODEL );
|
|
CAF_ASSERT( fractureActiveCellInfo );
|
|
fractureActiveCellInfo->clear();
|
|
|
|
RimReloadCaseTools::clearAllGridData( cornerPointCase.eclipseCaseData() );
|
|
|
|
// Clear the existing grid geometry before building the new grid
|
|
RigMainGrid* mainGrid = cornerPointCase.eclipseCaseData()->mainGrid();
|
|
CAF_ASSERT( mainGrid );
|
|
mainGrid->reservoirCells().clear();
|
|
mainGrid->clearNodes();
|
|
|
|
buildGrid( *cornerPointCase.eclipseCaseData(), nx, ny, nz, coord, zcorn, actnum );
|
|
cornerPointCase.computeCachedData();
|
|
computeDepthRelatedResults( cornerPointCase );
|
|
|
|
return {};
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::array<cvf::Vec3d, 8> RimCornerPointCase::getCorners( const RigMainGrid& grid,
|
|
const std::vector<float>& coord,
|
|
const std::vector<float>& zcorn,
|
|
const size_t cellIdx,
|
|
const cvf::Vec3d& offset,
|
|
const cvf::Vec3d& scale )
|
|
{
|
|
size_t i;
|
|
size_t j;
|
|
size_t k;
|
|
grid.ijkFromCellIndex( cellIdx, &i, &j, &k );
|
|
|
|
const size_t nx = grid.cellCountI();
|
|
const size_t ny = grid.cellCountJ();
|
|
|
|
// Get depths from zcorn
|
|
std::array<size_t, 8> zind;
|
|
zind[0] = ( k * nx * ny * 8 + j * nx * 4 + i * 2 );
|
|
zind[1] = ( zind[0] + 1 );
|
|
zind[2] = ( zind[0] + nx * 2 );
|
|
zind[3] = ( zind[2] + 1 );
|
|
|
|
for ( size_t n = 0; n < 4; n++ )
|
|
zind[n + 4] = ( zind[n] + nx * ny * 4 );
|
|
|
|
std::array<cvf::Vec3d, 8> corners;
|
|
for ( size_t n = 0; n < 8; n++ )
|
|
corners[n].z() = zcorn[zind[n]];
|
|
|
|
// calculate indices for grid pillars in COORD array
|
|
std::array<size_t, 4> pind;
|
|
pind[0] = j * ( nx + 1 ) * 6 + i * 6;
|
|
pind[1] = pind[0] + 6;
|
|
pind[2] = pind[0] + ( nx + 1 ) * 6;
|
|
pind[3] = pind[2] + 6;
|
|
|
|
for ( size_t n = 0; n < 4; n++ )
|
|
{
|
|
const double zt = coord[pind[n] + 2];
|
|
const double zb = coord[pind[n] + 5];
|
|
|
|
const double xt = coord[pind[n]];
|
|
const double yt = coord[pind[n] + 1];
|
|
const double xb = coord[pind[n] + 3];
|
|
const double yb = coord[pind[n] + 4];
|
|
|
|
const double diffZ = zt - zb;
|
|
|
|
if ( diffZ == 0.0 )
|
|
{
|
|
corners[n].x() = xt;
|
|
corners[n + 4].x() = xt;
|
|
corners[n].y() = yt;
|
|
corners[n + 4].y() = yt;
|
|
}
|
|
else
|
|
{
|
|
corners[n].x() = xt + ( xb - xt ) / diffZ * ( zt - corners[n].z() );
|
|
corners[n + 4].x() = xt + ( xb - xt ) / diffZ * ( zt - corners[n + 4].z() );
|
|
corners[n].y() = yt + ( yb - yt ) / diffZ * ( zt - corners[n].z() );
|
|
corners[n + 4].y() = yt + ( yb - yt ) / diffZ * ( zt - corners[n + 4].z() );
|
|
}
|
|
}
|
|
|
|
return corners;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimCornerPointCase::buildGrid( RigEclipseCaseData& eclipseCaseData,
|
|
const int nx,
|
|
const int ny,
|
|
const int nz,
|
|
const std::vector<float>& coord,
|
|
const std::vector<float>& zcorn,
|
|
const std::vector<float>& actnum )
|
|
{
|
|
auto startTime = high_resolution_clock::now();
|
|
|
|
RigActiveCellInfo* activeCellInfo = eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
|
|
CAF_ASSERT( activeCellInfo );
|
|
|
|
RigActiveCellInfo* fractureActiveCellInfo = eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::FRACTURE_MODEL );
|
|
CAF_ASSERT( fractureActiveCellInfo );
|
|
|
|
RigMainGrid* mainGrid = eclipseCaseData.mainGrid();
|
|
CAF_ASSERT( mainGrid );
|
|
|
|
mainGrid->setCellCounts( cvf::Vec3st( nx, ny, nz ) );
|
|
|
|
mainGrid->setGridName( "Main grid" );
|
|
|
|
size_t totalCellCount = nx * ny * nz;
|
|
|
|
activeCellInfo->setGridCount( 1 );
|
|
fractureActiveCellInfo->setGridCount( 1 );
|
|
|
|
activeCellInfo->setReservoirCellCount( totalCellCount );
|
|
fractureActiveCellInfo->setReservoirCellCount( totalCellCount );
|
|
|
|
// Reserve room for the cells and nodes and fill them with data
|
|
mainGrid->reservoirCells().reserve( totalCellCount );
|
|
|
|
int progTicks = 100;
|
|
caf::ProgressInfo progInfo( progTicks, "" );
|
|
|
|
int cellCount = static_cast<int>( totalCellCount );
|
|
size_t cellStartIndex = mainGrid->reservoirCells().size();
|
|
size_t nodeStartIndex = mainGrid->nodes().size();
|
|
|
|
RigCell defaultCell;
|
|
defaultCell.setHostGrid( mainGrid );
|
|
mainGrid->reservoirCells().resize( cellStartIndex + cellCount, defaultCell );
|
|
|
|
std::span<cvf::Vec3d> gridNodes = mainGrid->allocateNodes( static_cast<size_t>( cellCount ) * 8 );
|
|
|
|
const size_t cellMappingECLRi[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
|
|
|
|
cvf::Vec3d offset( 0.0, 0.0, 0.0 );
|
|
cvf::Vec3d scale( 1.0, 1.0, 1.0 );
|
|
|
|
// Convert actnum to integer
|
|
std::vector<int> activeCells( nx * ny * nz, 0 );
|
|
CAF_ASSERT( activeCells.size() == actnum.size() );
|
|
for ( size_t i = 0; i < activeCells.size(); i++ )
|
|
{
|
|
activeCells[i] = static_cast<int>( actnum[i] > 0.0 );
|
|
}
|
|
|
|
// Precompute the active cell matrix index
|
|
size_t numActiveCells = RifRoffFileTools::computeActiveCellMatrixIndex( activeCells );
|
|
|
|
// Loop over cells and fill them with data
|
|
#pragma omp parallel for
|
|
for ( int gridLocalCellIndex = 0; gridLocalCellIndex < cellCount; ++gridLocalCellIndex )
|
|
{
|
|
RigCell& cell = mainGrid->cell( cellStartIndex + gridLocalCellIndex );
|
|
|
|
cell.setGridLocalCellIndex( gridLocalCellIndex );
|
|
|
|
// Active cell index
|
|
int matrixActiveIndex = activeCells[gridLocalCellIndex];
|
|
if ( matrixActiveIndex != -1 )
|
|
{
|
|
activeCellInfo->setCellResultIndex( ReservoirCellIndex( cellStartIndex + gridLocalCellIndex ),
|
|
ActiveCellIndex( matrixActiveIndex ) );
|
|
}
|
|
|
|
cell.setParentCellIndex( cvf::UNDEFINED_SIZE_T );
|
|
|
|
std::array<cvf::Vec3d, 8> corners = getCorners( *mainGrid, coord, zcorn, gridLocalCellIndex, offset, scale );
|
|
|
|
// Corner coordinates
|
|
for ( int cIdx = 0; cIdx < 8; ++cIdx )
|
|
{
|
|
double* point = gridNodes[(size_t)gridLocalCellIndex * 8 + cellMappingECLRi[cIdx]].ptr();
|
|
|
|
point[0] = corners[cIdx].x();
|
|
point[1] = corners[cIdx].y();
|
|
point[2] = -corners[cIdx].z();
|
|
|
|
cell.cornerIndices()[cIdx] = nodeStartIndex + (size_t)gridLocalCellIndex * 8 + cIdx;
|
|
}
|
|
|
|
// Mark inactive long pyramid looking cells as invalid
|
|
cell.setInvalid( cell.isLongPyramidCell() );
|
|
}
|
|
|
|
activeCellInfo->setGridActiveCellCounts( 0, numActiveCells );
|
|
fractureActiveCellInfo->setGridActiveCellCounts( 0, 0 );
|
|
|
|
mainGrid->initAllSubGridsParentGridPointer();
|
|
activeCellInfo->computeDerivedData();
|
|
fractureActiveCellInfo->computeDerivedData();
|
|
|
|
createActnumResult( eclipseCaseData );
|
|
|
|
auto endTime = high_resolution_clock::now();
|
|
|
|
auto totalDuration = duration_cast<milliseconds>( endTime - startTime );
|
|
RiaLogging::info( std::format( "Total: {} ms", totalDuration.count() ) );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Materialize an ACTNUM static result so the created grid exposes the same queryable property as a
|
|
/// grid loaded from an Eclipse file. The corner point case has no reader interface, so the data
|
|
/// cannot be loaded lazily and must be created here.
|
|
///
|
|
/// Why every value is 1.0: a STATIC_NATIVE result is stored in the per-active-cell layout, i.e. the
|
|
/// array has exactly one slot per active cell, indexed by the active-cell result index. Inactive
|
|
/// cells (input actnum <= 0) are already excluded upstream in buildGrid(): they receive no result
|
|
/// index and are not counted by reservoirActiveCellCount(). Inactivity is therefore encoded by the
|
|
/// absence of a slot, not by storing a 0 here, so every slot that exists corresponds to an active
|
|
/// cell and must hold 1.0. This is the same thing the native readers do
|
|
/// (RifReaderEclipseOutput/RifReaderOpmCommon::staticResult: resize(reservoirActiveCellCount(), 1.0)).
|
|
/// Walking RigActiveCellInfo to set values per cell would only be needed for a per-reservoir-cell
|
|
/// layout, where inactive cells need an explicit 0 — that is not the layout used for ACTNUM.
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimCornerPointCase::createActnumResult( RigEclipseCaseData& eclipseCaseData )
|
|
{
|
|
RigActiveCellInfo* activeCellInfo = eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
|
|
CAF_ASSERT( activeCellInfo );
|
|
|
|
RigCaseCellResultsData* matrixResults = eclipseCaseData.results( RiaDefines::PorosityModelType::MATRIX_MODEL );
|
|
CAF_ASSERT( matrixResults );
|
|
|
|
RigEclipseResultAddress resAddr( RiaDefines::ResultCatType::STATIC_NATIVE, "ACTNUM" );
|
|
if ( matrixResults->hasResultEntry( resAddr ) ) return;
|
|
|
|
// One slot per active cell; all active by construction (see comment above), so fill with 1.0.
|
|
size_t activeCellCount = activeCellInfo->reservoirActiveCellCount();
|
|
matrixResults->addStaticScalarResult( RiaDefines::ResultCatType::STATIC_NATIVE, "ACTNUM", false, activeCellCount );
|
|
|
|
auto modifiableData = matrixResults->modifiableCellScalarResultTimesteps( resAddr );
|
|
CAF_ASSERT( modifiableData && !modifiableData->empty() );
|
|
std::fill( ( *modifiableData )[0].begin(), ( *modifiableData )[0].end(), 1.0 );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Compute the depth related geometry properties (DEPTH, DX, DY, DZ, TOPS, BOTTOM), matching the
|
|
/// behavior when a grid is imported from file.
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimCornerPointCase::computeDepthRelatedResults( RimCornerPointCase& cornerPointCase )
|
|
{
|
|
if ( RiaPreferencesGrid::current()->autoComputeDepthRelatedProperties() )
|
|
{
|
|
cornerPointCase.eclipseCaseData()->computeDepthRelatedResults();
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RimCornerPointCase::defineUiOrdering( QString uiConfigName, caf::PdmUiOrdering& uiOrdering )
|
|
{
|
|
uiOrdering.add( &m_caseUserDescription );
|
|
uiOrdering.add( &m_displayNameOption );
|
|
uiOrdering.add( &m_caseId );
|
|
uiOrdering.add( &m_caseFileName );
|
|
|
|
auto group = uiOrdering.addNewGroup( "Case Options" );
|
|
group->add( &m_activeFormationNames );
|
|
group->add( &m_flipXAxis );
|
|
group->add( &m_flipYAxis );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
QString RimCornerPointCase::locationOnDisc() const
|
|
{
|
|
if ( gridFileName().isEmpty() ) return QString();
|
|
|
|
QFileInfo fi( gridFileName() );
|
|
return fi.absolutePath();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RimCornerPointCase::importAsciiInputProperties( const QStringList& fileNames )
|
|
{
|
|
return true;
|
|
}
|