OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #646 from osae/satFunc-eclState-swl

Browse Source 
Initialization of SWL-family from EclipseState.  (Needs PR opm-parser/#278.)
This commit is contained in:
Atgeirr Flø Rasmussen 2014-11-06 10:51:51 +01:00
commit 665f08592c
12 changed files with 1909 additions and 392 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists_files.cmake
View File

@ -172,6 +172,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_parallel_linearsolver.cpp
tests/test_param.cpp
tests/test_blackoilfluid.cpp
tests/test_satfunc.cpp
tests/test_shadow.cpp
tests/test_equil.cpp
tests/test_regionmapping.cpp
@ -200,6 +201,12 @@ list (APPEND TEST_DATA_FILES
tests/equil_liveoil.DATA
tests/equil_rsvd_and_rvvd.DATA
tests/wetgas.DATA
tests/satfuncStandard.DATA
tests/satfuncEPSBase.DATA
tests/satfuncEPS_A.DATA
tests/satfuncEPS_B.DATA
tests/satfuncEPS_C.DATA
tests/satfuncEPS_D.DATA
tests/testBlackoilState1.DATA
tests/testBlackoilState2.DATA
tests/wells_manager_data.data

20
opm/core/props/satfunc/SatFuncGwseg.hpp
View File

@ -360,8 +360,8 @@ namespace Opm
const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
const double krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
@ -373,8 +373,8 @@ namespace Opm
// Derivatives.
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow));
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
@ -469,21 +469,21 @@ namespace Opm
if (ssow >= sat_hyst->sow_hyst) { // Drainage
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow));
krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
} else { // Imbibition
double ssow_shifted = ssow + sat_hyst->sow_shift;
double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst_hyst->watoil.scaleSatDeriv(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_);
dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(1.0-_ssow));
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_);
dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(std::max(1.0-_ssow,_swco)));
}
// Oil in gas and connate water - use drainage curve only
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_);
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps)));
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;

6
opm/core/props/satfunc/SatFuncSimple.cpp
View File

