diff --git a/ApplicationCode/Adm/LicenseInformation.txt b/ApplicationCode/Adm/LicenseInformation.txt index 69a1d1bdea..b7a312a144 100644 --- a/ApplicationCode/Adm/LicenseInformation.txt +++ b/ApplicationCode/Adm/LicenseInformation.txt @@ -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 . + +=============================================================================== + 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. diff --git a/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.cpp b/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.cpp index 64654400c7..9098e6f4e3 100644 --- a/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.cpp +++ b/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.cpp @@ -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; +} + diff --git a/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.h b/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.h index 6fae301645..592d8780fb 100644 --- a/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.h +++ b/ApplicationCode/Commands/FlowCommands/RicPlotProductionRateFeature.h @@ -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); }; diff --git a/ApplicationCode/Commands/FlowCommands/RicShowFlowCharacteristicsPlotFeature.cpp b/ApplicationCode/Commands/FlowCommands/RicShowFlowCharacteristicsPlotFeature.cpp index 3aa4b902d3..9177768572 100644 --- a/ApplicationCode/Commands/FlowCommands/RicShowFlowCharacteristicsPlotFeature.cpp +++ b/ApplicationCode/Commands/FlowCommands/RicShowFlowCharacteristicsPlotFeature.cpp @@ -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(); diff --git a/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp b/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp index 9727079aa2..47f37946ce 100644 --- a/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp +++ b/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp @@ -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; } diff --git a/ApplicationCode/ProjectDataModel/Flow/RimFlowCharacteristicsPlot.cpp b/ApplicationCode/ProjectDataModel/Flow/RimFlowCharacteristicsPlot.cpp index db3622ec06..90acb4e88b 100644 --- a/ApplicationCode/ProjectDataModel/Flow/RimFlowCharacteristicsPlot.cpp +++ b/ApplicationCode/ProjectDataModel/Flow/RimFlowCharacteristicsPlot.cpp @@ -141,14 +141,16 @@ QList RimFlowCharacteristicsPlot::calculateValueOptions( { std::vector 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 ) diff --git a/ApplicationCode/ProjectDataModel/Flow/RimWellAllocationPlot.cpp b/ApplicationCode/ProjectDataModel/Flow/RimWellAllocationPlot.cpp index 89d9ce4244..0d9a0d0861 100644 --- a/ApplicationCode/ProjectDataModel/Flow/RimWellAllocationPlot.cpp +++ b/ApplicationCode/ProjectDataModel/Flow/RimWellAllocationPlot.cpp @@ -357,29 +357,46 @@ std::map *> 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 = "[m3/day]"; + break; + case RigEclipseCaseData::UNITS_FIELD: + unitText = "[Brl/day]"; + break; + case RigEclipseCaseData::UNITS_LAB: + unitText = "[cm3/hr]"; + break; + default: + break; + + } plotTrack->setXAxisTitle("Reservoir Flow Rate " + unitText); } else { + QString unitText; + switch ( unitSet ) + { + case RigEclipseCaseData::UNITS_METRIC: + unitText = "[Liquid Sm3/day], [Gas kSm3/day]"; + break; + case RigEclipseCaseData::UNITS_FIELD: + unitText = "[Liquid BBL/day], [Gas BOE/day]"; + break; + case RigEclipseCaseData::UNITS_LAB: + unitText = "[cm3/hr]"; + break; + default: + break; + + } plotTrack->setXAxisTitle("Surface Flow Rate " + unitText); } } diff --git a/ApplicationCode/ProjectDataModel/Summary/RimSummaryCurveAppearanceCalculator.h b/ApplicationCode/ProjectDataModel/Summary/RimSummaryCurveAppearanceCalculator.h index 1845a97d07..d094221c9b 100644 --- a/ApplicationCode/ProjectDataModel/Summary/RimSummaryCurveAppearanceCalculator.h +++ b/ApplicationCode/ProjectDataModel/Summary/RimSummaryCurveAppearanceCalculator.h @@ -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 mapNameToAppearanceIndex(CurveAppearanceType & appearance, const std::set& 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); diff --git a/ApplicationCode/ReservoirDataModel/RigAccWellFlowCalculator.h b/ApplicationCode/ReservoirDataModel/RigAccWellFlowCalculator.h index ea7d3a6d9f..4c2027cb71 100644 --- a/ApplicationCode/ReservoirDataModel/RigAccWellFlowCalculator.h +++ b/ApplicationCode/ReservoirDataModel/RigAccWellFlowCalculator.h @@ -74,7 +74,6 @@ public: const std::vector& pseudoLengthFromTop(size_t branchIdx) const; const std::vector& trueVerticalDepth(size_t branchIdx) const; const std::vector& accumulatedTracerFlowPrPseudoLength(const QString& tracerName, size_t branchIdx) const; - const std::vector& flowPrPseudoLength( size_t branchIdx) const; const std::vector& tracerFlowPrPseudoLength(const QString& tracerName, size_t branchIdx) const; diff --git a/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h b/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h index b5292e26eb..d7856563fc 100644 --- a/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h +++ b/ApplicationCode/ReservoirDataModel/RigFlowDiagInterfaceTools.h @@ -21,7 +21,7 @@ namespace RigFlowDiagInterfaceTools { template 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); }; diff --git a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp index b9e34ea5c5..a69cdf5192 100644 --- a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp +++ b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp @@ -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)); diff --git a/ResInsightVersion.cmake b/ResInsightVersion.cmake index ca98751bcc..bdbe0b89f3 100644 --- a/ResInsightVersion.cmake +++ b/ResInsightVersion.cmake @@ -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 diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake index e30cdd67a1..18f440371b 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/CMakeLists_files.cmake @@ -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 ) diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp index f29df16e8b..58656e8e7e 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/computeFlowStorageCurve.cpp @@ -38,34 +38,27 @@ try { std::vector start; auto fwd = fdTool.computeInjectionDiagnostics(start); auto rev = fdTool.computeProductionDiagnostics(start); - - // Give disconnected cells zero pore volume. - std::vector 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 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); diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp index a536c2730e..5d79b5f760 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp @@ -31,6 +31,7 @@ #include #include +#include #include #include @@ -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 well_fluxes; Opm::FlowDiagnostics::Toolbox toolbox; bool compute_fluxes_ = false; + bool useEPS_ = false; }; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp new file mode 100644 index 0000000000..6795fdeb02 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp @@ -0,0 +1,1163 @@ +/* + 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 . +*/ + +#if HAVE_CONFIG_H +#include +#endif + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace { + std::vector + unscaledTwoPt(const std::vector& min, + const std::vector& max) + { + assert ((min.size() == max.size()) && "Internal Logic Error"); + assert ((! min.empty()) && "Internal Logic Error"); + + using TEP = Opm::SatFunc::EPSEvalInterface::TableEndPoints; + + auto tep = std::vector{}; + tep.reserve(min.size()); + + for (auto n = min.size(), i = 0*n; i < n; ++i) { + const auto smin = min[i]; + + // Ignore 'sdisp' in the two-point scaling. + tep.push_back(TEP{ smin, smin, max[i] }); + } + + return tep; + } + + std::vector + unscaledThreePt(const std::vector& min , + const std::vector& disp, + const std::vector& max ) + { + assert ((min.size() == max .size()) && "Internal Logic Error"); + assert ((min.size() == disp.size()) && "Internal Logic Error"); + assert ((! min.empty()) && "Internal Logic Error"); + + using TEP = Opm::SatFunc::EPSEvalInterface::TableEndPoints; + + auto tep = std::vector{}; + tep.reserve(min.size()); + + for (auto n = min.size(), i = 0*n; i < n; ++i) { + tep.push_back(TEP{ min[i], disp[i], max[i] }); + } + + return tep; + } +} + +// --------------------------------------------------------------------- +// Class Opm::TwoPointScaling::Impl +// --------------------------------------------------------------------- + +class Opm::SatFunc::TwoPointScaling::Impl +{ +public: + Impl(std::vector smin, + std::vector smax) + : smin_(std::move(smin)) + , smax_(std::move(smax)) + { + if (this->smin_.size() != this->smax_.size()) { + throw std::invalid_argument { + "Size Mismatch Between Minimum and " + "Maximum Saturation Arrays" + }; + } + } + + std::vector + eval(const TableEndPoints& tep, + const SaturationPoints& sp) const; + +private: + std::vector smin_; + std::vector smax_; +}; + +std::vector +Opm::SatFunc::TwoPointScaling:: +Impl::eval(const TableEndPoints& tep, + const SaturationPoints& sp) const +{ + const auto srng = tep.high - tep.low; + + auto effsat = std::vector{}; + effsat.reserve(sp.size()); + + for (const auto& eval_pt : sp) { + const auto cell = eval_pt.cell; + + const auto sLO = this->smin_[cell]; + const auto sHI = this->smax_[cell]; + + effsat.push_back(0.0); + auto& s_eff = effsat.back(); + + if (! (eval_pt.sat > sLO)) { + // s <= sLO + s_eff = tep.low; + } + else if (! (eval_pt.sat < sHI)) { + // s >= sHI + s_eff = tep.high; + } + else { + // s \in (sLO, sHI) + const auto scaled_sat = + ((eval_pt.sat - sLO) / (sHI - sLO)) * srng; + + s_eff = tep.low + scaled_sat; + } + } + + return effsat; +} + +// --------------------------------------------------------------------- +// Class Opm::ThreePointScaling::Impl +// --------------------------------------------------------------------- + +class Opm::SatFunc::ThreePointScaling::Impl +{ +public: + Impl(std::vector smin , + std::vector sdisp, + std::vector smax ) + : smin_ (std::move(smin )) + , sdisp_(std::move(sdisp)) + , smax_ (std::move(smax )) + { + if ((this->sdisp_.size() != this->smin_.size()) || + (this->sdisp_.size() != this->smax_.size())) + { + throw std::invalid_argument { + "Size Mismatch Between Minimum, Displacing " + "and Maximum Saturation Arrays" + }; + } + } + + std::vector + eval(const TableEndPoints& tep, + const SaturationPoints& sp) const; + +private: + std::vector smin_; + std::vector sdisp_; + std::vector smax_; +}; + +std::vector +Opm::SatFunc::ThreePointScaling:: +Impl::eval(const TableEndPoints& tep, + const SaturationPoints& sp) const +{ + auto effsat = std::vector{}; + effsat.reserve(sp.size()); + + for (const auto& eval_pt : sp) { + const auto cell = eval_pt.cell; + + effsat.push_back(0.0); + auto& s_eff = effsat.back(); + + const auto sLO = this->smin_ [cell]; + const auto sR = this->sdisp_[cell]; + const auto sHI = this->smax_ [cell]; + + if (! (eval_pt.sat > sLO)) { + // s <= sLO + s_eff = tep.low; + } + else if (! (eval_pt.sat < sHI)) { + // s >= sHI + s_eff = tep.high; + } + else if (eval_pt.sat < sR) { + // s \in (sLO, sR) + const auto t = (eval_pt.sat - sLO) / (sR - sLO); + + s_eff = tep.low + t*(tep.disp - tep.low); + } + else { + // s \in (sR, sHI) + const auto t = (eval_pt.sat - sR) / (sHI - sR); + + s_eff = tep.disp + t*(tep.high - tep.disp); + } + } + + return effsat; +} + +// --------------------------------------------------------------------- +// EPS factory functions for two-point and three-point scaling options +// --------------------------------------------------------------------- + +namespace Create { + using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions; + using RTEP = ::Opm::SatFunc::CreateEPS::RawTableEndPoints; + using TEP = ::Opm::SatFunc::EPSEvalInterface::TableEndPoints; + + namespace TwoPoint { + using EPS = ::Opm::SatFunc::TwoPointScaling; + using EPSPtr = std::unique_ptr; + + struct Kr { + static EPSPtr + G(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + + static EPSPtr + OG(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + + static EPSPtr + OW(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + + static EPSPtr + W(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + }; + + struct Pc { + static EPSPtr + GO(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + + static EPSPtr + OW(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + }; + + static EPSPtr + scalingFunction(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init, + const EPSOpt& opt); + + static std::vector + unscaledEndPoints(const RTEP& ep, const EPSOpt& opt); + } // namespace TwoPoint + + namespace ThreePoint { + using EPS = ::Opm::SatFunc::ThreePointScaling; + using EPSPtr = std::unique_ptr; + + struct Kr { + static EPSPtr + G(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + + static EPSPtr + OG(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + + static EPSPtr + OW(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + + static EPSPtr + W(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init); + }; + + static EPSPtr + scalingFunction(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init, + const ::Opm::SatFunc::CreateEPS::EPSOptions& opt); + + static std::vector + unscaledEndPoints(const RTEP& ep, const EPSOpt& opt); + } // namespace ThreePoint +} // namespace Create + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Implementation of Create::TwoPoint::scalingFunction() +Create::TwoPoint::EPSPtr +Create::TwoPoint::Kr::G(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto sgcr = G.rawLinearisedCellData(init, "SGCR"); + auto sgu = G.rawLinearisedCellData(init, "SGU"); + + if ((sgcr.size() != sgu.size()) || + (sgcr.size() != G.numCells())) + { + throw std::invalid_argument { + "Missing or Mismatching Gas End-Point " + "Specifications (SGCR and/or SGU)" + }; + } + + return EPSPtr { new EPS { std::move(sgcr), std::move(sgu) } }; +} + +Create::TwoPoint::EPSPtr +Create::TwoPoint::Kr::OG(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto sogcr = G.rawLinearisedCellData(init, "SOGCR"); + + if (sogcr.size() != G.numCells()) { + throw std::invalid_argument { + "Missing or Mismatching Critical Oil " + "Saturation in Oil/Gas System" + }; + } + + auto smax = std::vector(sogcr.size(), 1.0); + + // Adjust maximum S_o for scaled connate gas saturations. + { + const auto sgl = G.rawLinearisedCellData(init, "SGL"); + + if (sgl.size() != sogcr.size()) { + throw std::invalid_argument { + "Missing or Mismatching Connate Gas " + "Saturation in Oil/Gas System" + }; + } + + for (auto n = sgl.size(), i = 0*n; i < n; ++i) { + smax[i] -= sgl[i]; + } + } + + // Adjust maximum S_o for scaled connate water saturations (if relevant). + { + const auto swl = G.rawLinearisedCellData(init, "SWL"); + + if (swl.size() == sogcr.size()) { + for (auto n = swl.size(), i = 0*n; i < n; ++i) { + smax[i] -= swl[i]; + } + } + else if (! swl.empty()) { + throw std::invalid_argument { + "Mismatching Connate Water " + "Saturation in Oil/Gas System" + }; + } + } + + return EPSPtr { new EPS { std::move(sogcr), std::move(smax) } }; +} + +Create::TwoPoint::EPSPtr +Create::TwoPoint::Kr::OW(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto sowcr = G.rawLinearisedCellData(init, "SOWCR"); + + if (sowcr.size() != G.numCells()) { + throw std::invalid_argument { + "Missing or Mismatching Critical Oil " + "Saturation in Oil/Water System" + }; + } + + auto smax = std::vector(sowcr.size(), 1.0); + + // Adjust maximum S_o for scaled connate water saturations. + { + const auto swl = G.rawLinearisedCellData(init, "SWL"); + + if (swl.size() != sowcr.size()) { + throw std::invalid_argument { + "Missing or Mismatching Connate Water " + "Saturation in Oil/Water System" + }; + } + + for (auto n = swl.size(), i = 0*n; i < n; ++i) { + smax[i] -= swl[i]; + } + } + + // Adjust maximum S_o for scaled connate gas saturations (if relevant). + { + const auto sgl = G.rawLinearisedCellData(init, "SGL"); + + if (sgl.size() == sowcr.size()) { + for (auto n = sgl.size(), i = 0*n; i < n; ++i) { + smax[i] -= sgl[i]; + } + } + else if (! sgl.empty()) { + throw std::invalid_argument { + "Mismatching Connate Gas " + "Saturation in Oil/Water System" + }; + } + } + + return EPSPtr { new EPS { std::move(sowcr), std::move(smax) } }; +} + +Create::TwoPoint::EPSPtr +Create::TwoPoint::Kr::W(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto swcr = G.rawLinearisedCellData(init, "SWCR"); + auto swu = G.rawLinearisedCellData(init, "SWU"); + + if (swcr.empty() || swu.empty()) { + throw std::invalid_argument { + "Missing Water End-Point Specifications (SWCR and/or SWU)" + }; + } + + return EPSPtr { new EPS { std::move(swcr), std::move(swu) } }; +} + +Create::TwoPoint::EPSPtr +Create::TwoPoint::Pc::GO(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto sgl = G.rawLinearisedCellData(init, "SGL"); + auto sgu = G.rawLinearisedCellData(init, "SGU"); + + if ((sgl.size() != sgu.size()) || + (sgl.size() != G.numCells())) + { + throw std::invalid_argument { + "Missing or Mismatching Connate or Maximum Gas " + "Saturation in Pcgo EPS" + }; + } + + return EPSPtr { new EPS { std::move(sgl), std::move(sgu) } }; +} + +Create::TwoPoint::EPSPtr +Create::TwoPoint::Pc::OW(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto swl = G.rawLinearisedCellData(init, "SWL"); + auto swu = G.rawLinearisedCellData(init, "SWU"); + + if ((swl.size() != swu.size()) || + (swl.size() != G.numCells())) + { + throw std::invalid_argument { + "Missing or Mismatching Connate or Maximum Water " + "Saturation in Pcow EPS" + }; + } + + return EPSPtr { new EPS { std::move(swl), std::move(swu) } }; +} + +Create::TwoPoint::EPSPtr +Create::TwoPoint:: +scalingFunction(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init, + const ::Opm::SatFunc::CreateEPS::EPSOptions& opt) +{ + using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; + using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; + using PhIdx = ::Opm::ECLPhaseIndex; + + assert (((! opt.use3PtScaling) || (opt.curve == FCat::CapPress)) + && "Internal Error Selecting EPS Family"); + + if (opt.curve == FCat::Relperm) { + if (opt.subSys == SSys::OilWater) { + if (opt.thisPh == PhIdx::Vapour) { + throw std::invalid_argument { + "Cannot Create an EPS for Gas Relperm " + "in an Oil/Water System" + }; + } + + if (opt.thisPh == PhIdx::Aqua) { + return Create::TwoPoint::Kr::W(G, init); + } + + return Create::TwoPoint::Kr::OW(G, init); + } + + if (opt.subSys == SSys::OilGas) { + if (opt.thisPh == PhIdx::Aqua) { + throw std::invalid_argument { + "Cannot Create an EPS for Water Relperm " + "in an Oil/Gas System" + }; + } + + if (opt.thisPh == PhIdx::Vapour) { + return Create::TwoPoint::Kr::G(G, init); + } + + return Create::TwoPoint::Kr::OG(G, init); + } + } + + if (opt.curve == FCat::CapPress) { + if (opt.thisPh == PhIdx::Liquid) { + throw std::invalid_argument { + "Creating Capillary Pressure EPS as a Function " + "of Oil Saturation is not Supported" + }; + } + + if (opt.thisPh == PhIdx::Vapour) { + return Create::TwoPoint::Pc::GO(G, init); + } + + if (opt.thisPh == PhIdx::Aqua) { + return Create::TwoPoint::Pc::OW(G, init); + } + } + + // Invalid. + return EPSPtr{}; +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Implementation of Create::TwoPoint::unscaledEndPoints() +std::vector +Create::TwoPoint::unscaledEndPoints(const RTEP& ep, const EPSOpt& opt) +{ + using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; + using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; + using PhIdx = ::Opm::ECLPhaseIndex; + + assert (((opt.curve == FCat::CapPress) || + (! opt.use3PtScaling)) && "Internal Logic Error"); + + if (opt.curve == FCat::CapPress) { + // Left node is connate saturation, right node is max saturation. + + if (opt.thisPh == PhIdx::Liquid) { + throw std::invalid_argument { + "No Capillary Pressure Function for Oil" + }; + } + + if (opt.thisPh == PhIdx::Aqua) { + return unscaledTwoPt(ep.conn.water, ep.smax.water); + } + + if (opt.thisPh == PhIdx::Vapour) { + return unscaledTwoPt(ep.conn.gas, ep.smax.gas); + } + } + + if (opt.curve == FCat::Relperm) { + // Left node is critical saturation, right node is max saturation. + + if (opt.subSys == SSys::OilGas) { + if (opt.thisPh == PhIdx::Aqua) { + throw std::invalid_argument { + "Void Request for Unscaled Water Saturation " + "End-Points in Oil-Gas System" + }; + } + + if (opt.thisPh == PhIdx::Liquid) { + return unscaledTwoPt(ep.crit.oil_in_gas, ep.smax.oil); + } + + if (opt.thisPh == PhIdx::Vapour) { + return unscaledTwoPt(ep.crit.gas, ep.smax.gas); + } + } + + if (opt.subSys == SSys::OilWater) { + if (opt.thisPh == PhIdx::Aqua) { + return unscaledTwoPt(ep.crit.water, ep.smax.water); + } + + if (opt.thisPh == PhIdx::Liquid) { + return unscaledTwoPt(ep.crit.oil_in_water, ep.smax.oil); + } + + if (opt.thisPh == PhIdx::Vapour) { + throw std::invalid_argument { + "Void Request for Unscaled Gas Saturation " + "End-Points in Oil-Water System" + }; + } + } + } + + // Invalid. + return {}; +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Implementation of Create::ThreePoint::scalingFunction() + +Create::ThreePoint::EPSPtr +Create::ThreePoint::Kr::G(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto sgcr = G.rawLinearisedCellData(init, "SGCR"); + auto sgu = G.rawLinearisedCellData(init, "SGU"); + + if ((sgcr.size() != sgu.size()) || + (sgcr.size() != G.numCells())) + { + throw std::invalid_argument { + "Missing or Mismatching Gas End-Point " + "Specifications (SGCR and/or SGU)" + }; + } + + auto sr = std::vector(G.numCells(), 1.0); + + // Adjust displacing saturation for connate water. + { + const auto swl = G.rawLinearisedCellData(init, "SWL"); + + if (swl.size() == sgcr.size()) { + for (auto n = swl.size(), i = 0*n; i < n; ++i) { + sr[i] -= swl[i]; + } + } + else if (! swl.empty()) { + throw std::invalid_argument { + "Connate Water Saturation Array Mismatch " + "in Three-Point Scaling Option" + }; + } + } + + // Adjust displacing saturation for critical S_o in O/G system. + { + const auto sogcr = G.rawLinearisedCellData(init, "SOGCR"); + + if (sogcr.size() == sgcr.size()) { + for (auto n = sogcr.size(), i = 0*n; i < n; ++i) { + sr[i] -= sogcr[i]; + } + } + else if (! sogcr.empty()) { + throw std::invalid_argument { + "Critical Oil Saturation (O/G System) Array " + "Size Mismatch in Three-Point Scaling Option" + }; + } + } + + return EPSPtr { + new EPS { std::move(sgcr), std::move(sr), std::move(sgu) } + }; +} + +Create::ThreePoint::EPSPtr +Create::ThreePoint::Kr::OG(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto sogcr = G.rawLinearisedCellData(init, "SOGCR"); + + if (sogcr.size() != G.numCells()) { + throw std::invalid_argument { + "Missing or Mismatching Critical Oil " + "Saturation in Oil/Gas System" + }; + } + + auto smax = std::vector(sogcr.size(), 1.0); + + // Adjust maximum S_o for scaled connate gas saturations. + { + const auto sgl = G.rawLinearisedCellData(init, "SGL"); + + if (sgl.size() != sogcr.size()) { + throw std::invalid_argument { + "Missing or Mismatching Connate Gas " + "Saturation in Oil/Gas System" + }; + } + + for (auto n = sgl.size(), i = 0*n; i < n; ++i) { + smax[i] -= sgl[i]; + } + } + + auto sdisp = std::vector(sogcr.size(), 1.0); + + // Adjust displacing S_o for scaled critical gas saturation. + { + const auto sgcr = G.rawLinearisedCellData(init, "SGCR"); + + if (sgcr.size() != sogcr.size()) { + throw std::invalid_argument { + "Missing or Mismatching Scaled Critical Gas " + "Saturation in Oil/Gas System" + }; + } + + for (auto n = sgcr.size(), i = 0*n; i < n; ++i) { + sdisp[i] -= sgcr[i]; + } + } + + // Adjust displacing and maximum S_o for scaled connate water + // saturations (if relevant). + { + const auto swl = G.rawLinearisedCellData(init, "SWL"); + + if (swl.size() == sogcr.size()) { + for (auto n = swl.size(), i = 0*n; i < n; ++i) { + sdisp[i] -= swl[i]; + smax [i] -= swl[i]; + } + } + else if (! swl.empty()) { + throw std::invalid_argument { + "Mismatching Scaled Connate Water " + "Saturation in Oil/Gas System" + }; + } + } + + return EPSPtr { + new EPS { std::move(sogcr), std::move(sdisp), std::move(smax) } + }; +} + +Create::ThreePoint::EPSPtr +Create::ThreePoint::Kr::OW(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto sowcr = G.rawLinearisedCellData(init, "SOWCR"); + + if (sowcr.size() != G.numCells()) { + throw std::invalid_argument { + "Missing or Mismatching Critical Oil " + "Saturation in Oil/Water System" + }; + } + + auto smax = std::vector(sowcr.size(), 1.0); + + // Adjust maximum S_o for scaled connate water saturations. + { + const auto swl = G.rawLinearisedCellData(init, "SWL"); + + if (swl.size() != sowcr.size()) { + throw std::invalid_argument { + "Missing or Mismatching Connate Water " + "Saturation in Oil/Water System" + }; + } + + for (auto n = swl.size(), i = 0*n; i < n; ++i) { + smax[i] -= swl[i]; + } + } + + auto sdisp = std::vector(sowcr.size(), 1.0); + + // Adjust displacing S_o for scaled critical water saturations. + { + const auto swcr = G.rawLinearisedCellData(init, "SWCR"); + + if (swcr.size() != sowcr.size()) { + throw std::invalid_argument { + "Missing or Mismatching Scaled Critical Water " + "Saturation in Oil/Water System" + }; + } + + for (auto n = swcr.size(), i = 0*n; i < n; ++i) { + sdisp[i] -= swcr[i]; + } + } + + // Adjust displacing and maximum S_o for scaled connate gas saturations + // (if relevant). + { + const auto sgl = G.rawLinearisedCellData(init, "SGL"); + + if (sgl.size() == sowcr.size()) { + for (auto n = sgl.size(), i = 0*n; i < n; ++i) { + sdisp[i] -= sgl[i]; + smax [i] -= sgl[i]; + } + } + else if (! sgl.empty()) { + throw std::invalid_argument { + "Mismatching Connate Gas " + "Saturation in Oil/Water System" + }; + } + } + + return EPSPtr { + new EPS { std::move(sowcr), std::move(sdisp), std::move(smax) } + }; +} + +Create::ThreePoint::EPSPtr +Create::ThreePoint::Kr::W(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init) +{ + auto swcr = G.rawLinearisedCellData(init, "SWCR"); + auto swu = G.rawLinearisedCellData(init, "SWU"); + + if ((swcr.size() != G.numCells()) || + (swcr.size() != swu.size())) + { + throw std::invalid_argument { + "Missing Water End-Point Specifications (SWCR and/or SWU)" + }; + } + + auto sdisp = std::vector(swcr.size(), 1.0); + + // Adjust displacing S_w for scaled critical oil saturation. + { + const auto sowcr = G.rawLinearisedCellData(init, "SOWCR"); + + if (sowcr.size() == swcr.size()) { + for (auto n = sowcr.size(), i = 0*n; i < n; ++i) { + sdisp[i] -= sowcr[i]; + } + } + else if (! sowcr.empty()) { + throw std::invalid_argument { + "Missing or Mismatching Scaled Critical " + "Oil Saturation in Oil/Water System" + }; + } + } + + // Adjust displacing S_w for scaled connate gas saturation. + { + const auto sgl = G.rawLinearisedCellData(init, "SGL"); + + if (sgl.size() == swcr.size()) { + for (auto n = sgl.size(), i = 0*n; i < n; ++i) { + sdisp[i] -= sgl[i]; + } + } + else if (! sgl.empty()) { + throw std::invalid_argument { + "Missing or Mismatching Scaled Connate " + "Gas Saturation in Oil/Water System" + }; + } + } + + return EPSPtr { + new EPS { std::move(swcr), std::move(sdisp), std::move(swu) } + }; +} + +Create::ThreePoint::EPSPtr +Create::ThreePoint:: +scalingFunction(const ::Opm::ECLGraph& G, + const ::Opm::ECLInitFileData& init, + const ::Opm::SatFunc::CreateEPS::EPSOptions& opt) +{ + using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; + using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; + using PhIdx = ::Opm::ECLPhaseIndex; + + assert ((opt.use3PtScaling && (opt.curve == FCat::Relperm)) + && "Internal Error Selecting EPS Family"); + + if (opt.subSys == SSys::OilWater) { + if (opt.thisPh == PhIdx::Vapour) { + throw std::invalid_argument { + "Cannot Create a Three-Point EPS for " + "Gas Relperm in an Oil/Water System" + }; + } + + if (opt.thisPh == PhIdx::Aqua) { + return Create::ThreePoint::Kr::W(G, init); + } + + return Create::ThreePoint::Kr::OW(G, init); + } + + if (opt.subSys == SSys::OilGas) { + if (opt.thisPh == PhIdx::Aqua) { + throw std::invalid_argument { + "Cannot Create a Three-Point EPS for " + "Water Relperm in an Oil/Gas System" + }; + } + + if (opt.thisPh == PhIdx::Vapour) { + return Create::ThreePoint::Kr::G(G, init); + } + + return Create::ThreePoint::Kr::OG(G, init); + } + + // Invalid. + return EPSPtr{}; +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Implementation of Create::ThreePoint::unscaledEndPoints() +std::vector +Create::ThreePoint::unscaledEndPoints(const RTEP& ep, const EPSOpt& opt) +{ + using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; + using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; + using PhIdx = ::Opm::ECLPhaseIndex; + + assert ((opt.use3PtScaling && (opt.curve == FCat::Relperm)) + && "Internal Error Selecting EPS Family"); + + auto sdisp = [](const std::vector& s1, + const std::vector& s2) + -> std::vector + { + auto sr = std::vector(s1.size(), 1.0); + + for (auto n = s1.size(), i = 0*n; i < n; ++i) { + sr[i] -= s1[i] + s2[i]; + } + + return sr; + }; + + // Left node is critical saturation, middle node is displacing critical + // saturation, and right node is maximum saturation. + + if (opt.subSys == SSys::OilGas) { + if (opt.thisPh == PhIdx::Aqua) { + throw std::invalid_argument { + "Void Request for Unscaled Water Saturation " + "End-Points in Oil-Gas System" + }; + } + + if (opt.thisPh == PhIdx::Liquid) { + return unscaledThreePt(ep.crit.oil_in_gas, + sdisp(ep.crit.gas, ep.conn.water), + ep.smax.oil); + } + + if (opt.thisPh == PhIdx::Vapour) { + return unscaledThreePt(ep.crit.gas, + sdisp(ep.crit.oil_in_gas, ep.conn.water), + ep.smax.gas); + } + } + + if (opt.subSys == SSys::OilWater) { + if (opt.thisPh == PhIdx::Aqua) { + return unscaledThreePt(ep.crit.water, + sdisp(ep.crit.oil_in_water, ep.conn.gas), + ep.smax.water); + } + + if (opt.thisPh == PhIdx::Liquid) { + return unscaledThreePt(ep.crit.oil_in_water, + sdisp(ep.crit.water, ep.conn.gas), + ep.smax.oil); + } + + if (opt.thisPh == PhIdx::Vapour) { + throw std::invalid_argument { + "Void Request for Unscaled Gas Saturation " + "End-Points in Oil-Water System" + }; + } + } + + // Invalid. + return {}; +} + +// ##################################################################### +// ===================================================================== +// Public Interface Below Separator +// ===================================================================== +// ##################################################################### + +// Class Opm::SatFunc::EPSEvalInterface +Opm::SatFunc::EPSEvalInterface::~EPSEvalInterface() +{} + +// --------------------------------------------------------------------- + +// Class Opm::SatFunc::TwoPointScaling +Opm::SatFunc::TwoPointScaling:: +TwoPointScaling(std::vector smin, + std::vector smax) + : pImpl_(new Impl(std::move(smin), std::move(smax))) +{} + +Opm::SatFunc::TwoPointScaling::~TwoPointScaling() +{} + +Opm::SatFunc::TwoPointScaling:: +TwoPointScaling(const TwoPointScaling& rhs) + : pImpl_(new Impl(*rhs.pImpl_)) +{} + +Opm::SatFunc::TwoPointScaling:: +TwoPointScaling(TwoPointScaling&& rhs) + : pImpl_(std::move(rhs.pImpl_)) +{} + +Opm::SatFunc::TwoPointScaling& +Opm::SatFunc::TwoPointScaling::operator=(const TwoPointScaling& rhs) +{ + this->pImpl_.reset(new Impl(*rhs.pImpl_)); + + return *this; +} + +Opm::SatFunc::TwoPointScaling& +Opm::SatFunc::TwoPointScaling::operator=(TwoPointScaling&& rhs) +{ + this->pImpl_ = std::move(rhs.pImpl_); + + return *this; +} + +std::vector +Opm::SatFunc::TwoPointScaling::eval(const TableEndPoints& tep, + const SaturationPoints& sp) const +{ + return this->pImpl_->eval(tep, sp); +} + +std::unique_ptr +Opm::SatFunc::TwoPointScaling::clone() const +{ + return std::unique_ptr(new TwoPointScaling(*this)); +} + +// --------------------------------------------------------------------- + +// Class Opm::SatFunc::ThreePointScaling +Opm::SatFunc::ThreePointScaling:: +ThreePointScaling(std::vector smin, + std::vector sdisp, + std::vector smax) + : pImpl_(new Impl(std::move(smin), std::move(sdisp), std::move(smax))) +{} + +Opm::SatFunc::ThreePointScaling::~ThreePointScaling() +{} + +Opm::SatFunc::ThreePointScaling:: +ThreePointScaling(const ThreePointScaling& rhs) + : pImpl_(new Impl(*rhs.pImpl_)) +{} + +Opm::SatFunc::ThreePointScaling::ThreePointScaling(ThreePointScaling&& rhs) + : pImpl_(std::move(rhs.pImpl_)) +{} + +Opm::SatFunc::ThreePointScaling& +Opm::SatFunc::ThreePointScaling::operator=(const ThreePointScaling& rhs) +{ + this->pImpl_.reset(new Impl(*rhs.pImpl_)); + + return *this; +} + +Opm::SatFunc::ThreePointScaling& +Opm::SatFunc::ThreePointScaling::operator=(ThreePointScaling&& rhs) +{ + this->pImpl_ = std::move(rhs.pImpl_); + + return *this; +} + +std::vector +Opm::SatFunc::ThreePointScaling::eval(const TableEndPoints& tep, + const SaturationPoints& sp) const +{ + return this->pImpl_->eval(tep, sp); +} + +std::unique_ptr +Opm::SatFunc::ThreePointScaling::clone() const +{ + return std::unique_ptr(new ThreePointScaling(*this)); +} + +// --------------------------------------------------------------------- +// Factory function Opm::SatFunc::CreateEPS::fromECLOutput() + +std::unique_ptr +Opm::SatFunc::CreateEPS:: +fromECLOutput(const ECLGraph& G, + const ECLInitFileData& init, + const EPSOptions& opt) +{ + if ((opt.curve == FunctionCategory::CapPress) || + (! opt.use3PtScaling)) + { + return Create::TwoPoint::scalingFunction(G, init, opt); + } + + if ((opt.curve == FunctionCategory::Relperm) && opt.use3PtScaling) + { + return Create::ThreePoint::scalingFunction(G, init, opt); + } + + // Invalid + return std::unique_ptr{}; +} + +// --------------------------------------------------------------------- +// Factory function Opm::SatFunc::CreateEPS::unscaledEndPoints() + +std::vector +Opm::SatFunc::CreateEPS:: +unscaledEndPoints(const RawTableEndPoints& ep, + const EPSOptions& opt) +{ + if ((opt.curve == FunctionCategory::CapPress) || + (! opt.use3PtScaling)) + { + return Create::TwoPoint::unscaledEndPoints(ep, opt); + } + + if ((opt.curve == FunctionCategory::Relperm) && opt.use3PtScaling) + { + return Create::ThreePoint::unscaledEndPoints(ep, opt); + } + + // Invalid + return {}; +} diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp new file mode 100644 index 0000000000..149119b904 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp @@ -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 . +*/ + +#ifndef OPM_ECLENDPOINTSCALING_HEADER_INCLUDED +#define OPM_ECLENDPOINTSCALING_HEADER_INCLUDED + +#include + +#include +#include +#include +#include + +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::size_type cell; + + /// Saturation value. + double sat; + }; + + /// Convenience type alias. + using SaturationPoints = std::vector; + + /// 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 + eval(const TableEndPoints& tep, + const SaturationPoints& sp) const = 0; + + virtual std::unique_ptr 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 smin, + std::vector 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 + eval(const TableEndPoints& tep, + const SaturationPoints& sp) const override; + + virtual std::unique_ptr clone() const override; + + private: + /// Implementation class. + class Impl; + + /// Pimpl idiom. + std::unique_ptr 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 smin, + std::vector sdisp, + std::vector 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 + eval(const TableEndPoints& tep, + const SaturationPoints& sp) const override; + + virtual std::unique_ptr clone() const override; + + private: + /// Implementation class. + class Impl; + + /// Pimpl idiom. + std::unique_ptr 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 oil; + + /// Connate gas saturation for each table in total set of + /// tabulated saturation functions. + std::vector gas; + + /// Connate water saturation for each table in total set of + /// tabulated saturation functions. + std::vector 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 oil_in_gas; + + /// Critical oil saturation in 2p OW system from total set + /// of tabulated saturation functions. + std::vector oil_in_water; + + /// Critical gas saturation in 2p OG or 3p OGW system from + /// total set of tabulated saturation functions. + std::vector gas; + + /// Critical water saturation in 2p OW or 3p OGW system from + /// total set of tabulated saturation functions. + std::vector water; + }; + + /// Collection of maximum saturation end points. + struct Maximum { + /// Maximum oil saturation for each table in total set of + /// tabulated saturation functions. + std::vector oil; + + /// Maximum gas saturation for each table in total set of + /// tabulated saturation functions. + std::vector gas; + + /// Maximum water saturation for each table in total set of + /// tabulated saturation functions. + std::vector 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 + 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 + unscaledEndPoints(const RawTableEndPoints& ep, + const EPSOptions& opt); + }; +}} // namespace Opm::SatFunc + +#endif // OPM_ECLENDPOINTSCALING_HEADER_INCLUDED diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp index 1e597a4064..712c433e60 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp @@ -22,11 +22,15 @@ #include +#include + 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 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 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]); } diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp index 3f230ab23f..9e7b6ac6c7 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp @@ -21,6 +21,9 @@ #define OPM_ECLFLUXCALC_HEADER_INCLUDED #include +#include +#include + #include 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 flux(const ECLRestartData& rstrt, - const PhaseIndex phase) const; + const ECLPhaseIndex phase) const; private: struct DynamicData { std::vector pressure; + std::vector relperm; }; double singleFlux(const int connection, const DynamicData& dyn_data) const; const ECLGraph& graph_; + ECLSaturationFunc satfunc_; std::vector neighbours_; std::vector transmissibility_; }; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp index 52d94ced43..f56f880c35 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.cpp @@ -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& activePhases() const; + const std::vector& activePhases() const; /// Retrieve the simulation scenario's set of active grids. /// @@ -1322,7 +1322,7 @@ public: /// all). std::vector 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 activePhases_; + std::vector activePhases_; /// Set of active grids in result set. std::vector 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& +const std::vector& Opm::ECLGraph::Impl::activePhases() const { return this->activePhases_; @@ -2028,7 +2028,7 @@ Opm::ECLGraph::Impl::transmissibility() const std::vector 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(ih[INTEHEAD_PHASE_INDEX]); - using Check = std::pair; + using Check = std::pair; const auto allPhases = std::vector { - { 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& +const std::vector& Opm::ECLGraph::activePhases() const { return this->pImpl_->activePhases(); @@ -2351,7 +2351,7 @@ std::vector Opm::ECLGraph::transmissibility() const std::vector Opm::ECLGraph::flux(const ECLRestartData& rstrt, - const PhaseIndex phase) const + const ECLPhaseIndex phase) const { return this->pImpl_->flux(rstrt, phase); } diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp index ff60cf67dc..4ad8bc7490 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLGraph.hpp @@ -20,6 +20,7 @@ #ifndef OPM_ECLGRAPH_HEADER_INCLUDED #define OPM_ECLGRAPH_HEADER_INCLUDED +#include #include #include @@ -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& activePhases() const; + const std::vector& activePhases() const; /// Retrieve the simulation scenario's set of active grids. /// @@ -167,7 +165,7 @@ namespace Opm { /// (rm^3/s). std::vector 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. diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPhaseIndex.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPhaseIndex.hpp new file mode 100644 index 0000000000..0ffa899a78 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPhaseIndex.hpp @@ -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 . +*/ + +#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 diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.cpp new file mode 100644 index 0000000000..b49b3b2696 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.cpp @@ -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 . +*/ + +#include + +#include +#include +#include +#include +#include +#include +#include + +Opm::SatFuncInterpolant::SingleTable:: +SingleTable(ElmIt xBegin, + ElmIt xEnd, + std::vector& 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_castx_.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 +Opm::SatFuncInterpolant::SingleTable:: +interpolate(const ECLPropTableRawData::SizeType nCols, + const ResultColumn& c, + const std::vector& x) const +{ + auto y = std::vector{}; 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{ 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 +Opm::SatFuncInterpolant::interpolate(const InTable& t, + const ResultColumn& c, + const std::vector& 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 +Opm::SatFuncInterpolant::connateSat() const +{ + auto sconn = std::vector{}; + sconn.reserve(this->table_.size()); + + for (const auto& t : this->table_) { + sconn.push_back(t.connateSat()); + } + + return sconn; +} + +std::vector +Opm::SatFuncInterpolant::criticalSat(const ResultColumn& c) const +{ + auto scrit = std::vector{}; + scrit.reserve(this->table_.size()); + + for (const auto& t : this->table_) { + scrit.push_back(t.criticalSat(this->nResCols_, c)); + } + + return scrit; +} + +std::vector +Opm::SatFuncInterpolant::maximumSat() const +{ + auto smax = std::vector{}; + smax.reserve(this->table_.size()); + + for (const auto& t : this->table_) { + smax.push_back(t.maximumSat()); + } + + return smax; +} diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.hpp new file mode 100644 index 0000000000..32a7a81106 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPropTable.hpp @@ -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 . +*/ + +#ifndef OPM_ECLPROPTABLE_HEADER_INCLUDED +#define OPM_ECLPROPTABLE_HEADER_INCLUDED + +#include + +/// \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; + + /// 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 + interpolate(const InTable& t, + const ResultColumn& c, + const std::vector& x) const; + + /// Retrieve connate saturation from all tables. + std::vector connateSat() const; + + /// Retrieve critical saturation for particular result column in all + /// tables. + std::vector criticalSat(const ResultColumn& c) const; + + /// Retrieve maximum saturation in all tables. + std::vector 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& 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 + interpolate(const ECLPropTableRawData::SizeType nCols, + const ResultColumn& c, + const std::vector& 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 x_; + + /// Dependent variable (or variables). Row major (i.e., C) + /// ordering. Number of elements: x_.size() * host.nCols_. + std::vector 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 table_; + }; + +} // namespace Opm + +#endif // OPM_ECLPROPTABLE_HEADER_INCLUDED diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp index 2add62a4e6..8d2cff11bb 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp @@ -112,7 +112,9 @@ namespace { result.reserve(x.size()); for (const auto& xi : x) { - result.emplace_back(xi); + // push_back(T(xi)) because vector 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 + { + /// 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 struct ExtractKeywordElements; + /// Extract ERT keyword Boolean (LOGICAL) data. + template <> + struct ExtractKeywordElements + { + using EType = ElementType::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(x[i] == ECL_BOOL_TRUE_INT); + } + } + }; + /// Extract ERT keyword integer data. template <> struct ExtractKeywordElements @@ -489,6 +533,10 @@ namespace { return GetKeywordData:: as(kw, makeStringVector); + case ECL_BOOL_TYPE: + return GetKeywordData:: + as(kw, makeStringVector); + case ECL_INT_TYPE: return GetKeywordData:: as(kw, makeStringVector); @@ -1569,6 +1617,10 @@ namespace Opm { ECLRestartData::keywordData(const std::string& vector, const std::string& gridID) const; + template std::vector + ECLRestartData::keywordData(const std::string& vector, + const std::string& gridID) const; + template std::vector ECLRestartData::keywordData(const std::string& vector, const std::string& gridID) const; @@ -1647,6 +1699,10 @@ namespace Opm { ECLInitFileData::keywordData(const std::string& vector, const std::string& gridID) const; + template std::vector + ECLInitFileData::keywordData(const std::string& vector, + const std::string& gridID) const; + template std::vector ECLInitFileData::keywordData(const std::string& vector, const std::string& gridID) const; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp new file mode 100644 index 0000000000..9c01c77cc8 --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp @@ -0,0 +1,1272 @@ +/* + 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 . +*/ + +#if HAVE_CONFIG_H +#include +#endif + +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +namespace { + std::vector + oil_saturation(const std::vector& sg, + const std::vector& sw, + const ::Opm::ECLGraph& G, + const ::Opm::ECLRestartData& rstrt) + { + auto so = G.rawLinearisedCellData(rstrt, "SOIL"); + + if (so.size() == G.numCells()) { + // Use "SOIL" directly if available. + return so; + } + + // SOIL vector not provided. Compute from SWAT and/or SGAS. + + so.assign(G.numCells(), 1.0); + + auto adjust_So_for_other_phase = + [&so](const std::vector& s) + { + std::transform(std::begin(so), std::end(so), + std::begin(s) , + std::begin(so), std::minus()); + }; + + if (sg.size() == G.numCells()) { + adjust_So_for_other_phase(sg); + } + + if (sw.size() == G.numCells()) { + adjust_So_for_other_phase(sw); + } + + return so; + } +} // Anonymous + +class RegionMapping +{ +private: + using Map = ::Opm::AssembledConnections; + using NeighIt = Map::Neighbours::const_iterator; + using MapIdx = Map::Offset; + +public: + RegionMapping(const std::size_t numCells, + const std::vector& regIdx); + + using NeighRng = ::Opm::SimpleIteratorRange; + + MapIdx numRegions() const + { + return this->map_.numRows(); + } + + NeighRng cells(const MapIdx regId) const + { + assert (regId < this->numRegions()); + + const auto& start = this->map_.startPointers(); + const auto& neigh = this->map_.neighbourhood(); + + auto b = std::begin(neigh) + start[regId + 0]; + auto e = std::begin(neigh) + start[regId + 1]; + + return { b, e }; + } + +private: + Opm::AssembledConnections map_; +}; + +RegionMapping::RegionMapping(const std::size_t numCells, + const std::vector& regIdx) +{ + if (regIdx.empty()) { + // No explicit region mapping. Put all active cells in single + // region (region ID 0). This is somewhat roundabout since class + // AssembledConnections does not have a direct way of expressing + // this case. + const auto nc = static_cast(numCells); + + for (auto c = 0*nc; c < nc; ++c) { + this->map_.addConnection(0, c); + } + + this->map_.compress(1); + } + else if (regIdx.size() != numCells) { + throw std::invalid_argument { + "Region Array Size Does Not " + "Match Number of Active Cells" + }; + } + else { + // Caller provided explicit region mapping for all active cells. + // Assume that the region IDs themselves are one-based indices + // (e.g., SATNUM) and adjust accordingly. + + auto maxReg = -1; + auto c = 0; + + for (const auto& regId : regIdx) { + const auto regId_0based = regId - 1; + + this->map_.addConnection(regId_0based, c++); + + if (regId_0based > maxReg) { maxReg = regId_0based; } + } + + this->map_.compress(maxReg + 1); + } +} + +// ===================================================================== + +namespace Relperm { + namespace Gas { + namespace Details { + Opm::ECLPropTableRawData + tableData(const std::vector& tabdims, + const std::vector& tab); + } + + class KrFunction + { + public: + KrFunction(const std::vector& tabdims, + const std::vector& tab) + : func_(Details::tableData(tabdims, tab)) + {} + + std::vector sgco() const + { + return this->func_.connateSat(); + } + + std::vector sgcr() const + { + return this->func_.criticalSat(this->krcol()); + } + + std::vector sgmax() const + { + return this->func_.maximumSat(); + } + + std::vector + krg(const std::size_t regID, + const std::vector& sg) const + { + const auto t = this->table(regID); + const auto c = this->krcol(); + + return this->func_.interpolate(t, c, sg); + } + + private: + ::Opm::SatFuncInterpolant func_; + + Opm::SatFuncInterpolant::InTable + table(const Opm::ECLPropTableRawData::SizeType regID) const + { + return ::Opm::SatFuncInterpolant::InTable{regID}; + } + + Opm::SatFuncInterpolant::ResultColumn krcol() const + { + return ::Opm::SatFuncInterpolant::ResultColumn{0}; + } + }; + } // namespace Gas + + namespace Oil { + namespace Details { + Opm::ECLPropTableRawData + tableData(const std::vector& tabdims, + const bool isTwoP, + const std::vector& tab); + } // namespace Details + + class KrFunction + { + public: + KrFunction(const std::vector& tabdims, + const bool isTwoP, + const std::vector& tab) + : func_(Details::tableData(tabdims, isTwoP, tab)) + , twop_(isTwoP) + {} + + virtual ~KrFunction() {} + + std::vector soco() const + { + return this->func_.connateSat(); + } + + std::vector sogcr() const + { + return this->func_.criticalSat(this->gas_column()); + } + + std::vector sowcr() const + { + return this->func_.criticalSat(this->wat_column()); + } + + std::vector somax() const + { + return this->func_.maximumSat(); + } + + struct SGas { + std::vector data; + }; + + struct SWat { + std::vector data; + }; + + struct SOil { + std::vector data; + }; + + std::vector + kro(const std::size_t regID, + const SOil& so_g, + const SGas& sg, + const SOil& so_w, + const SWat& sw) const + { + return this->kroImpl(regID, so_g, sg, so_w, sw); + } + + virtual std::unique_ptr clone() const = 0; + + protected: + std::vector + krog(const std::size_t regID, + const std::vector& so) const + { + return this->eval(regID, this->gas_column(), so); + } + + std::vector + krow(const std::size_t regID, + const std::vector& so) const + { + return this->eval(regID, this->wat_column(), so); + } + + private: + using ResCol = ::Opm::SatFuncInterpolant::ResultColumn; + + Opm::SatFuncInterpolant func_; + bool twop_; + + Opm::SatFuncInterpolant::InTable + table(const Opm::ECLPropTableRawData::SizeType regID) const + { + return ::Opm::SatFuncInterpolant::InTable{regID}; + } + + ResCol gas_column() const + { + // Table format: + // + // So kro in two-phase runs + // So krow krog in three-phase runs + // + // Therefore kro(so) in the O/G subsystem is result column + // (dependent variable) zero in two-phase runs and result + // column 1 in three-phase runs. + const auto colID = + static_cast + (this->twop_ ? 0 : 1); + + return { colID }; + } + + ResCol wat_column() const + { + // Note: kro(so) in the O/W subsystem is dependent variable + // (result column) zero in two-phase runs and three-phase + // runs. + return ResCol{0}; + } + + std::vector + eval(const std::size_t regID, + const ResCol c, + const std::vector& so) const + { + return this->func_.interpolate(this->table(regID), c, so); + } + + virtual std::vector + kroImpl(const std::size_t regID, + const SOil& so_g, + const SGas& sg, + const SOil& so_w, + const SWat& sw) const = 0; + }; + + class TwoPhase : public KrFunction + { + public: + enum class SubSys { OilGas, OilWater }; + + TwoPhase(const SubSys subsys, + const std::vector& tabdims, + const std::vector& tab) + : KrFunction(tabdims, true, tab) + , subsys_ (subsys) + {} + + virtual std::unique_ptr clone() const override + { + return std::unique_ptr(new TwoPhase(*this)); + } + + private: + SubSys subsys_; + + virtual std::vector + kroImpl(const std::size_t regID, + const SOil& so_g, + const SGas& /* sg */, + const SOil& so_w, + const SWat& /* sw */) const override + { + switch (this->subsys_) { + case SubSys::OilGas: + return this->krog(regID, so_g.data); + + case SubSys::OilWater: + return this->krow(regID, so_w.data); + } + + return {}; + } + }; + + class ECLStdThreePhase : public KrFunction + { + public: + ECLStdThreePhase(const std::vector& tabdims, + const std::vector& tab, + std::vector swco) + : KrFunction(tabdims, false, tab) + , swco_ (std::move(swco)) + {} + + virtual std::unique_ptr clone() const override + { + return std::unique_ptr(new ECLStdThreePhase(*this)); + } + + private: + std::vector swco_; + + virtual std::vector + kroImpl(const std::size_t regID, + const SOil& so_g, + const SGas& sg, + const SOil& so_w, + const SWat& sw) const override + { + const auto kr_og = this->krog(regID, so_g.data); + const auto kr_ow = this->krow(regID, so_w.data); + + const auto swco = this->swco_[regID]; + + auto kro = std::vector{}; + kro.reserve(sw.data.size()); + + for (auto n = sw.data.size(), i = 0*n; i < n; ++i) { + const auto den = sg.data[i] + sw.data[i] - swco; + + const auto xg = sg.data[i] / den; + const auto xw = (sw.data[i] + swco) / den; + + kro.push_back(xg*kr_og[i] + xw*kr_ow[i]); + } + + return kro; + } + }; + } // namespace Oil + + namespace Water { + namespace Details { + Opm::ECLPropTableRawData + tableData(const std::vector& tabdims, + const std::vector& tab); + } // namespace Details + + class KrFunction + { + public: + KrFunction(const std::vector& tabdims, + const std::vector& tab) + : func_(Details::tableData(tabdims, tab)) + {} + + std::vector swco() const + { + return this->func_.connateSat(); + } + + std::vector swcr() const + { + return this->func_.criticalSat(this->krcol()); + } + + std::vector swmax() const + { + return this->func_.maximumSat(); + } + + std::vector + krw(const std::size_t regID, + const std::vector& sw) const + { + const auto t = this->table(regID); + const auto c = this->krcol(); + + return this->func_.interpolate(t, c, sw); + } + + private: + ::Opm::SatFuncInterpolant func_; + + Opm::SatFuncInterpolant::InTable + table(const Opm::ECLPropTableRawData::SizeType regID) const + { + return ::Opm::SatFuncInterpolant::InTable{regID}; + } + + Opm::SatFuncInterpolant::ResultColumn krcol() const + { + return ::Opm::SatFuncInterpolant::ResultColumn{0}; + } + }; + } // namespace Water +} // namespace Relperm + +Opm::ECLPropTableRawData +Relperm::Gas::Details::tableData(const std::vector& tabdims, + const std::vector& tab) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.numCols = 3; // Sg, Krg, Pcgo + t.numRows = tabdims[ TABDIMS_NSSGFN_ITEM ]; + t.numTables = tabdims[ TABDIMS_NTSGFN_ITEM ]; + + const auto nTabElems = t.numRows * t.numTables * t.numCols; + + // Subtract one to account for offset being one-based index. + const auto start = tabdims[ TABDIMS_IBSGFN_OFFSET_ITEM ] - 1; + + t.data.assign(&tab[start], &tab[start] + nTabElems); + + return t; +} + +Opm::ECLPropTableRawData +Relperm::Oil::Details::tableData(const std::vector& tabdims, + const bool isTwoP, + const std::vector& tab) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.numCols = isTwoP + ? 2 // So, Kro + : 3; // So, Krow, Krog + + t.numRows = tabdims[ TABDIMS_NSSOFN_ITEM ]; + t.numTables = tabdims[ TABDIMS_NTSOFN_ITEM ]; + + const auto nTabElems = t.numRows * t.numTables * t.numCols; + + // Subtract one to account for offset being one-based index. + const auto start = tabdims[ TABDIMS_IBSOFN_OFFSET_ITEM ] - 1; + + t.data.assign(&tab[start], &tab[start] + nTabElems); + + return t; +} + +Opm::ECLPropTableRawData +Relperm::Water::Details::tableData(const std::vector& tabdims, + const std::vector& tab) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.numCols = 3; // Sw, Krw, Pcow + t.numRows = tabdims[ TABDIMS_NSSWFN_ITEM ]; + t.numTables = tabdims[ TABDIMS_NTSWFN_ITEM ]; + + const auto nTabElems = t.numRows * t.numTables * t.numCols; + + // Subtract one to account for offset being one-based index. + const auto start = tabdims[ TABDIMS_IBSWFN_OFFSET_ITEM ] - 1; + + t.data.assign(&tab[start], &tab[start] + nTabElems); + + return t; +} + +// ===================================================================== + +class Opm::ECLSaturationFunc::Impl +{ +public: + Impl(const ECLGraph& G, + const ECLInitFileData& init); + + Impl(Impl&& rhs); + Impl(const Impl& rhs); + + void init(const ECLGraph& G, + const ECLInitFileData& init, + const bool useEPS); + + std::vector + relperm(const ECLGraph& G, + const ECLRestartData& rstrt, + const ECLPhaseIndex p) const; + +private: + class EPSEvaluator + { + public: + using RawTEP = ::Opm::SatFunc::CreateEPS::RawTableEndPoints; + using FuncCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; + + struct ActPh + { + ActPh(const unsigned int iphs) + : oil((iphs & (1u << 0)) != 0) + , gas((iphs & (1u << 2)) != 0) + , wat((iphs & (1u << 1)) != 0) + {} + + bool oil; + bool gas; + bool wat; + }; + + void define(const ECLGraph& G, + const ECLInitFileData& init, + const RawTEP& ep, + const bool use3PtScaling, + const FuncCat curve, + const ActPh& active) + { + auto opt = Create::EPSOptions{}; + opt.use3PtScaling = use3PtScaling; + opt.curve = curve; + + if (active.oil) { + this->create_oil_eps(G, init, ep, active, opt); + } + + if (active.gas) { + this->create_gas_eps(G, init, ep, opt); + } + + if (active.wat) { + this->create_wat_eps(G, init, ep, opt); + } + } + + void scaleOG(const RegionMapping& rmap, + std::vector& so) const + { + this->scale(this->oil_in_og_, rmap, so); + } + + void scaleOW(const RegionMapping& rmap, + std::vector& so) const + { + this->scale(this->oil_in_ow_, rmap, so); + } + + void scaleGas(const RegionMapping& rmap, + std::vector& sg) const + { + this->scale(this->gas_, rmap, sg); + } + + void scaleWat(const RegionMapping& rmap, + std::vector& sw) const + { + this->scale(this->wat_, rmap, sw); + } + + private: + using Create = ::Opm::SatFunc::CreateEPS; + using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; + using PhIdx = ::Opm::ECLPhaseIndex; + + using EPSInterface = ::Opm::SatFunc::EPSEvalInterface; + using EPSEndPts = ::Opm::SatFunc::EPSEvalInterface::TableEndPoints; + using EPSEndPtVec = std::vector; + + using EPSPtr = std::unique_ptr; + using EndPtsPtr = std::unique_ptr; + + struct EPS { + EPS() {} + + EPS(const EPS& rhs) + { + if (rhs.scaling) { + this->scaling = rhs.scaling->clone(); + } + + if (rhs.tep) { + this->tep = EndPtsPtr(new EPSEndPtVec(*rhs.tep)); + } + } + + EPS(EPS&& rhs) + : scaling(std::move(rhs.scaling)) + , tep (std::move(rhs.tep)) + {} + + EPS& operator=(const EPS& rhs) + { + if (rhs.scaling) { + this->scaling = rhs.scaling->clone(); + } + + if (rhs.tep) { + this->tep = EndPtsPtr(new EPSEndPtVec(*rhs.tep)); + } + + return *this; + } + + EPS& operator=(EPS&& rhs) + { + this->scaling = std::move(rhs.scaling); + this->tep = std::move(rhs.tep); + + return *this; + } + + EPSPtr scaling; + EndPtsPtr tep; + }; + + EPS oil_in_og_; + EPS oil_in_ow_; + EPS gas_; + EPS wat_; + + void scale(const EPS& eps, + const RegionMapping& rmap, + std::vector& s) const + { + if (! eps.scaling) { + // No end-point scaling defined for this curve. Return + // unchanged. + return; + } + + if (! eps.tep) { + throw std::logic_error { + "Cannot Perform EPS in Absence of " + "Table End-Point Data" + }; + } + + if (rmap.numRegions() == 1) { + this->scaleSingleRegion(eps, s); + } + else { + this->scaleMultiRegion(eps, rmap, s); + } + } + + void scaleSingleRegion(const EPS& eps, std::vector& s) const + { + assert (eps.tep->size() == 1); + + using Assoc = EPSInterface::SaturationAssoc; + using CellID = decltype(std::declval().cell); + + auto sp = EPSInterface::SaturationPoints{}; + + sp.reserve(s.size()); + + auto cell = static_cast(0); + for (const auto& si : s) { + sp.push_back(Assoc{ cell++, si }); + } + + s = eps.scaling->eval((*eps.tep)[0], sp); + } + + void scaleMultiRegion(const EPS& eps, + const RegionMapping& rmap, + std::vector& s) const + { + const auto nreg = rmap.numRegions(); + + assert (eps.tep->size() == nreg); + + using Assoc = EPSInterface::SaturationAssoc; + using CellID = decltype(std::declval().cell); + + for (auto reg = 0*nreg; reg < nreg; ++reg) { + auto sp = EPSInterface::SaturationPoints{}; + + for (const auto& cell : rmap.cells(reg)) { + sp.push_back(Assoc{CellID(cell), s[cell]}); + } + + const auto& sr = + eps.scaling->eval((*eps.tep)[reg], sp); + + auto i = static_cast(0); + for (const auto& cell : rmap.cells(reg)) { + s[cell] = sr[i++]; + } + } + } + + void create_oil_eps(const ECLGraph& G, + const ECLInitFileData& init, + const RawTEP& ep, + const ActPh& active, + Create::EPSOptions& opt) + { + opt.thisPh = PhIdx::Liquid; + + if (active.gas) { + opt.subSys = SSys::OilGas; + + this->oil_in_og_.scaling = + Create::fromECLOutput(G, init, opt); + + this->oil_in_og_.tep = this->endPoints(ep, opt); + } + + if (active.wat) { + opt.subSys = SSys::OilWater; + + this->oil_in_ow_.scaling = + Create::fromECLOutput(G, init, opt); + + this->oil_in_ow_.tep = this->endPoints(ep, opt); + } + } + + void create_gas_eps(const ECLGraph& G, + const ECLInitFileData& init, + const RawTEP& ep, + Create::EPSOptions& opt) + { + opt.thisPh = PhIdx::Vapour; + opt.subSys = SSys::OilGas; + + this->gas_.scaling = + Create::fromECLOutput(G, init, opt); + + this->gas_.tep = this->endPoints(ep, opt); + } + + void create_wat_eps(const ECLGraph& G, + const ECLInitFileData& init, + const RawTEP& ep, + Create::EPSOptions& opt) + { + opt.thisPh = PhIdx::Aqua; + opt.subSys = SSys::OilWater; + + this->wat_.scaling = + Create::fromECLOutput(G, init, opt); + + this->wat_.tep = this->endPoints(ep, opt); + } + + EndPtsPtr + endPoints(const RawTEP& ep, const Create::EPSOptions& opt) + { + return EndPtsPtr { + new EPSEndPtVec(Create::unscaledEndPoints(ep, opt)) + }; + } + }; + + using RegionID = + decltype(std::declval().numRegions()); + + RegionMapping rmap_; + + std::unique_ptr oil_; + std::unique_ptr gas_; + std::unique_ptr wat_; + + std::unique_ptr eps_; + + void initRelPermInterp(const EPSEvaluator::ActPh& active, + const ECLInitFileData& init); + + void initEPS(const EPSEvaluator::ActPh& active, + const bool use3PtScaling, + const ECLGraph& G, + const ECLInitFileData& init); + + std::vector + kro(const ECLGraph& G, + const ECLRestartData& rstrt) const; + + std::vector + krg(const ECLGraph& G, + const ECLRestartData& rstrt) const; + + std::vector + krw(const ECLGraph& G, + const ECLRestartData& rstrt) const; + + EPSEvaluator::RawTEP + extractRawTableEndPoints(const EPSEvaluator::ActPh& active) const; + + template + std::vector + gatherRegionSubset(const RegionID reg, + const std::vector& x) const + { + auto y = std::vector{}; + + if (x.empty()) { + return y; + } + + for (const auto& cell : this->rmap_.cells(reg)) { + y.push_back(x[cell]); + } + + return y; + } + + template + void scatterRegionResults(const RegionID reg, + const std::vector& x_reg, + std::vector& x) const + { + auto i = static_cast(0); + + for (const auto& cell : this->rmap_.cells(reg)) { + x[cell] = x_reg[i++]; + } + } + + template + void regionLoop(RegionOperation&& regOp) const + { + for (auto nreg = this->rmap_.numRegions(), + reg = 0*nreg; reg < nreg; ++reg) + { + regOp(reg); + } + } +}; + +Opm::ECLSaturationFunc::Impl::Impl(const ECLGraph& G, + const ECLInitFileData& init) + : rmap_(G.numCells(), G.rawLinearisedCellData(init, "SATNUM")) +{ +} + +Opm::ECLSaturationFunc::Impl::Impl(Impl&& rhs) + : rmap_(std::move(rhs.rmap_)) + , oil_ (std::move(rhs.oil_ )) + , gas_ (std::move(rhs.gas_ )) + , wat_ (std::move(rhs.wat_ )) +{} + +// --------------------------------------------------------------------- + +void +Opm::ECLSaturationFunc::Impl::init(const ECLGraph& G, + const ECLInitFileData& init, + const bool useEPS) +{ + // Extract INTEHEAD from main grid + const auto& ih = init.keywordData(INTEHEAD_KW); + const auto iphs = static_cast(ih[INTEHEAD_PHASE_INDEX]); + + const auto active = EPSEvaluator::ActPh{iphs}; + + this->initRelPermInterp(active, init); + + if (useEPS) { + const auto& lh = init.keywordData(LOGIHEAD_KW); + + const auto haveEPS = static_cast( + lh[LOGIHEAD_ENDPOINT_SCALING_INDEX]); + + if (haveEPS) { + const auto use3PtScaling = static_cast( + lh[LOGIHEAD_ALT_ENDPOINT_SCALING_INDEX]); + + // Must be called *after* initRelPermInterp(). + this->initEPS(active, use3PtScaling, G, init); + } + } +} + +void +Opm::ECLSaturationFunc:: +Impl::initRelPermInterp(const EPSEvaluator::ActPh& active, + const ECLInitFileData& init) +{ + const auto isThreePh = + active.oil && active.gas && active.wat; + + // Extract tabular data from main grid + const auto& tabdims = init.keywordData("TABDIMS"); + const auto& tab = init.keywordData("TAB"); + + if (active.gas) { + this->gas_.reset(new Relperm::Gas::KrFunction(tabdims, tab)); + } + + if (active.wat) { + this->wat_.reset(new Relperm::Water::KrFunction(tabdims, tab)); + } + + if (active.oil) { + if (! isThreePh) { + using KrModel = Relperm::Oil::TwoPhase; + + if (active.gas) { + const auto subsys = KrModel::SubSys::OilGas; + + this->oil_.reset(new KrModel(subsys, tabdims, tab)); + } + else if (active.wat) { + const auto subsys = KrModel::SubSys::OilWater; + + this->oil_.reset(new KrModel(subsys, tabdims, tab)); + } + else { + throw std::invalid_argument { + "Single-Phase Oil System Not Supported" + }; + } + } + + if (isThreePh) { + using KrModel = Relperm::Oil::ECLStdThreePhase; + + this->oil_.reset(new KrModel(tabdims, tab, this->wat_->swco())); + } + } +} + +void Opm::ECLSaturationFunc:: +Impl::initEPS(const EPSEvaluator::ActPh& active, + const bool use3PtScaling, + const ECLGraph& G, + const ECLInitFileData& init) +{ + const auto ep = this->extractRawTableEndPoints(active); + + const auto curve = ::Opm::SatFunc::CreateEPS:: + FunctionCategory::Relperm; + + this->eps_.reset(new EPSEvaluator()); + + this->eps_->define(G, init, ep, use3PtScaling, curve, active); +} + +// ##################################################################### + +Opm::ECLSaturationFunc::Impl::Impl(const Impl& rhs) + : rmap_(rhs.rmap_) +{ + if (rhs.oil_) { + // Polymorphic object must use clone(). + this->oil_ = rhs.oil_->clone(); + } + + if (rhs.gas_) { + this->gas_.reset(new Relperm::Gas::KrFunction(*rhs.gas_)); + } + + if (rhs.wat_) { + this->wat_.reset(new Relperm::Water::KrFunction(*rhs.wat_)); + } +} + +// ##################################################################### + +std::vector +Opm::ECLSaturationFunc::Impl:: +relperm(const ECLGraph& G, + const ECLRestartData& rstrt, + const ECLPhaseIndex p) const +{ + switch (p) { + case ECLPhaseIndex::Aqua: + return this->krw(G, rstrt); + + case ECLPhaseIndex::Liquid: + return this->kro(G, rstrt); + + case ECLPhaseIndex::Vapour: + return this->krg(G, rstrt); + } + + return {}; +} + +std::vector +Opm::ECLSaturationFunc::Impl:: +kro(const ECLGraph& G, + const ECLRestartData& rstrt) const +{ + auto kr = std::vector{}; + + if (! this->oil_) { + return kr; + } + + const auto& sg = G.rawLinearisedCellData(rstrt, "SGAS"); + const auto& sw = G.rawLinearisedCellData(rstrt, "SWAT"); + + auto so_g = oil_saturation(sg, sw, G, rstrt); + auto so_w = so_g; + + if (this->eps_) { + // Independent scaling of So in O/G and O/W sub-systems of an O/G/W + // run. Performs duplicate work in a two-phase case. Need to take + // action if this becomes a bottleneck. + this->eps_->scaleOG(this->rmap_, so_g); + this->eps_->scaleOW(this->rmap_, so_w); + } + + // Allocate result. Member function scatterRegionResult() depends on + // having an allocated result vector into which to write the values from + // a single region. + kr.resize(so_g.size(), 0.0); + + // Compute relative permeability per region. + this->regionLoop([this, &so_g, &so_w, &sg, &sw, &kr] + (const RegionID reg) + { + const auto So_g = Relperm::Oil::KrFunction::SOil { + this->gatherRegionSubset(reg, so_g) + }; + + const auto So_w = Relperm::Oil::KrFunction::SOil { + this->gatherRegionSubset(reg, so_w) + }; + + const auto Sg = Relperm::Oil::KrFunction::SGas { + // Empty in case of Oil/Water system + this->gatherRegionSubset(reg, sg) + }; + + const auto Sw = Relperm::Oil::KrFunction::SWat { + // Empty in case of Oil/Gas system + this->gatherRegionSubset(reg, sw) + }; + + const auto& kro_reg = + this->oil_->kro(reg, So_g, Sg, So_w, Sw); + + this->scatterRegionResults(reg, kro_reg, kr); + }); + + return kr; +} + +std::vector +Opm::ECLSaturationFunc::Impl:: +krg(const ECLGraph& G, + const ECLRestartData& rstrt) const +{ + auto kr = std::vector{}; + + if (! this->gas_) { + return kr; + } + + auto sg = G.rawLinearisedCellData(rstrt, "SGAS"); + + if (this->eps_) { + this->eps_->scaleGas(this->rmap_, sg); + } + + // Allocate result. Member function scatterRegionResult() depends on + // having an allocated result vector into which to write the values from + // a single region. + kr.resize(sg.size(), 0.0); + + // Compute relative permeability per region. + this->regionLoop([this, &sg, &kr](const RegionID reg) + { + const auto sg_reg = this->gatherRegionSubset(reg, sg); + + const auto krg_reg = + this->gas_->krg(reg, sg_reg); + + this->scatterRegionResults(reg, krg_reg, kr); + }); + + return kr; +} + +std::vector +Opm::ECLSaturationFunc::Impl:: +krw(const ECLGraph& G, + const ECLRestartData& rstrt) const +{ + auto kr = std::vector{}; + + if (! this->wat_) { + return kr; + } + + auto sw = G.rawLinearisedCellData(rstrt, "SWAT"); + + if (this->eps_) { + this->eps_->scaleWat(this->rmap_, sw); + } + + // Allocate result. Member function scatterRegionResult() depends on + // having an allocated result vector into which to write the values from + // a single region. + kr.resize(sw.size(), 0.0); + + // Compute relative permeability per region. + this->regionLoop([this, &sw, &kr](const RegionID reg) + { + const auto sw_reg = this->gatherRegionSubset(reg, sw); + + const auto krw_reg = + this->wat_->krw(reg, sw_reg); + + this->scatterRegionResults(reg, krw_reg, kr); + }); + + return kr; +} + +Opm::ECLSaturationFunc::Impl::EPSEvaluator::RawTEP +Opm::ECLSaturationFunc::Impl:: +extractRawTableEndPoints(const EPSEvaluator::ActPh& active) const +{ + auto ep = EPSEvaluator::RawTEP{}; + + if (active.oil) { + ep.conn.oil = this->oil_->soco(); + ep.crit.oil_in_gas = this->oil_->sogcr(); + ep.crit.oil_in_water = this->oil_->sowcr(); + ep.smax.oil = this->oil_->somax(); + } + + if (active.gas) { + ep.conn.gas = this->gas_->sgco(); + ep.crit.gas = this->gas_->sgcr(); + ep.smax.gas = this->gas_->sgmax(); + } + + if (active.wat) { + ep.conn.water = this->wat_->swco(); + ep.crit.water = this->wat_->swcr(); + ep.smax.water = this->wat_->swmax(); + } + + return ep; +} + +// ===================================================================== + +Opm::ECLSaturationFunc:: +ECLSaturationFunc(const ECLGraph& G, + const ECLInitFileData& init, + const bool useEPS) + : pImpl_(new Impl(G, init)) +{ + this->pImpl_->init(G, init, useEPS); +} + +Opm::ECLSaturationFunc::~ECLSaturationFunc() +{} + +Opm::ECLSaturationFunc::ECLSaturationFunc(ECLSaturationFunc&& rhs) + : pImpl_(std::move(rhs.pImpl_)) +{} + +Opm::ECLSaturationFunc::ECLSaturationFunc(const ECLSaturationFunc& rhs) + : pImpl_(new Impl(*rhs.pImpl_)) +{} + +Opm::ECLSaturationFunc& +Opm::ECLSaturationFunc::operator=(ECLSaturationFunc&& rhs) +{ + this->pImpl_ = std::move(rhs.pImpl_); + + return *this; +} + +Opm::ECLSaturationFunc& +Opm::ECLSaturationFunc::operator=(const ECLSaturationFunc& rhs) +{ + this->pImpl_.reset(new Impl(*rhs.pImpl_)); + + return *this; +} + +std::vector +Opm::ECLSaturationFunc:: +relperm(const ECLGraph& G, + const ECLRestartData& rstrt, + const ECLPhaseIndex p) const +{ + return this->pImpl_->relperm(G, rstrt, p); +} diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp new file mode 100644 index 0000000000..658b4048cd --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp @@ -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 . +*/ + +#ifndef OPM_ECLSATURATIONFUNC_HEADER_INCLUDED +#define OPM_ECLSATURATIONFUNC_HEADER_INCLUDED + +#include + +#include +#include + +/// \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 + 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 pImpl_; + }; + +} // namespace Opm + +#endif // OPM_ECLSATURATIONFUNC_HEADER_INCLUDED diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp index 1811f8a3ed..4fdd583e62 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runLinearisedCellDataTest.cpp @@ -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)) { diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp index 766a11048c..e42ee4e862 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/runTransTest.cpp @@ -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); diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp new file mode 100644 index 0000000000..1aff18e58c --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp @@ -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 . +*/ + +#if HAVE_CONFIG_H +#include +#endif // HAVE_CONFIG_H + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE + +#define BOOST_TEST_MODULE TEST_ECLENDPOINTSCALING + +#include +#include +#include + +#include + +#include +#include +#include + +namespace { + template + 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& 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{ 0.0 }; + const auto smax = std::vector{ 1.0 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.2 }; + const auto smax = std::vector{ 1.0 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.0, 0.0, 1.0 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.0 }; + const auto smax = std::vector{ 0.8 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.0, 0.0, 1.0 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.2 }; + const auto smax = std::vector{ 0.8 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.0, 0.0, 1.0 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.2 }; + const auto smax = std::vector{ 0.8 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.2, 0.0, 0.8 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.0 }; + const auto smax = std::vector{ 1.0 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.2, 0.0, 0.8 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.2 }; + const auto smax = std::vector{ 1.0 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.2, 0.0, 0.8 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.5 }; + const auto smax = std::vector{ 0.7 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.2, 0.0, 0.8 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.0 }; + const auto sdisp = std::vector{ 0.2 }; + const auto smax = std::vector{ 1.0 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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{ 0.1 }; + const auto sdisp = std::vector{ 0.4 }; + const auto smax = std::vector{ 1.0 }; + + const auto tep = SF::EPSEvalInterface:: + TableEndPoints { 0.0, 0.2, 1.0 }; + + const auto s = std::vector { + 0.0, + 0.2, + 0.4, + 0.6, + 0.8, + 1.0, + }; + + const auto sp = associate(s); + const auto expect = std::vector { + 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 () diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclproptable.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclproptable.cpp new file mode 100644 index 0000000000..f6e4f4f92c --- /dev/null +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclproptable.cpp @@ -0,0 +1,2417 @@ +/* + 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 . +*/ + +#if HAVE_CONFIG_H +#include +#endif // HAVE_CONFIG_H + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE + +#define BOOST_TEST_MODULE TEST_ECLPROPTABLE + +#include +#include +#include + +#include + +#include +#include + +namespace { + template + 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::ECLPropTableRawData + toRawTableFormat(Opm::ECLPropTableRawData t) + { + // Note: Raw table format is nTab*nRows consecutive values for one + // column followed by nTab*nRows consecutive values for the next + // column &c. + + const auto d = t.data; + const auto rTabStride = t.numRows * t.numCols; + const auto wColStride = t.numRows * t.numTables; + + for (auto c = 0*t.numCols; c < t.numCols; ++c) { + const auto wStart = c * wColStride; + + for (auto k = 0*t.numTables; k < t.numTables; ++k) { + const auto rStart = k * rTabStride; + const auto wOff = k * t.numRows; + + for (auto i = 0*t.numRows; i < t.numRows; ++i) { + t.data[wStart + wOff + i] = + d [rStart + i*t.numCols + c]; + } + } + } + + return t; + } +} // Namespace Anonymous + +// ===================================================================== +// Invalid tables (error handling/input validation) +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (InvalidTables) + +BOOST_AUTO_TEST_CASE (EmptyTable) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + }; + + t.numRows = 0; + t.numCols = 3; + t.numTables = 1; + + BOOST_CHECK_THROW(Opm::SatFuncInterpolant(toRawTableFormat(t)), + std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE (SingleNode) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.3 , 0.1 , 0.0, + }; + + t.numRows = 1; + t.numCols = 3; + t.numTables = 1; + + BOOST_CHECK_THROW(Opm::SatFuncInterpolant(toRawTableFormat(t)), + std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE (NoResultColumns) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s + 0.2, + 0.3, + 0.7, + 0.8, + }; + + t.numRows = 4; + t.numCols = 1; + t.numTables = 1; + + BOOST_CHECK_THROW(Opm::SatFuncInterpolant(toRawTableFormat(t)), + std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE (EmptyTableLargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , kr , pc + -1.0e+20, -1.0e+20, 0.0, + -1.0e+20, -1.0e+20, 0.0, + -1.0e+20, -1.0e+20, 0.0, + -1.0e+20, -1.0e+20, 0.0, + 1.0e+20 , 1.0e+20, 0.0, + 1.0e+20 , 1.0e+20, 0.0, + 1.0e+20 , 1.0e+20, 0.0, + 1.0e+20 , 1.0e+20, 0.0, + }; + + t.numRows = 8; + t.numCols = 3; + t.numTables = 1; + + BOOST_CHECK_THROW(Opm::SatFuncInterpolant(toRawTableFormat(t)), + std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE (SingleNodeLargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , kr , pc + -1.0e+20, -1.0e+20, 0.0, + -1.0e+20, -1.0e+20, 0.0, + -1.0e+20, -1.0e+20, 0.0, + -1.0e+20, -1.0e+20, 0.0, + 0.3 , 0.1 , 0.0, + 1.0e+20 , 1.0e+20, 0.0, + 1.0e+20 , 1.0e+20, 0.0, + 1.0e+20 , 1.0e+20, 0.0, + 1.0e+20 , 1.0e+20, 0.0, + }; + + t.numRows = 9; + t.numCols = 3; + t.numTables = 1; + + BOOST_CHECK_THROW(Opm::SatFuncInterpolant(toRawTableFormat(t)), + std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE (NoResultColumnsLargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s + -1.0e+20, + -1.0e+20, + -1.0e+20, + -1.0e+20, + 0.2, + 0.3, + 0.7, + 0.8, + 1.0e+20, + 1.0e+20, + 1.0e+20, + 1.0e+20, + }; + + t.numRows = 12; + t.numCols = 1; + t.numTables = 1; + + BOOST_CHECK_THROW(Opm::SatFuncInterpolant(toRawTableFormat(t)), + std::invalid_argument); +} + +BOOST_AUTO_TEST_SUITE_END () + +// ===================================================================== +// Single table (i.e., a single region). +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (InterpolationSingleTable) + +BOOST_AUTO_TEST_CASE (AtNodes) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 1; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ 0.8, 0.3, 0.3, 0.2 }; + const auto kr_expect = std::vector{ 0.5, 0.1, 0.1, 0.0 }; + const auto pc_expect = std::vector{ 0.0, 0.0, 0.0, 0.0 }; + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + // Check interpolation + + const auto kr = swfunc.interpolate(InTable{0}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{0}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); + + // Check error handling + + // Table ID out of range. + BOOST_CHECK_THROW(swfunc.interpolate(InTable{10}, ResultColumn{0}, s), + std::invalid_argument); + + // Result Column ID out of range. + BOOST_CHECK_THROW(swfunc.interpolate(InTable{0}, ResultColumn{2}, s), + std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE (AboveAndBelow) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 1; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ 0.80000001, 0.9, 0.199999999, 0.1 }; + const auto kr_expect = std::vector{ 0.5, 0.5, 0.0, 0.0 }; + const auto pc_expect = std::vector{ 0.0, 0.0, 0.0, 0.0 }; + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto kr = swfunc.interpolate(InTable{0}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{0}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); +} + +BOOST_AUTO_TEST_CASE (Interpolation) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 1; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ + 0.2000, + 0.2300, + 0.2600, + 0.2900, + 0.3200, + 0.3500, + 0.3800, + 0.4100, + 0.4400, + 0.4700, + 0.5000, + 0.5300, + 0.5600, + 0.5900, + 0.6200, + 0.6500, + 0.6800, + 0.7100, + 0.7400, + 0.7700, + 0.8000, + }; + + const auto kr_expect = std::vector{ + 0, + 0.0300, + 0.0600, + 0.0900, + 0.1160, + 0.1400, + 0.1640, + 0.1880, + 0.2120, + 0.2360, + 0.2600, + 0.2840, + 0.3080, + 0.3320, + 0.3560, + 0.3800, + 0.4040, + 0.4280, + 0.4520, + 0.4760, + 0.5000, + }; + + const auto pc_expect = std::vector(s.size(), 0.0); + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto kr = swfunc.interpolate(InTable{0}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{0}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); +} + +BOOST_AUTO_TEST_CASE (InterpolationLargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + // 1e20 is a sentinel value that counts as row "ignored". + t.data = std::vector{ + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.8 , 0.5 , 0.0, // 8 + 1.0e20 , 1.0e+100 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + }; + + t.numRows = 15; + t.numCols = 3; + t.numTables = 1; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ + 0.0000, + 0.1000, + 0.1500, + 0.1900, + 0.2000, + 0.2300, + 0.2600, + 0.2900, + 0.3200, + 0.3500, + 0.3800, + 0.4100, + 0.4400, + 0.4700, + 0.5000, + 0.5300, + 0.5600, + 0.5900, + 0.6200, + 0.6500, + 0.6800, + 0.7100, + 0.7400, + 0.7700, + 0.8000, + 0.8100, + 0.8500, + 0.9000, + 1.0000, + }; + + const auto kr_expect = std::vector{ + 0, + 0, + 0, + 0, + 0, + 0.0300, + 0.0600, + 0.0900, + 0.1160, + 0.1400, + 0.1640, + 0.1880, + 0.2120, + 0.2360, + 0.2600, + 0.2840, + 0.3080, + 0.3320, + 0.3560, + 0.3800, + 0.4040, + 0.4280, + 0.4520, + 0.4760, + 0.5000, + 0.5000, + 0.5000, + 0.5000, + 0.5000, + }; + + const auto pc_expect = std::vector(s.size(), 0.0); + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto kr = swfunc.interpolate(InTable{0}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{0}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); +} + +BOOST_AUTO_TEST_SUITE_END () + +// ===================================================================== +// Multiple tables (i.e., multiple regions). +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (InterpolationFourTables) + +BOOST_AUTO_TEST_CASE (AtNodes) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // Table 0 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 1 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 2 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 3 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 4; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ 0.8, 0.3, 0.3, 0.2 }; + const auto kr_expect = std::vector{ 0.5, 0.1, 0.1, 0.0 }; + const auto pc_expect = std::vector{ 0.0, 0.0, 0.0, 0.0 }; + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + for (auto ti = 0*t.numTables; ti < t.numTables; ++ti) { + const auto kr = swfunc.interpolate(InTable{ti}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{ti}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); + + // Check error handling + + // Table ID out of range. + BOOST_CHECK_THROW(swfunc.interpolate(InTable{10}, ResultColumn{0}, s), + std::invalid_argument); + + // Result Column ID out of range. + BOOST_CHECK_THROW(swfunc.interpolate(InTable{0}, ResultColumn{2}, s), + std::invalid_argument); + } +} + +BOOST_AUTO_TEST_CASE (AboveAndBelow) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // Table 0 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 1 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 2 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 3 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 4; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ 0.80000001, 0.9, 0.199999999, 0.1 }; + const auto kr_expect = std::vector{ 0.5, 0.5, 0.0, 0.0 }; + const auto pc_expect = std::vector{ 0.0, 0.0, 0.0, 0.0 }; + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + for (auto ti = 0*t.numTables; ti < t.numTables; ++ti) { + const auto kr = swfunc.interpolate(InTable{ti}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{ti}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); + } +} + +BOOST_AUTO_TEST_CASE (Interpolation) +{ + auto t = Opm::ECLPropTableRawData{}; + + // 1e20 is a sentinel value that counts as row "ignored". + t.data = std::vector{ + // Table 0 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 1 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 2 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // Table 3 + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 4; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ + 0.2000, + 0.2300, + 0.2600, + 0.2900, + 0.3200, + 0.3500, + 0.3800, + 0.4100, + 0.4400, + 0.4700, + 0.5000, + 0.5300, + 0.5600, + 0.5900, + 0.6200, + 0.6500, + 0.6800, + 0.7100, + 0.7400, + 0.7700, + 0.8000, + }; + + const auto kr_expect = std::vector{ + 0, + 0.0300, + 0.0600, + 0.0900, + 0.1160, + 0.1400, + 0.1640, + 0.1880, + 0.2120, + 0.2360, + 0.2600, + 0.2840, + 0.3080, + 0.3320, + 0.3560, + 0.3800, + 0.4040, + 0.4280, + 0.4520, + 0.4760, + 0.5000, + }; + + const auto pc_expect = std::vector(s.size(), 0.0); + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + for (auto ti = 0*t.numTables; ti < t.numTables; ++ti) { + const auto kr = swfunc.interpolate(InTable{ti}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{ti}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); + } +} + +BOOST_AUTO_TEST_CASE (InterpolationLargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + // 1e20 is a sentinel value that counts as row "ignored". + t.data = std::vector{ + // Table 0 + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.7 , 0.15 , 0.0, // 8 + 0.8 , 0.5 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + + // Table 1 + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.7 , 0.15 , 0.0, // 8 + 0.8 , 0.5 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + + // Table 2 + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.7 , 0.15 , 0.0, // 8 + 0.8 , 0.5 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + + // Table 3 + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.7 , 0.15 , 0.0, // 8 + 0.8 , 0.5 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + }; + + t.numRows = 15; + t.numCols = 3; + t.numTables = 4; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + const auto s = std::vector{ + 0.0000, + 0.1000, + 0.1500, + 0.1900, + 0.2000, + 0.2300, + 0.2600, + 0.2900, + 0.3200, + 0.3500, + 0.3800, + 0.4100, + 0.4400, + 0.4700, + 0.5000, + 0.5300, + 0.5600, + 0.5900, + 0.6200, + 0.6500, + 0.6800, + 0.7100, + 0.7400, + 0.7700, + 0.8000, + 0.8100, + 0.8500, + 0.9000, + 1.0000, + }; + + const auto kr_expect = std::vector{ + 0, + 0, + 0, + 0, + 0, + 3.0000e-02, + 6.0000e-02, + 9.0000e-02, + 1.0250e-01, + 1.0625e-01, + 1.1000e-01, + 1.1375e-01, + 1.1750e-01, + 1.2125e-01, + 1.2500e-01, + 1.2875e-01, + 1.3250e-01, + 1.3625e-01, + 1.4000e-01, + 1.4375e-01, + 1.4750e-01, + 1.8500e-01, + 2.9000e-01, + 3.9500e-01, + 5.0000e-01, + 5.0000e-01, + 5.0000e-01, + 5.0000e-01, + 5.0000e-01, + }; + + const auto pc_expect = std::vector(s.size(), 0.0); + + using InTable = Opm::SatFuncInterpolant::InTable; + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + for (auto ti = 0*t.numTables; ti < t.numTables; ++ti) { + const auto kr = swfunc.interpolate(InTable{ti}, ResultColumn{0}, s); + const auto pc = swfunc.interpolate(InTable{ti}, ResultColumn{1}, s); + + check_is_close(kr, kr_expect); + check_is_close(pc, pc_expect); + } +} + +BOOST_AUTO_TEST_SUITE_END () + +// ===================================================================== +// Single Table End-Points +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (SingleTableEndPoints) + +BOOST_AUTO_TEST_CASE (SWFN_CritIsConn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_expect = std::vector{ 0.2 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SWFN_CritIsConn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + // 1e20 is a sentinel value that counts as row "ignored". + t.data = std::vector{ + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.8 , 0.5 , 0.0, // 8 + 1.0e20 , 1.0e+100 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + }; + + t.numRows = 15; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_expect = std::vector{ 0.2 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SWFN) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 4; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_expect = std::vector{ 0.21 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SWFN_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + // 1e20 is a sentinel value that counts as row "ignored". + t.data = std::vector{ + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.0, // 7 + 0.3 , 0.1 , 0.0, // 8 + 0.8 , 0.5 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + 1.0e20 , 1.0e+100 , 0.0, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_expect = std::vector{ 0.21 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_CritIsConn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.2 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_CritIsConn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.5, // 7 + 0.8 , 0.5 , 0.8, // 8 + 1.0e20 , 1.0e+100, 1.0e+100, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + }; + + t.numRows = 15; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.2 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOGCR_is_Conn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.1, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + }; + + t.numRows = 4; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21 }; + const auto scrit_krog_expect = std::vector{ 0.2 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOGCR_is_Conn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.1, // 7 + 0.3 , 0.1 , 0.5, // 8 + 0.8 , 0.5 , 0.8, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21 }; + const auto scrit_krog_expect = std::vector{ 0.2 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOWCR_is_Conn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.1 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + }; + + t.numRows = 4; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.21 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOWCR_is_Conn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.1 , 0.0, // 7 + 0.3 , 0.1 , 0.5, // 8 + 0.8 , 0.5 , 0.8, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.21 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SCR_Not_Conn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.205, 0.0 , 0.0, + 0.21 , 0.0 , 0.1, + 0.25 , 0.1 , 0.2, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + }; + + t.numRows = 6; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21 }; + const auto scrit_krog_expect = std::vector{ 0.205 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SCR_Not_Conn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.2 , 0.0 , 0.0 , // 5 + 0.205 , 0.0 , 0.0 , // 6 + 0.21 , 0.0 , 0.1 , // 7 + 0.25 , 0.1 , 0.2 , // 8 + 0.3 , 0.1 , 0.5 , // 9 + 0.8 , 0.5 , 0.8 , // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 1; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21 }; + const auto scrit_krog_expect = std::vector{ 0.205 }; + const auto smax_expect = std::vector{ 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_SUITE_END () + +// ===================================================================== +// Table End-Points in Multiple Tables +// --------------------------------------------------------------------- + +BOOST_AUTO_TEST_SUITE (TableEndPointsMultiTable) + +BOOST_AUTO_TEST_CASE (SWFN_CritIsConn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto scrit_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto smax_expect = std::vector{ 0.8, 0.8, 0.8, 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SWFN_CritIsConn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + // 1e20 is a sentinel value that counts as row "ignored". + t.data = std::vector{ + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.8 , 0.5 , 0.0, // 8 + 1.0e20 , 1.0e+100 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.8 , 0.5 , 0.0, // 8 + 1.0e20 , 1.0e+100 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.8 , 0.5 , 0.0, // 8 + 1.0e20 , 1.0e+100 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.0, // 7 + 0.8 , 0.5 , 0.0, // 8 + 1.0e20 , 1.0e+100 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + }; + + t.numRows = 15; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto scrit_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto smax_expect = std::vector{ 0.8, 0.8, 0.8, 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SWFN) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + + // s, kr , pc + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.0, + 0.3 , 0.1 , 0.0, + 0.8 , 0.5 , 0.0, + }; + + t.numRows = 4; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.2 , 0.2 , 0.2 }; + const auto scrit_expect = std::vector{ 0.21, 0.21, 0.21, 0.21 }; + const auto smax_expect = std::vector{ 0.8 , 0.8 , 0.8 , 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SWFN_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + // 1e20 is a sentinel value that counts as row "ignored". + t.data = std::vector{ + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + -1.0e20, -1.0e+100, 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.0, // 7 + 0.3 , 0.1 , 0.0, // 8 + 0.7 , 0.5 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + 1.0e20 , 1.0e+100 , 0.0, // 16 + + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + 0.1 , 0.0 , 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.0, // 7 + 0.3 , 0.1 , 0.0, // 8 + 0.8 , 0.5 , 0.0, // 9 + 1.0e20 , 1.0e+100 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + 1.0e20 , 1.0e+100 , 0.0, // 16 + + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + -1.0e20, -1.0e+100, 0.0, // 4 + 0.1 , 0.0 , 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.0, // 7 + 0.3 , 0.1 , 0.0, // 8 + 0.5 , 0.35 , 0.0, // 9 + 0.9 , 0.5 , 0.0, // 10 + 1.0e20 , 1.0e+100 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + 1.0e20 , 1.0e+100 , 0.0, // 16 + + // s , kr , pc + -1.0e20, -1.0e+100, 0.0, // 1 + -1.0e20, -1.0e+100, 0.0, // 2 + -1.0e20, -1.0e+100, 0.0, // 3 + 0.0 , 0.0 , 0.0, // 4 + 0.1 , 0.0 , 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.0, // 7 + 0.3 , 0.1 , 0.0, // 8 + 0.5 , 0.35 , 0.0, // 9 + 0.8 , 0.5 , 0.0, // 10 + 0.95 , 0.5 , 0.0, // 11 + 1.0e20 , 1.0e+100 , 0.0, // 12 + 1.0e20 , 1.0e+100 , 0.0, // 13 + 1.0e20 , 1.0e+100 , 0.0, // 14 + 1.0e20 , 1.0e+100 , 0.0, // 15 + 1.0e20 , 1.0e+100 , 0.0, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.1 , 0.1 , 0.0 }; + const auto scrit_expect = std::vector{ 0.21, 0.21, 0.21, 0.21 }; + const auto smax_expect = std::vector{ 0.7 , 0.8 , 0.9 , 0.95 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit = swfunc.criticalSat(ResultColumn{0}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn, sconn_expect); + check_is_close(scrit, scrit_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_CritIsConn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + }; + + t.numRows = 3; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto smax_expect = std::vector{ 0.8, 0.8, 0.8, 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_CritIsConn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.5, // 7 + 0.8 , 0.5 , 0.8, // 8 + 1.0e20 , 1.0e+100, 1.0e+100, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.1 , 0.0 , 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.5, // 7 + 0.8 , 0.5 , 0.8, // 8 + 1.0e20 , 1.0e+100, 1.0e+100, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + 0.0 , 0.0 , 0.0, // 4 + 0.15 , 0.0 , 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.5, // 7 + 0.7 , 0.35 , 0.8, // 8 + 0.8 , 0.5 , 0.8, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.1 , 0.0 , 0.0, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.3 , 0.1 , 0.5, // 7 + 0.8 , 0.5 , 0.8, // 8 + 0.9 , 0.85 , 0.8, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + }; + + t.numRows = 15; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2, 0.1, 0.0, 0.1 }; + const auto scrit_krow_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.2, 0.2, 0.2, 0.2 }; + const auto smax_expect = std::vector{ 0.8, 0.8, 0.8, 0.9 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOGCR_is_Conn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.1, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.1, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.1, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.0 , 0.1, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + }; + + t.numRows = 4; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.2 , 0.2 , 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21, 0.21, 0.21, 0.21 }; + const auto scrit_krog_expect = std::vector{ 0.2 , 0.2 , 0.2 , 0.2 }; + const auto smax_expect = std::vector{ 0.8 , 0.8 , 0.8 , 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOGCR_is_Conn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.1, // 7 + 0.3 , 0.1 , 0.5, // 8 + 0.8 , 0.5 , 0.8, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.1, // 7 + 0.25 , 0.0 , 0.5, // 8 + 0.3 , 0.1 , 0.5, // 9 + 0.8 , 0.5 , 0.8, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.1 , 0.0 , 0.0, // 5 + 0.2 , 0.0 , 0.00001, // 6 + 0.21 , 0.0 , 0.1, // 7 + 0.3 , 0.1 , 0.5, // 8 + 0.8 , 0.5 , 0.8, // 9 + 0.9 , 0.75 , 0.9, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.0 , 0.1, // 7 + 0.3 , 0.0 , 0.25, // 8 + 0.4 , 0.00001 , 0.30, // 9 + 0.5 , 0.1 , 0.5, // 10 + 0.8 , 0.5 , 0.8, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.2 , 0.1 , 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21, 0.25, 0.21, 0.3 }; + const auto scrit_krog_expect = std::vector{ 0.2 , 0.2 , 0.1 , 0.2 }; + const auto smax_expect = std::vector{ 0.8 , 0.8 , 0.9 , 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOWCR_is_Conn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.1 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.1 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.1 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.21, 0.1 , 0.0, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + }; + + t.numRows = 4; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.2 , 0.2 , 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.2 , 0.2 , 0.2 , 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.21, 0.21, 0.21, 0.21 }; + const auto smax_expect = std::vector{ 0.8 , 0.8 , 0.8 , 0.8 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SOWCR_is_Conn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.1 , 0.0, // 7 + 0.3 , 0.1 , 0.5, // 8 + 0.8 , 0.5 , 0.8, // 9 + 1.0e20 , 1.0e+100, 1.0e+100, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.1 , 0.0, // 7 + 0.25 , 0.15 , 0.0, // 8 + 0.3 , 0.1 , 0.5, // 9 + 0.8 , 0.5 , 0.8, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.0 , 0.0 , 0.0, // 5 + 0.1 , 0.000001, 0.0, // 6 + 0.21 , 0.1 , 0.0, // 7 + 0.3 , 0.15 , 0.5, // 8 + 0.8 , 0.5 , 0.8, // 9 + 0.82 , 0.6 , 0.8, // 10 + 0.9 , 0.8 , 0.89, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + -1.0e20, -1.0e+100, -1.0e+100, // 5 + 0.2 , 0.0 , 0.0, // 6 + 0.21 , 0.1 , 0.0, // 7 + 0.3 , 0.1 , 0.5, // 8 + 0.8 , 0.5 , 0.8, // 9 + 0.80001, 0.500001, 0.8, // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.2 , 0.0 , 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.2 , 0.2 , 0.0 , 0.2 }; + const auto scrit_krog_expect = std::vector{ 0.21, 0.25, 0.21, 0.21 }; + const auto smax_expect = std::vector{ 0.8 , 0.8 , 0.9 , 0.80001 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SCR_Not_Conn) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.205, 0.0 , 0.0, + 0.21 , 0.0 , 0.1, + 0.25 , 0.1 , 0.2, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.206, 0.0 , 0.0, + 0.211, 0.0 , 0.1, + 0.25 , 0.1 , 0.2, + 0.3 , 0.1 , 0.5, + 0.8 , 0.5 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.207, 0.0 , 0.0, + 0.212, 0.0 , 0.1, + 0.25 , 0.1 , 0.2, + 0.3 , 0.1 , 0.5, + 0.85 , 0.6 , 0.8, + + // s, krow, krog + 0.2 , 0.0 , 0.0, + 0.209, 0.0 , 0.0, + 0.225, 0.0 , 0.1, + 0.25 , 0.1 , 0.2, + 0.3 , 0.1 , 0.5, + 0.9 , 0.8 , 0.9, + }; + + t.numRows = 6; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.2 , 0.2 , 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21 , 0.211, 0.212, 0.225 }; + const auto scrit_krog_expect = std::vector{ 0.205, 0.206, 0.207, 0.209 }; + const auto smax_expect = std::vector{ 0.8 , 0.8 , 0.85 , 0.9 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_CASE (SOF3_SCR_Not_Conn_LargeNodeAlloc) +{ + auto t = Opm::ECLPropTableRawData{}; + + t.data = std::vector{ + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.2 , 0.0 , 0.0 , // 5 + 0.205 , 0.0 , 0.0 , // 6 + 0.21 , 0.0 , 0.1 , // 7 + 0.25 , 0.1 , 0.2 , // 8 + 0.3 , 0.1 , 0.5 , // 9 + 0.8 , 0.5 , 0.8 , // 10 + 1.0e20 , 1.0e+100, 1.0e+100, // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + 0.1 , 0.0 , 0.0 , // 4 + 0.2 , 0.0 , 0.0 , // 5 + 0.205 , 0.0 , 0.0 , // 6 + 0.21 , 0.0 , 0.1 , // 7 + 0.25 , 0.1 , 0.2 , // 8 + 0.3 , 0.15 , 0.5 , // 9 + 0.5 , 0.25 , 0.6 , // 10 + 0.8 , 0.5 , 0.8 , // 11 + 0.9 , 0.75 , 0.9 , // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.2 , 0.0 , 0.0 , // 5 + 0.205 , 0.000001, 0.000001, // 6 + 0.21 , 0.000002, 0.1 , // 7 + 0.25 , 0.1 , 0.2 , // 8 + 0.3 , 0.1 , 0.5 , // 9 + 0.8 , 0.5 , 0.8 , // 10 + 0.825 , 0.75 , 0.825 , // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + + // s , krow , krog + -1.0e20, -1.0e+100, -1.0e+100, // 1 + -1.0e20, -1.0e+100, -1.0e+100, // 2 + -1.0e20, -1.0e+100, -1.0e+100, // 3 + -1.0e20, -1.0e+100, -1.0e+100, // 4 + 0.2 , 0.0 , 0.0 , // 5 + 0.205 , 0.0 , 0.0 , // 6 + 0.21 , 0.0 , 0.0 , // 7 + 0.25 , 0.0 , 0.0 , // 8 + 0.3 , 0.0 , 0.000001, // 9 + 0.8 , 0.0 , 0.8 , // 10 + 0.99 , 0.000001, 0.99 , // 11 + 1.0e20 , 1.0e+100, 1.0e+100, // 12 + 1.0e20 , 1.0e+100, 1.0e+100, // 13 + 1.0e20 , 1.0e+100, 1.0e+100, // 14 + 1.0e20 , 1.0e+100, 1.0e+100, // 15 + 1.0e20 , 1.0e+100, 1.0e+100, // 16 + }; + + t.numRows = 16; + t.numCols = 3; + t.numTables = 4; + + // Table end-points + const auto sconn_expect = std::vector{ 0.2 , 0.1 , 0.2 , 0.2 }; + const auto scrit_krow_expect = std::vector{ 0.21 , 0.21 , 0.2 , 0.8 }; + const auto scrit_krog_expect = std::vector{ 0.205, 0.205, 0.2 , 0.25 }; + const auto smax_expect = std::vector{ 0.8 , 0.9 , 0.825, 0.99 }; + + // Note: Need to convert input table to column major (Fortran) order + // because that is the format in which PropTable1D expects the tabular + // data. + const auto swfunc = Opm::SatFuncInterpolant(toRawTableFormat(t)); + + using ResultColumn = Opm::SatFuncInterpolant::ResultColumn; + + const auto sconn = swfunc.connateSat(); + const auto scrit_krow = swfunc.criticalSat(ResultColumn{0}); + const auto scrit_krog = swfunc.criticalSat(ResultColumn{1}); + const auto smax = swfunc.maximumSat(); + + check_is_close(sconn , sconn_expect ); + check_is_close(scrit_krow, scrit_krow_expect); + check_is_close(scrit_krog, scrit_krog_expect); + check_is_close(smax , smax_expect ); +} + +BOOST_AUTO_TEST_SUITE_END () diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp index dca1df6d21..71f416f233 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp +++ b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.cpp @@ -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& pv) + const std::vector& 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& injector_tof, const std::vector& producer_tof, - const std::vector& pv) + const std::vector& 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 D2; std::vector 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()); diff --git a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp index 22023ab708..86975628a4 100644 --- a/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp +++ b/ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/DerivedQuantities.hpp @@ -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& pore_volume); + const std::vector& 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& injector_tof, const std::vector& producer_tof, - const std::vector& pore_volume); + const std::vector& pore_volume, + const double max_pv_fraction = 1.0); /// The Lorenz coefficient from the F-Phi curve. ///