Skip to content

Commit

Permalink
Fault reactivation: Improve POR-Bar extraction.
Browse files Browse the repository at this point in the history
  • Loading branch information
@kriben
kriben committed Feb 9, 2024
1 parent 5037473 commit 5a3894b
Show file tree
Hide file tree
Showing 14 changed files with 163 additions and 108 deletions.
28 changes: 16 additions & 12 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@
#include "RimFaultReactivationDataAccessorVoidRatio.h"
#include "RimFaultReactivationEnums.h"
#include "RimFaultReactivationModel.h"
#include <limits>

//--------------------------------------------------------------------------------------------------
///
Expand Down Expand Up @@ -179,20 +180,24 @@ std::vector<double> RimFaultReactivationDataAccess::extractModelData( const RigF
CAF_ASSERT( it != borderSurfaceElements.end() && "Sea bed border surface does not exist" );
std::set<unsigned int> seabedElements( it->second.begin(), it->second.end() );

std::vector<double> values;

if ( nodeProperties.contains( property ) )
{
for ( auto& node : grid->dataNodes() )
int numNodes = static_cast<int>( grid->dataNodes().size() );
std::vector<double> values( numNodes, std::numeric_limits<double>::infinity() );

for ( int nodeIndex = 0; nodeIndex < numNodes; nodeIndex++ )
{
double value = accessor->valueAtPosition( node, model, gridPart );
values.push_back( value );
double value = accessor->valueAtPosition( grid->dataNodes()[nodeIndex], model, gridPart );
values[nodeIndex] = value;
}
return values;
}
else
{
size_t numElements = grid->elementIndices().size();
for ( size_t elementIndex = 0; elementIndex < numElements; elementIndex++ )
int numElements = static_cast<int>( grid->elementIndices().size() );
std::vector<double> values( numElements, std::numeric_limits<double>::infinity() );

for ( int elementIndex = 0; elementIndex < numElements; elementIndex++ )
{
std::vector<cvf::Vec3d> corners = grid->elementDataCorners( elementIndex );

Expand All @@ -203,13 +208,12 @@ std::vector<double> RimFaultReactivationDataAccess::extractModelData( const RigF
double topDepth = computeAverageDepth( corners, { 0, 1, 2, 3 } ) - topDepthAdjust;
double bottomDepth = computeAverageDepth( corners, { 4, 5, 6, 7 } );

cvf::Vec3d position = RigCaseToCaseCellMapperTools::calculateCellCenter( corners.data() );
double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth, elementIndex );
values.push_back( value );
cvf::Vec3d position = RigCaseToCaseCellMapperTools::calculateCellCenter( corners.data() );
double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth, elementIndex );
values[elementIndex] = value;
}
return values;
}

return values;
}

return {};
Expand Down
19 changes: 19 additions & 0 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include "RigFaultReactivationModel.h"
#include "RigMainGrid.h"

#include "RimFaultReactivationEnums.h"
#include "cafAssert.h"

//--------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -48,3 +49,21 @@ void RimFaultReactivationDataAccessor::setModelAndTimeStep( const RigFaultReacti
m_timeStep = timeStep;
updateResultAccessor();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<bool, RimFaultReactivation::ElementSets> RimFaultReactivationDataAccessor::findElementSetForElementIndex(
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets,
int elementIndex )
{
for ( auto [s, indexes] : elementSets )
{
if ( std::find( indexes.begin(), indexes.end(), elementIndex ) != indexes.end() )
{
return std::pair<bool, RimFaultReactivation::ElementSets>( true, s );
}
}

return std::pair<bool, RimFaultReactivation::ElementSets>( false, RimFaultReactivation::ElementSets::OverBurden );
}
5 changes: 5 additions & 0 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@

#include <limits>
#include <map>
#include <vector>

class RigMainGrid;
class RigFaultReactivationModel;
Expand Down Expand Up @@ -52,6 +53,10 @@ class RimFaultReactivationDataAccessor
protected:
virtual void updateResultAccessor() = 0;

static std::pair<bool, RimFaultReactivation::ElementSets>
findElementSetForElementIndex( const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets,
int elementIndex );

