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

Merge branch 'dev' into fishbones

Browse Source
This commit is contained in:
Bjørnar Grip Fjær 2017-06-01 09:02:38 +02:00
commit 1aedf92efc
33 changed files with 6868 additions and 177 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
ApplicationCode/Adm/LicenseInformation.txt
View File

@ -322,3 +322,46 @@ exceptions:
CRAVA is a software package for seismic inversion and conditioning of
geological reservoir models. CRAVA is copyrighted by the Norwegian
Computing Center and Statoil and licensed under GPLv3+.
===============================================================================
Notice for opm-flowdiagnostics and opm-flowdiagnostics-applications libraries
===============================================================================
Copyright 2016, 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
===============================================================================
Notice for the NightCharts code
===============================================================================
NightCharts
Copyright (C) 2010 by Alexander A. Avdonin, Artem N. Ivanov / ITGears Co.
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Please contact gordos.kund@gmail.com with any questions on this license.

183
ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.cpp
View File

@ -22,8 +22,10 @@
#include "RiaPreferences.h"
#include "RifEclipseSummaryAddress.h"
#include "RifReaderEclipseSummary.h"
#include "RigSingleWellResultsData.h"
#include "RigSummaryCaseData.h"
#include "RimEclipseResultCase.h"
#include "RimEclipseWell.h"
@ -33,7 +35,7 @@
#include "RimProject.h"
#include "RimSummaryCaseCollection.h"
#include "RimSummaryCurve.h"
#include "RimSummaryCurveFilter.h"
#include "RimSummaryCurveAppearanceCalculator.h"
#include "RimSummaryPlot.h"
#include "RimSummaryPlotCollection.h"
#include "RimView.h"
@ -94,67 +96,92 @@ void RicPlotProductionRateFeature::onActionTriggered(bool isChecked)
RimGridSummaryCase* gridSummaryCase = RicPlotProductionRateFeature::gridSummaryCaseForWell(well);
if (!gridSummaryCase) continue;
QString curveFilterText = "W*PR:";
QString description = "Well Production Rates : ";
RigSingleWellResultsData* wRes = well->wellResults();
if (wRes)
if (isInjector(well))
{
RimView* rimView = nullptr;
well->firstAncestorOrThisOfTypeAsserted(rimView);
int currentTimeStep = rimView->currentTimeStep();
if (wRes->hasWellResult(currentTimeStep))
{
const RigWellResultFrame& wrf = wRes->wellResultFrame(currentTimeStep);
if ( wrf.m_productionType == RigWellResultFrame::OIL_INJECTOR
|| wrf.m_productionType == RigWellResultFrame::GAS_INJECTOR
|| wrf.m_productionType == RigWellResultFrame::WATER_INJECTOR)
{
curveFilterText = "W*IR:";
description = "Well Injection Rates : ";
}
}
description = "Well Injection Rates : ";
}
curveFilterText += well->name();
description += well->name();
RimSummaryPlot* plot = new RimSummaryPlot();
summaryPlotColl->summaryPlots().push_back(plot);
description += well->name();
plot->setDescription(description);
if (isInjector(well))
{
RimSummaryCurveFilter* newCurveFilter = new RimSummaryCurveFilter();
plot->addCurveFilter(newCurveFilter);
// Left Axis
newCurveFilter->createCurves(gridSummaryCase, curveFilterText);
RimDefines::PlotAxis plotAxis = RimDefines::PLOT_AXIS_LEFT;
{
// Note : The parameter "WOIR" is probably never-existing, but we check for existence before creating curve
// Oil
QString parameterName = "WOIR";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledGreenColor(0));
}
{
// Water
QString parameterName = "WWIR";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledBlueColor(0));
}
{
// Gas
QString parameterName = "WGIR";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledRedColor(0));
}
}
else
{
// Left Axis
RimDefines::PlotAxis plotAxis = RimDefines::PLOT_AXIS_LEFT;
{
// Oil
QString parameterName = "WOPR";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledGreenColor(0));
}
{
// Water
QString parameterName = "WWPR";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledBlueColor(0));
}
{
// Gas
QString parameterName = "WGPR";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledRedColor(0));
}
}
// Right Axis
{
RimSummaryCurve* newCurve = new RimSummaryCurve();
plot->addCurve(newCurve);
RimDefines::PlotAxis plotAxis = RimDefines::PLOT_AXIS_RIGHT;
newCurve->setSummaryCase(gridSummaryCase);
{
QString parameterName = "WTHP";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledNoneRGBBrColor(0));
}
RifEclipseSummaryAddress addr( RifEclipseSummaryAddress::SUMMARY_WELL,
"WBHP",
-1,
-1,
"",
well->name().toStdString(),
-1,
"",
-1,
-1,
-1);
newCurve->setSummaryAddress(addr);
newCurve->setYAxis(RimDefines::PlotAxis::PLOT_AXIS_RIGHT);
{
QString parameterName = "WBHP";
RicPlotProductionRateFeature::addSummaryCurve(plot, well, gridSummaryCase, parameterName,
plotAxis, RimSummaryCurveAppearanceCalculator::cycledNoneRGBBrColor(1));
}
}
summaryPlotColl->updateConnectedEditors();
@ -210,3 +237,71 @@ RimGridSummaryCase* RicPlotProductionRateFeature::gridSummaryCaseForWell(RimEcli
return nullptr;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RicPlotProductionRateFeature::isInjector(RimEclipseWell* well)
{
RigSingleWellResultsData* wRes = well->wellResults();
if (wRes)
{
RimView* rimView = nullptr;
well->firstAncestorOrThisOfTypeAsserted(rimView);
int currentTimeStep = rimView->currentTimeStep();
if (wRes->hasWellResult(currentTimeStep))
{
const RigWellResultFrame& wrf = wRes->wellResultFrame(currentTimeStep);
if ( wrf.m_productionType == RigWellResultFrame::OIL_INJECTOR
|| wrf.m_productionType == RigWellResultFrame::GAS_INJECTOR
|| wrf.m_productionType == RigWellResultFrame::WATER_INJECTOR)
{
return true;
}
}
}
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RimSummaryCurve* RicPlotProductionRateFeature::addSummaryCurve( RimSummaryPlot* plot, const RimEclipseWell* well,
RimGridSummaryCase* gridSummaryCase, const QString& vectorName,
RimDefines::PlotAxis plotAxis, const cvf::Color3f& color)
{
CVF_ASSERT(plot);
CVF_ASSERT(gridSummaryCase);
CVF_ASSERT(well);
RifEclipseSummaryAddress addr(RifEclipseSummaryAddress::SUMMARY_WELL,
vectorName.toStdString(),
-1,
-1,
"",
well->name().toStdString(),
-1,
"",
-1,
-1,
-1);
if (!gridSummaryCase->caseData()->summaryReader()->hasAddress(addr))
{
return nullptr;
}
RimSummaryCurve* newCurve = new RimSummaryCurve();
plot->addCurve(newCurve);
newCurve->setSummaryCase(gridSummaryCase);
newCurve->setSummaryAddress(addr);
newCurve->setColor(color);
newCurve->setYAxis(plotAxis);
return newCurve;
}

9
ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.h
View File

@ -21,9 +21,12 @@
#include "cafCmdFeature.h"
#include "RimFlowDiagSolution.h"
#include "RimDefines.h"
class RimGridSummaryCase;
class RimEclipseWell;
class RimSummaryPlot;
class RimSummaryCurve;
//==================================================================================================
///
@ -39,7 +42,11 @@ protected:
virtual void setupActionLook( QAction* actionToSetup ) override;
private:
static RimGridSummaryCase* gridSummaryCaseForWell(RimEclipseWell* well);
static RimGridSummaryCase* gridSummaryCaseForWell(RimEclipseWell* well);
static bool isInjector(RimEclipseWell* well);
static RimSummaryCurve* addSummaryCurve(RimSummaryPlot* plot, const RimEclipseWell* well,
RimGridSummaryCase* gridSummaryCase, const QString& vectorName,
RimDefines::PlotAxis plotAxis, const cvf::Color3f& color);
};

12
ApplicationCode/Commands/FlowCommands/RicShowFlowCharacteristicsPlotFeature.cpp
View File

@ -23,6 +23,8 @@
#include "RimEclipseResultCase.h"
#include "RimEclipseView.h"
#include "RimFlowCharacteristicsPlot.h"
#include "RigFlowDiagResults.h"
#include "RimFlowDiagSolution.h"
#include "RimFlowPlotCollection.h"
#include "RimMainPlotCollection.h"
#include "RimProject.h"
@ -71,6 +73,16 @@ void RicShowFlowCharacteristicsPlotFeature::onActionTriggered(bool isChecked)
if (eclCase && eclCase->defaultFlowDiagSolution())
{
// Make sure flow results for the the active timestep is calculated, to avoid an empty plot
{
RimView * activeView = RiaApplication::instance()->activeReservoirView();
if (activeView && eclCase->defaultFlowDiagSolution()->flowDiagResults())
{
// Trigger calculation
eclCase->defaultFlowDiagSolution()->flowDiagResults()->maxAbsPairFlux(activeView->currentTimeStep());
}
}
if (RiaApplication::instance()->project())
{
RimFlowPlotCollection* flowPlotColl = RiaApplication::instance()->project()->mainPlotCollection->flowPlotCollection();

13
ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp
View File

@ -971,13 +971,22 @@ RigWellResultPoint RifReaderEclipseOutput::createWellResultPoint(const RigGridBa
resultPoint.m_oilRate = oilRate;
resultPoint.m_waterRate = waterRate;
/// Unit conversion for use with Well Allocation plots
// Convert Gas to oil equivalents
// If field unit, the Gas is in Mega ft^3 while the others are in [stb] (barrel)
// Unused Gas to Barrel conversion
// we convert gas to stb as well. Based on
// 1 [stb] = 0.15898729492800007 [m^3]
// 1 [ft] = 0.3048 [m]
// megaFt3ToStbFactor = 1.0 / (1.0e-6 * 0.15898729492800007 * ( 1.0 / 0.3048 )^3 )
double megaFt3ToStbFactor = 178107.60668;
if (m_eclipseCase->unitsType() == RigEclipseCaseData::UNITS_FIELD) gasRate = megaFt3ToStbFactor * gasRate;
// double megaFt3ToStbFactor = 178107.60668;
double fieldGasToOilEquivalent = 1.0e6/5800; // Mega ft^3 to BOE
double metricGasToOilEquivalent = 1.0/1.0e3; // Sm^3 Gas to Sm^3 oe
if (m_eclipseCase->unitsType() == RigEclipseCaseData::UNITS_FIELD) gasRate = fieldGasToOilEquivalent * gasRate;
if (m_eclipseCase->unitsType() == RigEclipseCaseData::UNITS_METRIC) gasRate = metricGasToOilEquivalent * gasRate;
resultPoint.m_gasRate = gasRate;
}

4
ApplicationCode/ProjectDataModel/Flow/RimFlowCharacteristicsPlot.cpp
View File

@ -141,14 +141,16 @@ QList<caf::PdmOptionItemInfo> RimFlowCharacteristicsPlot::calculateValueOptions(
{
std::vector<RimEclipseResultCase*> cases;
proj->descendantsIncludingThisOfType(cases);
RimEclipseResultCase* defaultCase = nullptr;
for ( RimEclipseResultCase* c : cases )
{
if ( c->defaultFlowDiagSolution() )
{
options.push_back(caf::PdmOptionItemInfo(c->caseUserDescription(), c, false, c->uiIcon()));
if (!defaultCase) defaultCase = c; // Select first
}
}
if (!m_case() && defaultCase) m_case = defaultCase;
}
}
else if ( fieldNeedingOptions == &m_flowDiagSolution )

47
ApplicationCode/ProjectDataModel/Flow/RimWellAllocationPlot.cpp
View File

@ -357,29 +357,46 @@ std::map<QString, const std::vector<double> *> RimWellAllocationPlot::findReleva
void RimWellAllocationPlot::updateWellFlowPlotXAxisTitle(RimWellLogTrack* plotTrack)
{
RigEclipseCaseData::UnitsType unitSet = m_case->eclipseCaseData()->unitsType();
QString unitText;
switch ( unitSet )
{
case RigEclipseCaseData::UNITS_METRIC:
unitText = "[m^3/day]";
break;
case RigEclipseCaseData::UNITS_FIELD:
unitText = "[Brl/day]";
break;
case RigEclipseCaseData::UNITS_LAB:
unitText = "[cm^3/hr]";
break;
default:
break;
}
if (m_flowDiagSolution)
{
QString unitText;
switch ( unitSet )
{
case RigEclipseCaseData::UNITS_METRIC:
unitText = "[m<sup>3</sup>/day]";
break;
case RigEclipseCaseData::UNITS_FIELD:
unitText = "[Brl/day]";
break;
case RigEclipseCaseData::UNITS_LAB:
unitText = "[cm<sup>3</sup>/hr]";
break;
default:
break;
}
plotTrack->setXAxisTitle("Reservoir Flow Rate " + unitText);
}
else
{
QString unitText;
switch ( unitSet )
{
case RigEclipseCaseData::UNITS_METRIC:
unitText = "[Liquid Sm<sup>3</sup>/day], [Gas kSm<sup>3</sup>/day]";
break;
case RigEclipseCaseData::UNITS_FIELD:
unitText = "[Liquid BBL/day], [Gas BOE/day]";
break;
case RigEclipseCaseData::UNITS_LAB:
unitText = "[cm<sup>3</sup>/hr]";
break;
default:
break;
}
plotTrack->setXAxisTitle("Surface Flow Rate " + unitText);
}
}

15
ApplicationCode/ProjectDataModel/Summary/RimSummaryCurveAppearanceCalculator.h
View File

@ -55,19 +55,18 @@ public:
void setupCurveLook(RimSummaryCurve* curve);
private:
static cvf::Color3f cycledPaletteColor(int colorIndex);
static cvf::Color3f cycledNoneRGBBrColor(int colorIndex);
static cvf::Color3f cycledGreenColor(int colorIndex);
static cvf::Color3f cycledBlueColor(int colorIndex);
static cvf::Color3f cycledRedColor(int colorIndex);
static cvf::Color3f cycledBrownColor(int colorIndex);
private:
void setOneCurveAppearance(CurveAppearanceType appeaType, size_t totalCount, int appeaIdx, RimSummaryCurve* curve);
void updateApperanceIndices();
std::map<std::string, size_t> mapNameToAppearanceIndex(CurveAppearanceType & appearance, const std::set<std::string>& names);
cvf::Color3f cycledPaletteColor(int colorIndex);
cvf::Color3f cycledNoneRGBBrColor(int colorIndex);
cvf::Color3f cycledGreenColor(int colorIndex);
cvf::Color3f cycledBlueColor(int colorIndex);
cvf::Color3f cycledRedColor(int colorIndex);
cvf::Color3f cycledBrownColor(int colorIndex);
RimPlotCurve::LineStyleEnum cycledLineStyle(int index);
RimPlotCurve::PointSymbolEnum cycledSymbol(int index);

1
ApplicationCode/ReservoirDataModel/RigAccWellFlowCalculator.h
View File

@ -74,7 +74,6 @@ public:
const std::vector<double>& pseudoLengthFromTop(size_t branchIdx) const;
const std::vector<double>& trueVerticalDepth(size_t branchIdx) const;
const std::vector<double>& accumulatedTracerFlowPrPseudoLength(const QString& tracerName, size_t branchIdx) const;
const std::vector<double>& flowPrPseudoLength( size_t branchIdx) const;
const std::vector<double>& tracerFlowPrPseudoLength(const QString& tracerName, size_t branchIdx) const;

21
ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h
View File

@ -21,7 +21,7 @@ namespace RigFlowDiagInterfaceTools {
template <class FluxCalc>
inline Opm::FlowDiagnostics::ConnectionValues
extractFluxField(const Opm::ECLGraph& G,
FluxCalc&& getFlux)
FluxCalc&& getFlux)
{
using ConnVals = Opm::FlowDiagnostics::ConnectionValues;
@ -52,24 +52,11 @@ namespace RigFlowDiagInterfaceTools {
}
inline Opm::FlowDiagnostics::ConnectionValues
extractFluxField(const Opm::ECLGraph& G,
const Opm::ECLRestartData& rstrt,
const bool compute_fluxes)
extractFluxFieldFromRestartFile(const Opm::ECLGraph& G,
const Opm::ECLRestartData& rstrt)
{
if (compute_fluxes) {
Opm::ECLFluxCalc calc(G);
auto getFlux = [&calc, &rstrt]
(const Opm::ECLGraph::PhaseIndex p)
{
return calc.flux(rstrt, p);
};
return extractFluxField(G, getFlux);
}
auto getFlux = [&G, &rstrt]
(const Opm::ECLGraph::PhaseIndex p)
(const Opm::ECLPhaseIndex p)
{
return G.flux(rstrt, p);
};

8
ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp
View File

@ -251,9 +251,8 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI
Opm::FlowDiagnostics::CellSetValues sumWellFluxPrCell;
{
Opm::FlowDiagnostics::ConnectionValues connectionsVals = RigFlowDiagInterfaceTools::extractFluxField(*(m_opmFlowDiagStaticData->m_eclGraph),
*currentRestartData,
false);
Opm::FlowDiagnostics::ConnectionValues connectionsVals = RigFlowDiagInterfaceTools::extractFluxFieldFromRestartFile(*(m_opmFlowDiagStaticData->m_eclGraph),
*currentRestartData);
m_opmFlowDiagStaticData->m_fldToolbox->assignConnectionFlux(connectionsVals);
@ -400,7 +399,8 @@ RigFlowDiagTimeStepResult RigFlowDiagSolverInterface::calculate(size_t timeStepI
{
Graph flowCapStorCapCurve = flowCapacityStorageCapacityCurve(*(injectorSolution.get()),
*(producerSolution.get()),
m_opmFlowDiagStaticData->m_poreVolume);
m_opmFlowDiagStaticData->m_poreVolume,
0.1);
result.setFlowCapStorageCapCurve(flowCapStorCapCurve);
result.setSweepEfficiencyCurve(sweepEfficiency(flowCapStorCapCurve));

10
ResInsightVersion.cmake
View File

@ -1,7 +1,7 @@
set(RESINSIGHT_MAJOR_VERSION 2016)
set(RESINSIGHT_MINOR_VERSION 11)
set(RESINSIGHT_INCREMENT_VERSION "flow.14")
set(RESINSIGHT_MAJOR_VERSION 2017)
set(RESINSIGHT_MINOR_VERSION 05)
set(RESINSIGHT_INCREMENT_VERSION "1-dev")
# https://github.com/CRAVA/crava/tree/master/libs/nrlib
@ -11,10 +11,10 @@ set(NRLIB_GITHUB_SHA "ba35d4359882f1c6f5e9dc30eb95fe52af50fd6f")
set(ERT_GITHUB_SHA "06a39878636af0bc52582430ad0431450e51139c")
# https://github.com/OPM/opm-flowdiagnostics
set(OPM_FLOWDIAGNOSTICS_SHA "2c5fb55db4c4ded49c14161dd16463e1207da049")
set(OPM_FLOWDIAGNOSTICS_SHA "b6e59ddcd2feba450c8612a7402c9239e442c0d4")
# https://github.com/OPM/opm-flowdiagnostics-applications
set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "570601718e7197b751bc3cba60c1e5fb7d842842")
set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "ccaaa4dd1b553e131a3051687fd615fe728b76ee")
# https://github.com/OPM/opm-parser/blob/master/opm/parser/eclipse/Units/Units.hpp
# This file was moved from opm-core to opm-parser october 2016

9
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake vendored
View File

@ -21,14 +21,19 @@
# the library needs it.
list (APPEND MAIN_SOURCE_FILES
opm/utility/ECLEndPointScaling.cpp
opm/utility/ECLFluxCalc.cpp
opm/utility/ECLGraph.cpp
opm/utility/ECLPropTable.cpp
opm/utility/ECLResultData.cpp
opm/utility/ECLSaturationFunc.cpp
opm/utility/ECLUnitHandling.cpp
opm/utility/ECLWellSolution.cpp
)
list (APPEND TEST_SOURCE_FILES
tests/test_eclendpointscaling.cpp
tests/test_eclproptable.cpp
tests/test_eclunithandling.cpp
)
@ -44,9 +49,13 @@ list (APPEND EXAMPLE_SOURCE_FILES
)
list (APPEND PUBLIC_HEADER_FILES
opm/utility/ECLEndPointScaling.hpp
opm/utility/ECLFluxCalc.hpp
opm/utility/ECLGraph.hpp
opm/utility/ECLPhaseIndex.hpp
opm/utility/ECLPropTable.hpp
opm/utility/ECLResultData.hpp
opm/utility/ECLSaturationFunc.hpp
opm/utility/ECLUnitHandling.hpp
opm/utility/ECLWellSolution.hpp
)

37
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp vendored
View File

@ -38,34 +38,27 @@ try {
std::vector<Opm::FlowDiagnostics::CellSet> start;
auto fwd = fdTool.computeInjectionDiagnostics(start);
auto rev = fdTool.computeProductionDiagnostics(start);
// Give disconnected cells zero pore volume.
std::vector<int> nbcount(setup.graph.numCells(), 0);
for (int nb : setup.graph.neighbours()) {
if (nb >= 0) {
++nbcount[nb];
}
}
auto pv = setup.graph.poreVolume();
for (size_t i = 0; i < pv.size(); ++i) {
if (nbcount[i] == 0) {
pv[i] = 0.0;
}
}
// Cells that have more than 100 times the average pore volume are
// probably aquifers, we ignore them (again by setting pore volume
// to zero). If this is the correct thing to do or not could
// depend on what you want to use the diagnostic for.
const double average_pv = std::accumulate(pv.begin(), pv.end(), 0.0) / double(pv.size());
for (double& pvval : pv) {
if (pvval > 100.0 * average_pv) {
pvval = 0.0;
const bool ignore_disconnected = setup.param.getDefault("ignore_disconnected", true);
if (ignore_disconnected) {
// Give disconnected cells zero pore volume.
std::vector<int> nbcount(setup.graph.numCells(), 0);
for (int nb : setup.graph.neighbours()) {
if (nb >= 0) {
++nbcount[nb];
}
}
for (size_t i = 0; i < pv.size(); ++i) {
if (nbcount[i] == 0) {
pv[i] = 0.0;
}
}
}
// Compute graph.
auto fphi = Opm::FlowDiagnostics::flowCapacityStorageCapacityCurve(fwd, rev, pv);
const double max_pv_fraction = setup.param.getDefault("max_pv_fraction", 0.1);
auto fphi = Opm::FlowDiagnostics::flowCapacityStorageCapacityCurve(fwd, rev, pv, max_pv_fraction);
// Write it to standard out.
std::cout.precision(16);

38
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp vendored
View File

@ -31,6 +31,7 @@
#include <opm/utility/ECLFluxCalc.hpp>
#include <opm/utility/ECLGraph.hpp>
#include <opm/utility/ECLPhaseIndex.hpp>
#include <opm/utility/ECLResultData.hpp>
#include <opm/utility/ECLWellSolution.hpp>
@ -115,15 +116,19 @@ namespace example {
}
inline Opm::FlowDiagnostics::ConnectionValues
extractFluxField(const Opm::ECLGraph& G,
const Opm::ECLRestartData& rstrt,
const bool compute_fluxes)
extractFluxField(const Opm::ECLGraph& G,
const Opm::ECLInitFileData& init,
const Opm::ECLRestartData& rstrt,
const bool compute_fluxes,
const bool useEPS)
{
if (compute_fluxes) {
Opm::ECLFluxCalc calc(G);
auto satfunc = ::Opm::ECLSaturationFunc(G, init, useEPS);
Opm::ECLFluxCalc calc(G, std::move(satfunc));
auto getFlux = [&calc, &rstrt]
(const Opm::ECLGraph::PhaseIndex p)
(const Opm::ECLPhaseIndex p)
{
return calc.flux(rstrt, p);
};
@ -132,7 +137,7 @@ namespace example {
}
auto getFlux = [&G, &rstrt]
(const Opm::ECLGraph::PhaseIndex p)
(const Opm::ECLPhaseIndex p)
{
return G.flux(rstrt, p);
};
@ -205,17 +210,6 @@ namespace example {
inline Opm::ECLGraph
initGraph(const FilePaths& file_paths)
{
const auto I = Opm::ECLInitFileData(file_paths.init);
return Opm::ECLGraph::load(file_paths.grid, I);
}
inline Opm::FlowDiagnostics::Toolbox
initToolbox(const Opm::ECLGraph& G)
{
@ -238,11 +232,13 @@ namespace example {
Setup(int argc, char** argv)
: param (initParam(argc, argv))
, file_paths (param)
, init (file_paths.init)
, rstrt (file_paths.restart)
, graph (initGraph(file_paths))
, graph (::Opm::ECLGraph::load(file_paths.grid, init))
, well_fluxes ()
, toolbox (initToolbox(graph))
, compute_fluxes_(param.getDefault("compute_fluxes", false))
, useEPS_ (param.getDefault("use_ep_scaling", false))
{
const int step = param.getDefault("step", 0);
@ -268,7 +264,9 @@ namespace example {
well_fluxes = wsol.solution(rstrt, graph.activeGrids());
}
toolbox.assignConnectionFlux(extractFluxField(graph, rstrt, compute_fluxes_));
toolbox.assignConnectionFlux(extractFluxField(graph, init, rstrt,
compute_fluxes_, useEPS_));
toolbox.assignInflowFlux(extractWellFlows(graph, well_fluxes));
return true;
@ -276,11 +274,13 @@ namespace example {
Opm::ParameterGroup param;
FilePaths file_paths;
Opm::ECLInitFileData init;
Opm::ECLRestartData rstrt;
Opm::ECLGraph graph;
std::vector<Opm::ECLWellSolution::WellData> well_fluxes;
Opm::FlowDiagnostics::Toolbox toolbox;
bool compute_fluxes_ = false;
bool useEPS_ = false;
};

1163
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp vendored Normal file
View File

File diff suppressed because it is too large Load Diff

431
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp vendored Normal file
View File

@ -0,0 +1,431 @@
/*
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLENDPOINTSCALING_HEADER_INCLUDED
#define OPM_ECLENDPOINTSCALING_HEADER_INCLUDED
#include <opm/utility/ECLPhaseIndex.hpp>
#include <functional>
#include <memory>
#include <string>
#include <vector>
namespace Opm {
class ECLGraph;
class ECLInitFileData;
} // namespace Opm
namespace Opm { namespace SatFunc {
/// Protocol for computing scaled saturation values.
class EPSEvalInterface
{
public:
/// Static (unscaled) end-points of a particular, tabulated
/// saturation function.
struct TableEndPoints {
/// Critical or connate/minimum saturation, depending on
/// subsequent function evaluation context.
///
/// Use critical saturation when computing look-up values for
/// relative permeability and connate/minimum saturation when
/// computing look-up values for capillary pressure.
double low;
/// Displacing saturation (3-pt option only).
double disp;
/// Maximum (high) saturation point.
double high;
};
/// Associate a saturation value to a specific cell.
struct SaturationAssoc {
/// Cell to which to connect a saturation value.
std::vector<int>::size_type cell;
/// Saturation value.
double sat;
};
/// Convenience type alias.
using SaturationPoints = std::vector<SaturationAssoc>;
/// Derive scaled saturations--inputs to subsequent evaluation of
/// saturation functions--corresponding to a sequence of unscaled
/// saturation values.
///
/// \param[in] tep Static end points that identify the saturation
/// scaling intervals of a particular tabulated saturation
/// function.
///
/// \param[in] sp Sequence of saturation points.
///
/// \return Sequence of scaled saturation values in order of the
/// input sequence. In particular the \c i-th element of this
/// result is the scaled version of \code sp[i].sat \endcode.
virtual std::vector<double>
eval(const TableEndPoints& tep,
const SaturationPoints& sp) const = 0;
virtual std::unique_ptr<EPSEvalInterface> clone() const = 0;
/// Destructor. Must be virtual.
virtual ~EPSEvalInterface();
};
/// Implementation of ECLIPSE's standard, two-point, saturation scaling
/// option.
class TwoPointScaling : public EPSEvalInterface
{
public:
/// Constructor.
///
/// Typically set up to define the end-point scaling of all active
/// cells in a model, but could alternatively be used as a means to
/// computing the effective saturation function of a single cell.
///
/// \param[in] smin Left end points for a set of cells.
///
/// \param[in] smax Right end points for a set of cells.
TwoPointScaling(std::vector<double> smin,
std::vector<double> smax);
/// Destructor.
~TwoPointScaling();
/// Copy constructor.
///
/// \param[in] rhs Existing object.
TwoPointScaling(const TwoPointScaling& rhs);
/// Move constructor.
///
/// Subsumes the implementation of an existing object.
///
/// \param[in] rhs Existing object.
TwoPointScaling(TwoPointScaling&& rhs);
/// Assignment operator.
///
/// Replaces current implementation with a copy of existing object's
/// implementation details.
///
/// \param[in] rhs Existing object.
///
/// \return \code *this \endcode.
TwoPointScaling&
operator=(const TwoPointScaling& rhs);
/// Move assingment operator.
///
/// Subsumes existing object's implementation details and uses those
/// to replace the current implementation.
///
/// \return \code *this \endcode.
TwoPointScaling&
operator=(TwoPointScaling&& rhs);
/// Derive scaled saturations using the two-point scaling definition
/// from a sequence of unscaled saturation values.
///
/// \param[in] tep Static end points that identify the saturation
/// scaling intervals of a particular tabulated saturation
/// function. The evaluation procedure considers only \code
/// tep.low \endcode and \code tep.high \endcode. The value of
/// \code tep.disp \endcode is never read.
///
/// \param[in] sp Sequence of saturation points. The maximum cell
/// index (\code sp[i].cell \endcode) must be strictly less than
/// the size of the input arrays that define the current
/// saturation regions.
///
/// \return Sequence of scaled saturation values in order of the
/// input sequence. In particular the \c i-th element of this
/// result is the scaled version of \code sp[i].sat \endcode.
virtual std::vector<double>
eval(const TableEndPoints& tep,
const SaturationPoints& sp) const override;
virtual std::unique_ptr<EPSEvalInterface> clone() const override;
private:
/// Implementation class.
class Impl;
/// Pimpl idiom.
std::unique_ptr<Impl> pImpl_;
};
/// Implementation of ECLIPSE's alternative, three-point, saturation
/// scaling option.
class ThreePointScaling : public EPSEvalInterface
{
public:
/// Constructor.
///
/// Typically set up to define the end-point scaling of all active
/// cells in a model, but could alternatively be used as a means to
/// computing the effective saturation function of a single cell.
///
/// \param[in] smin Left end points for a set of cells.
///
/// \param[in] sdisp Intermediate--displacing saturation--end points
/// for a set of cells.
///
/// \param[in] smax Right end points for a set of cells.
ThreePointScaling(std::vector<double> smin,
std::vector<double> sdisp,
std::vector<double> smax);
/// Destructor.
~ThreePointScaling();
/// Copy constructor.
///
/// \param[in] rhs Existing object.
ThreePointScaling(const ThreePointScaling& rhs);
/// Move constructor.
///
/// Subsumes the implementation of an existing object.
///
/// \param[in] rhs Existing object.
ThreePointScaling(ThreePointScaling&& rhs);
/// Assignment operator.
///
/// Replaces current implementation with a copy of existing object's
/// implementation details.
///
/// \param[in] rhs Existing object.
///
/// \return \code *this \endcode.
ThreePointScaling&
operator=(const ThreePointScaling& rhs);
/// Move assingment operator.
///
/// Subsumes existing object's implementation details and uses those
/// to replace the current implementation.
///
/// \return \code *this \endcode.
ThreePointScaling&
operator=(ThreePointScaling&& rhs);
/// Derive scaled saturations using the three-point scaling
/// definition from a sequence of unscaled saturation values.
///
/// \param[in] tep Static end points that identify the saturation
/// scaling intervals of a particular tabulated saturation
/// function. The evaluation procedure considers only \code
/// tep.low \endcode and \code tep.high \endcode. The value of
/// \code tep.disp \endcode is never read.
///
/// \param[in] sp Sequence of saturation points. The maximum cell
/// index (\code sp[i].cell \endcode) must be strictly less than
/// the size of the input arrays that define the current
/// saturation regions.
///
/// \return Sequence of scaled saturation values in order of the
/// input sequence. In particular the \c i-th element of this
/// result is the scaled version of \code sp[i].sat \endcode.
virtual std::vector<double>
eval(const TableEndPoints& tep,
const SaturationPoints& sp) const override;
virtual std::unique_ptr<EPSEvalInterface> clone() const override;
private:
/// Implementation class.
class Impl;
/// Pimpl idiom.
std::unique_ptr<Impl> pImpl_;
};
/// Named constructors for enabling saturation end-point scaling from an
/// ECL result set (see class \c ECLInitFileData).
struct CreateEPS
{
/// Category of function for which to create an EPS evaluator.
enum class FunctionCategory {
/// This EPS is for relative permeability. Possibly uses
/// three-point (alternative) formulation.
Relperm,
/// This EPS is for capillary pressure. Two-point option only.
CapPress,
};
/// Categories of saturation function systems for which to create an
/// EPS evaluator.
enum class SubSystem {
/// Create an EPS for a curve in the Oil-Water (sub-) system.
OilWater,
/// Create an EPS for a curve in the Oil-Gas (sub-) system.
OilGas,
};
/// Set of options that uniquely define a single EPS operation.
struct EPSOptions {
/// Whether or not to employ the alternative (i.e., 3-pt) scaling
/// procedure. Only applicable to FunctionCategory::Relperm and
/// ignored in the case of FunctionCategory::CapPress.
bool use3PtScaling;
/// Curve-type for which to create an EPS.
FunctionCategory curve;
/// Part of global fluid system for which to create an EPS.
SubSystem subSys;
/// Phase for whose \c curve in which \c subSys to create an
/// EPS.
///
/// Example: Create a standard (two-point) EPS for the relative
/// permeability of oil in the oil-gas subsystem of an
/// oil-gas-water active phase system.
///
/// \code
/// auto opt = EPSOptions{};
///
/// opt.use3PtScaling = false;
/// opt.curve = FunctionCategory::Relperm;
/// opt.subSys = SubSystem::OilGas;
/// opt.thisPh = ECLPhaseIndex::Oil;
///
/// auto eps = CreateEPS::fromECLOutput(G, init, opt);
/// \endcode
::Opm::ECLPhaseIndex thisPh;
};
/// Collection of raw saturation table end points.
struct RawTableEndPoints {
/// Collection of connate (minimum) saturation end points.
struct Connate {
/// Connate oil saturation for each table in total set of
/// tabulated saturation functions.
std::vector<double> oil;
/// Connate gas saturation for each table in total set of
/// tabulated saturation functions.
std::vector<double> gas;
/// Connate water saturation for each table in total set of
/// tabulated saturation functions.
std::vector<double> water;
};
/// Collection of critical saturations. Used in deriving scaled
/// displacing saturations in the alternative (three-point)
/// scaling procedure.
struct Critical {
/// Critical oil saturation in 2p OG system from total set
/// of tabulated saturation functions.
std::vector<double> oil_in_gas;
/// Critical oil saturation in 2p OW system from total set
/// of tabulated saturation functions.
std::vector<double> oil_in_water;
/// Critical gas saturation in 2p OG or 3p OGW system from
/// total set of tabulated saturation functions.
std::vector<double> gas;
/// Critical water saturation in 2p OW or 3p OGW system from
/// total set of tabulated saturation functions.
std::vector<double> water;
};
/// Collection of maximum saturation end points.
struct Maximum {
/// Maximum oil saturation for each table in total set of
/// tabulated saturation functions.
std::vector<double> oil;
/// Maximum gas saturation for each table in total set of
/// tabulated saturation functions.
std::vector<double> gas;
/// Maximum water saturation for each table in total set of
/// tabulated saturation functions.
std::vector<double> water;
};
/// Minimum saturation end points for all tabulated saturation
/// functions.
Connate conn;
/// Critical saturations for all tabulated saturation functions.
Critical crit;
/// Maximum saturation end points for all tabulated saturation
/// functions.
Maximum smax;
};
/// Construct an EPS evaluator from a particular ECL result set.
///
/// \param[in] G Connected topology of current model's active cells.
/// Needed to linearise table end-points that may be distributed
/// on local grids to all of the model's active cells (\code
/// member function G.rawLinearisedCellData() \endcode).
///
/// \param[in] init Container of tabulated saturation functions and
/// saturation table end points for all active cells.
///
/// \param[in] opt Options that identify a particular end-point
/// scaling behaviour of a particular saturation function curve.
///
/// \return EPS evaluator for the particular curve defined by the
/// input options.
static std::unique_ptr<EPSEvalInterface>
fromECLOutput(const ECLGraph& G,
const ECLInitFileData& init,
const EPSOptions& opt);
/// Extract table end points relevant to a particular EPS evaluator
/// from raw tabulated saturation functions.
///
/// \param[in] ep Collection of all raw table saturation end points
/// for all tabulated saturation functions. Typically computed
/// by direct calls to the \code connateSat() \endcode, \code
/// criticalSat() \endcode, and \code maximumSat() \endcode of
/// the currently active \code Opm::SatFuncInterpolant \code
/// objects.
///
/// \param[in] opt Options that identify a particular end-point
/// scaling behaviour of a particular saturation function curve.
///
/// \return Subset of the input end points in a format intended for
/// passing as the first argument of member function \code eval()
/// \endcode of the \code EPSEvalInterface \endcode that
/// corresponds to the input options.
static std::vector<EPSEvalInterface::TableEndPoints>
unscaledEndPoints(const RawTableEndPoints& ep,
const EPSOptions& opt);
};
}} // namespace Opm::SatFunc
#endif // OPM_ECLENDPOINTSCALING_HEADER_INCLUDED

17
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp vendored
View File

@ -22,11 +22,15 @@
#include <opm/parser/eclipse/Units/Units.hpp>
#include <utility>
namespace Opm
{
ECLFluxCalc::ECLFluxCalc(const ECLGraph& graph)
ECLFluxCalc::ECLFluxCalc(const ECLGraph& graph,
ECLSaturationFunc&& satfunc)
: graph_(graph)
, satfunc_(std::move(satfunc))
, neighbours_(graph.neighbours())
, transmissibility_(graph.transmissibility())
{
@ -38,7 +42,7 @@ namespace Opm
std::vector<double>
ECLFluxCalc::flux(const ECLRestartData& rstrt,
const PhaseIndex /* phase */) const
const ECLPhaseIndex phase) const
{
// Obtain dynamic data.
DynamicData dyn_data;
@ -46,6 +50,9 @@ namespace Opm
.linearisedCellData(rstrt, "PRESSURE",
&ECLUnits::UnitSystem::pressure);
dyn_data.relperm = this->satfunc_
.relperm(this->graph_, rstrt, phase);
// Compute fluxes per connection.
const int num_conn = transmissibility_.size();
std::vector<double> fluxvec(num_conn);
@ -66,8 +73,12 @@ namespace Opm
const int c2 = neighbours_[2*connection + 1];
const double transmissibility = transmissibility_[connection];
const double viscosity = 1.0 * prefix::centi * unit::Poise;
const double mobility = 1.0 / viscosity;
const auto& pressure = dyn_data.pressure;
const int upwind_cell = (pressure[c2] > pressure[c1]) ? c2 : c1;
const double kr = dyn_data.relperm[upwind_cell];
const double mobility = kr / viscosity;
return mobility * transmissibility * (pressure[c1] - pressure[c2]);
}

15
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp vendored
View File

@ -21,6 +21,9 @@
#define OPM_ECLFLUXCALC_HEADER_INCLUDED
#include <opm/utility/ECLGraph.hpp>
#include <opm/utility/ECLPhaseIndex.hpp>
#include <opm/utility/ECLSaturationFunc.hpp>
#include <vector>
namespace Opm
@ -34,13 +37,15 @@ namespace Opm
/// Construct from ECLGraph.
///
/// \param[in] graph Connectivity data, as well as providing a means to read data from the restart file.
explicit ECLFluxCalc(const ECLGraph& graph);
using PhaseIndex = ECLGraph::PhaseIndex;
explicit ECLFluxCalc(const ECLGraph& graph,
ECLSaturationFunc&& satfunc);
/// Retrive phase flux on all connections defined by \code
/// graph.neighbours() \endcode.
///
/// \param[in] rstrt ECL Restart data set from which to extract
/// relevant data per cell.
///
/// \param[in] phase Canonical phase for which to retrive flux.
///
/// \return Flux values corresponding to selected phase.
@ -48,18 +53,20 @@ namespace Opm
/// Numerical values in SI units (rm^3/s).
std::vector<double>
flux(const ECLRestartData& rstrt,
const PhaseIndex phase) const;
const ECLPhaseIndex phase) const;
private:
struct DynamicData
{
std::vector<double> pressure;
std::vector<double> relperm;
};
double singleFlux(const int connection,
const DynamicData& dyn_data) const;
const ECLGraph& graph_;
ECLSaturationFunc satfunc_;
std::vector<int> neighbours_;
std::vector<double> transmissibility_;
};

32
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp vendored
View File

@ -1278,7 +1278,7 @@ public:
///
/// Mostly useful to determine the set of \c PhaseIndex values for which
/// flux() may return non-zero values.
const std::vector<PhaseIndex>& activePhases() const;
const std::vector<ECLPhaseIndex>& activePhases() const;
/// Retrieve the simulation scenario's set of active grids.
///
@ -1322,7 +1322,7 @@ public:
/// all).
std::vector<double>
flux(const ECLRestartData& rstrt,
const PhaseIndex phase) const;
const ECLPhaseIndex phase) const;
/// Retrieve result set vector from current view linearised on active
/// cells.
@ -1609,7 +1609,7 @@ private:
/// Set of active phases in result set. Derived from .INIT on the
/// assumption that the set of active phases does not change throughout
/// the simulation run.
std::vector<PhaseIndex> activePhases_;
std::vector<ECLPhaseIndex> activePhases_;
/// Set of active grids in result set.
std::vector<std::string> activeGrids_;
@ -1639,7 +1639,7 @@ private:
/// \return Basename for ECl vector corresponding to particular phase
/// flux.
std::string
flowVector(const PhaseIndex phase) const;
flowVector(const ECLPhaseIndex phase) const;
/// Extract flux values corresponding to particular result set vector
/// for all identified non-neighbouring connections.
@ -1940,7 +1940,7 @@ Opm::ECLGraph::Impl::numConnections() const
return nconn + this->nnc_.numConnections();
}
const std::vector<Opm::ECLGraph::PhaseIndex>&
const std::vector<Opm::ECLPhaseIndex>&
Opm::ECLGraph::Impl::activePhases() const
{
return this->activePhases_;
@ -2028,7 +2028,7 @@ Opm::ECLGraph::Impl::transmissibility() const
std::vector<double>
Opm::ECLGraph::Impl::flux(const ECLRestartData& rstrt,
const PhaseIndex phase) const
const ECLPhaseIndex phase) const
{
auto fluxUnit = [&rstrt](const std::string& gridID)
{
@ -2226,12 +2226,12 @@ defineActivePhases(const ::Opm::ECLInitFileData& init)
const auto phaseMask =
static_cast<unsigned int>(ih[INTEHEAD_PHASE_INDEX]);
using Check = std::pair<PhaseIndex, unsigned int>;
using Check = std::pair<ECLPhaseIndex, unsigned int>;
const auto allPhases = std::vector<Check> {
{ PhaseIndex::Aqua , (1u << 1) },
{ PhaseIndex::Liquid, (1u << 0) },
{ PhaseIndex::Vapour, (1u << 2) },
{ ECLPhaseIndex::Aqua , (1u << 1) },
{ ECLPhaseIndex::Liquid, (1u << 0) },
{ ECLPhaseIndex::Vapour, (1u << 2) },
};
this->activePhases_.clear();
@ -2243,19 +2243,19 @@ defineActivePhases(const ::Opm::ECLInitFileData& init)
}
std::string
Opm::ECLGraph::Impl::flowVector(const PhaseIndex phase) const
Opm::ECLGraph::Impl::flowVector(const ECLPhaseIndex phase) const
{
const auto vector = std::string("FLR"); // Flow-rate, reservoir
if (phase == PhaseIndex::Aqua) {
if (phase == ECLPhaseIndex::Aqua) {
return vector + "WAT";
}
if (phase == PhaseIndex::Liquid) {
if (phase == ECLPhaseIndex::Liquid) {
return vector + "OIL";
}
if (phase == PhaseIndex::Vapour) {
if (phase == ECLPhaseIndex::Vapour) {
return vector + "GAS";
}
@ -2322,7 +2322,7 @@ std::size_t Opm::ECLGraph::numConnections() const
return this->pImpl_->numConnections();
}
const std::vector<Opm::ECLGraph::PhaseIndex>&
const std::vector<Opm::ECLPhaseIndex >&
Opm::ECLGraph::activePhases() const
{
return this->pImpl_->activePhases();
@ -2351,7 +2351,7 @@ std::vector<double> Opm::ECLGraph::transmissibility() const
std::vector<double>
Opm::ECLGraph::flux(const ECLRestartData& rstrt,
const PhaseIndex phase) const
const ECLPhaseIndex phase) const
{
return this->pImpl_->flux(rstrt, phase);
}

8
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp vendored
View File

@ -20,6 +20,7 @@
#ifndef OPM_ECLGRAPH_HEADER_INCLUDED
#define OPM_ECLGRAPH_HEADER_INCLUDED
#include <opm/utility/ECLPhaseIndex.hpp>
#include <opm/utility/ECLResultData.hpp>
#include <opm/utility/ECLUnitHandling.hpp>
@ -45,9 +46,6 @@ namespace Opm {
class ECLGraph
{
public:
/// Enum for indicating requested phase from the flux() method.
enum class PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
/// Disabled default constructor.
ECLGraph() = delete;
@ -125,7 +123,7 @@ namespace Opm {
///
/// Mostly useful to determine the set of \c PhaseIndex values for
/// which flux() will return non-zero values if data available.
const std::vector<PhaseIndex>& activePhases() const;
const std::vector<ECLPhaseIndex>& activePhases() const;
/// Retrieve the simulation scenario's set of active grids.
///
@ -167,7 +165,7 @@ namespace Opm {
/// (rm^3/s).
std::vector<double>
flux(const ECLRestartData& rstrt,
const PhaseIndex phase) const;
const ECLPhaseIndex phase) const;
/// Retrieve result set vector from current view (e.g., particular
/// report step) linearised on active cells.

32
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPhaseIndex.hpp vendored Normal file
View File

@ -0,0 +1,32 @@
/*
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLPHASEINDEX_HEADER_INCLUDED
#define OPM_ECLPHASEINDEX_HEADER_INCLUDED
namespace Opm {
/// Enum for indicating the phase--or set of phases--on which to apply a
/// phase-dependent operation (e.g., extracting flux data from a result
/// set or computing relative permeabilities from tabulated functions).
enum class ECLPhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
} // namespace Opm
#endif // OPM_ECLPHASEINDEX_HEADER_INCLUDED

305
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.cpp vendored Normal file
View File

@ -0,0 +1,305 @@
/*
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/utility/ECLPropTable.hpp>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <exception>
#include <iterator>
#include <stdexcept>
#include <utility>
Opm::SatFuncInterpolant::SingleTable::
SingleTable(ElmIt xBegin,
ElmIt xEnd,
std::vector<ElmIt>& colIt)
{
// There must be at least one dependent variable/result variable.
assert (colIt.size() >= 1);
const auto nRows = std::distance(xBegin, xEnd);
this->x_.reserve(nRows);
this->y_.reserve(nRows * colIt.size());
auto keyValid = [](const double xi)
{
// Indep. variable values <= -1.0e20 or >= 1.0e20 signal "unused"
// table nodes (rows). These nodes are in the table to fill out the
// allocated size if one particular sub-table does not use all
// nodes. The magic value 1.0e20 is documented in the Fileformats
// Reference Manual.
return std::abs(xi) < 1.0e20;
};
while (xBegin != xEnd) {
// Extract relevant portion of the table. Preallocated rows that
// are not actually part of the result set (i.e., those that are set
// to a sentinel value) are discarded.
if (keyValid(*xBegin)) {
this->x_.push_back(*xBegin);
for (auto ci : colIt) {
// Store 'y_' with column index cycling most rapidly.
this->y_.push_back(*ci);
}
}
// -------------------------------------------------------------
// Advance iterators.
// 1) Independent variable.
++xBegin;
// 2) Dependent/result/columns.
for (auto& ci : colIt) {
++ci;
}
}
// Dispose of any excess capacity.
if (this->x_.size() < static_cast<decltype(this->x_.size())>(nRows)) {
this->x_.shrink_to_fit();
this->y_.shrink_to_fit();
}
if (this->x_.size() < 2) {
// Table has no interval that supports interpolation. Either just a
// single node or no nodes at all. We can't do anything useful
// here, so don't pretend that this is okay.
throw std::invalid_argument {
"No Interpolation Intervals of Non-Zero Size"
};
}
}
double
Opm::SatFuncInterpolant::SingleTable::
y(const ECLPropTableRawData::SizeType nCols,
const ECLPropTableRawData::SizeType row,
const ResultColumn& c) const
{
assert (row * nCols < this->y_.size());
assert (c.i < nCols);
// Recall: 'y_' stored with column index cycling the most rapidly (row
// major ordering).
return this->y_[row*nCols + c.i];
}
std::vector<double>
Opm::SatFuncInterpolant::SingleTable::
interpolate(const ECLPropTableRawData::SizeType nCols,
const ResultColumn& c,
const std::vector<double>& x) const
{
auto y = std::vector<double>{}; y.reserve(x.size());
auto yval = [nCols, c, this]
(const ECLPropTableRawData::SizeType i)
{
return this->y(nCols, i, c);
};
const auto yfirst =
yval(ECLPropTableRawData::SizeType{ 0 });
const auto ylast =
yval(ECLPropTableRawData::SizeType{ this->x_.size() - 1 });
for (const auto& xi : x) {
y.push_back(0.0);
auto& yi = y.back();
if (! (xi > this->x_.front())) {
// Constant extrapolation to the left of range.
yi = yfirst;
}
else if (! (xi < this->x_.back())) {
// Constant extrapolation to the right of range.
yi = ylast;
}
else {
// Somewhere in [min(x_), max(x_)]. Primary key (indep. var) is
// sorted range. Recall: lower_bound() returns insertion point,
// which translates to the *upper* (right-hand) end-point of the
// interval in this context.
auto b = std::begin(this->x_);
auto p = std::lower_bound(b, std::end(this->x_), xi);
assert ((p != b) && "Logic Error Left End-Point");
assert ((p != std::end(this->x_)) &&
"Logic Error Right End-Point");
// p = lower_bound() => left == i-1, right == i-0.
const auto i = p - b;
const auto left = i - 1;
const auto right = i - 0;
const auto xl = this->x_[left];
const auto t = (xi - xl) / (this->x_[right] - xl);
yi = (1.0 - t)*yval(left) + t*yval(right);
}
}
return y;
}
double
Opm::SatFuncInterpolant::SingleTable::connateSat() const
{
return this->x_.front();
}
double
Opm::SatFuncInterpolant::SingleTable::
criticalSat(const ECLPropTableRawData::SizeType nCols,
const ResultColumn& c) const
{
// Note: Relative permeability functions are presented as non-decreasing
// functions of the corresponding phase saturation. The internal table
// format essentially mirrors that of input deck keywords SWFN, SGFN,
// and SOF* (i.e., saturation function family II). Extracting the
// critical saturation--even for oil--therefore amounts to a forward,
// linear scan from row=0 to row=n-1 irrespective of the input format of
// the current saturation function.
const auto nRows = this->x_.size();
auto row = 0 * nRows;
for (; row < nRows; ++row) {
if (this->y(nCols, row, c) > 0.0) { break; }
}
if (row == 0) {
throw std::invalid_argument {
"Table Does Not Define Critical Saturation"
};
}
return this->x_[row - 1];
}
double
Opm::SatFuncInterpolant::SingleTable::maximumSat() const
{
return this->x_.back();
}
// =====================================================================
Opm::SatFuncInterpolant::SatFuncInterpolant(const ECLPropTableRawData& raw)
: nResCols_(raw.numCols - 1)
{
if (raw.numCols < 2) {
throw std::invalid_argument {
"Malformed Property Table"
};
}
this->table_.reserve(raw.numTables);
// Table format: numRows*numTables values of first column (indep. var)
// followed by numCols-1 dependent variable (function value result)
// columns of numRows*numTables values each, one column at a time.
const auto colStride = raw.numRows * raw.numTables;
// Position column iterators (independent variable and results
// respectively) at beginning of each pertinent table column.
auto xBegin = std::begin(raw.data);
auto colIt = std::vector<decltype(xBegin)>{ xBegin + colStride };
for (auto col = 0*raw.numCols + 1; col < raw.numCols - 1; ++col) {
colIt.push_back(colIt.back() + colStride);
}
for (auto t = 0*raw.numTables;
t < raw.numTables;
++t, xBegin += raw.numRows)
{
auto xEnd = xBegin + raw.numRows;
// Note: The SingleTable ctor advances each 'colIt' across numRows
// entries. That is a bit of a layering violation, but helps in the
// implementation of this loop.
this->table_.push_back(SingleTable(xBegin, xEnd, colIt));
}
}
std::vector<double>
Opm::SatFuncInterpolant::interpolate(const InTable& t,
const ResultColumn& c,
const std::vector<double>& x) const
{
if (t.i >= this->table_.size()) {
throw std::invalid_argument {
"Invalid Table ID"
};
}
if (c.i >= this->nResCols_) {
throw std::invalid_argument {
"Invalid Result Column ID"
};
}
return this->table_[t.i].interpolate(this->nResCols_, c, x);
}
std::vector<double>
Opm::SatFuncInterpolant::connateSat() const
{
auto sconn = std::vector<double>{};
sconn.reserve(this->table_.size());
for (const auto& t : this->table_) {
sconn.push_back(t.connateSat());
}
return sconn;
}
std::vector<double>
Opm::SatFuncInterpolant::criticalSat(const ResultColumn& c) const
{
auto scrit = std::vector<double>{};
scrit.reserve(this->table_.size());
for (const auto& t : this->table_) {
scrit.push_back(t.criticalSat(this->nResCols_, c));
}
return scrit;
}
std::vector<double>
Opm::SatFuncInterpolant::maximumSat() const
{
auto smax = std::vector<double>{};
smax.reserve(this->table_.size());
for (const auto& t : this->table_) {
smax.push_back(t.maximumSat());
}
return smax;
}

182
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.hpp vendored Normal file
View File

@ -0,0 +1,182 @@
/*
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLPROPTABLE_HEADER_INCLUDED
#define OPM_ECLPROPTABLE_HEADER_INCLUDED
#include <vector>
/// \file
///
/// ECL Tabulated Functions (e.g., saturation functions).
namespace Opm {
/// Raw table data from which to construct collection of interpolants.
struct ECLPropTableRawData
{
/// Representation of the raw table data. 1D array with implicit
/// substructure.
using DataVector = std::vector<double>;
/// Size type for subsets of table data.
using SizeType = DataVector::size_type;
/// Iterator to table elements. Must be random access.
using ElementIterator = DataVector::const_iterator;
/// Raw table data. Column major (Fortran) order. Typically
/// copied/extracted directly from TAB vector of INIT result-set.
DataVector data;
/// Number of rows allocated in the result set for each individual
/// table. Typically corresponds to setting in one of the *DIMS
/// keywords. Should normally be at least two.
SizeType numRows;
/// Number of columns in this table. Varies by keyword/table.
SizeType numCols;
/// Number of tables of this type. Must match the corresponding
/// region keyword.
SizeType numTables;
};
/// Collection of 1D interpolants from tabulated functions (e.g., the
/// saturation functions).
class SatFuncInterpolant
{
public:
/// Constructor.
///
/// \param[in] raw Raw table data for this collection.
explicit SatFuncInterpolant(const ECLPropTableRawData& raw);
/// Wrapper type to disambiguate API usage. Represents a table ID.
struct InTable {
/// Table ID.
ECLPropTableRawData::SizeType i;
};
/// Wrapper type to disambiguate API usage. Represents a column ID.
struct ResultColumn {
/// Column ID.
ECLPropTableRawData::SizeType i;
};
/// Evaluate 1D interpolant in sequence of points.
///
/// \param[in] t ID of sub-table of interpolant.
///
/// \param[in] c ID of result column/dependent variable.
///
/// \param[in] x Points at which to evaluate interpolant.
///
/// \return Function values of dependent variable \p c evaluated at
/// points \p x in table \p t.
std::vector<double>
interpolate(const InTable& t,
const ResultColumn& c,
const std::vector<double>& x) const;
/// Retrieve connate saturation from all tables.
std::vector<double> connateSat() const;
/// Retrieve critical saturation for particular result column in all
/// tables.
std::vector<double> criticalSat(const ResultColumn& c) const;
/// Retrieve maximum saturation in all tables.
std::vector<double> maximumSat() const;
private:
/// Single tabulated 1D interpolant.
class SingleTable
{
public:
using ElmIt = ECLPropTableRawData::ElementIterator;
/// Constructor.
///
/// \param[in] xBegin Beginning (initial element) of linar range
/// of independent variable values.
///
/// \param[in] xEnd One past the end of linear range of
/// independent variable values.
///
/// \param[in,out] colIt Dependent/column range iterators. On
/// input, point to the beginnings of ranges of results
/// pertinent to a single table. On output, each iterator is
/// advanced across all rows of the SingleTable (including
/// sentinel/invalid nodes) which makes the pointers valid
/// for the next table if relevant (and called in a loop).
SingleTable(ElmIt xBegin,
ElmIt xEnd,
std::vector<ElmIt>& colIt);
/// Evaluate 1D interpolant in sequence of points.
///
/// \param[in] nCols Number of table columns.
///
/// \param[in] c ID of result column/dependent variable.
///
/// \param[in] x Points at which to evaluate interpolant.
///
/// \return Function values of dependent variable \p c evaluated
/// at points \p x.
std::vector<double>
interpolate(const ECLPropTableRawData::SizeType nCols,
const ResultColumn& c,
const std::vector<double>& x) const;
/// Retrieve connate saturation in table.
double connateSat() const;
/// Retrieve critical saturation for particular result column in
/// table.
double criticalSat(const ECLPropTableRawData::SizeType nCols,
const ResultColumn& c) const;
/// Retrieve maximum saturation in table.
double maximumSat() const;
private:
/// Independent variable.
std::vector<double> x_;
/// Dependent variable (or variables). Row major (i.e., C)
/// ordering. Number of elements: x_.size() * host.nCols_.
std::vector<double> y_;
/// Value of dependent variable at position (row,c).
double y(const ECLPropTableRawData::SizeType nCols,
const ECLPropTableRawData::SizeType row,
const ResultColumn& c) const;
};
/// Number of result/dependent variables (== #table cols - 1).
ECLPropTableRawData::SizeType nResCols_;
/// Sequence of individual tables, indexed by *NUM-type vectors.
std::vector<SingleTable> table_;
};
} // namespace Opm
#endif // OPM_ECLPROPTABLE_HEADER_INCLUDED

58
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp vendored
View File

@ -112,7 +112,9 @@ namespace {
result.reserve(x.size());
for (const auto& xi : x) {
result.emplace_back(xi);
// push_back(T(xi)) because vector<bool> does not
// support emplace_back until C++14.
result.push_back(Output(xi));
}
return result;
@ -188,6 +190,17 @@ namespace {
using type = void;
};
/// Translate ERT type class to keyword element type.
///
/// Actual element type of \code ECL_INT_TYPE \endcode.
template <>
struct ElementType<ECL_BOOL_TYPE>
{
/// Element type of ERT Boolean (LOGICAL) data. Stored
/// internally as 'int'.
using type = int;
};
/// Translate ERT type class to keyword element type.
///
/// Actual element type of \code ECL_INT_TYPE \endcode.
@ -226,6 +239,37 @@ namespace {
template <ecl_type_enum Input>
struct ExtractKeywordElements;
/// Extract ERT keyword Boolean (LOGICAL) data.
template <>
struct ExtractKeywordElements<ECL_BOOL_TYPE>
{
using EType = ElementType<ECL_BOOL_TYPE>::type;
/// Function call operator.
///
/// Retrieve actual data elements from ERT keyword of integer
/// (specifically, \c int) type.
///
/// \param[in] kw ERT keyword instance.
///
/// \param[in,out] x Linearised keyword data elements. On
/// input points to memory block of size \code
/// ecl_kw_get_size(kw) * sizeof *x \endcode bytes. On
/// output, those bytes are filled with the actual data
/// values of \p kw.
void operator()(const ecl_kw_type* kw, EType* x) const
{
// 1) Extract raw 'int' values.
ecl_kw_get_memcpy_int_data(kw, x);
// 2) Convert to 'bool'-like values by comparing to
// magic constant ECL_BOOL_TRUE_INT (ecl_util.h).
for (auto n = ecl_kw_get_size(kw), i = 0*n; i < n; ++i) {
x[i] = static_cast<EType>(x[i] == ECL_BOOL_TRUE_INT);
}
}
};
/// Extract ERT keyword integer data.
template <>
struct ExtractKeywordElements<ECL_INT_TYPE>
@ -489,6 +533,10 @@ namespace {
return GetKeywordData<ECL_CHAR_TYPE>::
as<T>(kw, makeStringVector);
case ECL_BOOL_TYPE:
return GetKeywordData<ECL_BOOL_TYPE>::
as<T>(kw, makeStringVector);
case ECL_INT_TYPE:
return GetKeywordData<ECL_INT_TYPE>::
as<T>(kw, makeStringVector);
@ -1569,6 +1617,10 @@ namespace Opm {
ECLRestartData::keywordData<std::string>(const std::string& vector,
const std::string& gridID) const;
template std::vector<bool>
ECLRestartData::keywordData<bool>(const std::string& vector,
const std::string& gridID) const;
template std::vector<int>
ECLRestartData::keywordData<int>(const std::string& vector,
const std::string& gridID) const;
@ -1647,6 +1699,10 @@ namespace Opm {
ECLInitFileData::keywordData<std::string>(const std::string& vector,
const std::string& gridID) const;
template std::vector<bool>
ECLInitFileData::keywordData<bool>(const std::string& vector,
const std::string& gridID) const;
template std::vector<int>
ECLInitFileData::keywordData<int>(const std::string& vector,
const std::string& gridID) const;

1272
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp vendored Normal file
View File

File diff suppressed because it is too large Load Diff

139
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp vendored Normal file
View File

@ -0,0 +1,139 @@
/*
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLSATURATIONFUNC_HEADER_INCLUDED
#define OPM_ECLSATURATIONFUNC_HEADER_INCLUDED
#include <opm/utility/ECLPhaseIndex.hpp>
#include <memory>
#include <vector>
/// \file
///
/// Public interface to relative permeability evaluation machinery. The
/// back-end is aware of ECLIPSE's standard three-phase model for relative
/// permeability of oil and the two- and three-point saturation end-point
/// scaling options. Vertical scaling of relative permeability is not
/// supported at present.
namespace Opm {
class ECLGraph;
class ECLRestartData;
class ECLInitFileData;
/// Gateway to engine for computing relative permeability values based
/// on tabulated saturation functions in ECL output.
class ECLSaturationFunc
{
public:
/// Constructor
///
/// \param[in] G Connected topology of current model's active cells.
/// Needed to linearise region mapping (e.g., SATNUM) that is
/// distributed on local grids to all of the model's active cells
/// (\code member function G.rawLinearisedCellData() \endcode).
///
/// \param[in] init Container of tabulated saturation functions and
/// saturation table end points, if applicable, for all active
/// cells in the model \p G.
///
/// \param[in] useEPS Whether or not to include effects of
/// saturation end-point scaling. No effect if the INIT result
/// set does not actually include saturation end-point scaling
/// data. Otherwise, enables turning EPS off even if associate
/// data is present in the INIT result set.
///
/// Default value (\c true) means that effects of EPS are
/// included if requisite data is present in the INIT result.
ECLSaturationFunc(const ECLGraph& G,
const ECLInitFileData& init,
const bool useEPS = true);
/// Destructor.
~ECLSaturationFunc();
/// Move constructor.
///
/// Subsumes the implementation of an existing object.
///
/// \param[in] rhs Existing engine for saturation function
/// evaluation. Does not have a valid implementation when the
/// constructor completes.
ECLSaturationFunc(ECLSaturationFunc&& rhs);
/// Copy constructor.
///
/// \param[in] rhs Existing engine for saturation function
/// evaluation.
ECLSaturationFunc(const ECLSaturationFunc& rhs);
/// Move assignment operator.
///
/// Subsumes the implementation of an existing object.
///
/// \param[in] rhs Existing engine for saturation function
/// evaluation. Does not have a valid implementation when the
/// constructor completes.
///
/// \return \code *this \endcode.
ECLSaturationFunc& operator=(ECLSaturationFunc&& rhs);
/// Assignment operator.
///
/// \param[in] rhs Existing engine for saturation function
/// evaluation.
///
/// \return \code *this \endcode.
ECLSaturationFunc& operator=(const ECLSaturationFunc& rhs);
/// Compute relative permeability values in all active cells for a
/// single phase.
///
/// \param[in] G Connected topology of current model's active cells.
/// Needed to linearise phase saturations (e.g., SOIL) that are
/// distributed on local grids to all of the model's active cells
/// (\code member function G.rawLinearisedCellData() \endcode).
///
/// \param[in] rstrt ECLIPSE restart vectors. Result set view
/// assumed to be positioned at a particular report step of
/// interest.
///
/// \param[in] p Phase for which to compute relative permeability
/// values.
///
/// \return Derived relative permeability values of active phase \p
/// p for all active cells in model \p G. Empty if phase \p p is
/// not actually active in the current result set.
std::vector<double>
relperm(const ECLGraph& G,
const ECLRestartData& rstrt,
const ECLPhaseIndex p) const;
private:
/// Implementation backend.
class Impl;
/// Pointer to actual backend/engine object.
std::unique_ptr<Impl> pImpl_;
};
} // namespace Opm
#endif // OPM_ECLSATURATIONFUNC_HEADER_INCLUDED

10
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp vendored
View File

@ -554,6 +554,14 @@ namespace {
return errorAcceptable(E.absolute, tol.absolute)
&& errorAcceptable(E.relative, tol.relative);
}
::Opm::ECLGraph
constructGraph(const example::FilePaths& pth)
{
const auto I = ::Opm::ECLInitFileData(pth.init);
return ::Opm::ECLGraph::load(pth.grid, I);
}
} // namespace Anonymous
int main(int argc, char* argv[])
@ -564,7 +572,7 @@ try {
const auto rstrt = ::Opm::ECLRestartData(pth.restart);
const auto steps = availableReportSteps(pth);
const auto graph = example::initGraph(pth);
const auto graph = constructGraph(pth);
auto all_ok = true;
for (const auto& quant : testQuantities(prm)) {

10
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp vendored
View File

@ -206,13 +206,21 @@ namespace {
return ! ((pointMetric(diff) > tol.absolute) ||
(pointMetric(rat) > tol.relative));
}
::Opm::ECLGraph
constructGraph(const example::FilePaths& pth)
{
const auto I = ::Opm::ECLInitFileData(pth.init);
return ::Opm::ECLGraph::load(pth.grid, I);
}
} // namespace Anonymous
int main(int argc, char* argv[])
try {
const auto prm = example::initParam(argc, argv);
const auto pth = example::FilePaths(prm);
const auto G = example::initGraph(pth);
const auto G = constructGraph(pth);
const auto T = G.transmissibility();
const auto ok = transfieldAcceptable(prm, T);

469
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp vendored Normal file
View File

@ -0,0 +1,469 @@
/*
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM 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.
OPM 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 for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include <config.h>
#endif // HAVE_CONFIG_H
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE
#define BOOST_TEST_MODULE TEST_ECLENDPOINTSCALING
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/utility/ECLEndPointScaling.hpp>
#include <exception>
#include <stdexcept>
#include <vector>
namespace {
template <class Collection1, class Collection2>
void check_is_close(const Collection1& c1, const Collection2& c2)
{
BOOST_REQUIRE_EQUAL(c1.size(), c2.size());
if (! c1.empty()) {
auto i1 = c1.begin(), e1 = c1.end();
auto i2 = c2.begin();
for (; i1 != e1; ++i1, ++i2) {
BOOST_CHECK_CLOSE(*i1, *i2, 1.0e-10);
}
}
}
::Opm::SatFunc::EPSEvalInterface::SaturationPoints
associate(const std::vector<double>& s)
{
using SatAssoc = ::Opm::SatFunc::
EPSEvalInterface::SaturationAssoc;
auto sp = ::Opm::SatFunc::
EPSEvalInterface::SaturationPoints{};
sp.reserve(s.size());
for (const auto& si : s) {
sp.push_back(SatAssoc{ 0, si });
}
return sp;
}
} // Namespace Anonymous
// =====================================================================
// Two-point scaling
// ---------------------------------------------------------------------
BOOST_AUTO_TEST_SUITE (TwoPointScaling_FullRange)
BOOST_AUTO_TEST_CASE (NoScaling)
{
namespace SF = ::Opm::SatFunc;
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.0, 0.0, 1.0 };
const auto smin = std::vector<double>{ 0.0 };
const auto smax = std::vector<double>{ 1.0 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_CASE (ScaledConnate)
{
namespace SF = ::Opm::SatFunc;
// Mobile Range: [0.2, 1.0] maps to [ 0.0, 1.0 ]
const auto smin = std::vector<double>{ 0.2 };
const auto smax = std::vector<double>{ 1.0 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.0, 0.0, 1.0 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0,
0,
0.25,
0.5,
0.75,
1.0,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_CASE (ScaledMax)
{
namespace SF = ::Opm::SatFunc;
// Mobile Range: [0.0, 0.8] maps to [ 0.0, 1.0 ]
const auto smin = std::vector<double>{ 0.0 };
const auto smax = std::vector<double>{ 0.8 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.0, 0.0, 1.0 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0,
0.25,
0.5,
0.75,
1.0,
1.0,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_CASE (ScaledBoth)
{
namespace SF = ::Opm::SatFunc;
// Mobile Range: [0.2, 0.8] maps to [ 0.0, 1.0 ]
const auto smin = std::vector<double>{ 0.2 };
const auto smax = std::vector<double>{ 0.8 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.0, 0.0, 1.0 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0,
0.0,
1.0 / 3.0,
2.0 / 3.0,
1.0,
1.0,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_SUITE_END ()
// =====================================================================
BOOST_AUTO_TEST_SUITE (TwoPointScaling_ReducedRange)
BOOST_AUTO_TEST_CASE (NoScaling)
{
namespace SF = ::Opm::SatFunc;
const auto smin = std::vector<double>{ 0.2 };
const auto smax = std::vector<double>{ 0.8 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.2, 0.0, 0.8 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0.2,
0.2,
0.4,
0.6,
0.8,
0.8,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_CASE (ScaledConnate)
{
namespace SF = ::Opm::SatFunc;
// Mobile Range: [0.0, 1.0] maps to [ 0.2, 0.8 ]
// s_eff = 0.6*s + 0.2
const auto smin = std::vector<double>{ 0.0 };
const auto smax = std::vector<double>{ 1.0 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.2, 0.0, 0.8 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0.20,
0.32,
0.44,
0.56,
0.68,
0.80,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_CASE (ScaledMax)
{
namespace SF = ::Opm::SatFunc;
// Mobile Range: [0.2, 0.8] maps to [ 0.0, 1.0 ]
// s_eff = max(0.75*s + 0.05, 0.2)
const auto smin = std::vector<double>{ 0.2 };
const auto smax = std::vector<double>{ 1.0 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.2, 0.0, 0.8 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0.20,
0.20,
0.35,
0.50,
0.65,
0.80,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_CASE (ScaledBoth)
{
namespace SF = ::Opm::SatFunc;
// Mobile Range: [0.2, 0.8] maps to [ 0.5, 0.7 ]
// s_eff = min(max(0.2, 3*s - 13/10), 0.8)
const auto smin = std::vector<double>{ 0.5 };
const auto smax = std::vector<double>{ 0.7 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.2, 0.0, 0.8 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0.2,
0.2,
0.2,
0.5,
0.8,
0.8,
};
const auto eps = SF::TwoPointScaling{ smin, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_SUITE_END ()
// =====================================================================
// Three-point (alternative) scaling, applicable to relperm only.
// ---------------------------------------------------------------------
BOOST_AUTO_TEST_SUITE (ThreePointScaling_FullRange)
BOOST_AUTO_TEST_CASE (NoScaling)
{
namespace SF = ::Opm::SatFunc;
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.0, 0.2, 1.0 };
const auto smin = std::vector<double>{ 0.0 };
const auto sdisp = std::vector<double>{ 0.2 };
const auto smax = std::vector<double>{ 1.0 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto eps = SF::ThreePointScaling{ smin, sdisp, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_CASE (ScaledConnate)
{
namespace SF = ::Opm::SatFunc;
// Mobile Range: [0.4, 1.0] maps to [ 0.0, 1.0 ]
const auto smin = std::vector<double>{ 0.1 };
const auto sdisp = std::vector<double>{ 0.4 };
const auto smax = std::vector<double>{ 1.0 };
const auto tep = SF::EPSEvalInterface::
TableEndPoints { 0.0, 0.2, 1.0 };
const auto s = std::vector<double> {
0.0,
0.2,
0.4,
0.6,
0.8,
1.0,
};
const auto sp = associate(s);
const auto expect = std::vector<double> {
0,
1.0 / 15,
0.2,
7.0 / 15,
11.0 / 15,
1.0,
};
const auto eps = SF::ThreePointScaling{ smin, sdisp, smax };
const auto s_eff = eps.eval(tep, sp);
check_is_close(s_eff, expect);
}
BOOST_AUTO_TEST_SUITE_END ()

2417
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclproptable.cpp vendored Normal file
View File

File diff suppressed because it is too large Load Diff

22
ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp vendored
View File

@ -69,11 +69,13 @@ namespace FlowDiagnostics
/// al. (SPE 146446), Shook and Mitchell (SPE 124625).
Graph flowCapacityStorageCapacityCurve(const Toolbox::Forward& injector_solution,
const Toolbox::Reverse& producer_solution,
const std::vector<double>& pv)
const std::vector<double>& pv,
const double max_pv_fraction)
{
return flowCapacityStorageCapacityCurve(injector_solution.fd.timeOfFlight(),
producer_solution.fd.timeOfFlight(),
pv);
pv,
max_pv_fraction);
}
@ -88,20 +90,30 @@ namespace FlowDiagnostics
/// al. (SPE 146446), Shook and Mitchell (SPE 124625).
Graph flowCapacityStorageCapacityCurve(const std::vector<double>& injector_tof,
const std::vector<double>& producer_tof,
const std::vector<double>& pv)
const std::vector<double>& pv,
const double max_pv_fraction)
{
if (pv.size() != injector_tof.size() || pv.size() != producer_tof.size()) {
throw std::runtime_error("flowCapacityStorageCapacityCurve(): "
"Input solutions must have same size.");
}
// Compute max pv cutoff.
const double total_pv = std::accumulate(pv.begin(), pv.end(), 0.0);
const double max_pv = max_pv_fraction * total_pv;
// Sort according to total travel time.
const int n = pv.size();
typedef std::pair<double, double> D2;
std::vector<D2> time_and_pv(n);
for (int ii = 0; ii < n; ++ii) {
time_and_pv[ii].first = injector_tof[ii] + producer_tof[ii]; // Total travel time.
time_and_pv[ii].second = pv[ii];
if (pv[ii] > max_pv) {
time_and_pv[ii].first = 1e100;
time_and_pv[ii].second = 0.0;
} else {
time_and_pv[ii].first = injector_tof[ii] + producer_tof[ii]; // Total travel time.
time_and_pv[ii].second = pv[ii];
}
}
std::sort(time_and_pv.begin(), time_and_pv.end());

13
ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp vendored
View File

@ -44,11 +44,19 @@ namespace FlowDiagnostics
/// coefficient. For a technical description see Shavali et
/// al. (SPE 146446), Shook and Mitchell (SPE 124625).
///
/// Single cells with a very large pore volume can be filtered out
/// before creating the curve. The 'max_pv_fraction' parameter
/// gives a fraction such that, if a cell's fraction of the total
/// pore volume is above this number, that cell will be
/// ignored. This can be used to disregard numerical aquifers for
/// example.
///
/// Returns F (flow capacity) as a function of Phi (storage capacity),
/// that is for the returned Graph g, g.first is Phi and g.second is F.
Graph flowCapacityStorageCapacityCurve(const Toolbox::Forward& injector_solution,
const Toolbox::Reverse& producer_solution,
const std::vector<double>& pore_volume);
const std::vector<double>& pore_volume,
const double max_pv_fraction = 1.0);
/// This overload gets the injector and producer time-of-flight
/// directly instead of extracting it from the solution
@ -56,7 +64,8 @@ namespace FlowDiagnostics
/// overload.
Graph flowCapacityStorageCapacityCurve(const std::vector<double>& injector_tof,
const std::vector<double>& producer_tof,
const std::vector<double>& pore_volume);
const std::vector<double>& pore_volume,
const double max_pv_fraction = 1.0);
/// The Lorenz coefficient from the F-Phi curve.
///