From 934c4fffd80a56b7f291a983fea5869c86fffa87 Mon Sep 17 00:00:00 2001 From: Gaute Lindkvist Date: Fri, 14 Sep 2018 11:29:43 +0200 Subject: [PATCH] #3193 First implementation of pressure differential depletion scaling. --- .../RicExportFractureCompletionsImpl.cpp | 101 +++++- .../RicExportFractureCompletionsImpl.h | 65 ++-- ...ellPathExportCompletionDataFeatureImpl.cpp | 8 +- .../RigFractureTransmissibilityEquations.cpp | 49 ++- .../RigFractureTransmissibilityEquations.h | 14 + .../RigTransmissibilityCondenser.cpp | 288 +++++++++++++++++- .../RigTransmissibilityCondenser.h | 33 +- 7 files changed, 502 insertions(+), 56 deletions(-) diff --git a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp index e0604c5f26..a6a37c0a75 100644 --- a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp +++ b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp @@ -61,7 +61,9 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat RimWellPath* wellPath, RimEclipseCase* caseToApply, std::vector* fractureDataForReport, - QTextStream* outputStreamForIntermediateResultsText) + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling, + PressureDepletionTransCorrection transCorrection) { std::vector fracturesAlongWellPath; @@ -77,7 +79,7 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat wellPath->wellPathGeometry(), fracturesAlongWellPath, fractureDataForReport, - outputStreamForIntermediateResultsText); + outputStreamForIntermediateResultsText, pressureDropScaling, transCorrection); } //-------------------------------------------------------------------------------------------------- @@ -86,7 +88,9 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat std::vector RicExportFractureCompletionsImpl::generateCompdatValuesForSimWell(RimEclipseCase* eclipseCase, const RimSimWellInView* well, - QTextStream* outputStreamForIntermediateResultsText) + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling, + PressureDepletionTransCorrection transCorrection) { std::vector completionData; @@ -104,7 +108,7 @@ std::vector } std::vector branchCompletions = generateCompdatValues( - eclipseCase, well->name(), branches[branchIndex], fractures, nullptr, outputStreamForIntermediateResultsText); + eclipseCase, well->name(), branches[branchIndex], fractures, nullptr, outputStreamForIntermediateResultsText, pressureDropScaling, transCorrection); completionData.insert(completionData.end(), branchCompletions.begin(), branchCompletions.end()); } @@ -121,7 +125,9 @@ std::vector const RigWellPath* wellPathGeometry, const std::vector& fractures, std::vector* fractureDataReportItems, - QTextStream* outputStreamForIntermediateResultsText) + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling, + PressureDepletionTransCorrection transCorrection) { std::vector fractureCompletions; @@ -166,8 +172,14 @@ std::vector caseToApply->loadStaticResultsByName(resultNames); } + if (pressureDropScaling != NO_SCALING) + { + RigCaseCellResultsData* results = caseToApply->results(RiaDefines::MATRIX_MODEL); + results->findOrLoadScalarResult(RiaDefines::DYNAMIC_NATIVE, "PRESSURE"); + } + return generateCompdatValuesConst( - caseToApply, wellPathName, wellPathGeometry, fractures, fractureDataReportItems, outputStreamForIntermediateResultsText); + caseToApply, wellPathName, wellPathGeometry, fractures, fractureDataReportItems, outputStreamForIntermediateResultsText, pressureDropScaling, transCorrection); } //-------------------------------------------------------------------------------------------------- @@ -179,7 +191,9 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat const RigWellPath* wellPathGeometry, const std::vector& fractures, std::vector* fractureDataReportItems, - QTextStream* outputStreamForIntermediateResultsText) + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling, + PressureDepletionTransCorrection transCorrection) { std::vector fractureCompletions; @@ -191,6 +205,17 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat double cDarcyInCorrectUnit = RiaEclipseUnitTools::darcysConstant(caseToApply->eclipseCaseData()->unitsType()); const RigMainGrid* mainGrid = caseToApply->eclipseCaseData()->mainGrid(); + const RigCaseCellResultsData* results = caseToApply->results(RiaDefines::MATRIX_MODEL); + size_t pressureResultIndex = results->findScalarResultIndex(RiaDefines::DYNAMIC_NATIVE, "PRESSURE"); + const RigActiveCellInfo* actCellInfo = caseToApply->eclipseCaseData()->activeCellInfo(RiaDefines::MATRIX_MODEL); + + const std::vector& originalMatrixPressures = results->cellScalarResults(pressureResultIndex).front(); + const std::vector& currentMatrixPressures = results->cellScalarResults(pressureResultIndex).back(); + + // TODO: extract well pressure + double originalWellPressure = 200; + double currentWellPressure = 200; + // To handle several fractures in the same eclipse cell we need to keep track of the transmissibility // to the well from each fracture intersecting the cell and sum these transmissibilities at the end. // std::map > @@ -239,7 +264,67 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat ///// // Insert total transmissibility from eclipse-cell to well for this fracture into the map - std::map matrixToWellTrans = calculateMatrixToWellTransmissibilities(transCondenser); + std::map matrixToWellTrans = calculateMatrixToWellTransmissibilities(transCondenser); + + if (pressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_MAX_DP || pressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_AVG_DP) + { + RigTransmissibilityCondenser scaledCondenser = transCondenser; + // 1. Scale matrix to fracture transmissibilities by matrix to fracture pressure + std::map originalLumpedMatrixToFractureTrans = + scaledCondenser.scaleMatrixTransmissibilitiesByPressureMatrixFracture(actCellInfo, + currentWellPressure, + currentMatrixPressures, + pressureDropScaling == + MATRIX_TO_FRACTURE_DP_OVER_AVG_DP); + // 2: Calculate new external transmissibilities + scaledCondenser.calculateCondensedTransmissibilities(); + + if (transCorrection == NO_CORRECTION) + { + // Calculate effective matrix to well transmissibilities. + std::map effectiveMatrixToWellTransBeforeCorrection = calculateMatrixToWellTransmissibilities(scaledCondenser); + matrixToWellTrans = effectiveMatrixToWellTransBeforeCorrection; + } + else if (transCorrection == HOGSTOL_CORRECTION) + { + // Høgstøl correction. + // 1. Calculate new effective fracture to well transmissiblities + std::map fictitiousFractureToWellTransmissibilities = scaledCondenser.calculateFicticiousFractureToWellTransmissibilities(); + // 2. Calculate new effective matrix to well transmissibilities + std::map effectiveMatrixToWellTrans = scaledCondenser.calculateEffectiveMatrixToWellTransmissibilities( + originalLumpedMatrixToFractureTrans, fictitiousFractureToWellTransmissibilities); + matrixToWellTrans = effectiveMatrixToWellTrans; + } + } + else if (pressureDropScaling == MATRIX_TO_WELL_DP_OVER_INITIAL_DP) + { + RigTransmissibilityCondenser scaledCondenser = transCondenser; + // From Høgstøl "Hydraulic Fracturing SoW 2.8 outside contract Fracture Transmissibility Calculations for Differential Depletion": + // 1. Scale matrix to fracture transmissibilities by matrix to well pressure + std::map originalLumpedMatrixToFractureTrans = + scaledCondenser.scaleMatrixTransmissibilitiesByPressureMatrixWell( + actCellInfo, originalWellPressure, currentWellPressure, originalMatrixPressures, currentMatrixPressures); + // 2: Calculate new external transmissibilities + scaledCondenser.calculateCondensedTransmissibilities(); + + if (transCorrection == NO_CORRECTION) + { + // Calculate effective matrix to well transmissibilities. + std::map effectiveMatrixToWellTransBeforeCorrection = calculateMatrixToWellTransmissibilities(scaledCondenser); + matrixToWellTrans = effectiveMatrixToWellTransBeforeCorrection; + } + else if (transCorrection == HOGSTOL_CORRECTION) + { + // Høgstøl correction. + // 1. Calculate new effective fracture to well transmissiblities + std::map fictitiousFractureToWellTransmissibilities = scaledCondenser.calculateFicticiousFractureToWellTransmissibilities(); + // 2. Calculate new effective matrix to well transmissibilities + std::map effectiveMatrixToWellTrans = scaledCondenser.calculateEffectiveMatrixToWellTransmissibilities( + originalLumpedMatrixToFractureTrans, fictitiousFractureToWellTransmissibilities); + matrixToWellTrans = effectiveMatrixToWellTrans; + } + } + std::vector allCompletionsForOneFracture = generateCompdatValuesForFracture(matrixToWellTrans, wellPathName, caseToApply, fracture, fracTemplate); diff --git a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h index 1392c5df8c..98b6266fbd 100644 --- a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h +++ b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h @@ -43,32 +43,51 @@ class QString; class RicExportFractureCompletionsImpl { public: - static std::vector - generateCompdatValuesForWellPath(RimWellPath* wellPath, - RimEclipseCase* caseToApply, - std::vector* fractureDataForReport, - QTextStream* outputStreamForIntermediateResultsText); + enum PressureDepletionTransScaling + { + NO_SCALING = 0, + MATRIX_TO_FRACTURE_DP_OVER_MAX_DP, + MATRIX_TO_FRACTURE_DP_OVER_AVG_DP, + MATRIX_TO_WELL_DP_OVER_INITIAL_DP + }; - static std::vector generateCompdatValuesForSimWell(RimEclipseCase* eclipseCase, + enum PressureDepletionTransCorrection + { + NO_CORRECTION, + HOGSTOL_CORRECTION + }; + + static std::vector generateCompdatValuesForWellPath(RimWellPath* wellPath, + RimEclipseCase* caseToApply, + std::vector* fractureDataForReport, + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling = MATRIX_TO_WELL_DP_OVER_INITIAL_DP, + PressureDepletionTransCorrection transCorrection = HOGSTOL_CORRECTION); + + static std::vector generateCompdatValuesForSimWell(RimEclipseCase* eclipseCase, const RimSimWellInView* well, - QTextStream* outputStreamForIntermediateResultsText); + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling = MATRIX_TO_WELL_DP_OVER_INITIAL_DP, + PressureDepletionTransCorrection transCorrection = HOGSTOL_CORRECTION); - static std::vector - generateCompdatValues(RimEclipseCase* caseToApply, - const QString& wellPathName, - const RigWellPath* wellPathGeometry, - const std::vector& fractures, - std::vector* fractureDataReportItems, - QTextStream* outputStreamForIntermediateResultsText); + static std::vector generateCompdatValues(RimEclipseCase* caseToApply, + const QString& wellPathName, + const RigWellPath* wellPathGeometry, + const std::vector& fractures, + std::vector* fractureDataReportItems, + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling, + PressureDepletionTransCorrection transCorrection); private: - static std::vector - generateCompdatValuesConst(const RimEclipseCase* caseToApply, - const QString& wellPathName, - const RigWellPath* wellPathGeometry, - const std::vector& fractures, - std::vector* fractureDataReportItems, - QTextStream* outputStreamForIntermediateResultsText); + static std::vector generateCompdatValuesConst(const RimEclipseCase* caseToApply, + const QString& wellPathName, + const RigWellPath* wellPathGeometry, + const std::vector& fractures, + std::vector* fractureDataReportItems, + QTextStream* outputStreamForIntermediateResultsText, + PressureDepletionTransScaling pressureDropScaling, + PressureDepletionTransCorrection transCorrection); static bool checkForStimPlanConductivity(const RimFractureTemplate* fracTemplate, const RimFracture* fracture); @@ -80,8 +99,8 @@ private: static void computeNonDarcyFlowParameters(const RimFracture* fracture, std::vector allCompletionsForOneFracture); - static double sumUpTransmissibilities(const std::vector& allCompletionsForOneFracture); - + static double sumUpTransmissibilities(const std::vector& allCompletionsForOneFracture); + static void calculateAndSetReportItemData(const std::vector& allCompletionsForOneFracture, const RigEclipseToStimPlanCalculator& calculator, RicWellPathFractureReportItem& reportItem); static void outputIntermediateResultsText(QTextStream* outputStreamForIntermediateResultsText, const RimFracture* fracture, RigTransmissibilityCondenser &transCondenser, const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid); diff --git a/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp b/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp index aba0fbc4cb..fe99f3f7ad 100644 --- a/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp +++ b/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp @@ -225,7 +225,9 @@ void RicWellPathExportCompletionDataFeatureImpl::exportCompletions(const std::ve wellPath, exportSettings.caseToApply(), reportItems, - fractureTransmissibilityExportInformationStream.get()); + fractureTransmissibilityExportInformationStream.get(), + exportSettings.transScalingType(), + exportSettings.transScalingCorrection()); appendCompletionData(&completionsPerEclipseCellAllCompletionTypes, fractureCompletionData); appendCompletionData(&completionsPerEclipseCellFracture, fractureCompletionData); @@ -1804,7 +1806,9 @@ RicMswExportInfo wellPath->wellPathGeometry(), { fracture }, nullptr, - nullptr); + nullptr, + RicExportFractureCompletionsImpl::NO_SCALING, + RicExportFractureCompletionsImpl::NO_CORRECTION); assignFractureIntersections(caseToApply, fracture, completionData, &location, &foundSubGridIntersections); } diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp index fd7aa599a6..4fc22fe255 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp +++ b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp @@ -23,6 +23,7 @@ #include +const double RigFractureTransmissibilityEquations::EPSILON = 1.0e-9; //-------------------------------------------------------------------------------------------------- /// @@ -44,7 +45,7 @@ double RigFractureTransmissibilityEquations::centerToCenterFractureCellTrans(dou } //-------------------------------------------------------------------------------------------------- -/// +/// T_16 //-------------------------------------------------------------------------------------------------- double RigFractureTransmissibilityEquations::fractureCellToWellRadialTrans(double fractureCellConductivity, double fractureCellSizeX, @@ -69,7 +70,7 @@ double RigFractureTransmissibilityEquations::fractureCellToWellRadialTrans(doubl } //-------------------------------------------------------------------------------------------------- -/// +/// T_16 //-------------------------------------------------------------------------------------------------- double RigFractureTransmissibilityEquations::fractureCellToWellLinearTrans(double fractureCellConductivity, double fractureCellSizeX, @@ -97,7 +98,7 @@ double RigFractureTransmissibilityEquations::fractureCellToWellLinearTrans(doubl //-------------------------------------------------------------------------------------------------- -/// +/// T_51 //-------------------------------------------------------------------------------------------------- double RigFractureTransmissibilityEquations::matrixToFractureTrans(double perm, double NTG, @@ -110,7 +111,7 @@ double RigFractureTransmissibilityEquations::matrixToFractureTrans(double perm, double transmissibility; double slDivPi = 0.0; - if ( cvf::Math::abs(skinfactor) > 1e-9) + if ( cvf::Math::abs(skinfactor) > EPSILON) { slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D; } @@ -121,6 +122,45 @@ double RigFractureTransmissibilityEquations::matrixToFractureTrans(double perm, return transmissibility; } +//-------------------------------------------------------------------------------------------------- +/// T'_51 +//-------------------------------------------------------------------------------------------------- +double RigFractureTransmissibilityEquations::pressureScalingMatrixToFractureTrans(double originalWellPressure, double wellPressure, double originalMatrixPressure, double matrixPressure) +{ + double pressureDelta = originalMatrixPressure - originalWellPressure; + if (cvf::Math::abs(pressureDelta) > EPSILON) + { + return (matrixPressure - wellPressure) / pressureDelta; + } + return 0.0; +} + +//-------------------------------------------------------------------------------------------------- +/// T'_16 +//-------------------------------------------------------------------------------------------------- +double RigFractureTransmissibilityEquations::effectiveInternalFractureToWellTrans(double scaledMatrixToFractureTrans, double scaledMatrixToWellTrans) +{ + double divisor = scaledMatrixToFractureTrans - scaledMatrixToWellTrans; + if (cvf::Math::abs(divisor) > EPSILON) + { + return (scaledMatrixToFractureTrans * scaledMatrixToWellTrans) / divisor; + } + return 0.0; +} + +//-------------------------------------------------------------------------------------------------- +///T^dp_56 +//-------------------------------------------------------------------------------------------------- +double RigFractureTransmissibilityEquations::effectiveMatrixToWellTrans(double originalMatrixToFractureTrans, double effectiveInternalFractureToWellTrans) +{ + double divisor = originalMatrixToFractureTrans + effectiveInternalFractureToWellTrans; + if (cvf::Math::abs(divisor) > EPSILON) + { + return (originalMatrixToFractureTrans * effectiveInternalFractureToWellTrans) / divisor; + } + return 0.0; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -133,3 +173,4 @@ double RigFractureTransmissibilityEquations::centerToEdgeFractureCellTrans(doubl return transmissibility; } + diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h index 48073997e0..93655f8664 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h +++ b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h @@ -57,11 +57,25 @@ public: double fractureAreaWeightedlength, double cDarcy); + static double pressureScalingMatrixToFractureTrans(double originalWellPressure, + double wellPressure, + double originalMatrixPressure, + double matrixPressure); + + + static double effectiveInternalFractureToWellTrans(double scaledMatrixToFractureTrans, + double scaledMatrixToWellTrans); + + static double effectiveMatrixToWellTrans(double originalMatrixToFractureTrans, + double effectiveInternalFractureToWellTrans); + private: static double centerToEdgeFractureCellTrans(double conductivity, double sideLengthParallellTrans, double sideLengthNormalTrans, double cDarcyForRelevantUnit); +private: + static const double EPSILON; }; diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp index 20aee6a618..ce116c3d95 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp +++ b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp @@ -18,12 +18,54 @@ #include "RigTransmissibilityCondenser.h" +#include "RigActiveCellInfo.h" +#include "RigFractureTransmissibilityEquations.h" #include "RiaLogging.h" +#include "RiaWeightedMeanCalculator.h" + +#include "cvfAssert.h" +#include "cvfBase.h" +#include "cvfMath.h" #include #include #include +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigTransmissibilityCondenser::RigTransmissibilityCondenser() +{ + +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigTransmissibilityCondenser::RigTransmissibilityCondenser(const RigTransmissibilityCondenser& copyFrom) + : m_neighborTransmissibilities(copyFrom.m_neighborTransmissibilities) + , m_condensedTransmissibilities(copyFrom.m_condensedTransmissibilities) + , m_externalCellAddrSet(copyFrom.m_externalCellAddrSet) + , m_TiiInv(copyFrom.m_TiiInv) + , m_Tie(copyFrom.m_Tie) +{ +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigTransmissibilityCondenser& RigTransmissibilityCondenser::operator=(const RigTransmissibilityCondenser& rhs) +{ + m_neighborTransmissibilities = rhs.m_neighborTransmissibilities; + m_condensedTransmissibilities = rhs.m_condensedTransmissibilities; + m_externalCellAddrSet = rhs.m_externalCellAddrSet; + m_TiiInv = rhs.m_TiiInv; + m_Tie = rhs.m_Tie; + return *this; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -44,7 +86,10 @@ void RigTransmissibilityCondenser::addNeighborTransmissibility(CellAddress cell1 //-------------------------------------------------------------------------------------------------- std::set RigTransmissibilityCondenser::externalCells() { - calculateCondensedTransmissibilitiesIfNeeded(); + if (m_externalCellAddrSet.empty()) + { + calculateCondensedTransmissibilities(); + } return m_externalCellAddrSet; } @@ -56,7 +101,10 @@ double RigTransmissibilityCondenser::condensedTransmissibility(CellAddress exter { CAF_ASSERT(!(externalCell1 == externalCell2)); - calculateCondensedTransmissibilitiesIfNeeded(); + if (m_condensedTransmissibilities.empty()) + { + calculateCondensedTransmissibilities(); + } if ( externalCell2 < externalCell1 ) std::swap(externalCell1, externalCell2); @@ -72,13 +120,228 @@ double RigTransmissibilityCondenser::condensedTransmissibility(CellAddress exter return 0.0; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::map RigTransmissibilityCondenser::scaleMatrixTransmissibilitiesByPressureMatrixWell( + const RigActiveCellInfo* actCellInfo, + double originalWellPressure, + double currentWellPressure, + const std::vector& originalMatrixPressures, + const std::vector& currentMatrixPressures) +{ + CVF_ASSERT(originalMatrixPressures.size() == currentMatrixPressures.size()); + + std::map originalLumpedMatrixToFractureTrans; // Sum(T_mf) + + for (auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it) + { + if (it->first.m_cellIndexSpace == CellAddress::STIMPLAN) + { + for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) + { + if (jt->first.m_cellIndexSpace == CellAddress::ECLIPSE) + { + size_t globalMatrixCellIdx = jt->first.m_globalCellIdx; + size_t eclipseResultIndex = actCellInfo->cellResultIndex(globalMatrixCellIdx); + CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); + + originalLumpedMatrixToFractureTrans[globalMatrixCellIdx] += jt->second; + + jt->second *= RigFractureTransmissibilityEquations::pressureScalingMatrixToFractureTrans( + originalWellPressure, + currentWellPressure, + originalMatrixPressures[eclipseResultIndex], + currentMatrixPressures[eclipseResultIndex]); + } + } + } + } + return originalLumpedMatrixToFractureTrans; +} + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::map RigTransmissibilityCondenser::scaleMatrixTransmissibilitiesByPressureMatrixFracture(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector& currentMatrixPressures, bool divideByAverageDP) +{ + // Solve for fracture pressures + Eigen::VectorXd matrixPressures(m_Tie.cols()); + { + size_t rowIndex = 0u; + for (const CellAddress& externalCell : m_externalCellAddrSet) + { + if (externalCell.m_cellIndexSpace == CellAddress::ECLIPSE) + { + size_t eclipseResultIndex = actCellInfo->cellResultIndex(externalCell.m_globalCellIdx); + CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); + matrixPressures[rowIndex++] = currentMatrixPressures[eclipseResultIndex]; + } + else + { + CVF_ASSERT(externalCell.m_cellIndexSpace == CellAddress::WELL); + matrixPressures[rowIndex++] = currentWellPressure; + } + } + } + Eigen::VectorXd fracturePressures = m_TiiInv * (m_Tie * matrixPressures * -1.0); + + // Extract fracture pressures into a map + std::map fracturePressureMap; + { + size_t rowIndex = 0u; + for (const ConnectionTransmissibility& connectionTrans : m_neighborTransmissibilities) + { + if (connectionTrans.first.m_cellIndexSpace == CellAddress::STIMPLAN) + { + fracturePressureMap[connectionTrans.first.m_globalCellIdx] = fracturePressures[rowIndex++]; + } + } + } + + // Calculate maximum and average pressure drop + double maxPressureDrop = 0.0; + RiaWeightedMeanCalculator meanCalculator; + for (auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it) + { + if (it->first.m_cellIndexSpace == CellAddress::STIMPLAN) + { + size_t globalFractureCellIdx = it->first.m_globalCellIdx; + double fracturePressure = fracturePressureMap[globalFractureCellIdx]; + + for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) + { + if (jt->first.m_cellIndexSpace == CellAddress::ECLIPSE) + { + size_t globalMatrixCellIdx = jt->first.m_globalCellIdx; + size_t eclipseResultIndex = actCellInfo->cellResultIndex(globalMatrixCellIdx); + CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); + double matrixPressure = currentMatrixPressures[eclipseResultIndex]; + double pressureDrop = std::abs(matrixPressure - fracturePressure); + meanCalculator.addValueAndWeight(pressureDrop, 1.0); + maxPressureDrop = std::max(maxPressureDrop, pressureDrop); + } + } + } + } + if (divideByAverageDP && !meanCalculator.validAggregatedWeight()) + { + return std::map(); + } + else if (!divideByAverageDP && maxPressureDrop < 1.0e-9) + { + return std::map(); + } + double averagePressureDrop = meanCalculator.weightedMean(); + + std::map originalLumpedMatrixToFractureTrans; // Sum(T_mf) + for (auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it) + { + if (it->first.m_cellIndexSpace == CellAddress::STIMPLAN) + { + size_t globalFractureCellIdx = it->first.m_globalCellIdx; + double fracturePressure = fracturePressureMap[globalFractureCellIdx]; + + for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) + { + if (jt->first.m_cellIndexSpace == CellAddress::ECLIPSE) + { + size_t globalMatrixCellIdx = jt->first.m_globalCellIdx; + size_t eclipseResultIndex = actCellInfo->cellResultIndex(globalMatrixCellIdx); + CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); + double matrixPressure = currentMatrixPressures[eclipseResultIndex]; + double pressureDrop = std::abs(matrixPressure - fracturePressure); + + // Add to Sum(T_mf) + originalLumpedMatrixToFractureTrans[globalMatrixCellIdx] += jt->second; + + double pressureScaling = pressureDrop / (divideByAverageDP ? averagePressureDrop : maxPressureDrop); + jt->second *= pressureScaling; + } + } + } + } + return originalLumpedMatrixToFractureTrans; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::map RigTransmissibilityCondenser::calculateFicticiousFractureToWellTransmissibilities() +{ + std::map matrixToAllFracturesTrans; + for (auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it) + { + if (it->first.m_cellIndexSpace == CellAddress::STIMPLAN) + { + for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) + { + if (jt->first.m_cellIndexSpace == CellAddress::ECLIPSE) + { + size_t globalMatrixCellIdx = jt->first.m_globalCellIdx; + // T'_mf + double matrixToFractureTrans = jt->second; + // Sum(T'_mf) + matrixToAllFracturesTrans[globalMatrixCellIdx] += matrixToFractureTrans; + } + } + } + } + + std::map fictitiousFractureToWellTrans; // T'_fjw + for (const CellAddress& externalCell : m_externalCellAddrSet) + { + if (externalCell.m_cellIndexSpace == CellAddress::ECLIPSE) + { + size_t globalMatrixCellIdx = externalCell.m_globalCellIdx; + // Sum(T'_mf) + double scaledMatrixToFractureTrans = matrixToAllFracturesTrans[globalMatrixCellIdx]; + // T'mw + double scaledMatrixToWellTrans = condensedTransmissibility(externalCell, { true, RigTransmissibilityCondenser::CellAddress::WELL, 1 }); + // T'_fjw + fictitiousFractureToWellTrans[globalMatrixCellIdx] = + RigFractureTransmissibilityEquations::effectiveInternalFractureToWellTrans(scaledMatrixToFractureTrans, scaledMatrixToWellTrans); + } + } + return fictitiousFractureToWellTrans; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::map RigTransmissibilityCondenser::calculateEffectiveMatrixToWellTransmissibilities( + const std::map& originalLumpedMatrixToFractureTrans, + const std::map& ficticuousFractureToWellTransMap) +{ + std::map effectiveMatrixToWellTrans; + for (const CellAddress& externalCell : m_externalCellAddrSet) + { + if (externalCell.m_cellIndexSpace == CellAddress::ECLIPSE) + { + size_t globalMatrixCellIdx = externalCell.m_globalCellIdx; + + auto matrixToFractureIt = originalLumpedMatrixToFractureTrans.find(globalMatrixCellIdx); + CVF_ASSERT(matrixToFractureIt != originalLumpedMatrixToFractureTrans.end()); + // Sum(T_mf) + double lumpedOriginalMatrixToFractureT = matrixToFractureIt->second; + // T'_fjw + auto fictitiousFractureToWellIt = ficticuousFractureToWellTransMap.find(globalMatrixCellIdx); + CVF_ASSERT(fictitiousFractureToWellIt != ficticuousFractureToWellTransMap.end()); + double fictitiousFractureToWellTrans = fictitiousFractureToWellIt->second; + // T^dp_mw + effectiveMatrixToWellTrans[globalMatrixCellIdx] = + RigFractureTransmissibilityEquations::effectiveMatrixToWellTrans(lumpedOriginalMatrixToFractureT, fictitiousFractureToWellTrans); + } + } + return effectiveMatrixToWellTrans; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -void RigTransmissibilityCondenser::calculateCondensedTransmissibilitiesIfNeeded() +void RigTransmissibilityCondenser::calculateCondensedTransmissibilities() { - if (m_condensedTransmissibilities.size()) return; - if (m_neighborTransmissibilities.empty()) return; // Find all equations, and their total ordering @@ -147,10 +410,13 @@ void RigTransmissibilityCondenser::calculateCondensedTransmissibilitiesIfNeeded( // std::cout << "T = " << std::endl << totalSystem << std::endl; int externalEquationCount = totalEquationCount - internalEquationCount; - MatrixXd condensedSystem = totalSystem.bottomRightCorner(externalEquationCount, externalEquationCount) - - totalSystem.bottomLeftCorner(externalEquationCount, internalEquationCount) - * totalSystem.topLeftCorner(internalEquationCount, internalEquationCount).inverse() - * totalSystem.topRightCorner(internalEquationCount, externalEquationCount ); + + MatrixXd Tee = totalSystem.bottomRightCorner(externalEquationCount, externalEquationCount); + MatrixXd Tei = totalSystem.bottomLeftCorner(externalEquationCount, internalEquationCount); + m_TiiInv = totalSystem.topLeftCorner(internalEquationCount, internalEquationCount).inverse(); + m_Tie = totalSystem.topRightCorner(internalEquationCount, externalEquationCount); + + MatrixXd condensedSystem = Tee - Tei * m_TiiInv * m_Tie; // std::cout << "Te = " << std::endl << condensedSystem << std::endl << std::endl; @@ -174,10 +440,6 @@ void RigTransmissibilityCondenser::calculateCondensedTransmissibilitiesIfNeeded( } } - - - - #include "RimStimPlanFractureTemplate.h" #include "RigMainGrid.h" #include "RigFractureCell.h" diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h index 7664b5f87c..f9650f8907 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h +++ b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h @@ -20,10 +20,13 @@ #include "cafAssert.h" +#include + #include #include #include +class RigActiveCellInfo; class RigMainGrid; class RimStimPlanFractureTemplate; class RigFractureGrid; @@ -34,6 +37,10 @@ class RigFractureGrid; class RigTransmissibilityCondenser { public: + RigTransmissibilityCondenser(); + RigTransmissibilityCondenser(const RigTransmissibilityCondenser& copyFrom); + RigTransmissibilityCondenser& operator=(const RigTransmissibilityCondenser& rhs); + class CellAddress { public: @@ -81,10 +88,24 @@ public: std::string neighborTransDebugOutput(const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid); std::string condensedTransDebugOutput(const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid); -private: - void calculateCondensedTransmissibilitiesIfNeeded(); + std::map scaleMatrixTransmissibilitiesByPressureMatrixWell(const RigActiveCellInfo* actCellInfo, double originalWellPressure, double currentWellPressure, const std::vector& originalMatrixPressures, const std::vector& currentMatrixPressures); + std::map scaleMatrixTransmissibilitiesByPressureMatrixFracture(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector& currentMatrixPressures, bool divideByAverageDP); + std::map calculateFicticiousFractureToWellTransmissibilities(); + std::map calculateEffectiveMatrixToWellTransmissibilities(const std::map& originalLumpedMatrixToFractureTrans, + const std::map& ficticuousFractureToWellTransMap); - std::map > m_neighborTransmissibilities; - std::map > m_condensedTransmissibilities; - std::set m_externalCellAddrSet; -}; + void calculateCondensedTransmissibilities(); + +protected: + + typedef std::pair> ConnectionTransmissibility; + typedef std::map> ConnectionTransmissibilities; + + ConnectionTransmissibilities m_neighborTransmissibilities; + ConnectionTransmissibilities m_condensedTransmissibilities; + + std::set m_externalCellAddrSet; + + Eigen::MatrixXd m_TiiInv; + Eigen::MatrixXd m_Tie; +};