mirror of
https://github.com/OPM/ResInsight.git
synced 2026-10-26 15:10:22 -05:00
Kristian Bendiksen
6357cef1ab
The sign of the projected formation direction vector flips with the fracture direction normal, which follows the well drilling direction, while the fracture azimuth defining the dip rotation in RimFracture::transformMatrix() is only defined mod 180 degrees. Canonicalize the formation direction against the fracture azimuth direction before computing the dip, so a positive dip means the formation descends toward the azimuth direction. Restore std::abs() in calculateFormationDipFromHorizontal() so the StimPlan model formation dip (BedDipDeg in the Asymmetric FRK export) keeps the non-negative convention from #8877.
323 lines
15 KiB
C++
323 lines
15 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (C) 2021- 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 "RigStimPlanModelTools.h"
|
|
|
|
#include "RiaLogging.h"
|
|
#include "RiaQStringFormatter.h"
|
|
|
|
#include "RigCell.h"
|
|
#include "RigEclipseCaseData.h"
|
|
#include "RigFault.h"
|
|
#include "RigMainGrid.h"
|
|
#include "Well/RigSimulationWellCoordsAndMD.h"
|
|
#include "Well/RigWellLogExtractor.h"
|
|
#include "Well/RigWellPathIntersectionTools.h"
|
|
|
|
#include "cvfBoundingBox.h"
|
|
#include "cvfGeometryTools.h"
|
|
#include "cvfMath.h"
|
|
#include "cvfPlane.h"
|
|
|
|
#include <QString>
|
|
|
|
#include <cmath>
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::Vec3d RigStimPlanModelTools::calculateTSTDirection( RigEclipseCaseData* eclipseCaseData,
|
|
const cvf::Vec3d& anchorPosition,
|
|
double boundingBoxHorizontal,
|
|
double boundingBoxVertical )
|
|
{
|
|
cvf::Vec3d defaultDirection = cvf::Vec3d( 0.0, 0.0, -1.0 );
|
|
|
|
if ( !eclipseCaseData ) return defaultDirection;
|
|
|
|
RigMainGrid* mainGrid = eclipseCaseData->mainGrid();
|
|
if ( !mainGrid ) return defaultDirection;
|
|
|
|
cvf::Vec3d boundingBoxSize( boundingBoxHorizontal, boundingBoxHorizontal, boundingBoxVertical );
|
|
|
|
// Find upper face of cells close to the anchor point
|
|
cvf::BoundingBox boundingBox( anchorPosition - boundingBoxSize, anchorPosition + boundingBoxSize );
|
|
std::vector<size_t> closeCells = mainGrid->findIntersectingCells( boundingBox );
|
|
|
|
// The stratigraphic thickness is the average of normals of the top face
|
|
cvf::Vec3d direction = cvf::Vec3d::ZERO;
|
|
|
|
int numContributingCells = 0;
|
|
for ( size_t globalCellIndex : closeCells )
|
|
{
|
|
const RigCell& cell = mainGrid->cell( globalCellIndex );
|
|
|
|
if ( !cell.isInvalid() )
|
|
{
|
|
direction += cell.faceNormalWithAreaLength( cvf::StructGridInterface::NEG_K ).getNormalized();
|
|
numContributingCells++;
|
|
}
|
|
}
|
|
|
|
RiaLogging::info( std::format( "TST contributing cells: {}/{}", numContributingCells, closeCells.size() ) );
|
|
if ( numContributingCells == 0 )
|
|
{
|
|
// No valid close cells found: just point straight up
|
|
return defaultDirection;
|
|
}
|
|
|
|
direction = ( direction / static_cast<double>( numContributingCells ) ).getNormalized();
|
|
|
|
// A surface has normals in both directions: if the normal points upwards we flip it to
|
|
// make it point downwards. This necessary when finding the TST start and end points later.
|
|
if ( direction.z() > 0.0 )
|
|
{
|
|
direction *= -1.0;
|
|
}
|
|
|
|
return direction;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigStimPlanModelTools::calculateAngleFromVertical( const cvf::Vec3d& direction )
|
|
{
|
|
// Angle from vertical is the inclination of a vector from the vertical axis.
|
|
return cvf::Math::toDegrees( cvf::GeometryTools::getAngle( direction, -cvf::Vec3d::Z_AXIS ) );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigStimPlanModelTools::calculateFormationDipFromHorizontal( const cvf::Vec3d& direction )
|
|
{
|
|
// Convert from angle-with-vertical to angle-from-horizontal, and take abs since
|
|
// direction sign depends on fracture normal orientation.
|
|
return std::abs( calculateAngleFromVertical( direction ) - 90.0 );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigStimPlanModelTools::calculateFormationDipAlignedToAzimuth( const cvf::Vec3d& formationDirection, double azimuthDegrees )
|
|
{
|
|
// The orientation of the formation direction vector depends on the fracture normal, which can
|
|
// point either way along the well. Canonicalize it to point along the fracture azimuth direction
|
|
// before computing the dip: a positive result then means the formation descends toward the
|
|
// azimuth direction, matching the sign convention of the fracture dip rotation in
|
|
// RimFracture::transformMatrix().
|
|
double azimuthRadians = cvf::Math::toRadians( azimuthDegrees );
|
|
cvf::Vec3d azimuthDirection( std::sin( azimuthRadians ), std::cos( azimuthRadians ), 0.0 );
|
|
|
|
cvf::Vec3d direction = formationDirection;
|
|
if ( direction * azimuthDirection < 0.0 ) direction = -direction;
|
|
|
|
return 90.0 - calculateAngleFromVertical( direction );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::tuple<const RigFault*, double, cvf::Vec3d, double> RigStimPlanModelTools::findClosestFaultBarrier( RigEclipseCaseData* eclipseCaseData,
|
|
const cvf::Vec3d& position,
|
|
const cvf::Vec3d& directionToBarrier )
|
|
{
|
|
std::vector<WellPathCellIntersectionInfo> intersections =
|
|
RigStimPlanModelTools::generateBarrierIntersections( eclipseCaseData, position, directionToBarrier );
|
|
|
|
RiaLogging::info( std::format( "Intersections: {}", intersections.size() ) );
|
|
|
|
double shortestDistance = std::numeric_limits<double>::max();
|
|
|
|
RigMainGrid* mainGrid = eclipseCaseData->mainGrid();
|
|
cvf::Vec3d barrierPosition;
|
|
double barrierDip = 0.0;
|
|
const RigFault* foundFault = nullptr;
|
|
for ( const WellPathCellIntersectionInfo& intersection : intersections )
|
|
{
|
|
// Find the closest cell face which is a fault
|
|
double distance = position.pointDistance( intersection.startPoint );
|
|
const RigFault* fault = mainGrid->findFaultFromCellIndexAndCellFace( intersection.globCellIndex, intersection.intersectedCellFaceIn );
|
|
if ( fault && distance < shortestDistance )
|
|
{
|
|
foundFault = fault;
|
|
shortestDistance = distance;
|
|
barrierPosition = intersection.startPoint;
|
|
|
|
const RigCell& cell = mainGrid->cell( intersection.globCellIndex );
|
|
cvf::Vec3d faceNormal = cell.faceNormalWithAreaLength( intersection.intersectedCellFaceIn );
|
|
barrierDip = RigStimPlanModelTools::calculateAngleFromVertical( faceNormal );
|
|
}
|
|
}
|
|
|
|
return std::make_tuple( foundFault, shortestDistance, barrierPosition, barrierDip );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<WellPathCellIntersectionInfo> RigStimPlanModelTools::generateBarrierIntersections( RigEclipseCaseData* eclipseCaseData,
|
|
const cvf::Vec3d& position,
|
|
const cvf::Vec3d& directionToBarrier )
|
|
{
|
|
double randoDistance = 10000.0;
|
|
cvf::Vec3d forwardPosition = position + ( directionToBarrier * randoDistance );
|
|
cvf::Vec3d backwardPosition = position + ( directionToBarrier * -randoDistance );
|
|
std::vector<WellPathCellIntersectionInfo> intersections =
|
|
generateBarrierIntersectionsBetweenPoints( eclipseCaseData, position, forwardPosition );
|
|
std::vector<WellPathCellIntersectionInfo> backwardIntersections =
|
|
generateBarrierIntersectionsBetweenPoints( eclipseCaseData, position, backwardPosition );
|
|
|
|
// Merge the intersections for the search for closest
|
|
intersections.insert( intersections.end(), backwardIntersections.begin(), backwardIntersections.end() );
|
|
return intersections;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<WellPathCellIntersectionInfo> RigStimPlanModelTools::generateBarrierIntersectionsBetweenPoints( RigEclipseCaseData* eclipseCaseData,
|
|
const cvf::Vec3d& startPosition,
|
|
const cvf::Vec3d& endPosition )
|
|
{
|
|
// Create a fake well path from the anchor point to
|
|
// a point far away in the direction barrier direction
|
|
std::vector<cvf::Vec3d> pathCoords;
|
|
pathCoords.push_back( startPosition );
|
|
pathCoords.push_back( endPosition );
|
|
|
|
RigSimulationWellCoordsAndMD helper( pathCoords );
|
|
return RigWellPathIntersectionTools::findCellIntersectionInfosAlongPath( eclipseCaseData,
|
|
"",
|
|
helper.wellPathPoints(),
|
|
helper.measuredDepths() );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigStimPlanModelTools::findThicknessTargetPoints( RigEclipseCaseData* eclipseCaseData,
|
|
const cvf::Vec3d& position,
|
|
const cvf::Vec3d& direction,
|
|
double extractionDepthTop,
|
|
double extractionDepthBottom,
|
|
cvf::Vec3d& topPosition,
|
|
cvf::Vec3d& bottomPosition )
|
|
{
|
|
if ( !eclipseCaseData ) return false;
|
|
|
|
RiaLogging::info( std::format( "Position: {}", RigStimPlanModelTools::vecToString( position ) ) );
|
|
RiaLogging::info( std::format( "Direction: {}", RigStimPlanModelTools::vecToString( direction ) ) );
|
|
|
|
// Create a "fake" well path which from top to bottom of formation
|
|
// passing through the point and with the given direction
|
|
|
|
const cvf::BoundingBox& allCellsBoundingBox = eclipseCaseData->mainGrid()->boundingBox();
|
|
|
|
RiaLogging::info( QString( "All cells bounding box: %1 %2" )
|
|
.arg( RigStimPlanModelTools::vecToString( allCellsBoundingBox.min() ) )
|
|
.arg( RigStimPlanModelTools::vecToString( allCellsBoundingBox.max() ) )
|
|
.toStdString() );
|
|
cvf::BoundingBox geometryBoundingBox( allCellsBoundingBox );
|
|
|
|
// Use smaller depth bounding box for extraction if configured
|
|
if ( extractionDepthTop > 0.0 && extractionDepthBottom > 0.0 && extractionDepthTop < extractionDepthBottom )
|
|
{
|
|
cvf::Vec3d bbMin( allCellsBoundingBox.min().x(), allCellsBoundingBox.min().y(), -extractionDepthBottom );
|
|
cvf::Vec3d bbMax( allCellsBoundingBox.max().x(), allCellsBoundingBox.max().y(), -extractionDepthTop );
|
|
geometryBoundingBox = cvf::BoundingBox( bbMin, bbMax );
|
|
}
|
|
|
|
if ( !geometryBoundingBox.contains( position ) )
|
|
{
|
|
//
|
|
RiaLogging::error( "Anchor position is outside the grid bounding box. Unable to compute direction." );
|
|
return false;
|
|
}
|
|
|
|
// Create plane on top and bottom of formation
|
|
cvf::Plane topPlane;
|
|
topPlane.setFromPointAndNormal( geometryBoundingBox.max(), cvf::Vec3d::Z_AXIS );
|
|
|
|
cvf::Plane bottomPlane;
|
|
bottomPlane.setFromPointAndNormal( geometryBoundingBox.min(), cvf::Vec3d::Z_AXIS );
|
|
|
|
// Find and add point on top plane
|
|
cvf::Vec3d abovePlane = position + ( direction * -10000.0 );
|
|
if ( !topPlane.intersect( position, abovePlane, &topPosition ) )
|
|
{
|
|
RiaLogging::error( "Unable to compute top position of thickness direction vector." );
|
|
return false;
|
|
}
|
|
|
|
RiaLogging::info( std::format( "Top: {}", RigStimPlanModelTools::vecToString( topPosition ) ) );
|
|
|
|
// Find and add point on bottom plane
|
|
cvf::Vec3d belowPlane = position + ( direction * 10000.0 );
|
|
if ( !bottomPlane.intersect( position, belowPlane, &bottomPosition ) )
|
|
{
|
|
RiaLogging::error( "Unable to compute bottom position of thickness direction vector." );
|
|
return false;
|
|
}
|
|
|
|
RiaLogging::info( std::format( "Bottom: {}", RigStimPlanModelTools::vecToString( bottomPosition ) ) );
|
|
|
|
return true;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
QString RigStimPlanModelTools::vecToString( const cvf::Vec3d& vec )
|
|
{
|
|
return QString( "[%1, %2, %3]" ).arg( vec.x() ).arg( vec.y() ).arg( vec.z() );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigStimPlanModelTools::calculatePerforationLength( const cvf::Vec3d& direction, double perforationLength )
|
|
{
|
|
// Deviation from vertical. Since well path is tending downwards we compare with negative z.
|
|
double inclination = cvf::GeometryTools::getAngle( direction, -cvf::Vec3d::Z_AXIS );
|
|
|
|
// Keep inclination in 0-90 degrees range
|
|
if ( inclination > cvf::PI_D / 2.0 )
|
|
{
|
|
inclination = cvf::PI_D - inclination;
|
|
}
|
|
|
|
double correctedPerforationLength = perforationLength * std::cos( inclination );
|
|
|
|
RiaLogging::info( QString( "Perforation length correction: original length: %1 inclination: %2 corrected length: %3" )
|
|
.arg( perforationLength )
|
|
.arg( cvf::Math::toDegrees( inclination ) )
|
|
.arg( correctedPerforationLength )
|
|
.toStdString() );
|
|
|
|
// Handle well inclination close to 90 dgr to ensure visual perforation interval in StimPlan model plot
|
|
if ( std::fabs( cvf::Math::toDegrees( inclination ) - 90.0 ) < 0.1 )
|
|
{
|
|
double minimumPerforationInterval = 0.5;
|
|
return std::max( minimumPerforationInterval, correctedPerforationLength );
|
|
}
|
|
|
|
return correctedPerforationLength;
|
|
}
|