From 301e05caa0124ac182ec91aee0f6587d4549d7ef Mon Sep 17 00:00:00 2001 From: Gaute Lindkvist Date: Tue, 25 Sep 2018 10:12:45 +0200 Subject: [PATCH] #3409 Add time dependent mat to frac dP option --- .../RicExportCompletionDataSettingsUi.cpp | 20 +- .../RicExportFractureCompletionsImpl.cpp | 103 ++++---- .../RicExportFractureCompletionsImpl.h | 10 +- ...ellPathExportCompletionDataFeatureImpl.cpp | 26 +- .../RimEclipseStatisticsCaseEvaluator.cpp | 4 +- .../RigFractureTransmissibilityEquations.cpp | 17 -- .../RigFractureTransmissibilityEquations.h | 6 - .../RigTransmissibilityCondenser.cpp | 235 +++++++++++------- .../RigTransmissibilityCondenser.h | 20 +- 9 files changed, 250 insertions(+), 191 deletions(-) diff --git a/ApplicationCode/Commands/CompletionExportCommands/RicExportCompletionDataSettingsUi.cpp b/ApplicationCode/Commands/CompletionExportCommands/RicExportCompletionDataSettingsUi.cpp index 31c9c64cb7..9f250dc278 100644 --- a/ApplicationCode/Commands/CompletionExportCommands/RicExportCompletionDataSettingsUi.cpp +++ b/ApplicationCode/Commands/CompletionExportCommands/RicExportCompletionDataSettingsUi.cpp @@ -55,11 +55,12 @@ namespace caf void RicExportCompletionDataSettingsUi::TransScalingType::setUp() { addItem(RicExportFractureCompletionsImpl::NO_SCALING, "NO_SCALING", "No scaling"); - addItem(RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_MAX_DP, "MATFRAC_DP_OVER_MAXDP", "Matrix to Fracture dP over max dP"); - addItem(RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_AVG_DP, "MATFRAC_DP_OVER_AVGDP", "Matrix to Fracture dP over avg dP"); - addItem(RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX, "MATFRAC_FLUX_OVER_MAXFLUX", "Matrix to Fracture Flux over max Flux"); - addItem(RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_FLUX_OVER_AVG_FLUX, "MATFRAC_FLUX_OVER_AVGFLUX", "Matrix to Fracture Flux over avg Flux"); + addItem(RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_INITIAL_DP, "MATFRAC_DP_OVER_INITIALDP", "Matrix to Fracture dP over initial dP"); + addItem(RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_MAX_INITIAL_DP, "MATFRAC_DP_OVER_MAX_INITIALDP", "Matrix to Fracture dP over max initial dP"); addItem(RicExportFractureCompletionsImpl::MATRIX_TO_WELL_DP_OVER_INITIAL_DP, "MATWELL_DP_OVER_INITIALDP", "Matrix to Well dP over initial dP"); + addItem(RicExportFractureCompletionsImpl::MATRIX_TO_WELL_DP_OVER_MAX_INITIAL_DP, "MATWELL_DP_OVER_MAX_INITIALDP", "Matrix to Well dP over max initial dP"); + addItem(RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX, "MATFRAC_FLUX_OVER_MAXFLUX", "Matrix to Fracture Flux over max Flux"); + setDefault(RicExportFractureCompletionsImpl::NO_SCALING); } @@ -211,17 +212,6 @@ void RicExportCompletionDataSettingsUi::fieldChangedByUi(const caf::PdmFieldHand includeFractures = true; } } - else if (changedField == &transScalingType) - { - if (transScalingType == RicExportFractureCompletionsImpl::MATRIX_TO_WELL_DP_OVER_INITIAL_DP) - { - transScalingCorrection = RicExportFractureCompletionsImpl::HOGSTOL_CORRECTION; - } - else - { - transScalingCorrection = RicExportFractureCompletionsImpl::NO_CORRECTION; - } - } } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp index ef1370dc5b..7d214fb4c8 100644 --- a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp +++ b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp @@ -325,16 +325,25 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat // Insert total transmissibility from eclipse-cell to well for this fracture into the map std::map matrixToWellTrans = calculateMatrixToWellTransmissibilities(transCondenser); - if (currentPressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_MAX_DP || - currentPressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_AVG_DP) + + //////////////////////////////////////////////////////// + // WARNING!!! // + // PROTOTYPE-CODE for Pressure Differential Depletion // + // MAY CHANGE A LOT // + //////////////////////////////////////////////////////// + if (currentPressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_INITIAL_DP || + currentPressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_MAX_INITIAL_DP) { RigTransmissibilityCondenser scaledCondenser = transCondenser; // 1. Scale matrix to fracture transmissibilities by matrix to fracture pressure - std::map originalLumpedMatrixToFractureTrans = scaledCondenser.scaleMatrixToFracTransByMatrixFracDP( - actCellInfo, - currentWellPressure, - *currentMatrixPressures, - currentPressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_AVG_DP); + std::map originalLumpedMatrixToFractureTrans = + scaledCondenser.scaleMatrixToFracTransByMatrixFracInitialDP(actCellInfo, + initialWellPressure, + currentWellPressure, + *initialMatrixPressures, + *currentMatrixPressures, + currentPressureDropScaling == + MATRIX_TO_FRACTURE_DP_OVER_MAX_INITIAL_DP); // 2: Calculate new external transmissibilities scaledCondenser.calculateCondensedTransmissibilities(); @@ -358,47 +367,18 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat matrixToWellTrans = effectiveMatrixToWellTrans; } } - else if (currentPressureDropScaling == MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX || - currentPressureDropScaling == MATRIX_TO_FRACTURE_FLUX_OVER_AVG_FLUX) + else if (currentPressureDropScaling == MATRIX_TO_WELL_DP_OVER_INITIAL_DP || + currentPressureDropScaling == MATRIX_TO_WELL_DP_OVER_MAX_INITIAL_DP) { RigTransmissibilityCondenser scaledCondenser = transCondenser; // 1. Scale matrix to fracture transmissibilities by matrix to fracture pressure - std::map originalLumpedMatrixToFractureTrans = scaledCondenser.scaleMatrixToFracTransByMatrixFracFlux( - actCellInfo, - currentWellPressure, - *currentMatrixPressures, - currentPressureDropScaling == MATRIX_TO_FRACTURE_FLUX_OVER_AVG_FLUX); - // 2: Calculate new external transmissibilities - scaledCondenser.calculateCondensedTransmissibilities(); - - if (pdParams.transCorrection == NO_CORRECTION) - { - // Calculate effective matrix to well transmissibilities. - std::map effectiveMatrixToWellTransBeforeCorrection = - calculateMatrixToWellTransmissibilities(scaledCondenser); - matrixToWellTrans = effectiveMatrixToWellTransBeforeCorrection; - } - else if (pdParams.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 (currentPressureDropScaling == 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.scaleMatrixToFracTransByMatrixWellDP( - actCellInfo, initialWellPressure, currentWellPressure, *initialMatrixPressures, *currentMatrixPressures); + actCellInfo, + initialWellPressure, + currentWellPressure, + *initialMatrixPressures, + *currentMatrixPressures, + currentPressureDropScaling == MATRIX_TO_FRACTURE_DP_OVER_MAX_INITIAL_DP); // 2: Calculate new external transmissibilities scaledCondenser.calculateCondensedTransmissibilities(); @@ -422,6 +402,41 @@ std::vector RicExportFractureCompletionsImpl::generateCompdat matrixToWellTrans = effectiveMatrixToWellTrans; } } + else if (currentPressureDropScaling == MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX) + { + RigTransmissibilityCondenser scaledCondenser = transCondenser; + // 1. Scale matrix to fracture transmissibilities by matrix to fracture pressure Depletion: + std::map originalLumpedMatrixToFractureTrans = + scaledCondenser.scaleMatrixToFracTransByMatrixFracFlux(actCellInfo, + currentWellPressure, + *currentMatrixPressures, + false); + // 2: Calculate new external transmissibilities + scaledCondenser.calculateCondensedTransmissibilities(); + + if (pdParams.transCorrection == NO_CORRECTION) + { + // Calculate effective matrix to well transmissibilities. + std::map effectiveMatrixToWellTransBeforeCorrection = + calculateMatrixToWellTransmissibilities(scaledCondenser); + matrixToWellTrans = effectiveMatrixToWellTransBeforeCorrection; + } + else if (pdParams.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; + } + } + //////////////////////////////////////////////////////////// + // END PROTOTYPE CODE FOR PRESSURE DIFFERENTIAL DEPLETION // + //////////////////////////////////////////////////////////// std::vector allCompletionsForOneFracture = generateCompdatValuesForFracture(matrixToWellTrans, wellPathName, caseToApply, fracture, fracTemplate); diff --git a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h index 7d755d36a3..115cba7ab6 100644 --- a/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h +++ b/ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h @@ -47,11 +47,11 @@ public: enum PressureDepletionTransScaling { NO_SCALING = 0, - MATRIX_TO_FRACTURE_DP_OVER_MAX_DP, - MATRIX_TO_FRACTURE_DP_OVER_AVG_DP, - MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX, - MATRIX_TO_FRACTURE_FLUX_OVER_AVG_FLUX, - MATRIX_TO_WELL_DP_OVER_INITIAL_DP + MATRIX_TO_FRACTURE_DP_OVER_INITIAL_DP, + MATRIX_TO_FRACTURE_DP_OVER_MAX_INITIAL_DP, + MATRIX_TO_WELL_DP_OVER_INITIAL_DP, + MATRIX_TO_WELL_DP_OVER_MAX_INITIAL_DP, + MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX }; enum PressureDepletionTransCorrection diff --git a/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp b/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp index ea6813aa2c..96bbd9e4d0 100644 --- a/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp +++ b/ApplicationCode/Commands/CompletionExportCommands/RicWellPathExportCompletionDataFeatureImpl.cpp @@ -2491,22 +2491,34 @@ QString RicWellPathExportCompletionDataFeatureImpl::createPressureDepletionFileN QString suffix; if (exportSettings.transScalingType() != RicExportFractureCompletionsImpl::NO_CORRECTION) { - if (exportSettings.transScalingType() == RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_AVG_DP) + if (exportSettings.transScalingType() == RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_INITIAL_DP) { - suffix += QString("_PAVG_"); + suffix += QString("_1"); } - else if (exportSettings.transScalingType() == RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_MAX_DP) + else if (exportSettings.transScalingType() == RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_DP_OVER_MAX_INITIAL_DP) { - suffix += QString("_PMAX_"); + suffix += QString("_2"); } - else + else if (exportSettings.transScalingType() == RicExportFractureCompletionsImpl::MATRIX_TO_WELL_DP_OVER_INITIAL_DP) { - suffix += QString("_PMWD_"); + suffix += QString("_3"); + } + else if (exportSettings.transScalingType() == RicExportFractureCompletionsImpl::MATRIX_TO_WELL_DP_OVER_MAX_INITIAL_DP) + { + suffix += QString("_4"); + } + else if (exportSettings.transScalingType() == RicExportFractureCompletionsImpl::MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX) + { + suffix += QString("_5"); } if (exportSettings.transScalingCorrection() == RicExportFractureCompletionsImpl::HOGSTOL_CORRECTION) { - suffix += QString("_HC_"); + suffix += QString("B_"); + } + else + { + suffix += QString("A_"); } if (exportSettings.transScalingSummaryWBHP()) diff --git a/ApplicationCode/ProjectDataModel/RimEclipseStatisticsCaseEvaluator.cpp b/ApplicationCode/ProjectDataModel/RimEclipseStatisticsCaseEvaluator.cpp index c4e0fbcc68..434fc62491 100644 --- a/ApplicationCode/ProjectDataModel/RimEclipseStatisticsCaseEvaluator.cpp +++ b/ApplicationCode/ProjectDataModel/RimEclipseStatisticsCaseEvaluator.cpp @@ -112,9 +112,9 @@ void RimEclipseStatisticsCaseEvaluator::evaluateForResults(const QList& if (activeCellCount > 0) { - for (size_t i = 0; i < statisticalResultNames.size(); ++i) + for (size_t j = 0; j < statisticalResultNames.size(); ++j) { - addNamedResult(destCellResultsData, resultType, statisticalResultNames[i], activeCellCount); + addNamedResult(destCellResultsData, resultType, statisticalResultNames[j], activeCellCount); } } } diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp index c10e03f720..6ad05545b0 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp +++ b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp @@ -122,23 +122,6 @@ double RigFractureTransmissibilityEquations::matrixToFractureTrans(double perm, return transmissibility; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RigFractureTransmissibilityEquations::pressureScalingMatrixToFractureTransPDDHC(double originalWellPressure, - double wellPressure, - double originalMatrixPressure, - double matrixPressure) -{ - double pressureDelta = originalMatrixPressure - originalWellPressure; - if (cvf::Math::abs(pressureDelta) > EPSILON) - { - return (matrixPressure - wellPressure) / pressureDelta; - } - CVF_ASSERT(false); - return 0.0; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h index 70172139ce..16741f7871 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h +++ b/ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h @@ -56,12 +56,6 @@ public: double fractureAreaWeightedlength, double cDarcy); - // Pressure Differential Depletion Høgstøl-correction (PDDHC) methods. - static double pressureScalingMatrixToFractureTransPDDHC(double originalWellPressure, - double wellPressure, - double originalMatrixPressure, - double matrixPressure); - static double effectiveInternalFractureToWellTransPDDHC(double sumScaledMatrixToFractureTrans, double scaledMatrixToWellTrans); diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp index d05891e2d5..8ff0d13ab8 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp +++ b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp @@ -123,17 +123,40 @@ double RigTransmissibilityCondenser::condensedTransmissibility(CellAddress exter //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixWellDP( - const RigActiveCellInfo* actCellInfo, - double originalWellPressure, - double currentWellPressure, - const std::vector& originalMatrixPressures, - const std::vector& currentMatrixPressures) +std::map + RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixWellDP(const RigActiveCellInfo* actCellInfo, + double initialWellPressure, + double currentWellPressure, + const std::vector& initialMatrixPressures, + const std::vector& currentMatrixPressures, + bool normalizeByMax) { - CVF_ASSERT(originalMatrixPressures.size() == currentMatrixPressures.size()); + CVF_ASSERT(initialMatrixPressures.size() == currentMatrixPressures.size()); std::map originalLumpedMatrixToFractureTrans; // Sum(T_mf) + double maxInitialDeltaPressure = 0.0; + if (normalizeByMax) + { + 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()); + + double initialDeltaPressure = initialMatrixPressures[eclipseResultIndex] - initialWellPressure; + maxInitialDeltaPressure = std::max(maxInitialDeltaPressure, initialDeltaPressure); + } + } + } + } + } for (auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it) { if (it->first.m_cellIndexSpace == CellAddress::STIMPLAN) @@ -143,16 +166,21 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM if (jt->first.m_cellIndexSpace == CellAddress::ECLIPSE) { size_t globalMatrixCellIdx = jt->first.m_globalCellIdx; - size_t eclipseResultIndex = actCellInfo->cellResultIndex(globalMatrixCellIdx); + size_t eclipseResultIndex = actCellInfo->cellResultIndex(globalMatrixCellIdx); CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); originalLumpedMatrixToFractureTrans[globalMatrixCellIdx] += jt->second; - jt->second *= RigFractureTransmissibilityEquations::pressureScalingMatrixToFractureTransPDDHC( - originalWellPressure, - currentWellPressure, - originalMatrixPressures[eclipseResultIndex], - currentMatrixPressures[eclipseResultIndex]); + double initialDeltaPressure = initialMatrixPressures[eclipseResultIndex] - initialWellPressure; + double currentDeltaPressure = currentMatrixPressures[eclipseResultIndex] - currentWellPressure; + if (normalizeByMax) + { + jt->second *= currentDeltaPressure / maxInitialDeltaPressure; + } + else + { + jt->second *= currentDeltaPressure / initialDeltaPressure; + } } } } @@ -160,90 +188,60 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM return originalLumpedMatrixToFractureTrans; } - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixFracDP(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector& currentMatrixPressures, bool divideByAverageDP) +std::map + RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixFracInitialDP(const RigActiveCellInfo* actCellInfo, + double initialWellPressure, + double currentWellPressure, + const std::vector& initialMatrixPressures, + const std::vector& currentMatrixPressures, + bool normalizeByMax) { - // 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); + // Solve for fracture pressures, current and initial + std::map initialFractureCellToPressureMap = + solveForFracturePressures(actCellInfo, initialMatrixPressures, initialWellPressure); + std::map currentFractureCellToPressureMap = + solveForFracturePressures(actCellInfo, currentMatrixPressures, currentWellPressure); - // Extract fracture pressures into a map - std::map fractureCellToPressureMap; + // Calculate maximum pressure drop + double maxInitialPressureDrop = 0.0; + if (normalizeByMax) { - size_t rowIndex = 0u; - for (const ConnectionTransmissibility& connectionTrans : m_neighborTransmissibilities) + for (auto it = m_neighborTransmissibilities.begin(); it != m_neighborTransmissibilities.end(); ++it) { - if (connectionTrans.first.m_cellIndexSpace == CellAddress::STIMPLAN) + if (it->first.m_cellIndexSpace == CellAddress::STIMPLAN) { - fractureCellToPressureMap[connectionTrans.first.m_globalCellIdx] = fracturePressures[rowIndex++]; - } - } - } + size_t globalFractureCellIdx = it->first.m_globalCellIdx; + double initialFracturePressure = initialFractureCellToPressureMap[globalFractureCellIdx]; - // 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 = fractureCellToPressureMap[globalFractureCellIdx]; - - for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) - { - if (jt->first.m_cellIndexSpace == CellAddress::ECLIPSE) + for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) { - size_t globalMatrixCellIdx = jt->first.m_globalCellIdx; + 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()); + size_t eclipseResultIndex = actCellInfo->cellResultIndex(globalMatrixCellIdx); + CVF_ASSERT(eclipseResultIndex < initialMatrixPressures.size()); - double matrixPressure = currentMatrixPressures[eclipseResultIndex]; - double pressureDrop = std::abs(matrixPressure - fracturePressure); - meanCalculator.addValueAndWeight(pressureDrop, 1.0); - maxPressureDrop = std::max(maxPressureDrop, pressureDrop); + double initialMatrixPressure = initialMatrixPressures[eclipseResultIndex]; + double initialPressureDrop = std::abs(initialMatrixPressure - initialFracturePressure); + maxInitialPressureDrop = std::max(maxInitialPressureDrop, initialPressureDrop); + } } } } } - 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 = fractureCellToPressureMap[globalFractureCellIdx]; + size_t globalFractureCellIdx = it->first.m_globalCellIdx; + double currentFracturePressure = currentFractureCellToPressureMap[globalFractureCellIdx]; + double initialFracturePressure = initialFractureCellToPressureMap[globalFractureCellIdx]; for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) { @@ -254,14 +252,24 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM size_t eclipseResultIndex = actCellInfo->cellResultIndex(globalMatrixCellIdx); CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); - double matrixPressure = currentMatrixPressures[eclipseResultIndex]; - double pressureDrop = std::abs(matrixPressure - fracturePressure); + double currentMatrixPressure = currentMatrixPressures[eclipseResultIndex]; + double pressureDrop = std::abs(currentMatrixPressure - currentFracturePressure); // Add to Sum(T_mf) originalLumpedMatrixToFractureTrans[globalMatrixCellIdx] += jt->second; - double pressureScaling = pressureDrop / (divideByAverageDP ? averagePressureDrop : maxPressureDrop); - jt->second *= pressureScaling; + if (normalizeByMax) + { + double pressureScaling = pressureDrop / maxInitialPressureDrop; + jt->second *= pressureScaling; + } + else + { + double initialMatrixPressure = initialMatrixPressures[eclipseResultIndex]; + double initialPressureDrop = std::abs(initialMatrixPressure - initialFracturePressure); + double pressureScaling = pressureDrop / initialPressureDrop; + jt->second *= pressureScaling; + } } } } @@ -272,7 +280,54 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixFracFlux(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector& currentMatrixPressures, bool divideByAverageFlux) +std::map + RigTransmissibilityCondenser::solveForFracturePressures(const RigActiveCellInfo* actCellInfo, + const std::vector& currentMatrixPressures, + double currentWellPressure) +{ + Eigen::VectorXd externalPressures(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()); + externalPressures[rowIndex++] = currentMatrixPressures[eclipseResultIndex]; + } + else + { + CVF_ASSERT(externalCell.m_cellIndexSpace == CellAddress::WELL); + externalPressures[rowIndex++] = currentWellPressure; + } + } + } + Eigen::VectorXd fracturePressures = m_TiiInv * (-m_Tie * externalPressures); + + // Extract fracture pressures into a map + std::map fractureCellToPressureMap; + { + size_t rowIndex = 0u; + for (const ConnectionTransmissibility& connectionTrans : m_neighborTransmissibilities) + { + if (connectionTrans.first.m_cellIndexSpace == CellAddress::STIMPLAN) + { + fractureCellToPressureMap[connectionTrans.first.m_globalCellIdx] = fracturePressures[rowIndex++]; + } + } + } + return fractureCellToPressureMap; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::map + RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixFracFlux(const RigActiveCellInfo* actCellInfo, + double currentWellPressure, + const std::vector& currentMatrixPressures, + bool divideByAverageFlux) { // Solve for fracture pressures Eigen::VectorXd matrixPressures(m_Tie.cols()); @@ -309,14 +364,14 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM } // Calculate maximum and average pressure drop - double maxFlux = 0.0; + double maxFlux = 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 = fractureCellToPressureMap[globalFractureCellIdx]; + double fracturePressure = fractureCellToPressureMap[globalFractureCellIdx]; for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) { @@ -328,8 +383,8 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); double matrixPressure = currentMatrixPressures[eclipseResultIndex]; - double pressureDrop = std::abs(matrixPressure - fracturePressure); - double flux = pressureDrop * jt->second; + double pressureDrop = std::abs(matrixPressure - fracturePressure); + double flux = pressureDrop * jt->second; meanCalculator.addValueAndWeight(flux, 1.0); maxFlux = std::max(maxFlux, flux); } @@ -352,7 +407,7 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM if (it->first.m_cellIndexSpace == CellAddress::STIMPLAN) { size_t globalFractureCellIdx = it->first.m_globalCellIdx; - double fracturePressure = fractureCellToPressureMap[globalFractureCellIdx]; + double fracturePressure = fractureCellToPressureMap[globalFractureCellIdx]; for (auto jt = it->second.begin(); jt != it->second.end(); ++jt) { @@ -364,8 +419,8 @@ std::map RigTransmissibilityCondenser::scaleMatrixToFracTransByM CVF_ASSERT(eclipseResultIndex < currentMatrixPressures.size()); double matrixPressure = currentMatrixPressures[eclipseResultIndex]; - double pressureDrop = std::abs(matrixPressure - fracturePressure); - double flux = jt->second * pressureDrop; + double pressureDrop = std::abs(matrixPressure - fracturePressure); + double flux = jt->second * pressureDrop; // Add to Sum(T_mf) originalLumpedMatrixToFractureTrans[globalMatrixCellIdx] += jt->second; @@ -411,10 +466,12 @@ std::map RigTransmissibilityCondenser::calculateFicticiousFractu // Sum(T'_mf) double scaledMatrixToFractureTrans = matrixToAllFracturesTrans[globalMatrixCellIdx]; // T'mw - double scaledMatrixToWellTrans = condensedTransmissibility(externalCell, { true, RigTransmissibilityCondenser::CellAddress::WELL, 1 }); + double scaledMatrixToWellTrans = + condensedTransmissibility(externalCell, {true, RigTransmissibilityCondenser::CellAddress::WELL, 1}); // T'_fjw fictitiousFractureToWellTrans[globalMatrixCellIdx] = - RigFractureTransmissibilityEquations::effectiveInternalFractureToWellTransPDDHC(scaledMatrixToFractureTrans, scaledMatrixToWellTrans); + RigFractureTransmissibilityEquations::effectiveInternalFractureToWellTransPDDHC(scaledMatrixToFractureTrans, + scaledMatrixToWellTrans); } } return fictitiousFractureToWellTrans; diff --git a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h index 6a30543ff5..21d5e6a8c9 100644 --- a/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h +++ b/ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h @@ -96,14 +96,22 @@ public: std::string condensedTransDebugOutput(const RigMainGrid* mainGrid, const RigFractureGrid* fractureGrid); std::map scaleMatrixToFracTransByMatrixWellDP(const RigActiveCellInfo* actCellInfo, - double originalWellPressure, - double currentWellPressure, - const std::vector& originalMatrixPressures, - const std::vector& currentMatrixPressures); - std::map scaleMatrixToFracTransByMatrixFracDP(const RigActiveCellInfo* actCellInfo, + double initialWellPressure, double currentWellPressure, + const std::vector& initialMatrixPressures, const std::vector& currentMatrixPressures, - bool divideByAverageDP); + bool normalizeByMax); + + std::map scaleMatrixToFracTransByMatrixFracInitialDP(const RigActiveCellInfo* actCellInfo, + double initialWellPressure, + double currentWellPressure, + const std::vector& initialMatrixPressures, + const std::vector& currentMatrixPressures, + bool normalizeByMax); + + std::map solveForFracturePressures(const RigActiveCellInfo* actCellInfo, const std::vector ¤tMatrixPressures, double currentWellPressure); + + std::map scaleMatrixToFracTransByMatrixFracFlux(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector& currentMatrixPressures,