Skip to content

Commit

Permalink
WIP: try to extract porbar
Browse files Browse the repository at this point in the history
  • Loading branch information
@kriben
kriben committed Dec 8, 2023
1 parent 7eb528a commit f16c60f
Show file tree
Hide file tree
Showing 12 changed files with 154 additions and 26 deletions.
2 changes: 1 addition & 1 deletion ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccess.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,7 @@ std::vector<double> RimFaultReactivationDataAccess::extractModelData( const RigF
double bottomDepth = computeAverageDepth( corners, { 4, 5, 6, 7 } );

cvf::Vec3d position = RigCaseToCaseCellMapperTools::calculateCellCenter( corners.data() );
double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth );
double value = accessor->valueAtPosition( position, model, gridPart, topDepth, bottomDepth, elementIndex );
values.push_back( value );
}
}
Expand Down
5 changes: 3 additions & 2 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessor.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,9 @@ class RimFaultReactivationDataAccessor
virtual double valueAtPosition( const cvf::Vec3d& position,
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity() ) const = 0;
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity(),
size_t elementIndex = std::numeric_limits<size_t>::max() ) const = 0;

protected:
virtual void updateResultAccessor() = 0;
Expand Down
4 changes: 2 additions & 2 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,8 +101,8 @@ double RimFaultReactivationDataAccessorGeoMech::valueAtPosition( const cvf::Vec3
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth,
double bottomDepth ) const

double bottomDepth,
size_t elementIndex ) const
{
if ( !m_resultFrames ) return std::numeric_limits<double>::infinity();

Expand Down
5 changes: 3 additions & 2 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorGeoMech.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,8 +43,9 @@ class RimFaultReactivationDataAccessorGeoMech : public RimFaultReactivationDataA
double valueAtPosition( const cvf::Vec3d& position,
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity() ) const override;
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity(),
size_t elementIndex = std::numeric_limits<size_t>::max() ) const override;