@ -159,7 +159,7 @@ namespace Opm
double EPSTransforms::Transform::scaleKrDeriv(double s, double krDeriv) const
{
if (doKrCrit) {
if (s <= scr) {
if (s < scr) {
return 0.0;
} else if (s <= sr) {
return krDeriv*krSlopeCrit;
@ -172,7 +172,7 @@ namespace Opm
return 0.0;
}
} else if (doKrMax) {
if (s <= scr) {
if (s < scr) {
return 0.0;
} else if (s <= smax) {
return krDeriv*krSlopeMax;
@ -180,7 +180,7 @@ namespace Opm
return 0.0;
}
} else {
if (s <= scr) {
if (s < scr) {
return 0.0;
} else if (s <= smax) {
return krDeriv;

19
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
View File

@ -60,7 +60,7 @@ namespace Opm
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
Opm::EclipseStateConstPtr eclipseState,
const UnstructuredGrid& grid,
const int samples);
@ -79,7 +79,7 @@ namespace Opm
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
template<class T>
void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
@ -165,14 +165,9 @@ namespace Opm
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSHyst(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
int dimensions,
const std::vector<std::string>& eps_kw,
std::vector<EPSTransforms>& eps_transf);
template<class T>
void initEPSKey(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
@ -204,8 +199,8 @@ namespace Opm
bool columnIsMasked_(Opm::DeckConstPtr deck,
const std::string& keywordName,
int /* columnIdx */)
{ return deck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; }
int columnIdx)
{ return deck->getKeyword(keywordName)->getRecord(columnIdx)->getItem(0)->getSIDouble(0) != -1.0; }
};

492
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
View File

@ -31,6 +31,7 @@
#include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
#include <iostream>
#include <map>
namespace Opm
{
@ -48,11 +49,11 @@ namespace Opm
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
Opm::EclipseStateConstPtr eclipseState,
const UnstructuredGrid& grid,
const int samples)
{
this->init(deck, eclState, grid.number_of_cells,
this->init(deck, eclipseState, grid.number_of_cells,
grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
}
@ -61,7 +62,7 @@ namespace Opm
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
@ -128,7 +129,7 @@ namespace Opm
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(eclState, table, phase_usage_, samples);
satfuncset_[table].init(eclipseState, table, phase_usage_, samples);
}
// Check EHYSTR status
@ -190,8 +191,11 @@ namespace Opm
// TODO: ENPTVD/ENKRVD: Too few tables gives a cryptical message from parser,
// superfluous tables are ignored by the parser without any warning ...
initEPS(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
const std::vector<std::string> eps_kw{"SWL", "SWU", "SWCR", "SGL", "SGU", "SGCR", "SOWCR",
"SOGCR", "KRW", "KRG", "KRO", "KRWR", "KRGR", "KRORW", "KRORG", "PCW", "PCG"};
eps_transf_.resize(number_of_cells);
initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids,
dimensions, eps_kw, eps_transf_);
if (do_hyst_) {
if (deck->hasKeyword("KRW")
@ -230,15 +234,20 @@ namespace Opm
// TODO: Make actual use of IMBNUM. For now we just consider the imbibition curve
// to be a scaled version of the drainage curve (confer Norne model).
}
initEPSHyst(deck, eclState, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
const std::vector<std::string> eps_i_kw{"ISWL", "ISWU", "ISWCR", "ISGL", "ISGU", "ISGCR", "ISOWCR",
"ISOGCR", "IKRW", "IKRG", "IKRO", "IKRWR", "IKRGR", "IKRORW", "IKRORG", "IPCW", "IPCG"};
eps_transf_hyst_.resize(number_of_cells);
sat_hyst_.resize(number_of_cells);
initEPS(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroids,
dimensions, eps_i_kw, eps_transf_hyst_);
}
}
}
/// \return P, the number of phases.
template <class SatFuncSet>
int SaturationPropsFromDeck<SatFuncSet>::numPhases() const
@ -462,204 +471,17 @@ namespace Opm
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions)
int dimensions,
const std::vector<std::string>& eps_kw,
std::vector<EPSTransforms>& eps_transf)
{
std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
std::vector<double> pcw, pcg;
std::vector<std::vector<double> > eps_vec(eps_kw.size());
const std::vector<double> dummy;
// Initialize saturation scaling parameter
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWL"), swl);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWU"), swu);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWCR"), swcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGL"), sgl);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGU"), sgu);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGCR"), sgcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOWCR"), sowcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOGCR"), sogcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRW"), krw);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRG"), krg);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRO"), kro);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRWR"), krwr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRGR"), krgr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORW"), krorw);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORG"), krorg);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("PCW"), pcw);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("PCG"), pcg);
eps_transf_.resize(number_of_cells);
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour];
const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
for (int cell = 0; cell < number_of_cells; ++cell) {
if (oilWater) {
// ### krw
initEPSParam(cell, eps_transf_[cell].wat, false,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_,
-1.0,
funcForCell(cell).krwr_,
funcForCell(cell).krwmax_,
funcForCell(cell).pcwmax_,
swl, swcr, swu, sowcr, sgl, krwr, krw, pcw);
// ### krow
initEPSParam(cell, eps_transf_[cell].watoil, true,
0.0,
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
-1.0,
funcForCell(cell).krorw_,
funcForCell(cell).kromax_,
0.0,
swl, sowcr, swl, swcr, sgl, krorw, kro, dummy);
} else if (oilGas) {
// ### krg
initEPSParam(cell, eps_transf_[cell].gas, false,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_,
-1.0,
funcForCell(cell).krgr_,
funcForCell(cell).krgmax_,
funcForCell(cell).pcgmax_,
sgl, sgcr, sgu, sogcr, swl, krgr, krg, pcg);
// ### krog
initEPSParam(cell, eps_transf_[cell].gasoil, true,
0.0,
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
-1.0,
funcForCell(cell).krorg_,
funcForCell(cell).kromax_,
0.0,
sgl, sogcr, sgl, sgcr, swl, krorg, kro, dummy);
} else if (threephase) {
// ### krw
initEPSParam(cell, eps_transf_[cell].wat, false,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).krwr_,
funcForCell(cell).krwmax_,
funcForCell(cell).pcwmax_,
swl, swcr, swu, sowcr, sgl, krwr, krw, pcw);
// ### krow
initEPSParam(cell, eps_transf_[cell].watoil, true,
0.0,
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).krorw_,
funcForCell(cell).kromax_,
0.0,
swl, sowcr, swl, swcr, sgl, krorw, kro, dummy);
// ### krg
initEPSParam(cell, eps_transf_[cell].gas, false,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).krgr_,
funcForCell(cell).krgmax_,
funcForCell(cell).pcgmax_,
sgl, sgcr, sgu, sogcr, swl, krgr, krg, pcg);
// ### krog
initEPSParam(cell, eps_transf_[cell].gasoil, true,
0.0,
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).krorg_,
funcForCell(cell).kromax_,
0.0,
sgl, sogcr, sgl, sgcr, swl, krorg, kro, dummy);
}
for (size_t i = 0; i < eps_kw.size(); ++i) {
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
eps_kw[i], eps_vec[i]);
}
}
// Initialize hysteresis saturation scaling parameters
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions)
{
std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
std::vector<double> ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg;
std::vector<double> ipcw, ipcg;
const std::vector<double> dummy;
// Initialize hysteresis saturation scaling parameters
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWL"), iswl);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWU"), iswu);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWCR"), iswcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGL"), isgl);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGU"), isgu);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGCR"), isgcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOWCR"), isowcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOGCR"), isogcr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRW"), ikrw);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRG"), ikrg);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRO"), ikro);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRWR"), ikrwr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRGR"), ikrgr);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORW"), ikrorw);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORG"), ikrorg);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IPCW"), ipcw);
initEPSKey(deck, eclipseState, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IPCG"), ipcg);
eps_transf_hyst_.resize(number_of_cells);
sat_hyst_.resize(number_of_cells);
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
@ -668,97 +490,53 @@ namespace Opm
const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
for (int cell = 0; cell < number_of_cells; ++cell) {
if (oilWater) {
// ### krw
initEPSParam(cell, eps_transf_hyst_[cell].wat, false,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_,
-1.0,
funcForCell(cell).krwr_,
funcForCell(cell).krwmax_,
funcForCell(cell).pcwmax_,
iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw, ipcw);
// ### krow
initEPSParam(cell, eps_transf_hyst_[cell].watoil, true,
0.0,
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
-1.0,
funcForCell(cell).krorw_,
funcForCell(cell).kromax_,
0.0,
iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro, dummy);
} else if (oilGas) {
// ### krg
initEPSParam(cell, eps_transf_hyst_[cell].gas, false,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_,
-1.0,
funcForCell(cell).krgr_,
funcForCell(cell).krgmax_,
funcForCell(cell).pcgmax_,
isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg, ipcg);
// ### krog
initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true,
0.0,
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
-1.0,
funcForCell(cell).krorg_,
funcForCell(cell).kromax_,
0.0,
isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro, dummy);
} else if (threephase) {
// ### krw
initEPSParam(cell, eps_transf_hyst_[cell].wat, false,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).krwr_,
funcForCell(cell).krwmax_,
funcForCell(cell).pcwmax_,
iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw, ipcw);
// ### krow
initEPSParam(cell, eps_transf_hyst_[cell].watoil, true,
0.0,
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).krorw_,
funcForCell(cell).kromax_,
0.0,
iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro, dummy);
// ### krg
initEPSParam(cell, eps_transf_hyst_[cell].gas, false,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).krgr_,
funcForCell(cell).krgmax_,
funcForCell(cell).pcgmax_,
isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg, ipcg);
// ### krog
initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true,
0.0,
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).krorg_,
funcForCell(cell).kromax_,
0.0,
isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro, dummy);
if (threephase || oilWater) {
// ### krw
initEPSParam(cell, eps_transf[cell].wat, false,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
funcForCell(cell).smax_[wpos],
funcForCell(cell).sowcr_,
oilWater ? -1.0 : funcForCell(cell).smin_[gpos],
funcForCell(cell).krwr_,
funcForCell(cell).krwmax_,
funcForCell(cell).pcwmax_,
eps_vec[0], eps_vec[2], eps_vec[1], eps_vec[6], eps_vec[3], eps_vec[11], eps_vec[8], eps_vec[15]);
// ### krow
initEPSParam(cell, eps_transf[cell].watoil, true,
0.0,
funcForCell(cell).sowcr_,
funcForCell(cell).smin_[wpos],
funcForCell(cell).swcr_,
oilWater ? -1.0 : funcForCell(cell).smin_[gpos],
funcForCell(cell).krorw_,
funcForCell(cell).kromax_,
0.0,
eps_vec[0], eps_vec[6], eps_vec[0], eps_vec[2], eps_vec[3], eps_vec[13], eps_vec[10], dummy);
}
if (threephase || oilGas) {
// ### krg
initEPSParam(cell, eps_transf[cell].gas, false,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
funcForCell(cell).smax_[gpos],
funcForCell(cell).sogcr_,
oilGas ? -1.0 : funcForCell(cell).smin_[wpos],
funcForCell(cell).krgr_,
funcForCell(cell).krgmax_,
funcForCell(cell).pcgmax_,
eps_vec[3], eps_vec[5], eps_vec[4], eps_vec[7], eps_vec[0], eps_vec[12], eps_vec[9], eps_vec[16]);
// ### krog
initEPSParam(cell, eps_transf[cell].gasoil, true,
0.0,
funcForCell(cell).sogcr_,
funcForCell(cell).smin_[gpos],
funcForCell(cell).sgcr_,
oilGas ? -1.0 : funcForCell(cell).smin_[wpos],
funcForCell(cell).krorg_,
funcForCell(cell).kromax_,
0.0,
eps_vec[3], eps_vec[7], eps_vec[3], eps_vec[5], eps_vec[0], eps_vec[14], eps_vec[10], dummy);
}
}
}
@ -779,6 +557,12 @@ namespace Opm
const bool useLiquid = phase_usage_.phase_used[Liquid];
const bool useVapour = phase_usage_.phase_used[Vapour];
bool useKeyword = deck->hasKeyword(keyword);
bool useStateKeyword = eclipseState->hasDoubleGridProperty(keyword);
const std::map<std::string, int> kw2tab = {
{"SWL", 1}, {"SWCR", 2}, {"SWU", 3}, {"SGL", 4},
{"SGCR", 5}, {"SGU", 6}, {"SOWCR", 7}, {"SOGCR", 8},
{"ISWL", 1}, {"ISWCR", 2}, {"ISWU", 3}, {"ISGL", 4},
{"ISGCR", 5}, {"ISGU", 6}, {"ISOWCR", 7}, {"ISOGCR", 8}};
bool hasENPTVD = deck->hasKeyword("ENPTVD");
bool hasENKRVD = deck->hasKeyword("ENKRVD");
int itab = 0;
@ -787,67 +571,10 @@ namespace Opm
std::vector<std::string> col_names;
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour];
if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) {
if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 0))) {
itab = 1;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
}
} else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 1))) {
itab = 2;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).swcr_;
}
} else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 2))) {
itab = 3;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
}
} else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 3))) {
itab = 4;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
}
} else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 4))) {
itab = 5;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sgcr_;
}
} else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 5))) {
itab = 6;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
}
} else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENPTVD", 6))) {
itab = 7;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sowcr_;
}
} else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 7))) {
itab = 8;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sogcr_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
if ((keyword[0] == 'S' && (useStateKeyword || hasENPTVD)) || (keyword[1] == 'S' && useStateKeyword) ) {
if (useAqua && (useStateKeyword || columnIsMasked_(deck, "ENPTVD", kw2tab.find(keyword)->second-1))) {
itab = kw2tab.find(keyword)->second;
scaleparam.resize(number_of_cells);
}
if (!useKeyword && itab > 0) {
const auto& enptvdTables = eclipseState->getEnptvdTables();
@ -926,7 +653,7 @@ namespace Opm
param_col[table_num] = enkrvdTable.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
}
}
} else if (useKeyword && (keyword[0] == 'P' || keyword[1] == 'P') ) {
} else if (useKeyword && (keyword[0] == 'P' || keyword[1] == 'P') ) {
if (useAqua && (keyword == std::string("PCW") || keyword == std::string("IPCW")) ) {
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
@ -940,16 +667,25 @@ namespace Opm
if (scaleparam.empty()) {
return;
} else if (useKeyword) {
}
if (useKeyword || useStateKeyword) {
// Keyword values from deck
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
const int* gc = global_cell;
const std::vector<double>& val = deck->getKeyword(keyword)->getSIDoubleData();
std::vector<double> val;
if (keyword[0] == 'S' || keyword[1] == 'S') { // Saturation from EclipseState
val = eclipseState->getDoubleGridProperty(keyword)->getData();
} else {
val = deck->getKeyword(keyword)->getSIDoubleData(); //KR and PC directly from deck.
}
for (int c = 0; c < int(scaleparam.size()); ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c];
scaleparam[c] = val[deck_pos];
}
} else {
}
if (itab > 0) {
const int dim = dimensions;
std::vector<int> endnum;
if ( deck->hasKeyword("ENDNUM")) {
@ -966,17 +702,39 @@ namespace Opm
// Default deck value is one
endnum.assign(number_of_cells, 0);
}
for (int cell = 0; cell < number_of_cells; ++cell) {
if (endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
double zc = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
dim-1);
if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
if (keyword[0] == 'S' || keyword[1] == 'S') { // From EclipseState
for (int cell = 0; cell < number_of_cells; ++cell) {
if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
double zc = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
dim-1);
if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
}
} else if (!std::isfinite(scaleparam[cell]) && endnum[cell] >= 0) {
// As of 1/9-2014: Reflects remaining work on opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp ...
OPM_THROW(std::runtime_error, " -- Inconsistent EclipseState: '" << keyword << "' (ENPTVD)");
}
}
} else { //KR and PC from deck.
for (int cell = 0; cell < number_of_cells; ++cell) {
if (endnum[cell] >= 0 && param_col[endnum[cell]][0] >= 0.0) {
double zc = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
dim-1);
if (zc >= depth_col[endnum[cell]].front() && zc <= depth_col[endnum[cell]].back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth_col[endnum[cell]], param_col[endnum[cell]], zc);
}
}
}
}
}
// std::cout << keyword << ":" << std::endl;
// for (int c = 0; c < int(scaleparam.size()); ++c) {
// std::cout << c << " " << scaleparam[c] << std::endl;
// }
}
// Saturation scaling

