Files
ResInsight/ApplicationLibCode/Commands/CompletionExportCommands/MswExport/RicMswBranchBuilder.cpp
T
Magne Sjaastad 6e8e658f6d #14980: Remove global max segment length from segment collection
The global max segment length no longer exists, so segment subdivision is defined only by the segment intervals.
2026-10-19 21:13:13 +02:00

1272 lines
63 KiB
C++

/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2026 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 "RicMswBranchBuilder.h"
#include "RicExportFractureCompletionsImpl.h"
#include "RicMswTableDataTools.h"
#include "RicPerforationCellFilterEvaluator.h"
#include "RimCellFilter.h"
#include "RigActiveCellInfo.h"
#include "RigEclipseCaseData.h"
#include "RigGridBase.h"
#include "RigMainGrid.h"
#include "Well/RigWellPathIntersectionTools.h"
#include "RimEclipseCase.h"
#include "RimFishbones.h"
#include "RimFishbonesCollection.h"
#include "RimFractureTemplate.h"
#include "RimPerforationCollection.h"
#include "RimPerforationInterval.h"
#include "RimSegmentCollection.h"
#include "RimWellPath.h"
#include "RimWellPathCompletions.h"
#include "RimWellPathFracture.h"
#include "RimWellPathFractureCollection.h"
#include "RimWellPathValve.h"
#include "RiaLogging.h"
#include "RiaQStringFormatter.h"
#include "Well/RigWellPath.h"
#include <algorithm>
#include <cmath>
#include <limits>
#include <map>
#include <set>
namespace RicMswBranchBuilder
{
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void applyEffectiveDiameters( const FishbonesExportContext& context, RigMswWellExportData& exportData )
{
if ( context.laterals.empty() ) return;
// Deff = sqrt(d1^2 + d2^2 + ..) for all cell intersections sharing a grid cell. A cell
// intersection contributes once, no matter how many WELSEGS rows it was split into.
std::map<size_t, double> sumOfSquaresPerCell;
for ( const auto& lateral : context.laterals )
{
for ( const auto& lateralSegment : lateral.segments )
{
sumOfSquaresPerCell[lateralSegment.globalCellIndex] += lateralSegment.equivalentDiameter * lateralSegment.equivalentDiameter;
}
}
std::map<int, double> effectiveDiameterPerSegment;
for ( const auto& lateral : context.laterals )
{
std::vector<double> effectiveDiameters;
for ( const auto& lateralSegment : lateral.segments )
{
effectiveDiameters.push_back( std::sqrt( sumOfSquaresPerCell[lateralSegment.globalCellIndex] ) );
}
// Reduce the diameter of the cell intersection shared with the main bore
if ( effectiveDiameters.size() > 1 ) effectiveDiameters[0] = effectiveDiameters[1];
for ( size_t i = 0; i < lateral.segments.size(); i++ )
{
for ( auto segmentNumber : lateral.segments[i].segmentNumbers )
{
effectiveDiameterPerSegment[segmentNumber] = effectiveDiameters[i];
}
}
}
for ( auto& branch : exportData.branches )
{
for ( auto& segment : branch.segments )
{
auto it = effectiveDiameterPerSegment.find( segment.segmentNumber );
if ( it != effectiveDiameterPerSegment.end() && it->second > 0.0 ) segment.diameter = it->second;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void applyIcdAreaPerCell( const FishbonesExportContext& context, RigMswWellExportData& exportData )
{
if ( context.icdSegments.empty() ) return;
std::set<int> icdSegmentNumbers;
for ( const auto& icdSegment : context.icdSegments )
{
icdSegmentNumbers.insert( icdSegment.segmentNumber );
}
std::map<int, double> originalAreaPerSegment;
for ( const auto& branch : exportData.branches )
{
for ( const auto& segment : branch.segments )
{
if ( !segment.wsegvalvData.has_value() ) continue;
if ( !icdSegmentNumbers.count( segment.segmentNumber ) ) continue;
originalAreaPerSegment[segment.segmentNumber] = segment.wsegvalvData->area;
}
}
// Every sum is computed from the original areas, so the result does not depend on the order the
// cells are visited in.
std::map<size_t, double> areaSumPerCell;
for ( const auto& icdSegment : context.icdSegments )
{
auto it = originalAreaPerSegment.find( icdSegment.segmentNumber );
if ( it == originalAreaPerSegment.end() ) continue;
for ( auto globalCellIndex : icdSegment.globalCellIndices )
{
areaSumPerCell[globalCellIndex] += it->second;
}
}
// An ICD sub spans a 0.1 m valve segment and normally connects to a single cell. When it reaches
// into more than one, it reports the largest of the sums it takes part in.
std::map<int, double> areaPerSegment;
for ( const auto& icdSegment : context.icdSegments )
{
double area = 0.0;
for ( auto globalCellIndex : icdSegment.globalCellIndices )
{
area = std::max( area, areaSumPerCell[globalCellIndex] );
}
if ( area > 0.0 ) areaPerSegment[icdSegment.segmentNumber] = area;
}
for ( auto& branch : exportData.branches )
{
for ( auto& segment : branch.segments )
{
if ( !segment.wsegvalvData.has_value() ) continue;
auto it = areaPerSegment.find( segment.segmentNumber );
if ( it != areaPerSegment.end() ) segment.wsegvalvData->area = it->second;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::optional<RigMswCellIntersection>
toMswCellIntersection( const WellPathCellIntersectionInfo& cellInfo, const RigMainGrid* mainGrid, double distanceStart, double distanceEnd )
{
if ( cellInfo.globCellIndex >= mainGrid->totalCellCount() ) return std::nullopt;
size_t localIdx = 0;
const RigGridBase* localGrid = mainGrid->gridAndGridLocalIdxFromGlobalCellIdx( cellInfo.globCellIndex, &localIdx );
size_t i, j, k;
localGrid->ijkFromCellIndex( localIdx, &i, &j, &k );
if ( mainGrid->isDualPorosity() ) k += mainGrid->cellCountK();
RigMswCellIntersection ci;
ci.i = i + 1; // 1-based
ci.j = j + 1;
ci.k = k + 1;
ci.distanceStart = distanceStart;
ci.distanceEnd = distanceEnd;
if ( localGrid != mainGrid ) ci.gridName = localGrid->gridName();
return ci;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double segmentNodeMD( const std::vector<CellSegmentEntry>& cellSegMap, int segmentNumber )
{
for ( const auto& entry : cellSegMap )
{
if ( entry.lastSubSegmentNumber == segmentNumber ) return entry.nodeMD;
}
return 0.0;
}
// Range-based lookup: find the segment whose [cellStartMD, cellEndMD] range contains md.
// Used for ICD valves where the valve overlap start always falls within the host main-bore cell.
// Fractures use findOutletSegmentForMD (midpoint-based) to match the tree approach.
int findContainingSegmentForMD( const std::vector<CellSegmentEntry>& cellSegMap, double md )
{
if ( cellSegMap.empty() ) return 1;
for ( const auto& entry : cellSegMap )
{
if ( md >= entry.cellStartMD && md < entry.cellEndMD ) return entry.lastSubSegmentNumber;
}
// Fallback: md is beyond all segments — return the last segment.
return cellSegMap.back().lastSubSegmentNumber;
}
int findOutletSegmentForMD( const std::vector<CellSegmentEntry>& cellSegMap, double md )
{
if ( cellSegMap.empty() ) return 1;
// Match the tree-based logic (findClosestSegmentWithLowerMD): find the segment whose node is
// closest to, but not greater than, the given MD. Range-based lookup is
// incorrect when a segment node lies above the completion MD (e.g. a fracture near the
// heel of a long cell), which would assign the wrong outlet segment.
const CellSegmentEntry* best = nullptr;
double bestDist = std::numeric_limits<double>::infinity();
for ( const auto& entry : cellSegMap )
{
const double dist = md - entry.nodeMD;
if ( dist >= 0.0 && dist < bestDist )
{
best = &entry;
bestDist = dist;
}
}
// Fallback: md is shallower than all segment nodes — connect to the first (shallowest) segment.
return best ? best->lastSubSegmentNumber : cellSegMap.front().lastSubSegmentNumber;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<NodeSpan> segmentSpans( const std::vector<std::pair<double, double>>& cellPieces, const std::vector<double>& nodes )
{
std::vector<NodeSpan> parts;
if ( nodes.empty() ) return parts;
// Nodes unchanged by the segmentation rules: each piece is served by its own centre node
bool nodesAreCandidates = nodes.size() == cellPieces.size();
for ( size_t i = 0; nodesAreCandidates && i < nodes.size(); ++i )
{
nodesAreCandidates = nodes[i] == 0.5 * ( cellPieces[i].first + cellPieces[i].second );
}
if ( nodesAreCandidates )
{
for ( size_t i = 0; i < cellPieces.size(); ++i )
{
parts.push_back( { cellPieces[i].first, cellPieces[i].second, i } );
}
return parts;
}
// Split the pieces at the nodes, each part belongs to the node upstream of it
for ( const auto& [pieceStart, pieceEnd] : cellPieces )
{
std::vector<double> cuts = { pieceStart };
for ( const double nodeMD : nodes )
{
if ( nodeMD > pieceStart && nodeMD < pieceEnd ) cuts.push_back( nodeMD );
}
std::ranges::sort( cuts );
cuts.push_back( pieceEnd );
for ( size_t k = 0; k + 1 < cuts.size(); ++k )
{
parts.push_back( { cuts[k], cuts[k + 1], *RicMswTableDataTools::upstreamNodeIndex( nodes, cuts[k] ) } );
}
}
// Merge consecutive parts served by the same node
std::vector<NodeSpan> spans;
for ( const auto& part : parts )
{
if ( !spans.empty() && spans.back().nodeIndex == part.nodeIndex )
{
spans.back().endMD = part.endMD;
}
else
{
spans.push_back( part );
}
}
return spans;
}
namespace
{
//--------------------------------------------------------------------------------------------------
/// Data shared by the COMPSEGS collection of all cells along a main bore branch
//--------------------------------------------------------------------------------------------------
struct PerforationCompsegsContext
{
const std::vector<const RimPerforationInterval*>& perforationIntervals;
const std::set<const RimPerforationInterval*>& valvedIntervals;
const std::optional<QDateTime>& exportDate;
const RigMainGrid* mainGrid;
const RigActiveCellInfo* activeCellInfo;
std::map<const RimPerforationInterval*, std::set<size_t>>& perforationCompletionCells;
std::map<const RimPerforationInterval*, RicPerforationCellFilterEvaluator>& filterEvaluators;
std::map<const RimPerforationInterval*, size_t>& overlapCountPerInterval;
std::map<const RimPerforationInterval*, size_t>& emittedCountPerInterval;
};
//--------------------------------------------------------------------------------------------------
/// Collect COMPSEGS only for bare perforations (no active valve) on the main bore.
/// Valved perforations get their COMPSEGS on the valve segment.
//--------------------------------------------------------------------------------------------------
std::vector<RigMswCellIntersection> collectCellCompsegs( const WellPathCellIntersectionInfo& cellInfo, PerforationCompsegsContext& context )
{
std::vector<RigMswCellIntersection> cellCompsegs;
for ( const auto* perf : context.perforationIntervals )
{
if ( context.valvedIntervals.count( perf ) ) continue; // valve handles this interval
if ( context.exportDate.has_value() && !perf->isActiveOnDate( context.exportDate.value() ) ) continue;
const double overlapStart = std::max( perf->startMD(), cellInfo.startMD );
const double overlapEnd = std::min( perf->endMD(), cellInfo.endMD );
if ( overlapEnd > overlapStart )
{
if ( context.activeCellInfo && cellInfo.globCellIndex < context.mainGrid->totalCellCount() &&
!context.activeCellInfo->isActive( ReservoirCellIndex( cellInfo.globCellIndex ) ) )
continue;
// Skip cells the perforation does not actually connect to (see perforationCompletionCells).
const auto& completionCells = context.perforationCompletionCells[perf];
if ( !completionCells.empty() && !completionCells.count( cellInfo.globCellIndex ) ) continue;
context.overlapCountPerInterval[perf]++;
auto& filterEval = context.filterEvaluators.at( perf );
if ( !filterEval.includesGlobalCell( cellInfo.globCellIndex ) ) continue;
if ( auto ci = toMswCellIntersection( cellInfo, context.mainGrid, overlapStart, overlapEnd ) )
{
// All candidates here describe the same grid cell, as they are derived from a single cellInfo.
// Several perforation intervals may overlap that cell, but the cell must be connected only once,
// matching the single COMPDAT connection. Widen the existing range instead of adding another row.
if ( cellCompsegs.empty() )
{
cellCompsegs.push_back( *ci );
}
else
{
auto& existing = cellCompsegs.front();
existing.distanceStart = std::min( existing.distanceStart, ci->distanceStart );
existing.distanceEnd = std::max( existing.distanceEnd, ci->distanceEnd );
}
context.emittedCountPerInterval[perf]++;
}
}
}
return cellCompsegs;
}
} // namespace
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigMswBranch buildMainBoreBranch( const RimWellPath* wellPath,
const std::vector<WellPathCellIntersectionInfo>& filteredIntersections,
const RigMainGrid* mainGrid,
const RimEclipseCase* eclipseCase,
const std::vector<const RimPerforationInterval*>& perforationIntervals,
const std::set<const RimPerforationInterval*>& valvedIntervals,
const std::string& infoType,
double heelMD,
double heelTVD,
int branchNumber,
int& segmentNumber,
int outletSegmentNumber,
const std::vector<std::pair<double, double>>& customSegmentIntervals,
const RicMswTableDataTools::SegmentationSettings& segmentationSettings,
const std::optional<QDateTime>& exportDate,
RiaDefines::EclipseUnitSystem unitSystem,
std::vector<CellSegmentEntry>* cellSegMap,
const RigActiveCellInfo* activeCellInfo )
{
std::vector<RigMswSegment> result;
const double segmentLengthThreshold = 1.0e-3;
double prevOutMD = heelMD;
double prevOutTVD = heelTVD;
int prevSegNum = outletSegmentNumber;
bool firstSeg = true;
std::map<const RimPerforationInterval*, RicPerforationCellFilterEvaluator> filterEvaluators;
std::map<const RimPerforationInterval*, size_t> overlapCountPerInterval;
std::map<const RimPerforationInterval*, size_t> emittedCountPerInterval;
for ( const auto* perf : perforationIntervals )
{
filterEvaluators.emplace( std::piecewise_construct,
std::forward_as_tuple( perf ),
std::forward_as_tuple( perf->cellFilter(), eclipseCase ) );
overlapCountPerInterval[perf] = 0;
emittedCountPerInterval[perf] = 0;
}
// For each perforation interval, determine the set of grid cells the perforation actually
// connects to. This mirrors the standard COMPDAT export (generatePerforationsCompdatValues):
// the well path is clipped to the perforation MD range and re-intersected with the grid.
// COMPSEGS rows are then emitted only for cells in this set.
//
// This matters for dual-porosity models, where the whole-well geometric intersection
// (filteredIntersections) may split a perforation across several matrix cells at cell
// boundaries, while the clipped perforation intersection collapses it to the cells that
// carry the connection. Restricting to the latter keeps the MSW COMPSEGS consistent with
// the COMPDAT connections. For single-porosity models the two sets are identical, so this
// is a no-op.
std::map<const RimPerforationInterval*, std::set<size_t>> perforationCompletionCells;
if ( const auto* wellPathGeometry = wellPath->wellPathGeometry() )
{
for ( const auto* perf : perforationIntervals )
{
const auto [clipCoords, clipMds] = wellPathGeometry->clippedPointSubset( perf->startMD(), perf->endMD() );
const auto perfCells = RigWellPathIntersectionTools::findCellIntersectionInfosAlongPath( eclipseCase->eclipseCaseData(),
wellPath->name(),
clipCoords,
clipMds );
std::set<size_t>& cellSet = perforationCompletionCells[perf];
for ( const auto& perfCell : perfCells )
{
const bool cellIsActive = !activeCellInfo || perfCell.globCellIndex >= mainGrid->totalCellCount() ||
activeCellInfo->isActive( ReservoirCellIndex( perfCell.globCellIndex ) );
if ( cellIsActive ) cellSet.insert( perfCell.globCellIndex );
}
}
}
PerforationCompsegsContext compsegsContext{ perforationIntervals,
valvedIntervals,
exportDate,
mainGrid,
activeCellInfo,
perforationCompletionCells,
filterEvaluators,
overlapCountPerInterval,
emittedCountPerInterval };
// Pass 1: the cell pieces give the candidate nodes at their centres, the segmentation rules place the nodes
std::vector<std::pair<double, double>> cellPieces;
std::vector<std::vector<RigMswCellIntersection>> compsegsPerCell;
std::vector<double> candidateNodes;
for ( const auto& cellInfo : filteredIntersections )
{
const double cellLength = std::fabs( cellInfo.endMD - cellInfo.startMD );
if ( cellLength <= segmentLengthThreshold )
{
RiaLogging::info( std::format( "Skipping segment, threshold = {}, length = {}", segmentLengthThreshold, cellLength ) );
continue;
}
compsegsPerCell.push_back( collectCellCompsegs( cellInfo, compsegsContext ) );
const auto subSegs = RicMswTableDataTools::createSubSegmentMDPairs( cellInfo.startMD, cellInfo.endMD, customSegmentIntervals );
for ( const auto& [subStartMD, subEndMD] : subSegs )
{
cellPieces.push_back( { subStartMD, subEndMD } );
candidateNodes.push_back( 0.5 * ( subStartMD + subEndMD ) );
}
}
std::vector<double> nodes = RicMswTableDataTools::placeSegmentNodes( heelMD,
candidateNodes,
segmentationSettings.intervals,
cellPieces.empty() ? heelMD : cellPieces.back().second,
cellPieces.empty() ? heelMD : cellPieces.front().first );
if ( !cellPieces.empty() && !perforationIntervals.empty() &&
( segmentationSettings.singleSegmentBeforeFirstPerforation || segmentationSettings.singleSegmentAfterLastPerforation ) )
{
double firstPerforationMD = std::numeric_limits<double>::infinity();
double lastPerforationMD = -std::numeric_limits<double>::infinity();
for ( const auto* perf : perforationIntervals )
{
if ( exportDate && !perf->isActiveOnDate( *exportDate ) ) continue;
firstPerforationMD = std::min( firstPerforationMD, perf->startMD() );
lastPerforationMD = std::max( lastPerforationMD, perf->endMD() );
}
nodes = RicMswTableDataTools::collapseNodesOutsidePerforations( nodes,
cellPieces.front().first,
cellPieces.back().second,
firstPerforationMD,
lastPerforationMD,
segmentationSettings.singleSegmentBeforeFirstPerforation,
segmentationSettings.singleSegmentAfterLastPerforation );
}
// Pass 2: one segment per node, the outlet is the previous node
const int firstSegmentNumber = segmentNumber;
for ( const double nodeMD : nodes )
{
const double nodeTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, nodeMD );
double curPrevOutMD = prevOutMD;
double curPrevOutTVD = prevOutTVD;
if ( nodeMD < curPrevOutMD )
{
curPrevOutMD = heelMD;
curPrevOutTVD = heelTVD;
}
double length = 0.0;
double depth = 0.0;
if ( infoType == "INC" )
{
length = nodeMD - curPrevOutMD;
depth = nodeTVD - curPrevOutTVD;
}
else
{
length = nodeMD;
depth = nodeTVD;
}
double diameter = 0.0;
double roughness = 0.0;
if ( exportDate.has_value() )
{
diameter = wellPath->segmentCollection()->getDiameterAtMD( nodeMD, unitSystem, *exportDate );
roughness = wellPath->segmentCollection()->getRoughnessAtMD( nodeMD, unitSystem, *exportDate );
}
else
{
diameter = wellPath->segmentCollection()->getDiameterAtMD( nodeMD, unitSystem );
roughness = wellPath->segmentCollection()->getRoughnessAtMD( nodeMD, unitSystem );
}
RigMswSegment seg;
seg.segmentNumber = segmentNumber;
seg.outletSegmentNumber = prevSegNum;
seg.length = length;
seg.depth = depth;
seg.diameter = diameter;
seg.roughness = roughness;
seg.sourceWellName = wellPath->name().toStdString();
if ( firstSeg )
{
seg.description = "Segments on main bore";
firstSeg = false;
}
prevSegNum = segmentNumber++;
prevOutMD = nodeMD;
prevOutTVD = nodeTVD;
result.push_back( std::move( seg ) );
}
// COMPSEGS rows go to the segment with the node nearest the row centre, as the simulator assigns them
for ( const auto& cellCompsegs : compsegsPerCell )
{
for ( const auto& ci : cellCompsegs )
{
if ( auto nodeIndex = RicMswTableDataTools::nearestNodeIndex( nodes, 0.5 * ( ci.distanceStart + ci.distanceEnd ) ) )
{
result[*nodeIndex].intersections.push_back( ci );
}
}
}
if ( cellSegMap )
{
for ( const auto& span : segmentSpans( cellPieces, nodes ) )
{
cellSegMap->push_back(
{ span.startMD, span.endMD, firstSegmentNumber + static_cast<int>( span.nodeIndex ), nodes[span.nodeIndex] } );
}
}
for ( const auto* perf : perforationIntervals )
{
const auto& filterEval = filterEvaluators.at( perf );
if ( !filterEval.isEnabled() ) continue;
if ( overlapCountPerInterval[perf] == 0 ) continue;
if ( emittedCountPerInterval[perf] > 0 ) continue;
RiaLogging::info( std::format( "Perforation '{}' on well '{}': all {} intersected cells suppressed by attached cell filter '{}'.",
perf->name(),
wellPath->name(),
filterEval.rejectedCellCount(),
perf->cellFilter() ? perf->cellFilter()->fullName() : QString() ) );
}
return RigMswBranch{ branchNumber, std::nullopt, std::move( result ) };
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RigMswBranch> buildValveBranches( const RimWellPath* wellPath,
const std::vector<WellPathCellIntersectionInfo>& filteredIntersections,
const RigMainGrid* mainGrid,
const RimEclipseCase* eclipseCase,
const std::vector<const RimPerforationInterval*>& perforationIntervals,
const std::vector<CellSegmentEntry>& cellSegMap,
const std::string& infoType,
const std::string& wellNameForExport,
int& segmentNumber,
int& branchNumber,
const std::vector<std::pair<double, double>>& customSegmentIntervals,
const std::optional<QDateTime>& exportDate,
RiaDefines::EclipseUnitSystem unitSystem )
{
std::vector<RigMswBranch> result;
const auto linerDiameter = wellPath->segmentCollection()->linerDiameter( unitSystem );
const auto roughnessFactor = wellPath->segmentCollection()->roughnessFactor( unitSystem );
for ( const auto* perf : perforationIntervals )
{
if ( exportDate.has_value() && !perf->isActiveOnDate( exportDate.value() ) ) continue;
RicPerforationCellFilterEvaluator filterEval( perf->cellFilter(), eclipseCase );
auto valves = perf->descendantsIncludingThisOfType<RimWellPathValve>();
for ( const auto* valve : valves )
{
if ( !valve->isChecked() ) continue;
if ( exportDate.has_value() && !valve->isActiveOnDate( exportDate.value() ) ) continue;
const auto valveType = valve->componentType();
const bool isIcdOrIcv =
( valveType == RiaDefines::WellPathComponentType::ICD || valveType == RiaDefines::WellPathComponentType::ICV );
const bool isAicdOrSicd =
( valveType == RiaDefines::WellPathComponentType::AICD || valveType == RiaDefines::WellPathComponentType::SICD );
if ( !isIcdOrIcv && !isAicdOrSicd ) continue;
const auto& valveLocations = valve->valveLocations();
const auto& valveSegs = valve->valveSegments(); // one (startMD, endMD) per location
// ICD: one branch per overlapping cell with proportional area, matching tree approach
if ( valveType == RiaDefines::WellPathComponentType::ICD )
{
for ( size_t vi = 0; vi < valveLocations.size(); ++vi )
{
const double valveSegStart = valveSegs[vi].first;
const double valveSegEnd = valveSegs[vi].second;
// Collect cells overlapping this ICD's coverage range and total overlap length
struct CellEntry
{
RigMswCellIntersection ci;
double overlapMD;
};
std::vector<CellEntry> cells;
double totalOverlapMD = 0.0;
for ( const auto& cellInfo : filteredIntersections )
{
const double overlapStart = std::max( valveSegStart, cellInfo.startMD );
const double overlapEnd = std::min( valveSegEnd, cellInfo.endMD );
if ( overlapEnd <= overlapStart ) continue;
if ( !filterEval.includesGlobalCell( cellInfo.globCellIndex ) ) continue;
if ( auto ci = toMswCellIntersection( cellInfo, mainGrid, overlapStart, overlapEnd ) )
{
double overlapMD = overlapEnd - overlapStart;
cells.push_back( { *ci, overlapMD } );
totalOverlapMD += overlapMD;
}
}
// One branch+segment per cell with area proportional to overlap length
for ( auto& entry : cells )
{
const int cellBranch = ++branchNumber;
const int outletSeg = findContainingSegmentForMD( cellSegMap, entry.ci.distanceStart );
const double valveMD = segmentNodeMD( cellSegMap, outletSeg );
const double valveEndMD = valveMD + RicMswTableDataTools::valveSegmentLength;
const double startTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, valveMD );
const double endTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, valveEndMD );
double length = 0.0, depth = 0.0;
if ( infoType == "INC" )
{
length = valveEndMD - valveMD;
depth = endTVD - startTVD;
}
else
{
length = valveEndMD;
depth = endTVD;
}
RigMswSegment seg;
seg.segmentNumber = segmentNumber;
seg.outletSegmentNumber = outletSeg;
seg.length = length;
seg.depth = depth;
seg.diameter = linerDiameter;
seg.roughness = roughnessFactor;
seg.sourceWellName = wellPath->name().toStdString();
seg.description = QString( "%1 #%2" ).arg( valve->name() ).arg( vi + 1 ).toStdString();
seg.intersections = { entry.ci };
WsegvalvRow wv;
wv.well = wellNameForExport;
wv.segmentNumber = segmentNumber;
wv.cv = valve->flowCoefficient();
wv.area = totalOverlapMD > 0.0 ? valve->area( unitSystem ) * entry.overlapMD / totalOverlapMD : 0.0;
wv.status = valve->isOpen() ? "OPEN" : "SHUT";
wv.description = valve->name().toStdString();
seg.wsegvalvData = wv;
segmentNumber++;
result.push_back( RigMswBranch{ cellBranch, std::nullopt, { std::move( seg ) } } );
}
}
continue;
}
// AICD: one branch per overlapping cell with proportional flow scaling factor, matching tree approach
if ( valveType == RiaDefines::WellPathComponentType::AICD )
{
const RimWellPathAicdParameters* aicdParams = valve->aicdParameters();
if ( aicdParams && aicdParams->isValid() )
{
const auto aicdValues = aicdParams->doubleValues();
auto setOptional = []( double value ) -> std::optional<double>
{ return std::isinf( value ) ? std::nullopt : std::optional<double>{ value }; };
for ( size_t vi = 0; vi < valveLocations.size(); ++vi )
{
const double valveSegStart = valveSegs[vi].first;
const double valveSegEnd = valveSegs[vi].second;
struct CellEntry
{
RigMswCellIntersection ci;
double overlapMD;
};
std::vector<CellEntry> cells;
double totalOverlapMD = 0.0;
for ( const auto& cellInfo : filteredIntersections )
{
const double overlapStart = std::max( valveSegStart, cellInfo.startMD );
const double overlapEnd = std::min( valveSegEnd, cellInfo.endMD );
if ( overlapEnd <= overlapStart ) continue;
if ( !filterEval.includesGlobalCell( cellInfo.globCellIndex ) ) continue;
if ( auto ci = toMswCellIntersection( cellInfo, mainGrid, overlapStart, overlapEnd ) )
{
cells.push_back( { *ci, overlapEnd - overlapStart } );
totalOverlapMD += overlapEnd - overlapStart;
}
}
for ( auto& entry : cells )
{
const int cellBranch = ++branchNumber;
const int outletSeg = findContainingSegmentForMD( cellSegMap, entry.ci.distanceStart );
const double valveMD = segmentNodeMD( cellSegMap, outletSeg );
const double valveEndMD = valveMD + RicMswTableDataTools::valveSegmentLength;
const double startTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, valveMD );
const double endTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, valveEndMD );
double length = 0.0, depth = 0.0;
if ( infoType == "INC" )
{
length = valveEndMD - valveMD;
depth = endTVD - startTVD;
}
else
{
length = valveEndMD;
depth = endTVD;
}
RigMswSegment seg;
seg.segmentNumber = segmentNumber;
seg.outletSegmentNumber = outletSeg;
seg.length = length;
seg.depth = depth;
seg.diameter = linerDiameter;
seg.roughness = roughnessFactor;
seg.sourceWellName = wellPath->name().toStdString();
seg.description = QString( "%1 #%2" ).arg( valve->name() ).arg( vi + 1 ).toStdString();
seg.intersections = { entry.ci };
WsegaicdRow wa;
wa.well = wellNameForExport;
wa.segment1 = segmentNumber;
wa.segment2 = segmentNumber;
wa.strength = aicdValues[AICD_STRENGTH];
wa.length = totalOverlapMD > 0.0 ? totalOverlapMD / entry.overlapMD : 0.0;
wa.densityCali = setOptional( aicdValues[AICD_DENSITY_CALIB_FLUID] );
wa.viscosityCali = setOptional( aicdValues[AICD_VISCOSITY_CALIB_FLUID] );
wa.criticalValue = setOptional( aicdValues[AICD_CRITICAL_WATER_IN_LIQUID_FRAC] );
wa.widthTrans = setOptional( aicdValues[AICD_EMULSION_VISC_TRANS_REGION] );
wa.maxViscRatio = setOptional( aicdValues[AICD_MAX_RATIO_EMULSION_VISC] );
wa.methodScalingFactor = 1;
wa.maxAbsRate = aicdValues[AICD_MAX_FLOW_RATE];
wa.flowRateExponent = aicdValues[AICD_VOL_FLOW_EXP];
wa.viscExponent = aicdValues[AICD_VISOSITY_FUNC_EXP];
wa.status = aicdParams->isOpen() ? "OPEN" : "SHUT";
wa.description = valve->name().toStdString();
seg.wsegaicdData = wa;
segmentNumber++;
result.push_back( RigMswBranch{ cellBranch, std::nullopt, { std::move( seg ) } } );
}
}
}
continue;
}
// ICV: one branch, one segment per valve location
const int valveBranch = ++branchNumber;
std::vector<RigMswSegment> branchSegs;
for ( size_t vi = 0; vi < valveLocations.size(); ++vi )
{
const double valveMD = valveLocations[vi];
const double valveEndMD = valveMD + RicMswTableDataTools::valveSegmentLength;
const double valveSegStart = valveSegs[vi].first;
const double valveSegEnd = valveSegs[vi].second;
const int outletSeg = findOutletSegmentForMD( cellSegMap, valveMD );
const double valveStartTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, valveMD );
const double valveEndTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, valveEndMD );
double length = 0.0;
double depth = 0.0;
if ( infoType == "INC" )
{
length = valveEndMD - valveMD;
depth = valveEndTVD - valveStartTVD;
}
else
{
length = valveEndMD;
depth = valveEndTVD;
}
RigMswSegment valveSeg;
valveSeg.segmentNumber = segmentNumber;
valveSeg.outletSegmentNumber = outletSeg;
valveSeg.length = length;
valveSeg.depth = depth;
valveSeg.diameter = linerDiameter;
valveSeg.roughness = roughnessFactor;
valveSeg.sourceWellName = wellPath->name().toStdString();
valveSeg.description = QString( "%1 #%2" ).arg( valve->name() ).arg( vi + 1 ).toStdString();
// COMPSEGS: cells that overlap this valve's coverage range (valveSegStart..valveSegEnd).
// The perforation intersections are redirected onto the ICV segment, so their COMPSEGS
// Start/End Length reflect the valve segment MD range (valveMD..valveEndMD) rather than
// the individual cell overlaps, matching the legacy tree-based export.
for ( const auto& cellInfo : filteredIntersections )
{
const double overlapStart = std::max( valveSegStart, cellInfo.startMD );
const double overlapEnd = std::min( valveSegEnd, cellInfo.endMD );
if ( overlapEnd > overlapStart )
{
if ( !filterEval.includesGlobalCell( cellInfo.globCellIndex ) ) continue;
if ( auto ci = toMswCellIntersection( cellInfo, mainGrid, valveMD, valveEndMD ) )
valveSeg.intersections.push_back( *ci );
}
}
// WSEGVALV for ICV
WsegvalvRow wv;
wv.well = wellNameForExport;
wv.segmentNumber = segmentNumber;
wv.cv = valve->flowCoefficient();
wv.area = valve->area( unitSystem );
wv.status = valve->isOpen() ? "OPEN" : "SHUT";
wv.description = valve->name().toStdString();
valveSeg.wsegvalvData = wv;
segmentNumber++;
branchSegs.push_back( std::move( valveSeg ) );
}
result.push_back( RigMswBranch{ valveBranch, std::nullopt, std::move( branchSegs ) } );
}
}
return result;
}
//--------------------------------------------------------------------------------------------------
/// Measured depth where the fracture segment starts. A fracture along the well path is centred on
/// the fracture MD, and starts half a perforation length above it.
//--------------------------------------------------------------------------------------------------
double fractureStartMD( const RimWellPathFracture* fracture )
{
double startMD = fracture->fractureMD();
if ( fracture->fractureTemplate()->orientationType() == RimFractureTemplate::ALONG_WELL_PATH )
{
startMD -= 0.5 * fracture->fractureTemplate()->perforationLength();
}
return startMD;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RigMswBranch> buildFractureBranches( RimEclipseCase* eclipseCase,
const RimWellPath* wellPath,
const RigMainGrid* mainGrid,
const std::vector<CellSegmentEntry>& cellSegMap,
const std::string& infoType,
int& segmentNumber,
int& branchNumber )
{
std::vector<RigMswBranch> result;
const QString wellNameForExport = wellPath->completionSettings()->wellNameForExport();
// Number the fracture branches by position along the well path, so that the branch numbers follow
// the well no matter in what order the fractures were created. The order also decides which
// fracture connects a grid cell shared by two of them.
auto fractures = wellPath->fractureCollection()->activeFractures();
std::stable_sort( fractures.begin(),
fractures.end(),
[]( const RimWellPathFracture* lhs, const RimWellPathFracture* rhs )
{ return fractureStartMD( lhs ) < fractureStartMD( rhs ); } );
for ( RimWellPathFracture* fracture : fractures )
{
fracture->ensureValidNonDarcyProperties();
const double position = fractureStartMD( fracture );
double width = fracture->fractureTemplate()->computeFractureWidth( fracture );
if ( fracture->fractureTemplate()->orientationType() == RimFractureTemplate::ALONG_WELL_PATH )
{
width = fracture->fractureTemplate()->perforationLength();
}
const double endMD = position + width;
const double startTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, position );
const double endTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, endMD );
double length = 0.0;
double depth = 0.0;
if ( infoType == "INC" )
{
length = width;
depth = endTVD - startTVD;
}
else
{
length = endMD;
depth = endTVD;
}
const int outletSeg = findOutletSegmentForMD( cellSegMap, position );
const int fracBranch = ++branchNumber;
// Get COMPSEGS cell intersections from fracture completion data
std::vector<RigCompletionData> completionData = RicExportFractureCompletionsImpl::generateCompdatValues( eclipseCase,
wellNameForExport,
wellPath->wellPathGeometry(),
{ fracture },
nullptr,
nullptr );
std::vector<RigMswCellIntersection> compsegs;
for ( const auto& compData : completionData )
{
const auto& cell = compData.completionDataGridCell();
// 0-based -> 1-based
size_t i = cell.localCellIndexI() + 1;
size_t j = cell.localCellIndexJ() + 1;
size_t k = cell.localCellIndexK() + 1;
// Shift K for dual porosity models
if ( mainGrid->isDualPorosity() ) k += mainGrid->cellCountK();
RigMswCellIntersection ci;
ci.i = i;
ci.j = j;
ci.k = k;
ci.distanceStart = position;
ci.distanceEnd = endMD;
ci.gridName = cell.lgrName().toStdString();
compsegs.push_back( ci );
}
RigMswSegment seg;
seg.segmentNumber = segmentNumber++;
seg.outletSegmentNumber = outletSeg;
seg.length = length;
seg.depth = depth;
seg.diameter = 0.15; // default segment diameter
seg.roughness = 5.0e-5; // default segment roughness
seg.sourceWellName = wellPath->name().toStdString();
seg.description = fracture->name().toStdString();
seg.intersections = std::move( compsegs );
result.push_back( RigMswBranch{ fracBranch, std::nullopt, { std::move( seg ) }, RigMswBranchSource::Fracture } );
}
return result;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<RigMswBranch> buildFishbonesBranches( const RimEclipseCase* eclipseCase,
const RimWellPath* wellPath,
const RigMainGrid* mainGrid,
const std::vector<WellPathCellIntersectionInfo>& filteredIntersections,
const std::vector<CellSegmentEntry>& cellSegMap,
const std::string& infoType,
const std::string& wellNameForExport,
int& segmentNumber,
int& branchNumber,
const std::vector<std::pair<double, double>>& customSegmentIntervals,
RiaDefines::EclipseUnitSystem unitSystem,
FishbonesExportContext& fishbonesContext )
{
std::vector<RigMswBranch> result;
auto fishbonesSubs = wellPath->completions()->fishbonesCollection()->activeFishbonesSubs();
for ( const RimFishbones* subs : fishbonesSubs )
{
// Group laterals by sub index
std::map<size_t, std::vector<size_t>> subAndLateralIndices;
for ( const auto& [subIndex, lateralIndex] : subs->installedLateralIndices() )
subAndLateralIndices[subIndex].push_back( lateralIndex );
// Identify which filtered intersections are closest to each sub, and which
// intersection contains the sub's MD. This mirrors appendFishbonesMswExportInfo.
const double fishboneSectionStart = subs->startMD();
const double fishboneSectionEnd = subs->endMD();
std::map<size_t, std::vector<size_t>> closestSubForCellIntersections;
std::map<size_t, size_t> cellIntersectionContainingSubIndex;
for ( size_t ii = 0; ii < filteredIntersections.size(); ++ii )
{
const auto& ci = filteredIntersections[ii];
if ( fishboneSectionEnd < ci.startMD || fishboneSectionStart > ci.endMD ) continue;
const double midpoint = 0.5 * ( ci.startMD + ci.endMD );
size_t closestSubIdx = 0;
double closestDist = std::numeric_limits<double>::infinity();
for ( const auto& [subIdx, lats] : subAndLateralIndices )
{
const double subMD = subs->measuredDepth( subIdx );
if ( ci.startMD <= subMD && subMD <= ci.endMD ) cellIntersectionContainingSubIndex[subIdx] = ii;
const double dist = std::abs( subMD - midpoint );
if ( dist < closestDist )
{
closestDist = dist;
closestSubIdx = subIdx;
}
}
closestSubForCellIntersections[closestSubIdx].push_back( ii );
}
// ICD parameters
const double icdOrificeRadius = subs->icdOrificeDiameter( unitSystem ) / 2.0;
const double icdArea = icdOrificeRadius * icdOrificeRadius * cvf::PI_D * static_cast<double>( subs->icdCount() );
const double icdCv = subs->icdFlowCoefficient();
// Track used global cell indices to deduplicate COMPSEGS, matching the tree's intersectedCells logic.
std::set<size_t> usedGlobalCellIndices;
for ( const auto& [subIndex, lateralIndices] : subAndLateralIndices )
{
const double subEndMd = subs->measuredDepth( subIndex );
const double startValveMd = subEndMd - RicMswTableDataTools::valveSegmentLength;
const double startValveTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, startValveMd );
const double endValveTVD = RicMswTableDataTools::tvdFromMeasuredDepth( wellPath, subEndMd );
const int outletSeg = findOutletSegmentForMD( cellSegMap, subEndMd );
const int icdBranch = ++branchNumber;
const int icdSegNum = segmentNumber++;
// COMPSEGS for the ICD segment: cells closest to this sub, clipped to [startValveMd, subEndMd]
std::set<size_t> icdCellIndices;
if ( closestSubForCellIntersections.count( subIndex ) )
{
for ( auto idx : closestSubForCellIntersections[subIndex] )
icdCellIndices.insert( idx );
}
if ( cellIntersectionContainingSubIndex.count( subIndex ) )
icdCellIndices.insert( cellIntersectionContainingSubIndex[subIndex] );
std::vector<RigMswCellIntersection> icdCompsegs;
std::set<size_t> icdGlobalCellIndices;
for ( auto idx : icdCellIndices )
{
const auto& ci = filteredIntersections[idx];
const double overlapStart = std::max( startValveMd, ci.startMD );
const double overlapEnd = std::min( subEndMd, ci.endMD );
if ( overlapEnd <= overlapStart ) continue;
if ( auto mci = toMswCellIntersection( ci, mainGrid, overlapStart, overlapEnd ) )
{
icdCompsegs.push_back( *mci );
icdGlobalCellIndices.insert( ci.globCellIndex );
usedGlobalCellIndices.insert( ci.globCellIndex );
}
}
double icdLength = 0.0;
double icdDepth = 0.0;
if ( infoType == "INC" )
{
icdLength = subEndMd - startValveMd;
icdDepth = endValveTVD - startValveTVD;
}
else
{
icdLength = subEndMd;
icdDepth = endValveTVD;
}
RigMswSegment icdSeg;
icdSeg.segmentNumber = icdSegNum;
icdSeg.outletSegmentNumber = outletSeg;
icdSeg.length = icdLength;
icdSeg.depth = icdDepth;
icdSeg.diameter = 0.15; // default segment diameter
icdSeg.roughness = 5.0e-5; // default segment roughness
icdSeg.sourceWellName = wellPath->name().toStdString();
icdSeg.description = QString( "ICD sub %1" ).arg( subIndex + 1 ).toStdString();
icdSeg.intersections = std::move( icdCompsegs );
WsegvalvRow wv;
wv.well = wellNameForExport;
wv.segmentNumber = icdSegNum;
wv.cv = icdCv;
wv.area = icdArea;
wv.status = "OPEN";
wv.description = QString( "ICD sub %1" ).arg( subIndex + 1 ).toStdString();
icdSeg.wsegvalvData = wv;
fishbonesContext.icdSegments.push_back( { icdSegNum, std::move( icdGlobalCellIndices ) } );
result.push_back( RigMswBranch{ icdBranch, std::nullopt, { std::move( icdSeg ) }, RigMswBranchSource::Fishbones } );
// Lateral sub-segments
for ( size_t lateralIndex : lateralIndices )
{
auto lateralCoordMDPairs = subs->coordsAndMDForLateral( subIndex, lateralIndex );
if ( lateralCoordMDPairs.empty() ) continue;
std::vector<cvf::Vec3d> lateralCoords;
std::vector<double> lateralMDs;
for ( auto& [coord, md] : lateralCoordMDPairs )
{
lateralCoords.push_back( coord );
lateralMDs.push_back( md );
}
auto lateralIntersections = RigWellPathIntersectionTools::findCellIntersectionInfosAlongPath( eclipseCase->eclipseCaseData(),
wellPath->name(),
lateralCoords,
lateralMDs );
if ( lateralIntersections.empty() ) continue;
const int latBranch = ++branchNumber;
double prevMD = lateralMDs.front();
double prevTVD = -lateralCoords.front().z();
int latOutletSeg = icdSegNum;
bool firstLatSeg = true;
const auto lateralLabel = QString( "Lateral %1" ).arg( lateralIndex + 1 ).toStdString();
std::vector<RigMswSegment> latSegs;
FishbonesLateral lateralContext;
for ( const auto& cellIntInfo : lateralIntersections )
{
if ( auto mci = toMswCellIntersection( cellIntInfo, mainGrid, cellIntInfo.startMD, cellIntInfo.endMD ) )
{
// Deduplicate COMPSEGS: skip cells already connected (matching tree intersectedCells logic)
const bool isNewCell = usedGlobalCellIndices.insert( cellIntInfo.globCellIndex ).second;
const double startMD = prevMD;
const double endMD = cellIntInfo.endMD;
const double startTVD = prevTVD;
const double endTVD = cellIntInfo.endTVD();
FishbonesLateralSegment lateralSegmentContext;
lateralSegmentContext.globalCellIndex = cellIntInfo.globCellIndex;
lateralSegmentContext.equivalentDiameter = subs->equivalentDiameter( unitSystem );
const auto subSegs = RicMswTableDataTools::createSubSegmentMDPairs( startMD, endMD, customSegmentIntervals );
bool firstSubSeg = true;
for ( const auto& [subStartMD, subEndMD] : subSegs )
{
// The lateral geometry is not part of the well path, so the TVD of a sub-segment
// cannot be interpolated along the well path. Interpolate linearly instead.
const double mdRange = endMD - startMD;
const double startWeight = mdRange > 0.0 ? ( subStartMD - startMD ) / mdRange : 0.0;
const double endWeight = mdRange > 0.0 ? ( subEndMD - startMD ) / mdRange : 1.0;
const double subStartTVD = startTVD * ( 1.0 - startWeight ) + endTVD * startWeight;
const double subEndTVD = startTVD * ( 1.0 - endWeight ) + endTVD * endWeight;
double len = 0.0;
double dep = 0.0;
if ( infoType == "INC" )
{
len = subEndMD - subStartMD;
dep = subEndTVD - subStartTVD;
}
else
{
len = subEndMD;
dep = subEndTVD;
}
RigMswSegment latSeg;
latSeg.segmentNumber = segmentNumber++;
latSeg.outletSegmentNumber = latOutletSeg;
latSeg.length = len;
latSeg.depth = dep;
latSeg.diameter = subs->equivalentDiameter( unitSystem );
latSeg.roughness = subs->openHoleRoughnessFactor( unitSystem );
latSeg.sourceWellName = wellPath->name().toStdString();
if ( firstLatSeg ) latSeg.description = lateralLabel;
// The cell is connected by the first sub-segment only
if ( firstSubSeg && isNewCell ) latSeg.intersections = { *mci };
lateralSegmentContext.segmentNumbers.push_back( latSeg.segmentNumber );
latOutletSeg = latSeg.segmentNumber;
firstLatSeg = false;
firstSubSeg = false;
latSegs.push_back( std::move( latSeg ) );
}
lateralContext.segments.push_back( std::move( lateralSegmentContext ) );
prevMD = cellIntInfo.endMD;
prevTVD = cellIntInfo.endTVD();
}
}
fishbonesContext.laterals.push_back( std::move( lateralContext ) );
result.push_back( RigMswBranch{ latBranch, std::nullopt, std::move( latSegs ), RigMswBranchSource::Fishbones } );
}
}
}
return result;
}
} // namespace RicMswBranchBuilder