Skip to content

Msj fault performance #93

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
magnesj wants to merge 7 commits into
base: dev
Choose a base branch
from
+189 −52
Draft

Msj fault performance #93

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f3b65f9
Add ScopeTimer to measure duration of a scope
magnesj Jul 26, 2024
2695f22
Add scope timer for calculation of faults
magnesj Jul 26, 2024
dee18b6
Avoid sqrt
magnesj Jul 26, 2024
3eae0d8
Avoid calling helper functions, major speedup
magnesj Jul 26, 2024
6954318
Use map for speedup
magnesj Jul 26, 2024
ee9ed42
Roll out ijkFromCellIndex
magnesj Jul 26, 2024
b93a483
openmp test
magnesj Jul 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 34 additions & 0 deletions ApplicationLibCode/Application/Tools/RiaLogging.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "RiaGuiApplication.h"
#include "RiaRegressionTestRunner.h"

#include <chrono>
#include <iostream>
#include <sstream>

Expand Down Expand Up @@ -263,6 +264,15 @@ void RiaLogging::debug( const QString& message )
}
}

//--------------------------------------------------------------------------------------------------
/// Usage: RiaLogging::debugTimer( "Message", [&]() { code(); } );
//--------------------------------------------------------------------------------------------------
void RiaLogging::debugTimer( const QString& message, std::function<void()> callable )
{
RiaScopeTimer timing( message );
callable();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -404,3 +414,27 @@ std::vector<QString> RiaThreadSafeLogger::messages() const
{
return m_messages;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RiaScopeTimer::RiaScopeTimer( const QString& message )
: m_startTime( std::chrono::high_resolution_clock::now() )
, m_message( message )
{
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RiaScopeTimer::~RiaScopeTimer()
{
auto end = std::chrono::high_resolution_clock::now();

// Default unit for duration is seconds
std::chrono::duration<double> duration = end - m_startTime;

auto text = m_message + QString( " (duration : %1 seconds)" ).arg( duration.count() );

RiaLogging::debug( text );
}
29 changes: 28 additions & 1 deletion ApplicationLibCode/Application/Tools/RiaLogging.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,14 @@

#pragma once

#include <chrono>
#include <functional>
#include <memory>
#include <string>
#include <vector>

class QString;
#include <QString>

class QWidget;

enum class RILogLevel
Expand Down Expand Up @@ -52,6 +55,28 @@ class RiaLogger
virtual void debug( const char* message ) = 0;
};

//==================================================================================================
//
// Timing class for measuring time spent in a scope. Will create a debug log message and the time spent.
//
// Usage:
// {
// RiaScopeTimer timer( "Message" );
// // Code to measure
// }
//
//==================================================================================================
class RiaScopeTimer
{
public:
RiaScopeTimer( const QString& message );
~RiaScopeTimer();

private:
std::chrono::time_point<std::chrono::high_resolution_clock> m_startTime;
QString m_message;
};

//==================================================================================================
//
//
Expand All @@ -70,6 +95,8 @@ class RiaLogging
static void info( const QString& message );
static void debug( const QString& message );

static void debugTimer( const QString& message, std::function<void()> callable );

static void errorInMessageBox( QWidget* parent, const QString& title, const QString& text );

private:
Expand Down
14 changes: 7 additions & 7 deletions ApplicationLibCode/ProjectDataModel/RimEclipseCase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -719,14 +719,14 @@ void RimEclipseCase::loadAndSynchronizeInputProperties( bool importGridOrFaultDa
//--------------------------------------------------------------------------------------------------
void RimEclipseCase::ensureFaultDataIsComputed()
{
RigEclipseCaseData* rigEclipseCase = eclipseCaseData();
if ( rigEclipseCase )
if ( !m_readerSettings.importFaults ) return;

if ( RigEclipseCaseData* rigEclipseCase = eclipseCaseData() )
{
if ( m_readerSettings.importFaults )
{
RigActiveCellInfo* actCellInfo = rigEclipseCase->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
rigEclipseCase->mainGrid()->calculateFaults( actCellInfo );
}
RiaScopeTimer timing( "Calculate faults" );

RigActiveCellInfo* actCellInfo = rigEclipseCase->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL );
rigEclipseCase->mainGrid()->calculateFaults( actCellInfo );
}
}

Expand Down
147 changes: 111 additions & 36 deletions ApplicationLibCode/ReservoirDataModel/RigMainGrid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -479,6 +479,48 @@ void RigMainGrid::computeBoundingBox()
}
}