168
tests/satfuncEPSBase.DATA Normal file
View File

@ -0,0 +1,168 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
ENDSCALE
--DIR REV NTENDP NSENDP
'NODIR' 'REVERS' 1 20 /
/
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
20*5.0
/
PORO
20*0.2
/
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SCALECRS
-- YES /
NO/
SWOF
0.1 0.0 1.0 0.9
0.2 0.0 0.8 0.8
0.3 0.1 0.6 0.7
0.4 0.2 0.4 0.6
0.7 0.5 0.1 0.3
0.8 0.6 0.0 0.2
0.9 0.7 0.0 0.1
/
SGOF
0.0 0.0 1.0 0.2
0.1 0.0 0.7 0.4
0.2 0.1 0.6 0.6
0.8 0.7 0.0 2.0
0.9 1.0 0.0 2.1
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
RSVD
0 0.0
100 100. /
RVVD
0. 0.
100. 0.0001 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

177
tests/satfuncEPS_A.DATA Normal file
View File

@ -0,0 +1,177 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 10
/
WELLDIMS
30 10 2 30 /
ENDSCALE
--DIR REV NTENDP NSENDP
'NODIR' 'REVERS' 1 20 /
/
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
10*5.0
/
PORO
10*0.2
/
PERMZ
10*1.0
/
PERMY
10*100.0
/
PERMX
10*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SCALECRS
-- YES /
NO/
SWL
4*0.1 4*0.2 2*0.1/
SWCR
0.2 0.2 0.4 0.4 0.2 0.2 0.4 0.4 2*0.2 /
SWU
0.9 0.7 0.9 0.7 0.9 0.7 0.9 0.7 2*0.9 /
SWOF
0.1 0.0 1.0 0.9
0.2 0.0 0.8 0.8
0.3 0.1 0.6 0.7
0.4 0.2 0.4 0.6
0.7 0.5 0.1 0.3
0.8 0.6 0.0 0.2
0.9 0.7 0.0 0.1
/
SGOF
0.0 0.0 1.0 0.2
0.1 0.0 0.7 0.4
0.2 0.1 0.6 0.6
0.8 0.7 0.0 2.0
0.9 1.0 0.0 2.1
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
RSVD
0 0.0
100 100. /
RVVD
0. 0.
100. 0.0001 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

