Support for loading only active cell geometry (#11624)

* Only load active cells for main grid, skip LGRs for now
* Handle wells with inactive cells
* Validate mapaxes transform before using it.
* Add log message
* Additional guarding when trying to find the geometrical location of a simulation cell
* Add extra safeguarding for init/restart file access in opm common. Only support unified restart files.
This commit is contained in:
jonjenssen
2024-08-28 18:22:57 +02:00
committed by GitHub
parent 27c46a65fd
commit 0572069511
24 changed files with 670 additions and 78 deletions

View File

@@ -86,6 +86,7 @@ set(SOURCE_GROUP_HEADER_FILES
${CMAKE_CURRENT_LIST_DIR}/RifRevealCsvSectionSummaryReader.h
${CMAKE_CURRENT_LIST_DIR}/RifStimPlanCsvSummaryReader.h
${CMAKE_CURRENT_LIST_DIR}/RifReaderOpmCommon.h
${CMAKE_CURRENT_LIST_DIR}/RifReaderOpmCommonActive.h
${CMAKE_CURRENT_LIST_DIR}/RifEclipseReportKeywords.h
${CMAKE_CURRENT_LIST_DIR}/RifInpExportTools.h
${CMAKE_CURRENT_LIST_DIR}/RifFaultReactivationModelExporter.h
@@ -188,6 +189,7 @@ set(SOURCE_GROUP_SOURCE_FILES
${CMAKE_CURRENT_LIST_DIR}/RifRevealCsvSectionSummaryReader.cpp
${CMAKE_CURRENT_LIST_DIR}/RifStimPlanCsvSummaryReader.cpp
${CMAKE_CURRENT_LIST_DIR}/RifReaderOpmCommon.cpp
${CMAKE_CURRENT_LIST_DIR}/RifReaderOpmCommonActive.cpp
${CMAKE_CURRENT_LIST_DIR}/RifInpExportTools.cpp
${CMAKE_CURRENT_LIST_DIR}/RifFaultReactivationModelExporter.cpp
${CMAKE_CURRENT_LIST_DIR}/RifThermalToStimPlanFractureXmlOutput.cpp

View File

@@ -132,6 +132,11 @@ size_t RifReaderEclipseWell::localGridCellIndexFromErtConnection( const RigGridB
return cvf::UNDEFINED_SIZE_T;
}
if ( ( cellI < 0 ) || ( cellJ < 0 ) )
{
return cvf::UNDEFINED_SIZE_T;
}
return grid->cellIndexFromIJK( cellI, cellJ, cellK );
}
@@ -158,6 +163,17 @@ RigWellResultPoint RifReaderEclipseWell::createWellResultPoint( const RigEclipse
RigWellResultPoint resultPoint;
if ( ( grid->cellCount() == 0 ) || ( gridCellIndex > grid->cellCount() - 1 ) )
{
return resultPoint;
}
const RigCell& c = grid->cell( gridCellIndex );
if ( c.isInvalid() )
{
return resultPoint;
}
if ( gridCellIndex != cvf::UNDEFINED_SIZE_T )
{
int branchId = -1, segmentId = -1, outletBranchId = -1, outletSegmentId = -1;
@@ -857,11 +873,17 @@ void RifReaderEclipseWell::readWellCells( RifEclipseRestartDataAccess* restartDa
{
prevResPoint = wellResFrame.wellHead();
}
else
else if ( rpIdx > 0 )
{
prevResPoint = wellResultBranch.branchResultPoints()[rpIdx - 1];
}
if ( !prevResPoint.isCell() )
{
// When importing only active cells, this situation can occur if the well head is a inactive cell.
continue;
}
cvf::Vec3d lastConnectionPos = grids[prevResPoint.gridIndex()]->cell( prevResPoint.cellIndex() ).center();
SegmentPositionContribution

View File

@@ -80,6 +80,14 @@ const QString RifReaderInterface::faultIncludeFileAbsolutePathPrefix() const
return m_readerSettings.includeFileAbsolutePathPrefix;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RifReaderInterface::onlyLoadActiveCells() const
{
return m_readerSettings.onlyLoadActiveCells;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

View File

@@ -57,6 +57,7 @@ public:
bool includeInactiveCellsInFaultGeometry() const;
bool loadWellDataEnabled() const;
const QString faultIncludeFileAbsolutePathPrefix() const;
bool onlyLoadActiveCells() const;
void setReaderSettings( RifReaderSettings readerSettings );

View File

@@ -29,6 +29,7 @@
#include "RifOpmRadialGridTools.h"
#include "RifReaderEclipseWell.h"
#include "RigActiveCellGrid.h"
#include "RigActiveCellInfo.h"
#include "RigCaseCellResultsData.h"
#include "RigEclipseCaseData.h"
@@ -87,18 +88,16 @@ bool RifReaderOpmCommon::open( const QString& fileName, RigEclipseCaseData* ecli
return false;
}
if ( isFaultImportEnabled() )
{
auto task = progress.task( "Reading faults", 25 );
if ( isFaultImportEnabled() )
{
cvf::Collection<RigFault> faults;
cvf::Collection<RigFault> faults;
importFaults( fileSet, &faults );
importFaults( fileSet, &faults );
RigMainGrid* mainGrid = eclipseCaseData->mainGrid();
mainGrid->setFaults( faults );
}
RigMainGrid* mainGrid = eclipseCaseData->mainGrid();
mainGrid->setFaults( faults );
}
m_eclipseCaseData = eclipseCaseData;
@@ -108,9 +107,10 @@ bool RifReaderOpmCommon::open( const QString& fileName, RigEclipseCaseData* ecli
buildMetaData( eclipseCaseData, progress );
}
auto task = progress.task( "Handling NCC Result data", 25 );
if ( isNNCsEnabled() )
{
auto task = progress.task( "Handling NCC Result data", 25 );
caf::ProgressInfo nncProgress( 10, "" );
RigMainGrid* mainGrid = eclipseCaseData->mainGrid();
@@ -278,9 +278,11 @@ bool RifReaderOpmCommon::importGrid( RigMainGrid* mainGrid, RigEclipseCaseData*
mapAxes[i] = opmMapAxes[i];
}
double norm_denominator = mapAxes[2] * mapAxes[5] - mapAxes[4] * mapAxes[3];
// Set the map axes transformation matrix on the main grid
mainGrid->setMapAxes( mapAxes );
mainGrid->setUseMapAxes( true );
mainGrid->setUseMapAxes( norm_denominator != 0.0 );
auto transform = mainGrid->mapAxisTransform();
@@ -414,6 +416,7 @@ void RifReaderOpmCommon::transferGeometry( Opm::EclIO::EGrid& opmMainGrid,
RigCell defaultCell;
defaultCell.setHostGrid( localGrid );
mainGrid->globalCellArray().resize( cellStartIndex + cellCount, defaultCell );
mainGrid->nodes().resize( nodeStartIndex + cellCount * 8, cvf::Vec3d( 0, 0, 0 ) );
@@ -724,12 +727,26 @@ void RifReaderOpmCommon::setupInitAndRestartAccess()
{
if ( ( m_initFile == nullptr ) && !m_initFileName.empty() )
{
m_initFile = std::make_unique<EclIO::EInit>( m_initFileName );
try
{
m_initFile = std::make_unique<EclIO::EInit>( m_initFileName );
}
catch ( ... )
{
m_initFile = nullptr;
}
}
if ( ( m_restartFile == nullptr ) && !m_restartFileName.empty() )
{
m_restartFile = std::make_unique<EclIO::ERst>( m_restartFileName );
try
{
m_restartFile = std::make_unique<EclIO::ERst>( m_restartFileName );
}
catch ( ... )
{
m_restartFile = nullptr;
}
}
}

View File

@@ -32,6 +32,7 @@ class EGrid;
} // namespace Opm::EclIO
class RigMainGrid;
class RigActiveCellGrid;
class RigGridBase;
class RigEclipseCaseData;
class RigEclipseTimeStepInfo;
@@ -58,36 +59,38 @@ public:
std::vector<QDateTime> timeStepsOnFile( QString gridFileName );
private:
void buildMetaData( RigEclipseCaseData* caseData, caf::ProgressInfo& progress );
bool importGrid( RigMainGrid* mainGrid, RigEclipseCaseData* caseData );
void transferGeometry( Opm::EclIO::EGrid& opmMainGrid,
Opm::EclIO::EGrid& opmGrid,
RigMainGrid* riMainGrid,
RigGridBase* riGrid,
RigEclipseCaseData* caseData );
protected:
virtual bool importGrid( RigMainGrid* mainGrid, RigEclipseCaseData* caseData );
void transferActiveCells( Opm::EclIO::EGrid& opmGrid,
size_t cellStartIndex,
RigEclipseCaseData* eclipseCaseData,
size_t matrixActiveStartIndex,
size_t fractureActiveStartIndex );
void transferStaticNNCData( Opm::EclIO::EGrid& opmMainGrid, std::vector<Opm::EclIO::EGrid>& lgrGrids, RigMainGrid* mainGrid );
bool verifyActiveCellInfo( int activeSizeMat, int activeSizeFrac );
void updateActiveCellInfo( RigEclipseCaseData* eclipseCaseData,
Opm::EclIO::EGrid& opmGrid,
std::vector<Opm::EclIO::EGrid>& lgrGrids,
RigMainGrid* mainGrid );
private:
void buildMetaData( RigEclipseCaseData* caseData, caf::ProgressInfo& progress );
std::vector<RigEclipseTimeStepInfo> createFilteredTimeStepInfos();
std::vector<std::vector<int>> readActiveCellInfoFromPorv( RigEclipseCaseData* eclipseCaseData, bool isDualPorosity );
void transferGeometry( Opm::EclIO::EGrid& opmMainGrid,
Opm::EclIO::EGrid& opmGrid,
RigMainGrid* riMainGrid,
RigGridBase* riGrid,
RigEclipseCaseData* caseData );
void transferDynamicNNCData( RigMainGrid* mainGrid );
void locateInitAndRestartFiles( QString gridFileName );
void setupInitAndRestartAccess();
std::vector<RigEclipseTimeStepInfo> createFilteredTimeStepInfos();
bool verifyActiveCellInfo( int activeSizeMat, int activeSizeFrac );
std::vector<std::vector<int>> readActiveCellInfoFromPorv( RigEclipseCaseData* eclipseCaseData, bool isDualPorosity );
void updateActiveCellInfo( RigEclipseCaseData* eclipseCaseData,
Opm::EclIO::EGrid& opmGrid,
std::vector<Opm::EclIO::EGrid>& lgrGrids,
RigMainGrid* mainGrid );
struct TimeDataFile
{
int sequenceNumber;
@@ -99,22 +102,22 @@ private:
std::vector<TimeDataFile> readTimeSteps();
private:
protected:
enum class ActiveType
{
ACTIVE_MATRIX_VALUE = 1,
ACTIVE_FRACTURE_VALUE = 2
};
std::string m_gridFileName;
std::string m_initFileName;
std::string m_restartFileName;
int m_gridUnit;
std::string m_gridFileName;
int m_gridUnit;
std::vector<std::string> m_gridNames;
RigEclipseCaseData* m_eclipseCaseData;
private:
std::string m_initFileName;
std::string m_restartFileName;
std::unique_ptr<Opm::EclIO::ERst> m_restartFile;
std::unique_ptr<Opm::EclIO::EInit> m_initFile;
std::vector<std::string> m_gridNames;
};

View File

@@ -0,0 +1,235 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2024 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 "RifReaderOpmCommonActive.h"
#include "RiaLogging.h"
#include "RiaStdStringTools.h"
#include "RifOpmRadialGridTools.h"
#include "RigActiveCellGrid.h"
#include "RigActiveCellInfo.h"
#include "RigEclipseCaseData.h"
#include "cafProgressInfo.h"
#include "opm/io/eclipse/EGrid.hpp"
#include <QStringList>
using namespace Opm;
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RifReaderOpmCommonActive::RifReaderOpmCommonActive()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RifReaderOpmCommonActive::~RifReaderOpmCommonActive()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RifReaderOpmCommonActive::importGrid( RigMainGrid* /* mainGrid*/, RigEclipseCaseData* eclipseCaseData )
{
RigActiveCellGrid* activeGrid = new RigActiveCellGrid();
eclipseCaseData->setMainGrid( activeGrid );
caf::ProgressInfo progInfo( 5, "Importing Eclipse Grid" );
Opm::EclIO::EGrid opmGrid( m_gridFileName );
const auto& dims = opmGrid.dimension();
activeGrid->setGridPointDimensions( cvf::Vec3st( dims[0] + 1, dims[1] + 1, dims[2] + 1 ) );
activeGrid->setGridName( "Main grid" );
activeGrid->setDualPorosity( opmGrid.porosity_mode() > 0 );
// assign grid unit, if found (1 = Metric, 2 = Field, 3 = Lab)
auto gridUnitStr = RiaStdStringTools::toUpper( opmGrid.grid_unit() );
if ( gridUnitStr.starts_with( 'M' ) )
m_gridUnit = 1;
else if ( gridUnitStr.starts_with( 'F' ) )
m_gridUnit = 2;
else if ( gridUnitStr.starts_with( 'C' ) )
m_gridUnit = 3;
auto globalMatrixActiveSize = opmGrid.activeCells();
auto globalFractureActiveSize = opmGrid.activeFracCells();
m_gridNames.clear();
m_gridNames.push_back( "global" );
std::vector<Opm::EclIO::EGrid> lgrGrids; // lgrs not supported here for now
// active cell information
{
auto task = progInfo.task( "Getting Active Cell Information", 1 );
// in case init file and grid file disagrees with number of active cells, read extra porv information from init file to correct this
if ( !verifyActiveCellInfo( globalMatrixActiveSize, globalFractureActiveSize ) )
{
updateActiveCellInfo( eclipseCaseData, opmGrid, lgrGrids, activeGrid );
}
activeGrid->transferActiveInformation( eclipseCaseData,
opmGrid.totalActiveCells(),
opmGrid.activeCells(),
opmGrid.activeFracCells(),
opmGrid.active_indexes(),
opmGrid.active_frac_indexes() );
}
// grid geometry
{
RiaLogging::info( QString( "Loading %0 active of %1 total cells." )
.arg( QString::fromStdString( RiaStdStringTools::formatThousandGrouping( opmGrid.totalActiveCells() ) ) )
.arg( QString::fromStdString( RiaStdStringTools::formatThousandGrouping( opmGrid.totalNumberOfCells() ) ) ) );
auto task = progInfo.task( "Loading Active Cell Main Grid Geometry", 1 );
transferActiveGeometry( opmGrid, activeGrid, eclipseCaseData );
}
activeGrid->initAllSubGridsParentGridPointer();
if ( isNNCsEnabled() )
{
auto task = progInfo.task( "Loading NNC data", 1 );
transferStaticNNCData( opmGrid, lgrGrids, activeGrid );
}
auto opmMapAxes = opmGrid.get_mapaxes();
if ( opmMapAxes.size() == 6 )
{
std::array<double, 6> mapAxes;
for ( size_t i = 0; i < opmMapAxes.size(); ++i )
{
mapAxes[i] = opmMapAxes[i];
}
double norm_denominator = mapAxes[2] * mapAxes[5] - mapAxes[4] * mapAxes[3];
// Set the map axes transformation matrix on the main grid
activeGrid->setMapAxes( mapAxes );
activeGrid->setUseMapAxes( norm_denominator != 0.0 );
auto transform = activeGrid->mapAxisTransform();
// Invert the transformation matrix to convert from file coordinates to domain coordinates
transform.invert();
#pragma omp parallel for
for ( long i = 0; i < static_cast<long>( activeGrid->nodes().size() ); i++ )
{
auto& n = activeGrid->nodes()[i];
n.transformPoint( transform );
}
}
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RifReaderOpmCommonActive::transferActiveGeometry( Opm::EclIO::EGrid& opmMainGrid,
RigActiveCellGrid* activeGrid,
RigEclipseCaseData* eclipseCaseData )
{
int cellCount = opmMainGrid.totalActiveCells();
RigCell defaultCell;
defaultCell.setHostGrid( activeGrid );
for ( size_t i = 0; i < 8; i++ )
defaultCell.cornerIndices()[i] = 0;
activeGrid->globalCellArray().resize( cellCount + 1, defaultCell );
activeGrid->globalCellArray()[cellCount].setInvalid( true );
activeGrid->nodes().resize( ( cellCount + 1 ) * 8, cvf::Vec3d( 0, 0, 0 ) );
auto& riNodes = activeGrid->nodes();
opmMainGrid.loadData();
opmMainGrid.load_grid_data();
const bool isRadialGrid = opmMainGrid.is_radial();
const auto& activeMatIndexes = opmMainGrid.active_indexes();
const auto& activeFracIndexes = opmMainGrid.active_frac_indexes();
// Compute the center of the LGR radial grid cells for each K layer
auto radialGridCenterTopLayerOpm = isRadialGrid
? RifOpmRadialGridTools::computeXyCenterForTopOfCells( opmMainGrid, opmMainGrid, activeGrid )
: std::map<int, std::pair<double, double>>();
// use same mapping as resdata
const size_t cellMappingECLRi[8] = { 0, 1, 3, 2, 4, 5, 7, 6 };
#pragma omp parallel for
for ( int opmCellIndex = 0; opmCellIndex < static_cast<int>( opmMainGrid.totalNumberOfCells() ); opmCellIndex++ )
{
if ( ( activeMatIndexes[opmCellIndex] < 0 ) && ( activeFracIndexes[opmCellIndex] < 0 ) ) continue;
auto opmIJK = opmMainGrid.ijk_from_global_index( opmCellIndex );
double xCenterCoordOpm = 0.0;
double yCenterCoordOpm = 0.0;
if ( isRadialGrid && radialGridCenterTopLayerOpm.contains( opmIJK[2] ) )
{
const auto& [xCenter, yCenter] = radialGridCenterTopLayerOpm[opmIJK[2]];
xCenterCoordOpm = xCenter;
yCenterCoordOpm = yCenter;
}
auto riReservoirIndex = activeGrid->cellIndexFromIJK( opmIJK[0], opmIJK[1], opmIJK[2] );
RigCell& cell = activeGrid->globalCellArray()[riReservoirIndex];
cell.setGridLocalCellIndex( riReservoirIndex );
cell.setParentCellIndex( cvf::UNDEFINED_SIZE_T );
// corner coordinates
std::array<double, 8> opmX{};
std::array<double, 8> opmY{};
std::array<double, 8> opmZ{};
opmMainGrid.getCellCorners( opmCellIndex, opmX, opmY, opmZ );
// Each cell has 8 nodes, use reservoir cell index and multiply to find first node index for cell
auto riNodeStartIndex = riReservoirIndex * 8;
for ( size_t opmNodeIndex = 0; opmNodeIndex < 8; opmNodeIndex++ )
{
auto riCornerIndex = cellMappingECLRi[opmNodeIndex];
size_t riNodeIndex = riNodeStartIndex + riCornerIndex;
auto& riNode = riNodes[riNodeIndex];
riNode.x() = opmX[opmNodeIndex] + xCenterCoordOpm;
riNode.y() = opmY[opmNodeIndex] + yCenterCoordOpm;
riNode.z() = -opmZ[opmNodeIndex];
cell.cornerIndices()[riCornerIndex] = riNodeIndex;
}
cell.setInvalid( cell.isLongPyramidCell() );
}
}

View File

@@ -0,0 +1,38 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2024- 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.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include "RifReaderOpmCommon.h"
class RigActiveCellGrid;
//==================================================================================================
//
//
//==================================================================================================
class RifReaderOpmCommonActive : public RifReaderOpmCommon
{
public:
RifReaderOpmCommonActive();
~RifReaderOpmCommonActive() override;
protected:
bool importGrid( RigMainGrid* mainGrid, RigEclipseCaseData* caseData ) override;
void transferActiveGeometry( Opm::EclIO::EGrid& opmMainGrid, RigActiveCellGrid* riMainGrid, RigEclipseCaseData* caseData );
};

View File

@@ -35,4 +35,5 @@ struct RifReaderSettings
bool skipWellData;
bool importSummaryData;
QString includeFileAbsolutePathPrefix;
bool onlyLoadActiveCells;
};