diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 9275f36aaf1..dba8c9fd0e9 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -573,6 +573,7 @@ list(APPEND TEST_SOURCE_FILES tests/test_Wells.cpp tests/test_WindowedArray.cpp tests/material/test_2dtables.cpp + tests/material/test_bakedsatfunctables.cpp tests/material/test_eclmateriallawmanager.cpp tests/material/test_hysteresis.cpp tests/material/test_spline.cpp @@ -1361,6 +1362,7 @@ list(APPEND PUBLIC_HEADER_FILES opm/material/fluidmatrixinteractions/EclHysteresisTwoPhaseLawParams.hpp opm/material/fluidmatrixinteractions/EclMaterialLawHystParams.hpp opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp + opm/material/fluidmatrixinteractions/BakedSatfuncTables.hpp opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp opm/material/fluidmatrixinteractions/EclMaterialLawReadEffectiveParams.hpp opm/material/fluidmatrixinteractions/EclMaterialLawTwoPhaseTypes.hpp diff --git a/opm/material/fluidmatrixinteractions/BakedSatfuncTables.hpp b/opm/material/fluidmatrixinteractions/BakedSatfuncTables.hpp new file mode 100644 index 00000000000..2940f253a25 --- /dev/null +++ b/opm/material/fluidmatrixinteractions/BakedSatfuncTables.hpp @@ -0,0 +1,528 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2026 Equinor 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 . +*/ +/*! + * \file + * \brief Baked-table evaluation of relative permeability and capillary + * pressure for the ECL-default three-phase material law. + * + * For each cell the installed material law (multiplexer params from the + * EclMaterialLawManager, endpoint scaling included) is sampled onto minimal + * exact piecewise-linear tables on three axes: + * Sw axis: [krw, pcow-diff] (pcow-diff = pC[w]-pC[o]) + * Sg axis: [krg, pcgo-diff] + * S axis: [kro_ow, kro_go] (S = max(Sw,Swco)+Sg; ECL default 3ph oil) + * Exactness comes from bake-then-refine: candidate nodes are densified (with + * kink snapping) until midpoint probes of every interval reproduce the true + * law, then collinear nodes are pruned. Tables are deduplicated across cells + * by the shape of the scaled EPS info; vertical endpoint scales (maxKr*, + * maxPc*, e.g. from SWATINIT) are factored into per-cell linear multipliers. + * Storage is a small table pool plus one index and a few scalars per cell. + * Correctness is checked at build time against the true law (spotValidate). + */ +#ifndef OPM_BAKED_SATFUNC_TABLES_HPP +#define OPM_BAKED_SATFUNC_TABLES_HPP + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + +// Minimal fluid state carrying only saturations (all the family-I curves read). +template +struct SatOnlyFluidState +{ + std::array sat{}; + Scalar saturation(unsigned phaseIdx) const + { return sat[phaseIdx]; } +}; + +// One two-column curve with precomputed slopes; x strictly increasing on [0,1]. +// Node values interleaved as {v0, s0, v1, s1} so an evaluation after the +// search touches a single cache line (the flowdiagnostics table layout). +struct BakedCurve +{ + std::vector x; + std::vector vals; // 4 doubles per node: v0, s0, v1, s1 + + // exact search acceleration: a uniform hint grid maps the abscissa to a + // starting node; the forward scan is bounded by nodes per hint cell. + static constexpr std::size_t hintCells = 64; + double x0 = 0.0; + double invH = 0.0; + std::array hint{}; + + void buildHint() + { + x0 = x.front(); + const double span = x.back() - x0; + invH = (span > 0.0) ? double(hintCells) / span : 0.0; + std::size_t i = 0; + for (std::size_t j = 0; j <= hintCells; ++j) { + const double xj = x0 + (invH > 0.0 ? double(j) / invH : 0.0); + while (i + 2 < x.size() && x[i + 1] <= xj) { + ++i; + } + hint[j] = static_cast(i); + } + } + + // one hint lookup + short scan serves both columns; derivatives analytic + template + void eval2(const Evaluation& xe, Evaluation& c0, Evaluation& c1) const + { + double xd; + if constexpr (std::is_floating_point_v) { + xd = xe; + } + else { + xd = xe.value(); + } + const double* xb = x.data(); + const double cell = std::min(double(hintCells), std::max(0.0, (xd - x0) * invH)); + std::size_t lo = hint[static_cast(cell)]; + while (lo + 2 < x.size() && xb[lo + 1] <= xd) { + ++lo; + } + const double dx = xd - xb[lo]; + const double* nd = vals.data() + 4 * lo; + if constexpr (std::is_floating_point_v) { + c0 = nd[0] + nd[1] * dx; + c1 = nd[2] + nd[3] * dx; + } + else { + c0.setValue(nd[0] + nd[1] * dx); + c1.setValue(nd[2] + nd[3] * dx); + for (int k = 0; k < xe.size(); ++k) { + c0.setDerivative(k, nd[1] * xe.derivative(k)); + c1.setDerivative(k, nd[3] * xe.derivative(k)); + } + } + } +}; + +struct BakedTable +{ + BakedCurve sw; // krw, pcow-diff + BakedCurve sg; // krg, pcgo-diff + BakedCurve so; // kro_ow, kro_go (abscissa is S = sw+sg) + double swco = 0.0; +}; + +template +class BakedSatfuncTables +{ + using MultiplexerParams = typename Manager::MaterialLawParams; + using DefaultMaterial = typename MaterialLaw::DefaultMaterial; + using DefaultParams = typename DefaultMaterial::Params; + static constexpr int numPhases = 3; + static constexpr unsigned waterIdx = MaterialLaw::Traits::wettingPhaseIdx; + static constexpr unsigned oilIdx = MaterialLaw::Traits::nonWettingPhaseIdx; + static constexpr unsigned gasIdx = MaterialLaw::Traits::gasPhaseIdx; + + using FS = SatOnlyFluidState; + +public: + // truth functions of one variable each (see header comment) + struct Truth + { + const DefaultParams* p; + double swco; + + std::array swAxis(double swv) const + { // krw, pcow-diff at (Sw=swv, Sg=0) + FS fs; fs.sat[waterIdx] = swv; fs.sat[oilIdx] = 1.0 - swv; fs.sat[gasIdx] = 0.0; + std::array pC{}; + DefaultMaterial::capillaryPressures(pC, *p, fs); + const double krw = DefaultMaterial::template krw(*p, fs); + return { krw, pC[waterIdx] - pC[oilIdx] }; + } + std::array sgAxis(double sgv) const + { // krg, pcgo-diff at (Sw=swco, Sg=sgv) + FS fs; fs.sat[waterIdx] = swco; fs.sat[gasIdx] = sgv; + fs.sat[oilIdx] = 1.0 - swco - sgv; + std::array pC{}; + DefaultMaterial::capillaryPressures(pC, *p, fs); + const double krg = DefaultMaterial::template krg(*p, fs); + return { krg, pC[gasIdx] - pC[oilIdx] }; + } + std::array soAxis(double S) const + { // kro_ow(S) at (Sw=S, Sg=0); kro_go(S) at (Sw=swco, Sg=S-swco) + FS a; a.sat[waterIdx] = S; a.sat[oilIdx] = 1.0 - S; a.sat[gasIdx] = 0.0; + FS b; b.sat[waterIdx] = swco; b.sat[gasIdx] = std::max(S - swco, 0.0); + b.sat[oilIdx] = 1.0 - swco - b.sat[gasIdx]; + const double kroOw = DefaultMaterial::template relpermOilInOilWaterSystem(*p, a); + const double kroGo = DefaultMaterial::template relpermOilInOilGasSystem(*p, b); + return { kroOw, kroGo }; + } + }; + + // per-cell factor slots: [krw, kro_ow, pcow, krg, kro_go, pcgo] + static constexpr int numFactors = 6; + + void build(const Manager& manager, std::size_t nc) + { + cellTable_.resize(nc); + cellSwco_.resize(nc); + cellFactors_.assign(numFactors * nc, 1.0); + std::map, std::uint32_t> dedup; + std::vector> tableRefs; + + for (std::size_t c = 0; c < nc; ++c) { + const auto& info = manager.oilWaterScaledEpsInfoDrainage(c); + const int satnum = manager.satnumRegionIdx(c); + // Vertical endpoint scales (maxKr*, maxPc*) are NOT part of the + // key: 2pt vertical scaling — and 3pt with equal Kr*r/maxKr* + // ratios, which the key preserves — is linear in the ordinate and + // factored into per-cell multipliers. Without this, per-cell pc + // (SWATINIT) or kr scaling degenerates to one table per cell. + const auto key = makeKey(satnum, info, + manager.oilWaterConfig(), manager.gasOilConfig()); + const std::array refs = + { info.maxKrw, info.maxKrow, info.maxPcow, + info.maxKrg, info.maxKrog, info.maxPcgo }; + auto [it, inserted] = dedup.try_emplace(key, static_cast(tables_.size())); + if (inserted) { + tables_.push_back(bakeOne(manager, c, info.Swl)); + tableRefs.push_back(refs); + } + cellTable_[c] = it->second; + cellSwco_[c] = tables_[it->second].swco; + for (int k = 0; k < numFactors; ++k) { + cellFactors_[numFactors * c + k] = linFactor(refs[k], tableRefs[it->second][k]); + } + } + + spotValidate(manager, nc); + } + + std::size_t numTables() const { return tables_.size(); } + std::size_t numNodes() const + { + std::size_t n = 0; + for (const auto& t : tables_) { + n += t.sw.x.size() + t.sg.x.size() + t.so.x.size(); + } + return n; + } + double maxBakeError() const { return maxBakeError_; } + + // Fused evaluation: kr[3], pcow-diff, pcgo-diff. One search per axis, + // analytic derivatives, ECL-default 3-phase oil combine (with the same + // epsilon regularization as EclDefaultMaterial::krn). + template + void evaluate(std::size_t cellIdx, + const Evaluation& sw, + const Evaluation& sg, + std::array& kr, // [water, oil, gas] canonical idx below + Evaluation& pcowDiff, + Evaluation& pcgoDiff) const + { + const BakedTable& t = tables_[cellTable_[cellIdx]]; + const double swco = cellSwco_[cellIdx]; + const double* f = &cellFactors_[numFactors * cellIdx]; + + Evaluation krw, krg, kroOw, kroGo; + t.sw.eval2(sw, krw, pcowDiff); + krw *= f[0]; + pcowDiff *= f[2]; + t.sg.eval2(sg, krg, pcgoDiff); + krg *= f[3]; + pcgoDiff *= f[5]; + + const Evaluation swClamped = max(Evaluation(swco), sw); + const Evaluation S = sg + swClamped; + t.so.eval2(S, kroOw, kroGo); + kroOw *= f[1]; // factors applied before the three-phase combine + kroGo *= f[4]; + + Evaluation kro; + constexpr double epsilon = 1e-5; + const double Sval = Opm::getValue(S); + if (Sval - swco < epsilon) { + const Evaluation kro2 = (kroOw + kroGo) / 2; + if (Sval - swco > epsilon / 2) { + const Evaluation kro1 = (sg * kroGo + (swClamped - swco) * kroOw) / (S - swco); + const Evaluation alpha = (epsilon - (S - swco)) / (epsilon / 2); + kro = kro2 * alpha + kro1 * (1 - alpha); + } + else { + kro = kro2; + } + } + else { + kro = (sg * kroGo + (swClamped - swco) * kroOw) / (S - swco); + } + + kr[waterIdx] = krw; + kr[oilIdx] = kro; + kr[gasIdx] = krg; + } + +private: + static double linFactor(double cellV, double tableV) + { + if (tableV > 0.0 && cellV > 0.0) { + const double fac = cellV / tableV; + if (std::isfinite(fac)) { + return fac; + } + } + return 1.0; + } + + template + static std::vector makeKey(int satnum, const Info& i, + const EpsConfig& owCfg, const EpsConfig& goCfg) + { + const auto q = [](double v) { return static_cast(std::llround(v * 1.0e9)); }; + std::vector key + { satnum, + q(i.Swl), q(i.Swcr), q(i.Swu), q(i.Sgl), q(i.Sgcr), q(i.Sgu), + q(i.Sowcr), q(i.Sogcr) }; + // Vertical scales enter only as shape information: with 3pt vertical + // scaling active the residual/max ratio shapes the curve; with plain + // 2pt scaling the magnitude is a pure per-cell factor and does not + // key. A non-positive max cannot be factored and keys the raw pair. + const auto shape = [&](double r, double m, bool threePt) { + if (m > 0.0) { + key.push_back(threePt ? q(r / m) : 0); + } + else { + key.push_back(std::numeric_limits::min()); + key.push_back(q(r)); + key.push_back(q(m)); + } + }; + // gas-oil system: wetting = oil (kro_go), non-wetting = gas (krg) + shape(i.Krwr, i.maxKrw, owCfg.enableThreePointKrwScaling()); + shape(i.Krorw, i.maxKrow, owCfg.enableThreePointKrnScaling()); + shape(0.0, i.maxPcow, false); + shape(i.Krgr, i.maxKrg, goCfg.enableThreePointKrnScaling()); + shape(i.Krorg, i.maxKrog, goCfg.enableThreePointKrwScaling()); + shape(0.0, i.maxPcgo, false); + return key; + } + + // Compare the baked evaluation (incl. the factored per-cell pc scale) + // against the true material law on a sample of cells and states. + void spotValidate(const Manager& manager, std::size_t nc) + { + const std::size_t stride = std::max(nc / 200, 1); + for (std::size_t c = 0; c < nc; c += stride) { + const auto& mp = manager.materialLawParams(c); + const DefaultParams& dp = + mp.template getRealParams(); + for (int a = 0; a <= 3; ++a) { + for (int b = 0; a + b <= 3; ++b) { + const double sw = 0.05 + 0.3 * a; + const double sg = 0.02 + 0.3 * b; + FS fs; + fs.sat[waterIdx] = sw; + fs.sat[gasIdx] = sg; + fs.sat[oilIdx] = 1.0 - sw - sg; + std::array krT{}, pcT{}; + DefaultMaterial::relativePermeabilities(krT, dp, fs); + DefaultMaterial::capillaryPressures(pcT, dp, fs); + std::array krB; + double pcow, pcgo; + evaluate(c, sw, sg, krB, pcow, pcgo); + double err = 0.0; + for (int ph = 0; ph < numPhases; ++ph) { + err = std::max(err, std::abs(krB[ph] - krT[ph]) + / std::max(std::abs(krT[ph]), 1.0)); + } + const double pcowT = pcT[waterIdx] - pcT[oilIdx]; + const double pcgoT = pcT[gasIdx] - pcT[oilIdx]; + err = std::max(err, std::abs(pcow - pcowT) / std::max(std::abs(pcowT), 1.0)); + err = std::max(err, std::abs(pcgo - pcgoT) / std::max(std::abs(pcgoT), 1.0)); + maxBakeError_ = std::max(maxBakeError_, err); + } + } + } + } + + template + BakedCurve bakeCurve(const Fn& truth, std::set& nodes) + { + // refine until every interval's midpoint matches the true law + constexpr double tol = 1.0e-12; + for (int pass = 0; pass < 12; ++pass) { + std::vector xs(nodes.begin(), nodes.end()); + std::vector insert; + for (std::size_t i = 0; i + 1 < xs.size(); ++i) { + const double xm = 0.5 * (xs[i] + xs[i + 1]); + if (xs[i + 1] - xs[i] < 1.0e-10) { + continue; + } + const auto lo = truth(xs[i]); + const auto hi = truth(xs[i + 1]); + const auto mid = truth(xm); + for (int col = 0; col < 2; ++col) { + const double lin = 0.5 * (lo[col] + hi[col]); + const double scale = std::max({ std::abs(lo[col]), std::abs(hi[col]), 1.0 }); + if (std::abs(mid[col] - lin) > tol * scale) { + insert.push_back(xm); + // kink snap: for piecewise-linear truth with one kink in + // the interval, intersecting the secants through the + // left and right quarters yields the kink exactly. + const double w = xs[i + 1] - xs[i]; + const double xq1 = xs[i] + 0.25 * w; + const double xq3 = xs[i] + 0.75 * w; + const auto q1 = truth(xq1); + const auto q3 = truth(xq3); + const double sl = (q1[col] - lo[col]) / (xq1 - xs[i]); + const double sr = (hi[col] - q3[col]) / (xs[i + 1] - xq3); + if (std::abs(sl - sr) > 1.0e-14 * std::max(std::abs(sl), std::abs(sr))) { + // lines: lo + sl (x - xi) and hi + sr (x - xi1) + const double xk = (hi[col] - lo[col] + sl * xs[i] - sr * xs[i + 1]) + / (sl - sr); + if (xk > xs[i] + 1.0e-12 && xk < xs[i + 1] - 1.0e-12) { + insert.push_back(xk); + } + } + break; + } + } + } + if (insert.empty()) { + break; + } + nodes.insert(insert.begin(), insert.end()); + } + + // sample final nodes, then prune collinear interior nodes + std::vector xs(nodes.begin(), nodes.end()); + std::vector> ys(xs.size()); + for (std::size_t i = 0; i < xs.size(); ++i) { + ys[i] = truth(xs[i]); + } + std::vector keep(xs.size(), 1); + for (std::size_t i = 1; i + 1 < xs.size(); ++i) { + bool collinear = true; + // find previous kept node + std::size_t p = i - 1; + while (!keep[p]) { --p; } + const double t = (xs[i] - xs[p]) / (xs[i + 1] - xs[p]); + for (int col = 0; col < 2 && collinear; ++col) { + const double lin = ys[p][col] * (1.0 - t) + ys[i + 1][col] * t; + const double scale = std::max(std::abs(ys[i][col]), 1.0); + collinear = std::abs(lin - ys[i][col]) <= 1.0e-13 * scale; + } + if (collinear) { + keep[i] = 0; + } + } + + BakedCurve c; + std::vector> kv; + for (std::size_t i = 0; i < xs.size(); ++i) { + if (keep[i]) { + c.x.push_back(xs[i]); + kv.push_back(ys[i]); + } + } + const std::size_t n = c.x.size(); + c.vals.assign(4 * n, 0.0); + for (std::size_t i = 0; i < n; ++i) { + const double idx = (i + 1 < n) ? 1.0 / (c.x[i + 1] - c.x[i]) : 0.0; + c.vals[4 * i + 0] = kv[i][0]; + c.vals[4 * i + 1] = (i + 1 < n) ? (kv[i + 1][0] - kv[i][0]) * idx : 0.0; + c.vals[4 * i + 2] = kv[i][1]; + c.vals[4 * i + 3] = (i + 1 < n) ? (kv[i + 1][1] - kv[i][1]) * idx : 0.0; + } + + // validation probe (golden-ratio points; belt and braces) + double maxErr = 0.0; + for (std::size_t i = 0; i + 1 < n; ++i) { + const double xp = c.x[i] + 0.618033988749895 * (c.x[i + 1] - c.x[i]); + const auto tv = truth(xp); + const double dx = xp - c.x[i]; + const double p0 = c.vals[4 * i + 0] + c.vals[4 * i + 1] * dx; + const double p1 = c.vals[4 * i + 2] + c.vals[4 * i + 3] * dx; + maxErr = std::max(maxErr, std::abs(p0 - tv[0]) / std::max(std::abs(tv[0]), 1.0)); + maxErr = std::max(maxErr, std::abs(p1 - tv[1]) / std::max(std::abs(tv[1]), 1.0)); + } + maxBakeError_ = std::max(maxBakeError_, maxErr); + c.buildHint(); + return c; + } + + BakedTable bakeOne(const Manager& manager, std::size_t cellIdx, double swco) + { + const auto& mp = manager.materialLawParams(cellIdx); + const DefaultParams& dp = + mp.template getRealParams(); + Truth truth{ &dp, swco }; + + // candidate nodes: endpoints + scaled EPS anchors + uniform fill; + // the refine loop discovers any interior table kinks. + const auto& info = manager.oilWaterScaledEpsInfoDrainage(cellIdx); + const auto seed = [&](std::initializer_list anchors) { + std::set s; + for (int i = 0; i <= 64; ++i) { + s.insert(double(i) / 64.0); + } + for (double a : anchors) { + if (a > 0.0 && a < 1.0) { + s.insert(a); + } + } + return s; + }; + + BakedTable t; + t.swco = swco; + { + auto nodes = seed({ info.Swl, info.Swcr, info.Swu, 1.0 - info.Sowcr, + 1.0 - info.Sowcr - info.Sgl }); + t.sw = bakeCurve([&](double x) { return truth.swAxis(x); }, nodes); + } + { + auto nodes = seed({ info.Sgl, info.Sgcr, info.Sgu, 1.0 - info.Sogcr - info.Swl }); + t.sg = bakeCurve([&](double x) { return truth.sgAxis(x); }, nodes); + } + { + auto nodes = seed({ info.Swl, info.Swcr, 1.0 - info.Sowcr, info.Swl + info.Sgcr, + 1.0 - info.Sogcr }); + t.so = bakeCurve([&](double x) { return truth.soAxis(x); }, nodes); + } + return t; + } + + std::vector tables_; + std::vector cellTable_; + std::vector cellSwco_; + std::vector cellFactors_; // numFactors per cell + double maxBakeError_ = 0.0; +}; + +} // namespace Opm + +#endif diff --git a/tests/material/test_bakedsatfunctables.cpp b/tests/material/test_bakedsatfunctables.cpp new file mode 100644 index 00000000000..80a948be48e --- /dev/null +++ b/tests/material/test_bakedsatfunctables.cpp @@ -0,0 +1,277 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2026 Equinor 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 . +*/ +/*! + * \file + * \brief Exactness and deduplication tests for BakedSatfuncTables: the baked + * piecewise-linear representation must reproduce the installed + * material law (endpoint scaling included) to near machine precision, + * and per-cell vertical scales must be factored, not multiply tables. + */ +#include "config.h" + +#define BOOST_TEST_MODULE BakedSatfuncTables +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace { + +const std::string commonDeck = + "DIMENS\n" + " 10 10 3 /\n" + "\n" + "TABDIMS\n" + "/\n" + "\n" + "OIL\n" + "GAS\n" + "WATER\n" + "\n" + "DISGAS\n" + "\n" + "FIELD\n" + "\n" + "GRID\n" + "\n" + "DX\n" + " 300*1000 /\n" + "DY\n" + " 300*1000 /\n" + "DZ\n" + " 100*20 100*30 100*50 /\n" + "\n" + "TOPS\n" + " 100*8325 /\n" + "\n" + "PORO\n" + " 300*0.15 /\n" + "PROPS\n" + "\n" + "SWOF\n" + "0.12 0 1 0\n" + "0.18 4.64876033057851E-008 1 0\n" + "0.24 0.000000186 0.997 0\n" + "0.3 4.18388429752066E-007 0.98 0.1\n" + "0.36 7.43801652892562E-007 0.7 0.2\n" + "0.42 1.16219008264463E-006 0.35 0.32\n" + "0.48 1.67355371900826E-006 0.2 0.6\n" + "0.54 2.27789256198347E-006 0.09 1.2\n" + "0.6 2.97520661157025E-006 0.021 2.5\n" + "0.66 3.7654958677686E-006 0.01 3.5\n" + "0.72 4.64876033057851E-006 0.001 4.0\n" + "0.78 0.000005625 0.0001 4.5\n" + "0.84 6.69421487603306E-006 0 5.0\n" + "0.91 8.05914256198347E-006 0 6.0\n" + "1 0.984 0 7.0 /\n" + "\n" + "SGOF\n" + "0 0 1 0\n" + "0.001 0 1 0\n" + "0.02 0 0.997 0.1\n" + "0.05 0.005 0.980 0.2\n" + "0.12 0.025 0.700 0.4\n" + "0.2 0.075 0.350 0.6\n" + "0.25 0.125 0.200 0.7\n" + "0.3 0.190 0.090 0.9\n" + "0.4 0.410 0.021 1.2\n" + "0.45 0.60 0.010 1.5\n" + "0.5 0.72 0.001 1.9\n" + "0.6 0.87 0.0001 2.4\n" + "0.7 0.94 0.000 3.0\n" + "0.85 0.98 0.000 3.5\n" + "0.88 0.984 0.000 4.0 /\n"; + +const std::string fam1DeckString = "RUNSPEC\n\n" + commonDeck; + +// Endpoint scaling with two SWL groups; KRW/PCW vary in four groups so the +// bake must factor them (linear vertical scaling) instead of keying them. +const std::string epsDeckString = + "RUNSPEC\n" + "\n" + "ENDSCALE\n" + "/\n\n" + + commonDeck + + "\n" + "SWL\n" + " 150*0.12 150*0.16 /\n" + "KRW\n" + " 75*0.984 75*0.7 75*0.5 75*0.35 /\n" + "PCW\n" + " 75*7.0 75*5.0 75*3.0 75*2.0 /\n"; + +template +struct Fixture +{ + enum { numPhases = 3 }; + enum { waterPhaseIdx = 0 }; + enum { oilPhaseIdx = 1 }; + enum { gasPhaseIdx = 2 }; + using MaterialTraits = Opm::ThreePhaseMaterialTraits; + using MaterialLawManager = Opm::EclMaterialLaw::Manager; + using MaterialLaw = typename MaterialLawManager::MaterialLaw; + using Baked = Opm::BakedSatfuncTables; +}; + +std::function(const Opm::FieldPropsManager&, const std::string&, bool)> lookup = + [](const Opm::FieldPropsManager& fp, const std::string& prop, bool translate) + { + std::vector dest; + const auto& raw = fp.get_int(prop); + dest.resize(raw.size()); + for (std::size_t i = 0; i < raw.size(); ++i) { + dest[i] = raw[i] - static_cast(translate); + } + return dest; + }; + +std::function identity = [](unsigned i) { return i; }; + +// Compare the baked evaluation against the installed law on a saturation +// sweep for the given cells. +template +void checkExactness(const typename Fixture::MaterialLawManager& manager, + const typename Fixture::Baked& baked, + const std::vector& cells) +{ + using Fix = Fixture; + using MaterialLaw = typename Fix::MaterialLaw; + using DefaultMaterial = typename MaterialLaw::DefaultMaterial; + using FS = Opm::SatOnlyFluidState; + + double maxErr = 0.0; + std::array worst{}; + for (const std::size_t c : cells) { + const auto& mp = manager.materialLawParams(c); + const auto& dp = mp.template getRealParams(); + for (int a = 0; a <= 20; ++a) { + for (int b = 0; a + b <= 20; ++b) { + const Scalar sw = Scalar(a) / 20.0; + const Scalar sg = Scalar(b) / 20.0; + FS fs; + fs.sat[Fix::waterPhaseIdx] = sw; + fs.sat[Fix::gasPhaseIdx] = sg; + fs.sat[Fix::oilPhaseIdx] = 1.0 - sw - sg; + + std::array krT{}, pcT{}; + DefaultMaterial::relativePermeabilities(krT, dp, fs); + DefaultMaterial::capillaryPressures(pcT, dp, fs); + + std::array krB; + Scalar pcow, pcgo; + baked.evaluate(c, sw, sg, krB, pcow, pcgo); + + const auto relErr = [](Scalar bakedV, Scalar trueV) { + if (!std::isfinite(double(trueV))) { + return 0.0; // outside the law's domain; nothing to match + } + return std::abs(double(bakedV - trueV)) + / std::max(std::abs(double(trueV)), 1.0); + }; + double err = 0.0; + for (int ph = 0; ph < Fix::numPhases; ++ph) { + err = std::max(err, relErr(krB[ph], krT[ph])); + } + const Scalar pcowT = pcT[Fix::waterPhaseIdx] - pcT[Fix::oilPhaseIdx]; + const Scalar pcgoT = pcT[Fix::gasPhaseIdx] - pcT[Fix::oilPhaseIdx]; + err = std::max(err, relErr(pcow, pcowT)); + err = std::max(err, relErr(pcgo, pcgoT)); + if (err > maxErr) { + maxErr = err; + worst = { double(sw), double(sg), double(c) }; + } + } + } + } + BOOST_TEST_MESSAGE("max rel err " << maxErr << " at sw=" << worst[0] + << " sg=" << worst[1] << " cell=" << worst[2]); + BOOST_CHECK_MESSAGE(maxErr < 1e-8, + "baked vs law mismatch: max rel err " << maxErr + << " at sw=" << worst[0] << " sg=" << worst[1] + << " cell=" << worst[2]); +} + +} // namespace + +BOOST_AUTO_TEST_CASE(UnscaledFamilyOneIsExactAndDeduplicates) +{ + using Fix = Fixture; + + Opm::Parser parser; + const auto deck = parser.parseString(fam1DeckString); + const Opm::EclipseState eclState(deck); + const std::size_t n = eclState.getInputGrid().getNumActive(); + + typename Fix::MaterialLawManager manager; + manager.initFromState(eclState); + manager.initParamsForElements(eclState, n, lookup, identity); + + typename Fix::Baked baked; + baked.build(manager, n); + + BOOST_CHECK_EQUAL(baked.numTables(), 1U); + BOOST_CHECK_LT(baked.maxBakeError(), 1e-9); + checkExactness(manager, baked, {0, n / 2, n - 1}); +} + +BOOST_AUTO_TEST_CASE(VerticalScalingIsFactoredNotKeyed) +{ + using Fix = Fixture; + + Opm::Parser parser; + const auto deck = parser.parseString(epsDeckString); + const Opm::EclipseState eclState(deck); + const std::size_t n = eclState.getInputGrid().getNumActive(); + + typename Fix::MaterialLawManager manager; + manager.initFromState(eclState); + manager.initParamsForElements(eclState, n, lookup, identity); + + typename Fix::Baked baked; + baked.build(manager, n); + + // two SWL groups shape the curves; the four KRW/PCW groups must be + // absorbed by the per-cell linear factors + BOOST_CHECK_EQUAL(baked.numTables(), 2U); + BOOST_CHECK_LT(baked.maxBakeError(), 1e-9); + + std::vector cells; + for (std::size_t c = 0; c < n; c += 37) { + cells.push_back(c); + } + checkExactness(manager, baked, cells); +}