190
tests/satfuncEPS_B.DATA Normal file
View File

@ -0,0 +1,190 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 10
/
WELLDIMS
30 10 2 30 /
ENDSCALE
--DIR REV NTENDP NSENDP
'NODIR' 'REVERS' 1 20 /
/
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
10*5.0
/
PORO
10*0.2
/
PERMZ
10*1.0
/
PERMY
10*100.0
/
PERMX
10*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SCALECRS
-- YES /
NO/
ENPTVD
-- Depth swl swcr swu sgl sgcr sgu sowcr sogcr
-- Sandstone
2.5 .1 0.2 0.9 .0 0.1 0.9 0.2 0.2
7.5 .1 0.2 0.7 .0 0.1 0.9 0.2 0.2
12.5 .1 0.4 0.9 .0 0.1 0.9 0.2 0.2
17.5 .1 0.4 0.7 .0 0.1 0.9 0.2 0.2
22.5 .2 0.2 0.9 .0 0.1 0.9 0.2 0.2
27.5 .2 0.2 0.7 .0 0.1 0.9 0.2 0.2
32.5 .2 0.4 0.9 .0 0.1 0.9 0.2 0.2
37.5 .2 0.4 0.7 .0 0.1 0.9 0.2 0.2
42.5 .1 0.2 0.9 .0 0.1 0.9 0.2 0.2
47.5 .1 0.2 0.9 .0 0.1 0.9 0.2 0.2
50.0 .1 0.2 0.9 .0 0.1 0.9 0.2 0.2 /
/
SWOF
0.1 0.0 1.0 0.9
0.2 0.0 0.8 0.8
0.3 0.1 0.6 0.7
0.4 0.2 0.4 0.6
0.7 0.5 0.1 0.3
0.8 0.6 0.0 0.2
0.9 0.7 0.0 0.1
/
SGOF
0.0 0.0 1.0 0.2
0.1 0.0 0.7 0.4
0.2 0.1 0.6 0.6
0.8 0.7 0.0 2.0
0.9 1.0 0.0 2.1
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
REGIONS ======
ENDNUM
10*1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
RSVD
0 0.0
100 100. /
RVVD
0. 0.
100. 0.0001 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

