OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2026-08-27 05:37:21 -05:00
Code Issues Packages Projects Releases Wiki Activity
Files
2ac1948abe49ebb7f95c3f550c710e0f56fff7f4
ResInsight/ApplicationLibCode/ProjectDataModel/RimCornerPointCase.cpp
T
Kristian Bendiksen a489db3e64 #14227 Compute depth-related geometry results once for both porosity models 
The depth-related geometry results (DEPTH/DX/DY/DZ/TOPS/BOTTOM) are derived
purely from the shared grid geometry and were computed independently for the
matrix and fracture porosity models. For dual-porosity cases this recomputed the
identical per-cell geometry twice.

Replace the per-model computeDepthRelatedResults() with a single static routine
that traverses the shared main grid once, computing each cell's geometry a single
time and writing it to every porosity model in which the cell is active. The
per-property already-computed guards and the temporary-grid recompute path are
preserved, so the stored values are unchanged. All six matrix+fracture call sites
now go through a thin RigEclipseCaseData::computeDepthRelatedResults() wrapper.
2026-06-19 12:25:55 +02:00

421 lines
18 KiB
C++
Raw Blame History

/////////////////////////////////////////////////////////////////////////////////
//
// 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 );
CVF_ASSERT( activeCellInfo );
activeCellInfo->clear();
RigActiveCellInfo* fractureActiveCellInfo =
cornerPointCase.eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::FRACTURE_MODEL );
CVF_ASSERT( fractureActiveCellInfo );
fractureActiveCellInfo->clear();
RimReloadCaseTools::clearAllGridData( cornerPointCase.eclipseCaseData() );
// Clear the existing grid geometry before building the new grid
RigMainGrid* mainGrid = cornerPointCase.eclipseCaseData()->mainGrid();
CVF_ASSERT( mainGrid );
mainGrid->reservoirCells().clear();
mainGrid->nodes().clear();
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 );
CVF_ASSERT( activeCellInfo );
RigActiveCellInfo* fractureActiveCellInfo = eclipseCaseData.activeCellInfo( RiaDefines::PorosityModelType::FRACTURE_MODEL );
CVF_ASSERT( fractureActiveCellInfo );
RigMainGrid* mainGrid = eclipseCaseData.mainGrid();
CVF_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 );
mainGrid->nodes().reserve( 8 * 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 );
mainGrid->nodes().resize( nodeStartIndex + static_cast<size_t>( cellCount ) * 8, cvf::Vec3d( 0, 0, 0 ) );
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 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 = mainGrid->nodes()[nodeStartIndex + (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 );
CVF_ASSERT( activeCellInfo );
RigCaseCellResultsData* matrixResults = eclipseCaseData.results( RiaDefines::PorosityModelType::MATRIX_MODEL );
CVF_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 );
CVF_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;
}
Reference in New Issue View Git Blame Copy Permalink