#3394 Differential Depletion: Flux as a scaling factor for transmissibility.

ApplicationCode/Commands/CompletionExportCommands/RicExportCompletionDataSettingsUi.cpp

@@ -57,6 +57,8 @@ namespace caf
         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_WELL_DP_OVER_INITIAL_DP, "MATWELL_DP_OVER_INITIALDP", "Matrix to Well dP over initial dP");
         setDefault(RicExportFractureCompletionsImpl::NO_SCALING);
     }

ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp

@@ -350,6 +350,36 @@ std::vector<RigCompletionData> RicExportFractureCompletionsImpl::generateCompdat
                 matrixToWellTrans = effectiveMatrixToWellTrans;
             }
         }
+        else if (currentPressureDropScaling == MATRIX_TO_FRACTURE_FLUX_OVER_MAX_FLUX ||
+            currentPressureDropScaling == MATRIX_TO_FRACTURE_FLUX_OVER_AVG_FLUX)
+        {
+            RigTransmissibilityCondenser scaledCondenser = transCondenser;
+            // 1. Scale matrix to fracture transmissibilities by matrix to fracture pressure
+            std::map<size_t, double> 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<size_t, double> 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<size_t, double> fictitiousFractureToWellTransmissibilities = scaledCondenser.calculateFicticiousFractureToWellTransmissibilities();
+                // 2. Calculate new effective matrix to well transmissibilities
+                std::map<size_t, double> effectiveMatrixToWellTrans = scaledCondenser.calculateEffectiveMatrixToWellTransmissibilities(
+                    originalLumpedMatrixToFractureTrans, fictitiousFractureToWellTransmissibilities);
+                matrixToWellTrans = effectiveMatrixToWellTrans;
+            }
+        }
         else if (currentPressureDropScaling == MATRIX_TO_WELL_DP_OVER_INITIAL_DP)
         {
             RigTransmissibilityCondenser scaledCondenser = transCondenser;

ApplicationCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h

@@ -49,6 +49,8 @@ public:
         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
     };
 

ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp

@@ -269,6 +269,115 @@ std::map<size_t, double> RigTransmissibilityCondenser::scaleMatrixToFracTransByM
     return originalLumpedMatrixToFractureTrans;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::map<size_t, double> RigTransmissibilityCondenser::scaleMatrixToFracTransByMatrixFracFlux(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector<double>& currentMatrixPressures, bool divideByAverageFlux)
+{
+    // 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<size_t, double> 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++];
+            }
+        }
+    }
+
+    // Calculate maximum and average pressure drop
+    double maxFlux = 0.0;
+    RiaWeightedMeanCalculator<double> 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)
+                {
+                    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);
+                    double flux = pressureDrop * jt->second;
+                    meanCalculator.addValueAndWeight(flux, 1.0);
+                    maxFlux = std::max(maxFlux, flux);
+                }
+            }
+        }
+    }
+    if (divideByAverageFlux && !meanCalculator.validAggregatedWeight())
+    {
+        return std::map<size_t, double>();
+    }
+    else if (!divideByAverageFlux && maxFlux < 1.0e-9)
+    {
+        return std::map<size_t, double>();
+    }
+    double averageFlux = meanCalculator.weightedMean();
+
+    std::map<size_t, double> 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];
+
+            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);
+                    double flux = jt->second * pressureDrop;
+                    // Add to Sum(T_mf)
+                    originalLumpedMatrixToFractureTrans[globalMatrixCellIdx] += jt->second;
+
+                    double pressureScaling = flux / (divideByAverageFlux ? averageFlux : maxFlux);
+                    jt->second *= pressureScaling;
+                }
+            }
+        }
+    }
+    return originalLumpedMatrixToFractureTrans;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h

@@ -104,6 +104,11 @@ public:
                                                                   double                     currentWellPressure,
                                                                   const std::vector<double>& currentMatrixPressures,
                                                                   bool                       divideByAverageDP);
+    std::map<size_t, double> scaleMatrixToFracTransByMatrixFracFlux(const RigActiveCellInfo*   actCellInfo,
+                                                                    double                     currentWellPressure,
+                                                                    const std::vector<double>& currentMatrixPressures,
+                                                                    bool                       divideByAverageFlux);
+
     std::map<size_t, double> calculateFicticiousFractureToWellTransmissibilities();
     std::map<size_t, double>
         calculateEffectiveMatrixToWellTransmissibilities(const std::map<size_t, double>& originalLumpedMatrixToFractureTrans,