OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-29 12:33:49 -06:00
Code Issues Packages Projects Releases Wiki Activity
58365df3ff
opm-simulators/tests/test_equil.cpp
Bård Skaflestad 2aac0d6cf5 Enable New Saturation Function Consistency Checks 
The new parameter CheckSatfuncConsistency, command line option
--check-satfunc-consistency, allows users to opt into running the
checks.  The option currently defaults to 'false' to reflect the
somewhat experimental nature of the new facility.

The new parameter --num-satfunc-consistency-sample-points allows the
user to select the maximum number of reported failures for each
individual consistency check.  By default, the simulator will report
at most five failures for each check.

We check the unscaled curves if the run does not activate the
end-point scaling option and the scaled curves otherwise.  At
present we're limited to reversible and non-directional saturation
functions for the drainage process, but those restrictions will be
lifted in due time.
2024-10-16 15:03:07 +02:00

1410 lines
64 KiB
C++
Raw Blame History

// -*- 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 <http://www.gnu.org/licenses/>.
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"
#define BOOST_TEST_MODULE Equil
#include <boost/test/unit_test.hpp>
#include <boost/version.hpp>
#if (BOOST_VERSION / 100000 == 1) && ((BOOST_VERSION / 100) % 1000 < 71)
#include <boost/test/floating_point_comparison.hpp>
#else
#include <boost/test/tools/floating_point_comparison.hpp>
#endif
#include <opm/grid/UnstructuredGrid.h>
#include <opm/grid/GridManager.hpp>
#include <opm/grid/cpgrid/GridHelpers.hpp>
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/models/utils/propertysystem.hh>
#include <opm/models/utils/start.hh>
#include <opm/simulators/flow/BlackoilModelParameters.hpp>
#include <opm/simulators/flow/FlowGenericVanguard.hpp>
#include <opm/simulators/flow/FlowProblemBlackoil.hpp>
#include <opm/simulators/flow/FlowProblemBlackoilProperties.hpp>
#include <opm/simulators/flow/equil/EquilibrationHelpers.hpp>
#include <opm/simulators/linalg/parallelbicgstabbackend.hh>
#include <opm/simulators/wells/BlackoilWellModel.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/misc/mpimanager.hh>
#else
#include <dune/common/parallel/mpihelper.hh>
#endif
#include <array>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <memory>
#include <numeric>
#include <string>
#include <vector>
namespace Opm::Properties {
namespace TTag {
struct TestEquilTypeTag {
using InheritsFrom = std::tuple<FlowBaseProblemBlackoil,
BlackOilModel>;
};
struct TestEquilVapwatTypeTag {
using InheritsFrom = std::tuple<FlowBaseProblemBlackoil,
BlackOilModel>;
};
} // namespace TTag
template<class TypeTag>
struct WellModel<TypeTag, TTag::TestEquilTypeTag> {
using type = BlackoilWellModel<TypeTag>;
};
template<class TypeTag>
struct EnableVapwat<TypeTag, TTag::TestEquilTypeTag> {
static constexpr bool value = true;
};
template<class TypeTag>
struct WellModel<TypeTag, TTag::TestEquilVapwatTypeTag> {
using type = BlackoilWellModel<TypeTag>;
};
template<class TypeTag>
struct EnableVapwat<TypeTag, TTag::TestEquilVapwatTypeTag> {
static constexpr bool value = true;
};
} // namespace Opm::Properties
namespace {
template <class TypeTag>
std::unique_ptr<Opm::GetPropType<TypeTag, Opm::Properties::Simulator>>
initSimulator(const char *filename)
{
using Simulator = Opm::GetPropType<TypeTag, Opm::Properties::Simulator>;
const auto filenameArg = std::string {"--ecl-deck-file-name="} + filename;
const char* argv[] = {
"test_equil",
filenameArg.c_str(),
"--check-satfunc-consistency=false",
};
Opm::setupParameters_<TypeTag>(/*argc=*/sizeof(argv)/sizeof(argv[0]), argv, /*registerParams=*/false);
Opm::FlowGenericVanguard::readDeck(filename);
return std::make_unique<Simulator>();
}
template <class GridView>
std::vector<std::pair<double,double>>
cellVerticalExtent(const GridView& gridView)
{
using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
int numElements = gridView.size(/*codim=*/0);
std::vector<std::pair<double,double>> cellZMinMax(numElements);
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
cellZMinMax[elemIdx] = Opm::EQUIL::Details::cellZMinMax<double>(element);
}
return cellZMinMax;
}
template <class TypeTag>
void initDefaultFluidSystem()
{
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
std::vector<std::pair<double, double> > Bo = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
};
std::vector<std::pair<double, double> > muo = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
};
std::vector<std::pair<double, double> > Bg = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
};
std::vector<std::pair<double, double> > 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<Opm::GasPvtMultiplexer<double>>();
gasPvt->setApproach(Opm::GasPvtApproach::DryGas);
auto& dryGasPvt = gasPvt->getRealPvt<Opm::GasPvtApproach::DryGas>();
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<Opm::OilPvtMultiplexer<double>>();
oilPvt->setApproach(Opm::OilPvtApproach::DeadOil);
auto& deadOilPvt = oilPvt->getRealPvt<Opm::OilPvtApproach::DeadOil>();
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<Opm::WaterPvtMultiplexer<double>>();
waterPvt->setApproach(Opm::WaterPvtApproach::ConstantCompressibilityWater);
auto& ccWaterPvt = waterPvt->getRealPvt<Opm::WaterPvtApproach::ConstantCompressibilityWater>();
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();
}
Opm::EquilRecord
mkEquilRecord(const double datd, const double datp,
const double zwoc, const double pcow_woc,
const double zgoc, const double pcgo_goc)
{
return { datd, datp, zwoc, pcow_woc, zgoc, pcgo_goc, true, true, 0, true};
}
template <typename Simulator>
double centerDepth(const Simulator& sim, const std::size_t cell)
{
return Opm::UgGridHelpers::cellCenterDepth(sim.vanguard().grid(), cell);
}
struct EquilFixture {
EquilFixture() {
int argc = boost::unit_test::framework::master_test_suite().argc;
char** argv = boost::unit_test::framework::master_test_suite().argv;
#if HAVE_DUNE_FEM
Dune::Fem::MPIManager::initialize(argc, argv);
#else
Dune::MPIHelper::instance(argc, argv);
#endif
using namespace Opm;
FlowGenericVanguard::setCommunication(std::make_unique<Opm::Parallel::Communication>());
Opm::ThreadManager::registerParameters();
BlackoilModelParameters<double>::registerParameters();
AdaptiveTimeStepping<TypeTag>::registerParameters();
Parameters::Register<Parameters::EnableTerminalOutput>("Dummy added for the well model to compile.");
registerAllParameters_<TypeTag>();
}
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
using Grid = Opm::GetPropType<TypeTag, Opm::Properties::Grid>;
using GridView = Opm::GetPropType<TypeTag, Opm::Properties::GridView>;
using ElementMapper = Opm::GetPropType<TypeTag, Opm::Properties::ElementMapper>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
using Initializer = Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem,
Grid,
GridView,
ElementMapper,
CartesianIndexMapper>;
};
} // Anonymous namespace
BOOST_GLOBAL_FIXTURE(EquilFixture);
BOOST_AUTO_TEST_CASE(PhasePressure)
{
const auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
std::vector<double> x = {0.0,100.0};
std::vector<double> y = {0.0,0.0};
Opm::Tabulated1DFunction<double> trivialSaltVdTable{2,x,y};
std::vector<double> yT = {298.15,298.15};
Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT};
auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>();
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const auto region = Opm::EQUIL::EquilReg<double> {
record,
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0
};
auto vspan = std::array<double, 2>{};
{
auto cells = std::vector<int>(simulator->vanguard().grid().size(0));
std::iota(cells.begin(), cells.end(), 0);
Opm::EQUIL::Details::verticalExtent(cells, cellVerticalExtent(simulator->vanguard().gridView()),
simulator->vanguard().gridView().comm(), vspan);
}
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav };
ptable.equilibrate(region, vspan);
const auto reltol = 1.0e-6;
const auto first = centerDepth(*simulator, 0);
const auto last = centerDepth(*simulator, simulator->vanguard().grid().size(0) - 1);
BOOST_CHECK_CLOSE(ptable.water(first), 90e3 , reltol);
BOOST_CHECK_CLOSE(ptable.water(last) , 180e3 , reltol);
BOOST_CHECK_CLOSE(ptable.oil (first), 103.5e3, reltol);
BOOST_CHECK_CLOSE(ptable.oil (last) , 166.5e3, reltol);
}
BOOST_AUTO_TEST_CASE(CellSubset)
{
using PVal = std::vector<double>;
using PPress = std::vector<PVal>;
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_base.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
initDefaultFluidSystem<TypeTag>();
const Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
std::vector<double> x = {0.0,100.0};
std::vector<double> y = {0.0,0.0};
Opm::Tabulated1DFunction<double> trivialSaltVdTable{2, x, y};
std::vector<double> yT = {298.15,298.15};
Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT};
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const Opm::EQUIL::EquilReg<double> region[] =
{
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
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<int> > cells(ncoarse);
for (int c = 0; c < simulator->vanguard().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);
}
auto vspan = std::array<double, 2>{};
{
auto vspancells = std::vector<int>(simulator->vanguard().grid().size(0));
std::iota(vspancells.begin(), vspancells.end(), 0);
Opm::EQUIL::Details::verticalExtent(vspancells, cellVerticalExtent(simulator->vanguard().gridView()),
simulator->vanguard().gridView().comm(), vspan);
}
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav };
auto ppress = PPress(2, PVal(simulator->vanguard().grid().size(0), 0.0));
for (auto r = cells.begin(), e = cells.end(); r != e; ++r) {
const int rno = int(r - cells.begin());
ptable.equilibrate(region[rno], vspan);
PVal::size_type i = 0;
for (auto c = r->begin(), ce = r->end(); c != ce; ++c, ++i) {
const auto depth = centerDepth(*simulator, *c);
ppress[0][*c] = ptable.water(depth);
ppress[1][*c] = ptable.oil (depth);
}
}
const int first = 0, last = simulator->vanguard().grid().size(0) - 1;
const double reltol = 1.0e-6;
BOOST_CHECK_CLOSE(ppress[FluidSystem::waterPhaseIdx][first] , 105e3 , reltol);
BOOST_CHECK_CLOSE(ppress[FluidSystem::waterPhaseIdx][last ] , 195e3 , reltol);
BOOST_CHECK_CLOSE(ppress[FluidSystem::oilPhaseIdx][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[FluidSystem::oilPhaseIdx][last ] , 166.5e3 , reltol);
}
BOOST_AUTO_TEST_CASE(RegMapping)
{
const Opm::EquilRecord record[] = {
mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ),
};
using PVal = std::vector<double>;
using PPress = std::vector<PVal>;
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_base.DATA");
initDefaultFluidSystem<TypeTag>();
std::vector<double> x = {0.0,100.0};
std::vector<double> y = {0.0,0.0};
Opm::Tabulated1DFunction<double> trivialSaltVdTable{2, x, y};
std::vector<double> yT = {298.15,298.15};
Opm::Tabulated1DFunction<double> trivialTempVdTable{2, x, yT};
using NoMix = Opm::EQUIL::Miscibility::NoMixing<double>;
const Opm::EQUIL::EquilReg<double> region[] =
{
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg<double>(record[0],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg<double>(record[1],
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
std::make_shared<NoMix>(),
trivialTempVdTable,
trivialSaltVdTable,
0)
};
auto vspan = std::array<double, 2>{};
{
auto cells = std::vector<int>(simulator->vanguard().grid().size(0));
std::iota(cells.begin(), cells.end(), 0);
Opm::EQUIL::Details::verticalExtent(cells, cellVerticalExtent(simulator->vanguard().gridView()),
simulator->vanguard().gridView().comm(), vspan);
}
const auto grav = 10.0;
auto ptable = Opm::EQUIL::Details::PressureTable<
FluidSystem, Opm::EQUIL::EquilReg<double>
>{ grav };
std::vector<int> eqlnum(simulator->vanguard().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;
}
}
}
const Opm::RegionMapping<> eqlmap(eqlnum);
auto ppress = PPress(2, PVal(simulator->vanguard().grid().size(0), 0.0));
for (const auto& r : eqlmap.activeRegions()) {
ptable.equilibrate(region[r], vspan);
for (const auto& c : eqlmap.cells(r)) {
const auto depth = centerDepth(*simulator, c);
ppress[0][c] = ptable.water(depth);
ppress[1][c] = ptable.oil (depth);
}
}
const int first = 0, last = simulator->vanguard().grid().size(0) - 1;
const double reltol = 1.0e-6;
BOOST_CHECK_CLOSE(ppress[FluidSystem::waterPhaseIdx][first] , 105e3 , reltol);
BOOST_CHECK_CLOSE(ppress[FluidSystem::waterPhaseIdx][last ] , 195e3 , reltol);
BOOST_CHECK_CLOSE(ppress[FluidSystem::oilPhaseIdx][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[FluidSystem::oilPhaseIdx][last ] , 166.5e3 , reltol);
}
BOOST_AUTO_TEST_CASE(DeckAllDead)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_deadfluids.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(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-1;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first] , 1.496329839e7 , reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last ] , 1.504526940e7 , reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last] , 1.504526940e7 , reltol);
}
BOOST_AUTO_TEST_CASE(CapillaryInversion)
{
// Test setup.
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_capillary.DATA");
// Test the capillary inversion for oil-water.
const int cell = 0;
const double reltol = 1.0e-5;
{
const int phase = FluidSystem::waterPhaseIdx;
const bool increasing = false;
const std::vector<double> pc = { 10.0e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.099e5, 0.0e5, -10.0e5 };
const std::vector<double> s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 };
BOOST_REQUIRE_EQUAL(pc.size(), s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem>(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
// Test the capillary inversion for gas-oil.
{
const int phase = FluidSystem::gasPhaseIdx;
const bool increasing = true;
const std::vector<double> pc = { 10.0e5, 0.6e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.0e5, -10.0e5 };
const std::vector<double> s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 };
BOOST_REQUIRE_EQUAL(pc.size(), s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem>(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
// Test the capillary inversion for gas-water.
{
const int water = FluidSystem::waterPhaseIdx;
const int gas = FluidSystem::gasPhaseIdx;
const std::vector<double> pc = { 0.9e5, 0.8e5, 0.6e5, 0.4e5, 0.3e5 };
const std::vector<double> s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 };
BOOST_REQUIRE_EQUAL(pc.size(), s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromSumOfPcs<FluidSystem>(*simulator->problem().materialLawManager(), water, gas, cell, pc[i]);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE(DeckWithCapillary)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_capillary.DATA");
auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(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-4;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.469769063e7 , reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last ], 15452880.328284413, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx] [last ], 15462880.328284413, reltol);
const auto& sats = comp.saturation();
std::vector<double> s[3];
s[FluidSystem::waterPhaseIdx] = { 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 };
s[FluidSystem::oilPhaseIdx] = { 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 };
s[FluidSystem::gasPhaseIdx] = { 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) {
BOOST_REQUIRE_EQUAL(sats[phase].size(), s[phase].size());
for (size_t i = 0; i < s[phase].size(); ++i) {
BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE(DeckWithCapillaryOverlap)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_capillary_overlap.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(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-4;
const double reltol_ecl = 100.0;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.48324e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.49224e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last], 1.54901e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first] , 14832467.14, reltol); // opm
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last ] , 15479883.47, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][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;
// }
std::vector<double> s_ecl[3]; // eclipse
s_ecl[FluidSystem::waterPhaseIdx] = { 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 };
s_ecl[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20039, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_ecl[FluidSystem::gasPhaseIdx] = { 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 };
std::vector<double> s_opm[3]; // opm
s_opm[FluidSystem::waterPhaseIdx] = { 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 };
s_opm[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20023624994125844, 0.084602875330224592, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_opm[FluidSystem::gasPhaseIdx] = { 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) {
BOOST_REQUIRE_EQUAL(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';
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE(DeckWithLiveOil)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_liveoil.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
// Initialize the fluid system
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(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-4;
const double reltol_ecl = 100.0;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.48324e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.49224e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last], 1.54901e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.483246714e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.547991652e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.492246714e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][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;
// }
std::vector<double> s_ecl[3]; // eclipse
s_ecl[FluidSystem::waterPhaseIdx] = { 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 };
s_ecl[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20073, 0.08469, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_ecl[FluidSystem::gasPhaseIdx] = { 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 };
std::vector<double> s_opm[3]; // opm
s_opm[FluidSystem::waterPhaseIdx] = { 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 };
s_opm[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20057438297017782, 0.084716756377890667, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_opm[FluidSystem::gasPhaseIdx] = { 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) {
BOOST_REQUIRE_EQUAL(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';
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
}
std::cout << std::endl;
}
const auto& rs = comp.rs();
const std::vector<double> 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<double> 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';
BOOST_CHECK_CLOSE(rs[i], rs_opm[i], reltol);
BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_CASE(DeckWithCO2STORE)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator1 = initSimulator<TypeTag>("equil_co2store_go.DATA");
EquilFixture::Initializer comp_go(*simulator1->problem().materialLawManager(),
simulator1->vanguard().eclState(),
simulator1->vanguard().grid(),
simulator1->vanguard().gridView(),
simulator1->vanguard().cartesianMapper(), 9.80665);
auto simulator2 = initSimulator<TypeTag>("equil_co2store_gw.DATA");
EquilFixture::Initializer comp_gw(*simulator2->problem().materialLawManager(),
simulator2->vanguard().eclState(),
simulator2->vanguard().grid(),
simulator2->vanguard().gridView(),
simulator2->vanguard().cartesianMapper(), 9.80665);
Opm::GridManager gm(simulator2->vanguard().eclState().getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
const double reltol = 1.0e-5;
const auto& pressures_go = comp_go.press();
BOOST_REQUIRE_EQUAL(pressures_go.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures_go[0].size()), grid.number_of_cells);
const auto& pressures_gw = comp_gw.press();
BOOST_REQUIRE_EQUAL(pressures_gw.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures_gw[0].size()), grid.number_of_cells);
const auto& sats_go = comp_go.saturation();
const auto& sats_gw = comp_gw.saturation();
for (int i = 0; i < grid.number_of_cells; ++i) {
BOOST_CHECK_CLOSE(pressures_go[FluidSystem::gasPhaseIdx][i], pressures_gw[FluidSystem::gasPhaseIdx][i], reltol);
BOOST_CHECK_CLOSE(pressures_go[FluidSystem::oilPhaseIdx][i], pressures_gw[FluidSystem::waterPhaseIdx][i], reltol);
BOOST_CHECK_CLOSE(sats_go[FluidSystem::gasPhaseIdx][i], sats_gw[FluidSystem::gasPhaseIdx][i], reltol);
BOOST_CHECK_CLOSE(sats_go[FluidSystem::oilPhaseIdx][i], sats_gw[FluidSystem::waterPhaseIdx][i], reltol);
}
}
BOOST_AUTO_TEST_CASE(DeckWithWetGas)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_wetgas.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(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-1;
const double reltol_ecl = 100.0;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.48215e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.49115e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last], 1.54901e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.482150311e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.547988347e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.491150311e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][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;
// }
std::vector<double> s_ecl[3]; // eclipse
s_ecl[FluidSystem::waterPhaseIdx] = { 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 };
s_ecl[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18311, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_ecl[FluidSystem::gasPhaseIdx] = { 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 };
std::vector<double> s_opm[3]; // opm
s_opm[FluidSystem::waterPhaseIdx] = { 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 };
s_opm[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18288667, 0.0846, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_opm[FluidSystem::gasPhaseIdx] = { 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) {
BOOST_REQUIRE_EQUAL(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';
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol);
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
}
std::cout << std::endl;
}
const auto& rv = comp.rv();
const std::vector<double> 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<double> 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) {
BOOST_CHECK_CLOSE(rv[i], rv_opm[i], reltol);
BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_CASE(DeckWithHumidWetGas)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
FluidSystem::setEnableVaporizedWater(true);
auto simulator = initSimulator<TypeTag>("equil_humidwetgas.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
const double reltol = 1.0e-1;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.480599988e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.549297524e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.489599988e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last], 1.550297524e7, reltol);
const auto& sats = comp.saturation();
std::vector<double> s_opm[3];
s_opm[FluidSystem::waterPhaseIdx] = { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.35838026, 0.64069098, 0.9154626, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
s_opm[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.02738364, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_opm[FluidSystem::gasPhaseIdx] = { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.64161973, 0.359309012, 0.057153701, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE_EQUAL(sats[phase].size(), s_opm[phase].size());
for (size_t i = 0; i < s_opm[phase].size(); ++i) {
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol);
}
std::cout << std::endl;
}
const auto& rv = comp.rv();
const std::vector<double> rv_opm {
0.00024837999651755729, 0.00024869285236692635, 0.00024900604366769004, 0.00024931957094322978,
0.00024963343471801471, 0.00024994763551760586, 0.00025026217386865733, 0.00025057705029892072,
0.00025089226533724643, 0.00025120780158539152, 0.00025105, 0.00025105,
0.00025105, 0.00025105, 0.00025105, 0.00025105,
0.00025105, 0.00025105, 0.00025105, 0.00025105};
const auto& rvw = comp.rvw();
const std::vector<double> rvw_opm {
0.00024837999651755729, 0.00024869285236692635, 0.00024900604366769004, 0.00024931957094322978,
0.00024963343471801471, 0.00024994763551760586, 0.00025026217386865733, 0.00025057705029892072,
0.00025089226533724643, 0.00025120780158539152, 0.00025236969680655122, 0.00025384953117447344,
0.00025532939474124625, 0.00025680928750801825, 0.00025828920947593858, 0.00025976916064615645,
0.00026124914101982041, 0.00026272915059807997, 0.0002642091893820838, 0.00026568925737298143};
for (size_t i = 0; i < rv_opm.size(); ++i) {
BOOST_CHECK_CLOSE(rv[i], rv_opm[i], reltol);
BOOST_CHECK_CLOSE(rvw[i], rvw_opm[i], reltol);
}
}
BOOST_AUTO_TEST_CASE(DeckWithRSVDAndRVVD)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_rsvd_and_rvvd.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(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-4;
const double reltol_ecl = 100.0;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.48350e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.54794e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.49250e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last], 1.54894e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.483499660e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.547924516e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.492499660e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][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;
// }
std::vector<double> s_ecl[3]; // eclipse
s_ecl[FluidSystem::waterPhaseIdx] = { 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 };
s_ecl[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19656529, 0.081805572, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_ecl[FluidSystem::gasPhaseIdx] = { 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 };
std::vector<double> s_opm[3]; // opm
s_opm[FluidSystem::waterPhaseIdx] = { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2223045711692897, 0.52882298575945874, 0.78152142505479982, 0.91816512259416283, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
s_opm[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19637607881498206, 0.08183487740583717, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_opm[FluidSystem::gasPhaseIdx] = { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.7776954288307103, 0.47117701424054126, 0.02210249613021811, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE_EQUAL(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';
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
}
std::cout << std::endl;
}
const auto& rs = comp.rs();
const std::vector<double> 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<double> 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<double> 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<double> 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';
BOOST_CHECK_CLOSE(rs[i], rs_opm[i], reltol);
BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl);
BOOST_CHECK_CLOSE(rv[i], rv_opm[i], reltol);
BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_CASE(DeckWithPBVDAndPDVD)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
auto simulator = initSimulator<TypeTag>("equil_pbvd_and_pdvd.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(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 = 5.0e-4;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 14821552.0, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 15479828.0, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 14911552.0, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last], 15489828.0, 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;
// }
std::vector<double> s_opm[3]; // opm
s_opm[FluidSystem::waterPhaseIdx] = { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24257337312592703, 0.53834824764362788, 0.7844998821510003, 0.9152832369551807, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
s_opm[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18185970596719522, 0.084716763044819343, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_opm[FluidSystem::gasPhaseIdx] = { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75742662687407303, 0.46165175235637212, 0.033640411881804465, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE_EQUAL(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';
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
}
}
const auto& rs = comp.rs();
const std::vector<double> rs_opm { // opm
74.55776480956456,
74.6008507125663,
74.6439680789467,
74.68711693934459,
74.73029732443825,
74.77350926494491,
74.81675279162118,
74.86802321984302,
74.96677993174352,
75.09034523640406,
75, 75, 75,75,75, 75, 75, 75, 75, 75 };
const auto& rv = comp.rv();
const std::vector<double> rv_opm {
0.0002488465888573874,
0.0002491051042753978,
0.0002493638084736803,
0.0002496227016360676,
0.0002498817839466295,
0.00025,
0.00025,
0.00025,
0.00025,
0.000251180039180951,
0.0002522295187440788,
0.0002275000000000001,
0.0002125,
0.0001975,
0.0001825,
0.0001675,
0.0001525,
0.0001375,
0.0001225,
0.0001075};
for (size_t i = 0; i < rv_opm.size(); ++i) {
BOOST_CHECK_CLOSE(rs[i], rs_opm[i], reltol);
BOOST_CHECK_CLOSE(rv[i], rv_opm[i], reltol);
}
}
BOOST_AUTO_TEST_CASE(DeckWithRSVDAndRVVDAndRVWVD)
{
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
using FluidSystem = Opm::GetPropType<TypeTag, Opm::Properties::FluidSystem>;
FluidSystem::setEnableVaporizedWater(true);
auto simulator = initSimulator<TypeTag>("equil_rsvd_and_rvvd_and_rvwvd.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
EquilFixture::Initializer comp(*simulator->problem().materialLawManager(),
eclipseState,
simulator->vanguard().grid(),
simulator->vanguard().gridView(),
simulator->vanguard().cartesianMapper(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
const double reltol = 1.0e-4;
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][first], 1.483359963e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::waterPhaseIdx][last], 1.549297524e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][first], 1.492359963e7, reltol);
BOOST_CHECK_CLOSE(pressures[FluidSystem::oilPhaseIdx][last], 1.550297524e7, reltol);
const auto& sats = comp.saturation();
std::vector<double> s_opm[3];
s_opm[FluidSystem::waterPhaseIdx] = { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.32527877578987319, 0.62976875867666171, 0.918795223850500588, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
s_opm[FluidSystem::oilPhaseIdx] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.054786199472198836, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
s_opm[FluidSystem::gasPhaseIdx] = { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.674721224210102681, 0.37023124132333829, 0.026418562022795279, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE_EQUAL(sats[phase].size(), s_opm[phase].size());
for (size_t i = 0; i < s_opm[phase].size(); ++i) {
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
}
std::cout << std::endl;
}
const auto& rs = comp.rs();
const std::vector<double> rs_opm { // opm
74.617998198796087,74.652774471604374, 74.687905898686935, 74.723393674854691,
74.759238999357947, 74.795443075905553, 74.832007112684167, 74.892422092838459,
74.986801564438935, 75.088917653469338, 52.5, 57.5,
62.5, 67.5, 72.5, 76.528193441026076,
76.774856836636729, 77.021525099679991, 77.268198230347295, 77.514876228830232};
const auto& rv = comp.rv();
const std::vector<double> rv_opm {
2.5000000000000002e-06, 7.5000000000000002e-06, 1.2500000000000001e-05, 1.7500000000000002e-05,
2.2500000000000001e-05, 2.7500000000000004e-05, 3.2500000000000004e-05, 3.7500000000000003e-05,
4.2500000000000003e-05, 0.00025116322680309166, 5.2500000000000002e-05, 5.7500000000000002e-05,
6.2500000000000001e-05, 6.7500000000000001e-05, 7.25e-05, 7.75e-05,
8.25e-05, 8.7500000000000013e-05, 9.2500000000000012e-05, 9.7499999999999998e-05};
const auto& rvw = comp.rvw();
const std::vector<double> rvw_opm {
0.00024920798919277656, 0.00024941664682962629, 0.00024962743539212165, 0.00024984036204912818,
0.00025005543399614773, 0.00025027265845543336, 0.000250492042676105, 0.00025071359393426718,
0.00025093731953312241, 0.00025116322680309166, 0.00025236969680655122, 0.00025384953117447344,
0.00025532939474124625, 0.00025680928750801825, 0.00025828920947593858, 0.00025976916064615645,
0.00026124914101982041, 0.00026272915059807997, 0.0002642091893820838, 0.00026568925737298143};
for (size_t i = 0; i < rv_opm.size(); ++i) {
BOOST_CHECK_CLOSE(rs[i], rs_opm[i], reltol);
BOOST_CHECK_CLOSE(rv[i], rv_opm[i], reltol);
BOOST_CHECK_CLOSE(rvw[i], rvw_opm[i], reltol);
}
}
BOOST_AUTO_TEST_CASE(DeckWithSwatinit)
{
#if 0
using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;
auto simulator = initSimulator<TypeTag>("equil_capillary_swatinit.DATA");
const auto& eclipseState = simulator->vanguard().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
// Create material law manager.
std::vector<int> 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<double> 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<double> 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<double,
/*numPhases=*/3,
/*numComponents=*/3,
FluidSystem,
/*storePressure=*/false,
/*storeTemperature=*/false,
/*storeComposition=*/false,
/*storeFugacity=*/false,
/*storeSaturation=*/true,
/*storeDensity=*/false,
/*storeViscosity=*/false,
/*storeEnthalpy=*/false> 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<double> pc_original(numCells * FluidSystem::numPhases);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = s[FluidSystem::waterPhaseIdx][c];
double so = s[FluidSystem::oilPhaseIdx][c];
double sg = s[FluidSystem::gasPhaseIdx][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<double> 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
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> compScaled(materialLawManagerScaled, eclipseState, simulator->vanguard().gridView(), 9.81, true);
// don't apply swatinit
Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> compUnscaled(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.81, false);
// compute pc
std::vector<double> pc_scaled(numCells * FluidSystem::numPhases);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = compScaled.saturation().data()[FluidSystem::waterPhaseIdx][c];
double so = compScaled.saturation().data()[FluidSystem::oilPhaseIdx][c];
double sg = compScaled.saturation().data()[FluidSystem::gasPhaseIdx][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<double> pc_unscaled(numCells * FluidSystem::numPhases);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = compUnscaled.saturation().data()[FluidSystem::waterPhaseIdx][c];
double so = compUnscaled.saturation().data()[FluidSystem::oilPhaseIdx][c];
double sg = compUnscaled.saturation().data()[FluidSystem::gasPhaseIdx][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-1;
for (int phase = 0; phase < 3; ++phase) {
for (size_t i = 0; i < 20; ++i) {
BOOST_CHECK_CLOSE( pc_original[3*i + phase ], pc_unscaled[3*i + phase ], reltol);
BOOST_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) {
BOOST_CHECK_CLOSE(compUnscaled.saturation()[phase][i], s[phase][i], reltol);
BOOST_CHECK_CLOSE(compScaled.saturation()[phase][i], swatinit[phase][i], reltol);
}
}
#endif
}
Reference in New Issue View Git Blame Copy Permalink