OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

#3193 First implementation of pressure differential depletion scaling.

Browse Source
This commit is contained in:
Gaute Lindkvist
2018-09-14 11:29:43 +02:00
parent 9953c56c35
commit 934c4fffd8
7 changed files with 502 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.cpp
View File

@@ -23,6 +23,7 @@
#include <cmath>
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;
}

14
ApplicationCode/ReservoirDataModel/Completions/RigFractureTransmissibilityEquations.h
View File

@@ -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;
};

288
ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.cpp
View File

@@ -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 <Eigen/Core>
#include <Eigen/LU>
#include <iomanip>
#include <QDebug>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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::CellAddress> 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<size_t, double> RigTransmissibilityCondenser::scaleMatrixTransmissibilitiesByPressureMatrixWell(
const RigActiveCellInfo* actCellInfo,
double originalWellPressure,
double currentWellPressure,
const std::vector<double>& originalMatrixPressures,
const std::vector<double>& currentMatrixPressures)
{
CVF_ASSERT(originalMatrixPressures.size() == currentMatrixPressures.size());
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)
{
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<size_t, double> RigTransmissibilityCondenser::scaleMatrixTransmissibilitiesByPressureMatrixFracture(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector<double>& 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<size_t, double> 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<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 = 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<size_t, double>();
}
else if (!divideByAverageDP && maxPressureDrop < 1.0e-9)
{
return std::map<size_t, double>();
}
double averagePressureDrop = 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 = 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<size_t, double> RigTransmissibilityCondenser::calculateFicticiousFractureToWellTransmissibilities()
{
std::map<size_t, double> 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<size_t, double> 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<size_t, double> RigTransmissibilityCondenser::calculateEffectiveMatrixToWellTransmissibilities(
const std::map<size_t, double>& originalLumpedMatrixToFractureTrans,
const std::map<size_t, double>& ficticuousFractureToWellTransMap)
{
std::map<size_t, double> 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"

33
ApplicationCode/ReservoirDataModel/Completions/RigTransmissibilityCondenser.h
View File

@@ -20,10 +20,13 @@
#include "cafAssert.h"
#include <Eigen/Core>
#include <map>
#include <vector>
#include <set>
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<size_t, double> scaleMatrixTransmissibilitiesByPressureMatrixWell(const RigActiveCellInfo* actCellInfo, double originalWellPressure, double currentWellPressure, const std::vector<double>& originalMatrixPressures, const std::vector<double>& currentMatrixPressures);
std::map<size_t, double> scaleMatrixTransmissibilitiesByPressureMatrixFracture(const RigActiveCellInfo* actCellInfo, double currentWellPressure, const std::vector<double>& currentMatrixPressures, bool divideByAverageDP);
std::map<size_t, double> calculateFicticiousFractureToWellTransmissibilities();
std::map<size_t, double> calculateEffectiveMatrixToWellTransmissibilities(const std::map<size_t, double>& originalLumpedMatrixToFractureTrans,
const std::map<size_t, double>& ficticuousFractureToWellTransMap);
std::map<CellAddress, std::map<CellAddress, double> > m_neighborTransmissibilities;
std::map<CellAddress, std::map<CellAddress, double> > m_condensedTransmissibilities;
std::set<CellAddress> m_externalCellAddrSet;
};
void calculateCondensedTransmissibilities();
protected:
typedef std::pair<CellAddress, std::map<CellAddress, double>> ConnectionTransmissibility;
typedef std::map<CellAddress, std::map<CellAddress, double>> ConnectionTransmissibilities;
ConnectionTransmissibilities m_neighborTransmissibilities;
ConnectionTransmissibilities m_condensedTransmissibilities;
std::set<CellAddress> m_externalCellAddrSet;
Eigen::MatrixXd m_TiiInv;
Eigen::MatrixXd m_Tie;
};