private:
void updateResultAccessor() override;
Expand Down
3 changes: 2 additions & 1 deletion ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,8 @@ double RimFaultReactivationDataAccessorPorePressure::valueAtPosition( const cvf:
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth,
double bottomDepth ) const
double bottomDepth,
size_t elementIndex ) const
{
if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
{
Expand Down
5 changes: 3 additions & 2 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorPorePressure.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorPorePressure : public RimFaultReactivation
double valueAtPosition( const cvf::Vec3d& position,
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity() ) const override;
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity(),
size_t elementIndex = std::numeric_limits<size_t>::max() ) const override;

private:
void updateResultAccessor() override;
Expand Down
117 changes: 112 additions & 5 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,15 +20,18 @@

#include "RiaEclipseUnitTools.h"

#include "RigFaultReactivationModel.h"
#include "RigFemAddressDefines.h"
#include "RigFemPartCollection.h"
#include "RigFemPartResultsCollection.h"
#include "RigFemResultAddress.h"
#include "RigFemScalarResultFrames.h"
#include "RigGeoMechCaseData.h"
#include "RigGriddedPart3d.h"
#include "RigResultAccessorFactory.h"

#include "RimFaultReactivationEnums.h"
#include "RimFracture.h"
#include "RimGeoMechCase.h"
#include "RimWellIADataAccess.h"

Expand Down Expand Up @@ -101,7 +104,8 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth,
double bottomDepth ) const
double bottomDepth,
size_t elementIndex ) const
{
if ( !m_porFrames || !m_s11Frames || !m_s22Frames || !m_s33Frames || !m_femPart ) return std::numeric_limits<double>::infinity();

Expand All @@ -123,16 +127,18 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
const std::vector<float>& s33Data = m_s33Frames->frameData( timeStepIndex, frameIndex );
const std::vector<float>& porData = m_porFrames->frameData( timeStepIndex, frameIndex );

const RigGriddedPart3d* grid = model.grid( gridPart );

if ( m_property == RimFaultReactivation::Property::StressTop )
{
double s33 = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, s33Data );
double porBar = interpolatedResultValue( iaDataAccess, m_femPart, topPosition, porData );
double porBar = getPorBar( iaDataAccess, m_femPart, topPosition, porData, *grid, elementIndex );
return RiaEclipseUnitTools::barToPascal( s33 - porBar );
}
else if ( m_property == RimFaultReactivation::Property::StressBottom )
{
double s33 = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, s33Data );
double porBar = interpolatedResultValue( iaDataAccess, m_femPart, bottomPosition, porData );
double porBar = getPorBar( iaDataAccess, m_femPart, bottomPosition, porData, *grid, elementIndex );
return RiaEclipseUnitTools::barToPascal( s33 - porBar );
}
else if ( m_property == RimFaultReactivation::Property::DepthTop )
Expand All @@ -147,14 +153,14 @@ double RimFaultReactivationDataAccessorStress::valueAtPosition( const cvf::Vec3d
{
double s11 = interpolatedResultValue( iaDataAccess, m_femPart, position, s11Data );
double s33 = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data );
double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData );
double porBar = getPorBar( iaDataAccess, m_femPart, position, porData, *grid, elementIndex );
return ( s11 - porBar ) / ( s33 - porBar );
}
else if ( m_property == RimFaultReactivation::Property::LateralStressComponentY )
{
double s22 = interpolatedResultValue( iaDataAccess, m_femPart, position, s22Data );
double s33 = interpolatedResultValue( iaDataAccess, m_femPart, position, s33Data );
double porBar = interpolatedResultValue( iaDataAccess, m_femPart, position, porData );
double porBar = getPorBar( iaDataAccess, m_femPart, position, porData, *grid, elementIndex );
return ( s22 - porBar ) / ( s33 - porBar );
}
}
Expand All @@ -172,3 +178,104 @@ double RimFaultReactivationDataAccessorStress::interpolatedResultValue( RimWellI
{
return iaDataAccess.interpolatedResultValue( femPart, scalarResults, RIG_ELEMENT_NODAL, position );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFaultReactivationDataAccessorStress::getPorBar( RimWellIADataAccess& iaDataAccess,
const RigFemPart* femPart,
const cvf::Vec3d& position,
const std::vector<float>& porData,
const RigGriddedPart3d& grid,
size_t elementIndex ) const
{
// Find element set for the position
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets = grid.elementSets();
auto [isOk, elementSet] = findElementSetContainingElement( elementSets, elementIndex );

if ( !isOk )
{
printf( "No data for pos: [%f, %f, %f] %zu.\n", position.x(), position.y(), position.z(), elementIndex );
return std::numeric_limits<double>::infinity();
}

//
std::vector<std::pair<double, double>> layers = grid.layers( elementSet );

if ( elementSet == RimFaultReactivation::ElementSets::OverBurden )
{
// Find top reservoir pore pressure
std::vector<std::pair<double, double>> reservoirLayers = grid.layers( RimFaultReactivation::ElementSets::Reservoir );
if ( reservoirLayers.empty() )
{
printf( "MISSING RESERVOIR LAYER FOR OVERBURDEN POS: %f\n", position.z() );
return std::numeric_limits<double>::infinity();
}

cvf::Vec3d reservoirTopPos( position.x(), position.y(), reservoirLayers[0].first );
double value = interpolatedResultValue( iaDataAccess, femPart, reservoirTopPos, porData );
printf( "Overburden: depth: %f. data depth: %f (distance: %f) Value: %g\n",
position.z(),
reservoirTopPos.z(),
position.z() - reservoirTopPos.z(),
value );
return value;
}
else if ( elementSet == RimFaultReactivation::ElementSets::Reservoir )
{
double value = interpolatedResultValue( iaDataAccess, femPart, position, porData );
printf( "Reservoir: depth: %f. Value: %g\n", position.z(), value );
return value;
}
else if ( elementSet == RimFaultReactivation::ElementSets::IntraReservoir )
{
double value = interpolatedResultValue( iaDataAccess, femPart, position, porData );
printf( "Intra-reservoir: depth: %f. Value: %g\n", position.z(), value );

return value;
}
else if ( elementSet == RimFaultReactivation::ElementSets::UnderBurden )
{
// Find top reservoir pore pressure
std::vector<std::pair<double, double>> reservoirLayers = grid.layers( RimFaultReactivation::ElementSets::Reservoir );
if ( reservoirLayers.empty() )
{
printf( "MISSING RESERVOIR LAYER FOR UNDERBURDEN POS: %f\n", position.z() );
return std::numeric_limits<double>::infinity();
}

cvf::Vec3d reservoirPos( position.x(), position.y(), reservoirLayers.back().second );
double value = interpolatedResultValue( iaDataAccess, femPart, reservoirPos, porData );
printf( "Underburden: depth: %f. data depth: %f (distance: %f) Value: %g\n",
position.z(),
reservoirPos.z(),
position.z() - reservoirPos.z(),
value );
return value;
}

return interpolatedResultValue( iaDataAccess, m_femPart, position, porData );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<bool, RimFaultReactivation::ElementSets> RimFaultReactivationDataAccessorStress::findElementSetContainingElement(
const std::map<RimFaultReactivation::ElementSets, std::vector<unsigned int>>& elementSets,
unsigned int elmIdx )
{
// Iterate through the map to find the ElementSets for the given unsigned int
for ( auto [element, elementSet] : elementSets )
{
// Check if the targetValue is in the current vector
auto it = std::find( elementSet.begin(), elementSet.end(), elmIdx );

// If found, print the ElementSets
if ( it != elementSet.end() )
{
return { true, element };
}
}

return { false, RimFaultReactivation::ElementSets::OverBurden };
}
23 changes: 18 additions & 5 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorStress.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,12 @@

#include <vector>

class RigFemPart;
class RigFemPartResultsCollection;
class RimGeoMechCase;
class RigGeoMechCaseData;
class RigFemScalarResultFrames;
class RigFemPart;
class RigGeoMechCaseData;
class RigGriddedPart3d;
class RimGeoMechCase;
class RimWellIADataAccess;

//==================================================================================================
Expand All @@ -47,8 +48,9 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc
double valueAtPosition( const cvf::Vec3d& position,
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity() ) const override;
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity(),
size_t elementIndex = std::numeric_limits<size_t>::max() ) const override;

private:
void updateResultAccessor() override;
Expand All @@ -60,6 +62,17 @@ class RimFaultReactivationDataAccessorStress : public RimFaultReactivationDataAc
const cvf::Vec3d& position,
const std::vector<float>& scalarResults ) const;

double getPorBar( RimWellIADataAccess& iaDataAccess,
const RigFemPart* femPart,
const cvf::Vec3d& position,
const std::vector<float>& porData,
const RigGriddedPart3d& grid,
size_t elementIndex ) const;

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

RimGeoMechCase* m_geoMechCase;
RimFaultReactivation::Property m_property;
RigGeoMechCaseData* m_geoMechCaseData;
Expand Down
3 changes: 2 additions & 1 deletion ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,8 @@ double RimFaultReactivationDataAccessorTemperature::valueAtPosition( const cvf::
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth,
double bottomDepth ) const
double bottomDepth,
size_t elementIndex ) const
{
if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
{
Expand Down
5 changes: 3 additions & 2 deletions ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorTemperature.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,8 +44,9 @@ class RimFaultReactivationDataAccessorTemperature : public RimFaultReactivationD
double valueAtPosition( const cvf::Vec3d& position,
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity() ) const override;
double topDepth = std::numeric_limits<double>::infinity(),
double bottomDepth = std::numeric_limits<double>::infinity(),
size_t elementIndex = std::numeric_limits<size_t>::max() ) const override;

private:
void updateResultAccessor() override;
Expand Down
3 changes: 2 additions & 1 deletion ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationDataAccessorVoidRatio.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,8 @@ double RimFaultReactivationDataAccessorVoidRatio::valueAtPosition( const cvf::Ve
const RigFaultReactivationModel& model,
RimFaultReactivation::GridPart gridPart,
double topDepth,
double bottomDepth ) const
double bottomDepth,
size_t elementIndex ) const
{
if ( ( m_mainGrid != nullptr ) && m_resultAccessor.notNull() )
{
Expand Down
Loading

0 comments on commit f16c60f

Please sign in to comment.