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

Merge pull request #583 from andlaus/remove_EclipseGridParser

Browse Source 
Remove eclipse grid parser
This commit is contained in:
Bård Skaflestad committed 2014-06-04 00:26:41 +02:00
2 parents 55427b2220 c1c3f9d2e0
commit 8afa34db4d
65 files changed
+491 -8022
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

CMakeLists_files.cmake
-9
View File
@@ -39,7 +39,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/grid/cpgpreprocess/preprocess.c
opm/core/grid/cpgpreprocess/uniquepoints.c
opm/core/io/eclipse/EclipseGridInspector.cpp
opm/core/io/eclipse/EclipseGridParser.cpp
opm/core/io/eclipse/EclipseWriter.cpp
opm/core/io/eclipse/writeECLData.cpp
opm/core/io/OutputWriter.cpp
@@ -223,15 +222,10 @@ list (APPEND EXAMPLE_SOURCE_FILES
# originally generated with the command:
# find attic -name '*.c*' -printf '\t%p\n' | sort
list (APPEND ATTIC_FILES
attic/bo_resprop_test.cpp
attic/pvt_test.cpp
attic/relperm_test.cpp
attic/test_cfs_tpfa.c
attic/test_ert.cpp
attic/test_jacsys.cpp
attic/test_lapack.cpp
attic/test_read_grid.c
attic/test_readpolymer.cpp
attic/test_read_vag.cpp
attic/test_writeVtkData.cpp
)
@@ -261,11 +255,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/grid/cpgpreprocess/uniquepoints.h
opm/core/io/eclipse/CornerpointChopper.hpp
opm/core/io/eclipse/EclipseGridInspector.hpp
opm/core/io/eclipse/EclipseGridParser.hpp
opm/core/io/eclipse/EclipseGridParserHelpers.hpp
opm/core/io/eclipse/EclipseUnits.hpp
opm/core/io/eclipse/EclipseWriter.hpp
opm/core/io/eclipse/SpecialEclipseFields.hpp
opm/core/io/eclipse/writeECLData.hpp
opm/core/io/OutputWriter.hpp
opm/core/io/vag/vag.hpp
attic/bo_resprop_test.cpp
-55
View File
@@ -1,55 +0,0 @@
/*
Copyright 2012 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"
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/io/eclipse/EclipseGridInspector.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/grid.h>
#include <algorithm>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>
int main(int argc, char** argv)
{
// Parameters.
Opm::parameter::ParameterGroup param(argc, argv);
// Parser.
std::string ecl_file = param.get<std::string>("filename");
Opm::EclipseGridParser deck(ecl_file);
UnstructuredGrid grid;
grid.number_of_cells = 1;
grid.global_cell = NULL;
Opm::BlackoilPropertiesFromDeck props(deck, grid, param);
const int n = 1;
double p[n] = { 150e5 };
double z[3*n] = { 1.0, 1.0, 1.0 };
int cells[n] = { 0 };
double A[9*n];
props.matrix(n, p, z, cells, A, 0);
std::copy(A, A + 9, std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
}
attic/pvt_test.cpp
-145
View File
@@ -1,145 +0,0 @@
/*
Copyright 2012 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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/io/eclipse/EclipseGridInspector.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/grid.h>
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <iterator>
#include <string>
#include <vector>
int main(int argc, char** argv)
{
using namespace std;
// Parameters.
Opm::parameter::ParameterGroup param(argc, argv);
// Parser.
std::string ecl_file = param.get<std::string>("deck_filename");
std::string input_file = param.get<std::string>("input_filename");
std::string matrix_output = param.get<std::string>("matrix_output");
std::string rs_output = param.get<std::string>("rs_output");
std::string b_output = param.get<std::string>("b_output");
std::string mu_output = param.get<std::string>("mu_output");
Opm::EclipseGridParser deck(ecl_file);
UnstructuredGrid grid;
grid.number_of_cells = 1;
grid.global_cell = NULL;
grid.dimensions = 3;
grid.cartdims[0] = 1;
grid.cartdims[1] = 1;
grid.cartdims[2] = 1;
Opm::BlackoilPropertiesFromDeck props(deck, grid, param);
Opm::BlackoilPvtProperties pvt;
int samples = param.getDefault("dead_tab_size", 1025);
pvt.init(deck, samples);
const int n = 1;
std::fstream inos(input_file.c_str());
if(!inos.good()){
std::cout << "Could not open :" << input_file << std::endl;
exit(3);
}
std::fstream aos(matrix_output.c_str(), std::fstream::out | std::fstream::trunc);
aos << setiosflags(ios::scientific) << setprecision(12);
if(!(aos.good())){
std::cout << "Could not open"<< matrix_output << std::endl;
exit(3);
}
std::fstream bos(b_output.c_str(), std::fstream::out | std::fstream::trunc);
bos << setiosflags(ios::scientific) << setprecision(12);
if(!(bos.good())){
std::cout << "Could not open"<< b_output << std::endl;
exit(3);
}
std::fstream rsos(rs_output.c_str(), std::fstream::out | std::fstream::trunc);
rsos << setiosflags(ios::scientific) << setprecision(12);
if(!(rsos.good())){
std::cout << "Could not open"<< rs_output << std::endl;
exit(3);
}
std::fstream muos(mu_output.c_str(), std::fstream::out | std::fstream::trunc);
if(!(muos.good())){
std::cout << "Could not open"<< rs_output << std::endl;
exit(3);
}
const int np = props.numPhases();
const int max_np = 3;
if (np > max_np) {
OPM_THROW(std::runtime_error, "Max #phases is 3.");
}
while((inos.good()) && (!inos.eof())){
double p[n];
double z[max_np*n];
int cells[n] = { 0 };
inos >> p[0];
for(int i=0; i < np; ++i){
inos >> z[i];
}
if(inos.good()){
double A[max_np*max_np*n];
double dA[max_np*max_np*n];
props.matrix(n, p, z, cells, A, dA);
std::copy(A, A + np*np*n, std::ostream_iterator<double>(aos, " "));
std::copy(dA, dA + np*np*n, std::ostream_iterator<double>(aos, " "));
aos << std::endl;
double mu[max_np];
//double dmu[max_np];//not implemented
props.viscosity(n, p, z, cells, mu, 0);
std::copy(mu, mu + np*n, std::ostream_iterator<double>(muos, " "));
//std::copy(dmu, dmu + np*n, std::ostream_iterator<double>(muos, " "));
aos << std::endl;
double b[max_np];
double dbdp[max_np];
pvt.dBdp(n, p, z, b, dbdp);
std::copy(b, b + np*n, std::ostream_iterator<double>(bos, " "));
std::copy(dbdp, dbdp + np*n, std::ostream_iterator<double>(bos, " "));
bos << std::endl;
double rs[max_np];
double drs[max_np];
//pvt.R(n, p, z, rs);
pvt.dRdp(n, p, z, rs,drs);
std::copy(rs, rs + np*n, std::ostream_iterator<double>(rsos, " "));
std::copy(drs, drs + np*n, std::ostream_iterator<double>(rsos, " "));
rsos << std::endl;
}
}
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
attic/relperm_test.cpp
-97
View File
@@ -1,97 +0,0 @@
/*
Copyright 2012 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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/io/eclipse/EclipseGridInspector.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/grid.h>
#include <algorithm>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>
#include <fstream>
int main(int argc, char** argv)
{
// Parameters.
Opm::parameter::ParameterGroup param(argc, argv);
// Parser.
std::string ecl_file = param.get<std::string>("deck_filename");
std::string input_file = param.get<std::string>("input_filename");
std::string relperm_output = param.get<std::string>("relperm_output");
//std::string relperm_output = param.get<std::string>("relperm_outout");
Opm::EclipseGridParser deck(ecl_file);
UnstructuredGrid grid;
grid.number_of_cells = 1;
grid.global_cell = NULL;
grid.dimensions = 3;
grid.cartdims[0] = 1;
grid.cartdims[1] = 1;
grid.cartdims[2] = 1;
Opm::BlackoilPropertiesFromDeck props(deck, grid, param);
std::fstream inos(input_file.c_str());//, std::fstream::in);
if(!inos.good()){
std::cout << "Could not open :" << input_file << std::endl;
exit(3);
}
std::fstream kros(relperm_output.c_str(), std::fstream::out | std::fstream::trunc);
if(!kros.good()){
std::cout << "Could not open :" << input_file << std::endl;
exit(3);
}
const int np = props.numPhases();
const int max_np = 3;
if (np > max_np) {
OPM_THROW(std::runtime_error, "Max #phases is 3.");
}
while((inos.good()) && (!inos.eof())){
double s[max_np];
for(int i=0; i < np; ++i){
inos >> s[i];
}
if(inos.good()){
double kr[max_np];
double dkr[max_np*max_np];
int cell[1];
cell[0]=1;
props.relperm(1,s, cell, kr, dkr);
std::copy(s, s + np, std::ostream_iterator<double>(kros, " "));
kros << " ";
std::copy(kr, kr + np, std::ostream_iterator<double>(kros, " "));
kros << " ";
std::copy(dkr, dkr + np*np, std::ostream_iterator<double>(kros, " "));
kros << "\n";
}
}
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
attic/test_ert.cpp
-262
View File
@@ -1,262 +0,0 @@
/*
Copyright 2012 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"
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid.h>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <ert/ecl/ecl_grid.h>
#include <cstdio>
#include <memory>
#include <iostream>
using namespace std;
#if 0
static
void cell_nodes(const UnstructuredGrid * c_grid , int cell , std::vector<int>& nodes) {
int face_offset = c_grid->cell_facepos[cell];
int num_faces = c_grid->cell_facepos[cell + 1] - face_offset;
nodes.clear();
//printf("cell: %d \n",cell);
for (int iface = 0; iface < num_faces; iface++) {
int face = c_grid->cell_faces[ face_offset + iface];
//printf("face[%d] = %d \n",iface , face );
{
int node_offset = c_grid->face_nodepos[ face ];
int num_nodes = c_grid->face_nodepos[ face + 1] - node_offset;
for (int inode = 0; inode < num_nodes; inode++) {
int node = c_grid->face_nodes[ inode + node_offset ];
//printf(" node[%d] = %d \n",inode , node);
nodes.push_back( node );
}
}
}
{
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
sort( nodes.begin() , nodes.end());
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
unique( nodes.begin() , nodes.end() );
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
nodes.resize( 8 );
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
}
}
#endif
/*
static
void eclExport(Opm::GridManager& grid) {
const UnstructuredGrid * c_grid = grid.c_grid();
printf("dimensions : %d \n",c_grid->dimensions);
printf("number of cells : %d \n",c_grid->number_of_cells);
printf("number of nodes : %d \n",c_grid->number_of_nodes);
printf("number of faces : %d \n",c_grid->number_of_faces);
printf("length(face_nodes) : %d \n",c_grid->face_nodepos[ c_grid->number_of_faces ]);
printf("cartdims : %d %d %d \n",
c_grid->cartdims[0] ,
c_grid->cartdims[1] ,
c_grid->cartdims[2]);
printf("global_cell : %d %d %d %d %d\n",
c_grid->global_cell[0] ,
c_grid->global_cell[1] ,
c_grid->global_cell[2] ,
c_grid->global_cell[3] ,
c_grid->global_cell[4]);
{
std::vector<int> nodes;
cell_nodes( c_grid , 10 , nodes );
cell_nodes( c_grid , 15 , nodes );
cell_nodes( c_grid , 20 , nodes );
cell_nodes( c_grid , 25 , nodes );
}
{
ecl_grid_type * ecl_grid;
int num_coords = c_grid->number_of_cells;
int coords_size = 4;
int nx = c_grid->cartdims[0];
int ny = c_grid->cartdims[1];
int nz = c_grid->cartdims[2];
int ** coords;
float ** corners;
// float * mapaxes = NULL;
std::vector<int> nodes;
corners = (float **) malloc( num_coords * sizeof * corners );
coords = (int **) malloc( num_coords * sizeof * coords );
{
int c;
for (c=0; c < num_coords; c++) {
corners[c] = (float *) malloc( 24 * sizeof * corners[c] );
coords[c] = (int *) malloc( coords_size * sizeof * coords[c] );
}
for (c=0; c < num_coords; c++) {
cell_nodes( c_grid , c , nodes );
for (int p=0; p < 8; p++) {
int n = nodes[p];
for (int d=0; d < c_grid->dimensions; d++)
corners[c][3*p + d] = c_grid->node_coordinates[ c_grid->dimensions * n + d ];
}
{
int i,j,k;
{
int g = c_grid->global_cell[ c ];
k = g / nx*ny; g -= k * nx*ny;
j = g / nx; g -= j * nx;
i = g;
}
coords[c][0] = i + 1;
coords[c][1] = j + 1;
coords[c][2] = k + 1;
coords[c][3] = c_grid->global_cell[ c ] + 1;
}
}
}
ecl_grid = ecl_grid_alloc( "/private/joaho/ERT/NR/libenkf/src/Gurbat/EXAMPLE_01_BASE.EGRID" );
printf("Grid loaded ... \n");
ecl_grid_free( ecl_grid );
printf("Grid discarded ... \n");
ecl_grid = ecl_grid_alloc_GRID_data( num_coords , nx , ny , nz , coords_size , coords , corners , NULL );
ecl_grid_fwrite_GRID( ecl_grid , "/tmp/test.GRID" );
{
FILE * stream = fopen( "/tmp/test.grdecl" , "w");
ecl_grid_fprintf_grdecl( ecl_grid , stream );
fclose( stream );
}
ecl_grid_free( ecl_grid );
{
for (int c=0; c < num_coords; c++) {
free(corners[c]);
free(coords[c]);
}
}
free( corners );
free( coords );
}
}
*/
/*
#ifdef HAVE_ERT
ecl_grid_type * create_ecl_grid( const struct UnstructuredGrid * g) {
int num_coords = g->number_of_cells;
int nx = g->cartdims[0];
int ny = g->cartdims[1];
int nz = g->cartdims[2];
int coords_size = 4;
int ** coords;
float ** corners;
float * mapaxes = NULL;
corners = malloc( num_coords * sizeof * corners );
coords = malloc( num_coords * sizeof * coords );
{
for (int c=0; c < num_coords; c++) {
corners[c] = malloc( 24 * sizeof * corners[0] );
coords[c] = malloc( coords_size * sizeof * coords[0] );
}
}
{
for (int k=0; k < nz; k++) {
for (int j=0; j < ny; j++) {
for (int i=0; i < nx; i++) {
int global_index = i + j*nx + k*nx*ny;
coords[global_index][0] = i;
coords[global_index][1] = j;
coords[global_index][2] = k;
coords[global_index][3] = 1;
}
}
}
}
{
for (int c=0; c < num_coords; c++) {
free(corners[c]);
free(coords[c]);
}
}
free( corners );
free( coords );
}
#endif
*/
// struct grdecl : opm/core/grid/cpgpreprocess/preprocess.h
// struct processes_grid : opm/core/grid/cpgpreprocess/preprocess.h
int main(int /*argc*/ , char **argv)
{
std::string filename( argv[1] );
const std::unique_ptr<Opm::GridManager> grid;
const std::unique_ptr<Opm::IncompPropertiesInterface> props;
Opm::EclipseGridParser eclParser(filename , false);
//eclParser.saveEGRID_INIT("/tmp" , "OPM" );
grid.reset(new Opm::GridManager(eclParser));
props.reset(new Opm::IncompPropertiesFromDeck(eclParser , *grid->c_grid()));
}
attic/test_readpolymer.cpp
-62
View File
@@ -1,62 +0,0 @@
//===========================================================================
//
// File: test_readpolymer.cpp
//
// Created: Thu Jan 12 15:18:46 2012
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <iostream>
// Test program for reading Eclipse Polymer keywords.
int main(int argc, char** argv)
{
using namespace std;
std::string ecl_filename;
if (argc == 2) {
ecl_filename = argv[1];
} else {
std::cout << "\nUsage: test_readpolymer filename.grdecl\n";
exit( 1 );
}
bool convert_to_SI = true;
Opm::EclipseGridParser parser(ecl_filename, convert_to_SI);
std::cout << "\n Polymer fields\n\n";
if (parser.hasField("PLYVISC")) {
parser.getPLYVISC().write(std::cout);
}
if (parser.hasField("PLYROCK")) {
parser.getPLYROCK().write(std::cout);
}
if (parser.hasField("PLYADS")) {
parser.getPLYADS().write(std::cout);
}
if (parser.hasField("TLMIXPAR")) {
parser.getTLMIXPAR().write(std::cout);
}
if (parser.hasField("PLYMAX")) {
parser.getPLYMAX().write(std::cout);
}
if (parser.hasField("WPOLYMER")) {
parser.getWPOLYMER().write(std::cout);
}
}
examples/compute_tof.cpp
+7 -7
View File
@@ -104,7 +104,7 @@ try
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
Opm::DeckConstPtr newParserDeck;
Opm::DeckConstPtr deck;
std::unique_ptr<GridManager> grid;
std::unique_ptr<IncompPropertiesInterface> props;
std::unique_ptr<Opm::WellsManager> wells;
@@ -115,22 +115,22 @@ try
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr parser(new Opm::Parser());
newParserDeck = parser->parseFile(deck_filename);
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParserDeck));
deck = parser->parseFile(deck_filename);
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
// Grid init
grid.reset(new GridManager(newParserDeck));
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(newParserDeck, *grid->c_grid()));
props.reset(new IncompPropertiesFromDeck(deck, *grid->c_grid()));
// Wells init.
wells.reset(new Opm::WellsManager(eclipseState , 0 , *grid->c_grid(), props->permeability()));
// Gravity.
gravity[2] = newParserDeck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, newParserDeck, gravity[2], state);
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
} else {
// Grid init.
examples/sim_2p_comp_reorder.cpp
+9 -9
View File
@@ -90,7 +90,7 @@ try
std::unique_ptr<RockCompressibility> rock_comp;
ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck;
Opm::DeckConstPtr deck;
BlackoilState state;
// bool check_well_controls = false;
@@ -98,24 +98,24 @@ try
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
newParserDeck = parser->parseFile(deck_filename);
deck = parser->parseFile(deck_filename);
eclipseState.reset(new EclipseState(newParserDeck));
eclipseState.reset(new EclipseState(deck));
// Grid init
grid.reset(new GridManager(newParserDeck));
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(newParserDeck, *grid->c_grid(), param));
props.reset(new BlackoilPropertiesFromDeck(deck, *grid->c_grid(), param));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(newParserDeck));
rock_comp.reset(new RockCompressibility(deck));
// Gravity.
gravity[2] = newParserDeck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, newParserDeck, gravity[2], state);
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
initBlackoilSurfvol(*grid->c_grid(), *props, state);
} else {
@@ -225,7 +225,7 @@ try
int step = 0;
SimulatorTimer simtimer;
// Use timer for last epoch to obtain total time.
Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck));
Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
simtimer.init(timeMap);
const double total_time = simtimer.totalTime();
for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
examples/sim_2p_incomp.cpp
+9 -9
View File
@@ -100,7 +100,7 @@ try
bool use_deck = param.has("deck_filename");
EclipseStateConstPtr eclipseState;
Opm::DeckConstPtr newParserDeck;
Opm::DeckConstPtr deck;
std::unique_ptr<GridManager> grid;
std::unique_ptr<IncompPropertiesInterface> props;
std::unique_ptr<RockCompressibility> rock_comp;
@@ -112,23 +112,23 @@ try
ParserPtr parser(new Opm::Parser());
std::string deck_filename = param.get<std::string>("deck_filename");
newParserDeck = parser->parseFile(deck_filename);
eclipseState.reset( new EclipseState(newParserDeck));
deck = parser->parseFile(deck_filename);
eclipseState.reset( new EclipseState(deck));
// Grid init
grid.reset(new GridManager(newParserDeck));
grid.reset(new GridManager(deck));
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(newParserDeck, *grid->c_grid()));
props.reset(new IncompPropertiesFromDeck(deck, *grid->c_grid()));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
rock_comp.reset(new RockCompressibility(newParserDeck));
rock_comp.reset(new RockCompressibility(deck));
// Gravity.
gravity[2] = newParserDeck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, newParserDeck, gravity[2], state);
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
} else {
// Grid init.
@@ -214,7 +214,7 @@ try
param.writeParam(output_dir + "/simulation.param");
}
Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck));
Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< " (number of report steps: "
examples/wells_example.cpp
+6 -6
View File
@@ -35,16 +35,16 @@ try
// Read input file
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck = parser->parseFile(file_name);
Opm::DeckConstPtr deck = parser->parseFile(file_name);
std::cout << "Done!" << std::endl;
// Setup grid
GridManager grid(newParserDeck);
GridManager grid(deck);
// Define rock and fluid properties
IncompPropertiesFromDeck incomp_properties(newParserDeck, *grid.c_grid());
RockCompressibility rock_comp(newParserDeck);
EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParserDeck));
IncompPropertiesFromDeck incomp_properties(deck, *grid.c_grid());
RockCompressibility rock_comp(deck);
EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
// Finally handle the wells
WellsManager wells(eclipseState , 0 , *grid.c_grid(), incomp_properties.permeability());
@@ -76,7 +76,7 @@ try
Opm::TwophaseState state;
initStateFromDeck(*grid.c_grid(), incomp_properties, newParserDeck, gravity[2], state);
initStateFromDeck(*grid.c_grid(), incomp_properties, deck, gravity[2], state);
Opm::WellState well_state;
well_state.init(wells.c_wells(), state);
opm/core/doxygen_main.hpp
+2 -2
View File
@@ -48,7 +48,7 @@ creation and destruction of an UnstructuredGrid. The method
Opm::GridManager::c_grid() provides access to the underlying
UnstructuredGrid structure. This class also provides an easy way
to initialize a grid from an ECLIPSE-format input deck, via the
constructor taking an Opm::EclipseGridParser.
constructor taking an Opm::DeckConstPtr.
<h3>Well handling</h3>
@@ -109,7 +109,7 @@ There are two simulator programs for two-phase immiscible flow in opm-core:
<h3>Various utilities</h3>
Utilities contained in opm-core include:
- IO utilities (Opm::EclipseGridParser, binary I/O via the ERT library, vtk output)
- IO utilities (binary output via the ERT library, vtk output)
- Interpolation utilities (Opm::MonotCubicInterpolator, Opm::VelocityInterpolationECVI)
- Support for SI and non-SI units (Opm::unit and Opm::prefix)
- Low-order quadratures for general geometries (Opm::CellQuadrature, Opm::FaceQuadrature)
opm/core/grid/GridManager.cpp
+31 -141
View File
@@ -20,10 +20,10 @@
#include "config.h"
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <array>
#include <algorithm>
@@ -32,30 +32,6 @@
namespace Opm
{
/// Construct a 3d corner-point grid from a deck.
GridManager::GridManager(const Opm::EclipseGridParser& deck)
{
// We accept two different ways to specify the grid.
// 1. Corner point format.
// Requires ZCORN, COORDS, DIMENS or SPECGRID, optionally
// ACTNUM, optionally MAPAXES.
// For this format, we will verify that DXV, DYV, DZV,
// DEPTHZ and TOPS are not present.
// 2. Tensor grid format.
// Requires DXV, DYV, DZV, optionally DEPTHZ or TOPS.
// For this format, we will verify that ZCORN, COORDS
// and ACTNUM are not present.
// Note that for TOPS, we only allow a uniform vector of values.
if (deck.hasField("ZCORN") && deck.hasField("COORD")) {
initFromDeckCornerpoint(deck);
} else if (deck.hasField("DXV") && deck.hasField("DYV") && deck.hasField("DZV")) {
initFromDeckTensorgrid(deck);
} else {
OPM_THROW(std::runtime_error, "Could not initialize grid from deck. "
"Need either ZCORN + COORD or DXV + DYV + DZV keywords.");
}
}
GridManager::GridManager(Opm::EclipseGridConstPtr eclipseGrid) {
struct grdecl g;
std::vector<int> actnum;
@@ -84,8 +60,7 @@ namespace Opm
}
/// Construct a 3d corner-point grid from a deck.
GridManager::GridManager(Opm::DeckConstPtr newParserDeck)
GridManager::GridManager(Opm::DeckConstPtr deck)
{
// We accept two different ways to specify the grid.
// 1. Corner point format.
@@ -99,10 +74,10 @@ namespace Opm
// and ACTNUM are not present.
// Note that for TOPS, we only allow a uniform vector of values.
if (newParserDeck->hasKeyword("ZCORN") && newParserDeck->hasKeyword("COORD")) {
initFromDeckCornerpoint(newParserDeck);
} else if (newParserDeck->hasKeyword("DXV") && newParserDeck->hasKeyword("DYV") && newParserDeck->hasKeyword("DZV")) {
initFromDeckTensorgrid(newParserDeck);
if (deck->hasKeyword("ZCORN") && deck->hasKeyword("COORD")) {
initFromDeckCornerpoint(deck);
} else if (deck->hasKeyword("DXV") && deck->hasKeyword("DYV") && deck->hasKeyword("DZV")) {
initFromDeckTensorgrid(deck);
} else {
OPM_THROW(std::runtime_error, "Could not initialize grid from deck. "
"Need either ZCORN + COORD or DXV + DYV + DZV keywords.");
@@ -165,11 +140,6 @@ namespace Opm
}
}
void GridManager::saveEGRID(const std::string& filename , const Opm::EclipseGridParser& deck) {
deck.saveEGRID( filename , ug_->number_of_cells , ug_->global_cell );
}
/// Destructor.
GridManager::~GridManager()
{
@@ -187,28 +157,13 @@ namespace Opm
return ug_;
}
// Construct corner-point grid from deck.
void GridManager::initFromDeckCornerpoint(const Opm::EclipseGridParser& deck)
{
// Extract data from deck.
// Collect in input struct for preprocessing.
struct grdecl grdecl = deck.get_grdecl();
// Process grid.
ug_ = create_grid_cornerpoint(&grdecl, 0.0);
if (!ug_) {
OPM_THROW(std::runtime_error, "Failed to construct grid.");
}
}
// Construct corner-point grid from deck.
void GridManager::initFromDeckCornerpoint(Opm::DeckConstPtr newParserDeck)
void GridManager::initFromDeckCornerpoint(Opm::DeckConstPtr deck)
{
// Extract data from deck.
// Collect in input struct for preprocessing.
struct grdecl grdecl;
createGrdecl(newParserDeck, grdecl);
createGrdecl(deck, grdecl);
// Process grid.
ug_ = create_grid_cornerpoint(&grdecl, 0.0);
@@ -217,24 +172,24 @@ namespace Opm
}
}
void GridManager::createGrdecl(Opm::DeckConstPtr newParserDeck, struct grdecl &grdecl)
void GridManager::createGrdecl(Opm::DeckConstPtr deck, struct grdecl &grdecl)
{
// Extract data from deck.
const std::vector<double>& zcorn = newParserDeck->getKeyword("ZCORN")->getSIDoubleData();
const std::vector<double>& coord = newParserDeck->getKeyword("COORD")->getSIDoubleData();
const std::vector<double>& zcorn = deck->getKeyword("ZCORN")->getSIDoubleData();
const std::vector<double>& coord = deck->getKeyword("COORD")->getSIDoubleData();
const int* actnum = NULL;
if (newParserDeck->hasKeyword("ACTNUM")) {
actnum = &(newParserDeck->getKeyword("ACTNUM")->getIntData()[0]);
if (deck->hasKeyword("ACTNUM")) {
actnum = &(deck->getKeyword("ACTNUM")->getIntData()[0]);
}
std::array<int, 3> dims;
if (newParserDeck->hasKeyword("DIMENS")) {
Opm::DeckKeywordConstPtr dimensKeyword = newParserDeck->getKeyword("DIMENS");
if (deck->hasKeyword("DIMENS")) {
Opm::DeckKeywordConstPtr dimensKeyword = deck->getKeyword("DIMENS");
dims[0] = dimensKeyword->getRecord(0)->getItem(0)->getInt(0);
dims[1] = dimensKeyword->getRecord(0)->getItem(1)->getInt(0);
dims[2] = dimensKeyword->getRecord(0)->getItem(2)->getInt(0);
} else if (newParserDeck->hasKeyword("SPECGRID")) {
Opm::DeckKeywordConstPtr specgridKeyword = newParserDeck->getKeyword("SPECGRID");
} else if (deck->hasKeyword("SPECGRID")) {
Opm::DeckKeywordConstPtr specgridKeyword = deck->getKeyword("SPECGRID");
dims[0] = specgridKeyword->getRecord(0)->getItem(0)->getInt(0);
dims[1] = specgridKeyword->getRecord(0)->getItem(1)->getInt(0);
dims[2] = specgridKeyword->getRecord(0)->getItem(2)->getInt(0);
@@ -251,8 +206,8 @@ namespace Opm
grdecl.dims[1] = dims[1];
grdecl.dims[2] = dims[2];
if (newParserDeck->hasKeyword("MAPAXES")) {
Opm::DeckKeywordConstPtr mapaxesKeyword = newParserDeck->getKeyword("MAPAXES");
if (deck->hasKeyword("MAPAXES")) {
Opm::DeckKeywordConstPtr mapaxesKeyword = deck->getKeyword("MAPAXES");
Opm::DeckRecordConstPtr mapaxesRecord = mapaxesKeyword->getRecord(0);
// memleak alert: here we need to make sure that C code
@@ -280,82 +235,17 @@ namespace Opm
} // anonymous namespace
// Construct tensor grid from deck.
void GridManager::initFromDeckTensorgrid(const Opm::EclipseGridParser& deck)
{
// Extract logical cartesian size.
std::vector<int> dims;
if (deck.hasField("DIMENS")) {
dims = deck.getIntegerValue("DIMENS");
} else if (deck.hasField("SPECGRID")) {
dims = deck.getSPECGRID().dimensions;
} else {
OPM_THROW(std::runtime_error, "Deck must have either DIMENS or SPECGRID.");
}
// Extract coordinates (or offsets from top, in case of z).
const std::vector<double>& dxv = deck.getFloatingPointValue("DXV");
const std::vector<double>& dyv = deck.getFloatingPointValue("DYV");
const std::vector<double>& dzv = deck.getFloatingPointValue("DZV");
std::vector<double> x = coordsFromDeltas(dxv);
std::vector<double> y = coordsFromDeltas(dyv);
std::vector<double> z = coordsFromDeltas(dzv);
// Check that number of cells given are consistent with DIMENS/SPECGRID.
if (dims[0] != int(dxv.size())) {
OPM_THROW(std::runtime_error, "Number of DXV data points do not match DIMENS or SPECGRID.");
}
if (dims[1] != int(dyv.size())) {
OPM_THROW(std::runtime_error, "Number of DYV data points do not match DIMENS or SPECGRID.");
}
if (dims[2] != int(dzv.size())) {
OPM_THROW(std::runtime_error, "Number of DZV data points do not match DIMENS or SPECGRID.");
}
// Extract top corner depths, if available.
const double* top_depths = 0;
std::vector<double> top_depths_vec;
if (deck.hasField("DEPTHZ")) {
const std::vector<double>& depthz = deck.getFloatingPointValue("DEPTHZ");
if (depthz.size() != x.size()*y.size()) {
OPM_THROW(std::runtime_error, "Incorrect size of DEPTHZ: " << depthz.size());
}
top_depths = &depthz[0];
} else if (deck.hasField("TOPS")) {
// We only support constant values for TOPS.
// It is not 100% clear how we best can deal with
// varying TOPS (stair-stepping grid, or not).
const std::vector<double>& tops = deck.getFloatingPointValue("TOPS");
if (std::count(tops.begin(), tops.end(), tops[0]) != int(tops.size())) {
OPM_THROW(std::runtime_error, "We do not support nonuniform TOPS, please use ZCORN/COORDS instead.");
}
top_depths_vec.resize(x.size()*y.size(), tops[0]);
top_depths = &top_depths_vec[0];
} else {
OPM_THROW(std::runtime_error,
"Could not find either TOPS or DEPTHZ keyword, "
"one of them is required for initialization with DXV/DYV/DZV.");
}
// Construct grid.
ug_ = create_grid_tensor3d(dxv.size(), dyv.size(), dzv.size(),
&x[0], &y[0], &z[0], top_depths);
if (!ug_) {
OPM_THROW(std::runtime_error, "Failed to construct grid.");
}
}
void GridManager::initFromDeckTensorgrid(Opm::DeckConstPtr newParserDeck)
void GridManager::initFromDeckTensorgrid(Opm::DeckConstPtr deck)
{
// Extract logical cartesian size.
std::array<int, 3> dims;
if (newParserDeck->hasKeyword("DIMENS")) {
Opm::DeckKeywordConstPtr dimensKeyword = newParserDeck->getKeyword("DIMENS");
if (deck->hasKeyword("DIMENS")) {
Opm::DeckKeywordConstPtr dimensKeyword = deck->getKeyword("DIMENS");
dims[0] = dimensKeyword->getRecord(0)->getItem(0)->getInt(0);
dims[1] = dimensKeyword->getRecord(0)->getItem(1)->getInt(0);
dims[2] = dimensKeyword->getRecord(0)->getItem(2)->getInt(0);
} else if (newParserDeck->hasKeyword("SPECGRID")) {
Opm::DeckKeywordConstPtr specgridKeyword = newParserDeck->getKeyword("SPECGRID");
} else if (deck->hasKeyword("SPECGRID")) {
Opm::DeckKeywordConstPtr specgridKeyword = deck->getKeyword("SPECGRID");
dims[0] = specgridKeyword->getRecord(0)->getItem(0)->getInt(0);
dims[1] = specgridKeyword->getRecord(0)->getItem(1)->getInt(0);
dims[2] = specgridKeyword->getRecord(0)->getItem(2)->getInt(0);
@@ -364,9 +254,9 @@ namespace Opm
}
// Extract coordinates (or offsets from top, in case of z).
const std::vector<double>& dxv = newParserDeck->getKeyword("DXV")->getSIDoubleData();
const std::vector<double>& dyv = newParserDeck->getKeyword("DYV")->getSIDoubleData();
const std::vector<double>& dzv = newParserDeck->getKeyword("DZV")->getSIDoubleData();
const std::vector<double>& dxv = deck->getKeyword("DXV")->getSIDoubleData();
const std::vector<double>& dyv = deck->getKeyword("DYV")->getSIDoubleData();
const std::vector<double>& dzv = deck->getKeyword("DZV")->getSIDoubleData();
std::vector<double> x = coordsFromDeltas(dxv);
std::vector<double> y = coordsFromDeltas(dyv);
std::vector<double> z = coordsFromDeltas(dzv);
@@ -385,17 +275,17 @@ namespace Opm
// Extract top corner depths, if available.
const double* top_depths = 0;
std::vector<double> top_depths_vec;
if (newParserDeck->hasKeyword("DEPTHZ")) {
const std::vector<double>& depthz = newParserDeck->getKeyword("DEPTHZ")->getSIDoubleData();
if (deck->hasKeyword("DEPTHZ")) {
const std::vector<double>& depthz = deck->getKeyword("DEPTHZ")->getSIDoubleData();
if (depthz.size() != x.size()*y.size()) {
OPM_THROW(std::runtime_error, "Incorrect size of DEPTHZ: " << depthz.size());
}
top_depths = &depthz[0];
} else if (newParserDeck->hasKeyword("TOPS")) {
} else if (deck->hasKeyword("TOPS")) {
// We only support constant values for TOPS.
// It is not 100% clear how we best can deal with
// varying TOPS (stair-stepping grid, or not).
const std::vector<double>& tops = newParserDeck->getKeyword("TOPS")->getSIDoubleData();
const std::vector<double>& tops = deck->getKeyword("TOPS")->getSIDoubleData();
if (std::count(tops.begin(), tops.end(), tops[0]) != int(tops.size())) {
OPM_THROW(std::runtime_error, "We do not support nonuniform TOPS, please use ZCORN/COORDS instead.");
}
opm/core/grid/GridManager.hpp
+4 -14
View File
@@ -30,9 +30,6 @@ struct grdecl;
namespace Opm
{
class EclipseGridParser;
/// This class manages an Opm::UnstructuredGrid in the sense that it
/// encapsulates creation and destruction of the grid.
/// The following grid types can be constructed:
@@ -45,10 +42,7 @@ namespace Opm
{
public:
/// Construct a 3d corner-point grid or tensor grid from a deck.
explicit GridManager(const Opm::EclipseGridParser& deck);
/// Construct a 3d corner-point grid or tensor grid from a deck.
explicit GridManager(Opm::DeckConstPtr newParserDeck);
explicit GridManager(Opm::DeckConstPtr deck);
/// Construct a grid from an EclipseState::EclipseGrid instance
explicit GridManager(Opm::EclipseGridConstPtr eclipseGrid);
@@ -74,14 +68,12 @@ namespace Opm
/// Destructor.
~GridManager();
void saveEGRID(const std::string& filename , const Opm::EclipseGridParser& deck);
/// Access the managed UnstructuredGrid.
/// The method is named similarly to c_str() in std::string,
/// to make it clear that we are returning a C-compatible struct.
const UnstructuredGrid* c_grid() const;
static void createGrdecl(Opm::DeckConstPtr newParserDeck, struct grdecl &grdecl);
static void createGrdecl(Opm::DeckConstPtr deck, struct grdecl &grdecl);
private:
// Disable copying and assignment.
@@ -89,11 +81,9 @@ namespace Opm
GridManager& operator=(const GridManager& other);
// Construct corner-point grid from deck.
void initFromDeckCornerpoint(const Opm::EclipseGridParser& deck);
void initFromDeckCornerpoint(Opm::DeckConstPtr newParserDeck);
void initFromDeckCornerpoint(Opm::DeckConstPtr deck);
// Construct tensor grid from deck.
void initFromDeckTensorgrid(const Opm::EclipseGridParser& deck);
void initFromDeckTensorgrid(Opm::DeckConstPtr newParserDeck);
void initFromDeckTensorgrid(Opm::DeckConstPtr deck);
// The managed UnstructuredGrid.
UnstructuredGrid* ug_;
opm/core/io/eclipse/CornerpointChopper.hpp
-38
View File
@@ -20,7 +20,6 @@
#ifndef OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
#define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
@@ -44,7 +43,6 @@ namespace Opm
{
public:
CornerPointChopper(const std::string& file)
: parser_(file, false)
{
Opm::ParserPtr parser(new Opm::Parser());
deck_ = parser->parseFile(file);
@@ -258,38 +256,6 @@ namespace Opm
filterIntegerField("SATNUM", new_SATNUM_);
}
/// Return a subparser with fields corresponding to the selected subset.
/// Note that the returned parser is NOT converted to SI, that must be done
/// by the user afterwards with the parser's convertToSI() method.
EclipseGridParser subparser()
{
if (parser_.hasField("FIELD") || parser_.hasField("LAB") || parser_.hasField("PVT-M")) {
OPM_THROW(std::runtime_error, "CornerPointChopper::subparser() cannot handle any eclipse unit system other than METRIC.");
}
EclipseGridParser sp;
std::shared_ptr<SPECGRID> sg(new SPECGRID);
for (int dd = 0; dd < 3; ++dd) {
sg->dimensions[dd] = new_dims_[dd];
}
sp.setSpecialField("SPECGRID", sg);
sp.setFloatingPointField("COORD", new_COORD_);
sp.setFloatingPointField("ZCORN", new_ZCORN_);
if (!new_ACTNUM_.empty()) sp.setIntegerField("ACTNUM", new_ACTNUM_);
if (!new_PORO_.empty()) sp.setFloatingPointField("PORO", new_PORO_);
if (!new_NTG_.empty()) sp.setFloatingPointField("NTG", new_NTG_);
if (!new_SWCR_.empty()) sp.setFloatingPointField("SWCR", new_SWCR_);
if (!new_SOWCR_.empty()) sp.setFloatingPointField("SOWCR", new_SOWCR_);
if (!new_PERMX_.empty()) sp.setFloatingPointField("PERMX", new_PERMX_);
if (!new_PERMY_.empty()) sp.setFloatingPointField("PERMY", new_PERMY_);
if (!new_PERMZ_.empty()) sp.setFloatingPointField("PERMZ", new_PERMZ_);
if (!new_SATNUM_.empty()) sp.setIntegerField("SATNUM", new_SATNUM_);
sp.computeUnits(); // Always METRIC, since that is default.
return sp;
}
/// Return a sub-deck with fields corresponding to the selected subset.
Opm::DeckConstPtr subDeck()
{
@@ -319,7 +285,6 @@ namespace Opm
specGridKw->addRecord(specGridRecord);
subDeck->addKeyword(specGridKw);
addDoubleKeyword_(subDeck, "COORD", /*dimension=*/"Length", new_COORD_);
addDoubleKeyword_(subDeck, "ZCORN", /*dimension=*/"Length", new_ZCORN_);
addIntKeyword_(subDeck, "ACTNUM", new_ACTNUM_);
@@ -331,10 +296,8 @@ namespace Opm
addDoubleKeyword_(subDeck, "PERMY", /*dimension=*/"Permeability", new_PERMY_);
addDoubleKeyword_(subDeck, "PERMZ", /*dimension=*/"Permeability", new_PERMZ_);
addIntKeyword_(subDeck, "SATNUM", new_SATNUM_);
return subDeck;
}
void writeGrdecl(const std::string& filename)
{
// Output new versions of SPECGRID, COORD, ZCORN, ACTNUM, PERMX, PORO, SATNUM.
@@ -363,7 +326,6 @@ namespace Opm
bool hasSOWCR() const {return !new_SOWCR_.empty(); }
private:
Opm::EclipseGridParser parser_;
Opm::DeckConstPtr deck_;
std::shared_ptr<Opm::UnitSystem> metricUnits_;
opm/core/io/eclipse/EclipseGridInspector.cpp
+16 -26
View File
@@ -51,40 +51,30 @@
namespace Opm
{
EclipseGridInspector::EclipseGridInspector(const EclipseGridParser& oldParser)
{
Opm::ParserConstPtr parser(new Opm::Parser());
newParserDeck_ = parser->parseFile(oldParser.deckFileName());
init_();
}
EclipseGridInspector::EclipseGridInspector(Opm::DeckConstPtr newParserDeck)
: newParserDeck_(newParserDeck)
EclipseGridInspector::EclipseGridInspector(Opm::DeckConstPtr deck)
: deck_(deck)
{
init_();
}
void EclipseGridInspector::init_()
{
if (!newParserDeck_->hasKeyword("COORD")) {
if (!deck_->hasKeyword("COORD")) {
OPM_THROW(std::runtime_error, "Needed field \"COORD\" is missing in file");
}
if (!newParserDeck_->hasKeyword("ZCORN")) {
if (!deck_->hasKeyword("ZCORN")) {
OPM_THROW(std::runtime_error, "Needed field \"ZCORN\" is missing in file");
}
if (newParserDeck_->hasKeyword("SPECGRID")) {
if (deck_->hasKeyword("SPECGRID")) {
Opm::DeckRecordConstPtr specgridRecord =
newParserDeck_->getKeyword("SPECGRID")->getRecord(0);
deck_->getKeyword("SPECGRID")->getRecord(0);
logical_gridsize_[0] = specgridRecord->getItem("NX")->getInt(0);
logical_gridsize_[1] = specgridRecord->getItem("NY")->getInt(0);
logical_gridsize_[2] = specgridRecord->getItem("NZ")->getInt(0);
} else if (newParserDeck_->hasKeyword("DIMENS")) {
} else if (deck_->hasKeyword("DIMENS")) {
Opm::DeckRecordConstPtr dimensRecord =
newParserDeck_->getKeyword("DIMENS")->getRecord(0);
deck_->getKeyword("DIMENS")->getRecord(0);
logical_gridsize_[0] = dimensRecord->getItem("NX")->getInt(0);
logical_gridsize_[1] = dimensRecord->getItem("NY")->getInt(0);
logical_gridsize_[2] = dimensRecord->getItem("NZ")->getInt(0);
@@ -107,13 +97,13 @@ std::pair<double,double> EclipseGridInspector::cellDips(int i, int j, int k) con
{
checkLogicalCoords(i, j, k);
const std::vector<double>& pillc =
newParserDeck_->getKeyword("COORD")->getSIDoubleData();
deck_->getKeyword("COORD")->getSIDoubleData();
int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1);
if (6*num_pillars != int(pillc.size())) {
throw std::runtime_error("Wrong size of COORD field.");
}
const std::vector<double>& z =
newParserDeck_->getKeyword("ZCORN")->getSIDoubleData();
deck_->getKeyword("ZCORN")->getSIDoubleData();
int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
if (8*num_cells != int(z.size())) {
throw std::runtime_error("Wrong size of ZCORN field");
@@ -217,13 +207,13 @@ double EclipseGridInspector::cellVolumeVerticalPillars(int i, int j, int k) cons
// Checking parameters and obtaining values from parser.
checkLogicalCoords(i, j, k);
const std::vector<double>& pillc =
newParserDeck_->getKeyword("COORD")->getSIDoubleData();
deck_->getKeyword("COORD")->getSIDoubleData();
int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1);
if (6*num_pillars != int(pillc.size())) {
throw std::runtime_error("Wrong size of COORD field.");
}
const std::vector<double>& z =
newParserDeck_->getKeyword("ZCORN")->getSIDoubleData();
deck_->getKeyword("ZCORN")->getSIDoubleData();
int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
if (8*num_cells != int(z.size())) {
throw std::runtime_error("Wrong size of ZCORN field");
@@ -281,12 +271,12 @@ void EclipseGridInspector::checkLogicalCoords(int i, int j, int k) const
std::array<double, 6> EclipseGridInspector::getGridLimits() const
{
if (! (newParserDeck_->hasKeyword("COORD") && newParserDeck_->hasKeyword("ZCORN") && newParserDeck_->hasKeyword("SPECGRID")) ) {
if (! (deck_->hasKeyword("COORD") && deck_->hasKeyword("ZCORN") && deck_->hasKeyword("SPECGRID")) ) {
throw std::runtime_error("EclipseGridInspector: Grid does not have SPECGRID, COORD, and ZCORN, can't find dimensions.");
}
std::vector<double> coord = newParserDeck_->getKeyword("COORD")->getSIDoubleData();
std::vector<double> zcorn = newParserDeck_->getKeyword("ZCORN")->getSIDoubleData();
std::vector<double> coord = deck_->getKeyword("COORD")->getSIDoubleData();
std::vector<double> zcorn = deck_->getKeyword("ZCORN")->getSIDoubleData();
double xmin = +DBL_MAX;
double xmax = -DBL_MAX;
@@ -335,7 +325,7 @@ std::array<int, 3> EclipseGridInspector::gridSize() const
std::array<double, 8> EclipseGridInspector::cellZvals(int i, int j, int k) const
{
// Get the zcorn field.
const std::vector<double>& z = newParserDeck_->getKeyword("ZCORN")->getSIDoubleData();
const std::vector<double>& z = deck_->getKeyword("ZCORN")->getSIDoubleData();
int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
if (8*num_cells != int(z.size())) {
throw std::runtime_error("Wrong size of ZCORN field");
opm/core/io/eclipse/EclipseGridInspector.hpp
+4 -11
View File
@@ -38,8 +38,6 @@
#include <vector>
#include <array>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/io/eclipse/SpecialEclipseFields.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
namespace Opm
@@ -48,9 +46,8 @@ namespace Opm
/**
@brief A class for inspecting the contents of an eclipse file.
Given an EclipseGridParser that has successfully read en Eclipse
.grdecl-file, this class may be used to answer certain queries about
its contents.
Given an Eclipse deck, this class may be used to answer certain
queries about its contents.
@author Atgeirr F. Rasmussen <atgeirr@sintef.no>
@date 2008/06/02 09:46:08
@@ -60,11 +57,7 @@ class EclipseGridInspector
public:
/// Constructor taking a parser as argument.
/// The parser must already have read an Eclipse file.
EclipseGridInspector(const EclipseGridParser& parser);
/// Constructor taking a parser as argument.
/// The parser must already have read an Eclipse file.
EclipseGridInspector(Opm::DeckConstPtr newParserDeck);
EclipseGridInspector(Opm::DeckConstPtr deck);
/// Assuming that the pillars are vertical, compute the
/// volume of the cell given by logical coordinates (i, j, k).
@@ -98,7 +91,7 @@ public:
std::array<double, 8> cellZvals(int i, int j, int k) const;
private:
Opm::DeckConstPtr newParserDeck_;
Opm::DeckConstPtr deck_;
int logical_gridsize_[3];
void init_();
void checkLogicalCoords(int i, int j, int k) const;
opm/core/io/eclipse/EclipseGridParser.cpp
-1215
View File
File diff suppressed because it is too large. Load diff
opm/core/io/eclipse/EclipseGridParser.hpp
-289
View File
@@ -1,289 +0,0 @@
//===========================================================================
//
// File: EclipseGridParser.h
//
// Created: Wed Dec 5 17:05:13 2007
//
// Author: Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
//
// Revision: $Id: EclipseGridParser.h,v 1.3 2008/08/18 14:16:13 atgeirr Exp $
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLIPSEGRIDPARSER_HEADER
#define OPM_ECLIPSEGRIDPARSER_HEADER
#include <iosfwd>
#include <string>
#include <vector>
#include <map>
#include <set>
#include <memory>
#include <opm/core/io/eclipse/SpecialEclipseFields.hpp>
#include <opm/core/io/eclipse/EclipseUnits.hpp>
#include <opm/core/utility/Factory.hpp>
#include <opm/core/grid/cornerpoint_grid.h>
#ifdef HAVE_ERT
#include <ert/ecl/ecl_kw.h>
#include <ert/ecl/ecl_grid.h>
#endif
namespace Opm
{
/**
@brief A class for reading and parsing all fields of an eclipse file.
This object is constructed using an Eclipse .grdecl-file. All data
fields are extracted upon construction and written to vector data
structures, which can then be read out afterwards via
convenience functions.
There is also a convenience function to easily check which fields
were successfully parsed.
@author Atgeirr F. Rasmussen <atgeirr@sintef.no>
@date 2007/12/06 13:00:03
*/
enum FieldType {
Integer,
FloatingPoint,
Timestepping,
SpecialField,
IgnoreWithData,
IgnoreNoData,
Include,
Import,
Unknown
};
class EclipseGridParser
{
public:
/// Default constructor.
EclipseGridParser();
/// Constructor taking an eclipse filename. Unless the second
/// argument 'convert_to_SI' is false, all fields will be
/// converted to SI units.
explicit EclipseGridParser(const std::string& filename, bool convert_to_SI = true);
static FieldType classifyKeyword(const std::string& keyword);
static bool readKeyword(std::istream& is, std::string& keyword);
// temporary measure to get the file name passed to the
// constuctor above
const std::string &deckFileName() const
{ return deckFileName_; };
/// Read the given stream, overwriting any previous data. Unless
/// the second argument 'convert_to_SI' is false, all fields will
/// be converted to SI units.
void read(std::istream& is, bool convert_to_SI = true);
/// Convert all data to SI units, according to unit category as
/// specified in the eclipse file (METRIC, FIELD etc.).
/// It is unnecessary to call this if constructed with
/// 'convert_to_SI' equal to true, but it is not an error.
void convertToSI();
/// Returns true if the given keyword corresponds to a field that
/// was found in the file.
bool hasField(const std::string& keyword) const;
/// Returns true if all the given keywords correspond to fields
/// that were found in the file.
bool hasFields(const std::vector<std::string>& keywords) const;
/// The keywords/fields found in the file.
std::vector<std::string> fieldNames() const;
/// Returns the number of distinct epochs found in the deck SCHEDULE.
int numberOfEpochs() const;
/// Sets the current epoch.
/// Valid arguments are in [0, ..., numberOfEpochs() - 1].
/// After reading, current epoch always starts at 0.
void setCurrentEpoch(int epoch);
/// Returns the start_date_
boost::gregorian::date getStartDate() const;
/// Returns a reference to a vector containing the values
/// corresponding to the given integer keyword.
const std::vector<int>& getIntegerValue(const std::string& keyword) const;
/// Returns a reference to a vector containing the values
/// corresponding to the given floating-point keyword.
const std::vector<double>& getFloatingPointValue(const std::string& keyword) const;
/// Returns a reference to a vector containing the values
/// corresponding to the given keyword of a type only known
/// indirectly (through a template)
///
/// \tparam T Type of the keyword's value. Currently only int
/// and double are supported.
template <typename T>
const std::vector<T>& getValue(const std::string& keyword) const;
typedef std::shared_ptr<SpecialBase> SpecialFieldPtr;
/// Returns a reference to a vector containing pointers to the values
/// corresponding to the given keyword when the values are not only integers
/// or floats.
const SpecialFieldPtr getSpecialValue(const std::string& keyword) const;
// This macro implements support for a special field keyword. It requires that a subclass
// of SpecialBase exists, that has the same name as the keyword.
// After using SPECIAL_FIELD(KEYWORD), the public member getKEYWORD will be available.
#define SPECIAL_FIELD(keyword) \
private: \
struct X##keyword { X##keyword() { Factory<SpecialBase>::addCreator<keyword>(#keyword); } }; \
X##keyword x##keyword; \
public: \
const keyword& get##keyword() const \
{ return dynamic_cast<const keyword&>(*getSpecialValue(#keyword)); }
// Support for special fields.
SPECIAL_FIELD(SPECGRID)
SPECIAL_FIELD(FAULTS)
SPECIAL_FIELD(MULTFLT)
SPECIAL_FIELD(TITLE)
SPECIAL_FIELD(START)
SPECIAL_FIELD(DATES)
SPECIAL_FIELD(DENSITY)
SPECIAL_FIELD(PVDG)
SPECIAL_FIELD(PVDO)
SPECIAL_FIELD(PVTG)
SPECIAL_FIELD(PVTO)
SPECIAL_FIELD(PVTW)
SPECIAL_FIELD(SGOF)
SPECIAL_FIELD(SWOF)
SPECIAL_FIELD(ROCK)
SPECIAL_FIELD(ROCKTAB)
SPECIAL_FIELD(WELSPECS)
SPECIAL_FIELD(COMPDAT)
SPECIAL_FIELD(WCONINJE)
SPECIAL_FIELD(WCONPROD)
SPECIAL_FIELD(WELTARG)
SPECIAL_FIELD(WELOPEN)
SPECIAL_FIELD(EQUIL)
SPECIAL_FIELD(PVCDO)
SPECIAL_FIELD(TSTEP)
SPECIAL_FIELD(PLYVISC)
SPECIAL_FIELD(PLYROCK)
SPECIAL_FIELD(PLYADS)
SPECIAL_FIELD(PLYMAX)
SPECIAL_FIELD(PLYSHEAR)
SPECIAL_FIELD(TLMIXPAR)
SPECIAL_FIELD(WPOLYMER)
SPECIAL_FIELD(GRUPTREE)
SPECIAL_FIELD(GCONINJE)
SPECIAL_FIELD(GCONPROD)
SPECIAL_FIELD(WGRUPCON)
SPECIAL_FIELD(ENDSCALE)
SPECIAL_FIELD(SCALECRS)
SPECIAL_FIELD(ENPTVD)
SPECIAL_FIELD(ENKRVD)
// The following fields only have a dummy implementation
// that allows us to ignore them.
SPECIAL_FIELD(SWFN)
SPECIAL_FIELD(SOF2)
SPECIAL_FIELD(TUNING)
#undef SPECIAL_FIELD
/// Sets an integer field to have a particular value.
void setIntegerField(const std::string& keyword, const std::vector<int>& field);
/// Sets a floating point field to have a particular value.
void setFloatingPointField(const std::string& keyword, const std::vector<double>& field);
/// Sets a special field to have a particular value.
void setSpecialField(const std::string& keyword, SpecialFieldPtr field);
/// Compute the units used by the deck, depending on the presence
/// of keywords such as METRIC, FIELD etc. It is an error to call
/// this after conversion to SI has taken place.
void computeUnits();
/// The units specified by the eclipse file read.
const EclipseUnits& units() const;
struct grdecl get_grdecl() const;
/// Save grid parts of deck in EGRID format.
void saveEGRID(const std::string & filename, int num_cells , const int * global_cell) const;
#ifdef HAVE_ERT
//void saveEGRID_INIT( const std::string& output_dir , const std::string& basename, bool fmt_file = false);
void saveINIT( const std::string & filename , const ecl_grid_type * ecl_grid);
void saveEGRID_INIT( const std::string& output_dir , const std::string& basename, bool fmt_file);
ecl_grid_type * newGrid( );
#endif
private:
#ifdef HAVE_ERT
ecl_kw_type * newEclKW(const std::string &keyword , ecl_type_enum ecl_type) const;
void save_kw( fortio_type * fortio , const std::string & kw , ecl_type_enum ecl_type);
#endif
SpecialFieldPtr createSpecialField(std::istream& is, const std::string& fieldname);
SpecialFieldPtr cloneSpecialField(const std::string& fieldname,
const std::shared_ptr<SpecialBase> original);
void readImpl(std::istream& is);
void getNumericErtFields(const std::string& filename);
std::string deckFileName_;
std::string directory_;
std::map<std::string, std::vector<int> > integer_field_map_;
std::map<std::string, std::vector<double> > floating_field_map_;
// std::map<std::string, SpecialFieldPtr> special_field_map_;
std::set<std::string> ignored_fields_;
EclipseUnits units_;
// For SCHEDULE handling.
enum ReadingMode { Regular, Timesteps };
ReadingMode current_reading_mode_;
boost::gregorian::date start_date_;
double current_time_days_;
int current_epoch_;
typedef std::map<std::string, SpecialFieldPtr> SpecialMap;
std::vector<SpecialMap> special_field_by_epoch_;
};
} // namespace Opm
#endif // OPM_ECLIPSEGRIDPARSER_HEADER
opm/core/io/eclipse/EclipseGridParserHelpers.hpp
-465
View File
@@ -1,465 +0,0 @@
//===========================================================================
//
// File: EclipseGridParserHelpers.hpp
//
// Created: Tue Dec 22 11:35:32 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B�rd Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLIPSEGRIDPARSERHELPERS_HEADER
#define OPM_ECLIPSEGRIDPARSERHELPERS_HEADER
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <boost/date_time/gregorian/gregorian.hpp>
#include <limits>
#include <string>
#include <istream>
#include <vector>
#include <iostream>
namespace Opm
{
namespace
{
inline std::istream& ignoreLine(std::istream& is)
{
is.ignore(std::numeric_limits<int>::max(), '\n');
return is;
}
inline std::istream& ignoreSlashLine(std::istream& is)
{
is.ignore(std::numeric_limits<int>::max(), '/');
is.ignore(std::numeric_limits<int>::max(), '\n');
return is;
}
inline std::istream& ignoreWhitespace(std::istream& is)
{
// Getting the character type facet for is()
// We use the classic (i.e. C) locale.
const std::ctype<char>& ct = std::use_facet< std::ctype<char> >(std::locale::classic());
char c;
while (is.get(c)) {
if (!ct.is(std::ctype_base::space, c)) {
is.putback(c);
break;
}
}
return is;
}
inline std::string readString(std::istream& is)
{
std::string string_candidate;
is >> string_candidate;
const char quote('\'');
int beg = string_candidate[0] == quote ? 1 : 0;
int len = string_candidate[0] == quote ? string_candidate.size() - 2 : string_candidate.size();
return string_candidate.substr(beg, len);
}
// Reads data until '/' or an error is encountered.
template<typename T>
inline void readVectorData(std::istream& is, std::vector<T>& data)
{
data.clear();
while (is) {
T candidate;
is >> candidate;
if (is.rdstate() & std::ios::failbit) {
is.clear(is.rdstate() & ~std::ios::failbit);
is >> ignoreWhitespace;
char dummy;
is >> dummy;
if (dummy == '/') {
is >> ignoreLine;
break;
} else if (dummy == '-') { // "comment test"
is >> ignoreLine; // This line is a comment
} else {
char buffer[1000];
is.getline(buffer, sizeof(buffer));
std::cout << buffer<<std::endl;
OPM_THROW(std::runtime_error, "Encountered format error while reading data values. Value = " << dummy);
}
} else {
if (is.peek() == int('*')) {
is.ignore(); // ignore the '*'
int multiplier = int(candidate);
is >> candidate;
data.insert(data.end(), multiplier, candidate);
} else {
data.push_back(candidate);
}
}
}
if (!is) {
OPM_THROW(std::runtime_error, "Encountered error while reading data values.");
}
}
// Reads data items of type T. Not more than 'max_values' items.
// Asterisks may be used to signify 'repeat counts'. 5*3.14 will
// insert 3.14 five times. Asterisk followed by a space is used to
// signify default values. n* will default n consecutive quantities.
template<class Vec>
inline int readDefaultedVectorData(std::istream& is, Vec& data, int max_values)
{
assert(int(data.size()) >= max_values);
const std::ctype<char>& ct = std::use_facet< std::ctype<char> >(std::locale::classic());
int num_values = 0;
while (is) {
typename Vec::value_type candidate;
is >> candidate;
if (is.rdstate() & std::ios::failbit) {
is.clear(is.rdstate() & ~std::ios::failbit);
std::string dummy;
is >> dummy;
if (dummy == "/") {
is >> ignoreLine;
break;
} else if (dummy[0] == '-') { // "comment test"
is >> ignoreLine; // This line is a comment
} else {
OPM_THROW(std::runtime_error, "Encountered format error while reading data values. Value = " << dummy);
}
} else {
if (is.peek() == int('*')) {
is.ignore(); // ignore the '*'
int multiplier = (int)candidate;
if (ct.is(std::ctype_base::space, is.peek())) {
num_values += multiplier; // Use default value(s)
} else {
is >> candidate; // Use candidate 'multipler' times
for (int i=0; i<multiplier; ++i, ++num_values) {
data[num_values] = candidate;
}
}
} else {
data[num_values] = candidate;
++num_values;
}
}
if (num_values >= max_values) {
//is >> ignoreLine;
break;
}
}
if (!is) {
OPM_THROW(std::runtime_error, "Encountered error while reading data values.");
}
return num_values;
}
// Keywords SGOF and SWOF. Reads data until '/' or an error is encountered.
// Default values represented by 1* is replaced by -1. Use linear interpolation
// outside this function to replace -1.
template<typename T>
inline void readRelPermTable(std::istream& is, std::vector<T>& data)
{
data.clear();
while (is) {
T candidate;
is >> candidate;
if (is.rdstate() & std::ios::failbit) {
is.clear(is.rdstate() & ~std::ios::failbit);
std::string dummy;
is >> dummy;
if (dummy == "/") {
is >> ignoreLine;
break;
} else if (dummy[0] == '-') { // "comment test"
is >> ignoreLine; // This line is a comment
} else {
OPM_THROW(std::runtime_error, "Encountered format error while reading data values. Value = " << dummy);
}
} else {
if (is.peek() == int('*')) {
is.ignore(); // ignore the '*'
assert(int(candidate) == 1);
data.push_back(-1); // Set new flag for interpolation.
} else {
data.push_back(candidate);
}
}
}
if (!is) {
OPM_THROW(std::runtime_error, "Encountered error while reading data values.");
}
}
// Returns month number 1-12. Returns 0 if illegal month name.
inline int getMonthNumber(const std::string& month_name)
{
const int num_months = 12;
std::string months[num_months] = {"JAN", "FEB", "MAR", "APR", "MAY", "JUN",
"JLY", "AUG", "SEP", "OCT", "NOV", "DEC"};
if (month_name == "JUL") {
return 7; // "JUL" is an acceptable alternative to 'JLY'
}
int m = 0;
for (int i=0; i<num_months; ++i) {
if (month_name == months[i]) {
m = i+1;
break;
}
}
return m;
}
inline boost::gregorian::date readDate(std::istream& is)
{
while (is.peek() == int('-')) {
is >> ignoreLine; // This line is a comment
}
int day, year;
std::string month_name;
is >> day;
month_name = readString(is);
is >> year;
ignoreSlashLine(is);
int month = getMonthNumber(month_name);
return boost::gregorian::date(boost::gregorian::greg_year(year),
boost::gregorian::greg_month(month),
boost::gregorian::greg_day(day));
}
// Get first character after whitespace and comments, and decide
// next action.
// NB! Will fail for negative number in column one.
inline int next_action(std::istream& is)
{
int task(0); // 0:continue 1:return 2:throw
const std::ctype<char>& ct =
std::use_facet< std::ctype<char> >(std::locale::classic());
while (!is.eof()) {
is >> ignoreWhitespace;
char c;
is.get(c);
if (is.eof()) {
return task;
}
is.putback(c);
if (ct.is(std::ctype_base::digit, c) || c== '.') {
task = 0; // Decimal digit. Read more data.
break;
}
if (ct.is(std::ctype_base::alpha, c)) {
task = 1; // Alphabetic char. Read next keyword.
break;
}
if (c == '-') {
is >> ignoreLine; // This line is a comment
} else {
task = 2; // Read error. Unexpected character.
break;
}
}
return task;
}
// Reads keywords PVTG, PVTO
typedef std::vector<std::vector<std::vector<double> > > table_t;
inline void readPvtTable(std::istream& is, table_t& pvt_table,
const std::string& field_name)
{
const std::ctype<char>& ct =
std::use_facet< std::ctype<char> >(std::locale::classic());
std::vector<double> record;
std::vector<std::vector<double> > table;
while (!is.eof()) {
record.clear();
readVectorData(is, record);
table.push_back(record);
while (!is.eof()) {
is >> ignoreWhitespace;
char c;
is.get(c);
if (is.eof()) {
// Reached end of file, we should have pushed
// the last table, and emptied it. If not,
// we have an error.
if (!table.empty()) {
OPM_THROW(std::runtime_error, "Reached EOF while still building PVT table. Missing end-of-table (slash)?");
}
return;
}
is.putback(c);
if (ct.is(std::ctype_base::digit, c) || c== '.') {
break; // Decimal digit. Read more records.
}
if (ct.is(std::ctype_base::alpha, c)) {
return; // Alphabetic char. Read next keyword.
}
if (c == '-') {
is >> ignoreLine; // This line is a comment
continue;
}
if (c == '/') {
is >> ignoreLine;
pvt_table.push_back(table);
table.clear();
} else {
std::ostringstream oss;
oss << "Error reading " << field_name
<< ". Next character is " << (char)is.peek();
OPM_THROW(std::runtime_error, oss.str());
}
}
}
}
// Reads keywords PVDG, PVDO, ROCKTAB
inline void readPvdTable(std::istream& is, table_t& pvd_table,
const std::string& field_name, int ncol)
{
std::vector<double> record;
std::vector<std::vector<double> > table(ncol);
while (!is.eof()) {
record.clear();
readVectorData(is, record);
const int rec_size = record.size()/ncol;
for (int k=0; k<ncol; ++k) {
table[k].resize(rec_size);
}
for (int i=0, n=-1; i<rec_size; ++i) {
for (int k=0; k<ncol; ++k) {
table[k][i] = record[++n];
}
}
pvd_table.push_back(table);
int action = next_action(is); // 0:continue 1:return 2:throw
if (action == 1) {
return; // Alphabetic char. Read next keyword.
} else if (action == 2) {
std::ostringstream oss;
oss << "Error reading " << field_name
<< ". Next character is " << (char)is.peek();
OPM_THROW(std::runtime_error, oss.str());
}
}
}
// Replace default values -1 by linear interpolation
inline void insertDefaultValues(std::vector<std::vector<double> >& table,
int ncol,
double defaultOut = 0.0,
bool defaultInterpolation = true)
{
const int sz = table[0].size();
for (int k=1; k<ncol; ++k) {
std::vector<int> indx;
std::vector<double> x;
for (int i=0; i<sz; ++i) {
if (table[k][i] == -1) {
indx.push_back(i);
x.push_back(table[0][i]);
}
}
if (!indx.empty()) {
std::vector<double> xv, yv;
for (int i=0; i<sz; ++i) {
if (table[k][i] != -1) {
xv.push_back(table[0][i]);
yv.push_back(table[k][i]);
}
}
if (xv.empty()) {
// Nothing specified, the entire column is defaulted.
// We insert default value.
for (int i=0; i<int(indx.size()); ++i) {
table[k][indx[i]] = defaultOut;
}
} else if (defaultInterpolation) {
// Interpolate
for (int i=0; i<int(indx.size()); ++i) {
table[k][indx[i]] = linearInterpolation(xv, yv, x[i]);
}
} else {
// Interpolate
for (int i=0; i<int(indx.size()); ++i) {
table[k][indx[i]] = linearInterpolation(xv, yv, x[i]);
}
}
}
}
}
// Reads keywords SGOF and SWOF
inline void readSGWOF(std::istream& is, table_t& relperm_table,
const std::string& field_name, int ncol)
{
std::vector<double> record;
std::vector<std::vector<double> > table(ncol);
while (!is.eof()) {
record.clear();
readRelPermTable(is, record);
const int rec_size = record.size()/ncol;
for (int k=0; k<ncol; ++k) {
table[k].resize(rec_size);
}
for (int i=0, n=-1; i<rec_size; ++i) {
for (int k=0; k<ncol; ++k) {
table[k][i] = record[++n];
}
}
insertDefaultValues(table, ncol);
relperm_table.push_back(table);
int action = next_action(is); // 0:continue 1:return 2:throw
if (action == 1) {
return; // Alphabetic char. Read next keyword.
} else if (action == 2) {
std::ostringstream oss;
oss << "Error reading " << field_name
<< ". Next character is " << (char)is.peek();
OPM_THROW(std::runtime_error, oss.str());
}
}
}
} // anon namespace
} // namespace Opm
#endif // OPM_ECLIPSEGRIDPARSERHELPERS_HEADER
opm/core/io/eclipse/EclipseWriter.cpp
+67 -65
View File
@@ -24,6 +24,8 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/cpgpreprocess/preprocess.h>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
@@ -413,20 +415,20 @@ private:
};
/// Get dimensions of the grid from the parse of the input file
std::vector <int> parserDim (Opm::DeckConstPtr newParserDeck) {
std::vector <int> parserDim (Opm::DeckConstPtr deck) {
std::vector<int> dim(/* n = */ 3);
// dimensions explicitly given
if (newParserDeck->hasKeyword("SPECGRID")) {
SpecgridWrapper specgrid(newParserDeck->getKeyword("SPECGRID"));
if (deck->hasKeyword("SPECGRID")) {
SpecgridWrapper specgrid(deck->getKeyword("SPECGRID"));
dim = specgrid.numBlocksVector();
}
// dimensions implicitly given by number of deltas
else if (newParserDeck->hasKeyword("DXV")) {
assert(newParserDeck->hasKeyword("DYV"));
assert(newParserDeck->hasKeyword("DZV"));
dim[0] = newParserDeck->getKeyword("DXV")->getRawDoubleData().size();
dim[1] = newParserDeck->getKeyword("DYV")->getRawDoubleData().size();
dim[2] = newParserDeck->getKeyword("DZV")->getRawDoubleData().size();
else if (deck->hasKeyword("DXV")) {
assert(deck->hasKeyword("DYV"));
assert(deck->hasKeyword("DZV"));
dim[0] = deck->getKeyword("DXV")->getRawDoubleData().size();
dim[1] = deck->getKeyword("DYV")->getRawDoubleData().size();
dim[2] = deck->getKeyword("DZV")->getRawDoubleData().size();
}
else {
OPM_THROW(std::runtime_error,
@@ -461,9 +463,9 @@ struct EclipseRestart : public EclipseHandle <ecl_rst_file_type> {
void writeHeader (const SimulatorTimer& timer,
int outputStepIdx,
const PhaseUsage uses,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
const int num_active_cells) {
const std::vector <int> dim = parserDim (newParserDeck);
const std::vector <int> dim = parserDim (deck);
ecl_rst_file_fwrite_header (*this,
outputStepIdx,
timer.currentPosixTime(),
@@ -505,29 +507,29 @@ private:
*/
struct EclipseWriterGrid : public EclipseHandle <ecl_grid_type> {
/// Create a grid based on the keywords available in input file
static EclipseWriterGrid make (Opm::DeckConstPtr newParserDeck,
static EclipseWriterGrid make (Opm::DeckConstPtr deck,
int number_of_cells,
const int* cart_dims,
const int* global_cell)
{
if (newParserDeck->hasKeyword("DXV")) {
if (deck->hasKeyword("DXV")) {
// make sure that the DYV and DZV keywords are present if the
// DXV keyword is used in the deck...
assert(newParserDeck->hasKeyword("DYV"));
assert(newParserDeck->hasKeyword("DZV"));
assert(deck->hasKeyword("DYV"));
assert(deck->hasKeyword("DZV"));
const auto& dxv = newParserDeck->getKeyword("DXV")->getSIDoubleData();
const auto& dyv = newParserDeck->getKeyword("DYV")->getSIDoubleData();
const auto& dzv = newParserDeck->getKeyword("DZV")->getSIDoubleData();
const auto& dxv = deck->getKeyword("DXV")->getSIDoubleData();
const auto& dyv = deck->getKeyword("DYV")->getSIDoubleData();
const auto& dzv = deck->getKeyword("DZV")->getSIDoubleData();
return EclipseWriterGrid (dxv, dyv, dzv);
}
else if (newParserDeck->hasKeyword("ZCORN")) {
else if (deck->hasKeyword("ZCORN")) {
struct grdecl g;
GridManager::createGrdecl(newParserDeck, g);
GridManager::createGrdecl(deck, g);
auto coordData = getAllSiDoubles_(newParserDeck->getKeyword(COORD_KW));
auto zcornData = getAllSiDoubles_(newParserDeck->getKeyword(ZCORN_KW));
auto coordData = getAllSiDoubles_(deck->getKeyword(COORD_KW));
auto zcornData = getAllSiDoubles_(deck->getKeyword(ZCORN_KW));
EclipseKeyword<float> coord_kw (COORD_KW, coordData);
EclipseKeyword<float> zcorn_kw (ZCORN_KW, zcornData);
@@ -538,7 +540,7 @@ struct EclipseWriterGrid : public EclipseHandle <ecl_grid_type> {
EclipseKeyword<float> mapaxes_kw (MAPAXES_KW);
if (g.mapaxes) {
auto mapaxesData = getAllSiDoubles_(newParserDeck->getKeyword(MAPAXES_KW));
auto mapaxesData = getAllSiDoubles_(deck->getKeyword(MAPAXES_KW));
mapaxes_kw = std::move (EclipseKeyword<float> (MAPAXES_KW, mapaxesData));
}
@@ -638,13 +640,13 @@ struct EclipseInit : public EclipseHandle <fortio_type> {
const int* cart_dims,
const int* global_cell,
const SimulatorTimer& timer,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
const PhaseUsage uses)
{
auto dataField = getAllSiDoubles_(newParserDeck->getKeyword(PORO_KW));
auto dataField = getAllSiDoubles_(deck->getKeyword(PORO_KW));
restrictToActiveCells_(dataField, number_of_cells, global_cell);
EclipseWriterGrid eclGrid = EclipseWriterGrid::make (newParserDeck, number_of_cells,
EclipseWriterGrid eclGrid = EclipseWriterGrid::make (deck, number_of_cells,
cart_dims, global_cell);
EclipseKeyword<float> poro (PORO_KW, dataField);
@@ -697,9 +699,9 @@ struct EclipseSummary : public EclipseHandle <ecl_sum_type> {
EclipseSummary (const std::string& outputDir,
const std::string& baseName,
const SimulatorTimer& timer,
Opm::DeckConstPtr newParserDeck)
Opm::DeckConstPtr deck)
: EclipseHandle <ecl_sum_type> (
alloc_writer (outputDir, baseName, timer, newParserDeck),
alloc_writer (outputDir, baseName, timer, deck),
ecl_sum_free) { }
typedef std::unique_ptr <EclipseWellReport> var_t;
@@ -718,7 +720,7 @@ struct EclipseSummary : public EclipseHandle <ecl_sum_type> {
}
// add rate variables for each of the well in the input file
void addWells (Opm::DeckConstPtr newParserDeck,
void addWells (Opm::DeckConstPtr deck,
const PhaseUsage& uses);
// no inline implementation of this since it depends on the
@@ -753,11 +755,11 @@ private:
static ecl_sum_type* alloc_writer (const std::string& outputDir,
const std::string& baseName,
const SimulatorTimer& timer,
Opm::DeckConstPtr newParserDeck) {
Opm::DeckConstPtr deck) {
boost::filesystem::path casePath (outputDir);
casePath /= boost::to_upper_copy (baseName);
const std::vector <int> dim = parserDim (newParserDeck);
const std::vector <int> dim = parserDim (deck);
return ecl_sum_alloc_writer (casePath.string ().c_str (),
false, /* formatted */
true, /* unified */
@@ -777,7 +779,7 @@ private:
struct EclipseWellReport : public EclipseHandle <smspec_node_type> {
protected:
EclipseWellReport (const EclipseSummary& summary, /* section to add to */
Opm::DeckConstPtr newParserDeck, /* well names */
Opm::DeckConstPtr deck, /* well names */
int whichWell, /* index of well line */
PhaseUsage uses, /* phases present */
BlackoilPhases::PhaseIndex phase, /* oil, water or gas */
@@ -789,7 +791,7 @@ protected:
varName (phase,
type,
aggregation).c_str (),
wellName (newParserDeck, whichWell).c_str (),
wellName (deck, whichWell).c_str (),
/* num = */ 0,
unit.c_str(),
/* defaultValue = */ 0.))
@@ -818,10 +820,10 @@ private:
const double sign_;
/// Get the name associated with this well
std::string wellName (Opm::DeckConstPtr newParserDeck,
std::string wellName (Opm::DeckConstPtr deck,
int whichWell)
{
Opm::WelspecsWrapper welspecs(newParserDeck->getKeyword("WELSPECS"));
Opm::WelspecsWrapper welspecs(deck->getKeyword("WELSPECS"));
return welspecs.wellName(whichWell);
}
@@ -873,13 +875,13 @@ protected:
/// Monitors the rate given by a well.
struct EclipseWellRate : public EclipseWellReport {
EclipseWellRate (const EclipseSummary& summary,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
newParserDeck,
deck,
whichWell,
uses,
phase,
@@ -897,13 +899,13 @@ struct EclipseWellRate : public EclipseWellReport {
/// Monitors the total production in a well.
struct EclipseWellTotal : public EclipseWellReport {
EclipseWellTotal (const EclipseSummary& summary,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
newParserDeck,
deck,
whichWell,
uses,
phase,
@@ -939,13 +941,13 @@ private:
/// Monitors the bottom hole pressure in a well.
struct EclipseWellBhp : public EclipseWellReport {
EclipseWellBhp (const EclipseSummary& summary,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
newParserDeck,
deck,
whichWell,
uses,
phase,
@@ -982,12 +984,12 @@ EclipseSummary::writeTimeStep (const SimulatorTimer& timer,
static WellType WELL_TYPES[] = { INJECTOR, PRODUCER };
inline void
EclipseSummary::addWells (Opm::DeckConstPtr newParserDeck,
EclipseSummary::addWells (Opm::DeckConstPtr deck,
const PhaseUsage& uses) {
// TODO: Only create report variables that are requested with keywords
// (e.g. "WOPR") in the input files, and only for those wells that are
// mentioned in those keywords
Opm::DeckKeywordConstPtr welspecsKeyword = newParserDeck->getKeyword("WELSPECS");
Opm::DeckKeywordConstPtr welspecsKeyword = deck->getKeyword("WELSPECS");
const int numWells = welspecsKeyword->size();
for (int phaseCounter = 0;
phaseCounter != BlackoilPhases::MaxNumPhases;
@@ -1006,7 +1008,7 @@ EclipseSummary::addWells (Opm::DeckConstPtr newParserDeck,
// W{O,G,W}{I,P}R
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellRate (*this,
newParserDeck,
deck,
whichWell,
uses,
phase,
@@ -1014,7 +1016,7 @@ EclipseSummary::addWells (Opm::DeckConstPtr newParserDeck,
// W{O,G,W}{I,P}T
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellTotal (*this,
newParserDeck,
deck,
whichWell,
uses,
phase,
@@ -1036,7 +1038,7 @@ EclipseSummary::addWells (Opm::DeckConstPtr newParserDeck,
}
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellBhp (*this,
newParserDeck,
deck,
whichWell,
uses,
phase,
@@ -1055,7 +1057,7 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer)
return;
}
/* Grid files */
EclipseWriterGrid eclGrid = EclipseWriterGrid::make (newParserDeck_, number_of_cells_,
EclipseWriterGrid eclGrid = EclipseWriterGrid::make (deck_, number_of_cells_,
cart_dims_, global_cell_);
eclGrid.write (outputDir_, baseName_, /*stepIdx=*/0);
@@ -1064,35 +1066,35 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer)
cart_dims_,
global_cell_,
timer,
newParserDeck_,
deck_,
uses_);
if (newParserDeck_->hasKeyword("PERM")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERM"));
if (deck_->hasKeyword("PERM")) {
auto data = getAllSiDoubles_(deck_->getKeyword("PERM"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERM", data);
}
if (newParserDeck_->hasKeyword("PERMX")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMX"));
if (deck_->hasKeyword("PERMX")) {
auto data = getAllSiDoubles_(deck_->getKeyword("PERMX"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMX", data);
}
if (newParserDeck_->hasKeyword("PERMY")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMY"));
if (deck_->hasKeyword("PERMY")) {
auto data = getAllSiDoubles_(deck_->getKeyword("PERMY"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMY", data);
}
if (newParserDeck_->hasKeyword("PERMZ")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMZ"));
if (deck_->hasKeyword("PERMZ")) {
auto data = getAllSiDoubles_(deck_->getKeyword("PERMZ"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMZ", data);
}
/* Create summary object (could not do it at construction time,
since it requires knowledge of the start time). */
summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, newParserDeck_));
summary_->addWells (newParserDeck_, uses_);
summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, deck_));
summary_->addWells (deck_, uses_);
}
void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
@@ -1112,7 +1114,7 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
// start writing to files
EclipseRestart rst(outputDir_, baseName_, timer, outputTimeStepIdx_);
rst.writeHeader (timer, outputTimeStepIdx_, uses_, newParserDeck_, reservoirState.pressure().size ());
rst.writeHeader (timer, outputTimeStepIdx_, uses_, deck_, reservoirState.pressure().size ());
EclipseSolution sol (rst);
// write out the pressure of the reference phase (whatever
@@ -1192,29 +1194,29 @@ void EclipseWriter::writeTimeStep(
EclipseWriter::EclipseWriter (
const ParameterGroup& params,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
int number_of_cells, const int* global_cell, const int* cart_dims,
int dimensions)
: newParserDeck_ (newParserDeck)
: deck_ (deck)
, number_of_cells_(number_of_cells)
, dimensions_(dimensions)
, cart_dims_(cart_dims)
, global_cell_(global_cell)
, uses_ (phaseUsageFromDeck (newParserDeck_)) {
, uses_ (phaseUsageFromDeck (deck_)) {
init(params);
}
EclipseWriter::EclipseWriter (
const ParameterGroup& params,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
std::shared_ptr<const UnstructuredGrid> grid)
: newParserDeck_ (newParserDeck)
: deck_ (deck)
, grid_(grid)
, number_of_cells_(grid->number_of_cells)
, dimensions_(grid->dimensions)
, cart_dims_(grid->cartdims)
, global_cell_(grid->global_cell)
, uses_ (phaseUsageFromDeck (newParserDeck_)) {
, uses_ (phaseUsageFromDeck (deck_)) {
init(params);
}
void EclipseWriter::init(const ParameterGroup& params)
opm/core/io/eclipse/EclipseWriter.hpp
+3 -3
View File
@@ -59,7 +59,7 @@ public:
* binary files using ERT.
*/
EclipseWriter(const parameter::ParameterGroup& params,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
std::shared_ptr <const UnstructuredGrid> grid);
/*!
@@ -67,7 +67,7 @@ public:
* binary files using ERT.
*/
EclipseWriter(const parameter::ParameterGroup& params,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
int number_of_cells, const int* global_cell, const int* cart_dims,
int dimension);
@@ -103,7 +103,7 @@ public:
const WellState& wellState);
private:
Opm::DeckConstPtr newParserDeck_;
Opm::DeckConstPtr deck_;
std::shared_ptr <const UnstructuredGrid> grid_;
int number_of_cells_;
int dimensions_;
opm/core/io/eclipse/SpecialEclipseFields.hpp
-2520
View File
File diff suppressed because it is too large. Load diff
opm/core/props/BlackoilPropertiesFromDeck.cpp
+3 -21
View File
@@ -23,24 +23,15 @@
namespace Opm
{
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid,
bool init_rock)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
grid.dimensions, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
bool init_rock)
{
init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims,
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims,
grid.cell_centroids, grid.dimensions, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock)
@@ -49,15 +40,6 @@ namespace Opm
grid.dimensions, param, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock)
{
init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
grid.dimensions, param, init_rock);
}
BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()
{
}
opm/core/props/BlackoilPropertiesFromDeck.hpp
+14 -76
View File
@@ -25,7 +25,6 @@
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -47,15 +46,7 @@ namespace Opm
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(const EclipseGridParser &deck,
const UnstructuredGrid& grid, bool init_rock=true );
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid, bool init_rock=true );
/// Initialize from deck, grid and parameters.
@@ -69,62 +60,27 @@ namespace Opm
/// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2").
/// For both size parameters, a 0 or negative value indicates that no spline fitting is to
/// be done, and the input fluid data used directly for linear interpolation.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock=true);
/// Initialize from deck, grid and parameters.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
/// \param[in] param Parameters. Accepted parameters include:
/// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables.
/// sat_tab_size (200) number of uniform sample points for saturation tables.
/// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2").
/// For both size parameters, a 0 or negative value indicates that no spline fitting is to
/// be done, and the input fluid data used directly for linear interpolation.
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock=true);
template<class T>
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
template<class CentroidIterator>
BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
bool init_rock=true);
template<class T>
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
template<class CentroidIterator>
BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock=true);
template<class T>
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock=true);
template<class T>
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock=true);
@@ -250,38 +206,20 @@ namespace Opm
double* smax) const;
private:
template<class T>
void init(const EclipseGridParser& deck,
template<class CentroidIterator>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
bool init_rock);
template<class T>
void init(const EclipseGridParser& deck,
template<class CentroidIterator>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock);
opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
+34 -163
View File
@@ -1,13 +1,11 @@
namespace Opm
{
template<class T>
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
template<class CentroidIterator>
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
bool init_rock)
{
@@ -15,13 +13,12 @@ namespace Opm
init_rock);
}
template<class T>
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
template<class CentroidIterator>
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
const CentroidIterator& begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock)
@@ -30,74 +27,23 @@ namespace Opm
param, init_rock);
}
template<class T>
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock)
{
init(newParserDeck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
init_rock);
}
template<class T>
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock)
{
init(newParserDeck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
param, init_rock);
}
template<class T>
void BlackoilPropertiesFromDeck::init(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock)
template<class CentroidIterator>
inline void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
const CentroidIterator& begin_cell_centroids,
int dimension,
bool init_rock)
{
if (init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
pvt_.init(deck, 0);
pvt_.init(deck, /*numSamples=*/0);
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension, 200);
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
template<class T>
void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock)
{
if (init_rock){
rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
}
pvt_.init(newParserDeck, /*numSamples=*/0);
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids, dimension,
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
/*numSamples=*/0);
if (pvt_.numPhases() != satprops_->numPhases()) {
@@ -106,103 +52,28 @@ namespace Opm
}
}
template<class T>
void BlackoilPropertiesFromDeck::init(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock)
template<class CentroidIterator>
inline void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
const CentroidIterator& begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock)
{
if(init_rock){
rock_.init(deck, number_of_cells, global_cell, cart_dims);
}
const int pvt_samples = param.getDefault("pvt_tab_size", 0);
pvt_.init(deck, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 0);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (deck.hasField("ENDSCALE") && threephase_model != "simple") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
}
} else {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
}
}
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
template<class T>
void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock)
{
if(init_rock){
rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
}
const int pvt_samples = param.getDefault("pvt_tab_size", 200);
pvt_.init(newParserDeck, pvt_samples);
pvt_.init(deck, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "gwseg") {
if (deck->hasKeyword("ENDSCALE") && threephase_model != "gwseg") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only.");
}
if (sat_samples > 1) {
@@ -210,19 +81,19 @@ namespace Opm
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
@@ -232,19 +103,19 @@ namespace Opm
SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else {
OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
opm/core/props/IncompPropertiesFromDeck.cpp
+2 -14
View File
@@ -26,10 +26,10 @@
namespace Opm
{
IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck,
IncompPropertiesFromDeck::IncompPropertiesFromDeck(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid)
{
rock_.init(deck, grid);
rock_.init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims);
pvt_.init(deck);
satprops_.init(deck, grid, 200);
if (pvt_.numPhases() != satprops_.numPhases()) {
@@ -38,18 +38,6 @@ namespace Opm
}
}
IncompPropertiesFromDeck::IncompPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid)
{
rock_.init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims);
pvt_.init(newParserDeck);
satprops_.init(newParserDeck, grid, 200);
if (pvt_.numPhases() != satprops_.numPhases()) {
OPM_THROW(std::runtime_error, "IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
}
}
IncompPropertiesFromDeck::~IncompPropertiesFromDeck()
{
}
opm/core/props/IncompPropertiesFromDeck.hpp
+1 -9
View File
@@ -50,15 +50,7 @@ namespace Opm
/// \param grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
IncompPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \param grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
IncompPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
IncompPropertiesFromDeck(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid);
/// Destructor.
opm/core/props/phaseUsageFromDeck.hpp
+5 -43
View File
@@ -20,9 +20,8 @@
#ifndef OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED
#define OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -70,56 +69,19 @@ namespace Opm
/// Looks at presence of WATER, OIL and GAS keywords in deck
/// to determine active phases.
inline PhaseUsage phaseUsageFromDeck(const EclipseGridParser& deck)
inline PhaseUsage phaseUsageFromDeck(Opm::DeckConstPtr deck)
{
PhaseUsage pu;
std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0);
// Discover phase usage.
if (deck.hasField("WATER")) {
if (deck->hasKeyword("WATER")) {
pu.phase_used[BlackoilPhases::Aqua] = 1;
}
if (deck.hasField("OIL")) {
if (deck->hasKeyword("OIL")) {
pu.phase_used[BlackoilPhases::Liquid] = 1;
}
if (deck.hasField("GAS")) {
pu.phase_used[BlackoilPhases::Vapour] = 1;
}
pu.num_phases = 0;
for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) {
pu.phase_pos[i] = pu.num_phases;
pu.num_phases += pu.phase_used[i];
}
// Only 2 or 3 phase systems handled.
if (pu.num_phases < 2 || pu.num_phases > 3) {
OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases.");
}
// We need oil systems, since we do not support the keywords needed for
// water-gas systems.
if (!pu.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems.");
}
return pu;
}
/// Looks at presence of WATER, OIL and GAS keywords in deck
/// to determine active phases.
inline PhaseUsage phaseUsageFromDeck(Opm::DeckConstPtr newParserDeck)
{
PhaseUsage pu;
std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0);
// Discover phase usage.
if (newParserDeck->hasKeyword("WATER")) {
pu.phase_used[BlackoilPhases::Aqua] = 1;
}
if (newParserDeck->hasKeyword("OIL")) {
pu.phase_used[BlackoilPhases::Liquid] = 1;
}
if (newParserDeck->hasKeyword("GAS")) {
if (deck->hasKeyword("GAS")) {
pu.phase_used[BlackoilPhases::Vapour] = 1;
}
pu.num_phases = 0;
opm/core/props/pvt/BlackoilPvtProperties.cpp
+15 -88
View File
@@ -27,7 +27,6 @@
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
@@ -43,7 +42,7 @@ namespace Opm
{
}
void BlackoilPvtProperties::init(const EclipseGridParser& deck, const int samples)
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, int samples)
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
region_number_ = 0;
@@ -51,80 +50,8 @@ namespace Opm
phase_usage_ = phaseUsageFromDeck(deck);
// Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number_];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
if (phase_usage_.phase_used[Aqua]) {
densities_[phase_usage_.phase_pos[Aqua]] = d[ECL_water];
}
if (phase_usage_.phase_used[Vapour]) {
densities_[phase_usage_.phase_pos[Vapour]] = d[ECL_gas];
}
if (phase_usage_.phase_used[Liquid]) {
densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil];
}
} else {
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
}
// Set the properties.
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (deck.hasField("PVTW")) {
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(deck.getPVTW().pvtw_));
} else {
// Eclipse 100 default.
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
}
}
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (deck.hasField("PVDO")) {
if (samples > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(deck.getPVDO().pvdo_, samples));
} else {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_));
}
} else if (deck.hasField("PVTO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_));
} else if (deck.hasField("PVCDO")) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(deck.getPVCDO().pvcdo_));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
}
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
if (deck.hasField("PVDG")) {
if (samples > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(deck.getPVDG().pvdg_, samples));
} else {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_));
}
} else if (deck.hasField("PVTG")) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
}
// Must inform pvt property objects of phase structure.
for (int i = 0; i < phase_usage_.num_phases; ++i) {
props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos);
}
}
void BlackoilPvtProperties::init(Opm::DeckConstPtr newParserDeck, int samples)
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
region_number_ = 0;
phase_usage_ = phaseUsageFromDeck(newParserDeck);
// Surface densities. Accounting for different orders in eclipse and our code.
if (newParserDeck->hasKeyword("DENSITY")) {
Opm::DeckKeywordConstPtr densityKeyword = newParserDeck->getKeyword("DENSITY");
if (deck->hasKeyword("DENSITY")) {
Opm::DeckKeywordConstPtr densityKeyword = deck->getKeyword("DENSITY");
if (phase_usage_.phase_used[Liquid]) {
densities_[phase_usage_.phase_pos[Liquid]]
= densityKeyword->getRecord(region_number_)->getItem("OIL")->getSIDouble(0);
@@ -145,8 +72,8 @@ namespace Opm
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (newParserDeck->hasKeyword("PVTW")) {
Opm::PvtwTable pvtwTable(newParserDeck->getKeyword("PVTW"));
if (deck->hasKeyword("PVTW")) {
Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
} else {
@@ -156,19 +83,19 @@ namespace Opm
}
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (newParserDeck->hasKeyword("PVDO")) {
Opm::PvdoTable pvdoTable(newParserDeck->getKeyword("PVDO"), region_number_);
if (deck->hasKeyword("PVDO")) {
Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"), region_number_);
if (samples > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
} else {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
}
} else if (newParserDeck->hasKeyword("PVTO")) {
Opm::PvtoTable pvtoTable(newParserDeck->getKeyword("PVTO"), /*tableIdx=*/0);
} else if (deck->hasKeyword("PVTO")) {
Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
} else if (newParserDeck->hasKeyword("PVCDO")) {
Opm::PvcdoTable pvcdoTable(newParserDeck->getKeyword("PVCDO"));
} else if (deck->hasKeyword("PVCDO")) {
Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
} else {
@@ -177,15 +104,15 @@ namespace Opm
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
if (newParserDeck->hasKeyword("PVDG")) {
Opm::PvdgTable pvdgTable(newParserDeck->getKeyword("PVDG"), region_number_);
if (deck->hasKeyword("PVDG")) {
Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"), region_number_);
if (samples > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
} else {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
}
} else if (newParserDeck->hasKeyword("PVTG")) {
Opm::PvtgTable pvtgTable(newParserDeck->getKeyword("PVTG"), /*tableIdx=*/0);
} else if (deck->hasKeyword("PVTG")) {
Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
} else {
opm/core/props/pvt/BlackoilPvtProperties.hpp
-11
View File
@@ -22,7 +22,6 @@
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -47,16 +46,6 @@ namespace Opm
/// Default constructor.
BlackoilPvtProperties();
/// Initialize from deck.
/// \param deck An input deck.
/// \param samples If greater than zero, indicates the number of
/// uniform samples to be taken from monotone spline
/// curves interpolating the fluid data.
/// Otherwise, interpolate linearly in the original
/// data without fitting a spline.
void init(const EclipseGridParser& deck, const int samples);
/// Initialize from deck.
/// \param deck An input deck from the opm-parser module.
void init(Opm::DeckConstPtr deck, int samples);
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
+7 -62
View File
@@ -21,7 +21,6 @@
#include "config.h"
#include <opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
@@ -34,8 +33,7 @@ namespace Opm
{
}
void PvtPropertiesIncompFromDeck::init(const EclipseGridParser& deck)
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck )
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
int region_number = 0;
@@ -48,61 +46,8 @@ namespace Opm
}
// Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasField("DENSITY")) {
const std::vector<double>& d = deck.getDENSITY().densities_[region_number];
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = d[ECL_water];
surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = d[ECL_oil];
} else {
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
}
// Make reservoir densities the same as surface densities initially.
// We will modify them with formation volume factors if found.
reservoir_density_ = surface_density_;
// Water viscosity.
if (deck.hasField("PVTW")) {
const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
OPM_MESSAGE("Compressibility effects in PVTW are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] /= pvtw[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
} else {
// Eclipse 100 default.
// viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = 0.5*Opm::prefix::centi*Opm::unit::Poise;
OPM_THROW(std::runtime_error, "Input is missing PVTW\n");
}
// Oil viscosity.
if (deck.hasField("PVCDO")) {
const std::vector<double>& pvcdo = deck.getPVCDO().pvcdo_[region_number];
if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
OPM_MESSAGE("Compressibility effects in PVCDO are ignored.");
}
reservoir_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] /= pvcdo[1];
viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
} else {
OPM_THROW(std::runtime_error, "Input is missing PVCDO\n");
}
}
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr newParserDeck )
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
int region_number = 0;
PhaseUsage phase_usage = phaseUsageFromDeck(newParserDeck);
if (phase_usage.phase_used[PhaseUsage::Vapour] ||
!phase_usage.phase_used[PhaseUsage::Aqua] ||
!phase_usage.phase_used[PhaseUsage::Liquid]) {
OPM_THROW(std::runtime_error, "PvtPropertiesIncompFromDeck::init() -- must have gas and oil phases (only) in deck input.\n");
}
// Surface densities. Accounting for different orders in eclipse and our code.
if (newParserDeck->hasKeyword("DENSITY")) {
Opm::DeckRecordConstPtr densityRecord = newParserDeck->getKeyword("DENSITY")->getRecord(region_number);
if (deck->hasKeyword("DENSITY")) {
Opm::DeckRecordConstPtr densityRecord = deck->getKeyword("DENSITY")->getRecord(region_number);
surface_density_[phase_usage.phase_pos[PhaseUsage::Aqua]] = densityRecord->getItem("OIL")->getSIDouble(0);
surface_density_[phase_usage.phase_pos[PhaseUsage::Liquid]] = densityRecord->getItem("WATER")->getSIDouble(0);
} else {
@@ -114,8 +59,8 @@ namespace Opm
reservoir_density_ = surface_density_;
// Water viscosity.
if (newParserDeck->hasKeyword("PVTW")) {
Opm::DeckRecordConstPtr pvtwRecord = newParserDeck->getKeyword("PVTW")->getRecord(region_number);
if (deck->hasKeyword("PVTW")) {
Opm::DeckRecordConstPtr pvtwRecord = deck->getKeyword("PVTW")->getRecord(region_number);
if (pvtwRecord->getItem("WATER_COMPRESSIBILITY")->getSIDouble(0) != 0.0 ||
pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0) != 0.0) {
OPM_MESSAGE("Compressibility effects in PVTW are ignored.");
@@ -129,8 +74,8 @@ namespace Opm
}
// Oil viscosity.
if (newParserDeck->hasKeyword("PVCDO")) {
Opm::DeckRecordConstPtr pvcdoRecord = newParserDeck->getKeyword("PVCDO")->getRecord(region_number);
if (deck->hasKeyword("PVCDO")) {
Opm::DeckRecordConstPtr pvcdoRecord = deck->getKeyword("PVCDO")->getRecord(region_number);
if (pvcdoRecord->getItem("OIL_COMPRESSIBILITY")->getSIDouble(0) != 0.0 ||
pvcdoRecord->getItem("OIL_VISCOSIBILITY")->getSIDouble(0) != 0.0) {
opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
+1 -5
View File
@@ -20,7 +20,6 @@
#ifndef OPM_PVTPROPERTIESINCOMPFROMDECK_HEADER_INCLUDED
#define OPM_PVTPROPERTIESINCOMPFROMDECK_HEADER_INCLUDED
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <array>
@@ -39,10 +38,7 @@ namespace Opm
PvtPropertiesIncompFromDeck();
/// Initialize from deck.
void init(const EclipseGridParser& deck);
/// Initialize from deck.
void init(Opm::DeckConstPtr newParserDeck);
void init(Opm::DeckConstPtr deck);
/// Number of active phases.
int numPhases() const;
opm/core/props/pvt/SinglePvtConstCompr.hpp
+8 -23
View File
@@ -43,21 +43,6 @@ namespace Opm
class SinglePvtConstCompr : public SinglePvtInterface
{
public:
typedef std::vector<std::vector<double> > table_t;
SinglePvtConstCompr(const table_t& pvtw)
{
const int region_number = 0;
if (pvtw.size() != 1) {
OPM_THROW(std::runtime_error, "More than one PVD-region");
}
ref_press_ = pvtw[region_number][0];
ref_B_ = pvtw[region_number][1];
comp_ = pvtw[region_number][2];
viscosity_ = pvtw[region_number][3];
visc_comp_ = pvtw[region_number][4];
}
SinglePvtConstCompr(const Opm::PvtwTable &pvtwTable)
{
if (pvtwTable.numRows() != 1)
@@ -71,17 +56,17 @@ namespace Opm
visc_comp_ = pvtwTable.getViscosibilityColumn()[0];
}
SinglePvtConstCompr(const Opm::PvcdoTable &pvdcoTable)
SinglePvtConstCompr(const Opm::PvcdoTable &pvcdoTable)
{
if (pvdcoTable.numRows() != 1)
if (pvcdoTable.numRows() != 1)
OPM_THROW(std::runtime_error,
"The table specified by the PVDCO keyword is required"
"The table specified by the PVCDO keyword is required"
"to exhibit exactly one row.");
ref_press_ = pvdcoTable.getPressureColumn()[0];
ref_B_ = pvdcoTable.getFormationFactorColumn()[0];
comp_ = pvdcoTable.getCompressibilityColumn()[0];
viscosity_ = pvdcoTable.getViscosityColumn()[0];
visc_comp_ = pvdcoTable.getViscosibilityColumn()[0];
ref_press_ = pvcdoTable.getPressureColumn()[0];
ref_B_ = pvcdoTable.getFormationFactorColumn()[0];
comp_ = pvcdoTable.getCompressibilityColumn()[0];
viscosity_ = pvcdoTable.getViscosityColumn()[0];
visc_comp_ = pvcdoTable.getViscosibilityColumn()[0];
}
SinglePvtConstCompr(double visc)
opm/core/props/pvt/SinglePvtDead.cpp
-29
View File
@@ -33,35 +33,6 @@ namespace Opm
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
SinglePvtDead::SinglePvtDead(const table_t& pvd_table)
{
const int region_number = 0;
if (pvd_table.size() != 1) {
OPM_THROW(std::runtime_error, "More than one PVT-region");
}
// Copy data
const int sz = pvd_table[region_number][0].size();
std::vector<double> press(sz);
std::vector<double> b(sz);
std::vector<double> visc(sz);
for (int i = 0; i < sz; ++i) {
press[i] = pvd_table[region_number][0][i];
b[i] = 1.0 / pvd_table[region_number][1][i];
visc[i] = pvd_table[region_number][2][i];
}
b_ = NonuniformTableLinear<double>(press, b);
viscosity_ = NonuniformTableLinear<double>(press, visc);
// Dumping the created tables.
// static int count = 0;
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
// os.precision(15);
// os << "1/B\n\n" << one_over_B_
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
}
/// Constructor
SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable)
{
opm/core/props/pvt/SinglePvtDead.hpp
-2
View File
@@ -42,8 +42,6 @@ namespace Opm
class SinglePvtDead : public SinglePvtInterface
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtDead(const table_t& pvd_table);
SinglePvtDead(const Opm::PvdoTable& pvdoTable);
SinglePvtDead(const Opm::PvdgTable& pvdgTable);
virtual ~SinglePvtDead();
opm/core/props/pvt/SinglePvtDeadSpline.cpp
-22
View File
@@ -36,28 +36,6 @@ namespace Opm
// Member functions
//-------------------------------------------------------------------------
/// Constructor
SinglePvtDeadSpline::SinglePvtDeadSpline(const table_t& pvd_table, const int samples)
{
const int region_number = 0;
if (pvd_table.size() != 1) {
OPM_THROW(std::runtime_error, "More than one PVT-region");
}
// Copy data
const int sz = pvd_table[region_number][0].size();
std::vector<double> press(sz);
std::vector<double> B_inv(sz);
std::vector<double> visc(sz);
for (int i = 0; i < sz; ++i) {
press[i] = pvd_table[region_number][0][i];
B_inv[i] = 1.0 / pvd_table[region_number][1][i];
visc[i] = pvd_table[region_number][2][i];
}
buildUniformMonotoneTable(press, B_inv, samples, b_);
buildUniformMonotoneTable(press, visc, samples, viscosity_);
}
SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples)
{
int numRows = pvdoTable.numRows();
opm/core/props/pvt/SinglePvtDeadSpline.hpp
-3
View File
@@ -41,9 +41,6 @@ namespace Opm
class SinglePvtDeadSpline : public SinglePvtInterface
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtDeadSpline(const table_t& pvd_table, const int samples);
SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples);
SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples);
virtual ~SinglePvtDeadSpline();
opm/core/props/pvt/SinglePvtLiveGas.cpp
-37
View File
@@ -46,43 +46,6 @@ namespace Opm
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
SinglePvtLiveGas::SinglePvtLiveGas(const table_t& pvtg)
{
// GAS, PVTG
const int region_number = 0;
if (pvtg.size() != 1) {
OPM_THROW(std::runtime_error, "More than one PVD-region");
}
saturated_gas_table_.resize(4);
const int sz = pvtg[region_number].size();
for (int k=0; k<4; ++k) {
saturated_gas_table_[k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_gas_table_[0][i] = pvtg[region_number][i][0]; // p
saturated_gas_table_[1][i] = 1.0/pvtg[region_number][i][2]; // 1/Bg
saturated_gas_table_[2][i] = pvtg[region_number][i][3]; // mu_g
saturated_gas_table_[3][i] = pvtg[region_number][i][1]; // Rv
}
undersat_gas_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
undersat_gas_tables_[i].resize(3);
int tsize = (pvtg[region_number][i].size() - 1)/3;
undersat_gas_tables_[i][0].resize(tsize);
undersat_gas_tables_[i][1].resize(tsize);
undersat_gas_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) {
undersat_gas_tables_[i][0][j] = pvtg[region_number][i][++k]; // Rv
undersat_gas_tables_[i][1][j] = 1.0/pvtg[region_number][i][++k]; // 1/Bg
undersat_gas_tables_[i][2][j] = pvtg[region_number][i][++k]; // mu_g
}
}
}
SinglePvtLiveGas::SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable)
{
// GAS, PVTG
opm/core/props/pvt/SinglePvtLiveGas.hpp
-3
View File
@@ -38,9 +38,6 @@ namespace Opm
class SinglePvtLiveGas : public SinglePvtInterface
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtLiveGas(const table_t& pvto);
SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable);
virtual ~SinglePvtLiveGas();
opm/core/props/pvt/SinglePvtLiveOil.cpp
-68
View File
@@ -35,74 +35,6 @@ namespace Opm
// Member functions
//-------------------------------------------------------------------------
/// Constructor
SinglePvtLiveOil::SinglePvtLiveOil(const table_t& pvto)
{
// OIL, PVTO
const int region_number = 0;
if (pvto.size() != 1) {
OPM_THROW(std::runtime_error, "More than one PVD-region");
}
saturated_oil_table_.resize(4);
const int sz = pvto[region_number].size();
for (int k=0; k<4; ++k) {
saturated_oil_table_[k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_oil_table_[0][i] = pvto[region_number][i][1]; // p
saturated_oil_table_[1][i] = 1.0/pvto[region_number][i][2]; // 1/Bo
saturated_oil_table_[2][i] = pvto[region_number][i][3]; // mu_o
saturated_oil_table_[3][i] = pvto[region_number][i][0]; // Rs
}
undersat_oil_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
undersat_oil_tables_[i].resize(3);
int tsize = (pvto[region_number][i].size() - 1)/3;
undersat_oil_tables_[i][0].resize(tsize);
undersat_oil_tables_[i][1].resize(tsize);
undersat_oil_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) {
undersat_oil_tables_[i][0][j] = pvto[region_number][i][++k]; // p
undersat_oil_tables_[i][1][j] = 1.0/pvto[region_number][i][++k]; // 1/Bo
undersat_oil_tables_[i][2][j] = pvto[region_number][i][++k]; // mu_o
}
}
// Complete undersaturated tables by extrapolating from existing data
// as is done in Eclipse and Mrst
int iNext = -1;
for (int i=0; i<sz; ++i) {
// Skip records already containing undersaturated data
if (undersat_oil_tables_[i][0].size() > 1) {
continue;
}
// Look ahead for next record containing undersaturated data
if (iNext < i) {
iNext = i+1;
while (iNext<sz && undersat_oil_tables_[iNext][0].size() < 2) {
++iNext;
}
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
}
// Add undersaturated data to current record while maintaining compressibility and viscosibility
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
for (sz_t j=1; j<undersat_oil_tables_[iNext][0].size(); ++j) {
double diffPressure = undersat_oil_tables_[iNext][0][j]-undersat_oil_tables_[iNext][0][j-1];
double pressure = undersat_oil_tables_[i][0].back()+diffPressure;
undersat_oil_tables_[i][0].push_back(pressure);
double compr = (1.0/undersat_oil_tables_[iNext][1][j]-1.0/undersat_oil_tables_[iNext][1][j-1])
/ (0.5*(1.0/undersat_oil_tables_[iNext][1][j]+1.0/undersat_oil_tables_[iNext][1][j-1]));
double B = (1.0/undersat_oil_tables_[i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
undersat_oil_tables_[i][1].push_back(1.0/B);
double visc = (undersat_oil_tables_[iNext][2][j]-undersat_oil_tables_[iNext][2][j-1])
/ (0.5*(undersat_oil_tables_[iNext][2][j]+undersat_oil_tables_[iNext][2][j-1]));
double mu = (undersat_oil_tables_[i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
undersat_oil_tables_[i][2].push_back(mu);
}
}
}
SinglePvtLiveOil::SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable)
{
const auto saturatedPvto = pvtoTable.getOuterTable();
opm/core/props/pvt/SinglePvtLiveOil.hpp
-3
View File
@@ -38,9 +38,6 @@ namespace Opm
class SinglePvtLiveOil : public SinglePvtInterface
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtLiveOil(const table_t& pvto);
SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable);
virtual ~SinglePvtLiveOil();
opm/core/props/rock/RockCompressibility.cpp
+6 -34
View File
@@ -19,7 +19,6 @@
#include "config.h"
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@@ -41,46 +40,19 @@ namespace Opm
rock_comp_ = param.getDefault("rock_compressibility", 0.0)/unit::barsa;
}
RockCompressibility::RockCompressibility(const EclipseGridParser& deck)
RockCompressibility::RockCompressibility(Opm::DeckConstPtr deck)
: pref_(0.0),
rock_comp_(0.0)
{
if (deck.hasField("ROCKTAB")) {
const table_t& rt = deck.getROCKTAB().rocktab_;
if (rt.size() != 1) {
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB.");
}
const int n = rt[0][0].size();
p_.resize(n);
poromult_.resize(n);
transmult_.resize(n);
for (int i = 0; i < n; ++i) {
p_[i] = rt[0][0][i];
poromult_[i] = rt[0][1][i];
transmult_[i] = rt[0][2][i];
}
} else if (deck.hasField("ROCK")) {
const ROCK& r = deck.getROCK();
pref_ = r.rock_compressibilities_[0][0];
rock_comp_ = r.rock_compressibilities_[0][1];
} else {
std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl;
}
}
RockCompressibility::RockCompressibility(Opm::DeckConstPtr newParserDeck)
: pref_(0.0),
rock_comp_(0.0)
{
if (newParserDeck->hasKeyword("ROCKTAB")) {
Opm::DeckKeywordConstPtr rtKeyword = newParserDeck->getKeyword("ROCKTAB");
if (deck->hasKeyword("ROCKTAB")) {
Opm::DeckKeywordConstPtr rtKeyword = deck->getKeyword("ROCKTAB");
if (rtKeyword->size() != 1)
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB.");
// the number of colums of the "ROCKTAB" keyword
// depends on the presence of the "RKTRMDIR"
// keyword. Messy stuff...
bool isDirectional = newParserDeck->hasKeyword("RKTRMDIR");
bool isDirectional = deck->hasKeyword("RKTRMDIR");
if (isDirectional)
{
// well, okay. we don't support non-isotropic
@@ -94,8 +66,8 @@ namespace Opm
p_ = rocktabTable.getPressureColumn();
poromult_ = rocktabTable.getPoreVolumeMultiplierColumn();
transmult_ = rocktabTable.getTransmissibilityMultiplierColumn();
} else if (newParserDeck->hasKeyword("ROCK")) {
Opm::RockTable rockTable(newParserDeck->getKeyword("ROCK"));
} else if (deck->hasKeyword("ROCK")) {
Opm::RockTable rockTable(deck->getKeyword("ROCK"));
if (rockTable.numRows() != 1)
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK.");
opm/core/props/rock/RockCompressibility.hpp
+1 -6
View File
@@ -27,7 +27,6 @@
namespace Opm
{
class EclipseGridParser;
namespace parameter { class ParameterGroup; }
class RockCompressibility
@@ -35,11 +34,7 @@ namespace Opm
public:
/// Construct from input deck.
/// Looks for the keywords ROCK and ROCKTAB.
RockCompressibility(const EclipseGridParser& deck);
/// Construct from input deck.
/// Looks for the keywords ROCK and ROCKTAB.
RockCompressibility(Opm::DeckConstPtr newParserDeck);
RockCompressibility(Opm::DeckConstPtr deck);
/// Construct from parameters.
/// Accepts the following parameters (with defaults).
opm/core/props/rock/RockFromDeck.cpp
+41 -309
View File
@@ -21,6 +21,7 @@
#include "config.h"
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -34,13 +35,9 @@ namespace Opm
{
enum PermeabilityKind { ScalarPerm, DiagonalPerm, TensorPerm, None, Invalid };
PermeabilityKind classifyPermeability(const EclipseGridParser& parser);
void setScalarPermIfNeeded(std::array<int,9>& kmap,
int i, int j, int k);
PermeabilityKind fillTensor(const EclipseGridParser& parser,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap);
PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
PermeabilityKind fillTensor(Opm::DeckConstPtr deck,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap);
@@ -56,65 +53,23 @@ namespace Opm
{
}
/// Initialize from deck and cell mapping.
/// \param deck Deck input parser
/// \param grid grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
void RockFromDeck::init(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims);
}
/// Initialize from deck and cell mapping.
/// \param deck Deck input parser
/// \param number_of_cells The number of cells in the grid.
/// \param global_cell The mapping fom local to global cell indices.
/// global_cell[i] is the corresponding global index of i.
/// \param cart_dims The size of the underlying cartesian grid.
void RockFromDeck::init(const EclipseGridParser& deck,
void RockFromDeck::init(Opm::DeckConstPtr deck,
int number_of_cells, const int* global_cell,
const int* cart_dims)
{
assignPorosity(deck, number_of_cells, global_cell);
permfield_valid_.assign(number_of_cells, false);
const double perm_threshold = 0.0; // Maybe turn into parameter?
assignPermeability(deck, number_of_cells, global_cell, cart_dims, perm_threshold);
}
void RockFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell,
const int* cart_dims)
{
assignPorosity(newParserDeck, number_of_cells, global_cell);
permfield_valid_.assign(number_of_cells, false);
const double perm_threshold = 0.0; // Maybe turn into parameter?
assignPermeability(newParserDeck, number_of_cells, global_cell, cart_dims,
assignPermeability(deck, number_of_cells, global_cell, cart_dims,
perm_threshold);
}
void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
void RockFromDeck::assignPorosity(Opm::DeckConstPtr deck,
int number_of_cells, const int* global_cell)
{
porosity_.assign(number_of_cells, 1.0);
if (parser.hasField("PORO")) {
const std::vector<double>& poro = parser.getFloatingPointValue("PORO");
for (int c = 0; c < int(porosity_.size()); ++c) {
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
porosity_[c] = poro[deck_pos];
}
}
}
void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell)
{
porosity_.assign(number_of_cells, 1.0);
if (newParserDeck->hasKeyword("PORO")) {
const std::vector<double>& poro = newParserDeck->getKeyword("PORO")->getSIDoubleData();
if (deck->hasKeyword("PORO")) {
const std::vector<double>& poro = deck->getKeyword("PORO")->getSIDoubleData();
for (int c = 0; c < int(porosity_.size()); ++c) {
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
assert(0 <= c && c < (int) porosity_.size());
@@ -124,62 +79,7 @@ namespace Opm
}
}
void RockFromDeck::assignPermeability(const EclipseGridParser& parser,
int number_of_cells,
const int* global_cell,
const int* cartdims,
double perm_threshold)
{
const int dim = 3;
const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2];
assert(num_global_cells > 0);
permeability_.assign(dim * dim * number_of_cells, 0.0);
std::vector<const std::vector<double>*> tensor;
tensor.reserve(10);
const std::vector<double> zero(num_global_cells, 0.0);
tensor.push_back(&zero);
std::array<int,9> kmap;
PermeabilityKind pkind = fillTensor(parser, tensor, kmap);
if (pkind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid permeability field.");
}
// Assign permeability values only if such values are
// given in the input deck represented by 'parser'. In
// other words: Don't set any (arbitrary) default values.
// It is infinitely better to experience a reproducible
// crash than subtle errors resulting from a (poorly
// chosen) default value...
//
if (tensor.size() > 1) {
int off = 0;
for (int c = 0; c < number_of_cells; ++c, off += dim*dim) {
// SharedPermTensor K(dim, dim, &permeability_[off]);
int kix = 0;
const int glob = (global_cell == NULL) ? c : global_cell[c];
for (int i = 0; i < dim; ++i) {
for (int j = 0; j < dim; ++j, ++kix) {
// K(i,j) = (*tensor[kmap[kix]])[glob];
permeability_[off + kix] = (*tensor[kmap[kix]])[glob];
}
// K(i,i) = std::max(K(i,i), perm_threshold);
permeability_[off + 3*i + i] = std::max(permeability_[off + 3*i + i], perm_threshold);
}
permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
}
}
}
void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck,
void RockFromDeck::assignPermeability(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cartdims,
@@ -200,13 +100,13 @@ namespace Opm
tensor.push_back(&zero);
std::array<int,9> kmap;
PermeabilityKind pkind = fillTensor(newParserDeck, tensor, kmap);
PermeabilityKind pkind = fillTensor(deck, tensor, kmap);
if (pkind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid permeability field.");
}
// Assign permeability values only if such values are
// given in the input deck represented by 'newParserDeck'. In
// given in the input deck represented by 'deck'. In
// other words: Don't set any (arbitrary) default values.
// It is infinitely better to experience a reproducible
// crash than subtle errors resulting from a (poorly
@@ -236,7 +136,6 @@ namespace Opm
}
namespace {
/// @brief
/// Classify and verify a given permeability specification
/// from a structural point of view. In particular, we
@@ -244,7 +143,7 @@ namespace Opm
/// components such as @f$k_{xy}@f$ unless the
/// corresponding diagonal components are known as well.
///
/// @param parser [in]
/// @param deck [in]
/// An Eclipse data parser capable of answering which
/// permeability components are present in a given input
/// deck.
@@ -256,85 +155,19 @@ namespace Opm
/// TensorPerm at least one cross-component given.
/// None no components given.
/// Invalid invalid set of components given.
PermeabilityKind classifyPermeability(const EclipseGridParser& parser)
PermeabilityKind classifyPermeability(Opm::DeckConstPtr deck)
{
const bool xx = parser.hasField("PERMX" );
const bool xy = parser.hasField("PERMXY");
const bool xz = parser.hasField("PERMXZ");
const bool xx = deck->hasKeyword("PERMX" );
const bool xy = deck->hasKeyword("PERMXY");
const bool xz = deck->hasKeyword("PERMXZ");
const bool yx = parser.hasField("PERMYX");
const bool yy = parser.hasField("PERMY" );
const bool yz = parser.hasField("PERMYZ");
const bool yx = deck->hasKeyword("PERMYX");
const bool yy = deck->hasKeyword("PERMY" );
const bool yz = deck->hasKeyword("PERMYZ");
const bool zx = parser.hasField("PERMZX");
const bool zy = parser.hasField("PERMZY");
const bool zz = parser.hasField("PERMZ" );
int num_cross_comp = xy + xz + yx + yz + zx + zy;
int num_comp = xx + yy + zz + num_cross_comp;
PermeabilityKind retval = None;
if (num_cross_comp > 0) {
retval = TensorPerm;
} else {
if (num_comp == 1) {
retval = ScalarPerm;
} else if (num_comp >= 2) {
retval = DiagonalPerm;
}
}
bool ok = true;
if (num_comp > 0) {
// At least one tensor component specified on input.
// Verify that any remaining components are OK from a
// structural point of view. In particular, there
// must not be any cross-components (e.g., k_{xy})
// unless the corresponding diagonal component (e.g.,
// k_{xx}) is present as well...
//
ok = xx || !(xy || xz || yx || zx) ;
ok = ok && (yy || !(yx || yz || xy || zy));
ok = ok && (zz || !(zx || zy || xz || yz));
}
if (!ok) {
retval = Invalid;
}
return retval;
}
/// @brief
/// Classify and verify a given permeability specification
/// from a structural point of view. In particular, we
/// verify that there are no off-diagonal permeability
/// components such as @f$k_{xy}@f$ unless the
/// corresponding diagonal components are known as well.
///
/// @param newParserDeck [in]
/// An Eclipse data parser capable of answering which
/// permeability components are present in a given input
/// deck.
///
/// @return
/// An enum value with the following possible values:
/// ScalarPerm only one component was given.
/// DiagonalPerm more than one component given.
/// TensorPerm at least one cross-component given.
/// None no components given.
/// Invalid invalid set of components given.
PermeabilityKind classifyPermeability(Opm::DeckConstPtr newParserDeck)
{
const bool xx = newParserDeck->hasKeyword("PERMX" );
const bool xy = newParserDeck->hasKeyword("PERMXY");
const bool xz = newParserDeck->hasKeyword("PERMXZ");
const bool yx = newParserDeck->hasKeyword("PERMYX");
const bool yy = newParserDeck->hasKeyword("PERMY" );
const bool yz = newParserDeck->hasKeyword("PERMYZ");
const bool zx = newParserDeck->hasKeyword("PERMZX");
const bool zy = newParserDeck->hasKeyword("PERMZY");
const bool zz = newParserDeck->hasKeyword("PERMZ" );
const bool zx = deck->hasKeyword("PERMZX");
const bool zy = deck->hasKeyword("PERMZY");
const bool zz = deck->hasKeyword("PERMZ" );
int num_cross_comp = xy + xz + yx + yz + zx + zy;
int num_comp = xx + yy + zz + num_cross_comp;
@@ -397,7 +230,6 @@ namespace Opm
if (kmap[k] == 0) { kmap[k] = kmap[i]; }
}
/// @brief
/// Extract pointers to appropriate tensor components from
/// input deck. The permeability tensor is, generally,
@@ -430,11 +262,11 @@ namespace Opm
///
/// @param [out] tensor
/// @param [out] kmap
PermeabilityKind fillTensor(const EclipseGridParser& parser,
PermeabilityKind fillTensor(Opm::DeckConstPtr deck,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap)
{
PermeabilityKind kind = classifyPermeability(parser);
PermeabilityKind kind = classifyPermeability(deck);
if (kind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid set of permeability fields given.");
}
@@ -447,151 +279,51 @@ namespace Opm
// -----------------------------------------------------------
// 1st row: [kxx, kxy, kxz]
if (parser.hasField("PERMX" )) {
if (deck->hasKeyword("PERMX" )) {
kmap[xx] = tensor.size();
tensor.push_back(&parser.getFloatingPointValue("PERMX" ));
tensor.push_back(&deck->getKeyword("PERMX")->getSIDoubleData());
setScalarPermIfNeeded(kmap, xx, yy, zz);
}
if (parser.hasField("PERMXY")) {
if (deck->hasKeyword("PERMXY")) {
kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMXY"));
tensor.push_back(&deck->getKeyword("PERMXY")->getSIDoubleData());
}
if (parser.hasField("PERMXZ")) {
if (deck->hasKeyword("PERMXZ")) {
kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMXZ"));
tensor.push_back(&deck->getKeyword("PERMXZ")->getSIDoubleData());
}
// -----------------------------------------------------------
// 2nd row: [kyx, kyy, kyz]
if (parser.hasField("PERMYX")) {
if (deck->hasKeyword("PERMYX")) {
kmap[yx] = kmap[xy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMYX"));
tensor.push_back(&deck->getKeyword("PERMYX")->getSIDoubleData());
}
if (parser.hasField("PERMY" )) {
if (deck->hasKeyword("PERMY" )) {
kmap[yy] = tensor.size();
tensor.push_back(&parser.getFloatingPointValue("PERMY" ));
tensor.push_back(&deck->getKeyword("PERMY")->getSIDoubleData());
setScalarPermIfNeeded(kmap, yy, zz, xx);
}
if (parser.hasField("PERMYZ")) {
if (deck->hasKeyword("PERMYZ")) {
kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMYZ"));
tensor.push_back(&deck->getKeyword("PERMYZ")->getSIDoubleData());
}
// -----------------------------------------------------------
// 3rd row: [kzx, kzy, kzz]
if (parser.hasField("PERMZX")) {
if (deck->hasKeyword("PERMZX")) {
kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMZX"));
tensor.push_back(&deck->getKeyword("PERMZX")->getSIDoubleData());
}
if (parser.hasField("PERMZY")) {
if (deck->hasKeyword("PERMZY")) {
kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&parser.getFloatingPointValue("PERMZY"));
tensor.push_back(&deck->getKeyword("PERMZY")->getSIDoubleData());
}
if (parser.hasField("PERMZ" )) {
if (deck->hasKeyword("PERMZ" )) {
kmap[zz] = tensor.size();
tensor.push_back(&parser.getFloatingPointValue("PERMZ" ));
setScalarPermIfNeeded(kmap, zz, xx, yy);
}
return kind;
}
/// @brief
/// Extract pointers to appropriate tensor components from
/// input deck. The permeability tensor is, generally,
/// @code
/// [ kxx kxy kxz ]
/// K = [ kyx kyy kyz ]
/// [ kzx kzy kzz ]
/// @endcode
/// We store these values in a linear array using natural
/// ordering with the column index cycling the most rapidly.
/// In particular we use the representation
/// @code
/// [ 0 1 2 3 4 5 6 7 8 ]
/// K = [ kxx, kxy, kxz, kyx, kyy, kyz, kzx, kzy, kzz ]
/// @endcode
/// Moreover, we explicitly enforce symmetric tensors by
/// assigning
/// @code
/// 3 1 6 2 7 5
/// kyx = kxy, kzx = kxz, kzy = kyz
/// @endcode
/// However, we make no attempt at enforcing positive
/// definite tensors.
///
/// @param [in] parser
/// An Eclipse data parser capable of answering which
/// permeability components are present in a given input
/// deck as well as retrieving the numerical value of each
/// permeability component in each grid cell.
///
/// @param [out] tensor
/// @param [out] kmap
PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap)
{
PermeabilityKind kind = classifyPermeability(newParserDeck);
if (kind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid set of permeability fields given.");
}
assert(tensor.size() == 1);
for (int i = 0; i < 9; ++i) { kmap[i] = 0; }
enum { xx, xy, xz, // 0, 1, 2
yx, yy, yz, // 3, 4, 5
zx, zy, zz }; // 6, 7, 8
// -----------------------------------------------------------
// 1st row: [kxx, kxy, kxz]
if (newParserDeck->hasKeyword("PERMX" )) {
kmap[xx] = tensor.size();
tensor.push_back(&newParserDeck->getKeyword("PERMX")->getSIDoubleData());
setScalarPermIfNeeded(kmap, xx, yy, zz);
}
if (newParserDeck->hasKeyword("PERMXY")) {
kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMXY")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMXZ")) {
kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMXZ")->getSIDoubleData());
}
// -----------------------------------------------------------
// 2nd row: [kyx, kyy, kyz]
if (newParserDeck->hasKeyword("PERMYX")) {
kmap[yx] = kmap[xy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMYX")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMY" )) {
kmap[yy] = tensor.size();
tensor.push_back(&newParserDeck->getKeyword("PERMY")->getSIDoubleData());
setScalarPermIfNeeded(kmap, yy, zz, xx);
}
if (newParserDeck->hasKeyword("PERMYZ")) {
kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMYZ")->getSIDoubleData());
}
// -----------------------------------------------------------
// 3rd row: [kzx, kzy, kzz]
if (newParserDeck->hasKeyword("PERMZX")) {
kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMZX")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMZY")) {
kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMZY")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMZ" )) {
kmap[zz] = tensor.size();
tensor.push_back(&newParserDeck->getKeyword("PERMZ")->getSIDoubleData());
tensor.push_back(&deck->getKeyword("PERMZ")->getSIDoubleData());
setScalarPermIfNeeded(kmap, zz, xx, yy);
}
opm/core/props/rock/RockFromDeck.hpp
+4 -32
View File
@@ -20,9 +20,6 @@
#ifndef OPM_ROCKFROMDECK_HEADER_INCLUDED
#define OPM_ROCKFROMDECK_HEADER_INCLUDED
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <vector>
@@ -38,30 +35,13 @@ namespace Opm
/// Default constructor.
RockFromDeck();
/// Initialize from deck and grid.
/// \param deck Deck input parser
/// \param grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
void init(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
/// Initialize from deck and cell mapping.
/// \param deck Deck input parser
/// \param deck Deck produced by the opm-parser code
/// \param number_of_cells The number of cells in the grid.
/// \param global_cell The mapping fom local to global cell indices.
/// global_cell[i] is the corresponding global index of i.
/// \param cart_dims The size of the underlying cartesian grid.
void init(const EclipseGridParser& deck,
int number_of_cells, const int* global_cell,
const int* cart_dims);
/// Initialize from deck and cell mapping.
/// \param newParserDeck Deck produced by the opm-parser code
/// \param number_of_cells The number of cells in the grid.
/// \param global_cell The mapping fom local to global cell indices.
/// global_cell[i] is the corresponding global index of i.
/// \param cart_dims The size of the underlying cartesian grid.
void init(Opm::DeckConstPtr newParserDeck,
void init(Opm::DeckConstPtr deck,
int number_of_cells, const int* global_cell,
const int* cart_dims);
@@ -92,18 +72,10 @@ namespace Opm
}
private:
void assignPorosity(const EclipseGridParser& parser,
void assignPorosity(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell);
void assignPorosity(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell);
void assignPermeability(const EclipseGridParser& parser,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
const double perm_threshold);
void assignPermeability(Opm::DeckConstPtr newParserDeck,
void assignPermeability(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
opm/core/props/satfunc/SatFuncBase.hpp
+4 -141
View File
@@ -17,7 +17,6 @@
#ifndef SATFUNCBASE_HPP
#define SATFUNCBASE_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
@@ -79,11 +78,7 @@ namespace Opm
class SatFuncBase : public BlackoilPhases
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void init(Opm::DeckConstPtr newParserDeck,
void init(Opm::DeckConstPtr deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
@@ -126,7 +121,7 @@ namespace Opm
};
template <class TableType>
void SatFuncBase<TableType>::init(const EclipseGridParser& deck,
void SatFuncBase<TableType>::init(Opm::DeckConstPtr deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
@@ -135,139 +130,7 @@ namespace Opm
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
if (krw.front() != 0.0 || krow.back() != 0.0) {
OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)");
}
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sw.size();
std::vector<double> sw_ex(n+2);
std::vector<double> krw_ex(n+2);
std::vector<double> krow_ex(n+2);
std::vector<double> pcow_ex(n+2);
extendTable(sw,sw_ex,1);
extendTable(krw,krw_ex,0);
extendTable(krow,krow_ex,0);
extendTable(pcow,pcow_ex,0);
initializeTableType(krw_,sw_ex, krw_ex, samples);
initializeTableType(krow_,sw_ex, krow_ex, samples);
initializeTableType(pcow_,sw_ex, pcow_ex, samples);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
krwmax_ = krw.back();
kromax_ = krow.front();
swcr_ = swmax;
sowcr_ = 1.0 - swco;
krwr_ = krw.back();
krorw_ = krow.front();
for (std::vector<double>::size_type i=1; i<sw.size(); ++i) {
if (krw[i]> 0.0) {
swcr_ = sw[i-1];
krorw_ = krow[i-1];
break;
}
}
for (std::vector<double>::size_type i=sw.size()-1; i>=1; --i) {
if (krow[i-1]> 0.0) {
sowcr_ = 1.0 - sw[i];
krwr_ = krw[i];
break;
}
}
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sg.size();
std::vector<double> sg_ex(n+2);
std::vector<double> krg_ex(n+2);
std::vector<double> krog_ex(n+2);
std::vector<double> pcog_ex(n+2);
extendTable(sg,sg_ex,1);
extendTable(krg,krg_ex,0);
extendTable(krog,krog_ex,0);
extendTable(pcog,pcog_ex,0);
initializeTableType(krg_,sg_ex, krg_ex, samples);
initializeTableType(krog_,sg_ex, krog_ex, samples);
initializeTableType(pcog_,sg_ex, pcog_ex, samples);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
smin_[phase_usage.phase_pos[Vapour]] = sg.front();
krgmax_ = krg.back();
sgcr_ = sg.front();
sogcr_ = 1.0 - sg.back();
krgr_ = krg.back();
krorg_ = krg.front();
for (std::vector<double>::size_type i=1; i<sg.size(); ++i) {
if (krg[i]> 0.0) {
sgcr_ = sg[i-1];
krorg_ = krog[i-1];
break;
}
}
for (std::vector<double>::size_type i=sg.size()-1; i>=1; --i) {
if (krog[i-1]> 0.0) {
sogcr_ = 1.0 - sg[i];
krgr_ = krg[i];
break;
}
}
}
if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) {
sowcr_ -= smin_[phase_usage.phase_pos[Vapour]];
sogcr_ -= smin_[phase_usage.phase_pos[Aqua]];
smin_[phase_usage.phase_pos[Liquid]] = 0.0;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]
- smin_[phase_usage.phase_pos[Vapour]]; // First entry in SGOF-table supposed to be zero anyway ...
} else if (phase_usage.phase_used[Aqua]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]];
} else if (phase_usage.phase_used[Vapour]) {
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]];
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]];
}
}
template <class TableType>
void SatFuncBase<TableType>::init(Opm::DeckConstPtr newParserDeck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
Opm::SwofTable swof(newParserDeck->getKeyword("SWOF"), table_num);
Opm::SwofTable swof(deck->getKeyword("SWOF"), table_num);
const std::vector<double>& sw = swof.getSwColumn();
const std::vector<double>& krw = swof.getKrwColumn();
const std::vector<double>& krow = swof.getKrowColumn();
@@ -321,7 +184,7 @@ namespace Opm
}
}
if (phase_usage.phase_used[Vapour]) {
Opm::SgofTable sgof(newParserDeck->getKeyword("SGOF"), table_num);
Opm::SgofTable sgof(deck->getKeyword("SGOF"), table_num);
const std::vector<double>& sg = sgof.getSgColumn();
const std::vector<double>& krg = sgof.getKrgColumn();
const std::vector<double>& krog = sgof.getKrogColumn();
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
+7 -58
View File
@@ -22,7 +22,6 @@
#include <opm/core/props/satfunc/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
@@ -59,44 +58,15 @@ namespace Opm
/// to logical cartesian indices consistent with the deck.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
void init(const EclipseGridParser& deck,
void init(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid,
const int samples);
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
/// \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 newParserDeck,
const UnstructuredGrid& grid,
const int samples);
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] number_of_cells The number of cells of the grid to which property
/// object applies, needed for the
/// mapping from cell indices (typically from a processed
/// grid) to logical cartesian indices consistent with the
/// deck.
/// \param[in] global_cell The mapping from local cell indices of the grid to
/// global cell indices used in the deck.
/// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid.
/// \param[in] dimensions The dimensions of the grid.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// \tparam T The iterator Type for the cell centroids.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
template<class T>
void init(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples);
/// grid) to logical cartesian indices consistent with the
/// deck.
/// \param[in] global_cell The mapping from local cell indices of the grid to
@@ -106,7 +76,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.
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
void init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
@@ -178,41 +148,20 @@ namespace Opm
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
template<class T>
void initEPS(const EclipseGridParser& deck,
void initEPS(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSHyst(const EclipseGridParser& deck,
void initEPSHyst(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
template<class T>
void initEPS(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSHyst(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(Opm::DeckConstPtr newParserDeck,
void initEPSKey(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
@@ -237,10 +186,10 @@ namespace Opm
const std::vector<double>& krsr,
const std::vector<double>& krmax);
bool columnIsMasked_(Opm::DeckConstPtr newParserDeck,
bool columnIsMasked_(Opm::DeckConstPtr deck,
const std::string& keywordName,
int /* columnIdx */)
{ return newParserDeck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; }
{ return deck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; }
};
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
+86 -544
View File
@@ -49,30 +49,19 @@ namespace Opm
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
const UnstructuredGrid& grid,
const int samples)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
}
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const int samples)
{
init(newParserDeck, grid.number_of_cells,
grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
this->init(deck, grid.number_of_cells,
grid.global_cell, grid.cell_centroids,
grid.dimensions, samples);
}
/// Initialize from deck.
template <class SatFuncSet>
template< class T>
void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
@@ -86,16 +75,31 @@ namespace Opm
if (!phase_usage_.phase_used[Liquid]) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active.");
}
// Check SATOPTS status
bool hysteresis_switch = false;
if (deck->hasKeyword("SATOPTS")) {
const std::vector<std::string>& satopts = deck->getKeyword("SATOPTS")->getStringData();
for (size_t i = 0; i < satopts.size(); ++i) {
if (satopts[i] == std::string("HYSTER")) {
hysteresis_switch = true;
} else {
OPM_THROW(std::runtime_error, "Keyword SATOPTS: Switch " << satopts[i] << " not supported. ");
}
}
}
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
if (deck->hasKeyword("SATNUM")) {
const std::vector<int>& satnum = deck->getKeyword("SATNUM")->getIntData();
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
cell_to_func_.resize(number_of_cells);
for (int cell = 0; cell < number_of_cells; ++cell) {
const int deck_pos = (global_cell == NULL) ? cell : global_cell[cell];
const int num_cells = number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
cell_to_func_[cell] = satnum[deck_pos] - 1;
}
}
@@ -104,15 +108,13 @@ namespace Opm
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
num_tables = swof_table.size();
num_tables = deck->getKeyword("SWOF")->size();
if (num_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
int num_sgof_tables = sgof_table.size();
int num_sgof_tables = deck->getKeyword("SGOF")->size();
if (num_sgof_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
@@ -128,147 +130,34 @@ namespace Opm
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, samples);
}
// Saturation table scaling
do_hyst_ = false;
do_eps_ = false;
do_3pt_ = false;
if (deck.hasField("ENDSCALE")) {
//if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) {
// OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas.");
//}
if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted.");
}
if (deck.getENDSCALE().revers_switch_ != std::string("REVERS")) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'REVERS' accepted.");
}
if (deck.hasField("SCALECRS")) {
if (deck.getSCALECRS().scalecrs_ == std::string("YES")) {
do_3pt_ = true;
}
}
do_eps_ = true;
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions);
// For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR
do_hyst_ = deck.hasField("ISWL") || deck.hasField("ISWU") || deck.hasField("ISWCR") || deck.hasField("ISGL") ||
deck.hasField("ISGU") || deck.hasField("ISGCR") || deck.hasField("ISOWCR") || deck.hasField("ISOGCR");
if (do_hyst_) {
if (deck.hasField("KRW") || deck.hasField("KRG") || deck.hasField("KRO") || deck.hasField("KRWR") ||
deck.hasField("KRGR") || deck.hasField("KRORW") || deck.hasField("KRORG") ||
deck.hasField("IKRW") || deck.hasField("IKRG") || deck.hasField("IKRO") || deck.hasField("IKRWR") ||
deck.hasField("IKRGR") || deck.hasField("IKRORW") || deck.hasField("IKRORG") ) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently hysteresis and relperm value scaling can not be combined.");
}
initEPSHyst(deck, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
}
//OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Under construction ...");
}
}
/// Initialize from deck.
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples)
{
phase_usage_ = phaseUsageFromDeck(newParserDeck);
// Extract input data.
// Oil phase should be active.
if (!phase_usage_.phase_used[Liquid]) {
OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active.");
}
// Check SATOPTS status
bool hysteresis_switch = false;
if (newParserDeck->hasKeyword("SATOPTS")) {
const std::vector<std::string>& satopts = newParserDeck->getKeyword("SATOPTS")->getStringData();
for (size_t i = 0; i < satopts.size(); ++i) {
if (satopts[i] == std::string("HYSTER")) {
hysteresis_switch = true;
} else {
OPM_THROW(std::runtime_error, "Keyword SATOPTS: Switch " << satopts[i] << " not supported. ");
}
}
}
// Obtain SATNUM, if it exists, and create cell_to_func_.
// Otherwise, let the cell_to_func_ mapping be just empty.
int satfuncs_expected = 1;
if (newParserDeck->hasKeyword("SATNUM")) {
const std::vector<int>& satnum = newParserDeck->getKeyword("SATNUM")->getIntData();
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
const int num_cells = number_of_cells;
cell_to_func_.resize(num_cells);
const int* gc = global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int newParserDeck_pos = (gc == NULL) ? cell : gc[cell];
cell_to_func_[cell] = satnum[newParserDeck_pos] - 1;
}
}
// Find number of tables, check for consistency.
enum { Uninitialized = -1 };
int num_tables = Uninitialized;
if (phase_usage_.phase_used[Aqua]) {
num_tables = newParserDeck->getKeyword("SWOF")->size();
if (num_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
}
}
if (phase_usage_.phase_used[Vapour]) {
int num_sgof_tables = newParserDeck->getKeyword("SGOF")->size();
if (num_sgof_tables < satfuncs_expected) {
OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
}
if (num_tables == Uninitialized) {
num_tables = num_sgof_tables;
} else if (num_tables != num_sgof_tables) {
OPM_THROW(std::runtime_error, "Inconsistent number of tables in SWOF and SGOF.");
}
}
// Initialize tables.
satfuncset_.resize(num_tables);
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(newParserDeck, table, phase_usage_, samples);
}
// Check EHYSTR status
do_hyst_ = false;
if (hysteresis_switch && newParserDeck->hasKeyword("EHYSTR")) {
const int& relative_perm_hyst = newParserDeck->getKeyword("EHYSTR")->getRecord(0)->getItem(1)->getInt(0);
const std::string& limiting_hyst_flag = newParserDeck->getKeyword("EHYSTR")->getRecord(0)->getItem(4)->getString(0);
if (hysteresis_switch && deck->hasKeyword("EHYSTR")) {
const int& relative_perm_hyst = deck->getKeyword("EHYSTR")->getRecord(0)->getItem(1)->getInt(0);
const std::string& limiting_hyst_flag = deck->getKeyword("EHYSTR")->getRecord(0)->getItem(4)->getString(0);
if (relative_perm_hyst != int(0)) {
OPM_THROW(std::runtime_error, "Keyword EHYSTR, item 2: Flag '" << relative_perm_hyst << "' found, only '0' is supported. ");
}
if (limiting_hyst_flag != std::string("KR")) {
OPM_THROW(std::runtime_error, "Keyword EHYSTR, item 5: Flag '" << limiting_hyst_flag << "' found, only 'KR' is supported. ");
}
if ( ! newParserDeck->hasKeyword("ENDSCALE")) {
if ( ! deck->hasKeyword("ENDSCALE")) {
// TODO When use of IMBNUM is implemented, this constraint will be lifted.
OPM_THROW(std::runtime_error, "Currently hysteris effects is only available through endpoint scaling.");
}
do_hyst_ = true;
} else if (hysteresis_switch) {
OPM_THROW(std::runtime_error, "Switch HYSTER of keyword SATOPTS is active, but keyword EHYSTR not found.");
} else if (newParserDeck->hasKeyword("EHYSTR")) {
} else if (deck->hasKeyword("EHYSTR")) {
OPM_THROW(std::runtime_error, "Found keyword EHYSTR, but switch HYSTER of keyword SATOPTS is not set.");
}
// Saturation table scaling
do_eps_ = false;
do_3pt_ = false;
if (newParserDeck->hasKeyword("ENDSCALE")) {
Opm::EndscaleWrapper endscale(newParserDeck->getKeyword("ENDSCALE"));
if (deck->hasKeyword("ENDSCALE")) {
Opm::EndscaleWrapper endscale(deck->getKeyword("ENDSCALE"));
if (endscale.directionSwitch() != std::string("NODIR")) {
OPM_THROW(std::runtime_error,
"SaturationPropsFromDeck::init() -- ENDSCALE: "
@@ -279,8 +168,8 @@ namespace Opm
"SaturationPropsFromDeck::init() -- ENDSCALE: "
"Currently only 'REVERS' accepted.");
}
if (newParserDeck->hasKeyword("SCALECRS")) {
Opm::ScalecrsWrapper scalecrs(newParserDeck->getKeyword("SCALECRS"));
if (deck->hasKeyword("SCALECRS")) {
Opm::ScalecrsWrapper scalecrs(deck->getKeyword("SCALECRS"));
if (scalecrs.isEnabled()) {
do_3pt_ = true;
}
@@ -288,8 +177,8 @@ namespace Opm
do_eps_ = true;
// Make a consistency check of ENDNUM: #regions = NTENDP (ENDSCALE::3, TABDIMS::8)...
if (newParserDeck->hasKeyword("ENDNUM")) {
const std::vector<int>& endnum = newParserDeck->getKeyword("ENDNUM")->getIntData();
if (deck->hasKeyword("ENDNUM")) {
const std::vector<int>& endnum = deck->getKeyword("ENDNUM")->getIntData();
int endnum_regions = *std::max_element(endnum.begin(), endnum.end());
if (endnum_regions > endscale.numEndscaleTables()) {
OPM_THROW(std::runtime_error,
@@ -301,30 +190,29 @@ 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(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
initEPS(deck, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
if (do_hyst_) {
if (newParserDeck->hasKeyword("KRW")
|| newParserDeck->hasKeyword("KRG")
|| newParserDeck->hasKeyword("KRO")
|| newParserDeck->hasKeyword("KRWR")
|| newParserDeck->hasKeyword("KRGR")
|| newParserDeck->hasKeyword("KRORW")
|| newParserDeck->hasKeyword("KRORG")
|| newParserDeck->hasKeyword("ENKRVD")
|| newParserDeck->hasKeyword("IKRW")
|| newParserDeck->hasKeyword("IKRG")
|| newParserDeck->hasKeyword("IKRO")
|| newParserDeck->hasKeyword("IKRWR")
|| newParserDeck->hasKeyword("IKRGR")
|| newParserDeck->hasKeyword("IKRORW")
|| newParserDeck->hasKeyword("IKRORG") ) {
if (deck->hasKeyword("KRW")
|| deck->hasKeyword("KRG")
|| deck->hasKeyword("KRO")
|| deck->hasKeyword("KRWR")
|| deck->hasKeyword("KRGR")
|| deck->hasKeyword("KRORW")
|| deck->hasKeyword("KRORG")
|| deck->hasKeyword("ENKRVD")
|| deck->hasKeyword("IKRG")
|| deck->hasKeyword("IKRO")
|| deck->hasKeyword("IKRWR")
|| deck->hasKeyword("IKRGR")
|| deck->hasKeyword("IKRORW")
|| deck->hasKeyword("IKRORG") ) {
OPM_THROW(std::runtime_error,"Currently hysteresis and relperm value scaling cannot be combined.");
}
if (newParserDeck->hasKeyword("IMBNUM")) {
const std::vector<int>& imbnum = newParserDeck->getKeyword("IMBNUM")->getIntData();
if (deck->hasKeyword("IMBNUM")) {
const std::vector<int>& imbnum = deck->getKeyword("IMBNUM")->getIntData();
int imbnum_regions = *std::max_element(imbnum.begin(), imbnum.end());
if (imbnum_regions > num_tables) {
OPM_THROW(std::runtime_error,
@@ -343,7 +231,7 @@ namespace Opm
// to be a scaled version of the drainage curve (confer Norne model).
}
initEPSHyst(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
initEPSHyst(deck, number_of_cells, global_cell, begin_cell_centroids,
dimensions);
}
}
@@ -534,8 +422,8 @@ namespace Opm
// Initialize saturation scaling parameters
template <class SatFuncSet>
template <class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
@@ -615,93 +503,10 @@ namespace Opm
}
}
// Initialize saturation scaling parameters
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions)
{
std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
// Initialize saturation scaling parameter
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWL"), swl);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWU"), swu);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SWCR"), swcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGL"), sgl);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGU"), sgu);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SGCR"), sgcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOWCR"), sowcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("SOGCR"), sogcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRW"), krw);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRG"), krg);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRO"), kro);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRWR"), krwr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRGR"), krgr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORW"), krorw);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("KRORG"), krorg);
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_, swl, swcr, swu, sowcr, sgl, krwr, krw);
// ### 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_, swl, sowcr, swl, swcr, sgl, krorw, kro);
} 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_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
// ### 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_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
} 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_, swl, swcr, swu, sowcr, sgl, krwr, krw);
// ### 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_, swl, sowcr, swl, swcr, sgl, krorw, kro);
// ### 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_, sgl, sgcr, sgu, sogcr, swl, krgr, krg);
// ### 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_, sgl, sogcr, sgl, sgcr, swl, krorg, kro);
}
}
}
// Initialize hysteresis saturation scaling parameters
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(const EclipseGridParser& deck,
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
@@ -740,91 +545,6 @@ namespace Opm
std::string("IKRORW"), ikrorw);
initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORG"), ikrorg);
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];
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_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_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
// ### 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_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
} 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_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
// ### 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_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
} 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_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw);
// ### 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_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro);
// ### 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_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg);
// ### 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_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro);
}
}
}
// Initialize hysteresis saturation scaling parameters
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr newParserDeck,
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;
// Initialize hysteresis saturation scaling parameters
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWL"), iswl);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWU"), iswu);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISWCR"), iswcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGL"), isgl);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGU"), isgu);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISGCR"), isgcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOWCR"), isowcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("ISOGCR"), isogcr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRW"), ikrw);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRG"), ikrg);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRO"), ikro);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRWR"), ikrwr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRGR"), ikrgr);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORW"), ikrorw);
initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
std::string("IKRORG"), ikrorg);
eps_transf_hyst_.resize(number_of_cells);
sat_hyst_.resize(number_of_cells);
@@ -870,185 +590,7 @@ namespace Opm
// Initialize saturation scaling parameter
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam)
{
const bool useAqua = phase_usage_.phase_used[Aqua];
const bool useLiquid = phase_usage_.phase_used[Liquid];
const bool useVapour = phase_usage_.phase_used[Vapour];
bool useKeyword = deck.hasField(keyword);
bool hasENPTVD = deck.hasField("ENPTVD");
bool hasENKRVD = deck.hasField("ENKRVD");
int itab = 0;
std::vector<std::vector<std::vector<double> > > table_dummy;
std::vector<std::vector<std::vector<double> > >& table = table_dummy;
// 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 || deck.getENPTVD().mask_[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 || deck.getENPTVD().mask_[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 || deck.getENPTVD().mask_[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 || deck.getENPTVD().mask_[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 || deck.getENPTVD().mask_[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 || deck.getENPTVD().mask_[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 || deck.getENPTVD().mask_[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 || deck.getENPTVD().mask_[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 (!useKeyword && itab > 0) {
table = deck.getENPTVD().table_;
}
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[0])) {
itab = 1;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
}
} else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[1])) {
itab = 2;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgmax_;
}
} else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) {
if (useLiquid && (useKeyword || deck.getENKRVD().mask_[2])) {
itab = 3;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
}
} else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[3])) {
itab = 4;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
}
} else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[4])) {
itab = 5;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgr_;
}
} else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
if (useAqua && (useKeyword || deck.getENKRVD().mask_[5])) {
itab = 6;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
}
} else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
if (useVapour && (useKeyword || deck.getENKRVD().mask_[6])) {
itab = 7;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorg_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
table = deck.getENKRVD().table_;
}
}
if (scaleparam.empty()) {
return;
} else if (useKeyword) {
// Keyword values from deck
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
const std::vector<double>& val = deck.getFloatingPointValue(keyword);
for (int c = 0; c < int(scaleparam.size()); ++c) {
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
scaleparam[c] = val[deck_pos];
}
} else {
std::cout << "--- Scaling parameter '" << keyword << "' assigned via ";
if (keyword[0] == 'S')
deck.getENPTVD().write(std::cout);
else
deck.getENKRVD().write(std::cout);
for (int cell = 0; cell < number_of_cells; ++cell) {
int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
if (table[itab][jtab][0] != -1.0) {
std::vector<double>& depth = table[0][jtab];
std::vector<double>& val = table[itab][jtab];
double zc = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell,
dimensions),
dimensions-1);
if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth, val, zc);
}
}
}
}
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
template<class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(Opm::DeckConstPtr newParserDeck,
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(Opm::DeckConstPtr deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroid,
@@ -1059,9 +601,9 @@ namespace Opm
const bool useAqua = phase_usage_.phase_used[Aqua];
const bool useLiquid = phase_usage_.phase_used[Liquid];
const bool useVapour = phase_usage_.phase_used[Vapour];
bool useKeyword = newParserDeck->hasKeyword(keyword);
bool hasENPTVD = newParserDeck->hasKeyword("ENPTVD");
bool hasENKRVD = newParserDeck->hasKeyword("ENKRVD");
bool useKeyword = deck->hasKeyword(keyword);
bool hasENPTVD = deck->hasKeyword("ENPTVD");
bool hasENKRVD = deck->hasKeyword("ENKRVD");
int itab = 0;
std::vector<std::vector<double> > param_col;
std::vector<std::vector<double> > depth_col;
@@ -1072,56 +614,56 @@ namespace Opm
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_(newParserDeck, "ENPTVD", 0))) {
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_(newParserDeck, "ENPTVD", 1))) {
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_(newParserDeck, "ENPTVD", 2))) {
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_(newParserDeck, "ENPTVD", 3))) {
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_(newParserDeck, "ENPTVD", 4))) {
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_(newParserDeck, "ENPTVD", 5))) {
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_(newParserDeck, "ENPTVD", 6))) {
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_(newParserDeck, "ENPTVD", 7))) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENPTVD", 7))) {
itab = 8;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
@@ -1131,61 +673,61 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
int num_tables = newParserDeck->getKeyword("ENPTVD")->size();
int num_tables = deck->getKeyword("ENPTVD")->size();
param_col.resize(num_tables);
depth_col.resize(num_tables);
col_names.resize(9);
for (int table_num=0; table_num<num_tables; ++table_num) {
Opm::SingleRecordTable enptvd(newParserDeck->getKeyword("ENPTVD"), col_names, table_num);
Opm::SingleRecordTable enptvd(deck->getKeyword("ENPTVD"), col_names, table_num);
depth_col[table_num] = enptvd.getColumn(0); // depth
param_col[table_num] = enptvd.getColumn(itab); // itab=[1-8]: swl swcr swu sgl sgcr sgu sowcr sogcr
}
}
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 0))) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 0))) {
itab = 1;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
}
} else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 1))) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 1))) {
itab = 2;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgmax_;
}
} else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) {
if (useLiquid && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 2))) {
if (useLiquid && (useKeyword || columnIsMasked_(deck, "ENKRVD", 2))) {
itab = 3;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
}
} else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 3))) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 3))) {
itab = 4;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
}
} else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 4))) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 4))) {
itab = 5;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgr_;
}
} else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 5))) {
if (useAqua && (useKeyword || columnIsMasked_(deck, "ENKRVD", 5))) {
itab = 6;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
}
} else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 6))) {
if (useVapour && (useKeyword || columnIsMasked_(deck, "ENKRVD", 6))) {
itab = 7;
scaleparam.resize(number_of_cells);
for (int i=0; i<number_of_cells; ++i)
@@ -1195,12 +737,12 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
int num_tables = newParserDeck->getKeyword("ENKRVD")->size();
int num_tables = deck->getKeyword("ENKRVD")->size();
param_col.resize(num_tables);
depth_col.resize(num_tables);
col_names.resize(8);
for (int table_num=0; table_num<num_tables; ++table_num) {
Opm::SingleRecordTable enkrvd(newParserDeck->getKeyword("ENKRVD"), col_names, table_num);
Opm::SingleRecordTable enkrvd(deck->getKeyword("ENKRVD"), col_names, table_num);
depth_col[table_num] = enkrvd.getColumn(0); // depth
param_col[table_num] = enkrvd.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
}
@@ -1213,7 +755,7 @@ namespace Opm
// Keyword values from deck
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
const int* gc = global_cell;
const std::vector<double>& val = newParserDeck->getKeyword(keyword)->getSIDoubleData();
const std::vector<double>& val = deck->getKeyword(keyword)->getSIDoubleData();
for (int c = 0; c < int(scaleparam.size()); ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c];
scaleparam[c] = val[deck_pos];
@@ -1221,9 +763,9 @@ namespace Opm
} else {
const int dim = dimensions;
std::vector<int> endnum;
if ( newParserDeck->hasKeyword("ENDNUM")) {
if ( deck->hasKeyword("ENDNUM")) {
const std::vector<int>& e =
newParserDeck->getKeyword("ENDNUM")->getIntData();
deck->getKeyword("ENDNUM")->getIntData();
endnum.resize(number_of_cells);
const int* gc = global_cell;
for (int cell = 0; cell < number_of_cells; ++cell) {
opm/core/simulator/SimulatorTimer.cpp
-9
View File
@@ -21,7 +21,6 @@
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <ostream>
#include <numeric>
@@ -49,14 +48,6 @@ namespace Opm
total_time_ = num_psteps*stepsize;
}
/// Initialize from TSTEP field.
void SimulatorTimer::init(const EclipseGridParser& deck)
{
timesteps_ = deck.getTSTEP().tstep_;
total_time_ = std::accumulate(timesteps_.begin(), timesteps_.end(), 0.0);
start_date_ = deck.getStartDate();
}
/// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap
void SimulatorTimer::init(Opm::TimeMapConstPtr timeMap)
{
opm/core/simulator/SimulatorTimer.hpp
-6
View File
@@ -31,8 +31,6 @@ namespace Opm
{
namespace parameter { class ParameterGroup; }
class EclipseGridParser;
class SimulatorTimer
{
@@ -45,10 +43,6 @@ namespace Opm
/// stepsize_days (default 1)
void init(const parameter::ParameterGroup& param);
/// Initialize from TSTEP field.
/// Note that DATES are folded into TSTEP by the parser.
void init(const EclipseGridParser& deck);
/// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap
void init(TimeMapConstPtr timeMap);
opm/core/simulator/initState.hpp
+4 -47
View File
@@ -20,7 +20,7 @@
#ifndef OPM_INITSTATE_HEADER_INCLUDED
#define OPM_INITSTATE_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp> // DeckConstPtr
#include <opm/parser/eclipse/Deck/Deck.hpp>
struct UnstructuredGrid;
@@ -28,7 +28,6 @@ namespace Opm
{
namespace parameter { class ParameterGroup; }
class EclipseGridParser;
class IncompPropertiesInterface;
class BlackoilPropertiesInterface;
@@ -168,27 +167,7 @@ namespace Opm
template <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state);
/// Initialize a two-phase state from input deck.
/// If EQUIL is present:
/// - saturation is set according to the water-oil contact,
/// - pressure is set to hydrostatic equilibrium.
/// Otherwise:
/// - saturation is set according to SWAT,
/// - pressure is set according to PRESSURE.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
Opm::DeckConstPtr deck,
const double gravity,
State& state);
@@ -203,31 +182,9 @@ namespace Opm
template <class Props, class State>
void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
Opm::DeckConstPtr deck,
const double gravity,
State& state);
/// Initialize a two-phase water-oil blackoil state from input deck.
/// If EQUIL is present:
/// - saturation is set according to the water-oil contact,
/// - pressure is set to hydrostatic equilibrium.
/// Otherwise:
/// - saturation is set according to SWAT,
/// - pressure is set according to PRESSURE.
/// In addition, this function sets surfacevol.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initBlackoilStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state);
/// Initialize a blackoil state from input deck.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initBlackoilStateFromDeck(int number_of_cells,
@@ -238,7 +195,7 @@ namespace Opm
CCI begin_cell_centroids,
int dimensions,
const Props& props,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
const double gravity,
State& state);
} // namespace Opm
opm/core/simulator/initStateEquil.hpp
+18 -270
View File
@@ -22,7 +22,6 @@
#include <opm/core/simulator/EquilibrationHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/RegionMapping.hpp>
@@ -62,14 +61,7 @@ namespace Opm
*/
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck,
const double gravity,
BlackoilState& state);
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
const Opm::DeckConstPtr deck,
const double gravity,
BlackoilState& state);
@@ -188,61 +180,13 @@ namespace Opm
const std::vector<double> gas_saturation);
namespace DeckDependent {
inline
std::vector<EquilRecord>
getEquil(const EclipseGridParser& deck)
getEquil(const Opm::DeckConstPtr deck)
{
if (deck.hasField("EQUIL")) {
const EQUIL& eql = deck.getEQUIL();
typedef std::vector<EquilLine>::size_type sz_t;
const sz_t nrec = eql.equil.size();
std::vector<EquilRecord> ret;
ret.reserve(nrec);
for (sz_t r = 0; r < nrec; ++r) {
const EquilLine& rec = eql.equil[r];
EquilRecord record =
{
{ rec.datum_depth_ ,
rec.datum_depth_pressure_ }
,
{ rec.water_oil_contact_depth_ ,
rec.oil_water_cap_pressure_ }
,
{ rec.gas_oil_contact_depth_ ,
rec.gas_oil_cap_pressure_ }
,
rec.live_oil_table_index_
,
rec.wet_gas_table_index_
,
rec.N_
};
if (record.N != 0) {
OPM_THROW(std::domain_error,
"kw EQUIL, item 9: Only N=0 supported.");
}
ret.push_back(record);
}
return ret;
}
else {
OPM_THROW(std::domain_error,
"Deck does not provide equilibration data.");
}
}
inline
std::vector<EquilRecord>
getEquil(const Opm::DeckConstPtr newParserDeck)
{
if (newParserDeck->hasKeyword("EQUIL")) {
if (deck->hasKeyword("EQUIL")) {
Opm::EquilWrapper eql(newParserDeck->getKeyword("EQUIL"));
Opm::EquilWrapper eql(deck->getKeyword("EQUIL"));
const int nrec = eql.numRegions();
@@ -282,39 +226,16 @@ namespace Opm
}
}
inline
std::vector<int>
equilnum(const EclipseGridParser& deck,
equilnum(const Opm::DeckConstPtr deck,
const UnstructuredGrid& G )
{
std::vector<int> eqlnum;
if (deck.hasField("EQLNUM")) {
const std::vector<int>& e = deck.getIntegerValue("EQLNUM");
eqlnum.reserve(e.size());
std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
std::bind2nd(std::minus<int>(), 1));
}
else {
// No explicit equilibration region.
// All cells in region zero.
eqlnum.assign(G.number_of_cells, 0);
}
return eqlnum;
}
inline
std::vector<int>
equilnum(const Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& G )
{
std::vector<int> eqlnum;
if (newParserDeck->hasKeyword("EQLNUM")) {
if (deck->hasKeyword("EQLNUM")) {
eqlnum.resize(G.number_of_cells);
const std::vector<int>& e =
newParserDeck->getKeyword("EQLNUM")->getIntData();
deck->getKeyword("EQLNUM")->getIntData();
const int* gc = G.global_cell;
for (int cell = 0; cell < G.number_of_cells; ++cell) {
const int deck_pos = (gc == NULL) ? cell : gc[cell];
@@ -331,14 +252,10 @@ namespace Opm
}
template <class InputDeck>
class InitialStateComputer;
template <>
class InitialStateComputer<Opm::EclipseGridParser> {
class InitialStateComputer {
public:
InitialStateComputer(const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck ,
const Opm::DeckConstPtr deck,
const UnstructuredGrid& G ,
const double grav = unit::gravity)
: pp_(props.numPhases(),
@@ -352,186 +269,17 @@ namespace Opm
const std::vector<EquilRecord> rec = getEquil(deck);
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(deck, G));
const RegionMapping<> eqlmap(equilnum(deck, G));
// Create Rs functions.
rs_func_.reserve(rec.size());
if (deck.hasField("DISGAS")) {
if (deck->hasKeyword("DISGAS")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].live_oil_table_index > 0) {
if (deck.hasField("RSVD")) {
// TODO When this kw is actually parsed, also check for proper number of available tables
// For now, just use dummy ...
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
std::vector<double> rs; rs.push_back(0.0); rs.push_back(100.0);
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, depth, rs));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RSVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press;
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
rv_func_.reserve(rec.size());
if (deck.hasField("VAPOIL")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].wet_gas_table_index > 0) {
if (deck.hasField("RVVD")) {
// TODO When this kw is actually parsed, also check for proper number of available tables
// For now, just use dummy ...
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
std::vector<double> rv; rv.push_back(0.0); rv.push_back(0.0001);
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, depth, rv));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RVVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press + rec[i].goc.press;
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, props, G, grav);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
}
typedef std::vector<double> Vec;
typedef std::vector<Vec> PVec; // One per phase.
const PVec& press() const { return pp_; }
const PVec& saturation() const { return sat_; }
const Vec& rs() const { return rs_; }
const Vec& rv() const { return rv_; }
private:
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
typedef EquilReg<RhoCalc> EqReg;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
PVec pp_;
PVec sat_;
Vec rs_;
Vec rv_;
template <class RMap>
void
calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec ,
const Opm::BlackoilPropertiesInterface& props,
const UnstructuredGrid& G ,
const double grav)
{
typedef Miscibility::NoMixing NoMix;
for (typename RMap::RegionId
r = 0, nr = reg.numRegions();
r < nr; ++r)
{
const typename RMap::CellRange cells = reg.cells(r);
const int repcell = *cells.begin();
const RhoCalc calc(props, repcell);
const EqReg eqreg(rec[r], calc,
rs_func_[r], rv_func_[r],
props.phaseUsage());
PVec press = phasePressures(G, eqreg, cells, grav);
const PVec sat = phaseSaturations(eqreg, cells, props, press);
const int np = props.numPhases();
for (int p = 0; p < np; ++p) {
copyFromRegion(press[p], cells, pp_[p]);
copyFromRegion(sat[p], cells, sat_[p]);
}
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs, cells, rs_);
copyFromRegion(rv, cells, rv_);
}
}
}
template <class CellRangeType>
void copyFromRegion(const Vec& source,
const CellRangeType& cells,
Vec& destination)
{
auto s = source.begin();
auto c = cells.begin();
const auto e = cells.end();
for (; c != e; ++c, ++s) {
destination[*c] = *s;
}
}
};
template <>
class InitialStateComputer<Opm::DeckConstPtr> {
public:
InitialStateComputer(const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& G ,
const double grav = unit::gravity)
: pp_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
sat_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
rs_(G.number_of_cells),
rv_(G.number_of_cells)
{
// Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(newParserDeck);
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(newParserDeck, G));
// Create Rs functions.
rs_func_.reserve(rec.size());
if (newParserDeck->hasKeyword("DISGAS")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].live_oil_table_index > 0) {
if (newParserDeck->hasKeyword("RSVD") &&
size_t(rec[i].live_oil_table_index) <= newParserDeck->getKeyword("RSVD")->size()) {
Opm::SingleRecordTable rsvd(newParserDeck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
if (deck->hasKeyword("RSVD") &&
size_t(rec[i].live_oil_table_index) <= deck->getKeyword("RSVD")->size()) {
Opm::SingleRecordTable rsvd(deck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
std::vector<double> vd(rsvd.getColumn("vd"));
std::vector<double> rs(rsvd.getColumn("rs"));
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, vd, rs));
@@ -556,13 +304,13 @@ namespace Opm
}
rv_func_.reserve(rec.size());
if (newParserDeck->hasKeyword("VAPOIL")) {
if (deck->hasKeyword("VAPOIL")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].wet_gas_table_index > 0) {
if (newParserDeck->hasKeyword("RVVD") &&
size_t(rec[i].wet_gas_table_index) <= newParserDeck->getKeyword("RVVD")->size()) {
Opm::SingleRecordTable rvvd(newParserDeck->getKeyword("RVVD"),std::vector<std::string>{"vd", "rv"},rec[i].wet_gas_table_index-1);
if (deck->hasKeyword("RVVD") &&
size_t(rec[i].wet_gas_table_index) <= deck->getKeyword("RVVD")->size()) {
Opm::SingleRecordTable rvvd(deck->getKeyword("RVVD"),std::vector<std::string>{"vd", "rv"},rec[i].wet_gas_table_index-1);
std::vector<double> vd(rvvd.getColumn("vd"));
std::vector<double> rv(rvvd.getColumn("rv"));
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, vd, rv));
opm/core/simulator/initStateEquil_impl.hpp
+2 -22
View File
@@ -737,11 +737,11 @@ namespace Opm
*/
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck,
const Opm::DeckConstPtr deck,
const double gravity,
BlackoilState& state)
{
typedef Equil::DeckDependent::InitialStateComputer<EclipseGridParser> ISC;
typedef Equil::DeckDependent::InitialStateComputer ISC;
ISC isc(props, deck, grid, gravity);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
@@ -753,26 +753,6 @@ namespace Opm
state.rv() = isc.rv();
// TODO: state.surfacevol() must be computed from s, rs, rv.
}
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
const double gravity,
BlackoilState& state)
{
typedef Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> ISC;
ISC isc(props, newParserDeck, grid, gravity);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
? pu.phase_pos[BlackoilPhases::Liquid]
: pu.phase_pos[BlackoilPhases::Aqua];
state.pressure() = isc.press()[ref_phase];
state.saturation() = convertSats(isc.saturation());
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
// TODO: state.surfacevol() must be computed from s, rs, rv.
}
opm/core/simulator/initState_impl.hpp
+23 -202
View File
@@ -24,7 +24,6 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/MonotCubicInterpolator.hpp>
#include <opm/core/utility/Units.hpp>
@@ -553,14 +552,14 @@ namespace Opm
template <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
const double gravity,
State& state)
{
initStateFromDeck(grid.number_of_cells, grid.global_cell,
grid.number_of_faces, UgGridHelpers::faceCells(grid),
grid.face_centroids, grid.cell_centroids, grid.dimensions,
props, newParserDeck, gravity, state);
props, deck, gravity, state);
}
/// Initialize a state from input deck.
template <class FaceCells, class FCI, class CCI, class Props, class State>
@@ -572,23 +571,23 @@ namespace Opm
CCI begin_cell_centroids,
int dimensions,
const Props& props,
Opm::DeckConstPtr newParserDeck,
Opm::DeckConstPtr deck,
const double gravity,
State& state)
{
const int num_phases = props.numPhases();
const PhaseUsage pu = phaseUsageFromDeck(newParserDeck);
const PhaseUsage pu = phaseUsageFromDeck(deck);
if (num_phases != pu.num_phases) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
state.init(number_of_cells, number_of_faces, num_phases);
if (newParserDeck->hasKeyword("EQUIL") && newParserDeck->hasKeyword("PRESSURE")) {
if (deck->hasKeyword("EQUIL") && deck->hasKeyword("PRESSURE")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
"condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
}
if (newParserDeck->hasKeyword("EQUIL")) {
if (deck->hasKeyword("EQUIL")) {
if (num_phases != 2) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
}
@@ -596,7 +595,7 @@ namespace Opm
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
}
// Set saturations depending on oil-water contact.
EquilWrapper equil(newParserDeck->getKeyword("EQUIL"));
EquilWrapper equil(deck->getKeyword("EQUIL"));
if (equil.numRegions() != 1) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
}
@@ -608,19 +607,19 @@ namespace Opm
const double datum_p = equil.datumDepthPressure(0);
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
props, woc, gravity, datum_z, datum_p, state);
} else if (newParserDeck->hasKeyword("PRESSURE")) {
} else if (deck->hasKeyword("PRESSURE")) {
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
std::vector<double>& s = state.saturation();
std::vector<double>& p = state.pressure();
const std::vector<double>& p_deck = newParserDeck->getKeyword("PRESSURE")->getSIDoubleData();
const std::vector<double>& p_deck = deck->getKeyword("PRESSURE")->getSIDoubleData();
const int num_cells = number_of_cells;
if (num_phases == 2) {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
// oil-gas: we require SGAS
if (!newParserDeck->hasKeyword("SGAS")) {
if (!deck->hasKeyword("SGAS")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init");
}
const std::vector<double>& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData();
const std::vector<double>& sg_deck = deck->getKeyword("SGAS")->getSIDoubleData();
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
for (int c = 0; c < num_cells; ++c) {
@@ -631,10 +630,10 @@ namespace Opm
}
} else {
// water-oil or water-gas: we require SWAT
if (!newParserDeck->hasKeyword("SWAT")) {
if (!deck->hasKeyword("SWAT")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
}
const std::vector<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
const std::vector<double>& sw_deck = deck->getKeyword("SWAT")->getSIDoubleData();
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int nwpos = (wpos + 1) % 2;
for (int c = 0; c < num_cells; ++c) {
@@ -645,15 +644,15 @@ namespace Opm
}
}
} else if (num_phases == 3) {
const bool has_swat_sgas = newParserDeck->hasKeyword("SWAT") && newParserDeck->hasKeyword("SGAS");
const bool has_swat_sgas = deck->hasKeyword("SWAT") && deck->hasKeyword("SGAS");
if (!has_swat_sgas) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
}
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const std::vector<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
const std::vector<double>& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData();
const std::vector<double>& sw_deck = deck->getKeyword("SWAT")->getSIDoubleData();
const std::vector<double>& sg_deck = deck->getKeyword("SGAS")->getSIDoubleData();
for (int c = 0; c < num_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
s[3*c + wpos] = sw_deck[c_deck];
@@ -673,129 +672,6 @@ namespace Opm
begin_cell_centroids, state);
}
/// Initialize a state from input deck.
template <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state)
{
initStateFromDeck(grid.number_of_cells, grid.global_cell,
grid.number_of_faces, UgGridHelpers::faceCells(grid),
grid.face_centroids, grid.cell_centroids,
grid.dimensions, props, deck,
gravity, state);
}
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
const double gravity,
State& state)
{
const int num_phases = props.numPhases();
const PhaseUsage pu = phaseUsageFromDeck(deck);
if (num_phases != pu.num_phases) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
if (deck.hasField("EQUIL") && deck.hasField("PRESSURE")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
"condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
}
state.init(number_of_cells, number_of_faces, num_phases);
if (deck.hasField("EQUIL")) {
if (num_phases != 2) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
}
// Set saturations depending on oil-water contact.
const EQUIL& equil= deck.getEQUIL();
if (equil.equil.size() != 1) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
}
const double woc = equil.equil[0].water_oil_contact_depth_;
initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
props, woc, WaterBelow, state);
// Set pressure depending on densities and depths.
const double datum_z = equil.equil[0].datum_depth_;
const double datum_p = equil.equil[0].datum_depth_pressure_;
initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
props, woc, gravity, datum_z, datum_p, state);
} else if (deck.hasField("PRESSURE")) {
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
std::vector<double>& s = state.saturation();
std::vector<double>& p = state.pressure();
const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
if (num_phases == 2) {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
// oil-gas: we require SGAS
if (!deck.hasField("SGAS")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init");
}
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
s[2*c + gpos] = sg_deck[c_deck];
s[2*c + opos] = 1.0 - sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
// water-oil or water-gas: we require SWAT
if (!deck.hasField("SWAT")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
}
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int nwpos = (wpos + 1) % 2;
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
s[2*c + wpos] = sw_deck[c_deck];
s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
p[c] = p_deck[c_deck];
}
}
} else if (num_phases == 3) {
const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS");
if (!has_swat_sgas) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
}
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
s[3*c + wpos] = sw_deck[c_deck];
s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
s[3*c + gpos] = sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
OPM_THROW(std::runtime_error, "initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
}
} else {
OPM_THROW(std::runtime_error, "initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
}
// Finally, init face pressures.
initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
begin_cell_centroids, state);
}
/// Initialize surface volume from pressure and saturation by z = As.
/// Here saturation is used as an intial guess for z in the
/// computation of A.
@@ -975,12 +851,11 @@ namespace Opm
}
}
/// Initialize a blackoil state from input deck.
template <class Props, class State>
void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
const EclipseGridParser& deck,
Opm::DeckConstPtr deck,
const double gravity,
State& state)
{
@@ -990,6 +865,7 @@ namespace Opm
props, deck, gravity, state);
}
/// Initialize a blackoil state from input deck.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initBlackoilStateFromDeck(int number_of_cells,
const int* global_cell,
@@ -999,69 +875,15 @@ namespace Opm
CCI begin_cell_centroids,
int dimensions,
const Props& props,
const EclipseGridParser& deck,
Opm::DeckConstPtr deck,
const double gravity,
State& state)
{
initStateFromDeck(number_of_cells, global_cell, number_of_faces,
face_cells, begin_face_centroids, begin_cell_centroids,
dimensions, props, deck, gravity, state);
if (deck.hasField("RS")) {
const std::vector<double>& rs_deck = deck.getFloatingPointValue("RS");
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
state.gasoilratio()[c] = rs_deck[c_deck];
}
initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
computeSaturation(props,state);
} else if (deck.hasField("RV")){
const std::vector<double>& rv_deck = deck.getFloatingPointValue("RV");
for (int c = 0; c < number_of_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
state.rv()[c] = rv_deck[c_deck];
}
initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
computeSaturation(props,state);
}
else {
OPM_THROW(std::runtime_error, "Temporarily, we require the RS or the RV field.");
}
}
/// Initialize a blackoil state from input deck.
template <class Props, class State>
void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
Opm::DeckConstPtr newParserDeck,
const double gravity,
State& state)
{
initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell,
grid.number_of_faces, UgGridHelpers::faceCells(grid),
grid.face_centroids, grid.cell_centroids,grid.dimensions,
props, newParserDeck, gravity, state);
}
/// Initialize a blackoil state from input deck.
template <class FaceCells, class FCI, class CCI, class Props, class State>
void initBlackoilStateFromDeck(int number_of_cells,
const int* global_cell,
int number_of_faces,
FaceCells face_cells,
FCI begin_face_centroids,
CCI begin_cell_centroids,
int dimensions,
const Props& props,
Opm::DeckConstPtr newParserDeck,
const double gravity,
State& state)
{
initStateFromDeck(number_of_cells, global_cell, number_of_faces,
face_cells, begin_face_centroids, begin_cell_centroids,
dimensions, props, newParserDeck, gravity, state);
if (newParserDeck->hasKeyword("RS")) {
const std::vector<double>& rs_deck = newParserDeck->getKeyword("RS")->getSIDoubleData();
if (deck->hasKeyword("RS")) {
const std::vector<double>& rs_deck = deck->getKeyword("RS")->getSIDoubleData();
const int num_cells = number_of_cells;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
@@ -1069,8 +891,8 @@ namespace Opm
}
initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
computeSaturation(props,state);
} else if (newParserDeck->hasKeyword("RV")){
const std::vector<double>& rv_deck = newParserDeck->getKeyword("RV")->getSIDoubleData();
} else if (deck->hasKeyword("RV")){
const std::vector<double>& rv_deck = deck->getKeyword("RV")->getSIDoubleData();
const int num_cells = number_of_cells;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (global_cell == NULL) ? c : global_cell[c];
@@ -1084,7 +906,6 @@ namespace Opm
}
}
} // namespace Opm
opm/core/wells/WellCollection.cpp
-44
View File
@@ -79,50 +79,6 @@ namespace Opm
child->setParent(parent);
}
void WellCollection::addChild(const std::string& child_name,
const std::string& parent_name,
const EclipseGridParser& deck)
{
WellsGroupInterface* parent = findNode(parent_name);
if (!parent) {
roots_.push_back(createWellsGroup(parent_name, deck));
parent = roots_[roots_.size() - 1].get();
}
std::shared_ptr<WellsGroupInterface> child;
for (size_t i = 0; i < roots_.size(); ++i) {
if (roots_[i]->name() == child_name) {
child = roots_[i];
// We've found a new parent to the previously thought root, need to remove it
for(size_t j = i; j < roots_.size() - 1; ++j) {
roots_[j] = roots_[j+1];
}
roots_.resize(roots_.size()-1);
break;
}
}
if (!child.get()) {
child = createWellsGroup(child_name, deck);
}
WellsGroup* parent_as_group = static_cast<WellsGroup*> (parent);
if (!parent_as_group) {
OPM_THROW(std::runtime_error, "Trying to add child to group named " << parent_name << ", but it's not a group.");
}
parent_as_group->addChild(child);
if(child->isLeafNode()) {
leaf_nodes_.push_back(static_cast<WellNode*>(child.get()));
}
child->setParent(parent);
}
const std::vector<WellNode*>& WellCollection::getLeafNodes() const {
return leaf_nodes_;
}
opm/core/wells/WellCollection.hpp
-14
View File
@@ -16,9 +16,6 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLCOLLECTION_HPP
#define OPM_WELLCOLLECTION_HPP
@@ -27,7 +24,6 @@
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/core/grid.h>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
@@ -47,16 +43,6 @@ namespace Opm
void addGroup(GroupConstPtr groupChild, std::string parent_name,
size_t timeStep, const PhaseUsage& phaseUsage);
/// Adds and creates if necessary the child to the collection
/// and appends it to parent's children. Also adds and creates the parent
/// if necessary.
/// \param[in] child name of child node
/// \param[in] parent name of parent node
/// \param[in] deck deck from which we will extract group control data
void addChild(const std::string& child,
const std::string& parent,
const EclipseGridParser& deck);
/// Adds the child to the collection
/// and appends it to parent's children.
/// \param[in] child the child node
opm/core/wells/WellsGroup.cpp
-95
View File
@@ -1045,101 +1045,6 @@ namespace Opm
}
} // anonymous namespace
std::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name,
const EclipseGridParser& deck)
{
PhaseUsage phase_usage = phaseUsageFromDeck(deck);
std::shared_ptr<WellsGroupInterface> return_value;
// First we need to determine whether it's a group or just a well:
bool isWell = false;
if (deck.hasField("WELSPECS")) {
WELSPECS wspecs = deck.getWELSPECS();
for (size_t i = 0; i < wspecs.welspecs.size(); i++) {
if (wspecs.welspecs[i].name_ == name) {
isWell = true;
break;
}
}
}
// For now, assume that if it isn't a well, it's a group
if (isWell) {
ProductionSpecification production_specification;
InjectionSpecification injection_specification;
if (deck.hasField("WCONINJE")) {
WCONINJE wconinje = deck.getWCONINJE();
for (size_t i = 0; i < wconinje.wconinje.size(); i++) {
if (wconinje.wconinje[i].well_ == name) {
WconinjeLine line = wconinje.wconinje[i];
injection_specification.BHP_limit_ = line.BHP_limit_;
injection_specification.injector_type_ = toInjectorType(line.injector_type_);
injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
injection_specification.reservoir_flow_max_rate_ = line.reservoir_flow_max_rate_;
production_specification.guide_rate_ = 0.0; // We know we're not a producer
}
}
}
if (deck.hasField("WCONPROD")) {
WCONPROD wconprod = deck.getWCONPROD();
for (size_t i = 0; i < wconprod.wconprod.size(); i++) {
if (wconprod.wconprod[i].well_ == name) {
WconprodLine line = wconprod.wconprod[i];
production_specification.BHP_limit_ = line.BHP_limit_;
production_specification.reservoir_flow_max_rate_ = line.reservoir_flow_max_rate_;
production_specification.oil_max_rate_ = line.oil_max_rate_;
production_specification.control_mode_ = toProductionControlMode(line.control_mode_);
production_specification.water_max_rate_ = line.water_max_rate_;
injection_specification.guide_rate_ = 0.0; // we know we're not an injector
}
}
}
return_value.reset(new WellNode(name, production_specification, injection_specification, phase_usage));
} else {
InjectionSpecification injection_specification;
if (deck.hasField("GCONINJE")) {
GCONINJE gconinje = deck.getGCONINJE();
for (size_t i = 0; i < gconinje.gconinje.size(); i++) {
if (gconinje.gconinje[i].group_ == name) {
GconinjeLine line = gconinje.gconinje[i];
injection_specification.injector_type_ = toInjectorType(line.injector_type_);
injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
injection_specification.reservoir_flow_max_rate_ = line.resv_flow_max_rate_;
injection_specification.reinjection_fraction_target_ = line.reinjection_fraction_target_;
injection_specification.voidage_replacment_fraction_ = line.voidage_replacement_fraction_;
}
}
}
ProductionSpecification production_specification;
if (deck.hasField("GCONPROD")) {
std::cout << "Searching in gconprod " << std::endl;
std::cout << "name= " << name << std::endl;
GCONPROD gconprod = deck.getGCONPROD();
for (size_t i = 0; i < gconprod.gconprod.size(); i++) {
if (gconprod.gconprod[i].group_ == name) {
GconprodLine line = gconprod.gconprod[i];
production_specification.oil_max_rate_ = line.oil_max_rate_;
std::cout << "control_mode = " << line.control_mode_ << std::endl;
production_specification.control_mode_ = toProductionControlMode(line.control_mode_);
production_specification.water_max_rate_ = line.water_max_rate_;
production_specification.gas_max_rate_ = line.gas_max_rate_;
production_specification.liquid_max_rate_ = line.liquid_max_rate_;
production_specification.procedure_ = toProductionProcedure(line.procedure_);
production_specification.reservoir_flow_max_rate_ = line.resv_max_rate_;
}
}
}
return_value.reset(new WellsGroup(name, production_specification, injection_specification, phase_usage));
}
return return_value;
}
std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep, const PhaseUsage& phase_usage )
{
InjectionSpecification injection_specification;
opm/core/wells/WellsGroup.hpp
-7
View File
@@ -22,7 +22,6 @@
#include <opm/core/wells/InjectionSpecification.hpp>
#include <opm/core/wells/ProductionSpecification.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/BlackoilPhases.hpp>
@@ -403,12 +402,6 @@ namespace Opm
bool shut_well_;
};
/// Creates the WellsGroupInterface for the given name
/// \param[in] name the name of the wells group.
/// \param[in] deck the deck from which to fetch information.
std::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name,
const EclipseGridParser& deck);
/// Creates the WellsGroupInterface for the given well
/// \param[in] well the Well to construct object for
/// \param[in] timeStep the time step in question
opm/core/wells/WellsManager.cpp
-1
View File
@@ -21,7 +21,6 @@
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
opm/core/wells/WellsManager.hpp
-2
View File
@@ -33,8 +33,6 @@ struct UnstructuredGrid;
namespace Opm
{
class EclipseGridParser;
struct WellData
{
WellType type;
tests/test_blackoilfluid.cpp
+19 -19
View File
@@ -36,7 +36,7 @@ using namespace Opm;
using namespace std;
std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr newParserDeck,PhaseUsage phase_usage_){
std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr deck,PhaseUsage phase_usage_){
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
@@ -49,8 +49,8 @@ std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr new
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (newParserDeck->hasKeyword("PVTW")) {
Opm::PvtwTable pvtwTable(newParserDeck->getKeyword("PVTW"));
if (deck->hasKeyword("PVTW")) {
Opm::PvtwTable pvtwTable(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
} else {
// Eclipse 100 default.
@@ -60,19 +60,19 @@ std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr new
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (newParserDeck->hasKeyword("PVDO")) {
Opm::PvdoTable pvdoTable(newParserDeck->getKeyword("PVDO"));
if (deck->hasKeyword("PVDO")) {
Opm::PvdoTable pvdoTable(deck->getKeyword("PVDO"));
if (samples > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
} else {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
}
} else if (newParserDeck->hasKeyword("PVTO")) {
Opm::PvtoTable pvtoTable(newParserDeck->getKeyword("PVTO"), /*tableIdx=*/0);
} else if (deck->hasKeyword("PVTO")) {
Opm::PvtoTable pvtoTable(deck->getKeyword("PVTO"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
} else if (newParserDeck->hasKeyword("PVCDO")) {
Opm::PvcdoTable pvcdoTable(newParserDeck->getKeyword("PVCDO"), /*tableIdx=*/0);
} else if (deck->hasKeyword("PVCDO")) {
Opm::PvcdoTable pvcdoTable(deck->getKeyword("PVCDO"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
} else {
@@ -81,16 +81,16 @@ std::vector<std::shared_ptr<SinglePvtInterface> > getProps(Opm::DeckConstPtr new
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
if (newParserDeck->hasKeyword("PVDG")) {
Opm::PvdgTable pvdgTable(newParserDeck->getKeyword("PVDG"));
if (deck->hasKeyword("PVDG")) {
Opm::PvdgTable pvdgTable(deck->getKeyword("PVDG"));
if (samples > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
} else {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
}
} else if (newParserDeck->hasKeyword("PVTG")) {
Opm::PvtgTable pvtgTable(newParserDeck->getKeyword("PVTG"), /*tableIdx=*/0);
} else if (deck->hasKeyword("PVTG")) {
Opm::PvtgTable pvtgTable(deck->getKeyword("PVTG"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
} else {
@@ -232,11 +232,11 @@ BOOST_AUTO_TEST_CASE(test_liveoil)
const std::string filename = "liveoil.DATA";
cout << "Reading deck: " << filename << endl;
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(newParserDeck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(newParserDeck,phase_usage_);
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
@@ -307,11 +307,11 @@ BOOST_AUTO_TEST_CASE(test_wetgas)
const std::string filename = "wetgas.DATA";
cout << "Reading deck: " << filename << endl;
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(newParserDeck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(newParserDeck,phase_usage_);
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<SinglePvtInterface> > props_ = getProps(deck,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
tests/test_blackoilstate.cpp
+8 -8
View File
@@ -32,12 +32,12 @@ BOOST_AUTO_TEST_CASE(EqualsDifferentDeckReturnFalse) {
const string filename1 = "testBlackoilState1.DATA";
const string filename2 = "testBlackoilState2.DATA";
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck1(parser->parseFile(filename1));
Opm::DeckConstPtr newParserDeck2(parser->parseFile(filename2));
Opm::DeckConstPtr deck1(parser->parseFile(filename1));
Opm::DeckConstPtr deck2(parser->parseFile(filename2));
GridManager gridManager1(newParserDeck1);
GridManager gridManager1(deck1);
const UnstructuredGrid* grid1 = gridManager1.c_grid();
GridManager gridManager2(newParserDeck2);
GridManager gridManager2(deck2);
const UnstructuredGrid* grid2 = gridManager2.c_grid();
BlackoilState state1;
@@ -55,9 +55,9 @@ BOOST_AUTO_TEST_CASE(EqualsDifferentNumPhasesReturnFalse) {
const string filename = "testBlackoilState1.DATA";
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
GridManager gridManager(newParserDeck);
GridManager gridManager(deck);
const UnstructuredGrid* grid = gridManager.c_grid();
BlackoilState state1;
@@ -75,9 +75,9 @@ BOOST_AUTO_TEST_CASE(EqualsNumericalDifferenceReturnFalse) {
const string filename = "testBlackoilState1.DATA";
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
GridManager gridManager(newParserDeck);
GridManager gridManager(deck);
const UnstructuredGrid* grid = gridManager.c_grid();
BlackoilState state1;
tests/test_column_extract.cpp
+2 -2
View File
@@ -131,9 +131,9 @@ BOOST_AUTO_TEST_CASE(DisjointColumn)
correct_answer[9].resize(1);
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseString(grdecl));
Opm::DeckConstPtr deck(parser->parseString(grdecl));
Opm::GridManager manager(newParserDeck);
Opm::GridManager manager(deck);
VVI columns;
Opm::extractColumn(*manager.c_grid(), columns);
tests/test_equil.cpp
+6 -6
View File
@@ -337,7 +337,7 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("deadfluids.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, *grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, *grid, 10.0);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, *grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
@@ -416,7 +416,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 10.0);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -455,7 +455,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
Opm::DeckConstPtr deck = parser->parseFile("capillary_overlap.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -516,7 +516,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
Opm::DeckConstPtr deck = parser->parseFile("equil_liveoil.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -594,7 +594,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
Opm::DeckConstPtr deck = parser->parseFile("equil_livegas.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -675,7 +675,7 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
Opm::DeckConstPtr deck = parser->parseFile("equil_rsvd_and_rvvd.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
tests/test_wellsmanager.cpp
+12 -12
View File
@@ -176,10 +176,10 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
const std::string filename = "wells_manager_data.data";
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParserDeck));
Opm::GridManager gridManager(newParserDeck);
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::GridManager gridManager(deck);
{
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
@@ -199,10 +199,10 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
BOOST_AUTO_TEST_CASE(WellsEqual) {
const std::string filename = "wells_manager_data.data";
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParserDeck));
Opm::GridManager gridManager(newParserDeck);
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::GridManager gridManager(deck);
Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
@@ -215,10 +215,10 @@ BOOST_AUTO_TEST_CASE(WellsEqual) {
BOOST_AUTO_TEST_CASE(ControlsEqual) {
const std::string filename = "wells_manager_data.data";
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParserDeck));
Opm::GridManager gridManager(newParserDeck);
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::GridManager gridManager(deck);
Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
@@ -239,10 +239,10 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
BOOST_AUTO_TEST_CASE(WellHasSTOP_ExceptionIsThrown) {
const std::string filename = "wells_manager_data_wellSTOP.data";
Opm::ParserPtr parser(new Opm::Parser());
Opm::DeckConstPtr newParserDeck(parser->parseFile(filename));
Opm::DeckConstPtr deck(parser->parseFile(filename));
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParserDeck));
Opm::GridManager gridManager(newParserDeck);
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck));
Opm::GridManager gridManager(deck);
BOOST_CHECK_THROW( new Opm::WellsManager(eclipseState, 0, *gridManager.c_grid(), NULL), std::runtime_error );
}