const RigFaultReactivationModel* m_model;
size_t m_timeStep;
};
24 changes: 8 additions & 16 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,18 +61,14 @@ void RimFaultReactivationDataAccessorGeoMech::updateResultAccessor()
{
const int partIndex = 0;

auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr ) -> RigFemScalarResultFrames*
auto loadFrameLambda = [&]( auto femParts, RigFemResultAddress addr ) -> std::vector<float>
{
auto result = femParts->findOrLoadScalarResult( partIndex, addr );
if ( result->frameData( 0, 0 ).empty() )
{
return nullptr;
}
return result;
return result->frameData( 0, 0 );
};

auto femParts = m_geoMechCaseData->femPartResults();
m_resultFrames = loadFrameLambda( femParts, getResultAddress( m_property ) );
auto femParts = m_geoMechCaseData->femPartResults();
m_data = loadFrameLambda( femParts, getResultAddress( m_property ) );
}

//--------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -104,23 +100,19 @@ double RimFaultReactivationDataAccessorGeoMech::valueAtPosition( const cvf::Vec3
double bottomDepth,
size_t elementIndex ) const
{
if ( !m_resultFrames ) return std::numeric_limits<double>::infinity();
if ( !m_data.empty() ) return std::numeric_limits<double>::infinity();

RimWellIADataAccess iaDataAccess( m_geoMechCase );
int elementIdx = iaDataAccess.elementIndex( position );
if ( elementIdx != -1 )
{
int timeStepIndex = 0;
int frameIndex = 0;

const std::vector<float>& data = m_resultFrames->frameData( timeStepIndex, frameIndex );
if ( elementIdx >= static_cast<int>( data.size() ) ) return std::numeric_limits<double>::infinity();
if ( elementIdx >= static_cast<int>( m_data.size() ) ) return std::numeric_limits<double>::infinity();

if ( m_property == RimFaultReactivation::Property::YoungsModulus )
{
return RiaEclipseUnitTools::gigaPascalToPascal( data[elementIdx] );
return RiaEclipseUnitTools::gigaPascalToPascal( m_data[elementIdx] );
}
return data[elementIdx];
return m_data[elementIdx];
}

return std::numeric_limits<double>::infinity();
Expand Down
2 changes: 1 addition & 1 deletion ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,5 +55,5 @@ class RimFaultReactivationDataAccessorGeoMech : public RimFaultReactivationDataA
RimGeoMechCase* m_geoMechCase;
RimFaultReactivation::Property m_property;
RigGeoMechCaseData* m_geoMechCaseData;
RigFemScalarResultFrames* m_resultFrames;
std::vector<float> m_data;
};
20 changes: 10 additions & 10 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,24 +109,24 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf:
CAF_ASSERT( m_wellPaths.find( gridPart ) != m_wellPaths.end() );
auto wellPath = m_wellPaths.find( gridPart )->second;

auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
if ( cellIdx != cvf::UNDEFINED_SIZE_T )
{
double valueFromEclipse = m_resultAccessor->cellScalar( cellIdx );
if ( !std::isinf( valueFromEclipse ) ) return RiaEclipseUnitTools::barToPascal( valueFromEclipse );
}

auto [values, intersections] =
RimFaultReactivationDataAccessorWellLogExtraction::extractValuesAndIntersections( *m_resultAccessor.p(), *extractor.p(), *wellPath );

auto [value, pos] = RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( model,
RimFaultReactivation::ElementSets elementSet = RimFaultReactivation::ElementSets::UnderBurden;
auto [value, pos] = RimFaultReactivationDataAccessorWellLogExtraction::calculatePorBar( model,
gridPart,
intersections,
values,
position,
elementSet,
m_defaultPorePressureGradient );
if ( pos.isUndefined() )
{
auto cellIdx = m_mainGrid->findReservoirCellIndexFromPoint( position );
if ( cellIdx != cvf::UNDEFINED_SIZE_T )
{
double valueFromEclipse = m_resultAccessor->cellScalar( cellIdx );
if ( !std::isinf( valueFromEclipse ) ) return RiaEclipseUnitTools::barToPascal( valueFromEclipse );
}
}

return RiaEclipseUnitTools::barToPascal( value );
}
Expand Down
23 changes: 15 additions & 8 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,19 +83,22 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d

cvf::Vec3d topPosition( position.x(), position.y(), topDepth );
cvf::Vec3d bottomPosition( position.x(), position.y(), bottomDepth );
auto part = model.grid( gridPart );

