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

restore the hydrostatic equilibration test from opm-core

Browse Source 
This involved quite a bit of kicking and screaming. The result
certainly is not pretty, but it works.
This commit is contained in:
Andreas Lauser
2018-01-02 14:28:52 +01:00
parent c2e9a7a518
commit adb2783a0c
11 changed files with 482 additions and 736 deletions
Download Patch File Download Diff File Expand all files Collapse all files

174
examples/compute_initial_state.cpp
View File

@@ -1,174 +0,0 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2017 IRIS
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/>.
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <boost/filesystem.hpp>
#include <fstream>
namespace
{
void warnIfUnusedParams(const Opm::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
void outputData(const std::string& output_dir,
const std::string& name,
const std::vector<double>& data)
{
std::ostringstream fname;
fname << output_dir << "/" << name;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
fname << "/" << "initial.txt";
std::ofstream file(fname.str().c_str());
if (!file) {
OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
}
std::copy(data.begin(), data.end(), std::ostream_iterator<double>(file, "\n"));
}
/// Convert saturations from a vector of individual phase saturation vectors
/// to an interleaved format where all values for a given cell come before all
/// values for the next cell, all in a single vector.
template <class FluidSystem>
void convertSats(std::vector<double>& sat_interleaved, const std::vector< std::vector<double> >& sat, const Opm::PhaseUsage& pu)
{
assert(sat.size() == 3);
const auto nc = sat[0].size();
const auto np = sat_interleaved.size() / nc;
for (size_t c = 0; c < nc; ++c) {
if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int opos = pu.phase_pos[Opm::BlackoilPhases::Liquid];
const std::vector<double>& sat_p = sat[ FluidSystem::oilPhaseIdx];
sat_interleaved[np*c + opos] = sat_p[c];
}
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int wpos = pu.phase_pos[Opm::BlackoilPhases::Aqua];
const std::vector<double>& sat_p = sat[ FluidSystem::waterPhaseIdx];
sat_interleaved[np*c + wpos] = sat_p[c];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gpos = pu.phase_pos[Opm::BlackoilPhases::Vapour];
const std::vector<double>& sat_p = sat[ FluidSystem::gasPhaseIdx];
sat_interleaved[np*c + gpos] = sat_p[c];
}
}
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
using namespace Opm;
// Setup.
ParameterGroup param(argc, argv);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
const std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParseContext parseContext;
Opm::Parser parser;
const Opm::Deck& deck = parser.parseFile(deck_filename , parseContext);
const Opm::EclipseState eclipseState(deck, parseContext);
const double grav = param.getDefault("gravity", unit::gravity);
GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *gm.c_grid();
warnIfUnusedParams(param);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
typedef FluidSystems::BlackOil<double> FluidSystem;
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
// Initialisation.
//initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
FluidSystem::initFromDeck(deck, eclipseState);
PhaseUsage pu = phaseUsageFromDeck(deck);
typedef EQUIL::DeckDependent::InitialStateComputer<FluidSystem> ISC;
ISC isc(materialLawManager, eclipseState, grid, grav);
const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
const int oilpos = FluidSystem::oilPhaseIdx;
const int waterpos = FluidSystem::waterPhaseIdx;
const int ref_phase = oil ? oilpos : waterpos;
state.pressure() = isc.press()[ref_phase];
convertSats<FluidSystem>(state.saturation(), isc.saturation(), pu);
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
// Output.
const std::string output_dir = param.getDefault<std::string>("output_dir", "output");
outputData(output_dir, "pressure", state.pressure());
outputData(output_dir, "saturation", state.saturation());
outputData(output_dir, "rs", state.gasoilratio());
outputData(output_dir, "rv", state.rv());
}
catch (const std::exception& e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

115
tests/data/equil_base.DATA Normal file
View File

@@ -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

16
tests/capillary.DATA → tests/data/equil_capillary.DATA
View File

@@ -11,7 +11,7 @@ OIL
GAS
DIMENS
40 40 40 /
1 1 20 /
TABDIMS
1 1 40 20 1 20 /
@@ -27,16 +27,22 @@ GRID
-- specify a fake one...
DXV
40*1 /
1*1 /
DYV
40*1 /
1*1 /
DZV
40*1 /
20*5 /
DEPTHZ
1681*123.456 /
4*0.0 /
PORO
20*0.3 /
PERMX
20*500 /
-------------------------------------
PROPS

0
tests/capillary_overlap.DATA → tests/data/equil_capillary_overlap.DATA
View File

0
tests/capillarySwatinit.DATA → tests/data/equil_capillary_swatinit.DATA
View File

23
tests/deadfluids.DATA → tests/data/equil_deadfluids.DATA
View File

@@ -7,11 +7,13 @@
RUNSPEC
WATER
OIL
GAS
OIL
METRIC
DIMENS
3 3 3 /
1 1 10 /
TABDIMS
1 1 40 20 1 20 /
@@ -20,23 +22,28 @@ EQLDIMS
-- NTEQUL
1 /
-------------------------------------
GRID
-- Opm::EclipseState assumes that _some_ grid gets defined, so let's
-- specify a fake one...
DXV
1 2 3 /
1*1 /
DYV
4 5 6 /
1*1 /
DZV
7 8 9 /
10*1 /
DEPTHZ
16*123.456 /
TOPS
1*0.0 /
PORO
10*0.3 /
PERMX
10*500 /
-------------------------------------
PROPS

0
tests/equil_livegas.DATA → tests/data/equil_livegas.DATA
View File

0
tests/equil_liveoil.DATA → tests/data/equil_liveoil.DATA
View File

0
tests/equil_rsvd_and_rvvd.DATA → tests/data/equil_rsvd_and_rvvd.DATA
View File

661
tests/test_equil.cpp → tests/test_equil.cc
View File

@@ -1,45 +1,38 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
*/
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"
/* --- Boost.Test boilerplate --- */
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE UnitsTest
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
/* --- our own headers --- */
#include <opm/core/simulator/initStateEquil.hpp>
#include <ebos/equil/equilibrationhelpers.hh>
#include <ebos/eclproblem.hh>
#include <ewoms/common/start.hh>
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Deck/DeckItem.hpp>
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
#include <opm/parser/eclipse/Units/Dimension.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <dune/common/parallel/mpihelper.hh>
#include <array>
#include <iostream>
#include <limits>
@@ -48,28 +41,56 @@
#include <sstream>
#include <string>
#include <vector>
#include <string.h>
#define CHECK(value, expected) \
{ \
if ((value) != (expected)) \
std::abort(); \
}
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
#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(); \
}
#define CHECK(value, expected, reltol) \
{ \
if (std::fabs((expected)) < 1.e-14) \
BOOST_CHECK_SMALL((value), (reltol)); \
else \
BOOST_CHECK_CLOSE((value), (expected), (reltol)); \
namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(TestEquilTypeTag, INHERITS_FROM(BlackOilModel, EclBaseProblem));
}}
template <class TypeTag>
std::unique_ptr<typename GET_PROP_TYPE(TypeTag, Simulator)>
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_<TypeTag>(/*argc=*/sizeof(argv)/sizeof(argv[0]), argv, /*registerParams=*/false);
return std::unique_ptr<Simulator>(new Simulator);
}
namespace
template <class TypeTag>
static void initDefaultFluidSystem()
{
static void initDefaultFluidSystem() {
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
std::vector<std::pair<double, double> > Bo = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
@@ -130,13 +151,11 @@ static void initDefaultFluidSystem() {
FluidSystem::setWaterPvt(std::move(waterPvt));
FluidSystem::initEnd();
}
}
BOOST_AUTO_TEST_SUITE ()
static Opm::EquilRecord mkEquilRecord( double datd, double datp,
double zwoc, double pcow_woc,
double zgoc, double pcgo_goc ) {
double zgoc, double pcgo_goc )
{
using namespace Opm;
DeckItem dd( "datdep", double() );
@@ -190,71 +209,74 @@ static Opm::EquilRecord mkEquilRecord( double datd, double datp,
return EquilRecord( rec );
}
BOOST_AUTO_TEST_CASE (PhasePressure)
void test_PhasePressure()
{
typedef std::vector<double> PVal;
typedef std::vector<PVal> PPress;
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
initDefaultFluidSystem();
typedef TTAG(TestEquilTypeTag) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
auto simulator = initSimulator<TypeTag>("data/equil_base.DATA");
initDefaultFluidSystem<TypeTag>();
Opm::EQUIL::EquilReg
Ewoms::EQUIL::EquilReg
region(record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0);
std::vector<int> cells(G->number_of_cells);
std::vector<int> cells(simulator->gridManager().grid().size(0));
std::iota(cells.begin(), cells.end(), 0);
const double grav = 10;
const PPress ppress = Opm::EQUIL::phasePressures<FluidSystem>(*G, region, cells, grav);
const PPress ppress = Ewoms::EQUIL::phasePressures<FluidSystem>(simulator->gridManager().grid(), region, cells, grav);
const int first = 0, last = G->number_of_cells - 1;
const int first = 0, last = simulator->gridManager().grid().size(0) - 1;
const double reltol = 1.0e-8;
BOOST_CHECK_CLOSE(ppress[0][first] , 90e3 , reltol);
BOOST_CHECK_CLOSE(ppress[0][last ] , 180e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
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);
}
BOOST_AUTO_TEST_CASE (CellSubset)
void test_CellSubset()
{
typedef std::vector<double> PVal;
typedef std::vector<PVal> PPress;
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
initDefaultFluidSystem();
typedef TTAG(TestEquilTypeTag) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
auto simulator = initSimulator<TypeTag>("data/equil_base.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
initDefaultFluidSystem<TypeTag>();
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg region[] =
Ewoms::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[0],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[0],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[1],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[1],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
};
@@ -263,22 +285,22 @@ BOOST_AUTO_TEST_CASE (CellSubset)
for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
std::vector< std::vector<int> > cells(ncoarse);
for (int c = 0; c < G->number_of_cells; ++c) {
for (int c = 0; c < simulator->gridManager().grid().size(0); ++c) {
int ci = c;
const int i = ci % G->cartdims[0]; ci /= G->cartdims[0];
const int j = ci % G->cartdims[1];
const int k = ci / G->cartdims[1];
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 / (G->cartdims[0] / cdim[0]));
const int jc = (j / (G->cartdims[1] / cdim[1]));
const int kc = (k / (G->cartdims[2] / cdim[2]));
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(G->number_of_cells, 0));
PPress ppress(2, PVal(simulator->gridManager().grid().size(0), 0));
for (std::vector< std::vector<int> >::const_iterator
r = cells.begin(), e = cells.end();
r != e; ++r)
@@ -286,7 +308,7 @@ BOOST_AUTO_TEST_CASE (CellSubset)
const int rno = int(r - cells.begin());
const double grav = 10;
const PPress p =
Opm::EQUIL::phasePressures<FluidSystem>(*G, region[rno], *r, grav);
Ewoms::EQUIL::phasePressures<FluidSystem>(simulator->gridManager().grid(), region[rno], *r, grav);
PVal::size_type i = 0;
for (std::vector<int>::const_iterator
@@ -300,54 +322,51 @@ BOOST_AUTO_TEST_CASE (CellSubset)
}
}
const int first = 0, last = G->number_of_cells - 1;
const int first = 0, last = simulator->gridManager().grid().size(0) - 1;
const double reltol = 1.0e-8;
BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol);
BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
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);
}
BOOST_AUTO_TEST_CASE (RegMapping)
void test_RegMapping()
{
typedef std::vector<double> PVal;
typedef std::vector<PVal> PPress;
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
initDefaultFluidSystem();
typedef TTAG(TestEquilTypeTag) TypeTag;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
auto simulator = initSimulator<TypeTag>("data/equil_base.DATA");
initDefaultFluidSystem<TypeTag>();
Opm::EQUIL::EquilReg region[] =
Ewoms::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[0],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[0],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[1],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
Ewoms::EQUIL::EquilReg(record[1],
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Ewoms::EQUIL::Miscibility::NoMixing>(),
0)
};
std::vector<int> eqlnum(G->number_of_cells);
std::vector<int> eqlnum(simulator->gridManager().grid().size(0));
// [ 0 1; 2 3]
{
for (int i = 0; i < 5; ++i) {
@@ -368,16 +387,16 @@ BOOST_AUTO_TEST_CASE (RegMapping)
}
}
Opm::RegionMapping<> eqlmap(eqlnum);
Ewoms::RegionMapping<> eqlmap(eqlnum);
PPress ppress(2, PVal(G->number_of_cells, 0));
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 =
Opm::EQUIL::phasePressures<FluidSystem>(*G, region[rno], rng, grav);
Ewoms::EQUIL::phasePressures<FluidSystem>(simulator->gridManager().grid(), region[rno], rng, grav);
PVal::size_type i = 0;
for (const auto& c : rng) {
@@ -390,76 +409,47 @@ BOOST_AUTO_TEST_CASE (RegMapping)
}
}
const int first = 0, last = G->number_of_cells - 1;
const int first = 0, last = simulator->gridManager().grid().size(0) - 1;
const double reltol = 1.0e-8;
BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol);
BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
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);
}
BOOST_AUTO_TEST_CASE (DeckAllDead)
void test_DeckAllDead()
{
std::shared_ptr<UnstructuredGrid>
grid(create_grid_cart3d(1, 1, 10), destroy_grid);
Opm::ParseContext parseContext;
Opm::Parser parser;
Opm::Deck deck = parser.parseFile("deadfluids.DATA" , parseContext);
Opm::EclipseState eclipseState(deck, parseContext); // Create material law manager.
typedef TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_deadfluids.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid->number_of_cells, grid->global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, *grid, 10.0);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
REQUIRE(pressures.size() == 3);
REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
const int first = 0, last = grid->number_of_cells - 1;
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;
BOOST_CHECK_CLOSE(pressures[0][first] , 1.496329839e7 , reltol);
BOOST_CHECK_CLOSE(pressures[0][last ] , 1.504526940e7 , reltol);
BOOST_CHECK_CLOSE(pressures[1][last] , 1.504526940e7 , reltol);
CHECK_CLOSE(pressures[0][first] , 1.496329839e7 , reltol);
CHECK_CLOSE(pressures[0][last ] , 1.504526940e7 , reltol);
CHECK_CLOSE(pressures[1][last] , 1.504526940e7 , reltol);
}
BOOST_AUTO_TEST_CASE (CapillaryInversion)
void test_CapillaryInversion()
{
// Test setup.
Opm::GridManager gm(1, 1, 40, 1.0, 1.0, 2.5);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::Parser parser;
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
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<TypeTag>("data/equil_capillary.DATA");
// 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);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
typedef MaterialLawManager::MaterialLaw MaterialLaw;
if (!FluidSystem::isInitialized()) {
// make sure that we don't initialize the fluid system twice
FluidSystem::initFromDeck(deck, eclipseState);
}
// Test the capillary inversion for oil-water.
const int cell = 0;
const double reltol = 1.0e-7;
@@ -468,10 +458,10 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
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(pc.size() == s.size());
REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
const double s_computed = Ewoms::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing);
CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@@ -481,10 +471,10 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
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(pc.size() == s.size());
REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
const double s_computed = Ewoms::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(*simulator->problem().materialLawManager(), phase, cell, pc[i], increasing);
CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@@ -494,41 +484,27 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
const int gas = 2;
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(pc.size() == s.size());
REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromSumOfPcs<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, water, gas, cell, pc[i]);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
const double s_computed = Ewoms::EQUIL::satFromSumOfPcs<FluidSystem, MaterialLaw, MaterialLawManager>(*simulator->problem().materialLawManager(), water, gas, cell, pc[i]);
CHECK_CLOSE(s_computed, s[i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE (DeckWithCapillary)
void test_DeckWithCapillary()
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
typedef typename TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_capillary.DATA");
auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::Parser parser;
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
// 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);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 10.0);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
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,
@@ -536,9 +512,9 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-6;
BOOST_CHECK_CLOSE(pressures[0][first] , 1.469769063e7 , reltol);
BOOST_CHECK_CLOSE(pressures[0][last ] , 15452880.328284413 , reltol);
BOOST_CHECK_CLOSE(pressures[1][last] , 15462880.328284413 , reltol);
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<double> s[3]{
@@ -547,39 +523,25 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
{ 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(sats[phase].size() == s[phase].size());
REQUIRE(sats[phase].size() == s[phase].size());
for (size_t i = 0; i < s[phase].size(); ++i) {
CHECK(sats[phase][i], s[phase][i], reltol);
CHECK_CLOSE(sats[phase][i], s[phase][i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
void test_DeckWithCapillaryOverlap()
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
typedef typename TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_capillary_overlap.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::Parser parser;
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary_overlap.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
// 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);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
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,
@@ -588,14 +550,14 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
// the true answer or something else.
const double reltol = 1.0e-6;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl);
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);
BOOST_CHECK_CLOSE(pressures[0][first] , 14832467.14, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last ] , 15479883.47, reltol);
BOOST_CHECK_CLOSE(pressures[1][last ] , 15489883.47, reltol);
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";
@@ -618,39 +580,28 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
{ 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(sats[phase].size() == s_opm[phase].size());
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(sats[phase][i], s_ecl[phase][i], reltol_ecl);
CHECK(sats[phase][i], s_opm[phase][i], reltol);
CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
void test_DeckWithLiveOil()
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
typedef typename TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_liveoil.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::Parser parser;
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_liveoil.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
// 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);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
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,
@@ -659,15 +610,15 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
// the true answer or something else.
const double reltol = 1.0e-6;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl);
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);
BOOST_CHECK_CLOSE(pressures[0][first], 1.483246714e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last], 1.547991652e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][first], 1.492246714e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][last], 1.548991652e7, reltol);
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";
@@ -688,11 +639,11 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
{ 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(sats[phase].size() == s_opm[phase].size());
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(sats[phase][i], s_opm[phase][i], reltol);
CHECK(sats[phase][i], s_ecl[phase][i], reltol_ecl);
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;
}
@@ -710,36 +661,23 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
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);
CHECK_CLOSE(rs[i], rs_opm[i], reltol);
CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
void test_DeckWithLiveGas()
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
typedef typename TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_livegas.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::Parser parser;
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_livegas.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
// 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);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
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,
@@ -748,15 +686,15 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
// the true answer or something else.
const double reltol = 1.0e-3;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48215e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49115e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl);
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);
BOOST_CHECK_CLOSE(pressures[0][first], 1.482150311e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last], 1.547988347e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][first], 1.491150311e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][last], 1.548988347e7, reltol);
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";
@@ -778,11 +716,11 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
{ 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(sats[phase].size() == s_opm[phase].size());
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(sats[phase][i], s_opm[phase][i], 100.*reltol);
CHECK(sats[phase][i], s_ecl[phase][i], reltol_ecl);
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;
}
@@ -803,34 +741,23 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
0.25104999E-03, 0.25104999E-03, 0.25104999E-03, 0.25104999E-03};
for (size_t i = 0; i < rv_opm.size(); ++i) {
CHECK(rv[i], rv_opm[i], reltol);
CHECK(rv[i], rv_ecl[i], reltol_ecl);
CHECK_CLOSE(rv[i], rv_opm[i], reltol);
CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
void test_DeckWithRSVDAndRVVD()
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
typedef typename TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_rsvd_and_rvvd.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParseContext parseContext;
Opm::Parser parser;
Opm::Deck deck = parser.parseFile("equil_rsvd_and_rvvd.DATA", parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
// 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);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
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,
@@ -839,15 +766,15 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
// the true answer or something else.
const double reltol = 1.0e-6;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48350e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54794e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49250e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54894e+07, reltol_ecl);
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);
BOOST_CHECK_CLOSE(pressures[0][first], 1.483499660e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last], 1.547924516e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][first], 1.492499660e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][last], 1.548924516e7, reltol);
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";
@@ -870,11 +797,11 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
};
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
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(sats[phase][i], s_opm[phase][i], reltol);
CHECK(sats[phase][i], s_ecl[phase][i], reltol_ecl);
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;
}
@@ -911,20 +838,19 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
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);
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);
}
}
BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
void test_DeckWithSwatinit()
{
//Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
Opm::Parser parser;
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillarySwatinit.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
#if 0
typedef typename TTAG(TestEquilTypeTag) TypeTag;
auto simulator = initSimulator<TypeTag>("data/equil_capillary_swatinit.DATA");
const auto& eclipseState = simulator->gridManager().eclState();
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
@@ -1008,9 +934,9 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
// compute the initial state
// apply swatinit
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> compScaled(materialLawManagerScaled, eclipseState, grid, 9.81, true);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> compScaled(materialLawManagerScaled, eclipseState, simulator->gridManager().grid(), 9.81, true);
// don't apply swatinit
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> compUnscaled(materialLawManager, eclipseState, grid, 9.81, false);
Ewoms::EQUIL::DeckDependent::InitialStateComputer<TypeTag> compUnscaled(*simulator->problem().materialLawManager(), eclipseState, simulator->gridManager().grid(), 9.81, false);
// compute pc
std::vector<double> pc_scaled(numCells * FluidSystem::numPhases);
@@ -1051,17 +977,38 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
const double reltol = 1.0e-3;
for (int phase = 0; phase < 3; ++phase) {
for (size_t i = 0; i < 20; ++i) {
CHECK( pc_original[3*i + phase ], pc_unscaled[3*i + phase ], reltol);
CHECK( pc_scaled_truth[3*i + phase], pc_scaled[3*i + phase ], reltol);
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(compUnscaled.saturation()[phase][i], s[phase][i], reltol);
CHECK(compScaled.saturation()[phase][i], swatinit[phase][i], reltol);
CHECK_CLOSE(compUnscaled.saturation()[phase][i], s[phase][i], reltol);
CHECK_CLOSE(compScaled.saturation()[phase][i], swatinit[phase][i], reltol);
}
}
#endif
}
BOOST_AUTO_TEST_SUITE_END()
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
typedef TTAG(TestEquilTypeTag) TypeTag;
Ewoms::registerAllParameters_<TypeTag>();
test_PhasePressure();
test_CellSubset();
test_RegMapping();
test_DeckAllDead();
test_CapillaryInversion();
test_DeckWithCapillary();
test_DeckWithCapillaryOverlap();
test_DeckWithLiveOil();
test_DeckWithLiveGas();
test_DeckWithRSVDAndRVVD();
//test_DeckWithSwatinit();
return 0;
}

