#1091, #1092 - pre-proto - Simplifying calulation of projected areas and updating export functionality to include intermediate results of transmissibility calculation

This commit is contained in:
astridkbjorke 2017-01-27 10:46:45 +01:00
parent fccd7d9314
commit 035107b48b
4 changed files with 142 additions and 75 deletions

View File

@ -34,6 +34,7 @@
#include <QTextStream>
#include <QFile>
#include <QString>
@ -59,6 +60,15 @@ RifEclipseExportTools::~RifEclipseExportTools()
//--------------------------------------------------------------------------------------------------
bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, const std::vector< RimFracture*>& fractures)
{
RiaApplication* app = RiaApplication::instance();
RimView* activeView = RiaApplication::instance()->activeReservoirView();
if (!activeView) return false;
RimEclipseView* activeRiv = dynamic_cast<RimEclipseView*>(activeView);
if (!activeRiv) return false;
const RigMainGrid* mainGrid = activeRiv->mainGrid();
if (!mainGrid) return false;
QFile file(fileName);
if (!file.open(QIODevice::WriteOnly | QIODevice::Text))
@ -66,11 +76,6 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c
return false;
}
QTextStream out(&file);
out << "\n";
out << "-- Exported from ResInsight" << "\n\n";
out << "COMPDAT" << "\n" << right << qSetFieldWidth(8);
caf::ProgressInfo pi(fractures.size(), QString("Writing data to file %1").arg(fileName));
RimEclipseWell* simWell = nullptr;
@ -79,6 +84,89 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c
size_t progress =0;
std::vector<size_t> ijk;
QTextStream out(&file);
out << "\n";
out << "-- Exported from ResInsight" << "\n\n";
out << "-- Background data for calculation" << "\n\n";
//Write header line
out << qSetFieldWidth(4);
out << "--";
out << qSetFieldWidth(12);
out << "Well name";
out << qSetFieldWidth(16);
out << "Fracture Name";
out << qSetFieldWidth(5);
out << "i";
out << "j";
out << "k";
out << qSetFieldWidth(8);
out << "Ax";
out << " ";
out << "Ay";
out << " ";
out << "Az";
out << " ";
out << "TotalArea";
out << "\n";
//Write background numbers (with -- first to comment out line)
for (RimFracture* fracture : fractures)
{
fracture->computeTransmissibility();
std::vector<RigFractureData> fracDataVector = fracture->attachedRigFracture()->fractureData();
for (RigFractureData fracData : fracDataVector)
{
out << qSetFieldWidth(4);
out << "--";
out << qSetFieldWidth(12);
wellPath, simWell = nullptr;
fracture->firstAncestorOrThisOfType(simWell);
if (simWell) out << simWell->name; // 1. Well name
fracture->firstAncestorOrThisOfType(wellPath);
if (wellPath) out << wellPath->name; // 1. Well name
out << qSetFieldWidth(16);
out << fracture->name();
out << qSetFieldWidth(5);
size_t i, j, k;
mainGrid->ijkFromCellIndex(fracData.reservoirCellIndex, &i, &j, &k);
out << i + 1; // 2. I location grid block, adding 1 to go to eclipse 1-based grid definition
out << j + 1; // 3. J location grid block, adding 1 to go to eclipse 1-based grid definition
out << k + 1; // 4. K location of upper connecting grid block, adding 1 to go to eclipse 1-based grid definition
out << " ";
out << fracData.projectedAreas.x();
out << " ";
out << fracData.projectedAreas.y();
out << " ";
out << fracData.projectedAreas.z();
out << " ";
out << fracData.totalArea;
out << "\n";
}
}
out << "\n";
//Write COMPDAT keyword and entries
out << "COMPDAT" << "\n" << right << qSetFieldWidth(8);
for (RimFracture* fracture : fractures)
{
fracture->computeTransmissibility();
@ -97,15 +185,6 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c
out << qSetFieldWidth(5);
RiaApplication* app = RiaApplication::instance();
RimView* activeView = RiaApplication::instance()->activeReservoirView();
if (!activeView) return false;
RimEclipseView* activeRiv = dynamic_cast<RimEclipseView*>(activeView);
if (!activeRiv) return false;
const RigMainGrid* mainGrid = activeRiv->mainGrid();
if (!mainGrid) return false;
size_t i, j, k;
mainGrid->ijkFromCellIndex(fracData.reservoirCellIndex, &i, &j, &k);
out << i+1; // 2. I location grid block, adding 1 to go to eclipse 1-based grid definition
@ -113,8 +192,9 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c
out << k+1; // 4. K location of upper connecting grid block, adding 1 to go to eclipse 1-based grid definition
out << k+1; // 5. K location of lower connecting grid block, adding 1 to go to eclipse 1-based grid definition
out << "OPEN"; // 6. Open / Shut flag of connection
out << "1* "; // 7. Saturation table number for connection rel perm. Default value
out << "2*"; // Default value for
//6. Open / Shut flag of connection
// 7. Saturation table number for connection rel perm. Default value
out << qSetFieldWidth(8);
// 8. Transmissibility
@ -122,29 +202,27 @@ bool RifEclipseExportTools::writeFracturesToTextFile(const QString& fileName, c
else out << "UNDEF";
out << qSetFieldWidth(4);
out << "1* "; // 9. Well bore diameter. Set to default
out << qSetFieldWidth(8);
out << "2*"; // Default value for
// 9. Well bore diameter. Set to default
// 10. Effective Kh (perm times width)
if (fracture->attachedFractureDefinition())
{
out << fracture->attachedFractureDefinition()->effectiveKh(); // 10. Effective Kh (perm times width)
out << qSetFieldWidth(4);
out << fracture->attachedFractureDefinition()->skinFactor; // 11. Skin factor
}
else //If no attached fracture definition these parameters are set to UNDEF
{
out << "UNDEF";
out << qSetFieldWidth(4);
out << "UNDEF";
}
out << "1*"; // 12. D-factor for handling non-Darcy flow of free gas. Default value.
out << "Z"; // 13. Direction well is penetrating the grid block. Z is default.
out << "1*"; // 14. Pressure equivalent radius, Default
out << "/" << "\n";
out << "/";
out << " " << fracture->name(); //Fracture name as comment
out << "\n"; // Terminating entry
}
//TODO: If same cell is used for multiple fractures, the sum of contributions should be added to table.
progress++;
pi.setProgress(progress);
}

