mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
#1487 - pre-proto - Moving RigStimPlanUpscallingCalc to separate file
This commit is contained in:
astridkbjorke
@@ -45,6 +45,7 @@
|
||||
#include <QFile>
|
||||
#include <QString>
|
||||
#include <QTextStream>
|
||||
#include "RigStimPlanUpscalingCalc.h"
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -60,6 +60,8 @@ ${CEE_CURRENT_LIST_DIR}RigTesselatorTools.h
|
||||
${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.h
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.h
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanFracTemplateCell.h
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanUpscalingCalc.h
|
||||
|
||||
|
||||
|
||||
)
|
||||
@@ -116,6 +118,7 @@ ${CEE_CURRENT_LIST_DIR}RigTesselatorTools.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanFracTemplateCell.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanUpscalingCalc.cpp
|
||||
|
||||
)
|
||||
|
||||
|
||||
@@ -278,312 +278,6 @@ bool RigFractureTransCalc::planeCellIntersectionPolygons(cvf::Vec3d hexCorners[8
|
||||
return isCellIntersected;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::pair<double, double> RigStimPlanUpscalingCalc::flowAcrossLayersUpscaling(QString resultName, QString resultUnit, size_t timeStepIndex, RimDefines::UnitSystem unitSystem, size_t eclipseCellIndex)
|
||||
{
|
||||
RimStimPlanFractureTemplate* fracTemplateStimPlan;
|
||||
if (dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition()))
|
||||
{
|
||||
fracTemplateStimPlan = dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition());
|
||||
}
|
||||
else return std::make_pair(cvf::UNDEFINED_DOUBLE, cvf::UNDEFINED_DOUBLE);
|
||||
|
||||
std::vector<RigStimPlanFracTemplateCell> stimPlanCells = fracTemplateStimPlan->getStimPlanCells();
|
||||
|
||||
|
||||
|
||||
cvf::Vec3d localX, localY, localZ; //Not used in calculation here, but needed for function to find planCellPolygons
|
||||
std::vector<std::vector<cvf::Vec3d> > planeCellPolygons;
|
||||
|
||||
const RigMainGrid* mainGrid = m_case->eclipseCaseData()->mainGrid();
|
||||
cvf::Vec3d hexCorners[8];
|
||||
mainGrid->cellCornerVertices(eclipseCellIndex, hexCorners);
|
||||
|
||||
|
||||
|
||||
bool isPlanIntersected = RigFractureTransCalc::planeCellIntersectionPolygons(hexCorners, m_fracture->transformMatrix(), planeCellPolygons);
|
||||
if (!isPlanIntersected || planeCellPolygons.size() == 0) return std::make_pair(cvf::UNDEFINED_DOUBLE, cvf::UNDEFINED_DOUBLE);
|
||||
|
||||
//Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon/stimPlan mesh already is located)
|
||||
cvf::Mat4f invertedTransMatrix = m_fracture->transformMatrix().getInverted();
|
||||
for (std::vector<cvf::Vec3d> & planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
for (cvf::Vec3d& v : planeCellPolygon)
|
||||
{
|
||||
v.transformPoint(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
}
|
||||
}
|
||||
|
||||
cvf::Vec3d directionAcrossLayers;
|
||||
cvf::Vec3d directionAlongLayers;
|
||||
directionAcrossLayers = cvf::Vec3d(0.0, -1.0, 0.0);
|
||||
directionAlongLayers = cvf::Vec3d(1.0, 0.0, 0.0);
|
||||
|
||||
std::vector<cvf::Vec3f> fracPolygon = m_fracture->attachedFractureDefinition()->fracturePolygon(unitSystem);
|
||||
std::vector<std::vector<cvf::Vec3d> > polygonsDescribingFractureInCell;
|
||||
|
||||
std::vector<double> upscaledConductivitiesHA;
|
||||
std::vector<double> upscaledConductivitiesAH;
|
||||
|
||||
for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
double condHA = computeHAupscale(fracTemplateStimPlan, stimPlanCells, planeCellPolygon, directionAlongLayers, directionAcrossLayers);
|
||||
upscaledConductivitiesHA.push_back(condHA);
|
||||
|
||||
double condAH = computeAHupscale(fracTemplateStimPlan, stimPlanCells, planeCellPolygon, directionAlongLayers, directionAcrossLayers);
|
||||
upscaledConductivitiesAH.push_back(condAH);
|
||||
}
|
||||
|
||||
return std::make_pair(arithmeticAverage(upscaledConductivitiesHA), arithmeticAverage(upscaledConductivitiesAH));
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigStimPlanUpscalingCalc::computeHAupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers)
|
||||
{
|
||||
std::vector<double> DcolSum;
|
||||
std::vector<double> lavgCol;
|
||||
std::vector<double> CondHarmCol;
|
||||
|
||||
for (size_t j = 0; j < fracTemplateStimPlan->stimPlanGridNumberOfColums(); j++)
|
||||
{
|
||||
std::vector<double> conductivitiesInStimPlanCells;
|
||||
std::vector<double> lengthsLiOfStimPlanCol;
|
||||
std::vector<double> heightsDiOfStimPlanCells;
|
||||
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellsCol = getColOfStimPlanCells(stimPlanCells, j);
|
||||
for (RigStimPlanFracTemplateCell* stimPlanCell : stimPlanCellsCol)
|
||||
{
|
||||
if (stimPlanCell->getConductivtyValue() > 1e-7)
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
|
||||
if (clippedStimPlanPolygons.size() > 0)
|
||||
{
|
||||
for (auto clippedStimPlanPolygon : clippedStimPlanPolygons)
|
||||
{
|
||||
conductivitiesInStimPlanCells.push_back(stimPlanCell->getConductivtyValue());
|
||||
lengthsLiOfStimPlanCol.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAlongLayers, clippedStimPlanPolygon));
|
||||
heightsDiOfStimPlanCells.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAcrossLayers, clippedStimPlanPolygon));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
//Regne ut average
|
||||
double sumDiDivCondLi = 0.0;
|
||||
double sumDi = 0.0;
|
||||
double sumLiDi = 0.0;
|
||||
for (int i = 0; i < conductivitiesInStimPlanCells.size(); i++)
|
||||
{
|
||||
sumDiDivCondLi += heightsDiOfStimPlanCells[i] / (conductivitiesInStimPlanCells[i] * lengthsLiOfStimPlanCol[i]);
|
||||
sumDi += heightsDiOfStimPlanCells[i];
|
||||
sumLiDi += heightsDiOfStimPlanCells[i] * lengthsLiOfStimPlanCol[i];
|
||||
}
|
||||
|
||||
if (sumDiDivCondLi != 0)
|
||||
{
|
||||
DcolSum.push_back(sumDi);
|
||||
double lAvgValue = sumLiDi / sumDi;
|
||||
lavgCol.push_back(lAvgValue);
|
||||
double harmMeanForCol = (sumDi / lAvgValue) * (1 / sumDiDivCondLi);
|
||||
CondHarmCol.push_back(harmMeanForCol);
|
||||
}
|
||||
}
|
||||
|
||||
//Do arithmetic upscaling based on harmonic upscaled values for coloums
|
||||
double sumCondHLiDivDi = 0.0;
|
||||
double sumLi = 0.0;
|
||||
double sumDiLi = 0.0;
|
||||
for (int i = 0; i < CondHarmCol.size(); i++)
|
||||
{
|
||||
sumLi += lavgCol[i];
|
||||
sumDiLi += DcolSum[i] * lavgCol[i];
|
||||
sumCondHLiDivDi += CondHarmCol[i] * lavgCol[i] / DcolSum[i];
|
||||
}
|
||||
double Davg = sumDiLi / sumLi;
|
||||
double condHA = (Davg / sumLi) * sumCondHLiDivDi;
|
||||
|
||||
return condHA;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigStimPlanUpscalingCalc::computeAHupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers)
|
||||
{
|
||||
std::vector<double> DrowAvg;
|
||||
std::vector<double> liRowSum;
|
||||
std::vector<double> CondAritRow;
|
||||
|
||||
for (size_t j = 0; j < fracTemplateStimPlan->stimPlanGridNumberOfRows(); j++)
|
||||
{
|
||||
std::vector<double> conductivitiesInStimPlanCells;
|
||||
std::vector<double> lengthsLiOfStimPlanCol;
|
||||
std::vector<double> heightsDiOfStimPlanCells;
|
||||
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellsCol = getRowOfStimPlanCells(stimPlanCells, j);
|
||||
for (RigStimPlanFracTemplateCell* stimPlanCell : stimPlanCellsCol)
|
||||
{
|
||||
if (stimPlanCell->getConductivtyValue() > 1e-7)
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
|
||||
if (clippedStimPlanPolygons.size() > 0)
|
||||
{
|
||||
for (auto clippedStimPlanPolygon : clippedStimPlanPolygons)
|
||||
{
|
||||
conductivitiesInStimPlanCells.push_back(stimPlanCell->getConductivtyValue());
|
||||
lengthsLiOfStimPlanCol.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAlongLayers, clippedStimPlanPolygon));
|
||||
heightsDiOfStimPlanCells.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAcrossLayers, clippedStimPlanPolygon));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
//Calculate sums needed for (arithmetic) average for coloum
|
||||
double sumCondLiDivDi = 0.0;
|
||||
double sumDi = 0.0;
|
||||
double sumLiDi = 0.0;
|
||||
double sumLi = 0.0;
|
||||
for (int i = 0; i < conductivitiesInStimPlanCells.size(); i++)
|
||||
{
|
||||
sumCondLiDivDi += (conductivitiesInStimPlanCells[i] * lengthsLiOfStimPlanCol[i]) / heightsDiOfStimPlanCells[i];
|
||||
sumDi += heightsDiOfStimPlanCells[i];
|
||||
sumLiDi += heightsDiOfStimPlanCells[i] * lengthsLiOfStimPlanCol[i];
|
||||
sumLi += lengthsLiOfStimPlanCol[i];
|
||||
}
|
||||
|
||||
if (sumCondLiDivDi != 0)
|
||||
{
|
||||
//Calculating art avg
|
||||
double dAvg = sumLiDi / sumLi;
|
||||
DrowAvg.push_back(dAvg);
|
||||
liRowSum.push_back(sumLi);
|
||||
CondAritRow.push_back(dAvg / sumLi * sumCondLiDivDi);
|
||||
}
|
||||
}
|
||||
|
||||
//Do harmonic upscaling based on arithmetric upscaled values for coloums
|
||||
double sumDiDivCondALi = 0.0;
|
||||
double sumDi = 0.0;
|
||||
double sumDiLi = 0.0;
|
||||
for (int i = 0; i < CondAritRow.size(); i++)
|
||||
{
|
||||
sumDi += DrowAvg[i];
|
||||
sumDiLi += DrowAvg[i] * liRowSum[i];
|
||||
sumDiDivCondALi += DrowAvg[i] / (CondAritRow[i] * liRowSum[i]);
|
||||
|
||||
}
|
||||
double Lavg = sumDiLi / sumDi;
|
||||
double condAH = (sumDi / Lavg) * (1 / sumDiDivCondALi);
|
||||
|
||||
return condAH;
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigStimPlanUpscalingCalc::arithmeticAverage(std::vector<double> values)
|
||||
{
|
||||
if (values.size() == 0) return cvf::UNDEFINED_DOUBLE;
|
||||
|
||||
double sumValue = 0.0;
|
||||
size_t numberOfValues = 0;
|
||||
for (double value : values)
|
||||
{
|
||||
sumValue += value;
|
||||
numberOfValues++;
|
||||
}
|
||||
|
||||
return sumValue / numberOfValues;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RigStimPlanUpscalingCalc::RigStimPlanUpscalingCalc(RimEclipseCase* caseToApply, RimFracture* fracture)
|
||||
{
|
||||
m_case = caseToApply;
|
||||
m_fracture = fracture;
|
||||
|
||||
|
||||
//Set correct unit system:
|
||||
RigEclipseCaseData::UnitsType caseUnit = m_case->eclipseCaseData()->unitsType();
|
||||
|
||||
if (caseUnit == RigEclipseCaseData::UNITS_METRIC)
|
||||
{
|
||||
RiaLogging::debug(QString("Calculating transmissibilities in metric units"));
|
||||
m_unitForCalculation = RimDefines::UNITS_METRIC;
|
||||
}
|
||||
else if (caseUnit == RigEclipseCaseData::UNITS_FIELD)
|
||||
{
|
||||
RiaLogging::debug(QString("Calculating transmissibilities in field units"));
|
||||
m_unitForCalculation = RimDefines::UNITS_FIELD;
|
||||
}
|
||||
else
|
||||
{
|
||||
//TODO: How to handle lab units for eclipse case?
|
||||
RiaLogging::error(QString("Unit system for case not supported for fracture export."));
|
||||
RiaLogging::error(QString("Export will be in metric units, but results might be wrong."));
|
||||
m_unitForCalculation = RimDefines::UNITS_METRIC;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<RigFracturedEclipseCellExportData> RigStimPlanUpscalingCalc::computeUpscaledPropertyFromStimPlan( QString resultName, QString resultUnit, size_t timeStepIndex)
|
||||
{
|
||||
std::vector<RigFracturedEclipseCellExportData> fracDataVec;
|
||||
|
||||
RimStimPlanFractureTemplate* fracTemplateStimPlan;
|
||||
if (dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition()))
|
||||
{
|
||||
fracTemplateStimPlan = dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition());
|
||||
}
|
||||
else
|
||||
{
|
||||
return fracDataVec;
|
||||
}
|
||||
|
||||
std::vector<std::vector<cvf::Vec3d> > stimPlanCellsAsPolygons;
|
||||
std::vector<double> stimPlanParameterValues;
|
||||
fracTemplateStimPlan->getStimPlanDataAsPolygonsAndValues(stimPlanCellsAsPolygons, stimPlanParameterValues, resultName, resultUnit, timeStepIndex);
|
||||
|
||||
//TODO: A lot of common code with function above... Can be cleaned up...?
|
||||
|
||||
std::vector<size_t> fracCells = m_fracture->getPotentiallyFracturedCells(m_case->eclipseCaseData()->mainGrid());
|
||||
|
||||
RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS;
|
||||
RimReservoirCellResultsStorage* gridCellResults = m_case->results(porosityModel);
|
||||
RigActiveCellInfo* activeCellInfo = m_case->eclipseCaseData()->activeCellInfo(porosityModel);
|
||||
|
||||
for (size_t fracCell : fracCells)
|
||||
{
|
||||
RigFracturedEclipseCellExportData fracData;
|
||||
fracData.reservoirCellIndex = fracCell;
|
||||
|
||||
std::pair<double, double> upscaledCondFlowAcrossLayers = flowAcrossLayersUpscaling(resultName, resultUnit, timeStepIndex, m_unitForCalculation, fracCell);
|
||||
|
||||
double upscaledStimPlanValueHA = upscaledCondFlowAcrossLayers.first;
|
||||
double upscaledStimPlanValueAH = upscaledCondFlowAcrossLayers.second;
|
||||
|
||||
if (upscaledStimPlanValueHA != cvf::UNDEFINED_DOUBLE)
|
||||
{
|
||||
fracData.upscaledStimPlanValueHA = upscaledStimPlanValueHA;
|
||||
fracData.upscaledStimPlanValueAH = upscaledStimPlanValueAH;
|
||||
|
||||
fracDataVec.push_back(fracData);
|
||||
}
|
||||
}
|
||||
|
||||
return fracDataVec;
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
@@ -682,41 +376,6 @@ double RigFractureTransCalc::cDarcy()
|
||||
return c;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<RigStimPlanFracTemplateCell*> RigStimPlanUpscalingCalc::getRowOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t i)
|
||||
{
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellRow;
|
||||
|
||||
for (RigStimPlanFracTemplateCell stimPlanCell : allStimPlanCells)
|
||||
{
|
||||
if (stimPlanCell.getI() == i)
|
||||
{
|
||||
stimPlanCellRow.push_back(&stimPlanCell);
|
||||
}
|
||||
}
|
||||
|
||||
return stimPlanCellRow;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<RigStimPlanFracTemplateCell*> RigStimPlanUpscalingCalc::getColOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t j)
|
||||
{
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellCol;
|
||||
|
||||
for (RigStimPlanFracTemplateCell stimPlanCell : allStimPlanCells)
|
||||
{
|
||||
if (stimPlanCell.getJ() == j)
|
||||
{
|
||||
stimPlanCellCol.push_back(&stimPlanCell);
|
||||
}
|
||||
}
|
||||
|
||||
return stimPlanCellCol;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
|
||||
@@ -89,26 +89,3 @@ private:
|
||||
RimDefines::UnitSystem m_unitForCalculation;
|
||||
|
||||
};
|
||||
|
||||
class RigStimPlanUpscalingCalc
|
||||
{
|
||||
public:
|
||||
explicit RigStimPlanUpscalingCalc(RimEclipseCase* caseToApply, RimFracture* fracture);
|
||||
|
||||
std::vector<RigFracturedEclipseCellExportData> computeUpscaledPropertyFromStimPlan(QString resultName, QString resultUnit, size_t timeStepIndex);
|
||||
|
||||
private:
|
||||
std::pair<double, double> flowAcrossLayersUpscaling(QString resultName, QString resultUnit, size_t timeStepIndex, RimDefines::UnitSystem unitSystem, size_t eclipseCellIndex);
|
||||
double computeHAupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers);
|
||||
double computeAHupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers);
|
||||
static double arithmeticAverage(std::vector<double> values);
|
||||
|
||||
static std::vector<RigStimPlanFracTemplateCell*> getRowOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t i);
|
||||
static std::vector<RigStimPlanFracTemplateCell*> getColOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t j);
|
||||
|
||||
private:
|
||||
RimEclipseCase* m_case;
|
||||
RimFracture* m_fracture;
|
||||
RimDefines::UnitSystem m_unitForCalculation;
|
||||
|
||||
};
|
||||
353
ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.cpp
Normal file
353
ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.cpp
Normal file
@@ -0,0 +1,353 @@
|
||||
#include "RigStimPlanUpscalingCalc.h"
|
||||
#include "RigMainGrid.h"
|
||||
#include "RimEclipseCase.h"
|
||||
#include "RigFractureTransCalc.h"
|
||||
#include "RigEclipseCaseData.h"
|
||||
#include "RigCellGeometryTools.h"
|
||||
#include "RigStimPlanFracTemplateCell.h"
|
||||
#include "RiaLogging.h"
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RigStimPlanUpscalingCalc::RigStimPlanUpscalingCalc(RimEclipseCase* caseToApply, RimFracture* fracture)
|
||||
{
|
||||
m_case = caseToApply;
|
||||
m_fracture = fracture;
|
||||
|
||||
|
||||
//Set correct unit system:
|
||||
RigEclipseCaseData::UnitsType caseUnit = m_case->eclipseCaseData()->unitsType();
|
||||
|
||||
if (caseUnit == RigEclipseCaseData::UNITS_METRIC)
|
||||
{
|
||||
RiaLogging::debug(QString("Calculating transmissibilities in metric units"));
|
||||
m_unitForCalculation = RimDefines::UNITS_METRIC;
|
||||
}
|
||||
else if (caseUnit == RigEclipseCaseData::UNITS_FIELD)
|
||||
{
|
||||
RiaLogging::debug(QString("Calculating transmissibilities in field units"));
|
||||
m_unitForCalculation = RimDefines::UNITS_FIELD;
|
||||
}
|
||||
else
|
||||
{
|
||||
//TODO: How to handle lab units for eclipse case?
|
||||
RiaLogging::error(QString("Unit system for case not supported for fracture export."));
|
||||
RiaLogging::error(QString("Export will be in metric units, but results might be wrong."));
|
||||
m_unitForCalculation = RimDefines::UNITS_METRIC;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::pair<double, double> RigStimPlanUpscalingCalc::flowAcrossLayersUpscaling(QString resultName, QString resultUnit, size_t timeStepIndex, RimDefines::UnitSystem unitSystem, size_t eclipseCellIndex)
|
||||
{
|
||||
RimStimPlanFractureTemplate* fracTemplateStimPlan;
|
||||
if (dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition()))
|
||||
{
|
||||
fracTemplateStimPlan = dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition());
|
||||
}
|
||||
else return std::make_pair(cvf::UNDEFINED_DOUBLE, cvf::UNDEFINED_DOUBLE);
|
||||
|
||||
std::vector<RigStimPlanFracTemplateCell> stimPlanCells = fracTemplateStimPlan->getStimPlanCells();
|
||||
|
||||
|
||||
|
||||
cvf::Vec3d localX, localY, localZ; //Not used in calculation here, but needed for function to find planCellPolygons
|
||||
std::vector<std::vector<cvf::Vec3d> > planeCellPolygons;
|
||||
|
||||
const RigMainGrid* mainGrid = m_case->eclipseCaseData()->mainGrid();
|
||||
cvf::Vec3d hexCorners[8];
|
||||
mainGrid->cellCornerVertices(eclipseCellIndex, hexCorners);
|
||||
|
||||
|
||||
|
||||
bool isPlanIntersected = RigFractureTransCalc::planeCellIntersectionPolygons(hexCorners, m_fracture->transformMatrix(), planeCellPolygons);
|
||||
if (!isPlanIntersected || planeCellPolygons.size() == 0) return std::make_pair(cvf::UNDEFINED_DOUBLE, cvf::UNDEFINED_DOUBLE);
|
||||
|
||||
//Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon/stimPlan mesh already is located)
|
||||
cvf::Mat4f invertedTransMatrix = m_fracture->transformMatrix().getInverted();
|
||||
for (std::vector<cvf::Vec3d> & planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
for (cvf::Vec3d& v : planeCellPolygon)
|
||||
{
|
||||
v.transformPoint(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
}
|
||||
}
|
||||
|
||||
cvf::Vec3d directionAcrossLayers;
|
||||
cvf::Vec3d directionAlongLayers;
|
||||
directionAcrossLayers = cvf::Vec3d(0.0, -1.0, 0.0);
|
||||
directionAlongLayers = cvf::Vec3d(1.0, 0.0, 0.0);
|
||||
|
||||
std::vector<cvf::Vec3f> fracPolygon = m_fracture->attachedFractureDefinition()->fracturePolygon(unitSystem);
|
||||
std::vector<std::vector<cvf::Vec3d> > polygonsDescribingFractureInCell;
|
||||
|
||||
std::vector<double> upscaledConductivitiesHA;
|
||||
std::vector<double> upscaledConductivitiesAH;
|
||||
|
||||
for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
double condHA = computeHAupscale(fracTemplateStimPlan, stimPlanCells, planeCellPolygon, directionAlongLayers, directionAcrossLayers);
|
||||
upscaledConductivitiesHA.push_back(condHA);
|
||||
|
||||
double condAH = computeAHupscale(fracTemplateStimPlan, stimPlanCells, planeCellPolygon, directionAlongLayers, directionAcrossLayers);
|
||||
upscaledConductivitiesAH.push_back(condAH);
|
||||
}
|
||||
|
||||
return std::make_pair(arithmeticAverage(upscaledConductivitiesHA), arithmeticAverage(upscaledConductivitiesAH));
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigStimPlanUpscalingCalc::computeHAupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers)
|
||||
{
|
||||
std::vector<double> DcolSum;
|
||||
std::vector<double> lavgCol;
|
||||
std::vector<double> CondHarmCol;
|
||||
|
||||
for (size_t j = 0; j < fracTemplateStimPlan->stimPlanGridNumberOfColums(); j++)
|
||||
{
|
||||
std::vector<double> conductivitiesInStimPlanCells;
|
||||
std::vector<double> lengthsLiOfStimPlanCol;
|
||||
std::vector<double> heightsDiOfStimPlanCells;
|
||||
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellsCol = getColOfStimPlanCells(stimPlanCells, j);
|
||||
for (RigStimPlanFracTemplateCell* stimPlanCell : stimPlanCellsCol)
|
||||
{
|
||||
if (stimPlanCell->getConductivtyValue() > 1e-7)
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
|
||||
if (clippedStimPlanPolygons.size() > 0)
|
||||
{
|
||||
for (auto clippedStimPlanPolygon : clippedStimPlanPolygons)
|
||||
{
|
||||
conductivitiesInStimPlanCells.push_back(stimPlanCell->getConductivtyValue());
|
||||
lengthsLiOfStimPlanCol.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAlongLayers, clippedStimPlanPolygon));
|
||||
heightsDiOfStimPlanCells.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAcrossLayers, clippedStimPlanPolygon));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
//Regne ut average
|
||||
double sumDiDivCondLi = 0.0;
|
||||
double sumDi = 0.0;
|
||||
double sumLiDi = 0.0;
|
||||
for (int i = 0; i < conductivitiesInStimPlanCells.size(); i++)
|
||||
{
|
||||
sumDiDivCondLi += heightsDiOfStimPlanCells[i] / (conductivitiesInStimPlanCells[i] * lengthsLiOfStimPlanCol[i]);
|
||||
sumDi += heightsDiOfStimPlanCells[i];
|
||||
sumLiDi += heightsDiOfStimPlanCells[i] * lengthsLiOfStimPlanCol[i];
|
||||
}
|
||||
|
||||
if (sumDiDivCondLi != 0)
|
||||
{
|
||||
DcolSum.push_back(sumDi);
|
||||
double lAvgValue = sumLiDi / sumDi;
|
||||
lavgCol.push_back(lAvgValue);
|
||||
double harmMeanForCol = (sumDi / lAvgValue) * (1 / sumDiDivCondLi);
|
||||
CondHarmCol.push_back(harmMeanForCol);
|
||||
}
|
||||
}
|
||||
|
||||
//Do arithmetic upscaling based on harmonic upscaled values for coloums
|
||||
double sumCondHLiDivDi = 0.0;
|
||||
double sumLi = 0.0;
|
||||
double sumDiLi = 0.0;
|
||||
for (int i = 0; i < CondHarmCol.size(); i++)
|
||||
{
|
||||
sumLi += lavgCol[i];
|
||||
sumDiLi += DcolSum[i] * lavgCol[i];
|
||||
sumCondHLiDivDi += CondHarmCol[i] * lavgCol[i] / DcolSum[i];
|
||||
}
|
||||
double Davg = sumDiLi / sumLi;
|
||||
double condHA = (Davg / sumLi) * sumCondHLiDivDi;
|
||||
|
||||
return condHA;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigStimPlanUpscalingCalc::computeAHupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers)
|
||||
{
|
||||
std::vector<double> DrowAvg;
|
||||
std::vector<double> liRowSum;
|
||||
std::vector<double> CondAritRow;
|
||||
|
||||
for (size_t j = 0; j < fracTemplateStimPlan->stimPlanGridNumberOfRows(); j++)
|
||||
{
|
||||
std::vector<double> conductivitiesInStimPlanCells;
|
||||
std::vector<double> lengthsLiOfStimPlanCol;
|
||||
std::vector<double> heightsDiOfStimPlanCells;
|
||||
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellsCol = getRowOfStimPlanCells(stimPlanCells, j);
|
||||
for (RigStimPlanFracTemplateCell* stimPlanCell : stimPlanCellsCol)
|
||||
{
|
||||
if (stimPlanCell->getConductivtyValue() > 1e-7)
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d> >clippedStimPlanPolygons = RigCellGeometryTools::clipPolygons(stimPlanCell->getPolygon(), planeCellPolygon);
|
||||
if (clippedStimPlanPolygons.size() > 0)
|
||||
{
|
||||
for (auto clippedStimPlanPolygon : clippedStimPlanPolygons)
|
||||
{
|
||||
conductivitiesInStimPlanCells.push_back(stimPlanCell->getConductivtyValue());
|
||||
lengthsLiOfStimPlanCol.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAlongLayers, clippedStimPlanPolygon));
|
||||
heightsDiOfStimPlanCells.push_back(RigCellGeometryTools::polygonAreaWeightedLength(directionAcrossLayers, clippedStimPlanPolygon));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
//Calculate sums needed for (arithmetic) average for coloum
|
||||
double sumCondLiDivDi = 0.0;
|
||||
double sumDi = 0.0;
|
||||
double sumLiDi = 0.0;
|
||||
double sumLi = 0.0;
|
||||
for (int i = 0; i < conductivitiesInStimPlanCells.size(); i++)
|
||||
{
|
||||
sumCondLiDivDi += (conductivitiesInStimPlanCells[i] * lengthsLiOfStimPlanCol[i]) / heightsDiOfStimPlanCells[i];
|
||||
sumDi += heightsDiOfStimPlanCells[i];
|
||||
sumLiDi += heightsDiOfStimPlanCells[i] * lengthsLiOfStimPlanCol[i];
|
||||
sumLi += lengthsLiOfStimPlanCol[i];
|
||||
}
|
||||
|
||||
if (sumCondLiDivDi != 0)
|
||||
{
|
||||
//Calculating art avg
|
||||
double dAvg = sumLiDi / sumLi;
|
||||
DrowAvg.push_back(dAvg);
|
||||
liRowSum.push_back(sumLi);
|
||||
CondAritRow.push_back(dAvg / sumLi * sumCondLiDivDi);
|
||||
}
|
||||
}
|
||||
|
||||
//Do harmonic upscaling based on arithmetric upscaled values for coloums
|
||||
double sumDiDivCondALi = 0.0;
|
||||
double sumDi = 0.0;
|
||||
double sumDiLi = 0.0;
|
||||
for (int i = 0; i < CondAritRow.size(); i++)
|
||||
{
|
||||
sumDi += DrowAvg[i];
|
||||
sumDiLi += DrowAvg[i] * liRowSum[i];
|
||||
sumDiDivCondALi += DrowAvg[i] / (CondAritRow[i] * liRowSum[i]);
|
||||
|
||||
}
|
||||
double Lavg = sumDiLi / sumDi;
|
||||
double condAH = (sumDi / Lavg) * (1 / sumDiDivCondALi);
|
||||
|
||||
return condAH;
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigStimPlanUpscalingCalc::arithmeticAverage(std::vector<double> values)
|
||||
{
|
||||
if (values.size() == 0) return cvf::UNDEFINED_DOUBLE;
|
||||
|
||||
double sumValue = 0.0;
|
||||
size_t numberOfValues = 0;
|
||||
for (double value : values)
|
||||
{
|
||||
sumValue += value;
|
||||
numberOfValues++;
|
||||
}
|
||||
|
||||
return sumValue / numberOfValues;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<RigFracturedEclipseCellExportData> RigStimPlanUpscalingCalc::computeUpscaledPropertyFromStimPlan( QString resultName, QString resultUnit, size_t timeStepIndex)
|
||||
{
|
||||
std::vector<RigFracturedEclipseCellExportData> fracDataVec;
|
||||
|
||||
RimStimPlanFractureTemplate* fracTemplateStimPlan;
|
||||
if (dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition()))
|
||||
{
|
||||
fracTemplateStimPlan = dynamic_cast<RimStimPlanFractureTemplate*>(m_fracture->attachedFractureDefinition());
|
||||
}
|
||||
else
|
||||
{
|
||||
return fracDataVec;
|
||||
}
|
||||
|
||||
std::vector<std::vector<cvf::Vec3d> > stimPlanCellsAsPolygons;
|
||||
std::vector<double> stimPlanParameterValues;
|
||||
fracTemplateStimPlan->getStimPlanDataAsPolygonsAndValues(stimPlanCellsAsPolygons, stimPlanParameterValues, resultName, resultUnit, timeStepIndex);
|
||||
|
||||
//TODO: A lot of common code with function above... Can be cleaned up...?
|
||||
|
||||
std::vector<size_t> fracCells = m_fracture->getPotentiallyFracturedCells(m_case->eclipseCaseData()->mainGrid());
|
||||
|
||||
RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS;
|
||||
RimReservoirCellResultsStorage* gridCellResults = m_case->results(porosityModel);
|
||||
RigActiveCellInfo* activeCellInfo = m_case->eclipseCaseData()->activeCellInfo(porosityModel);
|
||||
|
||||
for (size_t fracCell : fracCells)
|
||||
{
|
||||
RigFracturedEclipseCellExportData fracData;
|
||||
fracData.reservoirCellIndex = fracCell;
|
||||
|
||||
std::pair<double, double> upscaledCondFlowAcrossLayers = flowAcrossLayersUpscaling(resultName, resultUnit, timeStepIndex, m_unitForCalculation, fracCell);
|
||||
|
||||
double upscaledStimPlanValueHA = upscaledCondFlowAcrossLayers.first;
|
||||
double upscaledStimPlanValueAH = upscaledCondFlowAcrossLayers.second;
|
||||
|
||||
if (upscaledStimPlanValueHA != cvf::UNDEFINED_DOUBLE)
|
||||
{
|
||||
fracData.upscaledStimPlanValueHA = upscaledStimPlanValueHA;
|
||||
fracData.upscaledStimPlanValueAH = upscaledStimPlanValueAH;
|
||||
|
||||
fracDataVec.push_back(fracData);
|
||||
}
|
||||
}
|
||||
|
||||
return fracDataVec;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<RigStimPlanFracTemplateCell*> RigStimPlanUpscalingCalc::getRowOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t i)
|
||||
{
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellRow;
|
||||
|
||||
for (RigStimPlanFracTemplateCell stimPlanCell : allStimPlanCells)
|
||||
{
|
||||
if (stimPlanCell.getI() == i)
|
||||
{
|
||||
stimPlanCellRow.push_back(&stimPlanCell);
|
||||
}
|
||||
}
|
||||
|
||||
return stimPlanCellRow;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<RigStimPlanFracTemplateCell*> RigStimPlanUpscalingCalc::getColOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t j)
|
||||
{
|
||||
std::vector<RigStimPlanFracTemplateCell*> stimPlanCellCol;
|
||||
|
||||
for (RigStimPlanFracTemplateCell stimPlanCell : allStimPlanCells)
|
||||
{
|
||||
if (stimPlanCell.getJ() == j)
|
||||
{
|
||||
stimPlanCellCol.push_back(&stimPlanCell);
|
||||
}
|
||||
}
|
||||
|
||||
return stimPlanCellCol;
|
||||
}
|
||||
|
||||
@@ -0,0 +1,30 @@
|
||||
#pragma once
|
||||
|
||||
#include "RimEclipseCase.h"
|
||||
#include "RimFracture.h"
|
||||
#include "RimStimPlanFractureTemplate.h"
|
||||
|
||||
|
||||
class RigStimPlanUpscalingCalc
|
||||
{
|
||||
public:
|
||||
explicit RigStimPlanUpscalingCalc(RimEclipseCase* caseToApply, RimFracture* fracture);
|
||||
|
||||
std::vector<RigFracturedEclipseCellExportData> computeUpscaledPropertyFromStimPlan(QString resultName, QString resultUnit, size_t timeStepIndex);
|
||||
|
||||
private:
|
||||
std::pair<double, double> flowAcrossLayersUpscaling(QString resultName, QString resultUnit, size_t timeStepIndex, RimDefines::UnitSystem unitSystem, size_t eclipseCellIndex);
|
||||
double computeHAupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers);
|
||||
double computeAHupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector<RigStimPlanFracTemplateCell> stimPlanCells, std::vector<cvf::Vec3d> planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers);
|
||||
static double arithmeticAverage(std::vector<double> values);
|
||||
|
||||
static std::vector<RigStimPlanFracTemplateCell*> getRowOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t i);
|
||||
static std::vector<RigStimPlanFracTemplateCell*> getColOfStimPlanCells(std::vector<RigStimPlanFracTemplateCell>& allStimPlanCells, size_t j);
|
||||
|
||||
private:
|
||||
RimEclipseCase* m_case;
|
||||
RimFracture* m_fracture;
|
||||
RimDefines::UnitSystem m_unitForCalculation;
|
||||
|
||||
};
|
||||
|
||||
Reference in New Issue
Block a user