diff --git a/ApplicationCode/FileInterface/CMakeLists_files.cmake b/ApplicationCode/FileInterface/CMakeLists_files.cmake index a72f417f28..f486dae921 100644 --- a/ApplicationCode/FileInterface/CMakeLists_files.cmake +++ b/ApplicationCode/FileInterface/CMakeLists_files.cmake @@ -22,7 +22,6 @@ ${CEE_CURRENT_LIST_DIR}RifReaderMockModel.h ${CEE_CURRENT_LIST_DIR}RifReaderSettings.h ${CEE_CURRENT_LIST_DIR}RifEclipseSummaryAddress.h ${CEE_CURRENT_LIST_DIR}RifWellPathImporter.h -${CEE_CURRENT_LIST_DIR}RifFractureExportTools.h ${CEE_CURRENT_LIST_DIR}RifStimPlanXmlReader.h ) @@ -44,7 +43,6 @@ ${CEE_CURRENT_LIST_DIR}RifReaderMockModel.cpp ${CEE_CURRENT_LIST_DIR}RifReaderSettings.cpp ${CEE_CURRENT_LIST_DIR}RifEclipseSummaryAddress.cpp ${CEE_CURRENT_LIST_DIR}RifWellPathImporter.cpp -${CEE_CURRENT_LIST_DIR}RifFractureExportTools.cpp ${CEE_CURRENT_LIST_DIR}RifStimPlanXmlReader.cpp ) diff --git a/ApplicationCode/FileInterface/RifFractureExportTools.cpp b/ApplicationCode/FileInterface/RifFractureExportTools.cpp deleted file mode 100644 index cb30968d20..0000000000 --- a/ApplicationCode/FileInterface/RifFractureExportTools.cpp +++ /dev/null @@ -1,650 +0,0 @@ -///////////////////////////////////////////////////////////////////////////////// -// -// Copyright (C) 2017- Statoil ASA -// -// ResInsight is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. -// -// See the GNU General Public License at -// for more details. -// -///////////////////////////////////////////////////////////////////////////////// - -#include "RifFractureExportTools.h" - -#include "RiaApplication.h" -#include "RiaLogging.h" - -#include "RigEclipseCaseData.h" -#include "RigFractureTransCalc.h" -#include "RigFractureCell.h" -#include "RigMainGrid.h" -#include "RigEclipseToStimPlanCellTransmissibilityCalculator.h" -#include "RigFractureTransmissibilityEquations.h" - -#include "RimEclipseCase.h" -#include "RimEclipseResultDefinition.h" -#include "RimEclipseView.h" -#include "RimEclipseWell.h" -#include "RimEllipseFractureTemplate.h" -#include "RimFracture.h" -#include "RimFractureTemplate.h" -#include "RimSimWellFracture.h" -#include "RimStimPlanFractureTemplate.h" -#include "RimWellPath.h" - -#include "cafProgressInfo.h" - -#include -#include -#include -#include "RigStimPlanUpscalingCalc.h" -#include "RigTransmissibilityCondenser.h" -#include "RigWellPathStimplanIntersector.h" -#include "RigFractureGrid.h" - - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -bool RifFractureExportTools::exportFracturesToEclipseDataInputFile(const QString& fileName, const std::vector< RimFracture*>& fractures, RimEclipseCase* caseToApply) -{ - RiaLogging::info(QString("Computing and writing COMPDAT values to file %1").arg(fileName)); - - const RigMainGrid* mainGrid = caseToApply->eclipseCaseData()->mainGrid(); - if (!mainGrid) return false; - - QFile file(fileName); - if (!file.open(QIODevice::WriteOnly | QIODevice::Text)) - { - return false; - } - - caf::ProgressInfo pi(fractures.size(), QString("Writing data to file %1").arg(fileName)); - - size_t progress = 0; - std::vector ijk; - - QTextStream out(&file); - out << "\n"; - out << "-- Exported from ResInsight" << "\n"; - - QString wellName; - { - RimEclipseWell* simWell = nullptr; - fractures[0]->firstAncestorOrThisOfType(simWell); - if ( simWell ) wellName = simWell->name; - - RimWellPath* wellPath = nullptr; - fractures[0]->firstAncestorOrThisOfType(wellPath); - if ( wellPath ) wellName = wellPath->name; - } - - RigEclipseCaseData::UnitsType caseUnit = caseToApply->eclipseCaseData()->unitsType(); - if (caseUnit == RigEclipseCaseData::UNITS_METRIC) out << "-- Using metric unit system" << "\n"; - if (caseUnit == RigEclipseCaseData::UNITS_FIELD) out << "-- Using field unit system" << "\n"; - out << "\n"; - - //Included for debug / prototyping only - printTransmissibilityFractureToWell(fractures, out, caseToApply); - - printStimPlanFractureTrans(fractures, caseToApply, out); - - printStimPlanCellsMatrixTransContributions(fractures, caseToApply, out, wellName, mainGrid); - - printBackgroundDataHeaderLine(out); - - RiaLogging::debug(QString("Writing intermediate results from COMPDAT calculation")); - - std::map > exportDataPrFracture; - - for (RimFracture* fracture : fractures) - { - RigFractureTransCalc transmissibilityCalculator(caseToApply, fracture); - - //TODO: Check that there is a fracture template available for given fracture.... - - std::vector fracDataVector = transmissibilityCalculator.computeTransmissibilityFromPolygonWithInfiniteConductivityInFracture(); - exportDataPrFracture[fracture] = fracDataVector; - - for (RigFracturedEclipseCellExportData fracData : fracDataVector) - { - printBackgroundData(out, wellName, fracture, mainGrid, fracData); - } - } - - out << "\n"; - - out << qSetFieldWidth(7) << "COMPDAT" << "\n" << right << qSetFieldWidth(8); - for (RimFracture* fracture : fractures) - { - RiaLogging::debug(QString("Writing COMPDAT values for fracture %1").arg(fracture->name())); - std::vector fracDataVector = exportDataPrFracture[fracture]; - - double skinFactor = cvf::UNDEFINED_DOUBLE; - if (fracture->fractureTemplate()) skinFactor = fracture->fractureTemplate()->skinFactor(); - QString fractureName = fracture->name(); - - for (RigFracturedEclipseCellExportData fracData : fracDataVector) - { - if ( fracData.transmissibility > 0 ) - { - size_t i, j, k; - mainGrid->ijkFromCellIndex(fracData.reservoirCellIndex, &i, &j, &k); - printCOMPDATvalues(out, fracData.transmissibility, i, j, k, fractureName, skinFactor, wellName); - } - } - - //TODO: If same cell is used for multiple fractures, the sum of contributions should be added to table. - - progress++; - pi.setProgress(progress); - } - - out << "/ \n"; - - RiaLogging::info(QString("Competed writing COMPDAT data to file %1").arg(fileName)); - return true; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RifFractureExportTools::printCOMPDATvalues(QTextStream & out, - double transmissibility, - size_t i, size_t j, size_t k, - const QString& fractureName, - double skinFactor, - const QString& wellName) -{ - out << qSetFieldWidth(8); - - if (transmissibility == cvf::UNDEFINED_DOUBLE || skinFactor == cvf::UNDEFINED_DOUBLE) - { - out << "--"; //Commenting out line in output file - } - - out << wellName; - out << qSetFieldWidth(5); - - - 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 << k + 1; // 5. K location of lower connecting grid block, adding 1 to go to eclipse 1-based grid definition - - out << "2* "; // Default value for - //6. Open / Shut flag of connection - // 7. Saturation table number for connection rel perm. Default value - - out << qSetFieldWidth(12); - - // 8. Transmissibility - if (transmissibility != cvf::UNDEFINED_DOUBLE) - { - out << QString::number(transmissibility, 'e', 4); - } - else - { - out << "UNDEF"; - } - - out << qSetFieldWidth(4); - out << "2* "; // Default value for - // 9. Well bore diameter. Set to default - // 10. Effective Kh (perm times width) - - - if (skinFactor != cvf::UNDEFINED_DOUBLE) - { - out << skinFactor; // 11. Skin factor - } - else //If no attached fracture definition these parameters are set to UNDEF - { - out << "UNDEF"; - } - - out << "/"; - out << " " << fractureName; //Fracture name as comment - out << "\n"; // Terminating entry - -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RifFractureExportTools::printStimPlanCellsMatrixTransContributions(const std::vector& fractures, - RimEclipseCase* caseToApply, - QTextStream &out, - const QString& wellName, - const RigMainGrid* mainGrid) -{ - out << "StimPlan cells' matrix transmissibility and Eclipse Cell contributions \n"; - - out << qSetFieldWidth(4); - out << "-- "; - - out << qSetFieldWidth(12); - out << "Well name "; // 1. Well name - - out << qSetFieldWidth(16); - out << "Fracture name "; - - out << qSetFieldWidth(5); - out << "Ec i"; - out << "Ec j"; - out << "Ec k"; - - out << qSetFieldWidth(10); - out << "Ecl cell"; - - out << qSetFieldWidth(5); - out << "SP i"; - out << "SP j"; - - out << qSetFieldWidth(10); - out << "Tm contr"; - - out << "\n"; - - - - for (RimFracture* fracture : fractures) - { - RimStimPlanFractureTemplate* fracTemplateStimPlan; - if (dynamic_cast(fracture->fractureTemplate())) - { - fracTemplateStimPlan = dynamic_cast(fracture->fractureTemplate()); - } - else continue; - - double cDarcyInCorrectUnit = caseToApply->eclipseCaseData()->darchysValue(); - - std::vector stimPlanCells = fracTemplateStimPlan->fractureGrid()->fractureCells(); - - for (RigFractureCell stimPlanCell : stimPlanCells) - { - if (stimPlanCell.getConductivtyValue() < 1e-7) - { - continue; - } - - RigEclipseToStimPlanCellTransmissibilityCalculator eclToStimPlanTransCalc(caseToApply, - fracture->transformMatrix(), - fracture->fractureTemplate()->skinFactor, - cDarcyInCorrectUnit, - stimPlanCell); - - std::vector stimPlanContributingEclipseCells = eclToStimPlanTransCalc.globalIndeciesToContributingEclipseCells(); - std::vector stimPlanContributingEclipseCellTransmissibilities = eclToStimPlanTransCalc.contributingEclipseCellTransmissibilities(); - - for (size_t i = 0; i < stimPlanContributingEclipseCells.size(); i++) - { - out << qSetFieldWidth(4); - out << "-- "; - - out << qSetFieldWidth(12); - out << wellName + " "; - out << qSetFieldWidth(16); - out << fracture->name().left(15) + " "; - - out << qSetFieldWidth(5); - size_t ii, jj, kk; - mainGrid->ijkFromCellIndex(stimPlanContributingEclipseCells[i], &ii, &jj, &kk); - out << ii + 1; - out << jj + 1; - out << kk + 1; - - out << qSetFieldWidth(10); - out << stimPlanContributingEclipseCells[i]; - - out << qSetFieldWidth(5); - size_t spi = stimPlanCell.getI(); - size_t spj = stimPlanCell.getJ(); - - out << spi; - out << spj; - - out << qSetFieldWidth(10); - out << QString::number(stimPlanContributingEclipseCellTransmissibilities[i], 'e', 3); - - out << "\n"; - } - - //TODO: add RigFractureStimPlanCellData to m_StimPlanCellsFractureData i RigFracture??? - } - } - return; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RifFractureExportTools::printStimPlanFractureTrans(const std::vector& fractures, RimEclipseCase* caseToApply, QTextStream &out) -{ - - double cDarcyInCorrectUnit = caseToApply->eclipseCaseData()->darchysValue(); - - out << "StimPlan cells' fracture transmissibility \n"; - - out << qSetFieldWidth(4); - out << "-- "; - - out << qSetFieldWidth(5); - out << "SP i"; - out << "SP j"; - - out << qSetFieldWidth(10); - out << "Tf_hor"; - out << "Tf_vert"; - - out << "\n"; - - if (fractures.size() < 1) return; - RimFracture* fracture = fractures[0]; - - RimStimPlanFractureTemplate* fracTemplateStimPlan; - if (dynamic_cast(fracture->fractureTemplate())) - { - fracTemplateStimPlan = dynamic_cast(fracture->fractureTemplate()); - } - else return; - - std::vector stimPlanCells = fracTemplateStimPlan->fractureGrid()->fractureCells(); - - for (RigFractureCell stimPlanCell : stimPlanCells) - { - if (stimPlanCell.getConductivtyValue() < 1e-7) - { - //If conductivity in stimPlanCell is 0, contributions might not be relevant... - continue; - } - - double verticalTrans = RigFractureTransmissibilityEquations::centerToEdgeFractureCellTrans(stimPlanCell.getConductivtyValue(), - stimPlanCell.cellSizeX(), - stimPlanCell.cellSizeZ(), - cDarcyInCorrectUnit); - double horizontalTrans = RigFractureTransmissibilityEquations::centerToEdgeFractureCellTrans(stimPlanCell.getConductivtyValue(), - stimPlanCell.cellSizeZ(), - stimPlanCell.cellSizeX(), - cDarcyInCorrectUnit); - - out << qSetFieldWidth(5); - size_t spi = stimPlanCell.getI(); - size_t spj = stimPlanCell.getJ(); - - out << spi; - out << spj; - - out << qSetFieldWidth(10); - out << QString::number(verticalTrans, 'e', 3); - out << QString::number(horizontalTrans, 'e', 3); - - out << "\n"; - } - -return; - -} -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RifFractureExportTools::printBackgroundDataHeaderLine(QTextStream & out) -{ - out << "-- Background data for calculation" << "\n\n"; - - - //Write header line - out << qSetFieldWidth(4); - out << "--"; - out << qSetFieldWidth(12); - out << "Well "; - - out << qSetFieldWidth(16); - out << "Fracture "; - - out << qSetFieldWidth(5); - out << "i"; - out << "j"; - out << "k"; - - out << qSetFieldWidth(12); - out << "Ax"; - out << "Ay"; - out << "Az"; - out << "TotArea"; - - out << "skinfac"; - out << "FracLen"; - - out << qSetFieldWidth(10); - out << "DX"; - out << "DY"; - out << "DZ"; - - out << qSetFieldWidth(12); - out << "PermX"; - out << "PermY"; - out << "PermZ"; - - out << qSetFieldWidth(8); - out << "NTG"; - - out << qSetFieldWidth(12); - out << "T_x"; - out << "T_y"; - out << "T_z"; - - out << qSetFieldWidth(15); - out << "Transm"; - - out << qSetFieldWidth(20); - out << "Status"; - - out << "\n"; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RifFractureExportTools::printBackgroundData(QTextStream & out, const QString& wellName, RimFracture* fracture, const RigMainGrid* mainGrid, RigFracturedEclipseCellExportData &fracData) -{ - out << qSetFieldWidth(4); - out << "-- "; - - out << qSetFieldWidth(12); - - out << wellName + " "; - out << qSetFieldWidth(16); - out << fracture->name().left(15) + " "; - - - 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 << qSetFieldWidth(12); - //Use f for float, e for exponent float and g for best choice of these two. - out << QString::number(fracData.projectedAreas.x(), 'g', 4); - out << QString::number(fracData.projectedAreas.y(), 'g', 4); - out << QString::number(fracData.projectedAreas.z(), 'g', 4); - out << QString::number(fracData.totalArea, 'g', 4); - - out << QString::number(fracData.skinFactor, 'f', 2); - out << QString::number(fracData.fractureLenght, 'g', 3); - - out << qSetFieldWidth(10); - out << QString::number(fracData.cellSizes.x(), 'f', 2); - out << QString::number(fracData.cellSizes.y(), 'f', 2); - out << QString::number(fracData.cellSizes.z(), 'f', 2); - - out << qSetFieldWidth(12); - out << QString::number(fracData.permeabilities.x(), 'e', 3); - out << QString::number(fracData.permeabilities.y(), 'e', 3); - out << QString::number(fracData.permeabilities.z(), 'e', 3); - - out << qSetFieldWidth(8); - out << QString::number(fracData.NTG, 'f', 2); - - out << qSetFieldWidth(12); - out << QString::number(fracData.transmissibilities.x(), 'e', 3); - out << QString::number(fracData.transmissibilities.y(), 'e', 3); - out << QString::number(fracData.transmissibilities.z(), 'e', 3); - - out << qSetFieldWidth(15); - out << QString::number(fracData.transmissibility, 'e', 3); - - if (!fracData.cellIsActive) - { - out << qSetFieldWidth(20); - out << " INACTIVE CELL "; - } - - else if (fracData.cellIsActive && fracData.transmissibility > 0) - { - out << qSetFieldWidth(20); - out << " ACTIVE CELL "; - } - else - { - out << qSetFieldWidth(20); - out << " INVALID DATA "; - } - - out << "\n"; - -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RifFractureExportTools::printTransmissibilityFractureToWell(const std::vector& fractures, QTextStream &out, RimEclipseCase* caseToApply) -{ - out << "-- Transmissibility From Fracture To Well \n"; - - out << qSetFieldWidth(12); - out << "Well name "; - - out << qSetFieldWidth(16); - out << "Fracture name "; - out << "Inflow type "; - - out << qSetFieldWidth(5); - out << " i "; - out << " j "; - - out << "Tw"; - out << "\n"; - - for (RimFracture* fracture : fractures) - { - out << qSetFieldWidth(12); - RimEclipseWell* simWell = nullptr; - RimWellPath* wellPath = 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().left(15) + " "; - - - if (fracture->fractureTemplate()->orientationType == RimFractureTemplate::ALONG_WELL_PATH) - { - out << "Linear inflow"; - out << qSetFieldWidth(5); - - RimStimPlanFractureTemplate* fracTemplateStimPlan; - if (dynamic_cast(fracture->fractureTemplate())) - { - fracTemplateStimPlan = dynamic_cast(fracture->fractureTemplate()); - } - else continue; - - //TODO: Can be removed when implementation of dip angle is more general: - RimSimWellFracture* simWellFrac; - if (dynamic_cast(fracture)) - { - simWellFrac = dynamic_cast(fracture); - } - else continue; - - double wellDip = simWellFrac->wellDipAtFracturePosition(); - - double perforationLengthVert = fracture->perforationLength * cos(wellDip); - double perforationLengthHor = fracture->perforationLength * sin(wellDip); - - std::pair wellCenterStimPlanCellIJ = fracTemplateStimPlan->fractureGrid()->fractureCellAtWellCenter(); - out << qSetFieldWidth(5); - out << wellCenterStimPlanCellIJ.first; - out << wellCenterStimPlanCellIJ.second; - - - //RigStimPlanCell* stimPlanCell = fracTemplateStimPlan->getStimPlanCellAtIJ(wellCenterStimPlanCellIJ.first, wellCenterStimPlanCellIJ.second); - const RigFractureCell& stimPlanCell = fracTemplateStimPlan->fractureGrid()->cellFromIndex(fracTemplateStimPlan->fractureGrid()->getGlobalIndexFromIJ(wellCenterStimPlanCellIJ.first, wellCenterStimPlanCellIJ.second)); - - double linTransInStimPlanCell = RigFractureTransmissibilityEquations::fractureCellToWellLinearTrans(stimPlanCell.getConductivtyValue(), - stimPlanCell.cellSizeX(), - stimPlanCell.cellSizeZ(), - perforationLengthVert, - perforationLengthHor, - fracture->perforationEfficiency, - fracture->fractureTemplate()->skinFactor(), - caseToApply->eclipseCaseData()->darchysValue()); - - out << qSetFieldWidth(10); - out << QString::number(linTransInStimPlanCell, 'f', 2); - out << "\n"; - } - - - if (fracture->fractureTemplate()->orientationType == RimFractureTemplate::TRANSVERSE_WELL_PATH - || fracture->fractureTemplate()->orientationType == RimFractureTemplate::AZIMUTH) - { - out << "Radial inflow"; - - RimStimPlanFractureTemplate* fracTemplateStimPlan; - if (dynamic_cast(fracture->fractureTemplate())) - { - fracTemplateStimPlan = dynamic_cast(fracture->fractureTemplate()); - } - else continue; - - std::pair wellCenterStimPlanCellIJ = fracTemplateStimPlan->fractureGrid()->fractureCellAtWellCenter(); - out << qSetFieldWidth(5); - out << wellCenterStimPlanCellIJ.first; - out << wellCenterStimPlanCellIJ.second; - - const RigFractureCell& stimPlanCell = fracTemplateStimPlan->fractureGrid()->cellFromIndex(fracTemplateStimPlan->fractureGrid()->getGlobalIndexFromIJ(wellCenterStimPlanCellIJ.first, wellCenterStimPlanCellIJ.second)); - - double radTransInStimPlanCell = RigFractureTransmissibilityEquations::fractureCellToWellRadialTrans(stimPlanCell.getConductivtyValue(), - stimPlanCell.cellSizeX(), - stimPlanCell.cellSizeZ(), - fracture->wellRadius(), - fracture->fractureTemplate()->skinFactor(), - caseToApply->eclipseCaseData()->darchysValue()); - - out << qSetFieldWidth(10); - out << QString::number(radTransInStimPlanCell, 'f', 2); - out << "\n"; - - - } - - - - } - - out << "\n"; -} diff --git a/ApplicationCode/FileInterface/RifFractureExportTools.h b/ApplicationCode/FileInterface/RifFractureExportTools.h deleted file mode 100644 index ddbeb702dc..0000000000 --- a/ApplicationCode/FileInterface/RifFractureExportTools.h +++ /dev/null @@ -1,82 +0,0 @@ -///////////////////////////////////////////////////////////////////////////////// -// -// Copyright (C) 2017- Statoil ASA -// -// ResInsight is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. -// -// See the GNU General Public License at -// for more details. -// -///////////////////////////////////////////////////////////////////////////////// - -#pragma once - -#include "cvfBase.h" -#include "cvfObject.h" -#include "cvfLibCore.h" - -#include "ert/ecl/ecl_kw.h" - -#include -#include - - -class QFile; -class QTextStream; -class RigFracturedEclipseCellExportData; -class RigMainGrid; -class RimEclipseCase; -class RimEclipseWell; -class RimFracture; -class RimFracture; -class RimWellPath; - -//================================================================================================== -// -// Class for access to Eclipse "keyword" files using libecl -// -//================================================================================================== -class RifFractureExportTools -{ -public: - static bool exportFracturesToEclipseDataInputFile(const QString& fileName, - const std::vector& fractures, - RimEclipseCase* caseToApply); - -private: - - static void printCOMPDATvalues(QTextStream & out, - double transmissibility, - size_t i, size_t j, size_t k, - const QString& fractureName, - double skinFactor, - const QString& wellName); - - static void printStimPlanCellsMatrixTransContributions(const std::vector& fractures, - RimEclipseCase* caseToApply, - QTextStream &out, - const QString& wellName, - const RigMainGrid* mainGrid); - static void printStimPlanFractureTrans(const std::vector& fractures, - RimEclipseCase* caseToApply, - QTextStream &out); - static void printTransmissibilityFractureToWell(const std::vector& fractures, - QTextStream &out, - RimEclipseCase* caseToApply); - - - static void printBackgroundDataHeaderLine(QTextStream & out); - - static void printBackgroundData(QTextStream & out, - const QString& wellName, - RimFracture* fracture, - const RigMainGrid* mainGrid, - RigFracturedEclipseCellExportData &fracData); -}; diff --git a/ApplicationCode/ProjectDataModel/Completions/RimStimPlanFractureTemplate.cpp b/ApplicationCode/ProjectDataModel/Completions/RimStimPlanFractureTemplate.cpp index f12a801d74..7851e91f9e 100644 --- a/ApplicationCode/ProjectDataModel/Completions/RimStimPlanFractureTemplate.cpp +++ b/ApplicationCode/ProjectDataModel/Completions/RimStimPlanFractureTemplate.cpp @@ -22,7 +22,6 @@ #include "RiaLogging.h" #include "RigStimPlanFractureDefinition.h" -#include "RigFractureTransCalc.h" #include "RigFractureGrid.h" #include "RimEclipseView.h" diff --git a/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake b/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake index 414c88b647..12c8bf304e 100644 --- a/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake +++ b/ApplicationCode/ReservoirDataModel/CMakeLists_files.cmake @@ -57,13 +57,11 @@ ${CEE_CURRENT_LIST_DIR}RigSummaryCaseData.h ${CEE_CURRENT_LIST_DIR}RigLasFileExporter.h ${CEE_CURRENT_LIST_DIR}RigSimulationWellCoordsAndMD.h ${CEE_CURRENT_LIST_DIR}RigFishbonesGeometry.h -${CEE_CURRENT_LIST_DIR}RigFractureTransCalc.h ${CEE_CURRENT_LIST_DIR}RigFractureTransmissibilityEquations.h ${CEE_CURRENT_LIST_DIR}RigWellPathStimplanIntersector.h ${CEE_CURRENT_LIST_DIR}RigTesselatorTools.h ${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.h ${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.h -${CEE_CURRENT_LIST_DIR}RigStimPlanUpscalingCalc.h ${CEE_CURRENT_LIST_DIR}RigFractureGrid.h ${CEE_CURRENT_LIST_DIR}RigFractureCell.h ${CEE_CURRENT_LIST_DIR}RigWellPathIntersectionTools.h @@ -119,13 +117,11 @@ ${CEE_CURRENT_LIST_DIR}RigSummaryCaseData.cpp ${CEE_CURRENT_LIST_DIR}RigLasFileExporter.cpp ${CEE_CURRENT_LIST_DIR}RigSimulationWellCoordsAndMD.cpp ${CEE_CURRENT_LIST_DIR}RigFishbonesGeometry.cpp -${CEE_CURRENT_LIST_DIR}RigFractureTransCalc.cpp ${CEE_CURRENT_LIST_DIR}RigFractureTransmissibilityEquations.cpp ${CEE_CURRENT_LIST_DIR}RigWellPathStimplanIntersector.cpp ${CEE_CURRENT_LIST_DIR}RigTesselatorTools.cpp ${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.cpp ${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.cpp -${CEE_CURRENT_LIST_DIR}RigStimPlanUpscalingCalc.cpp ${CEE_CURRENT_LIST_DIR}RigFractureGrid.cpp ${CEE_CURRENT_LIST_DIR}RigFractureCell.cpp ${CEE_CURRENT_LIST_DIR}RigWellPathIntersectionTools.cpp diff --git a/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.cpp b/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.cpp index 5570ae382c..76f330752c 100644 --- a/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.cpp +++ b/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.cpp @@ -22,7 +22,6 @@ #include "RigCellGeometryTools.h" #include "RigEclipseCaseData.h" #include "RigFractureCell.h" -#include "RigFractureTransCalc.h" #include "RigFractureTransmissibilityEquations.h" #include "RigMainGrid.h" #include "RigResultAccessorFactory.h" @@ -135,7 +134,7 @@ void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsM mainGrid->cellCornerVertices(fracCell, hexCorners.data()); std::vector > planeCellPolygons; - bool isPlanIntersected = RigFractureTransCalc::planeCellIntersectionPolygons(hexCorners.data(), m_fractureTransform, planeCellPolygons); + bool isPlanIntersected = planeCellIntersectionPolygons(hexCorners.data(), m_fractureTransform, planeCellPolygons); if (!isPlanIntersected || planeCellPolygons.size() == 0) continue; cvf::Vec3d localX; @@ -242,3 +241,28 @@ std::vector RigEclipseToStimPlanCellTransmissibilityCalculator::getPoten return cellIndices; } + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigEclipseToStimPlanCellTransmissibilityCalculator::planeCellIntersectionPolygons(cvf::Vec3d hexCorners[8], + cvf::Mat4f transformMatrixForPlane, + std::vector > & polygons) +{ + bool isCellIntersected = false; + + cvf::Plane fracturePlane; + fracturePlane.setFromPointAndNormal(static_cast(transformMatrixForPlane.translation()), + static_cast(transformMatrixForPlane.col(2))); + + //Find line-segments where cell and fracture plane intersects + std::list > intersectionLineSegments; + + isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments); + + RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons); + + return isCellIntersected; +} + diff --git a/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.h b/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.h index 66250adf1e..f9f8e8264c 100644 --- a/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.h +++ b/ApplicationCode/ReservoirDataModel/RigEclipseToStimPlanCellTransmissibilityCalculator.h @@ -46,6 +46,11 @@ private: void calculateStimPlanCellsMatrixTransmissibility(); std::vector getPotentiallyFracturedCellsForPolygon(std::vector polygon); + + static bool planeCellIntersectionPolygons(cvf::Vec3d hexCorners[8], + cvf::Mat4f transformMatrixForPlane, + std::vector > & polygons); + private: const RimEclipseCase* m_case; double m_cDarcy; diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp deleted file mode 100644 index 9eadd793f3..0000000000 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp +++ /dev/null @@ -1,301 +0,0 @@ -///////////////////////////////////////////////////////////////////////////////// -// -// Copyright (C) 2017 Statoil ASA -// -// ResInsight is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. -// -// See the GNU General Public License at -// for more details. -// -///////////////////////////////////////////////////////////////////////////////// - -#include "RigFractureTransCalc.h" -#include "RimFractureTemplate.h" -#include "RigEclipseCaseData.h" -#include "RimEclipseCase.h" -#include "RiaLogging.h" -#include "QString" -#include "RimReservoirCellResultsStorage.h" -#include "RigResultAccessorFactory.h" -#include "RimFracture.h" -#include "cvfGeometryTools.h" -#include "RigCellGeometryTools.h" -#include "RigActiveCellInfo.h" -#include "RimStimPlanFractureTemplate.h" - -#include -#include "RimEllipseFractureTemplate.h" -#include "cafAppEnum.h" -#include "RigCell.h" -#include "RigMainGrid.h" -#include "cvfMath.h" -#include "RiaEclipseUnitTools.h" -#include "RigFractureCell.h" -#include //Used for log -#include "RiaApplication.h" -#include "RimEclipseView.h" -#include "RigFractureTransmissibilityEquations.h" - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -RigFractureTransCalc::RigFractureTransCalc(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 = RiaEclipseUnitTools::UNITS_METRIC; - } - else if (caseUnit == RigEclipseCaseData::UNITS_FIELD) - { - RiaLogging::debug(QString("Calculating transmissibilities in field units")); - m_unitForCalculation = RiaEclipseUnitTools::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 = RiaEclipseUnitTools::UNITS_METRIC; - } - - -} - - - -//-------------------------------------------------------------------------------------------------- -/// TODO: Document equation -//-------------------------------------------------------------------------------------------------- -std::vector RigFractureTransCalc::computeTransmissibilityFromPolygonWithInfiniteConductivityInFracture() -{ - if (m_fracture->fractureTemplate()->conductivityType == RimFractureTemplate::FINITE_CONDUCTIVITY) - { - RiaLogging::warning(QString("Transimssibility for finite conductity in fracture not yet implemented.")); - RiaLogging::warning(QString("Performing calculation for infinite conductivity instead.")); - } - - RigEclipseCaseData* eclipseCaseData = m_case->eclipseCaseData(); - double cDarchy = eclipseCaseData->darchysValue(); - - RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS; - RimReservoirCellResultsStorage* gridCellResults = m_case->results(porosityModel); - - size_t scalarSetIndex; - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "DX"); - cvf::ref dataAccessObjectDx = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DX"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "DY"); - cvf::ref dataAccessObjectDy = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DY"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "DZ"); - cvf::ref dataAccessObjectDz = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DZ"); //assuming 0 time step and main grid (so grid index =0) - - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "PERMX"); - cvf::ref dataAccessObjectPermX = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMX"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "PERMY"); - cvf::ref dataAccessObjectPermY = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMY"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "PERMZ"); - cvf::ref dataAccessObjectPermZ = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMZ"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RiaDefines::STATIC_NATIVE, "NTG"); - cvf::ref dataAccessObjectNTG = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "NTG"); //assuming 0 time step and main grid (so grid index =0) - - RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); - - std::vector fracDataVec; - std::vector fracCells = m_fracture->getPotentiallyFracturedCells(eclipseCaseData->mainGrid()); - - for (size_t fracCell : fracCells) - { - bool cellIsActive = activeCellInfo->isActive(fracCell); - - double permX = dataAccessObjectPermX->cellScalarGlobIdx(fracCell); - double permY = dataAccessObjectPermY->cellScalarGlobIdx(fracCell); - double permZ = dataAccessObjectPermZ->cellScalarGlobIdx(fracCell); - - double dx = dataAccessObjectDx->cellScalarGlobIdx(fracCell); - double dy = dataAccessObjectDy->cellScalarGlobIdx(fracCell); - double dz = dataAccessObjectDz->cellScalarGlobIdx(fracCell); - - double NTG = dataAccessObjectNTG->cellScalarGlobIdx(fracCell); - - const RigMainGrid* mainGrid = m_case->eclipseCaseData()->mainGrid(); - std::array hexCorners; - mainGrid->cellCornerVertices(fracCell, hexCorners.data()); - - std::vector > planeCellPolygons; - bool isPlanIntersected = planeCellIntersectionPolygons(hexCorners.data(), m_fracture->transformMatrix(), planeCellPolygons); - if (!isPlanIntersected || planeCellPolygons.size() == 0) continue; - - cvf::Vec3d localX; - cvf::Vec3d localY; - cvf::Vec3d localZ; - RigCellGeometryTools::findCellLocalXYZ(hexCorners, localX, localY, localZ); - - //Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon already is located) - cvf::Mat4f invertedTransMatrix = m_fracture->transformMatrix().getInverted(); - for (std::vector & planeCellPolygon : planeCellPolygons) - { - for (cvf::Vec3d& v : planeCellPolygon) - { - v.transformPoint(static_cast(invertedTransMatrix)); - } - } - - cvf::Vec3d localZinFracPlane; - localZinFracPlane = localZ; - localZinFracPlane.transformVector(static_cast(invertedTransMatrix)); - cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO; - directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1)); - directionOfLength.normalize(); - - RigFracturedEclipseCellExportData fracData; - fracData.reservoirCellIndex = fracCell; - - std::vector fracPolygon = m_fracture->fractureTemplate()->fractureBorderPolygon(m_unitForCalculation); - - std::vector fracPolygonDouble; - for (auto v : fracPolygon) fracPolygonDouble.push_back(static_cast(v)); - - std::vector > polygonsDescribingFractureInCell; - cvf::Vec3d areaVector; - - for (std::vector planeCellPolygon : planeCellPolygons) - { - std::vector >clippedPolygons = RigCellGeometryTools::intersectPolygons(planeCellPolygon, fracPolygonDouble); - for (std::vector clippedPolygon : clippedPolygons) - { - polygonsDescribingFractureInCell.push_back(clippedPolygon); - } - } - - double area; - std::vector areaOfFractureParts; - double length; - std::vector lengthXareaOfFractureParts; - double Ax = 0.0; - double Ay = 0.0; - double Az = 0.0; - - for (std::vector fracturePartPolygon : polygonsDescribingFractureInCell) - { - areaVector = cvf::GeometryTools::polygonAreaNormal3D(fracturePartPolygon); - area = areaVector.length(); - areaOfFractureParts.push_back(area); - length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon); - lengthXareaOfFractureParts.push_back(length * area); - - cvf::Plane fracturePlane; - cvf::Mat4f m = m_fracture->transformMatrix(); - bool isCellIntersected = false; - - fracturePlane.setFromPointAndNormal(static_cast(m.translation()), - static_cast(m.col(2))); - - Ax += abs(area*(fracturePlane.normal().dot(localY))); - Ay += abs(area*(fracturePlane.normal().dot(localX))); - Az += abs(area*(fracturePlane.normal().dot(localZ))); - //TODO: resulting values have only been checked for vertical fracture... - } - - double fractureArea = 0.0; - for (double area : areaOfFractureParts) fractureArea += area; - - double totalAreaXLength = 0.0; - for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea; - - double fractureAreaWeightedlength = totalAreaXLength / fractureArea; - double skinfactor = m_fracture->fractureTemplate()->skinFactor; - - double transmissibility_X = RigFractureTransmissibilityEquations::matrixToFractureTrans(permY, NTG, Ay, dx, skinfactor, fractureAreaWeightedlength, cDarchy); - double transmissibility_Y = RigFractureTransmissibilityEquations::matrixToFractureTrans(permX, NTG, Ax, dy, skinfactor, fractureAreaWeightedlength, cDarchy); - double transmissibility_Z = RigFractureTransmissibilityEquations::matrixToFractureTrans(permZ, 1.0, Az, dz, skinfactor, fractureAreaWeightedlength, cDarchy); - - double transmissibility = sqrt(transmissibility_X * transmissibility_X - + transmissibility_Y * transmissibility_Y - + transmissibility_Z * transmissibility_Z); - - fracData.transmissibility = transmissibility; - fracData.transmissibilities = cvf::Vec3d(transmissibility_X, transmissibility_Y, transmissibility_Z); - - fracData.totalArea = fractureArea; - fracData.projectedAreas = cvf::Vec3d(Ax, Ay, Az); - fracData.fractureLenght = fractureAreaWeightedlength; - - fracData.cellSizes = cvf::Vec3d(dx, dy, dz); - fracData.permeabilities = cvf::Vec3d(permX, permY, permZ); - fracData.NTG = NTG; - fracData.skinFactor = skinfactor; - fracData.cellIsActive = cellIsActive; - - //Since we loop over all potentially fractured cells, we only keep FractureData for cells where fracture have an non-zero area. - if (fractureArea > 1e-5) - { - fracDataVec.push_back(fracData); - } - - } - - return fracDataVec; -} - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -bool RigFractureTransCalc::planeCellIntersectionPolygons(cvf::Vec3d hexCorners[8], - cvf::Mat4f transformMatrixForPlane, - std::vector > & polygons) -{ - - //Lage static func - input: transform-matrix for plan, hexcorners for celle - - bool isCellIntersected = false; - - cvf::Plane fracturePlane; - fracturePlane.setFromPointAndNormal(static_cast(transformMatrixForPlane.translation()), - static_cast(transformMatrixForPlane.col(2))); - - //Find line-segments where cell and fracture plane intersects - std::list > intersectionLineSegments; - - isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments); - - RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons); - - return isCellIntersected; -} - - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RigFractureTransCalc::convertConductivtyValue(double Kw, RiaEclipseUnitTools::UnitSystem fromUnit, RiaEclipseUnitTools::UnitSystem toUnit) -{ - - if (fromUnit == toUnit) return Kw; - - else if (fromUnit == RiaEclipseUnitTools::UNITS_METRIC && toUnit == RiaEclipseUnitTools::UNITS_FIELD) - { - return RiaEclipseUnitTools::meterToFeet(Kw); - } - else if (fromUnit == RiaEclipseUnitTools::UNITS_METRIC && toUnit == RiaEclipseUnitTools::UNITS_FIELD) - { - return RiaEclipseUnitTools::feetToMeter(Kw); - } - - return cvf::UNDEFINED_DOUBLE; -} - diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h deleted file mode 100644 index be625ebab1..0000000000 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h +++ /dev/null @@ -1,96 +0,0 @@ -///////////////////////////////////////////////////////////////////////////////// -// -// Copyright (C) 2017 Statoil ASA -// -// ResInsight is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. -// -// See the GNU General Public License at -// for more details. -// -///////////////////////////////////////////////////////////////////////////////// - -#pragma once - - -#include "RiaEclipseUnitTools.h" - -#include "cafAppEnum.h" - -#include "cvfBase.h" -#include "cvfMath.h" -#include "cvfVector3.h" -#include "cvfMatrix4.h" -#include "cvfPlane.h" - -#include -#include - -//================================================================================================== -/// -//================================================================================================== - -class RigFracturedEclipseCellExportData -{ -public: - RigFracturedEclipseCellExportData() {}; - - // Compdat export data - size_t reservoirCellIndex; - double transmissibility; // Total cell to well transmissibility finally used in COMPDAT keyword - bool cellIsActive; - - // General intermediate results - double NTG; - cvf::Vec3d permeabilities; - double skinFactor; - - // Elipse fracture related values - cvf::Vec3d transmissibilities; //matrixToFractureTransmissibilitiesXYZ - double totalArea; // Elipse cell overlap area - double fractureLenght; - cvf::Vec3d projectedAreas; - - cvf::Vec3d cellSizes; - - //TODO: Used for upscaling - should be moved? - double upscaledStimPlanValueHA; - double upscaledStimPlanValueAH; -}; - - -class RimFracture; -class RimEclipseCase; -class RigFractureCell; -class RimStimPlanFractureTemplate; - -//================================================================================================== -/// -//================================================================================================== -class RigFractureTransCalc -{ -public: - explicit RigFractureTransCalc(RimEclipseCase* caseToApply, RimFracture* fracture); - - // Calculations based on fracture polygon and eclipse grid cells - std::vector computeTransmissibilityFromPolygonWithInfiniteConductivityInFracture(); - - static bool planeCellIntersectionPolygons(cvf::Vec3d hexCorners[8], - cvf::Mat4f transformMatrixForPlane, - std::vector > & polygons); - -private: - static double convertConductivtyValue(double Kw, RiaEclipseUnitTools::UnitSystem fromUnit, RiaEclipseUnitTools::UnitSystem toUnit); - -private: - RimEclipseCase* m_case; - RimFracture* m_fracture; - RiaEclipseUnitTools::UnitSystem m_unitForCalculation; - -}; diff --git a/ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.cpp b/ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.cpp deleted file mode 100644 index 448daf8f6c..0000000000 --- a/ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.cpp +++ /dev/null @@ -1,400 +0,0 @@ -#include "RiaLogging.h" - -#include "RigCellGeometryTools.h" -#include "RigEclipseCaseData.h" -#include "RigFractureCell.h" -#include "RigFractureGrid.h" -#include "RigFractureTransCalc.h" -#include "RigMainGrid.h" -#include "RigStimPlanUpscalingCalc.h" - -#include "RimEclipseCase.h" -#include "RimStimPlanFractureTemplate.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 = RiaEclipseUnitTools::UNITS_METRIC; - } - else if (caseUnit == RigEclipseCaseData::UNITS_FIELD) - { - RiaLogging::debug(QString("Calculating transmissibilities in field units")); - m_unitForCalculation = RiaEclipseUnitTools::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 = RiaEclipseUnitTools::UNITS_METRIC; - } - -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::pair RigStimPlanUpscalingCalc::flowAcrossLayersUpscaling(QString resultName, QString resultUnit, size_t timeStepIndex, RiaEclipseUnitTools::UnitSystem unitSystem, size_t eclipseCellIndex) -{ - RimStimPlanFractureTemplate* fracTemplateStimPlan; - if (dynamic_cast(m_fracture->fractureTemplate())) - { - fracTemplateStimPlan = dynamic_cast(m_fracture->fractureTemplate()); - } - else return std::make_pair(cvf::UNDEFINED_DOUBLE, cvf::UNDEFINED_DOUBLE); - - std::vector stimPlanCells = fracTemplateStimPlan->fractureGrid()->fractureCells(); - - - - cvf::Vec3d localX, localY, localZ; //Not used in calculation here, but needed for function to find planCellPolygons - std::vector > 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 & planeCellPolygon : planeCellPolygons) - { - for (cvf::Vec3d& v : planeCellPolygon) - { - v.transformPoint(static_cast(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 fracPolygon = m_fracture->fractureTemplate()->fractureBorderPolygon(unitSystem); - std::vector > polygonsDescribingFractureInCell; - - std::vector upscaledConductivitiesHA; - std::vector upscaledConductivitiesAH; - - for (std::vector 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 stimPlanCells, std::vector planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers) -{ - std::vector DcolSum; - std::vector lavgCol; - std::vector CondHarmCol; - - for (size_t j = 0; j < fracTemplateStimPlan->fractureGrid()->iCellCount(); j++) - { - std::vector conductivitiesInStimPlanCells; - std::vector lengthsLiOfStimPlanCol; - std::vector heightsDiOfStimPlanCells; - - std::vector stimPlanCellsCol = getColOfStimPlanCells(stimPlanCells, j); - for (RigFractureCell* stimPlanCell : stimPlanCellsCol) - { - if (stimPlanCell->getConductivtyValue() > 1e-7) - { - std::vector >clippedStimPlanPolygons = RigCellGeometryTools::intersectPolygons(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 (size_t 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 (size_t 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 stimPlanCells, std::vector planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers) -{ - std::vector DrowAvg; - std::vector liRowSum; - std::vector CondAritRow; - - for (size_t j = 0; j < fracTemplateStimPlan->fractureGrid()->jCellCount(); j++) - { - std::vector conductivitiesInStimPlanCells; - std::vector lengthsLiOfStimPlanCol; - std::vector heightsDiOfStimPlanCells; - - std::vector stimPlanCellsCol = getRowOfStimPlanCells(stimPlanCells, j); - for (RigFractureCell* stimPlanCell : stimPlanCellsCol) - { - if (stimPlanCell->getConductivtyValue() > 1e-7) - { - std::vector >clippedStimPlanPolygons = RigCellGeometryTools::intersectPolygons(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 (size_t 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 (size_t 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 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 RigStimPlanUpscalingCalc::computeUpscaledPropertyFromStimPlan( QString resultName, QString resultUnit, size_t timeStepIndex) -{ - std::vector fracDataVec; - - RimStimPlanFractureTemplate* fracTemplateStimPlan; - if (dynamic_cast(m_fracture->fractureTemplate())) - { - fracTemplateStimPlan = dynamic_cast(m_fracture->fractureTemplate()); - } - else - { - return fracDataVec; - } - - std::vector > stimPlanCellsAsPolygons; - std::vector stimPlanParameterValues; - fracTemplateStimPlan->getStimPlanDataAsPolygonsAndValues(stimPlanCellsAsPolygons, stimPlanParameterValues, resultName, resultUnit, timeStepIndex); - - //TODO: A lot of common code with function above... Can be cleaned up...? - - std::vector 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) - { - //TODO: Lage ny classe for å holde upscaledData - RigFracturedEclipseCellExportData fracData; - fracData.reservoirCellIndex = fracCell; - - std::pair 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 RigStimPlanUpscalingCalc::getRowOfStimPlanCells(std::vector& allStimPlanCells, size_t i) -{ - std::vector stimPlanCellRow; - - for (RigFractureCell stimPlanCell : allStimPlanCells) - { - if (stimPlanCell.getI() == i) - { - stimPlanCellRow.push_back(&stimPlanCell); - } - } - - return stimPlanCellRow; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RigStimPlanUpscalingCalc::getColOfStimPlanCells(std::vector& allStimPlanCells, size_t j) -{ - std::vector stimPlanCellCol; - - for (RigFractureCell stimPlanCell : allStimPlanCells) - { - if (stimPlanCell.getJ() == j) - { - stimPlanCellCol.push_back(&stimPlanCell); - } - } - - return stimPlanCellCol; -} - -#include "RigStimPlanFractureDefinition.h" - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RimStimPlanFractureTemplate::getStimPlanDataAsPolygonsAndValues(std::vector > &cellsAsPolygons, - std::vector ¶meterValues, - const QString& resultName, - const QString& unitName, - size_t timeStepIndex) -{ - std::vector< std::vector > propertyValuesAtTimeStep = m_stimPlanFractureDefinitionData->getMirroredDataAtTimeIndex(resultName, unitName, timeStepIndex); - - cellsAsPolygons.clear(); - parameterValues.clear(); - - //TODO: Code partly copied from RivWellFracturePartMgr - can this be combined in some function? - std::vector depthCoordsAtNodes = m_stimPlanFractureDefinitionData->adjustedDepthCoordsAroundWellPathPosition(m_wellPathDepthAtFracture()); - std::vector xCoordsAtNodes = m_stimPlanFractureDefinitionData->getNegAndPosXcoords(); - - //Cells are around nodes instead of between nodes - std::vector xCoords; - for (int i = 0; i < xCoordsAtNodes.size() - 1; i++) xCoords.push_back((xCoordsAtNodes[i] + xCoordsAtNodes[i + 1]) / 2); - std::vector depthCoords; - for (int i = 0; i < depthCoordsAtNodes.size() - 1; i++) depthCoords.push_back((depthCoordsAtNodes[i] + depthCoordsAtNodes[i + 1]) / 2); - - for (int i = 0; i < xCoords.size() - 1; i++) - { - for (int j = 0; j < depthCoords.size() - 1; j++) - { - std::vector cellAsPolygon; - cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i]), static_cast(depthCoords[j]), 0.0)); - cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i + 1]), static_cast(depthCoords[j]), 0.0)); - cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i + 1]), static_cast(depthCoords[j + 1]), 0.0)); - cellAsPolygon.push_back(cvf::Vec3d(static_cast(xCoords[i]), static_cast(depthCoords[j + 1]), 0.0)); - cellsAsPolygons.push_back(cellAsPolygon); - //TODO: Values for both neg and pos x values... - parameterValues.push_back(propertyValuesAtTimeStep[j+1][i+1]); //TODO test that this value exsist... - - } - } -} diff --git a/ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.h b/ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.h deleted file mode 100644 index b025ff5c6c..0000000000 --- a/ApplicationCode/ReservoirDataModel/RigStimPlanUpscalingCalc.h +++ /dev/null @@ -1,29 +0,0 @@ -#pragma once - -#include "RimEclipseCase.h" -#include "RimFracture.h" -#include "RimStimPlanFractureTemplate.h" - - -class RigStimPlanUpscalingCalc -{ -public: - explicit RigStimPlanUpscalingCalc(RimEclipseCase* caseToApply, RimFracture* fracture); - - std::vector computeUpscaledPropertyFromStimPlan(QString resultName, QString resultUnit, size_t timeStepIndex); - -private: - std::pair flowAcrossLayersUpscaling(QString resultName, QString resultUnit, size_t timeStepIndex, RiaEclipseUnitTools::UnitSystem unitSystem, size_t eclipseCellIndex); - double computeHAupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector stimPlanCells, std::vector planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers); - double computeAHupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector stimPlanCells, std::vector planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers); - static double arithmeticAverage(std::vector values); - - static std::vector getRowOfStimPlanCells(std::vector& allStimPlanCells, size_t i); - static std::vector getColOfStimPlanCells(std::vector& allStimPlanCells, size_t j); - -private: - RimEclipseCase* m_case; - RimFracture* m_fracture; - RiaEclipseUnitTools::UnitSystem m_unitForCalculation; -}; -