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

Update Flow Diagnostics Application Library

Browse Source 
Fixes incorrect horizontal and vertical end-point scaling of model's
saturation functions.

API Change: No longer supports user-selected behaviour for treating
scaled end-points with a sentinel value (-1.0E+20).  That option was
introduced due to incomplete understanding of the semantics of the
sentinel value.  Now that we understand the meaning (use actual,
unscaled end-point value from input table), we no longer need the
option.  Update the calling code in RigFlowDiagSolverInterface.cpp
accordingly.
This commit is contained in:
Bård Skaflestad
2018-10-23 15:58:30 +02:00
parent 87622d2e85
commit dc8931f782
11 changed files with 135 additions and 543 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp
View File

@@ -682,7 +682,6 @@ std::vector<RigFlowDiagSolverInterface::RelPermCurve> RigFlowDiagSolverInterface
if (!useEps) { if (!useEps) {
scaling.enable = static_cast<unsigned char>(0); scaling.enable = static_cast<unsigned char>(0);
} }
scaling.invalid = Opm::SatFunc::EPSEvalInterface::InvalidEndpointBehaviour::IgnorePoint;
std::vector<Opm::FlowDiagnostics::Graph> graphArr = m_opmFlowDiagStaticData->m_eclSaturationFunc->getSatFuncCurve(satFuncRequests, static_cast<int>(activeCellIndex), scaling); std::vector<Opm::FlowDiagnostics::Graph> graphArr = m_opmFlowDiagStaticData->m_eclSaturationFunc->getSatFuncCurve(satFuncRequests, static_cast<int>(activeCellIndex), scaling);
for (size_t i = 0; i < graphArr.size(); i++) for (size_t i = 0; i < graphArr.size(); i++)
{ {

1
ThirdParty/custom-opm-flowdiag-app/CMakeLists.txt vendored
View File

@@ -5,7 +5,6 @@ project (custom-opm-flowdiag-app)
include_directories( include_directories(
../custom-opm-flowdiagnostics/opm-flowdiagnostics ../custom-opm-flowdiagnostics/opm-flowdiagnostics
opm-flowdiagnostics-applications opm-flowdiagnostics-applications
opmCore
) )
include (opm-flowdiagnostics-applications/CMakeLists_files.cmake) include (opm-flowdiagnostics-applications/CMakeLists_files.cmake)

65
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractPropCurves.cpp vendored
View File

@@ -54,13 +54,13 @@ namespace {
const auto& x = graph.first; const auto& x = graph.first;
const auto& y = graph.second; const auto& y = graph.second;
os << name << '{' << k << "} = [\n"; os << name << '{' << k << "} = extendTab([\n";
for (auto n = x.size(), i = 0*n; i < n; ++i) { for (auto n = x.size(), i = 0*n; i < n; ++i) {
os << x[i] << ' ' << y[i] << '\n'; os << x[i] << ' ' << y[i] << '\n';
} }
os << "];\n\n"; os << "]);\n\n";
k += 1; k += 1;
} }
@@ -118,7 +118,7 @@ namespace {
const auto graph = const auto graph =
sfunc.getSatFuncCurve(func, activeCell, scaling); sfunc.getSatFuncCurve(func, activeCell, scaling);
printGraph(std::cout, "krg", graph); printGraph(std::cout, "crv.krg", graph);
} }
void krog(const Opm::ECLSaturationFunc& sfunc, void krog(const Opm::ECLSaturationFunc& sfunc,
@@ -140,7 +140,7 @@ namespace {
const auto graph = const auto graph =
sfunc.getSatFuncCurve(func, activeCell, scaling); sfunc.getSatFuncCurve(func, activeCell, scaling);
printGraph(std::cout, "krog", graph); printGraph(std::cout, "crv.krog", graph);
} }
void krow(const Opm::ECLSaturationFunc& sfunc, void krow(const Opm::ECLSaturationFunc& sfunc,
@@ -162,7 +162,7 @@ namespace {
const auto graph = const auto graph =
sfunc.getSatFuncCurve(func, activeCell, scaling); sfunc.getSatFuncCurve(func, activeCell, scaling);
printGraph(std::cout, "krow", graph); printGraph(std::cout, "crv.krow", graph);
} }
void krw(const Opm::ECLSaturationFunc& sfunc, void krw(const Opm::ECLSaturationFunc& sfunc,
@@ -184,7 +184,7 @@ namespace {
const auto graph = const auto graph =
sfunc.getSatFuncCurve(func, activeCell, scaling); sfunc.getSatFuncCurve(func, activeCell, scaling);
printGraph(std::cout, "krw", graph); printGraph(std::cout, "crv.krw", graph);
} }
// ----------------------------------------------------------------- // -----------------------------------------------------------------
@@ -209,7 +209,7 @@ namespace {
const auto graph = const auto graph =
sfunc.getSatFuncCurve(func, activeCell, scaling); sfunc.getSatFuncCurve(func, activeCell, scaling);
printGraph(std::cout, "pcgo", graph); printGraph(std::cout, "crv.pcgo", graph);
} }
void pcow(const Opm::ECLSaturationFunc& sfunc, void pcow(const Opm::ECLSaturationFunc& sfunc,
@@ -231,7 +231,7 @@ namespace {
const auto graph = const auto graph =
sfunc.getSatFuncCurve(func, activeCell, scaling); sfunc.getSatFuncCurve(func, activeCell, scaling);
printGraph(std::cout, "pcow", graph); printGraph(std::cout, "crv.pcow", graph);
} }
// ----------------------------------------------------------------- // -----------------------------------------------------------------
@@ -245,7 +245,7 @@ namespace {
const auto graph = pvtCurves const auto graph = pvtCurves
.getPvtCurve(RC::FVF, Opm::ECLPhaseIndex::Vapour, activeCell); .getPvtCurve(RC::FVF, Opm::ECLPhaseIndex::Vapour, activeCell);
printGraph(std::cout, "Bg", graph); printGraph(std::cout, "crv.Bg", graph);
} }
void mu_g(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, void mu_g(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves,
@@ -256,7 +256,7 @@ namespace {
const auto graph = pvtCurves const auto graph = pvtCurves
.getPvtCurve(RC::Viscosity, Opm::ECLPhaseIndex::Vapour, activeCell); .getPvtCurve(RC::Viscosity, Opm::ECLPhaseIndex::Vapour, activeCell);
printGraph(std::cout, "mu_g", graph); printGraph(std::cout, "crv.mu_g", graph);
} }
void Bo(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, void Bo(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves,
@@ -267,7 +267,7 @@ namespace {
const auto graph = pvtCurves const auto graph = pvtCurves
.getPvtCurve(RC::FVF, Opm::ECLPhaseIndex::Liquid, activeCell); .getPvtCurve(RC::FVF, Opm::ECLPhaseIndex::Liquid, activeCell);
printGraph(std::cout, "Bo", graph); printGraph(std::cout, "crv.Bo", graph);
} }
void mu_o(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, void mu_o(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves,
@@ -278,7 +278,7 @@ namespace {
const auto graph = pvtCurves const auto graph = pvtCurves
.getPvtCurve(RC::Viscosity, Opm::ECLPhaseIndex::Liquid, activeCell); .getPvtCurve(RC::Viscosity, Opm::ECLPhaseIndex::Liquid, activeCell);
printGraph(std::cout, "mu_o", graph); printGraph(std::cout, "crv.mu_o", graph);
} }
// ----------------------------------------------------------------- // -----------------------------------------------------------------
@@ -293,7 +293,7 @@ namespace {
const auto graph = pvtCurves const auto graph = pvtCurves
.getPvtCurve(RC::SaturatedState, PI::Vapour, activeCell); .getPvtCurve(RC::SaturatedState, PI::Vapour, activeCell);
printGraph(std::cout, "rvSat", graph); printGraph(std::cout, "crv.rvSat", graph);
} }
void rsSat(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, void rsSat(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves,
@@ -305,7 +305,7 @@ namespace {
const auto graph = pvtCurves const auto graph = pvtCurves
.getPvtCurve(RC::SaturatedState, PI::Liquid, activeCell); .getPvtCurve(RC::SaturatedState, PI::Liquid, activeCell);
printGraph(std::cout, "rsSat", graph); printGraph(std::cout, "crv.rsSat", graph);
} }
// ----------------------------------------------------------------- // -----------------------------------------------------------------
@@ -340,37 +340,6 @@ namespace {
return Opm::ECLUnits::serialisedUnitConventions(init); return Opm::ECLUnits::serialisedUnitConventions(init);
} }
Opm::SatFunc::EPSEvalInterface::InvalidEndpointBehaviour
handleInvalid(const std::string& behaviour)
{
using IEB = Opm::SatFunc::
EPSEvalInterface::InvalidEndpointBehaviour;
if ((behaviour == "ignore") ||
(behaviour == "ignore-point") ||
(behaviour == "ignore_point") ||
(behaviour == "ignorepoint"))
{
return IEB::IgnorePoint;
}
return IEB::UseUnscaled;
}
auto handleInvalid(const Opm::ParameterGroup& prm)
-> decltype(handleInvalid("ignore"))
{
for (const auto* param : { "handle_invalid" ,
"hInv", "handleInv" })
{
if (prm.has(param)) {
return handleInvalid(prm.get<std::string>(param));
}
}
return handleInvalid("ignore");
}
int getActiveCell(const Opm::ECLGraph& G, int getActiveCell(const Opm::ECLGraph& G,
const Opm::ParameterGroup& prm) const Opm::ParameterGroup& prm)
{ {
@@ -425,7 +394,6 @@ namespace {
auto scaling = Opm::ECLSaturationFunc::SatFuncScaling{}; auto scaling = Opm::ECLSaturationFunc::SatFuncScaling{};
scaling.enable = static_cast<unsigned char>(0); scaling.enable = static_cast<unsigned char>(0);
scaling.invalid = handleInvalid(prm);
if (std::regex_search(useEPS, horiz)) { if (std::regex_search(useEPS, horiz)) {
scaling.enable |= T::Horizontal; scaling.enable |= T::Horizontal;
@@ -460,6 +428,8 @@ try {
pvtCC.setOutputUnits(std::move(units)); pvtCC.setOutputUnits(std::move(units));
} }
std::cout << "function crv = pcurves()\n";
// ----------------------------------------------------------------- // -----------------------------------------------------------------
// Relative permeability // Relative permeability
@@ -470,7 +440,8 @@ try {
// ----------------------------------------------------------------- // -----------------------------------------------------------------
// Capillary pressure // Capillary pressure
if (prm.getDefault("pcgo", false)) { pcgo(sfunc, cellID, scaling); } if (prm.getDefault("pcog", false) || // Alias pcog -> pcgo
prm.getDefault("pcgo", false)) { pcgo(sfunc, cellID, scaling); }
if (prm.getDefault("pcow", false)) { pcow(sfunc, cellID, scaling); } if (prm.getDefault("pcow", false)) { pcow(sfunc, cellID, scaling); }
// ----------------------------------------------------------------- // -----------------------------------------------------------------

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

@@ -1,4 +1,5 @@
/* /*
Copyright 2018 Equinor ASA.
Copyright 2017 Statoil ASA. Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media Project (OPM). This file is part of the Open Porous Media Project (OPM).
@@ -34,8 +35,8 @@
#include <cassert> #include <cassert>
#include <cmath> #include <cmath>
#include <exception> #include <exception>
#include <initializer_list>
#include <iterator> #include <iterator>
#include <limits>
#include <numeric> #include <numeric>
#include <stdexcept> #include <stdexcept>
#include <utility> #include <utility>
@@ -217,21 +218,6 @@ namespace {
return (std::abs(s) < 1.0e+20) ? s : dflt; return (std::abs(s) < 1.0e+20) ? s : dflt;
} }
bool validSaturation(const double s)
{
return (! (s < 0.0)) && (! (s > 1.0));
}
bool validSaturations(std::initializer_list<double> sats)
{
return std::accumulate(std::begin(sats),
std::end (sats), true,
[](const bool result, const double s) -> bool
{
return result && validSaturation(s);
});
}
bool bool
haveScaledRelPermAtCritSat(const ::Opm::ECLGraph& G, haveScaledRelPermAtCritSat(const ::Opm::ECLGraph& G,
const ::Opm::ECLInitFileData& init, const ::Opm::ECLInitFileData& init,
@@ -263,9 +249,7 @@ class Opm::SatFunc::TwoPointScaling::Impl
public: public:
Impl(std::vector<double> smin, Impl(std::vector<double> smin,
std::vector<double> smax) std::vector<double> smax)
: smin_ (std::move(smin)) : smin_(std::move(smin)), smax_(std::move(smax))
, smax_ (std::move(smax))
, handle_invalid_(InvalidEndpointBehaviour::UseUnscaled)
{ {
if (this->smin_.size() != this->smax_.size()) { if (this->smin_.size() != this->smax_.size()) {
throw std::invalid_argument { throw std::invalid_argument {
@@ -287,11 +271,6 @@ private:
std::vector<double> smin_; std::vector<double> smin_;
std::vector<double> smax_; std::vector<double> smax_;
InvalidEndpointBehaviour handle_invalid_;
void handleInvalidEndpoint(const SaturationAssoc& sp,
std::vector<double>& effsat) const;
double sMin(const std::vector<int>::size_type cell, double sMin(const std::vector<int>::size_type cell,
const TableEndPoints& tep) const const TableEndPoints& tep) const
{ {
@@ -327,12 +306,6 @@ Impl::eval(const TableEndPoints& tep,
const auto sLO = this->sMin(cell, tep); const auto sLO = this->sMin(cell, tep);
const auto sHI = this->sMax(cell, tep); const auto sHI = this->sMax(cell, tep);
if (! validSaturations({ sLO, sHI })) {
this->handleInvalidEndpoint(eval_pt, effsat);
continue;
}
effsat.push_back(0.0); effsat.push_back(0.0);
auto& s_eff = effsat.back(); auto& s_eff = effsat.back();
@@ -372,12 +345,6 @@ Impl::reverse(const TableEndPoints& tep,
const auto sLO = this->sMin(cell, tep); const auto sLO = this->sMin(cell, tep);
const auto sHI = this->sMax(cell, tep); const auto sHI = this->sMax(cell, tep);
if (! validSaturations({ sLO, sHI })) {
this->handleInvalidEndpoint(eval_pt, unscaledsat);
continue;
}
unscaledsat.push_back(0.0); unscaledsat.push_back(0.0);
auto& s_unsc = unscaledsat.back(); auto& s_unsc = unscaledsat.back();
@@ -404,26 +371,6 @@ Impl::reverse(const TableEndPoints& tep,
return unscaledsat; return unscaledsat;
} }
void
Opm::SatFunc::TwoPointScaling::Impl::
handleInvalidEndpoint(const SaturationAssoc& sp,
std::vector<double>& effsat) const
{
if (this->handle_invalid_ == InvalidEndpointBehaviour::UseUnscaled) {
// User requests that invalid scaling be treated as unscaled
// saturations. Pick that.
effsat.push_back(sp.sat);
return;
}
if (this->handle_invalid_ == InvalidEndpointBehaviour::IgnorePoint) {
// User requests that invalid scaling be ignored. Signal invalid
// scaled saturation to caller as NaN.
effsat.push_back(std::nan(""));
return;
}
}
// --------------------------------------------------------------------- // ---------------------------------------------------------------------
// Class Opm::SatFunc::PureVerticalScaling::Impl // Class Opm::SatFunc::PureVerticalScaling::Impl
// --------------------------------------------------------------------- // ---------------------------------------------------------------------
@@ -433,6 +380,9 @@ class Opm::SatFunc::PureVerticalScaling::Impl
public: public:
explicit Impl(std::vector<double> fmax) explicit Impl(std::vector<double> fmax)
: fmax_(std::move(fmax)) : fmax_(std::move(fmax))
, inf_ (std::numeric_limits<double>::has_infinity
? std::numeric_limits<double>::infinity()
: std::numeric_limits<double>::max())
{} {}
std::vector<double> std::vector<double>
@@ -442,6 +392,16 @@ public:
private: private:
std::vector<double> fmax_; std::vector<double> fmax_;
double inf_;
std::vector<double>
zeroMaxVal(const SaturationPoints& sp) const;
std::vector<double>
nonZeroMaxVal(const SaturationPoints& sp,
const double maxVal,
std::vector<double> val) const;
}; };
std::vector<double> std::vector<double>
@@ -452,10 +412,40 @@ vertScale(const FunctionValues& f,
{ {
assert (sp.size() == val.size() && "Internal Error in Vertical Scaling"); assert (sp.size() == val.size() && "Internal Error in Vertical Scaling");
auto ret = std::move(val);
const auto maxVal = f.max.val; const auto maxVal = f.max.val;
if (! (std::abs(maxVal) > 0.0)) {
return this->zeroMaxVal(sp);
}
return this->nonZeroMaxVal(sp, maxVal, std::move(val));
}
std::vector<double>
Opm::SatFunc::PureVerticalScaling::Impl::
zeroMaxVal(const SaturationPoints& sp) const
{
auto ret = std::vector<double>(sp.size(), 0.0);
for (auto n = sp.size(), i = 0*n; i < n; ++i) {
const auto fmax = this->fmax_[ sp[i].cell ];
ret[i] = (std::abs(fmax) > 0.0)
? (std::signbit(fmax) ? -this->inf_ : this->inf_)
: 0.0;
}
return ret;
}
std::vector<double>
Opm::SatFunc::PureVerticalScaling::Impl::
nonZeroMaxVal(const SaturationPoints& sp,
const double maxVal,
std::vector<double> val) const
{
auto ret = std::move(val);
for (auto n = sp.size(), i = 0*n; i < n; ++i) { for (auto n = sp.size(), i = 0*n; i < n; ++i) {
ret[i] *= this->fmax_[ sp[i].cell ] / maxVal; ret[i] *= this->fmax_[ sp[i].cell ] / maxVal;
} }
@@ -513,24 +503,29 @@ vertScale(const FunctionValues& f,
const auto c = sp[i].cell; const auto c = sp[i].cell;
const auto s = sp[i].sat; const auto s = sp[i].sat;
const auto sr = this->sdisp_[c]; const auto sr = std::min(this->sdisp_[c], this->smax_[c]);
const auto fr = this->fdisp_[c]; const auto fr = this->fdisp_[c];
const auto sm = this->smax_ [c]; const auto sm = this->smax_ [c];
const auto fm = this->fmax_ [c]; const auto fm = this->fmax_ [c];
if (! (s > sr)) { if (s < sr) {
// s <= sr: Pure vertical scaling in left interval. // s < sr: Pure vertical scaling in left interval.
y *= fr / fdisp; y *= fr / fdisp;
} }
else if (sepfv) { else if (sepfv) {
// s > sr; normal case: Kr(Smax) > Kr(Sr) // s \in [sr, sm), sm > sr; normal case: Kr(Smax) > Kr(Sr).
//
// Linear function between (sr,fr) and (sm,fm) in terms of
// function value 'y'.
const auto t = (y - fdisp) / (fmax - fdisp); const auto t = (y - fdisp) / (fmax - fdisp);
y = fr + t*(fm - fr); y = fr + t*(fm - fr);
} }
else if (s < sm) { else if (s < sm) {
// s > sr; special case: Kr(Smax) == Kr(Sr) in table. Use // s \in [sr, sm), sm > sr; special case: Kr(Smax) == Kr(Sr).
// linear function between (sr,fr) and (sm,fm). //
// Use linear function between (sr,fr) and (sm,fm) in terms of
// saturation value 's'.
const auto t = (s - sr) / (sm - sr); const auto t = (s - sr) / (sm - sr);
y = fr + t*(fm - fr); y = fr + t*(fm - fr);
@@ -554,10 +549,7 @@ public:
Impl(std::vector<double> smin , Impl(std::vector<double> smin ,
std::vector<double> sdisp, std::vector<double> sdisp,
std::vector<double> smax ) std::vector<double> smax )
: smin_ (std::move(smin )) : smin_(std::move(smin)), sdisp_(std::move(sdisp)), smax_(std::move(smax))
, sdisp_ (std::move(sdisp))
, smax_ (std::move(smax ))
, handle_invalid_(InvalidEndpointBehaviour::UseUnscaled)
{ {
if ((this->sdisp_.size() != this->smin_.size()) || if ((this->sdisp_.size() != this->smin_.size()) ||
(this->sdisp_.size() != this->smax_.size())) (this->sdisp_.size() != this->smax_.size()))
@@ -582,11 +574,6 @@ private:
std::vector<double> sdisp_; std::vector<double> sdisp_;
std::vector<double> smax_; std::vector<double> smax_;
InvalidEndpointBehaviour handle_invalid_;
void handleInvalidEndpoint(const SaturationAssoc& sp,
std::vector<double>& effsat) const;
double sMin(const std::vector<int>::size_type cell, double sMin(const std::vector<int>::size_type cell,
const TableEndPoints& tep) const const TableEndPoints& tep) const
{ {
@@ -630,12 +617,6 @@ Impl::eval(const TableEndPoints& tep,
const auto sR = this->sDisp(cell, tep); const auto sR = this->sDisp(cell, tep);
const auto sHI = this->sMax (cell, tep); const auto sHI = this->sMax (cell, tep);
if (! validSaturations({ sLO, sR, sHI })) {
this->handleInvalidEndpoint(eval_pt, effsat);
continue;
}
effsat.push_back(0.0); effsat.push_back(0.0);
auto& s_eff = effsat.back(); auto& s_eff = effsat.back();
@@ -643,22 +624,26 @@ Impl::eval(const TableEndPoints& tep,
// s <= sLO // s <= sLO
s_eff = tep.low; s_eff = tep.low;
} }
else if (! (eval_pt.sat < sHI)) { else if (eval_pt.sat < std::min(sR, sHI)) {
// s >= sHI // s in scaled interval [sLO, sR)
s_eff = tep.high; // Map to tabulated saturation in [tep.low, tep.disp)
}
else if (eval_pt.sat < sR) {
// s \in (sLO, sR)
const auto t = (eval_pt.sat - sLO) / (sR - sLO); const auto t = (eval_pt.sat - sLO) / (sR - sLO);
s_eff = tep.low + t*(tep.disp - tep.low); s_eff = tep.low + t*(tep.disp - tep.low);
} }
else { else if (eval_pt.sat < sHI) {
// s \in (sR, sHI) // s in scaled interval [sR, sHI)
// Map to tabulated saturation in [tep.disp, tep.high)
assert (sHI > sR);
const auto t = (eval_pt.sat - sR) / (sHI - sR); const auto t = (eval_pt.sat - sR) / (sHI - sR);
s_eff = tep.disp + t*(tep.high - tep.disp); s_eff = tep.disp + t*(tep.high - tep.disp);
} }
else {
// s >= sHI
s_eff = tep.high;
}
} }
return effsat; return effsat;
@@ -679,12 +664,6 @@ Impl::reverse(const TableEndPoints& tep,
const auto sR = this->sDisp(cell, tep); const auto sR = this->sDisp(cell, tep);
const auto sHI = this->sMax (cell, tep); const auto sHI = this->sMax (cell, tep);
if (! validSaturations({ sLO, sR, sHI })) {
this->handleInvalidEndpoint(eval_pt, unscaledsat);
continue;
}
unscaledsat.push_back(0.0); unscaledsat.push_back(0.0);
auto& s_unsc = unscaledsat.back(); auto& s_unsc = unscaledsat.back();
@@ -693,54 +672,45 @@ Impl::reverse(const TableEndPoints& tep,
// Map to Minimum Input Saturation in cell (sLO). // Map to Minimum Input Saturation in cell (sLO).
s_unsc = sLO; s_unsc = sLO;
} }
else if (! (eval_pt.sat < tep.high)) {
// s >= maximum tabulated saturation.
// Map to Maximum Input Saturation in cell (sHI).
s_unsc = sHI;
}
else if (eval_pt.sat < tep.disp) { else if (eval_pt.sat < tep.disp) {
// s in tabulated interval (tep.low, tep.disp) // s in tabulated interval [tep.low, tep.disp)
// Map to Input Saturation in (sLO, sR) // Map to Input Saturation in (sLO, sR)
const auto t = const auto t =
(eval_pt.sat - tep.low) (eval_pt.sat - tep.low)
/ (tep.disp - tep.low); / (tep.disp - tep.low);
s_unsc = sLO + t*(sR - sLO); s_unsc = std::min(sLO + t*(sR - sLO), sHI);
} }
else { else if (eval_pt.sat < tep.high) {
// s in tabulated interval (tep.disp, tep.high) // s in tabulated interval [tep.disp, tep.high)
// Map to Input Saturation in (sR, sHI) // Map to Input Saturation in [sR, sHI)
assert (tep.high > tep.disp);
const auto t = const auto t =
(eval_pt.sat - tep.disp) (eval_pt.sat - tep.disp)
/ (tep.high - tep.disp); / (tep.high - tep.disp);
s_unsc = sR + t*(sHI - sR); s_unsc = std::min(sR + t*std::max(sHI - sR, 0.0), sHI);
}
else {
// s >= maximum tabulated saturation.
//
// Map to Maximum Input Saturation in cell (sHI) if maximum
// tabulated saturation is strictly greater than critical
// displacing saturation--otherwise map to critical displacing
// saturation.
//
// Needed to handle cases in which \code tep.disp==tep.high
// \endcode but scaled versions of these might differ, i.e. when
// sR < sHI, but the corresponding saturation points in the
// underlying input table coincide.
s_unsc = (tep.high > tep.disp) ? sHI : sR;
} }
} }
return unscaledsat; return unscaledsat;
} }
void
Opm::SatFunc::ThreePointScaling::Impl::
handleInvalidEndpoint(const SaturationAssoc& sp,
std::vector<double>& effsat) const
{
if (this->handle_invalid_ == InvalidEndpointBehaviour::UseUnscaled) {
// User requests that invalid scaling be treated as unscaled
// saturations. Pick that.
effsat.push_back(sp.sat);
return;
}
if (this->handle_invalid_ == InvalidEndpointBehaviour::IgnorePoint) {
// User requests that invalid scaling be ignored. Signal invalid
// scaled saturation to caller as NaN.
effsat.push_back(std::nan(""));
return;
}
}
// --------------------------------------------------------------------- // ---------------------------------------------------------------------
// EPS factory functions for two-point and three-point scaling options // EPS factory functions for two-point and three-point scaling options
// --------------------------------------------------------------------- // ---------------------------------------------------------------------
@@ -749,7 +719,6 @@ namespace Create {
using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions; using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions;
using RTEP = ::Opm::SatFunc::CreateEPS::RawTableEndPoints; using RTEP = ::Opm::SatFunc::CreateEPS::RawTableEndPoints;
using TEP = ::Opm::SatFunc::EPSEvalInterface::TableEndPoints; using TEP = ::Opm::SatFunc::EPSEvalInterface::TableEndPoints;
using InvBeh = ::Opm::SatFunc::EPSEvalInterface::InvalidEndpointBehaviour;
namespace TwoPoint { namespace TwoPoint {
using EPS = ::Opm::SatFunc::TwoPointScaling; using EPS = ::Opm::SatFunc::TwoPointScaling;
@@ -1742,7 +1711,7 @@ KrGO::twoPointMethod(const ::Opm::ECLGraph& G,
auto t = std::vector<double>(sogcr.size(), 0.0); auto t = std::vector<double>(sogcr.size(), 0.0);
{ {
const auto& sgco = tep.conn.gas; const auto& sgco = tep.conn.gas;
const auto& swco = tep.conn.gas; const auto& swco = tep.conn.water;
const auto& sgcr = tep.crit.gas; const auto& sgcr = tep.crit.gas;
for (auto n = sgcr.size(), i = 0*n; i < n; ++i) { for (auto n = sgcr.size(), i = 0*n; i < n; ++i) {
@@ -1790,7 +1759,7 @@ KrOW::twoPointMethod(const ::Opm::ECLGraph& G,
auto t = std::vector<double>(sowcr.size(), 0.0); auto t = std::vector<double>(sowcr.size(), 0.0);
{ {
const auto& sgco = tep.conn.gas; const auto& sgco = tep.conn.gas;
const auto& swco = tep.conn.gas; const auto& swco = tep.conn.water;
const auto& swcr = tep.crit.water; const auto& swcr = tep.crit.water;
for (auto n = swcr.size(), i = 0*n; i < n; ++i) { for (auto n = swcr.size(), i = 0*n; i < n; ++i) {

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

@@ -67,19 +67,6 @@ namespace Opm { namespace SatFunc {
double sat; double sat;
}; };
/// Policy for how to handle an invalid end-point scaling (e.g., if
/// lower and/or upper scaled saturations have nonsensical values
/// like -1.0E+20).
enum class InvalidEndpointBehaviour {
/// Use the unscaled value for this point.
UseUnscaled,
/// Ignore this scaling request (e.g., produce no graphical
/// output for a scaled saturation function when the scaled
/// end-points are meaningless).
IgnorePoint,
};
/// Convenience type alias. /// Convenience type alias.
using SaturationPoints = std::vector<SaturationAssoc>; using SaturationPoints = std::vector<SaturationAssoc>;
@@ -574,13 +561,6 @@ namespace Opm { namespace SatFunc {
/// auto eps = CreateEPS::fromECLOutput(G, init, opt); /// auto eps = CreateEPS::fromECLOutput(G, init, opt);
/// \endcode /// \endcode
::Opm::ECLPhaseIndex thisPh; ::Opm::ECLPhaseIndex thisPh;
/// How to handle an invalid end-point scaling (e.g., if lower
/// and/or upper scaled saturations have nonsensical values like
/// -1.0E+20).
EPSEvalInterface::InvalidEndpointBehaviour handle_invalid {
EPSEvalInterface::InvalidEndpointBehaviour::UseUnscaled
};
}; };
/// Collection of raw saturation table end points. /// Collection of raw saturation table end points.

2
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp vendored
View File

@@ -90,6 +90,8 @@ namespace {
class PVxGBase class PVxGBase
{ {
public: public:
virtual ~PVxGBase() {}
virtual std::vector<double> virtual std::vector<double>
formationVolumeFactor(const std::vector<double>& rv, formationVolumeFactor(const std::vector<double>& rv,
const std::vector<double>& pg) const = 0; const std::vector<double>& pg) const = 0;

2
ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp vendored
View File

@@ -129,6 +129,8 @@ namespace {
class PVxOBase class PVxOBase
{ {
public: public:
virtual ~PVxOBase() {}
virtual std::vector<double> virtual std::vector<double>
formationVolumeFactor(const std::vector<double>& rs, formationVolumeFactor(const std::vector<double>& rs,
const std::vector<double>& po) const = 0; const std::vector<double>& po) const = 0;

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

@@ -473,12 +473,12 @@ namespace {
explicit InitFileSections(const ecl_file_type* init); explicit InitFileSections(const ecl_file_type* init);
struct Section { struct Section {
Section(const ecl_file_view_type* blk) Section(ecl_file_view_type* blk)
: block (blk) : block (blk)
, first_kw(Details::firstBlockKeyword(block)) , first_kw(Details::firstBlockKeyword(block))
{} {}
const ecl_file_view_type* block; ecl_file_view_type* block;
std::string first_kw; std::string first_kw;
}; };
@@ -502,7 +502,7 @@ namespace {
} }
private: private:
const ecl_file_view_type* init_; ecl_file_view_type* init_;
std::vector<Section> sect_; std::vector<Section> sect_;
}; };
@@ -796,7 +796,7 @@ ECLImpl::InitFileSections::InitFileSections(const ecl_file_type* init)
? sectID - 1 : 0*sectID; ? sectID - 1 : 0*sectID;
// Main section 'sectID': [ start_kw, LGRSGONE ] // Main section 'sectID': [ start_kw, LGRSGONE ]
const auto* sect = auto* sect =
ecl_file_view_add_blockview2(this->init_, start_kw, ecl_file_view_add_blockview2(this->init_, start_kw,
end_kw, start_kw_occurrence); end_kw, start_kw_occurrence);
@@ -808,12 +808,12 @@ ECLImpl::InitFileSections::InitFileSections(const ecl_file_type* init)
const auto occurrence = 0; const auto occurrence = 0;
// Main grid sub-section of 'sectID': [ start_kw, LGR ] // Main grid sub-section of 'sectID': [ start_kw, LGR ]
const auto* main_grid_sect = auto* main_grid_sect =
ecl_file_view_add_blockview2(sect, firstkw.c_str(), ecl_file_view_add_blockview2(sect, firstkw.c_str(),
LGR_KW, occurrence); LGR_KW, occurrence);
// LGR sub-section of 'sectID': [ LGR, LGRSGONE ] // LGR sub-section of 'sectID': [ LGR, LGRSGONE ]
const auto* lgr_sect = auto* lgr_sect =
ecl_file_view_add_blockview2(sect, LGR_KW, ecl_file_view_add_blockview2(sect, LGR_KW,
end_kw, occurrence); end_kw, occurrence);
@@ -965,7 +965,7 @@ private:
/// Saved original active view from host (prior to restricting view /// Saved original active view from host (prior to restricting view
/// to single grid ID). /// to single grid ID).
const ecl_file_view_type* save_; ecl_file_view_type* save_;
}; };
/// Casename prefix. Mostly to implement copy ctor. /// Casename prefix. Mostly to implement copy ctor.
@@ -986,7 +986,7 @@ private:
std::unique_ptr<ECLImpl::GridIDCache> gridIDCache_; std::unique_ptr<ECLImpl::GridIDCache> gridIDCache_;
/// Current active result-set view. /// Current active result-set view.
mutable const ecl_file_view_type* activeBlock_{ nullptr }; mutable ecl_file_view_type* activeBlock_{ nullptr };
/// Support for passing \code *this \endcode to ERT functions that /// Support for passing \code *this \endcode to ERT functions that
/// require an \c ecl_file_type, particularly the function that selects /// require an \c ecl_file_type, particularly the function that selects
@@ -1276,7 +1276,7 @@ private:
/// Sections of the INIT result set. /// Sections of the INIT result set.
ECLImpl::InitFileSections sections_; ECLImpl::InitFileSections sections_;
mutable const ecl_file_view_type* activeBlock_{ nullptr }; mutable ecl_file_view_type* activeBlock_{ nullptr };
/// Negative look-up cache for haveKeywordData() queries. /// Negative look-up cache for haveKeywordData() queries.
mutable MissingKW missing_kw_; mutable MissingKW missing_kw_;
@@ -1415,6 +1415,10 @@ lookup(const std::string& vector, const std::string& gridName) const
} }
} }
// Status of 'vector' unknown for 'gridName'. Actually look for the
// vector in gridName's sections (main grid if gridName.empty(),
// otherwise named local grid).
if (gridName.empty()) { if (gridName.empty()) {
return this->lookupMainGrid(key); return this->lookupMainGrid(key);
} }

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

@@ -759,7 +759,6 @@ private:
{ {
auto opt = Create::EPSOptions{}; auto opt = Create::EPSOptions{};
opt.use3PtScaling = use3PtScaling; opt.use3PtScaling = use3PtScaling;
opt.handle_invalid = SatFuncScaling::IEB::UseUnscaled;
if (active.oil) { if (active.oil) {
this->create_oil_eps(host, G, init, ep, active, opt); this->create_oil_eps(host, G, init, ep, active, opt);

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

@@ -85,19 +85,12 @@ namespace Opm {
Vertical = 1 << 1u, Vertical = 1 << 1u,
}; };
using IEB = SatFunc::EPSEvalInterface::InvalidEndpointBehaviour;
SatFuncScaling() SatFuncScaling()
: enable (Type::Horizontal | Type::Vertical) : enable(Type::Horizontal | Type::Vertical)
, invalid(IEB::UseUnscaled)
{} {}
// Default: Use both Horizontal and Vertical if specified. // Default: Use both Horizontal and Vertical if specified.
unsigned char enable; unsigned char enable;
// Default: Use unscaled values in case of invalid scaled
// saturations occurring in the result set.
IEB invalid;
}; };
/// Constructor /// Constructor

326
ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp vendored
View File

@@ -1,326 +0,0 @@
//===========================================================================
//
// File: Units.hpp
//
// Created: Thu Jul 2 09:19:08 2009
//
// Author(s): Halvor M Nilsen <hnil@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010, 2011, 2012 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_UNITS_HEADER
#define OPM_UNITS_HEADER
/**
The unit sets emplyed in ECLIPSE, in particular the FIELD units,
are quite inconsistent. Ideally one should choose units for a set
of base quantities like Mass,Time and Length and then derive the
units for e.g. pressure and flowrate in a consisten manner. However
that is not the case; for instance in the metric system we have:
[Length] = meters
[time] = days
[mass] = kg
This should give:
[Pressure] = [mass] / ([length] * [time]^2) = kg / (m * days * days)
Instead pressure is given in Bars. When it comes to FIELD units the
number of such examples is long.
*/
namespace Opm {
namespace prefix
/// Conversion prefix for units.
{
constexpr const double micro = 1.0e-6; /**< Unit prefix [\f$\mu\f$] */
constexpr const double milli = 1.0e-3; /**< Unit prefix [m] */
constexpr const double centi = 1.0e-2; /**< Non-standard unit prefix [c] */
constexpr const double deci = 1.0e-1; /**< Non-standard unit prefix [d] */
constexpr const double kilo = 1.0e3; /**< Unit prefix [k] */
constexpr const double mega = 1.0e6; /**< Unit prefix [M] */
constexpr const double giga = 1.0e9; /**< Unit prefix [G] */
} // namespace prefix
namespace unit
/// Definition of various units.
/// All the units are defined in terms of international standard
/// units (SI). Example of use: We define a variable \c k which
/// gives a permeability. We want to set \c k to \f$1\,mD\f$.
/// \code
/// using namespace Opm::unit
/// double k = 0.001*darcy;
/// \endcode
/// We can also use one of the prefixes defined in Opm::prefix
/// \code
/// using namespace Opm::unit
/// using namespace Opm::prefix
/// double k = 1.0*milli*darcy;
/// \endcode
{
///\name Common powers
/// @{
constexpr double square(double v) { return v * v; }
constexpr double cubic (double v) { return v * v * v; }
/// @}
// --------------------------------------------------------------
// Basic (fundamental) units and conversions
// --------------------------------------------------------------
/// \name Length
/// @{
constexpr const double meter = 1;
constexpr const double inch = 2.54 * prefix::centi*meter;
constexpr const double feet = 12 * inch;
/// @}
/// \name Time
/// @{
constexpr const double second = 1;
constexpr const double minute = 60 * second;
constexpr const double hour = 60 * minute;
constexpr const double day = 24 * hour;
constexpr const double year = 365 * day;
/// @}
/// \name Volume
/// @{
constexpr const double gallon = 231 * cubic(inch);
constexpr const double stb = 42 * gallon;
constexpr const double liter = 1 * cubic(prefix::deci*meter);
/// @}
/// \name Mass
/// @{
constexpr const double kilogram = 1;
constexpr const double gram = 1.0e-3 * kilogram;
// http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound
constexpr const double pound = 0.45359237 * kilogram;
/// @}
// --------------------------------------------------------------
// Standardised constants
// --------------------------------------------------------------
/// \name Standardised constant
/// @{
constexpr const double gravity = 9.80665 * meter/square(second);
/// @}
// --------------------------------------------------------------
// Derived units and conversions
// --------------------------------------------------------------
/// \name Force
/// @{
constexpr const double Newton = kilogram*meter / square(second); // == 1
constexpr const double dyne = 1e-5*Newton;
constexpr const double lbf = pound * gravity; // Pound-force
/// @}
/// \name Pressure
/// @{
constexpr const double Pascal = Newton / square(meter); // == 1
constexpr const double barsa = 100000 * Pascal;
constexpr const double atm = 101325 * Pascal;
constexpr const double psia = lbf / square(inch);
/// @}
/// \name Temperature. This one is more complicated
/// because the unit systems used by Eclipse (i.e. degrees
/// Celsius and degrees Fahrenheit require to add or
/// subtract an offset for the conversion between from/to
/// Kelvin
/// @{
constexpr const double degCelsius = 1.0; // scaling factor <20>C -> K
constexpr const double degCelsiusOffset = 273.15; // offset for the <20>C -> K conversion
constexpr const double degFahrenheit = 5.0/9; // scaling factor <20>F -> K
constexpr const double degFahrenheitOffset = 255.37; // offset for the <20>C -> K conversion
/// @}
/// \name Viscosity
/// @{
constexpr const double Pas = Pascal * second; // == 1
constexpr const double Poise = prefix::deci*Pas;
/// @}
namespace perm_details {
constexpr const double p_grad = atm / (prefix::centi*meter);
constexpr const double area = square(prefix::centi*meter);
constexpr const double flux = cubic (prefix::centi*meter) / second;
constexpr const double velocity = flux / area;
constexpr const double visc = prefix::centi*Poise;
constexpr const double darcy = (velocity * visc) / p_grad;
// == 1e-7 [m^2] / 101325
// == 9.869232667160130e-13 [m^2]
}
/// \name Permeability
/// @{
///
/// A porous medium with a permeability of 1 darcy permits a flow (flux)
/// of \f$1\,\mathit{cm}^3/s\f$ of a fluid with viscosity
/// \f$1\,\mathit{cP}\f$ (\f$1\,mPa\cdot s\f$) under a pressure gradient
/// of \f$1\,\mathit{atm}/\mathit{cm}\f$ acting across an area of
/// \f$1\,\mathit{cm}^2\f$.
///
constexpr const double darcy = perm_details::darcy;
/// @}
/**
* Unit conversion routines.
*/
namespace convert {
/**
* Convert from external units of measurements to equivalent
* internal units of measurements. Note: The internal units of
* measurements are *ALWAYS*, and exclusively, SI.
*
* Example: Convert a double @c kx, containing a permeability value
* in units of milli-darcy (mD) to the equivalent value in SI units
* (i.e., \f$m^2\f$).
* \code
* using namespace Opm::unit;
* using namespace Opm::prefix;
* convert::from(kx, milli*darcy);
* \endcode
*
* @param[in] q Physical quantity.
* @param[in] unit Physical unit of measurement.
* @return Value of @c q in equivalent SI units of measurements.
*/
constexpr double from(const double q, const double unit)
{
return q * unit;
}
/**
* Convert from internal units of measurements to equivalent
* external units of measurements. Note: The internal units of
* measurements are *ALWAYS*, and exclusively, SI.
*
* Example: Convert a <CODE>std::vector<double> p</CODE>, containing
* pressure values in the SI unit Pascal (i.e., unit::Pascal) to the
* equivalent values in Psi (unit::psia).
* \code
* using namespace Opm::unit;
* std::transform(p.begin(), p.end(), p.begin(),
* boost::bind(convert::to, _1, psia));
* \endcode
*
* @param[in] q Physical quantity, measured in SI units.
* @param[in] unit Physical unit of measurement.
* @return Value of @c q in unit <CODE>unit</CODE>.
*/
constexpr double to(const double q, const double unit)
{
return q / unit;
}
} // namespace convert
}
namespace Metric {
using namespace prefix;
using namespace unit;
constexpr const double Pressure = barsa;
constexpr const double Temperature = degCelsius;
constexpr const double TemperatureOffset = degCelsiusOffset;
constexpr const double AbsoluteTemperature = degCelsius; // actually [K], but the these two are identical
constexpr const double Length = meter;
constexpr const double Time = day;
constexpr const double Mass = kilogram;
constexpr const double Permeability = milli*darcy;
constexpr const double Transmissibility = centi*Poise*cubic(meter)/(day*barsa);
constexpr const double LiquidSurfaceVolume = cubic(meter);
constexpr const double GasSurfaceVolume = cubic(meter);
constexpr const double ReservoirVolume = cubic(meter);
constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume;
constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume;
constexpr const double Density = kilogram/cubic(meter);
constexpr const double PolymerDensity = kilogram/cubic(meter);
constexpr const double Salinity = kilogram/cubic(meter);
constexpr const double Viscosity = centi*Poise;
constexpr const double Timestep = day;
constexpr const double SurfaceTension = dyne/(centi*meter);
}
namespace Field {
using namespace prefix;
using namespace unit;
constexpr const double Pressure = psia;
constexpr const double Temperature = degFahrenheit;
constexpr const double TemperatureOffset = degFahrenheitOffset;
constexpr const double AbsoluteTemperature = degFahrenheit; // actually [<5B>R], but the these two are identical
constexpr const double Length = feet;
constexpr const double Time = day;
constexpr const double Mass = pound;
constexpr const double Permeability = milli*darcy;
constexpr const double Transmissibility = centi*Poise*stb/(day*psia);
constexpr const double LiquidSurfaceVolume = stb;
constexpr const double GasSurfaceVolume = 1000*cubic(feet);
constexpr const double ReservoirVolume = stb;
constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume;
constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume;
constexpr const double Density = pound/cubic(feet);
constexpr const double PolymerDensity = pound/stb;
constexpr const double Salinity = pound/stb;
constexpr const double Viscosity = centi*Poise;
constexpr const double Timestep = day;
constexpr const double SurfaceTension = dyne/(centi*meter);
}
namespace Lab {
using namespace prefix;
using namespace unit;
constexpr const double Pressure = atm;
constexpr const double Temperature = degCelsius;
constexpr const double TemperatureOffset = degCelsiusOffset;
constexpr const double AbsoluteTemperature = degCelsius; // actually [K], but the these two are identical
constexpr const double Length = centi*meter;
constexpr const double Time = hour;
constexpr const double Mass = gram;
constexpr const double Permeability = milli*darcy;
constexpr const double Transmissibility = centi*Poise*cubic(centi*meter)/(hour*atm);
constexpr const double LiquidSurfaceVolume = cubic(centi*meter);
constexpr const double GasSurfaceVolume = cubic(centi*meter);
constexpr const double ReservoirVolume = cubic(centi*meter);
constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume;
constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume;
constexpr const double Density = gram/cubic(centi*meter);
constexpr const double PolymerDensity = gram/cubic(centi*meter);
constexpr const double Salinity = gram/cubic(centi*meter);
constexpr const double Viscosity = centi*Poise;
constexpr const double Timestep = hour;
constexpr const double SurfaceTension = dyne/(centi*meter);
}
}
#endif // OPM_UNITS_HEADER