namespace StructGridDefines
{
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
consteval auto cubeFaceIndices()
{
//
// 7---------6
// /| /| |k
// / | / | | /j
// 4---------5 | |/
// | 3------|--2 *---i
// | / | /
// |/ |/
// 0---------1

return std::array<std::pair<cvf::StructGridInterface::FaceType, std::array<cvf::ubyte, 4>>, 6>{
{ { cvf::StructGridInterface::FaceType::NEG_K, { 0, 3, 2, 1 } },
{ cvf::StructGridInterface::FaceType::POS_K, { 4, 5, 6, 7 } },
{ cvf::StructGridInterface::FaceType::NEG_J, { 0, 1, 5, 4 } },
{ cvf::StructGridInterface::FaceType::POS_J, { 3, 7, 6, 2 } },
{ cvf::StructGridInterface::FaceType::NEG_I, { 0, 4, 7, 3 } },
{ cvf::StructGridInterface::FaceType::POS_I, { 1, 2, 6, 5 } } } };
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
auto createFaceIndicesMap()
{
std::map<cvf::StructGridInterface::FaceType, std::array<cvf::ubyte, 4>> faultFaceToFaceIdxs;

for ( const auto& [key, value] : cubeFaceIndices() )
{
faultFaceToFaceIdxs[key] = value;
}

return faultFaceToFaceIdxs;
}
}; // namespace StructGridDefines

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -515,8 +557,12 @@ void RigMainGrid::calculateFaults( const RigActiveCellInfo* activeCellInfo )

const std::vector<cvf::Vec3d>& vxs = m_mainGrid->nodes();

auto faultFaceToFaceIdxs = StructGridDefines::createFaceIndicesMap();

std::vector<RigFault::FaultFace>& unNamedFaultFaces = unNamedFault->faultFaces();
std::vector<RigFault::FaultFace>& unNamedFaultFacesInactive = unNamedFaultWithInactive->faultFaces();

#pragma omp parallel for
for ( int gcIdx = 0; gcIdx < static_cast<int>( m_cells.size() ); ++gcIdx )
{
addUnNamedFaultFaces( gcIdx,
Expand All @@ -526,7 +572,8 @@ void RigMainGrid::calculateFaults( const RigActiveCellInfo* activeCellInfo )
unNamedFaultWithInactiveIdx,
unNamedFaultFaces,
unNamedFaultFacesInactive,
m_faultsPrCellAcc.p() );
m_faultsPrCellAcc.p(),
faultFaceToFaceIdxs );
}
}

