OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Result Divided by Area: Establish concept used to compute flow velocity and normalized trans (#7349)

Browse Source 
* Geometry Tools : Add convenience functions for polygon area

* #7232 Result Divided by Area: Add cell face result and show in GUI

Native support for flow rate is given by mass rate (mass per time) over a cell face. Add a derived result that takes flow rate divided by cell face area to get velocity (distance per time).

Add support for this concept on relevant native results, and indicate this result type in UI using a "/A" postfix

* Speed up divided-by-area calculations by using openmp

* Some refactoring of result data access.

* Make sure NNC data is scaled correctly in vector flow viz.

Co-authored-by: jonjenssen <jon@soundsoft.no>
This commit is contained in:
Magne Sjaastad
2021-02-11 03:01:17 +01:00
committed by GitHub
parent 424f66fd4e
commit cc292b197a
30 changed files with 631 additions and 321 deletions
Download Patch File Download Diff File Expand all files Collapse all files

134
ApplicationLibCode/ReservoirDataModel/RigCaseCellResultCalculator.cpp
View File

@@ -43,23 +43,21 @@ bool RigCaseCellResultCalculator::computeDifference( RigEclipseCaseData*
const RigEclipseResultAddress& address )
{
CVF_ASSERT( address.isValid() );
CVF_ASSERT( address.hasDifferenceCase() || address.isTimeLapse() );
CVF_ASSERT( address.isDeltaCaseActive() || address.isDeltaTimeStepActive() );
// Assume at this stage that data for the case is available
// It is up to the caller to make sure the case is read from file
RigEclipseCaseData* baseCase = sourceCase;
if ( address.hasDifferenceCase() )
if ( address.isDeltaCaseActive() )
{
auto eclipseCases = RimProject::current()->eclipseCases();
for ( RimEclipseCase* c : eclipseCases )
{
auto eclipseCases = RimProject::current()->eclipseCases();
for ( RimEclipseCase* c : eclipseCases )
if ( c && c->caseId() == address.deltaCaseId() && c->eclipseCaseData() )
{
if ( c && c->caseId() == address.m_differenceCaseId && c->eclipseCaseData() )
{
baseCase = c->eclipseCaseData();
}
baseCase = c->eclipseCaseData();
}
}
}
@@ -92,8 +90,8 @@ bool RigCaseCellResultCalculator::computeDifference( RigEclipseCaseData*
}
RigEclipseResultAddress nativeAddress( address );
nativeAddress.m_differenceCaseId = RigEclipseResultAddress::noCaseDiffValue();
nativeAddress.m_timeLapseBaseFrameIdx = RigEclipseResultAddress::noTimeLapseValue();
nativeAddress.setDeltaCaseId( RigEclipseResultAddress::noCaseDiffValue() );
nativeAddress.setDeltaTimeStepIndex( RigEclipseResultAddress::noTimeLapseValue() );
if ( !sourceCaseResults->ensureKnownResultLoaded( nativeAddress ) )
{
RiaLogging::error( "Failed to load destination diff result" );
@@ -128,7 +126,7 @@ bool RigCaseCellResultCalculator::computeDifference( RigEclipseCaseData*
size_t baseFrameCount = baseCaseResults->cellScalarResults( nativeAddress ).size();
size_t sourceFrameCount = sourceCaseResults->cellScalarResults( nativeAddress ).size();
size_t maxFrameCount = 0;
if ( address.isTimeLapse() )
if ( address.isDeltaTimeStepActive() )
{
// We have one defined time step for base case, loop over all source time steps
maxFrameCount = sourceFrameCount;
@@ -155,9 +153,9 @@ bool RigCaseCellResultCalculator::computeDifference( RigEclipseCaseData*
RigResultModifierFactory::createResultModifier( sourceCase, gridIdx, porosityModel, fIdx, address );
size_t baseFrameIdx = fIdx;
if ( address.isTimeLapse() )
if ( address.isDeltaTimeStepActive() )
{
baseFrameIdx = address.m_timeLapseBaseFrameIdx;
baseFrameIdx = address.deltaTimeStepIndex();
}
cvf::ref<RigResultAccessor> baseResultAccessor =
@@ -181,3 +179,113 @@ bool RigCaseCellResultCalculator::computeDifference( RigEclipseCaseData*
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCaseCellResultCalculator::computeDivideByCellFaceArea( RigMainGrid* mainGrid,
RigEclipseCaseData* destination,
RiaDefines::PorosityModelType porosityModel,
const RigEclipseResultAddress& address )
{
if ( !destination )
{
RiaLogging::error( "Missing input case for divide by area calculator" );
return false;
}
CVF_ASSERT( address.isValid() );
CVF_ASSERT( address.isDivideByCellFaceAreaActive() );
RigCaseCellResultsData* baseCaseResults = destination->results( porosityModel );
if ( !baseCaseResults )
{
RiaLogging::error( "Missing result data for divide by area calculator" );
return false;
}
RigEclipseResultAddress nativeAddress( address );
nativeAddress.enableDivideByCellFaceArea( false );
if ( !baseCaseResults->ensureKnownResultLoaded( nativeAddress ) )
{
RiaLogging::error( "Failed to load source case for divide by area result" );
return false;
}
// Initialize difference result with infinity for correct number of time steps and values per time step
{
const std::vector<std::vector<double>>& srcFrames = baseCaseResults->cellScalarResults( nativeAddress );
std::vector<std::vector<double>>* diffResultFrames =
baseCaseResults->modifiableCellScalarResultTimesteps( address );
diffResultFrames->resize( srcFrames.size() );
for ( size_t fIdx = 0; fIdx < srcFrames.size(); ++fIdx )
{
const std::vector<double>& srcVals = srcFrames[fIdx];
std::vector<double>& dstVals = diffResultFrames->at( fIdx );
// Clear the values, and resize with infinity as default value
dstVals.clear();
dstVals.resize( srcVals.size(), std::numeric_limits<double>::infinity() );
}
}
size_t baseFrameCount = baseCaseResults->cellScalarResults( nativeAddress ).size();
size_t maxGridCount = mainGrid->gridCount();
cvf::StructGridInterface::FaceType cellFace = cvf::StructGridInterface::NO_FACE;
QString resultName = address.resultName();
if ( resultName.contains( "I+" ) )
cellFace = cvf::StructGridInterface::POS_I;
else if ( resultName.contains( "J+" ) )
cellFace = cvf::StructGridInterface::POS_J;
else if ( resultName.contains( "K+" ) )
cellFace = cvf::StructGridInterface::POS_K;
else if ( resultName.contains( "TRANX" ) )
cellFace = cvf::StructGridInterface::POS_I;
else if ( resultName.contains( "TRANY" ) )
cellFace = cvf::StructGridInterface::POS_J;
else if ( resultName.contains( "TRANZ" ) )
cellFace = cvf::StructGridInterface::POS_K;
for ( size_t gridIdx = 0; gridIdx < maxGridCount; ++gridIdx )
{
auto grid = mainGrid->gridByIndex( gridIdx );
const RigActiveCellInfo* activeCellInfo = baseCaseResults->activeCellInfo();
for ( size_t fIdx = 0; fIdx < baseFrameCount; ++fIdx )
{
cvf::ref<RigResultAccessor> sourceResultAccessor =
RigResultAccessorFactory::createFromResultAddress( destination, gridIdx, porosityModel, fIdx, nativeAddress );
cvf::ref<RigResultModifier> resultModifier =
RigResultModifierFactory::createResultModifier( destination, gridIdx, porosityModel, fIdx, address );
#pragma omp parallel for
for ( int localGridCellIdx = 0; localGridCellIdx < static_cast<int>( grid->cellCount() ); localGridCellIdx++ )
{
const size_t reservoirCellIndex = grid->reservoirCellIndex( localGridCellIdx );
if ( activeCellInfo->isActive( reservoirCellIndex ) )
{
double sourceVal = sourceResultAccessor->cellScalar( localGridCellIdx );
const auto faceNormal = grid->cell( localGridCellIdx ).faceNormalWithAreaLength( cellFace );
const auto divisor = faceNormal.length();
const double epsilon = 1e-12;
if ( divisor > epsilon )
{
sourceVal /= divisor;
}
resultModifier->setCellScalar( localGridCellIdx, sourceVal );
}
}
}
}
return true;
}