View File

@ -253,8 +253,8 @@ void RimFracture::computeTransmissibility()
std::vector<RigFractureData> fracDataVec;
std::vector<size_t> fracCells = getPotentiallyFracturedCells();
for (auto fracCell : fracCells)
for (size_t fracCell : fracCells)
{
//TODO: Remove - only for simplifying debugging...
RiaApplication* app = RiaApplication::instance();
@ -274,6 +274,7 @@ void RimFracture::computeTransmissibility()
//Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon already is located)
cvf::Mat4f invertedTransMatrix = transformMatrix().getInverted();
for (std::vector<cvf::Vec3d> & planeCellPolygon : planeCellPolygons)
{
for (cvf::Vec3d& v : planeCellPolygon)
@ -287,11 +288,18 @@ void RimFracture::computeTransmissibility()
localZinFracPlane.transformVector(static_cast<cvf::Mat4d>(invertedTransMatrix));
cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO;
directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1));
directionOfLength.normalize();
RigFractureData fracData;
fracData.reservoirCellIndex = fracCell;
double transmissibility;
double fractureArea = 0.0;
double fractureAreaWeightedlength = 0.0;
double Ax = 0.0;
double Ay = 0.0;
double Az = 0.0;
if (attachedFractureDefinition())
{
@ -326,34 +334,26 @@ void RimFracture::computeTransmissibility()
areaOfFractureParts.push_back(area);
length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon);
lengthXareaOfFractureParts.push_back(length * area);
double AreaX = calculateProjectedArea(fracturePartPolygon, localX);
//Calculating area in x, y and z direction
cvf::Vec3d planeNormal = cvf::Vec3d::ZERO;
planeNormal.cross(localX, -localZ);
double Ax = calculateProjectedArea(fracturePartPolygon, planeNormal);
planeNormal.cross(localY, localZ);
double Ay = calculateProjectedArea(fracturePartPolygon, planeNormal);
planeNormal.cross(localX, localY);
double Az = calculateProjectedArea(fracturePartPolygon, planeNormal);
cvf::Plane fracturePlane;
cvf::Mat4f m = transformMatrix();
bool isCellIntersected = false;
fracturePlane.setFromPointAndNormal(static_cast<cvf::Vec3d>(m.translation()),
static_cast<cvf::Vec3d>(m.col(2)));
Ax += abs(area*(fracturePlane.normal().dot(localX)));
Ay += abs(area*(fracturePlane.normal().dot(localY)));
Az += abs(area*(fracturePlane.normal().dot(localZ)));
}
double totalArea = 0.0;
for (double area : areaOfFractureParts) totalArea += area;
fractureArea = 0.0;
for (double area : areaOfFractureParts) fractureArea += area;
double totalAreaXLength = 0.0;
for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea;
double fractureAreaWeightedlength = totalAreaXLength / totalArea;
double fractureAreaWeightedlength = totalAreaXLength / fractureArea;
//TODO: FInd direction for length calculation (normal to z, in fracture plane)
double flowLength = 2.718281828;
@ -361,13 +361,17 @@ void RimFracture::computeTransmissibility()
//TODO: read permeability from file (should use matrix permeability and not fracture perm)...
transmissibility = totalArea;
//transmissibility = 8 * c * attachedFractureDefinition()->permeability * areaOfCellPlaneFractureOverlap /
transmissibility = 8 * c * attachedFractureDefinition()->permeability * fractureArea /
( flowLength + (attachedFractureDefinition()->skinFactor * fractureAreaWeightedlength) / cvf::PI_D);
}
else transmissibility = cvf::UNDEFINED_DOUBLE;
else transmissibility = cvf::UNDEFINED_DOUBLE;
fracData.transmissibility = transmissibility;
fracData.fractureLenght = fractureAreaWeightedlength;
fracData.totalArea = fractureArea;
fracData.projectedAreas = cvf::Vec3d(Ax, Ay, Az);
//only keep fracData if transmissibility is non-zero
if (transmissibility > 0)
{
@ -438,24 +442,6 @@ bool RimFracture::planeCellIntersectionPolygons(size_t cellindex, std::vector<st
return isCellIntersected;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RimFracture::calculateProjectedArea(std::vector<cvf::Vec3d> polygon, cvf::Vec3d planeNormal)
{
//Set up plane
cvf::Plane plane;
plane.setFromPointAndNormal(polygon[0], planeNormal);
//Project points
for (cvf::Vec3d v : polygon) plane.projectPoint(v);
//calculate Area
cvf::Vec3d areaVector = cvf::GeometryTools::polygonAreaNormal3D(polygon);
double AreaInPlane = areaVector.length();
return AreaInPlane;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

View File

@ -93,8 +93,6 @@ private:
//Functions for area calculations - should these be in separate class
bool planeCellIntersectionPolygons(size_t cellindex, std::vector<std::vector<cvf::Vec3d> > & polygons, cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ);
double calculateProjectedArea(std::vector<cvf::Vec3d> fracturePartPolygon, cvf::Vec3d localX);
protected:
caf::PdmPtrField<RimEllipseFractureTemplate*> m_fractureTemplate;

View File

@ -33,6 +33,11 @@ public:
size_t reservoirCellIndex;
double transmissibility;
double totalArea;
double fractureLenght;
cvf::Vec3d projectedAreas;
};