Expand All @@ -540,15 +587,14 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx,
int unNamedFaultWithInactiveIdx,
std::vector<RigFault::FaultFace>& unNamedFaultFaces,
std::vector<RigFault::FaultFace>& unNamedFaultFacesInactive,
RigFaultsPrCellAccumulator* faultsPrCellAcc ) const
RigFaultsPrCellAccumulator* faultsPrCellAcc,
const std::map<cvf::StructGridInterface::FaceType, std::array<cvf::ubyte, 4>>& faceMap ) const
{
if ( m_cells[gcIdx].isInvalid() )
{
return;
}

size_t neighborReservoirCellIdx;
size_t neighborGridCellIdx;
size_t i = 0;
size_t j = 0;
size_t k = 0;
Expand All @@ -564,65 +610,98 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx,
{
cvf::StructGridInterface::FaceType face = cvf::StructGridInterface::FaceType( faceIdx );

// For faces that has no used defined Fault assigned:
int faultState = -1;
// #pragma omp critical
{
faultState = faultsPrCellAcc->faultIdx( gcIdx, face );
}

if ( m_faultsPrCellAcc->faultIdx( gcIdx, face ) == RigFaultsPrCellAccumulator::NO_FAULT )
// For faces that has no used defined Fault assigned:
if ( faultState == RigFaultsPrCellAccumulator::NO_FAULT )
{
// Find neighbor cell
if ( firstNO_FAULTFaceForCell ) // To avoid doing this for every face, and only when detecting a NO_FAULT
{
size_t gridLocalCellIndex;
hostGrid = gridAndGridLocalIdxFromGlobalCellIdx( gcIdx, &gridLocalCellIndex );
size_t index = gcIdx;

i = index % m_cellCount.x();
index /= m_cellCount.x();
j = index % m_cellCount.y();
k = index / m_cellCount.y();

hostGrid->ijkFromCellIndex( gridLocalCellIndex, &i, &j, &k );
isCellActive = activeCellInfo->isActive( gcIdx );

firstNO_FAULTFaceForCell = false;
}

if ( !hostGrid->cellIJKNeighbor( i, j, k, face, &neighborGridCellIdx ) )
{
continue;
}
size_t ni, nj, nk;
hostGrid->neighborIJKAtCellFace( i, j, k, face, &ni, &nj, &nk );

neighborReservoirCellIdx = hostGrid->reservoirCellIndex( neighborGridCellIdx );
if ( m_cells[neighborReservoirCellIdx].isInvalid() )
if ( ni >= m_cellCount.x() || nj >= m_cellCount.y() || nk >= m_cellCount.z() )
{
continue;
}

bool isNeighborCellActive = activeCellInfo->isActive( neighborReservoirCellIdx );
const size_t nbGridLocalCellIndex = ni + nj * m_cellCount.x() + nk * m_cellCount.x() * m_cellCount.y();
const size_t neighborReservoirCellIdx = m_indexToStartOfCells + nbGridLocalCellIndex;
const auto& nbCell = m_mainGrid->globalCellArray()[neighborReservoirCellIdx];

if ( nbCell.isInvalid() ) continue;

const double threshold = 1e-6;
const double tol = threshold * threshold;

double tolerance = 1e-6;
// cvf::ubyte faceNodeIndices[4];
size_t faceIdxs[4];

std::array<size_t, 4> faceIdxs;
m_cells[gcIdx].faceIndices( face, &faceIdxs );
std::array<size_t, 4> nbFaceIdxs;
m_cells[neighborReservoirCellIdx].faceIndices( StructGridInterface::oppositeFace( face ), &nbFaceIdxs );
// m_mainGrid->cellFaceVertexIndices( face, faceNodeIndices );
auto faceNodeIndices = faceMap.at( face );

const auto& cornerIndices = m_cells[gcIdx].cornerIndices();
faceIdxs[0] = cornerIndices[faceNodeIndices[0]];
faceIdxs[1] = cornerIndices[faceNodeIndices[1]];
faceIdxs[2] = cornerIndices[faceNodeIndices[2]];
faceIdxs[3] = cornerIndices[faceNodeIndices[3]];

// cvf::ubyte nbFaceNodeIndices[4];
size_t nbFaceIdxs[4];

// m_mainGrid->cellFaceVertexIndices( StructGridInterface::oppositeFace( face ), nbFaceNodeIndices );
auto nbFaceNodeIndices = faceMap.at( StructGridInterface::oppositeFace( face ) );

const auto& nbCornerIndices = m_cells[neighborReservoirCellIdx].cornerIndices();
nbFaceIdxs[0] = nbCornerIndices[nbFaceNodeIndices[0]];
nbFaceIdxs[1] = nbCornerIndices[nbFaceNodeIndices[1]];
nbFaceIdxs[2] = nbCornerIndices[nbFaceNodeIndices[2]];
nbFaceIdxs[3] = nbCornerIndices[nbFaceNodeIndices[3]];

bool sharedFaceVertices = true;
if ( sharedFaceVertices && vxs[faceIdxs[0]].pointDistance( vxs[nbFaceIdxs[0]] ) > tolerance ) sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[1]].pointDistance( vxs[nbFaceIdxs[3]] ) > tolerance ) sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[2]].pointDistance( vxs[nbFaceIdxs[2]] ) > tolerance ) sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[3]].pointDistance( vxs[nbFaceIdxs[1]] ) > tolerance ) sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[0]].pointDistanceSquared( vxs[nbFaceIdxs[0]] ) > tol ) sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[1]].pointDistanceSquared( vxs[nbFaceIdxs[3]] ) > tol ) sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[2]].pointDistanceSquared( vxs[nbFaceIdxs[2]] ) > tol ) sharedFaceVertices = false;
if ( sharedFaceVertices && vxs[faceIdxs[3]].pointDistanceSquared( vxs[nbFaceIdxs[1]] ) > tol ) sharedFaceVertices = false;

