From 8eda5bcbb979a7b6f449efaf391bf5c85f31781b Mon Sep 17 00:00:00 2001 From: Gaute Lindkvist Date: Thu, 6 Sep 2018 19:59:47 +0200 Subject: [PATCH] #3313 Implement weighted average of WBS data. --- .../RigGeoMechWellLogExtractor.cpp | 138 +++++++++++------- .../RigGeoMechWellLogExtractor.h | 6 +- Fwk/AppFwk/cafTensor/cafTensor3.h | 2 +- Fwk/AppFwk/cafTensor/cafTensor3.inl | 13 ++ 4 files changed, 100 insertions(+), 59 deletions(-) diff --git a/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp b/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp index c3d06573b7..9df849155b 100644 --- a/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp +++ b/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.cpp @@ -23,6 +23,7 @@ #include "RigGeoMechWellLogExtractor.h" #include "RiaDefines.h" +#include "RiaWeightedAverageCalculator.h" #include "RigFemTypes.h" #include "RigGeoMechBoreHoleStressCalculator.h" #include "RigFemPart.h" @@ -108,18 +109,18 @@ void RigGeoMechWellLogExtractor::curveData(const RigFemResultAddress& resAddr, i //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -float RigGeoMechWellLogExtractor::calculatePorePressureInSegment(int64_t intersectionIdx, float averageSegmentPorePressureBars, double hydroStaticPorePressureBars, double effectiveDepthMeters, const std::vector& poreElementPressuresPascal) const +float RigGeoMechWellLogExtractor::calculatePorePressureInSegment(int64_t intersectionIdx, float averageSegmentPorePressureBar, double hydroStaticPorePressureBar, double effectiveDepthMeters, const std::vector& poreElementPressuresPascal) const { - double porePressure = hydroStaticPorePressureBars; + double porePressure = hydroStaticPorePressureBar; // 1: Try pore pressure from the grid - if (porePressure == hydroStaticPorePressureBars && averageSegmentPorePressureBars != std::numeric_limits::infinity()) + if (porePressure == hydroStaticPorePressureBar && averageSegmentPorePressureBar != std::numeric_limits::infinity()) { - porePressure = averageSegmentPorePressureBars; + porePressure = averageSegmentPorePressureBar; } // 2: Try mud weight from LAS-file to generate pore pressure - if (porePressure == hydroStaticPorePressureBars && !m_wellLogMdAndMudWeightKgPerM3.empty()) + if (porePressure == hydroStaticPorePressureBar && !m_wellLogMdAndMudWeightKgPerM3.empty()) { double lasMudWeightKgPerM3 = getWellLogSegmentValue(intersectionIdx, m_wellLogMdAndMudWeightKgPerM3); if (lasMudWeightKgPerM3 != std::numeric_limits::infinity()) @@ -131,7 +132,7 @@ float RigGeoMechWellLogExtractor::calculatePorePressureInSegment(int64_t interse } size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx]; // 3: Try pore pressure from element property tables - if (porePressure == hydroStaticPorePressureBars && elmIdx < poreElementPressuresPascal.size()) + if (porePressure == hydroStaticPorePressureBar && elmIdx < poreElementPressuresPascal.size()) { // Pore pressure from element property tables are in pascal. porePressure = pascalToBar(poreElementPressuresPascal[elmIdx]); @@ -174,25 +175,25 @@ float RigGeoMechWellLogExtractor::calculateUcs(int64_t intersectionIdx, const st { // Typical UCS: http://ceae.colorado.edu/~amadei/CVEN5768/PDF/NOTES8.pdf // Typical UCS for Shale is 5 - 100 MPa -> 50 - 1000 bar. - const double defaultUniaxialStrengthInBars = 100.0; + const double defaultUniaxialStrengthInBar = 100.0; - double uniaxialStrengthInBars = defaultUniaxialStrengthInBars; + double uniaxialStrengthInBar = defaultUniaxialStrengthInBar; if (!m_wellLogMdAndUcsBar.empty()) { - double lasUniaxialStrengthInBars = getWellLogSegmentValue(intersectionIdx, m_wellLogMdAndUcsBar); - if (lasUniaxialStrengthInBars != std::numeric_limits::infinity()) + double lasUniaxialStrengthInBar = getWellLogSegmentValue(intersectionIdx, m_wellLogMdAndUcsBar); + if (lasUniaxialStrengthInBar != std::numeric_limits::infinity()) { - uniaxialStrengthInBars = lasUniaxialStrengthInBars; + uniaxialStrengthInBar = lasUniaxialStrengthInBar; } } size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx]; - if (uniaxialStrengthInBars == defaultUniaxialStrengthInBars && elmIdx < ucsValuesPascal.size()) + if (uniaxialStrengthInBar == defaultUniaxialStrengthInBar && elmIdx < ucsValuesPascal.size()) { // Read UCS from element table in Pascal - uniaxialStrengthInBars = pascalToBar(ucsValuesPascal[elmIdx]); + uniaxialStrengthInBar = pascalToBar(ucsValuesPascal[elmIdx]); } - return uniaxialStrengthInBars; + return uniaxialStrengthInBar; } @@ -310,13 +311,12 @@ void RigGeoMechWellLogExtractor::wellPathScaledCurveData(const RigFemResultAddre if (!(elmType == HEX8 || elmType == HEX8P)) continue; - const int* elmNodeIndices = femPart->connectivities(elmIdx); - cvf::Vec3f centroid = cellCentroid(elmNodeIndices, nodeCoords); + cvf::Vec3f centroid = cellCentroid(intersectionIdx); double trueVerticalDepth = -centroid.z(); double effectiveDepthMeters = trueVerticalDepth + m_rkbDiff; - double hydroStaticPorePressureBars = pascalToBar(effectiveDepthMeters * UNIT_WEIGHT_OF_WATER); + double hydroStaticPorePressureBar = pascalToBar(effectiveDepthMeters * UNIT_WEIGHT_OF_WATER); float averageUnscaledValue = std::numeric_limits::infinity(); bool validAverage = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfaceValues, std::numeric_limits::infinity(), &averageUnscaledValue); @@ -324,11 +324,11 @@ void RigGeoMechWellLogExtractor::wellPathScaledCurveData(const RigFemResultAddre if (resAddr.fieldName == "PP") { double segmentPorePressureFromGrid = averageUnscaledValue; - averageUnscaledValue = calculatePorePressureInSegment(intersectionIdx, segmentPorePressureFromGrid, hydroStaticPorePressureBars, effectiveDepthMeters, poreElementPressuresPascal); + averageUnscaledValue = calculatePorePressureInSegment(intersectionIdx, segmentPorePressureFromGrid, hydroStaticPorePressureBar, effectiveDepthMeters, poreElementPressuresPascal); } - (*values)[intersectionIdx] = static_cast(averageUnscaledValue) / hydroStaticPorePressureBars; + (*values)[intersectionIdx] = static_cast(averageUnscaledValue) / hydroStaticPorePressureBar; } } @@ -366,11 +366,11 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData(const RigFemResultAddress std::vector poissonRatios = resultCollection->resultValues(poissonResAddr, 0, frameIndex); std::vector ucsValuesPascal = resultCollection->resultValues(ucsResAddr, 0, frameIndex); - std::vector interpolatedInterfacePP; - interpolatedInterfacePP.resize(m_intersections.size(), std::numeric_limits::infinity()); + std::vector interpolatedInterfacePorePressureBar; + interpolatedInterfacePorePressureBar.resize(m_intersections.size(), std::numeric_limits::infinity()); - std::vector interpolatedInterfaceStress; - interpolatedInterfaceStress.resize(m_intersections.size()); + std::vector interpolatedInterfaceStressBar; + interpolatedInterfaceStressBar.resize(m_intersections.size()); #pragma omp parallel for for (int64_t intersectionIdx = 0; intersectionIdx < (int64_t)m_intersections.size(); ++intersectionIdx) { @@ -378,8 +378,8 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData(const RigFemResultAddress RigElementType elmType = femPart->elementType(elmIdx); if (!(elmType == HEX8 || elmType == HEX8P)) continue; - interpolatedInterfacePP[intersectionIdx] = interpolateGridResultValue(porBarResAddr.resultPosType, porePressures, intersectionIdx, false); - interpolatedInterfaceStress[intersectionIdx] = interpolateGridResultValue(stressResAddr.resultPosType, vertexStresses, intersectionIdx, false); + interpolatedInterfacePorePressureBar[intersectionIdx] = interpolateGridResultValue(porBarResAddr.resultPosType, porePressures, intersectionIdx, false); + interpolatedInterfaceStressBar[intersectionIdx] = interpolateGridResultValue(stressResAddr.resultPosType, vertexStresses, intersectionIdx, false); } values->resize(m_intersections.size(), 0.0f); @@ -392,32 +392,31 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData(const RigFemResultAddress if (!(elmType == HEX8 || elmType == HEX8P)) continue; - const int* elmNodeIndices = femPart->connectivities(elmIdx); - cvf::Vec3f centroid = cellCentroid(elmNodeIndices, nodeCoords); + cvf::Vec3f centroid = cellCentroid(intersectionIdx); double trueVerticalDepth = -centroid.z(); double effectiveDepthMeters = trueVerticalDepth + m_rkbDiff; - double hydroStaticPorePressureBars = pascalToBar(effectiveDepthMeters * UNIT_WEIGHT_OF_WATER); + double hydroStaticPorePressureBar = pascalToBar(effectiveDepthMeters * UNIT_WEIGHT_OF_WATER); - float averageUnscaledPP = std::numeric_limits::infinity(); - bool validGridPP = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfacePP, std::numeric_limits::infinity(), &averageUnscaledPP); + float averagePorePressureBar = std::numeric_limits::infinity(); + bool validGridPorePressure = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfacePorePressureBar, std::numeric_limits::infinity(), &averagePorePressureBar); - double porePressureBars = calculatePorePressureInSegment(intersectionIdx, averageUnscaledPP, hydroStaticPorePressureBars, effectiveDepthMeters, poreElementPressuresPascal); - double poissonRatio = calculatePoissonRatio(intersectionIdx, poissonRatios); - double ucsBars = calculateUcs(intersectionIdx, ucsValuesPascal); + double porePressureBar = calculatePorePressureInSegment(intersectionIdx, averagePorePressureBar, hydroStaticPorePressureBar, effectiveDepthMeters, poreElementPressuresPascal); + double poissonRatio = calculatePoissonRatio(intersectionIdx, poissonRatios); + double ucsBar = calculateUcs(intersectionIdx, ucsValuesPascal); caf::Ten3d segmentStress; - bool validSegmentStress = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfaceStress, caf::Ten3d::invalid(), &segmentStress); + bool validSegmentStress = averageIntersectionValuesToSegmentValue(intersectionIdx, interpolatedInterfaceStressBar, caf::Ten3d::invalid(), &segmentStress); cvf::Vec3d wellPathTangent = calculateWellPathTangent(intersectionIdx, TangentConstantWithinCell); caf::Ten3d wellPathStressFloat = transformTensorToWellPathOrientation(wellPathTangent, segmentStress); caf::Ten3d wellPathStressDouble(wellPathStressFloat); - RigGeoMechBoreHoleStressCalculator sigmaCalculator(wellPathStressDouble, porePressureBars, poissonRatio, ucsBars, 32); + RigGeoMechBoreHoleStressCalculator sigmaCalculator(wellPathStressDouble, porePressureBar, poissonRatio, ucsBar, 32); double resultValue = std::numeric_limits::infinity(); if (resAddr.fieldName == RiaDefines::wellPathFGResultName().toStdString()) { - if (validGridPP) + if (validGridPorePressure) { resultValue = sigmaCalculator.solveFractureGradient(); } @@ -425,20 +424,16 @@ void RigGeoMechWellLogExtractor::wellBoreWallCurveData(const RigFemResultAddress else { CVF_ASSERT(resAddr.fieldName == RiaDefines::wellPathSFGResultName().toStdString()); - if (!validGridPP) + if (!validGridPorePressure) { resultValue = sigmaCalculator.solveStassiDalia(); } } if (resultValue != std::numeric_limits::infinity()) { - if (hydroStaticPorePressureBars > 1.0e-8) + if (hydroStaticPorePressureBar > 1.0e-8) { - resultValue /= hydroStaticPorePressureBars; - } - else - { - resultValue = std::numeric_limits::infinity(); + resultValue /= hydroStaticPorePressureBar; } } (*values)[intersectionIdx] = resultValue; @@ -750,14 +745,23 @@ caf::Ten3d RigGeoMechWellLogExtractor::transformTensorToWellPathOrientation(cons //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -cvf::Vec3f RigGeoMechWellLogExtractor::cellCentroid(const int* elmNodeIndices, const std::vector& nodeCoords) +cvf::Vec3f RigGeoMechWellLogExtractor::cellCentroid(size_t intersectionIdx) const { + const RigFemPart* femPart = m_caseData->femParts()->part(0); + const std::vector& nodeCoords = femPart->nodes().coordinates; + + size_t elmIdx = m_intersectedCellsGlobIdx[intersectionIdx]; + RigElementType elmType = femPart->elementType(elmIdx); + int elementNodeCount = RigFemTypes::elmentNodeCount(elmType); + + const int* elmNodeIndices = femPart->connectivities(elmIdx); + cvf::Vec3f centroid(0.0, 0.0, 0.0); - for (int i = 0; i < 8; ++i) + for (int i = 0; i < elementNodeCount; ++i) { centroid += nodeCoords[elmNodeIndices[i]]; } - return centroid / 8.0; + return centroid / elementNodeCount; } //-------------------------------------------------------------------------------------------------- @@ -765,6 +769,8 @@ cvf::Vec3f RigGeoMechWellLogExtractor::cellCentroid(const int* elmNodeIndices, c //-------------------------------------------------------------------------------------------------- double RigGeoMechWellLogExtractor::getWellLogSegmentValue(size_t intersectionIdx, const std::vector>& wellLogValues) const { + const RigFemPart* femPart = m_caseData->femParts()->part(0); + double startMD, endMD; if (intersectionIdx % 2 == 0) { @@ -777,20 +783,25 @@ double RigGeoMechWellLogExtractor::getWellLogSegmentValue(size_t intersectionIdx endMD = m_intersectionMeasuredDepths[intersectionIdx]; } - int nValuesWithinSegment = 0; - double sumValuesWithinSegment = 0.0; + RiaWeightedAverageCalculator averageCalc; for (auto& depthAndValue : wellLogValues) { if (cvf::Math::valueInRange(depthAndValue.first, startMD, endMD)) { - sumValuesWithinSegment += depthAndValue.second; - nValuesWithinSegment++; + cvf::Vec3d position = m_wellPath->interpolatedPointAlongWellPath(depthAndValue.first); + cvf::Vec3d centroid(cellCentroid(intersectionIdx)); + double weight = 1.0; + double dist = (position - centroid).length(); + if (dist > 1.0) + { + weight = 1.0 / dist; + } + averageCalc.addValueAndWeight(depthAndValue.second, weight); } } - if (nValuesWithinSegment) + if (averageCalc.validAggregatedWeight()) { - // TODO: Could do average weighted by inverse distance to center of element. - return sumValuesWithinSegment / nValuesWithinSegment; + return averageCalc.weightedAverage(); } return std::numeric_limits::infinity(); @@ -808,18 +819,29 @@ double RigGeoMechWellLogExtractor::pascalToBar(double pascalValue) /// //-------------------------------------------------------------------------------------------------- template -bool RigGeoMechWellLogExtractor::averageIntersectionValuesToSegmentValue(size_t intersectionIdx, const std::vector& values, const T& invalidValue, T* averagedCellValue) +bool RigGeoMechWellLogExtractor::averageIntersectionValuesToSegmentValue(size_t intersectionIdx, const std::vector& values, const T& invalidValue, T* averagedCellValue) const { CVF_ASSERT(values.size() >= 2); + + *averagedCellValue = invalidValue; + T value1, value2; + cvf::Vec3d centroid(cellCentroid(intersectionIdx)); + double dist1 = 0.0, dist2 = 0.0; if (intersectionIdx % 2 == 0) { value1 = values[intersectionIdx]; value2 = values[intersectionIdx + 1]; + + dist1 = (centroid - m_intersections[intersectionIdx]).length(); + dist2 = (centroid - m_intersections[intersectionIdx + 1]).length(); } else { value1 = values[intersectionIdx - 1]; value2 = values[intersectionIdx]; + + dist1 = (centroid - m_intersections[intersectionIdx - 1]).length(); + dist2 = (centroid - m_intersections[intersectionIdx]).length(); } if (invalidValue == value1 || invalidValue == value2) @@ -827,6 +849,12 @@ bool RigGeoMechWellLogExtractor::averageIntersectionValuesToSegmentValue(size_t return false; } - *averagedCellValue = (value1 + value2) * 0.5; + RiaWeightedAverageCalculator averageCalc; + averageCalc.addValueAndWeight(value1, dist2); + averageCalc.addValueAndWeight(value2, dist1); + if (averageCalc.validAggregatedWeight()) + { + *averagedCellValue = averageCalc.weightedAverage(); + } return true; } diff --git a/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h b/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h index 6cf6492a35..120341ebc2 100644 --- a/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h +++ b/ApplicationCode/ReservoirDataModel/RigGeoMechWellLogExtractor.h @@ -66,7 +66,7 @@ private: TangentConstantWithinCell }; - float calculatePorePressureInSegment(int64_t intersectionIdx, float averageSegmentPorePressureBars, double hydroStaticPorePressureBars, double effectiveDepthMeters, const std::vector& poreElementPressuresPascal) const; + float calculatePorePressureInSegment(int64_t intersectionIdx, float averageSegmentPorePressureBar, double hydroStaticPorePressureBar, double effectiveDepthMeters, const std::vector& poreElementPressuresPascal) const; float calculatePoissonRatio(int64_t intersectionIdx, const std::vector& poissonRatios) const; float calculateUcs(int64_t intersectionIdx, const std::vector& ucsValuesPascal) const; @@ -88,11 +88,11 @@ private: static caf::Ten3d transformTensorToWellPathOrientation(const cvf::Vec3d& wellPathTangent, const caf::Ten3d& wellPathTensor); - static cvf::Vec3f cellCentroid(const int* elmNodeIndices, const std::vector& nodeCoords); + cvf::Vec3f cellCentroid(size_t intersectionIdx) const; double getWellLogSegmentValue(size_t intersectionIdx, const std::vector>& wellLogValues) const; template - static bool averageIntersectionValuesToSegmentValue(size_t intersectionIdx, const std::vector& intersectionValues, const T& invalidValue, T* averagedSegmentValue); + bool averageIntersectionValuesToSegmentValue(size_t intersectionIdx, const std::vector& intersectionValues, const T& invalidValue, T* averagedSegmentValue) const; static double pascalToBar(double pascalValue); private: cvf::ref m_caseData; diff --git a/Fwk/AppFwk/cafTensor/cafTensor3.h b/Fwk/AppFwk/cafTensor/cafTensor3.h index 7999ed4290..5e475a63b1 100644 --- a/Fwk/AppFwk/cafTensor/cafTensor3.h +++ b/Fwk/AppFwk/cafTensor/cafTensor3.h @@ -29,7 +29,7 @@ class Tensor3 { S m_tensor[6]; // SXX, SYY, SZZ, SXY, SYZ, SZX public: - Tensor3() {} + Tensor3(); Tensor3(S sxx, S syy, S szz, S sxy, S syz, S szx); Tensor3(const Tensor3& other); template diff --git a/Fwk/AppFwk/cafTensor/cafTensor3.inl b/Fwk/AppFwk/cafTensor/cafTensor3.inl index 2d0e0027ce..2bfe705b08 100644 --- a/Fwk/AppFwk/cafTensor/cafTensor3.inl +++ b/Fwk/AppFwk/cafTensor/cafTensor3.inl @@ -27,6 +27,19 @@ namespace caf { +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +template< typename S> +caf::Tensor3::Tensor3() +{ + m_tensor[0] = (S) 0.0; + m_tensor[1] = (S) 0.0; + m_tensor[2] = (S) 0.0; + m_tensor[3] = (S) 0.0; + m_tensor[4] = (S) 0.0; + m_tensor[5] = (S) 0.0; +} //---------------------------------------------------------------------------------------------------- /// Copy constructor