auto [isOk, elementSet] = findElementSetForElementIndex( part->elementSets(), static_cast<int>( elementIndex ) );

if ( isPositionValid( position, topPosition, bottomPosition, gridPart ) )
{
if ( m_property == RimFaultReactivation::Property::StressTop )
{
auto [porBar, extractionPos] = calculatePorBar( topPosition, m_gradient, gridPart );
auto [porBar, extractionPos] = calculatePorBar( topPosition, elementSet, m_gradient, gridPart );
if ( std::isinf( porBar ) ) return porBar;
double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
return -RiaEclipseUnitTools::barToPascal( s33 - porBar );
}
else if ( m_property == RimFaultReactivation::Property::StressBottom )
{
auto [porBar, extractionPos] = calculatePorBar( bottomPosition, m_gradient, gridPart );
auto [porBar, extractionPos] = calculatePorBar( bottomPosition, elementSet, m_gradient, gridPart );
if ( std::isinf( porBar ) ) return porBar;
double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
return -RiaEclipseUnitTools::barToPascal( s33 - porBar );
Expand All @@ -110,11 +113,11 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
}
else if ( m_property == RimFaultReactivation::Property::LateralStressComponentX )
{
return lateralStressComponentX( position, gridPart );
return lateralStressComponentX( position, elementSet, gridPart );
}
else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY )
{
return lateralStressComponentY( position, gridPart );
return lateralStressComponentY( position, elementSet, gridPart );
}
}

Expand All @@ -124,9 +127,11 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorStress::lateralStressComponentX( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const
double RimFaultReactivationDataAccessorStress::lateralStressComponentX( const cvf::Vec3d& position,
RimFaultReactivation::ElementSets elementSet,
RimFaultReactivation::GridPart gridPart ) const
{
auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, gridPart );
auto [porBar, extractionPos] = calculatePorBar( position, elementSet, m_gradient, gridPart );
if ( std::isinf( porBar ) ) return porBar;
double s11 = extractStressValue( StressType::S11, extractionPos, gridPart );
double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
Expand All @@ -136,9 +141,11 @@ double RimFaultReactivationDataAccessorStress::lateralStressComponentX( const cv
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorStress::lateralStressComponentY( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const
double RimFaultReactivationDataAccessorStress::lateralStressComponentY( const cvf::Vec3d& position,
RimFaultReactivation::ElementSets elementSet,
RimFaultReactivation::GridPart gridPart ) const
{
auto [porBar, extractionPos] = calculatePorBar( position, m_gradient, gridPart );
auto [porBar, extractionPos] = calculatePorBar( position, elementSet, m_gradient, gridPart );
if ( std::isinf( porBar ) ) return porBar;
double s22 = extractStressValue( StressType::S22, extractionPos, gridPart );
double s33 = extractStressValue( StressType::S33, extractionPos, gridPart );
Expand Down
14 changes: 10 additions & 4 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,16 +63,22 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc

virtual double extractStressValue( StressType stressType, const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const = 0;

virtual std::pair<double, cvf::Vec3d>
calculatePorBar( const cvf::Vec3d& position, double gradient, RimFaultReactivation::GridPart gridPart ) const = 0;
virtual std::pair<double, cvf::Vec3d> calculatePorBar( const cvf::Vec3d& position,
RimFaultReactivation::ElementSets elementSet,
double gradient,
RimFaultReactivation::GridPart gridPart ) const = 0;

virtual bool isPositionValid( const cvf::Vec3d& position,
const cvf::Vec3d& topPosition,
const cvf::Vec3d& bottomPosition,
RimFaultReactivation::GridPart gridPart ) const = 0;

virtual double lateralStressComponentX( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const;
virtual double lateralStressComponentY( const cvf::Vec3d& position, RimFaultReactivation::GridPart gridPart ) const;
virtual double lateralStressComponentX( const cvf::Vec3d& position,
RimFaultReactivation::ElementSets elementSet,
RimFaultReactivation::GridPart gridPart ) const;
virtual double lateralStressComponentY( const cvf::Vec3d& position,
RimFaultReactivation::ElementSets elementSet,
RimFaultReactivation::GridPart gridPart ) const;

RimFaultReactivation::Property m_property;
double m_gradient;
Expand Down
Loading

0 comments on commit 5a3894b

Please sign in to comment.