225
tests/satfuncEPS_C.DATA Normal file
View File

@ -0,0 +1,225 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 10
/
WELLDIMS
30 10 2 30 /
ENDSCALE
--DIR REV NTENDP NSENDP
'NODIR' 'REVERS' 1 20 /
/
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
10*5.0
/
PORO
10*0.2
/
PERMZ
10*1.0
/
PERMY
10*100.0
/
PERMX
10*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SCALECRS
-- YES /
NO/
--SWL
--4*0.1 4*0.2 2*0.1/
EQUALS
SWL 0.1 1 1 1 1 1 1 /
SWL 0.1 1 1 1 1 2 2 /
SWL 0.1 1 1 1 1 3 3 /
SWL 0.1 1 1 1 1 4 4 /
SWL 0.2 1 1 1 1 5 5 /
SWL 0.2 1 1 1 1 6 6 /
SWL 0.2 1 1 1 1 7 7 /
SWL 0.2 1 1 1 1 8 8 /
SWL 0.1 1 1 1 1 9 9 /
SWL 0.1 1 1 1 1 10 10 /
/
--SWCR
--0.2 0.2 0.4 0.4 0.2 0.2 0.4 0.4 2*0.2 /
COPY
SWL SWCR /
SWL SWU /
/
ADD
SWCR 0.1 1 1 1 1 1 2 /
SWCR 0.3 1 1 1 1 3 4 /
SWCR 0.2 1 1 1 1 7 8 /
SWCR 0.1 1 1 1 1 9 10 /
/
--SWU
--0.9 0.7 0.9 0.7 0.9 0.7 0.9 0.7 2*0.9 /
MULTIPLY
SWU -1.0 1 1 1 1 1 10 /
/
ADD
SWU 1.0 1 1 1 1 1 10 /
/
ADD
SWU -0.2 1 1 1 1 2 2 /
SWU -0.2 1 1 1 1 4 4 /
SWU 0.1 1 1 1 1 5 5 /
SWU -0.1 1 1 1 1 6 6 /
SWU 0.1 1 1 1 1 7 7 /
SWU -0.1 1 1 1 1 8 8 /
/
SWOF
0.1 0.0 1.0 0.9
0.2 0.0 0.8 0.8
0.3 0.1 0.6 0.7
0.4 0.2 0.4 0.6
0.7 0.5 0.1 0.3
0.8 0.6 0.0 0.2
0.9 0.7 0.0 0.1
/
SGOF
0.0 0.0 1.0 0.2
0.1 0.0 0.7 0.4
0.2 0.1 0.6 0.6
0.8 0.7 0.0 2.0
0.9 1.0 0.0 2.1
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
REGIONS ======
ENDNUM
10*1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
RSVD
0 0.0
100 100. /
RVVD
0. 0.
100. 0.0001 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

181
tests/satfuncEPS_D.DATA Normal file
View File

@ -0,0 +1,181 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 10
/
WELLDIMS
30 10 2 30 /
ENDSCALE
--DIR REV NTENDP NSENDP
'NODIR' 'REVERS' 1 20 /
/
SATOPTS
HYSTER /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
10*5.0
/
PORO
10*0.2
/
PERMZ
10*1.0
/
PERMY
10*100.0
/
PERMX
10*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
-- Hysteresis input
EHYSTR
0.1 0 0.1 1* KR /
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SCALECRS
-- YES /
NO/
SWOF
0.1 0.0 1.0 0.9
0.2 0.0 0.8 0.8
0.3 0.1 0.6 0.7
0.4 0.2 0.4 0.6
0.7 0.5 0.1 0.3
0.8 0.6 0.0 0.2
0.9 0.7 0.0 0.1
/
SGOF
0.0 0.0 1.0 0.2
0.1 0.0 0.7 0.4
0.2 0.1 0.6 0.6
0.8 0.7 0.0 2.0
0.9 1.0 0.0 2.1
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
REGIONS ======
ENDNUM
10*1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
RSVD
0 0.0
100 100. /
RVVD
0. 0.
100. 0.0001 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

168
tests/satfuncStandard.DATA Normal file
View File