if ( sharedFaceVertices )
{
continue;
}

// To avoid doing this calculation for the opposite face
int faultIdx = unNamedFaultIdx;
bool isNeighborCellActive = activeCellInfo->isActive( m_indexToStartOfCells + nbGridLocalCellIndex );
int faultIdx = unNamedFaultIdx;
if ( !( isCellActive && isNeighborCellActive ) ) faultIdx = unNamedFaultWithInactiveIdx;

faultsPrCellAcc->setFaultIdx( gcIdx, face, faultIdx );
faultsPrCellAcc->setFaultIdx( neighborReservoirCellIdx, StructGridInterface::oppositeFace( face ), faultIdx );
#pragma omp critical
{
faultsPrCellAcc->setFaultIdx( gcIdx, face, faultIdx );
faultsPrCellAcc->setFaultIdx( neighborReservoirCellIdx, StructGridInterface::oppositeFace( face ), faultIdx );

// Add as fault face only if the grid index is less than the neighbors

// Add as fault face only if the grid index is less than the neighbors
size_t cellIdxLow = std::min( static_cast<size_t>( gcIdx ), neighborReservoirCellIdx );
size_t cellIdxHigh = std::max( static_cast<size_t>( gcIdx ), neighborReservoirCellIdx );

if ( static_cast<size_t>( gcIdx ) < neighborReservoirCellIdx )
{
RigFault::FaultFace ff( gcIdx, cvf::StructGridInterface::FaceType( faceIdx ), neighborReservoirCellIdx );
RigFault::FaultFace ff( cellIdxLow, cvf::StructGridInterface::FaceType( faceIdx ), cellIdxHigh );
if ( isCellActive && isNeighborCellActive )
{
unNamedFaultFaces.push_back( ff );
Expand All @@ -632,10 +711,6 @@ void RigMainGrid::addUnNamedFaultFaces( int gcIdx,
unNamedFaultFacesInactive.push_back( ff );
}
}
else
{
CVF_FAIL_MSG( "Found fault with global neighbor index less than the native index. " );
}
}
}
}
Expand Down
17 changes: 9 additions & 8 deletions ApplicationLibCode/ReservoirDataModel/RigMainGrid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,14 +74,15 @@ class RigMainGrid : public RigGridBase
cvf::Collection<RigFault>& faults();
void calculateFaults( const RigActiveCellInfo* activeCellInfo );

void addUnNamedFaultFaces( int gcIdx,
const RigActiveCellInfo* activeCellInfo,
const std::vector<cvf::Vec3d>& vxs,
int unNamedFaultIdx,
int unNamedFaultWithInactiveIdx,
std::vector<RigFault::FaultFace>& unNamedFaultFaces,
std::vector<RigFault::FaultFace>& unNamedFaultFacesInactive,
RigFaultsPrCellAccumulator* faultsPrCellAcc ) const;
void addUnNamedFaultFaces( int gcIdx,
const RigActiveCellInfo* activeCellInfo,
const std::vector<cvf::Vec3d>& vxs,
int unNamedFaultIdx,
int unNamedFaultWithInactiveIdx,
std::vector<RigFault::FaultFace>& unNamedFaultFaces,
std::vector<RigFault::FaultFace>& unNamedFaultFacesInactive,
RigFaultsPrCellAccumulator* faultsPrCellAcc,
const std::map<cvf::StructGridInterface::FaceType, std::array<cvf::ubyte, 4>>& faceMap ) const;

void distributeNNCsToFaults();

Expand Down