155
tests/test_regionmapping.cpp
View File

@@ -1,155 +0,0 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
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/>.
*/
#include "config.h"
/* --- Boost.Test boilerplate --- */
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE RegionMapping
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
/* --- our own headers --- */
#include <opm/core/utility/RegionMapping.hpp>
#include <algorithm>
#include <map>
BOOST_AUTO_TEST_SUITE ()
BOOST_AUTO_TEST_CASE (Forward)
{
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
for (decltype(regions.size()) i = 0, n = regions.size(); i < n; ++i) {
BOOST_CHECK_EQUAL(rm.region(i), regions[i]);
}
}
BOOST_AUTO_TEST_CASE (ActiveRegions)
{
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
auto active = [&region_ids](const int reg)
{
auto b = region_ids.begin();
auto e = region_ids.end();
return std::find(b, e, reg) != e;
};
BOOST_CHECK_EQUAL(rm.activeRegions().size(), region_ids.size());
for (const auto& reg : rm.activeRegions()) {
BOOST_CHECK(active(reg));
}
}
BOOST_AUTO_TEST_CASE (Consecutive)
{
using RegionCells = std::map<int, std::vector<int>>;
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
RegionCells region_cells;
{
using VT = RegionCells::value_type;
region_cells.insert(VT(2, { 0, 2, 4 }));
region_cells.insert(VT(3, { 7 }));
region_cells.insert(VT(4, { 3 }));
region_cells.insert(VT(5, { 1 }));
region_cells.insert(VT(6, { 6, 8 }));
region_cells.insert(VT(7, { 5 }));
}
for (const auto& reg : region_ids) {
const auto& cells = rm.cells(reg);
const auto& expect = region_cells[reg];
BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(),
expect.begin(), expect.end());
}
// Verify that there are no cells in unused regions 0 and 1.
for (const auto& r : { 0, 1 }) {
BOOST_CHECK(rm.cells(r).empty());
}
}
BOOST_AUTO_TEST_CASE (NonConsecutive)
{
using RegionCells = std::map<int, std::vector<int>>;
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 4, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
std::vector<int> region_ids = { 2, 3, 4, 6, 7 };
RegionCells region_cells;
{
using VT = RegionCells::value_type;
region_cells.insert(VT(2, { 0, 2, 4 }));
region_cells.insert(VT(3, { 7 }));
region_cells.insert(VT(4, { 1, 3 }));
region_cells.insert(VT(6, { 6, 8 }));
region_cells.insert(VT(7, { 5 }));
}
for (const auto& reg : region_ids) {
const auto& cells = rm.cells(reg);
const auto& expect = region_cells[reg];
BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(),
expect.begin(), expect.end());
}
// Verify that there are no cells in unused regions 0, 1, and 5.
for (const auto& r : { 0, 1, 5 }) {
BOOST_CHECK(rm.cells(r).empty());
}
}
BOOST_AUTO_TEST_SUITE_END()