@ -0,0 +1,168 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 10
/
WELLDIMS
30 10 2 30 /
--ENDSCALE
--DIR REV NTENDP NSENDP
--'NODIR' 'REVERS' 1 20 /
--/
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
10*5.0
/
PORO
10*0.2
/
PERMZ
10*1.0
/
PERMY
10*100.0
/
PERMX
10*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
--SCALECRS
-- YES /
-- NO/
SWOF
0.1 0.0 1.0 0.9
0.2 0.0 0.8 0.8
0.3 0.1 0.6 0.7
0.4 0.2 0.4 0.6
0.7 0.5 0.1 0.3
0.8 0.6 0.0 0.2
0.9 0.7 0.0 0.1
/
SGOF
0.0 0.0 1.0 0.2
0.1 0.0 0.7 0.4
0.2 0.1 0.6 0.6
0.8 0.7 0.0 2.0
0.9 1.0 0.0 2.1
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
RSVD
0 0.0
100 100. /
RVVD
0. 0.
100. 0.0001 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

648
tests/test_satfunc.cpp Normal file
View File

@ -0,0 +1,648 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
*/
#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 <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/core/pressure/msmfem/partition.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <array>
#include <iostream>
#include <limits>
#include <memory>
#include <numeric>
#include <sstream>
#include <string>
#include <vector>
BOOST_AUTO_TEST_SUITE ()
BOOST_AUTO_TEST_CASE (GwsegStandard)
{
// This is the basic (no eps and hysteris) version of
// the Gwseg model.
//std::cout << "==================================== GwsegStandard ====================================" << std::endl;
Opm::parameter::ParameterGroup param;
Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("satfuncStandard.DATA");
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
const int np = props.numPhases();
const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
BOOST_REQUIRE(np == 3);
BOOST_REQUIRE(wpos == 0);
BOOST_REQUIRE(opos == 1);
BOOST_REQUIRE(gpos == 2);
const int n=11;
double s[n*np];
int cells[n];
double kr[n*np];
double dkrds[n*np*np];
for (int i=0; i<n; ++i) {
cells[i] = 0;
s[i*np+wpos] = i*0.1;
s[i*np+opos] = 1.0-s[i*np+wpos];
s[i*np+gpos] = 0.0;
}
props.relperm(n, s, cells, kr, dkrds);
double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0, -0.66666666666666741, -0.75, -0.8,
-0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
const double reltol = 1.0e-6;
for (int i=0; i<n; ++i) {
BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[i], reltol);
BOOST_CHECK_CLOSE(kr[i*np+opos], kro[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
}
/*
std::cout << std::setw(12) << "sw";
std::cout << std::setw(12) << "so";
std::cout << std::setw(12) << "sg";
std::cout << std::setw(12) << "krw";
std::cout << std::setw(12) << "kro";
std::cout << std::setw(12) << "krg";
std::cout << std::setw(12) << "DkrwDsw";
std::cout << std::setw(12) << "DkroDsw";
std::cout << std::setw(12) << "DkrgDsw";
std::cout << std::setw(12) << "DkrwDso";
std::cout << std::setw(12) << "DkroDso";
std::cout << std::setw(12) << "DkrgDso";
std::cout << std::setw(12) << "DkrwDsg";
std::cout << std::setw(12) << "DkroDsg";
std::cout << std::setw(12) << "DkrgDsg";
std::cout << std::endl;
for (int i=0; i<n; ++i) {
std::cout << std::setw(12) << s[i*np+wpos] << std::setw(12) << s[i*np+opos] << std::setw(12) << s[i*np+gpos]
<< std::setw(12) << kr[i*np+wpos] << std::setw(12) << kr[i*np+opos] << std::setw(12) << kr[i*np+gpos];
for (int j=0; j<np*np; ++j) {
std::cout << std::setw(12) << dkrds[i*np*np+j];
}
std::cout << std::endl;
}
*/
}
BOOST_AUTO_TEST_CASE (GwsegEPSBase)
{
// This is the eps (but no hysteris) version of the Gwseg model.
// However, only default scaling parameters, i.e no scaling.
//std::cout << "==================================== GwsegEPSBase ====================================" << std::endl;
Opm::parameter::ParameterGroup param;
Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("satfuncEPSBase.DATA");
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
const int np = props.numPhases();
const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
BOOST_REQUIRE(np == 3);
BOOST_REQUIRE(wpos == 0);
BOOST_REQUIRE(opos == 1);
BOOST_REQUIRE(gpos == 2);
const int n=11;
double s[n*np];
int cells[n];
double kr[n*np];
double dkrds[n*np*np];
for (int i=0; i<n; ++i) {
cells[i] = 0;
s[i*np+wpos] = i*0.1;
s[i*np+opos] = 1.0-s[i*np+wpos];
s[i*np+gpos] = 0.0;
}
props.relperm(n, s, cells, kr, dkrds);
double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0,-0.66666666666666741, -0.75, -0.8, -0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
const double reltol = 1.0e-6;
for (int i=0; i<n; ++i) {
BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[i], reltol);
BOOST_CHECK_CLOSE(kr[i*np+opos], kro[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
}
/*
std::cout << std::setw(12) << "sw";
std::cout << std::setw(12) << "so";
std::cout << std::setw(12) << "sg";
std::cout << std::setw(12) << "krw";
std::cout << std::setw(12) << "kro";
std::cout << std::setw(12) << "krg";
std::cout << std::setw(12) << "DkrwDsw";
std::cout << std::setw(12) << "DkroDsw";
std::cout << std::setw(12) << "DkrgDsw";
std::cout << std::setw(12) << "DkrwDso";
std::cout << std::setw(12) << "DkroDso";
std::cout << std::setw(12) << "DkrgDso";
std::cout << std::setw(12) << "DkrwDsg";
std::cout << std::setw(12) << "DkroDsg";
std::cout << std::setw(12) << "DkrgDsg";
std::cout << std::endl;
for (int i=0; i<n; ++i) {
std::cout << std::setw(12) << s[i*np+wpos] << std::setw(12) << s[i*np+opos] << std::setw(12) << s[i*np+gpos]
<< std::setw(12) << kr[i*np+wpos] << std::setw(12) << kr[i*np+opos] << std::setw(12) << kr[i*np+gpos];
for (int j=0; j<np*np; ++j) {
std::cout << std::setw(12) << dkrds[i*np*np+j];
}
std::cout << std::endl;
}
*/
}
BOOST_AUTO_TEST_CASE (GwsegEPS_A)
{
// This is the eps (but no hysteris) version of the Gwseg model.
// Scaling parameters from keyword SWL etc.
//std::cout << "==================================== GwsegEPS_A ====================================" << std::endl;
Opm::parameter::ParameterGroup param;
Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_A.DATA");
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
const int np = props.numPhases();
const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
BOOST_REQUIRE(np == 3);
BOOST_REQUIRE(wpos == 0);
BOOST_REQUIRE(opos == 1);
BOOST_REQUIRE(gpos == 2);
const int n=11;
double s[n*np];
int cells[n];
double kr[n*np];
double dkrds[n*np*np];
const int ncell = 8;
double krw[ncell][n] = {{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
{0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7},
{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
{0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7}};
double kro[ncell][n] = {{1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
double DkroDsg[ncell][n] = {{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}};
for (int icell=0; icell<ncell; ++icell) {
for (int i=0; i<n; ++i) {
cells[i] = icell;
s[i*np+wpos] = i*0.1;
s[i*np+opos] = 1.0-s[i*np+wpos];
s[i*np+gpos] = 0.0;
}
props.relperm(n, s, cells, kr, dkrds);
const double reltol = 1.0e-3;
for (int i=0; i<n; ++i) {
BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[icell][i], reltol);
BOOST_CHECK_CLOSE(kr[i*np+opos], kro[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[icell][i], reltol);
}
/*
std::cout << std::setw(12) << "sw: ";
for (int i=0; i<n; ++i) std::cout << s[i*np+wpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "so: ";
for (int i=0; i<n; ++i) std::cout << s[i*np+opos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "sg: ";
for (int i=0; i<n; ++i) std::cout << s[i*np+gpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "krw: ";
for (int i=0; i<n; ++i) std::cout << kr[i*np+wpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "kro: ";
for (int i=0; i<n; ++i) std::cout << kr[i*np+opos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "krg: ";
for (int i=0; i<n; ++i) std::cout << kr[i*np+gpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkrwDsw: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+wpos*np+wpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkroDsw: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+wpos*np+opos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkrgDsw: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+wpos*np+gpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkrwDso: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+opos*np+wpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkroDso: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+opos*np+opos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkrgDso: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+opos*np+gpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkrwDsg: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+gpos*np+wpos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkroDsg: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+gpos*np+opos] << ", ";
std::cout << std::endl;
std::cout << std::setw(12) << "DkrgDsg: ";
for (int i=0; i<n; ++i) std::cout << dkrds[i*np*np+gpos*np+gpos] << ", ";
std::cout << std::endl;
std::cout << std::endl;
*/
}
}
BOOST_AUTO_TEST_CASE (GwsegEPS_B)
{
// This is the eps (but no hysteris) version of the Gwseg model.
// Scaling parameters from ENPTVD table.
//std::cout << "==================================== GwsegEPS_B ====================================" << std::endl;
/*
Opm::parameter::ParameterGroup param;
Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_B.DATA");
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
const int np = props.numPhases();
const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
BOOST_REQUIRE(np == 3);
BOOST_REQUIRE(wpos == 0);
BOOST_REQUIRE(opos == 1);
BOOST_REQUIRE(gpos == 2);
const int n=11;
double s[n*np];
int cells[n];
double kr[n*np];
double dkrds[n*np*np];
const int ncell = 8;
double krw[ncell][n] = {{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
{0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7},
{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
{0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7}};
double kro[ncell][n] = {{1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
double DkroDsg[ncell][n] = {{-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-2.32831e-10, -2.32831e-10, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-2.32831e-10, -2.32831e-10, -2.32831e-10, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}};
for (int icell=0; icell<ncell; ++icell) {
for (int i=0; i<n; ++i) {
cells[i] = icell;
s[i*np+wpos] = i*0.1;
s[i*np+opos] = 1.0-s[i*np+wpos];
s[i*np+gpos] = 0.0;
}
props.relperm(n, s, cells, kr, dkrds);
const double reltol = 1.0e-3;
for (int i=0; i<n; ++i) {
BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[icell][i], reltol);
BOOST_CHECK_CLOSE(kr[i*np+opos], kro[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[icell][i], reltol);
}
}
*/
}
BOOST_AUTO_TEST_CASE (GwsegEPS_C)
{
// This is the eps (but no hysteris) version of the Gwseg model.
// Scaling parameters given the "norne-way", i.e EQUALS, COPY, ADD, MULTIPLY.
//std::cout << "==================================== GwsegEPS_C ====================================" << std::endl;
Opm::parameter::ParameterGroup param;
Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_C.DATA");
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
const int np = props.numPhases();
const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
BOOST_REQUIRE(np == 3);
BOOST_REQUIRE(wpos == 0);
BOOST_REQUIRE(opos == 1);
BOOST_REQUIRE(gpos == 2);
const int n=11;
double s[n*np];
int cells[n];
double kr[n*np];
double dkrds[n*np*np];
const int ncell = 8;
double krw[ncell][n] = {{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
{0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7},
{0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7},
{0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.14, 0.28, 0.42, 0.56, 0.7, 0.7},
{0, 0, 0, 0, 0, 0.233333, 0.466667, 0.7, 0.7, 0.7, 0.7}};
double kro[ncell][n] = {{1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0},
{1, 1, 1, 0.766667, 0.533333, 0.35, 0.233333, 0.116667, 0, 0, 0}};
double DkrwDsw[ncell][n] = {{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0},
{0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0},
{0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1.4, 1.4, 1.4, 1.4, 0, 0},
{0, 0, 0, 0, 0, 2.33333, 2.33333, 0, 0, 0, 0}};
double DkroDsw[ncell][n] = {{0.0, 0.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, -2, -2, -1, -1, -1, -1, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0},
{0, 0, 0, -2.33333, -2.33333, -1.16667, -1.16667, -1.16667, 0, 0, 0}};
double DkroDsg[ncell][n] = {{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3, -2, -0.666667, -0.75, -0.8, -0.833333, 0.142857, 0, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0},
{-3.0, -3.0, -3.0, -3.14286, -2.14286, -0.809524, -0.892857, -0.942857, 0.190476, 3.17207e-16, 0}};
for (int icell=0; icell<ncell; ++icell) {
for (int i=0; i<n; ++i) {
cells[i] = icell;
s[i*np+wpos] = i*0.1;
s[i*np+opos] = 1.0-s[i*np+wpos];
s[i*np+gpos] = 0.0;
}
props.relperm(n, s, cells, kr, dkrds);
const double reltol = 1.0e-3;
for (int i=0; i<n; ++i) {
BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[icell][i], reltol);
BOOST_CHECK_CLOSE(kr[i*np+opos], kro[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[icell][i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[icell][i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE (GwsegEPS_D)
{
// This is the eps and hysteris version of the Gwseg model.
// However, only default scaling parameters, i.e no scaling.
//std::cout << "==================================== GwsegEPS_D ====================================" << std::endl;
Opm::parameter::ParameterGroup param;
Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_D.DATA");
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false);
const int np = props.numPhases();
const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua];
const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid];
const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour];
BOOST_REQUIRE(np == 3);
BOOST_REQUIRE(wpos == 0);
BOOST_REQUIRE(opos == 1);
BOOST_REQUIRE(gpos == 2);
const int n=11;
double s[n*np];
int cells[n];
double kr[n*np];
double dkrds[n*np*np];
for (int i=0; i<n; ++i) {
cells[i] = 0;
s[i*np+wpos] = i*0.1;
s[i*np+opos] = 1.0-s[i*np+wpos];
s[i*np+gpos] = 0.0;
}
props.relperm(n, s, cells, kr, dkrds);
double krw[11] = {0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7};
double kro[11] = {1.0, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1, 0.0, 0.0, 0.0};
double DkrwDsw[11] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
double DkroDsw[11] = {-2.0, -2.0, -2.0, -2.0, -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 0.0};
double DkroDsg[11] = {-5.0, -5.0, -3.0, -2.0, -0.66666666666666741, -0.75, -0.8,
-0.83333333333333237, 0.14285714285714296, 0.0, 0.0};
const double reltol = 1.0e-6;
for (int i=0; i<n; ++i) {
BOOST_CHECK_CLOSE(kr[i*np+wpos], krw[i], reltol);
BOOST_CHECK_CLOSE(kr[i*np+opos], kro[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+wpos], DkrwDsw[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+opos], DkroDsw[i], reltol);
BOOST_CHECK_CLOSE(dkrds[i*np*np+np*gpos+opos], DkroDsg[i], reltol);
}
/*
std::cout << std::setw(12) << "sw";
std::cout << std::setw(12) << "so";
std::cout << std::setw(12) << "sg";
std::cout << std::setw(12) << "krw";
std::cout << std::setw(12) << "kro";
std::cout << std::setw(12) << "krg";
std::cout << std::setw(12) << "DkrwDsw";
std::cout << std::setw(12) << "DkroDsw";
std::cout << std::setw(12) << "DkrgDsw";
std::cout << std::setw(12) << "DkrwDso";
std::cout << std::setw(12) << "DkroDso";
std::cout << std::setw(12) << "DkrgDso";
std::cout << std::setw(12) << "DkrwDsg";
std::cout << std::setw(12) << "DkroDsg";
std::cout << std::setw(12) << "DkrgDsg";
std::cout << std::endl;
for (int i=0; i<n; ++i) {
std::cout << std::setw(12) << s[i*np+wpos] << std::setw(12) << s[i*np+opos] << std::setw(12) << s[i*np+gpos]
<< std::setw(12) << kr[i*np+wpos] << std::setw(12) << kr[i*np+opos] << std::setw(12) << kr[i*np+gpos];
for (int j=0; j<np*np; ++j) {
std::cout << std::setw(12) << dkrds[i*np*np+j];
}
std::cout << std::endl;
}
*/
}
BOOST_AUTO_TEST_SUITE_END()