diff --git a/tests/data/equil_base.DATA b/tests/data/equil_base.DATA new file mode 100644 index 000000000..740ae0850 --- /dev/null +++ b/tests/data/equil_base.DATA @@ -0,0 +1,115 @@ +RUNSPEC + +WATER +GAS +OIL + +METRIC + +DIMENS + 10 1 10 / + +GRID + +DX + 100*1 / +DY + 100*1 / +DZ + 100*1 / + +TOPS + 10*0. / + +PORO + 100*0.3 / + +PERMX + 100*500 / + +PROPS + +PVTW + 4017.55 1.038 3.22E-6 0.318 0.0 / + +ROCK + 14.7 3E-6 / + +SWOF +0.12 0 1 0 +0.18 4.64876033057851E-008 1 0 +0.24 0.000000186 0.997 0 +0.3 4.18388429752066E-007 0.98 0 +0.36 7.43801652892562E-007 0.7 0 +0.42 1.16219008264463E-006 0.35 0 +0.48 1.67355371900826E-006 0.2 0 +0.54 2.27789256198347E-006 0.09 0 +0.6 2.97520661157025E-006 0.021 0 +0.66 3.7654958677686E-006 0.01 0 +0.72 4.64876033057851E-006 0.001 0 +0.78 0.000005625 0.0001 0 +0.84 6.69421487603306E-006 0 0 +0.91 8.05914256198347E-006 0 0 +1 0.00001 0 0 / + + +SGOF +0 0 1 0 +0.001 0 1 0 +0.02 0 0.997 0 +0.05 0.005 0.980 0 +0.12 0.025 0.700 0 +0.2 0.075 0.350 0 +0.25 0.125 0.200 0 +0.3 0.190 0.090 0 +0.4 0.410 0.021 0 +0.45 0.60 0.010 0 +0.5 0.72 0.001 0 +0.6 0.87 0.0001 0 +0.7 0.94 0.000 0 +0.85 0.98 0.000 0 +0.88 0.984 0.000 0 / + +DENSITY + 53.66 64.49 0.0533 / + +PVDG +14.700 166.666 0.008000 +264.70 12.0930 0.009600 +514.70 6.27400 0.011200 +1014.7 3.19700 0.014000 +2014.7 1.61400 0.018900 +2514.7 1.29400 0.020800 +3014.7 1.08000 0.022800 +4014.7 0.81100 0.026800 +5014.7 0.64900 0.030900 +9014.7 0.38600 0.047000 / + +PVTO +0.0010 14.7 1.0620 1.0400 / +0.0905 264.7 1.1500 0.9750 / +0.1800 514.7 1.2070 0.9100 / +0.3710 1014.7 1.2950 0.8300 / +0.6360 2014.7 1.4350 0.6950 / +0.7750 2514.7 1.5000 0.6410 / +0.9300 3014.7 1.5650 0.5940 / +1.2700 4014.7 1.6950 0.5100 + 9014.7 1.5790 0.7400 / +1.6180 5014.7 1.8270 0.4490 + 9014.7 1.7370 0.6310 / +/ + +SOLUTION + +SWAT + 100*0.0 / + +SGAS + 100*0.0 / + +PRESSURE + 100*300.0 / + +SUMMARY + +SCHEDULE diff --git a/tests/data/equil_capillary.DATA b/tests/data/equil_capillary.DATA new file mode 100644 index 000000000..4289102e5 --- /dev/null +++ b/tests/data/equil_capillary.DATA @@ -0,0 +1,87 @@ +-- Most of the following sections are not actually needed by the test, +-- but it is required by the Eclipse reference manual that they are +-- present. Also, the higher level opm-parser classes +-- (i.e. Opm::EclipseState et al.) assume that they are present. + +------------------------------------- +RUNSPEC + +WATER +OIL +GAS + +DIMENS +1 1 20 / + +TABDIMS + 1 1 40 20 1 20 / + +EQLDIMS +-- NTEQUL + 1 / + +------------------------------------- +GRID + +-- Opm::EclipseState assumes that _some_ grid gets defined, so let's +-- specify a fake one... + +DXV +1*1 / + +DYV +1*1 / + +DZV +20*5 / + +DEPTHZ +4*0.0 / + +PORO + 20*0.3 / + +PERMX + 20*500 / + +------------------------------------- +PROPS + +PVDO +100 1.0 1.0 +200 0.9 1.0 +/ + +PVDG +100 0.010 0.1 +200 0.005 0.2 +/ + +PVTW +1.0 1.0 4.0E-5 0.96 0.0 +/ + +SWOF +0.2 0 1 0.4 +1 1 0 0.1 +/ + +SGOF +0 0 1 0.2 +0.8 1 0 0.5 +/ + +DENSITY +700 1000 1 +/ + +------------------------------------- +SOLUTION + +EQUIL +50 150 50 0.25 20 0.35 1* 1* 0 +/ + +------------------------------------- +SCHEDULE +-- empty section diff --git a/tests/data/equil_capillary_overlap.DATA b/tests/data/equil_capillary_overlap.DATA new file mode 100644 index 000000000..b59e42202 --- /dev/null +++ b/tests/data/equil_capillary_overlap.DATA @@ -0,0 +1,128 @@ +NOECHO + +RUNSPEC ====== + +WATER +OIL +GAS + +TABDIMS + 1 1 40 20 1 20 / + +DIMENS +1 1 20 +/ + +WELLDIMS + 30 10 2 30 / + +START + 1 'JAN' 1990 / + +NSTACK + 25 / + +EQLDIMS +-- NTEQUL + 1 / + + +FMTOUT +FMTIN + +GRID ====== + +DXV +1.0 +/ + +DYV +1.0 +/ + +DZV +20*5.0 +/ + + +PORO +20*0.2 +/ + + +PERMZ + 20*1.0 +/ + +PERMY +20*100.0 +/ + +PERMX +20*100.0 +/ + +BOX + 1 1 1 1 1 1 / + +TOPS +0.0 +/ + +PROPS ====== + + +PVDO +100 1.0 1.0 +200 0.9 1.0 +/ + +PVDG +100 0.010 0.1 +200 0.005 0.2 +/ + +SWOF +0.2 0 1 0.9 +1 1 0 0.1 +/ + +SGOF +0 0 1 0.2 +0.8 1 0 0.5 +/ + +PVTW +--RefPres Bw Comp Vw Cv + 1. 1.0 4.0E-5 0.96 0.0 / + + +ROCK +--RefPres Comp + 1. 5.0E-5 / + +DENSITY +700 1000 1 +/ + +SOLUTION ====== + +EQUIL +45 150 50 0.25 45 0.35 1* 1* 0 +/ + +RPTSOL +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=2' / + +SUMMARY ====== +RUNSUM + +SEPARATE + +SCHEDULE ====== + +RPTSCHED +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=3' 'NEWTON=2' / + + +END diff --git a/tests/data/equil_capillary_swatinit.DATA b/tests/data/equil_capillary_swatinit.DATA new file mode 100644 index 000000000..f4e693a58 --- /dev/null +++ b/tests/data/equil_capillary_swatinit.DATA @@ -0,0 +1,102 @@ +-- Most of the following sections are not actually needed by the test, +-- but it is required by the Eclipse reference manual that they are +-- present. Also, the higher level opm-parser classes +-- (i.e. Opm::EclipseState et al.) assume that they are present. + +------------------------------------- +RUNSPEC + +WATER +OIL +GAS + +DIMENS +1 1 20 / + +TABDIMS + 1 1 40 20 1 20 / + +EQLDIMS +-- NTEQUL + 1 / + +START + 1 'JAN' 2015 / +------------------------------------- +GRID + +-- Opm::EclipseState assumes that _some_ grid gets defined, so let's +-- specify a fake one... + +DXV +1 / + +DYV +1 / + +DZ +20*5 / + +TOPS +0 / + +PORO +20*0.3 / + +PERMX +20*500 / + +PERMZ +20*500 / +------------------------------------- +PROPS + +ROCK + 14.7 3E-6 / + +PVDO +100 1.0 1.0 +200 0.9 1.0 +/ + +PVDG +100 0.010 0.1 +200 0.005 0.2 +/ + +PVTW +1.0 1.0 4.0E-5 0.96 0.0 +/ + +SWOF +0.2 0 1 0.4 +1 1 0 0.1 +/ + +SGOF +0 0 1 0.2 +0.8 1 0 0.5 +/ + +DENSITY +700 1000 1 +/ + +SWATINIT + 5*0 + 10*0.5 +5*1 / + +------------------------------------- +SOLUTION + +EQUIL +50 150 50 0.25 20 0.35 1* 1* 0 +/ + +RPTSOL + SWATINIT SWAT SGAS SOIL PCOG PCOW +/ +------------------------------------- +SCHEDULE +-- empty section diff --git a/tests/data/equil_deadfluids.DATA b/tests/data/equil_deadfluids.DATA new file mode 100644 index 000000000..e0d63cb37 --- /dev/null +++ b/tests/data/equil_deadfluids.DATA @@ -0,0 +1,88 @@ +-- Most of the following sections are not actually needed by the test, +-- but it is required by the Eclipse reference manual that they are +-- present. Also, the higher level opm-parser classes +-- (i.e. Opm::EclipseState et al.) assume that they are present. + +------------------------------------- +RUNSPEC + +WATER +GAS +OIL + +METRIC + +DIMENS +1 1 10 / + +TABDIMS + 1 1 40 20 1 20 / + +EQLDIMS +-- NTEQUL + 1 / + +GRID + +-- Opm::EclipseState assumes that _some_ grid gets defined, so let's +-- specify a fake one... + +DXV +1*1 / + +DYV +1*1 / + +DZV +10*1 / + +TOPS +1*0.0 / + +PORO + 10*0.3 / + +PERMX + 10*500 / + +------------------------------------- +PROPS + +PVDO +100 1.0 1.0 +200 0.5 1.0 +/ + +PVDG +100 0.05 0.1 +200 0.02 0.2 +/ + +PVTW +1.0 1.0 4.0E-5 0.96 0.0 +/ + +SWOF +0 0 1 0 +1 1 0 0 +/ + +SGOF +0 0 1 0 +1 1 0 0 +/ + +DENSITY +700 1000 10 +/ + +------------------------------------- +SOLUTION + +EQUIL +5 150 5 0 2 0 1* 1* 0 +/ + +------------------------------------- +SCHEDULE +-- empty section diff --git a/tests/data/equil_livegas.DATA b/tests/data/equil_livegas.DATA new file mode 100644 index 000000000..a0676f52c --- /dev/null +++ b/tests/data/equil_livegas.DATA @@ -0,0 +1,131 @@ +NOECHO + +RUNSPEC ====== + +WATER +OIL +GAS +VAPOIL + +TABDIMS + 1 1 40 20 1 20 / + +DIMENS +1 1 20 +/ + +WELLDIMS + 30 10 2 30 / + +START + 1 'JAN' 1990 / + +NSTACK + 25 / + +EQLDIMS +-- NTEQUL + 1 / + + +FMTOUT +FMTIN + +GRID ====== + +DXV +1.0 +/ + +DYV +1.0 +/ + +DZV +20*5.0 +/ + + +PORO +20*0.2 +/ + + +PERMZ + 20*1.0 +/ + +PERMY +20*100.0 +/ + +PERMX +20*100.0 +/ + +BOX + 1 1 1 1 1 1 / + +TOPS +0.0 +/ + +PROPS ====== + +PVDO +100 1.0 1.0 +200 0.9 1.0 +/ + +PVTG +-- Pg Rv Bg Vg + 100 0.0001 0.010 0.1 + 0.0 0.0104 0.1 / + 200 0.0004 0.005 0.2 + 0.0 0.0054 0.2 / +/ + +SWOF +0.2 0 1 0.9 +1 1 0 0.1 +/ + +SGOF +0 0 1 0.2 +0.8 1 0 0.5 +/ + +PVTW +--RefPres Bw Comp Vw Cv + 1. 1.0 4.0E-5 0.96 0.0 / + + +ROCK +--RefPres Comp + 1. 5.0E-5 / + +DENSITY +700 1000 1 +/ + +SOLUTION ====== + +EQUIL +45 150 50 0.25 45 0.35 1* 1* 0 +/ + +RPTSOL +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' / + +SUMMARY ====== +RUNSUM + +SEPARATE + +SCHEDULE ====== + +RPTSCHED +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' / + + +END diff --git a/tests/data/equil_liveoil.DATA b/tests/data/equil_liveoil.DATA new file mode 100644 index 000000000..e1b417c46 --- /dev/null +++ b/tests/data/equil_liveoil.DATA @@ -0,0 +1,140 @@ +NOECHO + +RUNSPEC ====== + +WATER +OIL +GAS +DISGAS + +TABDIMS + 1 1 40 20 1 20 / + +DIMENS +1 1 20 +/ + +WELLDIMS + 30 10 2 30 / + +START + 1 'JAN' 1990 / + +NSTACK + 25 / + +EQLDIMS +-- NTEQUL + 1 / + + +FMTOUT +FMTIN + +GRID ====== + +DXV +1.0 +/ + +DYV +1.0 +/ + +DZV +20*5.0 +/ + + +PORO +20*0.2 +/ + + +PERMZ + 20*1.0 +/ + +PERMY +20*100.0 +/ + +PERMX +20*100.0 +/ + +BOX + 1 1 1 1 1 1 / + +TOPS +0.0 +/ + +PROPS ====== + + +PVTO +-- Rs Pbub Bo Vo + 0 1. 1.0000 1.20 / + 20 40. 1.0120 1.17 / + 40 80. 1.0255 1.14 / + 60 120. 1.0380 1.11 / + 80 160. 1.0510 1.08 / + 100 200. 1.0630 1.06 / + 120 240. 1.0750 1.03 / + 140 280. 1.0870 1.00 / + 160 320. 1.0985 .98 / + 180 360. 1.1100 .95 / + 200 400. 1.1200 .94 + 500. 1.1189 .94 / + / + +PVDG +100 0.010 0.1 +200 0.005 0.2 +/ + +SWOF +0.2 0 1 0.9 +1 1 0 0.1 +/ + +SGOF +0 0 1 0.2 +0.8 1 0 0.5 +/ + +PVTW +--RefPres Bw Comp Vw Cv + 1. 1.0 4.0E-5 0.96 0.0 / + + +ROCK +--RefPres Comp + 1. 5.0E-5 / + +DENSITY +700 1000 1 +/ + +SOLUTION ====== + +EQUIL +45 150 50 0.25 45 0.35 1* 1* 0 +/ + +RPTSOL +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' / + +SUMMARY ====== +RUNSUM + +SEPARATE + +SCHEDULE ====== + +RPTSCHED +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' / + + +END diff --git a/tests/data/equil_rsvd_and_rvvd.DATA b/tests/data/equil_rsvd_and_rvvd.DATA new file mode 100644 index 000000000..3076524e3 --- /dev/null +++ b/tests/data/equil_rsvd_and_rvvd.DATA @@ -0,0 +1,151 @@ +NOECHO + +RUNSPEC ====== + +WATER +OIL +GAS +DISGAS +VAPOIL + +TABDIMS + 1 1 40 20 1 20 / + +DIMENS +1 1 20 +/ + +WELLDIMS + 30 10 2 30 / + +START + 1 'JAN' 1990 / + +NSTACK + 25 / + +EQLDIMS +-- NTEQUL + 1 / + + +FMTOUT +FMTIN + +GRID ====== + +DXV +1.0 +/ + +DYV +1.0 +/ + +DZV +20*5.0 +/ + + +PORO +20*0.2 +/ + + +PERMZ + 20*1.0 +/ + +PERMY +20*100.0 +/ + +PERMX +20*100.0 +/ + +BOX + 1 1 1 1 1 1 / + +TOPS +0.0 +/ + +PROPS ====== + +PVTO +-- Rs Pbub Bo Vo + 0 1. 1.0000 1.20 / + 20 40. 1.0120 1.17 / + 40 80. 1.0255 1.14 / + 60 120. 1.0380 1.11 / + 80 160. 1.0510 1.08 / + 100 200. 1.0630 1.06 / + 120 240. 1.0750 1.03 / + 140 280. 1.0870 1.00 / + 160 320. 1.0985 .98 / + 180 360. 1.1100 .95 / + 200 400. 1.1200 .94 + 500. 1.1189 .94 / +/ + +PVTG +-- Pg Rv Bg Vg + 100 0.0001 0.010 0.1 + 0.0 0.0104 0.1 / + 200 0.0004 0.005 0.2 + 0.0 0.0054 0.2 / +/ + +SWOF +0.2 0 1 0.9 +1 1 0 0.1 +/ + +SGOF +0 0 1 0.2 +0.8 1 0 0.5 +/ + +PVTW +--RefPres Bw Comp Vw Cv + 1. 1.0 4.0E-5 0.96 0.0 / + + +ROCK +--RefPres Comp + 1. 5.0E-5 / + +DENSITY +700 1000 1 +/ + +SOLUTION ====== + +EQUIL +45 150 50 0.25 45 0.35 1 1 0 +/ + +RSVD + 0 0.0 + 100 100. / + +RVVD + 0. 0. + 100. 0.0001 / + +RPTSOL +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' / + +SUMMARY ====== +RUNSUM + +SEPARATE + +SCHEDULE ====== + +RPTSCHED +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' / + + +END diff --git a/tests/test_equil.cc b/tests/test_equil.cc new file mode 100644 index 000000000..80a7616d3 --- /dev/null +++ b/tests/test_equil.cc @@ -0,0 +1,1014 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#include "config.h" + +#include +#include +#include + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define CHECK(value, expected) \ + { \ + if ((value) != (expected)) \ + std::abort(); \ + } + +#define CHECK_CLOSE(value, expected, reltol) \ + { \ + if (std::fabs((expected) - (value)) > 1e-14 && \ + std::fabs(((expected) - (value))/((expected) + (value))) > reltol) \ + std::abort(); \ + } + +#define REQUIRE(cond) \ + { \ + if (!(cond)) \ + std::abort(); \ + } + +namespace Ewoms { +namespace Properties { +NEW_TYPE_TAG(TestEquilTypeTag, INHERITS_FROM(BlackOilModel, EclBaseProblem)); +}} + +template +std::unique_ptr +initSimulator(const char *filename) +{ + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + + std::string filenameArg = "--ecl-deck-file-name="; + filenameArg += filename; + + const char* argv[] = { + "test_equil", + filenameArg.c_str() + }; + + Ewoms::setupParameters_(/*argc=*/sizeof(argv)/sizeof(argv[0]), argv, /*registerParams=*/false); + + return std::unique_ptr(new Simulator); +} + +template +static void initDefaultFluidSystem() +{ + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + std::vector > Bo = { + { 101353, 1. }, + { 6.21542e+07, 1 } + }; + std::vector > muo = { + { 101353, 1. }, + { 6.21542e+07, 1 } + }; + + std::vector > Bg = { + { 101353, 1. }, + { 6.21542e+07, 1 } + }; + std::vector > mug = { + { 101353, 1. }, + { 6.21542e+07, 1 } + }; + + double rhoRefO = 700; // [kg/m3] + double rhoRefG = 1000; // [kg/m3] + double rhoRefW = 1000; // [kg/m3] + + FluidSystem::initBegin(/*numPvtRegions=*/1); + FluidSystem::setEnableDissolvedGas(false); + FluidSystem::setEnableVaporizedOil(false); + FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, /*regionIdx=*/0); + + auto gasPvt = std::make_shared>(); + gasPvt->setApproach(Opm::GasPvtMultiplexer::DryGasPvt); + auto& dryGasPvt = gasPvt->getRealPvt::DryGasPvt>(); + dryGasPvt.setNumRegions(/*numPvtRegion=*/1); + dryGasPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW); + dryGasPvt.setGasFormationVolumeFactor(/*regionIdx=*/0, Bg); + dryGasPvt.setGasViscosity(/*regionIdx=*/0, mug); + + auto oilPvt = std::make_shared>(); + oilPvt->setApproach(Opm::OilPvtMultiplexer::DeadOilPvt); + auto& deadOilPvt = oilPvt->getRealPvt::DeadOilPvt>(); + deadOilPvt.setNumRegions(/*numPvtRegion=*/1); + deadOilPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW); + deadOilPvt.setInverseOilFormationVolumeFactor(/*regionIdx=*/0, Bo); + deadOilPvt.setOilViscosity(/*regionIdx=*/0, muo); + + auto waterPvt = std::make_shared>(); + waterPvt->setApproach(Opm::WaterPvtMultiplexer::ConstantCompressibilityWaterPvt); + auto& ccWaterPvt = waterPvt->getRealPvt::ConstantCompressibilityWaterPvt>(); + ccWaterPvt.setNumRegions(/*numPvtRegions=*/1); + ccWaterPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW); + ccWaterPvt.setViscosity(/*regionIdx=*/0, 1); + ccWaterPvt.setCompressibility(/*regionIdx=*/0, 0); + + gasPvt->initEnd(); + oilPvt->initEnd(); + waterPvt->initEnd(); + + FluidSystem::setGasPvt(std::move(gasPvt)); + FluidSystem::setOilPvt(std::move(oilPvt)); + FluidSystem::setWaterPvt(std::move(waterPvt)); + FluidSystem::initEnd(); +} + +static Opm::EquilRecord mkEquilRecord( double datd, double datp, + double zwoc, double pcow_woc, + double zgoc, double pcgo_goc ) +{ + using namespace Opm; + + DeckItem dd( "datdep", double() ); + dd.push_back( datd ); + Opm::Dimension dd_dim( "dddim", 1 ); + dd.push_backDimension( dd_dim, dd_dim ); + + DeckItem dp( "datps", double() ); + dp.push_back( datp ); + Opm::Dimension dp_dim( "dpdim", 1 ); + dp.push_backDimension( dp_dim, dp_dim ); + + DeckItem zw( "zwoc", double() ); + zw.push_back( zwoc ); + Opm::Dimension zw_dim( "zwdim", 1 ); + zw.push_backDimension( zw_dim, zw_dim ); + + DeckItem pcow( "pcow", double() ); + pcow.push_back( pcow_woc ); + Opm::Dimension pcow_dim( "pcowdim", 1 ); + pcow.push_backDimension( pcow_dim, pcow_dim ); + + DeckItem zg( "zgoc", double() ); + zg.push_back( zgoc ); + Opm::Dimension zg_dim( "zgdim", 1 ); + zg.push_backDimension( zg_dim, zg_dim ); + + DeckItem pcgo( "pcgo", double() ); + pcgo.push_back( pcgo_goc ); + Opm::Dimension pcgo_dim( "pcgodim", 1 ); + pcgo.push_backDimension( pcgo_dim, pcgo_dim ); + + DeckItem i1( "i1", int() ); + DeckItem i2( "i2", int() ); + DeckItem i3( "i3", int() ); + i1.push_back( 0 ); + i2.push_back( 0 ); + i3.push_back( 0 ); + + DeckRecord rec; + rec.addItem( std::move( dd ) ); + rec.addItem( std::move( dp ) ); + rec.addItem( std::move( zw ) ); + rec.addItem( std::move( pcow ) ); + rec.addItem( std::move( zg ) ); + rec.addItem( std::move( pcgo ) ); + rec.addItem( std::move( i1 ) ); + rec.addItem( std::move( i2 ) ); + rec.addItem( std::move( i3 ) ); + + return EquilRecord( rec ); +} + +void test_PhasePressure() +{ + typedef std::vector PVal; + typedef std::vector PPress; + + auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 ); + + typedef TTAG(TestEquilTypeTag) TypeTag; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + auto simulator = initSimulator("data/equil_base.DATA"); + initDefaultFluidSystem(); + + Ewoms::EQUIL::EquilReg + region(record, + std::make_shared(), + std::make_shared(), + 0); + + std::vector cells(simulator->gridManager().grid().size(0)); + std::iota(cells.begin(), cells.end(), 0); + + const double grav = 10; + const PPress ppress = Ewoms::EQUIL::phasePressures(simulator->gridManager().grid(), region, cells, grav); + + const int first = 0, last = simulator->gridManager().grid().size(0) - 1; + const double reltol = 1.0e-8; + CHECK_CLOSE(ppress[0][first] , 90e3 , reltol); + CHECK_CLOSE(ppress[0][last ] , 180e3 , reltol); + CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol); + CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); +} + +void test_CellSubset() +{ + typedef std::vector PVal; + typedef std::vector PPress; + + typedef TTAG(TestEquilTypeTag) TypeTag; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + auto simulator = initSimulator("data/equil_base.DATA"); + const auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + initDefaultFluidSystem(); + + Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ), + mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) }; + + Ewoms::EQUIL::EquilReg region[] = + { + Ewoms::EQUIL::EquilReg(record[0], + std::make_shared(), + std::make_shared(), + 0) + , + Ewoms::EQUIL::EquilReg(record[0], + std::make_shared(), + std::make_shared(), + 0) + , + Ewoms::EQUIL::EquilReg(record[1], + std::make_shared(), + std::make_shared(), + 0) + , + Ewoms::EQUIL::EquilReg(record[1], + std::make_shared(), + std::make_shared(), + 0) + }; + + const int cdim[] = { 2, 1, 2 }; + int ncoarse = cdim[0]; + for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; } + + std::vector< std::vector > cells(ncoarse); + for (int c = 0; c < simulator->gridManager().grid().size(0); ++c) { + int ci = c; + const int i = ci % grid.cartdims[0]; ci /= grid.cartdims[0]; + const int j = ci % grid.cartdims[1]; + const int k = ci / grid.cartdims[1]; + + const int ic = (i / (grid.cartdims[0] / cdim[0])); + const int jc = (j / (grid.cartdims[1] / cdim[1])); + const int kc = (k / (grid.cartdims[2] / cdim[2])); + const int ix = ic + cdim[0]*(jc + cdim[1]*kc); + + assert ((0 <= ix) && (ix < ncoarse)); + cells[ix].push_back(c); + } + + PPress ppress(2, PVal(simulator->gridManager().grid().size(0), 0)); + for (std::vector< std::vector >::const_iterator + r = cells.begin(), e = cells.end(); + r != e; ++r) + { + const int rno = int(r - cells.begin()); + const double grav = 10; + const PPress p = + Ewoms::EQUIL::phasePressures(simulator->gridManager().grid(), region[rno], *r, grav); + + PVal::size_type i = 0; + for (std::vector::const_iterator + c = r->begin(), ce = r->end(); + c != ce; ++c, ++i) + { + assert (i < p[0].size()); + + ppress[0][*c] = p[0][i]; + ppress[1][*c] = p[1][i]; + } + } + + const int first = 0, last = simulator->gridManager().grid().size(0) - 1; + const double reltol = 1.0e-8; + CHECK_CLOSE(ppress[0][first] , 105e3 , reltol); + CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol); + CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol); + CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); +} + +void test_RegMapping() +{ + typedef std::vector PVal; + typedef std::vector PPress; + + Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ), + mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) }; + + typedef TTAG(TestEquilTypeTag) TypeTag; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + auto simulator = initSimulator("data/equil_base.DATA"); + initDefaultFluidSystem(); + + Ewoms::EQUIL::EquilReg region[] = + { + Ewoms::EQUIL::EquilReg(record[0], + std::make_shared(), + std::make_shared(), + 0) + , + Ewoms::EQUIL::EquilReg(record[0], + std::make_shared(), + std::make_shared(), + 0) + , + Ewoms::EQUIL::EquilReg(record[1], + std::make_shared(), + std::make_shared(), + 0) + , + Ewoms::EQUIL::EquilReg(record[1], + std::make_shared(), + std::make_shared(), + 0) + }; + + std::vector eqlnum(simulator->gridManager().grid().size(0)); + // [ 0 1; 2 3] + { + for (int i = 0; i < 5; ++i) { + for (int j = 0; j < 5; ++j) { + eqlnum[i*10 + j] = 0; + } + for (int j = 5; j < 10; ++j) { + eqlnum[i*10 + j] = 1; + } + } + for (int i = 5; i < 10; ++i) { + for (int j = 0; j < 5; ++j) { + eqlnum[i*10 + j] = 2; + } + for (int j = 5; j < 10; ++j) { + eqlnum[i*10 + j] = 3; + } + } + } + + Ewoms::RegionMapping<> eqlmap(eqlnum); + + PPress ppress(2, PVal(simulator->gridManager().grid().size(0), 0)); + for (const auto& r : eqlmap.activeRegions()) { + const auto& rng = eqlmap.cells(r); + + const int rno = r; + const double grav = 10; + const PPress p = + Ewoms::EQUIL::phasePressures(simulator->gridManager().grid(), region[rno], rng, grav); + + PVal::size_type i = 0; + for (const auto& c : rng) { + assert (i < p[0].size()); + + ppress[0][c] = p[0][i]; + ppress[1][c] = p[1][i]; + + ++i; + } + } + + const int first = 0, last = simulator->gridManager().grid().size(0) - 1; + const double reltol = 1.0e-8; + CHECK_CLOSE(ppress[0][first] , 105e3 , reltol); + CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol); + CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol); + CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); +} + +void test_DeckAllDead() +{ + typedef TTAG(TestEquilTypeTag) TypeTag; + auto simulator = initSimulator("data/equil_deadfluids.DATA"); + const auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + Ewoms::EQUIL::DeckDependent::InitialStateComputer comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 10.0); + const auto& pressures = comp.press(); + REQUIRE(pressures.size() == 3); + REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-3; + CHECK_CLOSE(pressures[0][first] , 1.496329839e7 , reltol); + CHECK_CLOSE(pressures[0][last ] , 1.504526940e7 , reltol); + CHECK_CLOSE(pressures[1][last] , 1.504526940e7 , reltol); +} + +void test_CapillaryInversion() +{ + // Test setup. + typedef typename TTAG(TestEquilTypeTag) TypeTag; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager MaterialLawManager; + auto simulator = initSimulator("data/equil_capillary.DATA"); + + // Test the capillary inversion for oil-water. + const int cell = 0; + const double reltol = 1.0e-7; + { + const int phase = 0; + const bool increasing = false; + const std::vector pc = { 10.0e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.099e5, 0.0e5, -10.0e5 }; + const std::vector s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 }; + REQUIRE(pc.size() == s.size()); + for (size_t i = 0; i < pc.size(); ++i) { + const double s_computed = Ewoms::EQUIL::satFromPc(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing); + CHECK_CLOSE(s_computed, s[i], reltol); + } + } + + // Test the capillary inversion for gas-oil. + { + const int phase = 2; + const bool increasing = true; + const std::vector pc = { 10.0e5, 0.6e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.0e5, -10.0e5 }; + const std::vector s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 }; + REQUIRE(pc.size() == s.size()); + for (size_t i = 0; i < pc.size(); ++i) { + const double s_computed = Ewoms::EQUIL::satFromPc(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing); + CHECK_CLOSE(s_computed, s[i], reltol); + } + } + + // Test the capillary inversion for gas-water. + { + const int water = 0; + const int gas = 2; + const std::vector pc = { 0.9e5, 0.8e5, 0.6e5, 0.4e5, 0.3e5 }; + const std::vector s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 }; + REQUIRE(pc.size() == s.size()); + for (size_t i = 0; i < pc.size(); ++i) { + const double s_computed = Ewoms::EQUIL::satFromSumOfPcs(*simulator->problem().materialLawManager(), water, gas, cell, pc[i]); + CHECK_CLOSE(s_computed, s[i], reltol); + } + } +} + +void test_DeckWithCapillary() +{ + typedef typename TTAG(TestEquilTypeTag) TypeTag; + auto simulator = initSimulator("data/equil_capillary.DATA"); + auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + Ewoms::EQUIL::DeckDependent::InitialStateComputer comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 10.0); + + const auto& pressures = comp.press(); + REQUIRE(pressures.size() == 3); + REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + CHECK_CLOSE(pressures[0][first] , 1.469769063e7 , reltol); + CHECK_CLOSE(pressures[0][last ] , 15452880.328284413 , reltol); + CHECK_CLOSE(pressures[1][last] , 15462880.328284413 , reltol); + + const auto& sats = comp.saturation(); + const std::vector s[3]{ + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.42190294373815257, 0.77800802072306474, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0.0073481611123183965, 0.79272270823081337, 0.8, 0.8, 0.8, 0.8, 0.57809705626184749, 0.22199197927693526, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.79265183888768165, 0.0072772917691866562, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + REQUIRE(sats[phase].size() == s[phase].size()); + for (size_t i = 0; i < s[phase].size(); ++i) { + CHECK_CLOSE(sats[phase][i], s[phase][i], reltol); + } + } +} + +void test_DeckWithCapillaryOverlap() +{ + typedef typename TTAG(TestEquilTypeTag) TypeTag; + auto simulator = initSimulator("data/equil_capillary_overlap.DATA"); + const auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + Ewoms::EQUIL::DeckDependent::InitialStateComputer comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665); + const auto& pressures = comp.press(); + REQUIRE(pressures.size() == 3); + REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + const double reltol_ecl = 1.0; + CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse + CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl); + + CHECK_CLOSE(pressures[0][first] , 14832467.14, reltol); // opm + CHECK_CLOSE(pressures[0][last ] , 15479883.47, reltol); + CHECK_CLOSE(pressures[1][last ] , 15489883.47, reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + + const std::vector s_ecl[3]{// eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22874042, 0.53397995, 0.78454906, 0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20039, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77125955, 0.46602005, 0.015063271, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22892931226886132, 0.53406457830052489, 0.78457075254244724, 0.91539712466977541, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20023624994125844, 0.084602875330224592, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77107068773113863, 0.46593542169947511, 0.015192997516294321, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); + CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol); + } + } +} + +void test_DeckWithLiveOil() +{ + typedef typename TTAG(TestEquilTypeTag) TypeTag; + auto simulator = initSimulator("data/equil_liveoil.DATA"); + const auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + // Initialize the fluid system + Ewoms::EQUIL::DeckDependent::InitialStateComputer comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665); + const auto& pressures = comp.press(); + REQUIRE(pressures.size() == 3); + REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + const double reltol_ecl = 1.0; + CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse + CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl); + + CHECK_CLOSE(pressures[0][first], 1.483246714e7, reltol); // opm + CHECK_CLOSE(pressures[0][last], 1.547991652e7, reltol); + CHECK_CLOSE(pressures[1][first], 1.492246714e7, reltol); + CHECK_CLOSE(pressures[1][last], 1.548991652e7, reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + const std::vector s_ecl[3]{ // eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22898, 0.53422, 0.78470, 0.91531, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20073, 0.08469, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77102, 0.46578, 0.01458, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22916963446461344, 0.53430490523774521, 0.78471886612242092, 0.91528324362210933, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20057438297017782, 0.084716756377890667, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77083036553538653, 0.46569509476225479, 0.014706750907401245, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol); + CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); + } + std::cout << std::endl; + } + + const auto& rs = comp.rs(); + const std::vector rs_opm {74.61233568, 74.64905212, 74.68578656, 74.72253902, // opm + 74.75930951, 74.79609803, 74.83290459, 74.87519876, + 74.96925416, 75.09067512, 75.0, 75.0, + 75.0, 75.0, 75.0, 75.0, + 75.0, 75.0, 75.0, 75.0}; + const std::vector rs_ecl {74.612228, 74.648956, 74.685707, 74.722473, // eclipse + 74.759254, 74.796051, 74.832870, 74.875145, + 74.969231, 75.090706, 75.000000, 75.000000, + 75.000000, 75.000000, 75.000000, 75.000000, + 75.000000, 75.000000, 75.000000, 75.000000}; + for (size_t i = 0; i < rs_opm.size(); ++i) { + //std::cout << std::setprecision(10) << rs[i] << '\n'; + CHECK_CLOSE(rs[i], rs_opm[i], reltol); + CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl); + } +} + +void test_DeckWithLiveGas() +{ + typedef typename TTAG(TestEquilTypeTag) TypeTag; + auto simulator = initSimulator("data/equil_livegas.DATA"); + const auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + Ewoms::EQUIL::DeckDependent::InitialStateComputer comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665); + const auto& pressures = comp.press(); + REQUIRE(pressures.size() == 3); + REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-3; + const double reltol_ecl = 1.0; + CHECK_CLOSE(pressures[0][first], 1.48215e+07, reltol_ecl); // eclipse + CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][first], 1.49115e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl); + + CHECK_CLOSE(pressures[0][first], 1.482150311e7, reltol); // opm + CHECK_CLOSE(pressures[0][last], 1.547988347e7, reltol); + CHECK_CLOSE(pressures[1][first], 1.491150311e7, reltol); + CHECK_CLOSE(pressures[1][last], 1.548988347e7, reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + const std::vector s_ecl[3]{ // eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24285614, 0.53869015, 0.78454906, 0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18311, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75714386, 0.46130988, 0.032345835, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24310545, 0.5388, 0.78458, 0.91540, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18288667, 0.0846, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75689455, 0.4612, 0.03253333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + for (int phase = 0; phase < 3; ++phase) { + REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol); + CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); + } + std::cout << std::endl; + } + + const auto& rv = comp.rv(); + const std::vector rv_opm { // opm + 2.4884509e-4, 2.4910378e-4, 2.4936267e-4, 2.4962174e-4, + 2.4988100e-4, 2.5014044e-4, 2.5040008e-4, 2.5065990e-4, + 2.5091992e-4, 2.5118012e-4, 2.5223082e-4, 2.5105e-4, + 2.5105e-4, 2.5105e-4, 2.5105e-4, 2.5105e-4, + 2.5105e-4, 2.5105e-4, 2.5105e-4, 2.5105e-4}; + + const std::vector rv_ecl { // eclipse + 0.24884584E-03, 0.24910446E-03, 0.24936325E-03, 0.24962222E-03, + 0.24988138E-03, 0.25014076E-03, 0.25040031E-03, 0.25066003E-03, + 0.25091995E-03, 0.25118008E-03, 0.25223137E-03, 0.25104999E-03, + 0.25104999E-03, 0.25104999E-03, 0.25104999E-03, 0.25104999E-03, + 0.25104999E-03, 0.25104999E-03, 0.25104999E-03, 0.25104999E-03}; + + for (size_t i = 0; i < rv_opm.size(); ++i) { + CHECK_CLOSE(rv[i], rv_opm[i], reltol); + CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl); + } +} + +void test_DeckWithRSVDAndRVVD() +{ + typedef typename TTAG(TestEquilTypeTag) TypeTag; + auto simulator = initSimulator("data/equil_rsvd_and_rvvd.DATA"); + const auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + Ewoms::EQUIL::DeckDependent::InitialStateComputer comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665); + const auto& pressures = comp.press(); + REQUIRE(pressures.size() == 3); + REQUIRE(int(pressures[0].size()) == grid.number_of_cells); + + const int first = 0, last = grid.number_of_cells - 1; + // The relative tolerance is too loose to be very useful, + // but the answer we are checking is the result of an ODE + // solver, and it is unclear if we should check it against + // the true answer or something else. + const double reltol = 1.0e-6; + const double reltol_ecl = 1.0; + CHECK_CLOSE(pressures[0][first], 1.48350e+07, reltol_ecl); // eclipse + CHECK_CLOSE(pressures[0][last], 1.54794e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][first], 1.49250e+07, reltol_ecl); + CHECK_CLOSE(pressures[1][last], 1.54894e+07, reltol_ecl); + + CHECK_CLOSE(pressures[0][first], 1.483499660e7, reltol); // opm + CHECK_CLOSE(pressures[0][last], 1.547924516e7, reltol); + CHECK_CLOSE(pressures[1][first], 1.492499660e7, reltol); + CHECK_CLOSE(pressures[1][last], 1.548924516e7, reltol); + + const auto& sats = comp.saturation(); + // std::cout << "Saturations:\n"; + // for (const auto& sat : sats) { + // for (const double s : sat) { + // std::cout << s << ' '; + // } + // std::cout << std::endl; + // } + const std::vector s_ecl[3]{ // eclipse + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22206347, 0.52871972, 0.78150368, 0.91819441, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19656529, 0.081805572, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77793652, 0.47128031, 0.021931054, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + const std::vector s_opm[3]{ // opm + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22231423543119974, 0.52882640735211706, 0.78152142505479982, 0.91816512259416283, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19636279642563928, 0.08183487740583717, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77768576456880023, 0.47117359264788294, 0.022115778519560897, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + for (int phase = 0; phase < 3; ++phase) { + REQUIRE(sats[phase].size() == s_opm[phase].size()); + for (size_t i = 0; i < s_opm[phase].size(); ++i) { + //std::cout << std::setprecision(10) << sats[phase][i] << '\n'; + CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol); + CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl); + } + std::cout << std::endl; + } + + const auto& rs = comp.rs(); + const std::vector rs_opm { // opm + 74.62498302, 74.65959041, 74.69438035, 74.72935336, + 74.76450995, 74.79985061, 74.83537588, 74.87527065, + 74.96863769, 75.08891765, 52.5, 57.5, + 62.5, 67.5, 72.5, 76.45954841, + 76.70621045, 76.95287736, 77.19954913, 77.44622578}; + + const std::vector rs_ecl { // eclipse + 74.625114, 74.659706, 74.694481, 74.729439, + 74.764580, 74.799904, 74.835419, 74.875252, + 74.968628, 75.088951, 52.500000, 57.500000, + 62.500000, 67.500000, 72.500000, 76.168388, + 76.349953, 76.531532, 76.713142, 76.894775,}; + + const auto& rv = comp.rv(); + const std::vector rv_opm { // opm + 2.50e-6, 7.50e-6, 1.25e-5, 1.75e-5, + 2.25e-5, 2.75e-5, 3.25e-5, 3.75e-5, + 4.25e-5, 2.51158386e-4, 2.52203372e-4, 5.75e-5, + 6.25e-5, 6.75e-5, 7.25e-5, 7.75e-5, + 8.25e-5, 8.75e-5, 9.25e-5, 9.75e-5}; + + const std::vector rv_ecl { // eclipse + 0.24999999E-05, 0.74999998E-05, 0.12500000E-04, 0.17500000E-04, + 0.22500000E-04, 0.27500000E-04, 0.32500000E-04, 0.37500002E-04, + 0.42500000E-04, 0.25115837E-03, 0.25220393E-03, 0.57500001E-04, + 0.62500003E-04, 0.67499997E-04, 0.72499999E-04, 0.77500001E-04, + 0.82500002E-04, 0.87499997E-04, 0.92499999E-04, 0.97500000E-04}; + + for (size_t i = 0; i < rv_opm.size(); ++i) { + //std::cout << std::setprecision(10) << rs[i] << '\n'; + CHECK_CLOSE(rs[i], rs_opm[i], reltol); + CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl); + CHECK_CLOSE(rv[i], rv_opm[i], reltol); + CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl); + } +} + +void test_DeckWithSwatinit() +{ +#if 0 + typedef typename TTAG(TestEquilTypeTag) TypeTag; + auto simulator = initSimulator("data/equil_capillary_swatinit.DATA"); + const auto& eclipseState = simulator->gridManager().eclState(); + Opm::GridManager gm(eclipseState.getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + // Create material law manager. + std::vector compressedToCartesianIdx + = Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell); + MaterialLawManager materialLawManager = MaterialLawManager(); + materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx); + + MaterialLawManager materialLawManagerScaled = MaterialLawManager(); + materialLawManagerScaled.initFromDeck(deck, eclipseState, compressedToCartesianIdx); + + // reference saturations + const std::vector s[3]{ + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.42528761746004229, 0.77462669821009045, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0.014813991154779993, 0.78525420807446045, 0.8, 0.8, 0.8, 0.8, 0.57471238253995771, 0.22537330178990955, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.78518600884522005, 0.014745791925539575, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + // sw in cell 1-5 is forced to be 0.2 since swl=0.2 + // sw in cell 13 and 14 is forced to be swu=1 since P_oil - P_wat < 0. + const std::vector swatinit[3]{ + { 0.2, 0.2, 0.2, 0.2, 0.2, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 0, 0, 0, 0.014813991154779993, 0.78525420807446045, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0.8, 0.8, 0.8, 0.78518600884522005, 0.014745791925539575, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } + }; + + // Adjust oil pressure according to gas saturation and cap pressure + typedef Opm::SimpleModularFluidState SatOnlyFluidState; + + SatOnlyFluidState fluidState; + typedef MaterialLawManager::MaterialLaw MaterialLaw; + + // Initialize the fluid system + FluidSystem::initFromDeck(deck, eclipseState); + + // reference pcs + const int numCells = Opm::UgGridHelpers::numCells(grid); + std::vector pc_original(numCells * FluidSystem::numPhases); + for (int c = 0; c < numCells; ++c) { + std::vector pc = {0,0,0}; + double sw = s[0][c]; + double so = s[1][c]; + double sg = s[2][c]; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, so); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg); + const auto& matParams = materialLawManager.materialLawParams(c); + MaterialLaw::capillaryPressures(pc, matParams, fluidState); + pc_original[3*c + 0] = pc[FluidSystem::oilPhaseIdx] - pc[FluidSystem::waterPhaseIdx]; + pc_original[3*c + 1] = 0.0; + pc_original[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx]; + } + + std::vector pc_scaled_truth = pc_original; + + // modify pcow for cell 1 - 12 (where sw is changed due to swatinit) + // for the reference scaled pc. + pc_scaled_truth[3*0 + 0] = 150031.3; + pc_scaled_truth[3*1 + 0] = 136815.6; + pc_scaled_truth[3*2 + 0] = 123612.7; + pc_scaled_truth[3*3 + 0] = 110422.7; + pc_scaled_truth[3*4 + 0] = 97245.4; + pc_scaled_truth[3*5 + 0] = 84081; + pc_scaled_truth[3*6 + 0] = 70929; + pc_scaled_truth[3*7 + 0] = 57791; + pc_scaled_truth[3*8 + 0] = 44665; + pc_scaled_truth[3*9 + 0] = 31552; + pc_scaled_truth[3*10 + 0] = 18451.5; + pc_scaled_truth[3*11 + 0] = 5364.1; + + // compute the initial state + // apply swatinit + Ewoms::EQUIL::DeckDependent::InitialStateComputer compScaled(materialLawManagerScaled, eclipseState, simulator->gridManager().grid(), 9.81, true); + // don't apply swatinit + Ewoms::EQUIL::DeckDependent::InitialStateComputer compUnscaled(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.81, false); + + // compute pc + std::vector pc_scaled(numCells * FluidSystem::numPhases); + for (int c = 0; c < numCells; ++c) { + std::vector pc = {0,0,0}; + double sw = compScaled.saturation().data()[0][c]; + double so = compScaled.saturation().data()[1][c]; + double sg = compScaled.saturation().data()[2][c]; + + fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, so); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg); + const auto& matParams = materialLawManagerScaled.materialLawParams(c); + MaterialLaw::capillaryPressures(pc, matParams, fluidState); + pc_scaled[3*c + 0] = pc[FluidSystem::oilPhaseIdx] - pc[FluidSystem::waterPhaseIdx]; + pc_scaled[3*c + 1] = 0.0; + pc_scaled[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx]; + } + std::vector pc_unscaled(numCells * FluidSystem::numPhases); + for (int c = 0; c < numCells; ++c) { + std::vector pc = {0,0,0}; + double sw = compUnscaled.saturation().data()[0][c]; + double so = compUnscaled.saturation().data()[1][c]; + double sg = compUnscaled.saturation().data()[2][c]; + + fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, so); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg); + + const auto& matParams = materialLawManager.materialLawParams(c); + MaterialLaw::capillaryPressures(pc, matParams, fluidState); + pc_unscaled[3*c + 0] = pc[FluidSystem::oilPhaseIdx] - pc[FluidSystem::waterPhaseIdx]; + pc_unscaled[3*c + 1] = 0.0; + pc_unscaled[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx]; + } + + // test + const double reltol = 1.0e-3; + for (int phase = 0; phase < 3; ++phase) { + for (size_t i = 0; i < 20; ++i) { + CHECK_CLOSE( pc_original[3*i + phase ], pc_unscaled[3*i + phase ], reltol); + CHECK_CLOSE( pc_scaled_truth[3*i + phase], pc_scaled[3*i + phase ], reltol); + } + } + + for (int phase = 0; phase < 3; ++phase) { + for (size_t i = 0; i < 20; ++i) { + CHECK_CLOSE(compUnscaled.saturation()[phase][i], s[phase][i], reltol); + CHECK_CLOSE(compScaled.saturation()[phase][i], swatinit[phase][i], reltol); + } + } +#endif +} + +int main(int argc, char** argv) +{ + Dune::MPIHelper::instance(argc, argv); + + typedef TTAG(TestEquilTypeTag) TypeTag; + Ewoms::registerAllParameters_(); + + test_PhasePressure(); + test_CellSubset(); + test_RegMapping(); + test_DeckAllDead(); + test_CapillaryInversion(); + test_DeckWithCapillary(); + test_DeckWithCapillaryOverlap(); + test_DeckWithLiveOil(); + test_DeckWithLiveGas(); + test_DeckWithRSVDAndRVVD(); + //test_DeckWithSwatinit(); + + return 0; +}