diff --git a/examples/sim_2p_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp index bc688b3f..4b42abb1 100644 --- a/examples/sim_2p_incomp_reorder.cpp +++ b/examples/sim_2p_incomp_reorder.cpp @@ -81,9 +81,7 @@ main(int argc, char** argv) // Grid init grid.reset(new GridManager(*deck)); // Rock and fluid init - const int* gc = grid->c_grid()->global_cell; - std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); - props.reset(new IncompPropertiesFromDeck(*deck, global_cell)); + 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. diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 1ffdd2b0..479594b6 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -178,9 +178,7 @@ main(int argc, char** argv) // Grid init grid.reset(new Opm::GridManager(deck)); // Rock and fluid init - const int* gc = grid->c_grid()->global_cell; - std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); - props.reset(new Opm::BlackoilPropertiesFromDeck(deck, global_cell)); + props.reset(new Opm::BlackoilPropertiesFromDeck(deck, *grid->c_grid())); // Wells init. wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); check_well_controls = param.getDefault("check_well_controls", false); @@ -423,8 +421,10 @@ main(int argc, char** argv) std::cout << "Making " << num_transport_substeps << " transport substeps." << std::endl; } for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) { + // Note that for now we do not handle rock compressibility, + // although the transport solver should be able to. reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0], - &porevol[0], &reorder_src[0], stepsize, state.saturation()); + &porevol[0], &porevol[0], &reorder_src[0], stepsize, state.saturation()); // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); if (use_segregation_split) { THROW("Segregation not implemented yet."); diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 69697378..2c6eea43 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -309,9 +309,7 @@ main(int argc, char** argv) // Grid init grid.reset(new Opm::GridManager(deck)); // Rock and fluid init - const int* gc = grid->c_grid()->global_cell; - std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); - props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell)); + props.reset(new Opm::IncompPropertiesFromDeck(deck, *grid->c_grid())); // Wells init. wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); check_well_controls = param.getDefault("check_well_controls", false); diff --git a/examples/wells_example.cpp b/examples/wells_example.cpp index 9f7be7ec..1c523420 100644 --- a/examples/wells_example.cpp +++ b/examples/wells_example.cpp @@ -38,10 +38,8 @@ int main(int argc, char** argv) // Finally handle the wells WellsManager wells(parser, *grid.c_grid(), NULL); - std::vector global_cells(grid.c_grid()->global_cell, grid.c_grid()->global_cell + grid.c_grid()->number_of_cells); - double gravity[3] = {0.0, 0.0, parameters.getDefault("gravity", 0.0)}; - IncompPropertiesFromDeck incomp_properties(parser, global_cells); + IncompPropertiesFromDeck incomp_properties(parser, *grid.c_grid()); RockCompressibility rock_comp(parser); diff --git a/opm/core/GridManager.cpp b/opm/core/GridManager.cpp index 16140a07..f84d1c55 100644 --- a/opm/core/GridManager.cpp +++ b/opm/core/GridManager.cpp @@ -22,6 +22,8 @@ #include #include #include +#include +#include @@ -33,36 +35,25 @@ namespace Opm /// Construct a 3d corner-point grid from a deck. GridManager::GridManager(const Opm::EclipseGridParser& deck) { - // Extract data from deck. - const std::vector& zcorn = deck.getFloatingPointValue("ZCORN"); - const std::vector& coord = deck.getFloatingPointValue("COORD"); - const int* actnum = 0; - if (deck.hasField("ACTNUM")) { - actnum = &(deck.getIntegerValue("ACTNUM")[0]); + // We accept two different ways to specify the grid. + // 1. Corner point format. + // Requires ZCORN, COORDS, DIMENS or SPECGRID, optionally ACTNUM. + // 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 { + THROW("Could not initialize grid from deck. " + "Need either ZCORN + COORD or DXV + DYV + DZV keywords."); } - std::vector dims; - if (deck.hasField("DIMENS")) { - dims = deck.getIntegerValue("DIMENS"); - } else if (deck.hasField("SPECGRID")) { - dims = deck.getSPECGRID().dimensions; - } else { - THROW("Deck must have either DIMENS or SPECGRID."); - } - - // Collect in input struct for preprocessing. - struct grdecl grdecl; - grdecl.zcorn = &zcorn[0]; - grdecl.coord = &coord[0]; - grdecl.actnum = actnum; - grdecl.dims[0] = dims[0]; - grdecl.dims[1] = dims[1]; - grdecl.dims[2] = dims[2]; - - // Process grid. - ug_ = create_grid_cornerpoint(&grdecl, 0.0); - if (!ug_) { - THROW("Failed to construct grid."); - } } @@ -71,10 +62,10 @@ namespace Opm /// Construct a 2d cartesian grid with cells of unit size. GridManager::GridManager(int nx, int ny) { - ug_ = create_grid_cart2d(nx, ny); - if (!ug_) { - THROW("Failed to construct grid."); - } + ug_ = create_grid_cart2d(nx, ny); + if (!ug_) { + THROW("Failed to construct grid."); + } } @@ -83,10 +74,10 @@ namespace Opm /// Construct a 3d cartesian grid with cells of unit size. GridManager::GridManager(int nx, int ny, int nz) { - ug_ = create_grid_cart3d(nx, ny, nz); - if (!ug_) { - THROW("Failed to construct grid."); - } + ug_ = create_grid_cart3d(nx, ny, nz); + if (!ug_) { + THROW("Failed to construct grid."); + } } @@ -96,10 +87,10 @@ namespace Opm GridManager::GridManager(int nx, int ny, int nz, double dx, double dy, double dz) { - ug_ = create_grid_hexa3d(nx, ny, nz, dx, dy, dz); - if (!ug_) { - THROW("Failed to construct grid."); - } + ug_ = create_grid_hexa3d(nx, ny, nz, dx, dy, dz); + if (!ug_) { + THROW("Failed to construct grid."); + } } @@ -108,7 +99,7 @@ namespace Opm /// Destructor. GridManager::~GridManager() { - destroy_grid(ug_); + destroy_grid(ug_); } @@ -119,10 +110,99 @@ namespace Opm /// to make it clear that we are returning a C-compatible struct. const UnstructuredGrid* GridManager::c_grid() const { - return ug_; + return ug_; } + // Construct corner-point grid from deck. + void GridManager::initFromDeckCornerpoint(const Opm::EclipseGridParser& deck) + { + // Extract data from deck. + const std::vector& zcorn = deck.getFloatingPointValue("ZCORN"); + const std::vector& coord = deck.getFloatingPointValue("COORD"); + const int* actnum = 0; + if (deck.hasField("ACTNUM")) { + actnum = &(deck.getIntegerValue("ACTNUM")[0]); + } + std::vector dims; + if (deck.hasField("DIMENS")) { + dims = deck.getIntegerValue("DIMENS"); + } else if (deck.hasField("SPECGRID")) { + dims = deck.getSPECGRID().dimensions; + } else { + THROW("Deck must have either DIMENS or SPECGRID."); + } + + // Collect in input struct for preprocessing. + struct grdecl grdecl; + grdecl.zcorn = &zcorn[0]; + grdecl.coord = &coord[0]; + grdecl.actnum = actnum; + grdecl.dims[0] = dims[0]; + grdecl.dims[1] = dims[1]; + grdecl.dims[2] = dims[2]; + + // Process grid. + ug_ = create_grid_cornerpoint(&grdecl, 0.0); + if (!ug_) { + THROW("Failed to construct grid."); + } + } + + + namespace + { + std::vector coordsFromDeltas(const std::vector& deltas) + { + std::vector coords(deltas.size() + 1); + coords[0] = 0.0; + std::partial_sum(deltas.begin(), deltas.end(), coords.begin() + 1); + return coords; + } + } // anonymous namespace + + + // Construct tensor grid from deck. + void GridManager::initFromDeckTensorgrid(const Opm::EclipseGridParser& deck) + { + // Extract coordinates (or offsets from top, in case of z). + const std::vector& dxv = deck.getFloatingPointValue("DXV"); + const std::vector& dyv = deck.getFloatingPointValue("DYV"); + const std::vector& dzv = deck.getFloatingPointValue("DZV"); + std::vector x = coordsFromDeltas(dxv); + std::vector y = coordsFromDeltas(dyv); + std::vector z = coordsFromDeltas(dzv); + + // Extract top corner depths, if available. + const double* top_depths = 0; + std::vector top_depths_vec; + if (deck.hasField("DEPTHZ")) { + const std::vector& depthz = deck.getFloatingPointValue("DEPTHZ"); + if (depthz.size() != x.size()*y.size()) { + THROW("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& tops = deck.getFloatingPointValue("TOPS"); + if (std::count(tops.begin(), tops.end(), tops[0]) != int(tops.size())) { + THROW("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]; + } + + // Construct grid. + ug_ = create_grid_tensor3d(dxv.size(), dyv.size(), dzv.size(), + &x[0], &y[0], &z[0], top_depths); + if (!ug_) { + THROW("Failed to construct grid."); + } + } + + } // namespace Opm diff --git a/opm/core/GridManager.hpp b/opm/core/GridManager.hpp index 97dd84fd..26e335f2 100644 --- a/opm/core/GridManager.hpp +++ b/opm/core/GridManager.hpp @@ -33,39 +33,46 @@ namespace Opm /// encapsulates creation and destruction of the grid. /// The following grid types can be constructed: /// - 3d corner-point grids (from deck input) + /// - 3d tensor grids (from deck input) /// - 2d cartesian grids /// - 3d cartesian grids /// The resulting UnstructuredGrid is available through the c_grid() method. class GridManager { public: - /// Construct a 3d corner-point grid from a deck. - GridManager(const Opm::EclipseGridParser& deck); + /// Construct a 3d corner-point grid or tensor grid from a deck. + GridManager(const Opm::EclipseGridParser& deck); - /// Construct a 2d cartesian grid with cells of unit size. - GridManager(int nx, int ny); + /// Construct a 2d cartesian grid with cells of unit size. + GridManager(int nx, int ny); - /// Construct a 3d cartesian grid with cells of unit size. - GridManager(int nx, int ny, int nz); + /// Construct a 3d cartesian grid with cells of unit size. + GridManager(int nx, int ny, int nz); - /// Construct a 3d cartesian grid with cells of size [dx, dy, dz]. - GridManager(int nx, int ny, int nz, - double dx, double dy, double dz); + /// Construct a 3d cartesian grid with cells of size [dx, dy, dz]. + GridManager(int nx, int ny, int nz, + double dx, double dy, double dz); - /// Destructor. - ~GridManager(); + /// Destructor. + ~GridManager(); - /// 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; + /// 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; private: - // Disable copying and assignment. - GridManager(const GridManager& other); - GridManager& operator=(const GridManager& other); - // The managed UnstructuredGrid. - UnstructuredGrid* ug_; + // Disable copying and assignment. + GridManager(const GridManager& other); + GridManager& operator=(const GridManager& other); + + // Construct corner-point grid from deck. + void initFromDeckCornerpoint(const Opm::EclipseGridParser& deck); + // Construct tensor grid from deck. + void initFromDeckTensorgrid(const Opm::EclipseGridParser& deck); + + // The managed UnstructuredGrid. + UnstructuredGrid* ug_; }; } // namespace Opm diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index 45017a46..b4ea9ee3 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -86,7 +86,7 @@ namespace EclipseKeywords string("MULTPV"), string("PRESSURE"), string("SGAS"), string("SWAT"), string("SOIL"), string("RS"), string("DXV"), string("DYV"), string("DZV"), - string("DEPTHZ") + string("DEPTHZ"), string("TOPS") }; const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); @@ -525,9 +525,9 @@ void EclipseGridParser::convertToSI() // Find the right unit. double unit = 1e100; bool do_convert = true; - if (key == "COORD" || key == "ZCORN" || - key == "DXV" || key == "DYV" || key == "DZV" || - key == "DEPTHZ") { + if (key == "COORD" || key == "ZCORN" || + key == "DXV" || key == "DYV" || key == "DZV" || + key == "DEPTHZ" || key == "TOPS") { unit = units_.length; } else if (key == "PERMX" || key == "PERMY" || key == "PERMZ" || key == "PERMXX" || key == "PERMYY" || key == "PERMZZ" || diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 0eb23cb0..004e6f88 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -23,11 +23,11 @@ namespace Opm { BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const std::vector& global_cell) + const UnstructuredGrid& grid) { - rock_.init(deck, global_cell); + rock_.init(deck, grid); pvt_.init(deck); - satprops_.init(deck, global_cell); + satprops_.init(deck, grid); if (pvt_.numPhases() != satprops_.numPhases()) { THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index 9c20b568..a2a6c3de 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -27,6 +27,8 @@ #include #include +struct UnstructuredGrid; + namespace Opm { @@ -35,12 +37,13 @@ namespace Opm class BlackoilPropertiesFromDeck : public BlackoilPropertiesInterface { public: - /// Construct from deck and cell mapping. - /// \param deck eclipse input parser - /// \param global_cell mapping from cell indices (typically from a processed 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. BlackoilPropertiesFromDeck(const EclipseGridParser& deck, - const std::vector& global_cell); + const UnstructuredGrid& grid); /// Destructor. virtual ~BlackoilPropertiesFromDeck(); diff --git a/opm/core/fluid/IncompPropertiesFromDeck.cpp b/opm/core/fluid/IncompPropertiesFromDeck.cpp index e61e13c8..570154bf 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.cpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.cpp @@ -27,11 +27,11 @@ namespace Opm { IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck, - const std::vector& global_cell) + const UnstructuredGrid& grid) { - rock_.init(deck, global_cell); + rock_.init(deck, grid); pvt_.init(deck); - satprops_.init(deck, global_cell); + satprops_.init(deck, grid); if (pvt_.numPhases() != satprops_.numPhases()) { THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); diff --git a/opm/core/fluid/IncompPropertiesFromDeck.hpp b/opm/core/fluid/IncompPropertiesFromDeck.hpp index d17dd1b7..68623e5c 100644 --- a/opm/core/fluid/IncompPropertiesFromDeck.hpp +++ b/opm/core/fluid/IncompPropertiesFromDeck.hpp @@ -26,6 +26,8 @@ #include #include +struct UnstructuredGrid; + namespace Opm { @@ -43,12 +45,13 @@ namespace Opm class IncompPropertiesFromDeck : public IncompPropertiesInterface { public: - /// Construct from deck and cell mapping. - /// \param deck eclipse input parser - /// \param global_cell mapping from cell indices (typically from a processed 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(const EclipseGridParser& deck, - const std::vector& global_cell); + const UnstructuredGrid& grid); /// Destructor. virtual ~IncompPropertiesFromDeck(); diff --git a/opm/core/fluid/RockFromDeck.cpp b/opm/core/fluid/RockFromDeck.cpp index 94777904..6568e8cc 100644 --- a/opm/core/fluid/RockFromDeck.cpp +++ b/opm/core/fluid/RockFromDeck.cpp @@ -19,7 +19,7 @@ #include - +#include #include namespace Opm @@ -53,28 +53,29 @@ namespace Opm /// Initialize from deck and cell mapping. /// \param deck Deck input parser - /// \param global_cell mapping from cell indices (typically from a processed grid) + /// \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 std::vector& global_cell) + const UnstructuredGrid& grid) { - assignPorosity(deck, global_cell); - permfield_valid_.assign(global_cell.size(), false); + assignPorosity(deck, grid); + permfield_valid_.assign(grid.number_of_cells, false); const double perm_threshold = 0.0; // Maybe turn into parameter? - assignPermeability(deck, global_cell, perm_threshold); + assignPermeability(deck, grid, perm_threshold); } void RockFromDeck::assignPorosity(const EclipseGridParser& parser, - const std::vector& global_cell) + const UnstructuredGrid& grid) { - porosity_.assign(global_cell.size(), 1.0); - + porosity_.assign(grid.number_of_cells, 1.0); + const int* gc = grid.global_cell; if (parser.hasField("PORO")) { const std::vector& poro = parser.getFloatingPointValue("PORO"); - for (int c = 0; c < int(porosity_.size()); ++c) { - porosity_[c] = poro[global_cell[c]]; + const int deck_pos = (gc == NULL) ? c : gc[c]; + porosity_[c] = poro[deck_pos]; } } } @@ -82,14 +83,16 @@ namespace Opm void RockFromDeck::assignPermeability(const EclipseGridParser& parser, - const std::vector& global_cell, + const UnstructuredGrid& grid, double perm_threshold) { const int dim = 3; const int num_global_cells = numGlobalCells(parser); + const int nc = grid.number_of_cells; + ASSERT (num_global_cells > 0); - permeability_.assign(dim * dim * global_cell.size(), 0.0); + permeability_.assign(dim * dim * nc, 0.0); std::vector*> tensor; tensor.reserve(10); @@ -111,13 +114,13 @@ namespace Opm // chosen) default value... // if (tensor.size() > 1) { - const int nc = global_cell.size(); - int off = 0; + const int* gc = grid.global_cell; + int off = 0; for (int c = 0; c < nc; ++c, off += dim*dim) { // SharedPermTensor K(dim, dim, &permeability_[off]); int kix = 0; - const int glob = global_cell[c]; + const int glob = (gc == NULL) ? c : gc[c]; for (int i = 0; i < dim; ++i) { for (int j = 0; j < dim; ++j, ++kix) { diff --git a/opm/core/fluid/RockFromDeck.hpp b/opm/core/fluid/RockFromDeck.hpp index abd9c14c..65027e42 100644 --- a/opm/core/fluid/RockFromDeck.hpp +++ b/opm/core/fluid/RockFromDeck.hpp @@ -24,6 +24,7 @@ #include #include +struct UnstructuredGrid; namespace Opm { @@ -34,12 +35,13 @@ namespace Opm /// Default constructor. RockFromDeck(); - /// Initialize from deck and cell mapping. + /// Initialize from deck and grid. /// \param deck Deck input parser - /// \param global_cell mapping from cell indices (typically from a processed grid) + /// \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 std::vector& global_cell); + const UnstructuredGrid& grid); /// \return D, the number of spatial dimensions. Always 3 for deck input. int numDimensions() const @@ -69,9 +71,9 @@ namespace Opm private: void assignPorosity(const EclipseGridParser& parser, - const std::vector& global_cell); + const UnstructuredGrid& grid); void assignPermeability(const EclipseGridParser& parser, - const std::vector& global_cell, + const UnstructuredGrid& grid, const double perm_threshold); std::vector porosity_; diff --git a/opm/core/fluid/SaturationPropsFromDeck.cpp b/opm/core/fluid/SaturationPropsFromDeck.cpp index 89f81278..36cf8364 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.cpp +++ b/opm/core/fluid/SaturationPropsFromDeck.cpp @@ -18,6 +18,7 @@ */ #include +#include #include #include #include @@ -33,7 +34,7 @@ namespace Opm /// Initialize from deck. void SaturationPropsFromDeck::init(const EclipseGridParser& deck, - const std::vector& global_cell) + const UnstructuredGrid& grid) { phase_usage_ = phaseUsageFromDeck(deck); @@ -49,10 +50,12 @@ namespace Opm if (deck.hasField("SATNUM")) { const std::vector& satnum = deck.getIntegerValue("SATNUM"); satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); - int num_cells = global_cell.size(); + const int num_cells = grid.number_of_cells; cell_to_func_.resize(num_cells); + const int* gc = grid.global_cell; for (int cell = 0; cell < num_cells; ++cell) { - cell_to_func_[cell] = satnum[global_cell[cell]] - 1; + const int deck_pos = (gc == NULL) ? cell : gc[cell]; + cell_to_func_[cell] = satnum[deck_pos] - 1; } } diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index 7c8be8d4..b83a1271 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -25,6 +25,8 @@ #include #include +struct UnstructuredGrid; + namespace Opm { @@ -34,10 +36,13 @@ namespace Opm /// Default constructor. SaturationPropsFromDeck(); - /// Initialize from deck. - /// global_cell maps from grid cells to their original logical Cartesian indices. + /// 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 std::vector& global_cell); + const UnstructuredGrid& grid); /// \return P, the number of phases. int numPhases() const; diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 7d936d6e..154f0633 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -341,7 +341,7 @@ namespace Opm const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1]; props_.density(1, &cell_A_[np*np*c[j]], &gravcontrib[j][0]); for (int p = 0; p < np; ++p) { - gravcontrib[j][p] *= depth_diff; + gravcontrib[j][p] *= depth_diff*grav; } } else { std::fill(gravcontrib[j].begin(), gravcontrib[j].end(), 0.0); diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp index c04aac0a..ce8ebb97 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp @@ -39,71 +39,73 @@ namespace Opm TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid, - const Opm::BlackoilPropertiesInterface& props, - const double tol, - const int maxit) - : grid_(grid), - props_(props), - tol_(tol), - maxit_(maxit), - darcyflux_(0), - source_(0), - dt_(0.0), - saturation_(grid.number_of_cells, -1.0), - fractionalflow_(grid.number_of_cells, -1.0), - mob_(2*grid.number_of_cells, -1.0), + const Opm::BlackoilPropertiesInterface& props, + const double tol, + const int maxit) + : grid_(grid), + props_(props), + tol_(tol), + maxit_(maxit), + darcyflux_(0), + source_(0), + dt_(0.0), + saturation_(grid.number_of_cells, -1.0), + fractionalflow_(grid.number_of_cells, -1.0), + mob_(2*grid.number_of_cells, -1.0), ia_upw_(grid.number_of_cells + 1, -1), - ja_upw_(grid.number_of_faces, -1), - ia_downw_(grid.number_of_cells + 1, -1), - ja_downw_(grid.number_of_faces, -1) + ja_upw_(grid.number_of_faces, -1), + ia_downw_(grid.number_of_cells + 1, -1), + ja_downw_(grid.number_of_faces, -1) { - if (props.numPhases() != 2) { - THROW("Property object must have 2 phases"); - } + if (props.numPhases() != 2) { + THROW("Property object must have 2 phases"); + } int np = props.numPhases(); - int num_cells = props.numCells(); - visc_.resize(np*num_cells); + int num_cells = props.numCells(); + visc_.resize(np*num_cells); A_.resize(np*np*num_cells); - smin_.resize(np*num_cells); - smax_.resize(np*num_cells); - allcells_.resize(num_cells); - for (int i = 0; i < num_cells; ++i) { - allcells_[i] = i; - } - props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]); + smin_.resize(np*num_cells); + smax_.resize(np*num_cells); + allcells_.resize(num_cells); + for (int i = 0; i < num_cells; ++i) { + allcells_[i] = i; + } + props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]); } void TransportModelCompressibleTwophase::solve(const double* darcyflux, const double* pressure, const double* surfacevol0, + const double* porevolume0, const double* porevolume, const double* source, const double dt, std::vector& saturation) { - darcyflux_ = darcyflux; + darcyflux_ = darcyflux; surfacevol0_ = surfacevol0; + porevolume0_ = porevolume0; porevolume_ = porevolume; - source_ = source; - dt_ = dt; + source_ = source; + dt_ = dt; toWaterSat(saturation, saturation_); props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); - std::vector seq(grid_.number_of_cells); - std::vector comp(grid_.number_of_cells + 1); - int ncomp; - compute_sequence_graph(&grid_, darcyflux_, - &seq[0], &comp[0], &ncomp, - &ia_upw_[0], &ja_upw_[0]); - const int nf = grid_.number_of_faces; - std::vector neg_darcyflux(nf); - std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate()); - compute_sequence_graph(&grid_, &neg_darcyflux[0], - &seq[0], &comp[0], &ncomp, - &ia_downw_[0], &ja_downw_[0]); - reorderAndTransport(grid_, darcyflux); + std::vector seq(grid_.number_of_cells); + std::vector comp(grid_.number_of_cells + 1); + int ncomp; + compute_sequence_graph(&grid_, darcyflux_, + &seq[0], &comp[0], &ncomp, + &ia_upw_[0], &ja_upw_[0]); + const int nf = grid_.number_of_faces; + std::vector neg_darcyflux(nf); + std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate()); + compute_sequence_graph(&grid_, &neg_darcyflux[0], + &seq[0], &comp[0], &ncomp, + &ia_downw_[0], &ja_downw_[0]); + reorderAndTransport(grid_, darcyflux); toBothSat(saturation_, saturation); } @@ -111,190 +113,190 @@ namespace Opm // // [[ incompressible was: r(s) = s - s0 + dt/pv*( influx + outflux*f(s) ) ]] // - // r(s) = s - B*z0 + dt/pv*( influx + outflux*f(s)) + // r(s) = s - B*z0 + dt/pv*( influx + outflux*f(s) ) // // @@@ What about the source term - // + // // where influx is water influx, outflux is total outflux. - // We need the formula influx = B_i sum_{j->i} b_j v_{ij} + q_w. - // outflux = B_i sum_{i->j} b_i v_{ij} - q = sum_{i->j} v_{ij} - q (as before) + // We need the formula influx = B_i sum_{j->i} b_j v_{ij} - B_i q_w. + // outflux = B_i sum_{i->j} b_i v_{ij} - B_i q = sum_{i->j} v_{ij} - B_i q // Influxes are negative, outfluxes positive. struct TransportModelCompressibleTwophase::Residual { - int cell; - double s0; - double influx; // sum_j min(v_ij, 0)*f(s_j) + q_w // TODO: fix comment. - double outflux; // sum_j max(v_ij, 0) - q + int cell; + double B_cell; + double z0; + double influx; // B_i sum_j b_j min(v_ij, 0)*f(s_j) - B_i q_w + double outflux; // sum_j max(v_ij, 0) - B_i q // @@@ TODO: figure out change to rock-comp. terms with fluid compr. - // double comp_term; // q - sum_j v_ij - double dtpv; // dt/pv(i) - const TransportModelCompressibleTwophase& tm; - explicit Residual(const TransportModelCompressibleTwophase& tmodel, int cell_index) - : tm(tmodel) - { - cell = cell_index; - s0 = tm.saturation_[cell]; + double comp_term; // Now: used to be: q - sum_j v_ij + double dtpv; // dt/pv(i) + const TransportModelCompressibleTwophase& tm; + explicit Residual(const TransportModelCompressibleTwophase& tmodel, int cell_index) + : tm(tmodel) + { + cell = cell_index; const int np = tm.props_.numPhases(); - const double B_cell = 1.0/tm.A_[np*np*cell + 0]; + z0 = tm.surfacevol0_[np*cell + 0]; // I.e. water surface volume + B_cell = 1.0/tm.A_[np*np*cell + 0]; double src_flux = -tm.source_[cell]; bool src_is_inflow = src_flux < 0.0; - influx = src_is_inflow ? src_flux : 0.0; - outflux = !src_is_inflow ? src_flux : 0.0; - // comp_term = tm.source_[cell]; // Note: this assumes that all source flux is water. - dtpv = tm.dt_/tm.porevolume_[cell]; - for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { - const int f = tm.grid_.cell_faces[i]; - double flux; - int other; - // Compute cell flux - if (cell == tm.grid_.face_cells[2*f]) { - flux = tm.darcyflux_[f]; - other = tm.grid_.face_cells[2*f+1]; - } else { - flux =-tm.darcyflux_[f]; - other = tm.grid_.face_cells[2*f]; - } - // Add flux to influx or outflux, if interior. - if (other != -1) { - if (flux < 0.0) { + influx = src_is_inflow ? B_cell*src_flux : 0.0; + outflux = !src_is_inflow ? B_cell*src_flux : 0.0; + comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell]; + dtpv = tm.dt_/tm.porevolume0_[cell]; + for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { + const int f = tm.grid_.cell_faces[i]; + double flux; + int other; + // Compute cell flux + if (cell == tm.grid_.face_cells[2*f]) { + flux = tm.darcyflux_[f]; + other = tm.grid_.face_cells[2*f+1]; + } else { + flux =-tm.darcyflux_[f]; + other = tm.grid_.face_cells[2*f]; + } + // Add flux to influx or outflux, if interior. + if (other != -1) { + if (flux < 0.0) { const double b_face = tm.A_[np*np*other + 0]; - influx += B_cell*b_face*flux*tm.fractionalflow_[other]; - } else { - outflux += flux; - } - // comp_term -= flux; - } - } - } - double operator()(double s) const - { - // return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx + s*comp_term); - return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx); - } + influx += B_cell*b_face*flux*tm.fractionalflow_[other]; + } else { + outflux += flux; // Because B_cell*b_face = 1 for outflow faces + } + } + } + } + double operator()(double s) const + { + // return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx + s*comp_term); + return s - B_cell*z0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx) + s*comp_term; + } }; void TransportModelCompressibleTwophase::solveSingleCell(const int cell) { - Residual res(*this, cell); - int iters_used; - saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); - fractionalflow_[cell] = fracFlow(saturation_[cell], cell); + Residual res(*this, cell); + int iters_used; + saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); + fractionalflow_[cell] = fracFlow(saturation_[cell], cell); } void TransportModelCompressibleTwophase::solveMultiCell(const int num_cells, const int* cells) { - // Experiment: when a cell changes more than the tolerance, - // mark all downwind cells as needing updates. After - // computing a single update in each cell, use marks - // to guide further updating. Clear mark in cell when - // its solution gets updated. - // Verdict: this is a good one! Approx. halved total time. - std::vector needs_update(num_cells, 1); - // This one also needs the mapping from all cells to - // the strongly connected subset to filter out connections - std::vector pos(grid_.number_of_cells, -1); - for (int i = 0; i < num_cells; ++i) { - const int cell = cells[i]; - pos[cell] = i; - } + // Experiment: when a cell changes more than the tolerance, + // mark all downwind cells as needing updates. After + // computing a single update in each cell, use marks + // to guide further updating. Clear mark in cell when + // its solution gets updated. + // Verdict: this is a good one! Approx. halved total time. + std::vector needs_update(num_cells, 1); + // This one also needs the mapping from all cells to + // the strongly connected subset to filter out connections + std::vector pos(grid_.number_of_cells, -1); + for (int i = 0; i < num_cells; ++i) { + const int cell = cells[i]; + pos[cell] = i; + } - // Note: partially copied from below. - const double tol = 1e-9; - const int max_iters = 300; - // Must store s0 before we start. - std::vector s0(num_cells); - // Must set initial fractional flows before we start. - // Also, we compute the # of upstream neighbours. - // std::vector num_upstream(num_cells); - for (int i = 0; i < num_cells; ++i) { - const int cell = cells[i]; - fractionalflow_[cell] = fracFlow(saturation_[cell], cell); - s0[i] = saturation_[cell]; - // num_upstream[i] = ia_upw_[cell + 1] - ia_upw_[cell]; - } - // Solve once in each cell. - // std::vector fully_marked_stack; - // fully_marked_stack.reserve(num_cells); - int num_iters = 0; - int update_count = 0; // Change name/meaning to cells_updated? - do { - update_count = 0; // Must reset count for every iteration. - for (int i = 0; i < num_cells; ++i) { - // while (!fully_marked_stack.empty()) { - // // std::cout << "# fully marked cells = " << fully_marked_stack.size() << std::endl; - // const int fully_marked_ci = fully_marked_stack.back(); - // fully_marked_stack.pop_back(); - // ++update_count; - // const int cell = cells[fully_marked_ci]; - // const double old_s = saturation_[cell]; - // saturation_[cell] = s0[fully_marked_ci]; - // solveSingleCell(cell); - // const double s_change = std::fabs(saturation_[cell] - old_s); - // if (s_change > tol) { - // // Mark downwind cells. - // for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) { - // const int downwind_cell = ja_downw_[j]; - // int ci = pos[downwind_cell]; - // ++needs_update[ci]; - // if (needs_update[ci] == num_upstream[ci]) { - // fully_marked_stack.push_back(ci); - // } - // } - // } - // // Unmark this cell. - // needs_update[fully_marked_ci] = 0; - // } - if (!needs_update[i]) { - continue; - } - ++update_count; - const int cell = cells[i]; - const double old_s = saturation_[cell]; - saturation_[cell] = s0[i]; - solveSingleCell(cell); - const double s_change = std::fabs(saturation_[cell] - old_s); - if (s_change > tol) { - // Mark downwind cells. - for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) { - const int downwind_cell = ja_downw_[j]; - int ci = pos[downwind_cell]; + // Note: partially copied from below. + const double tol = 1e-9; + const int max_iters = 300; + // Must store s0 before we start. + std::vector s0(num_cells); + // Must set initial fractional flows before we start. + // Also, we compute the # of upstream neighbours. + // std::vector num_upstream(num_cells); + for (int i = 0; i < num_cells; ++i) { + const int cell = cells[i]; + fractionalflow_[cell] = fracFlow(saturation_[cell], cell); + s0[i] = saturation_[cell]; + // num_upstream[i] = ia_upw_[cell + 1] - ia_upw_[cell]; + } + // Solve once in each cell. + // std::vector fully_marked_stack; + // fully_marked_stack.reserve(num_cells); + int num_iters = 0; + int update_count = 0; // Change name/meaning to cells_updated? + do { + update_count = 0; // Must reset count for every iteration. + for (int i = 0; i < num_cells; ++i) { + // while (!fully_marked_stack.empty()) { + // // std::cout << "# fully marked cells = " << fully_marked_stack.size() << std::endl; + // const int fully_marked_ci = fully_marked_stack.back(); + // fully_marked_stack.pop_back(); + // ++update_count; + // const int cell = cells[fully_marked_ci]; + // const double old_s = saturation_[cell]; + // saturation_[cell] = s0[fully_marked_ci]; + // solveSingleCell(cell); + // const double s_change = std::fabs(saturation_[cell] - old_s); + // if (s_change > tol) { + // // Mark downwind cells. + // for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) { + // const int downwind_cell = ja_downw_[j]; + // int ci = pos[downwind_cell]; + // ++needs_update[ci]; + // if (needs_update[ci] == num_upstream[ci]) { + // fully_marked_stack.push_back(ci); + // } + // } + // } + // // Unmark this cell. + // needs_update[fully_marked_ci] = 0; + // } + if (!needs_update[i]) { + continue; + } + ++update_count; + const int cell = cells[i]; + const double old_s = saturation_[cell]; + saturation_[cell] = s0[i]; + solveSingleCell(cell); + const double s_change = std::fabs(saturation_[cell] - old_s); + if (s_change > tol) { + // Mark downwind cells. + for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) { + const int downwind_cell = ja_downw_[j]; + int ci = pos[downwind_cell]; if (ci != -1) { needs_update[ci] = 1; } - // ++needs_update[ci]; - // if (needs_update[ci] == num_upstream[ci]) { - // fully_marked_stack.push_back(ci); - // } - } - } - // Unmark this cell. - needs_update[i] = 0; - } - // std::cout << "Iter = " << num_iters << " update_count = " << update_count - // << " # marked cells = " - // << std::accumulate(needs_update.begin(), needs_update.end(), 0) << std::endl; - } while (update_count > 0 && ++num_iters < max_iters); + // ++needs_update[ci]; + // if (needs_update[ci] == num_upstream[ci]) { + // fully_marked_stack.push_back(ci); + // } + } + } + // Unmark this cell. + needs_update[i] = 0; + } + // std::cout << "Iter = " << num_iters << " update_count = " << update_count + // << " # marked cells = " + // << std::accumulate(needs_update.begin(), needs_update.end(), 0) << std::endl; + } while (update_count > 0 && ++num_iters < max_iters); - // Done with iterations, check if we succeeded. - if (update_count > 0) { - THROW("In solveMultiCell(), we did not converge after " - << num_iters << " iterations. Remaining update count = " << update_count); - } - std::cout << "Solved " << num_cells << " cell multicell problem in " - << num_iters << " iterations." << std::endl; + // Done with iterations, check if we succeeded. + if (update_count > 0) { + THROW("In solveMultiCell(), we did not converge after " + << num_iters << " iterations. Remaining update count = " << update_count); + } + std::cout << "Solved " << num_cells << " cell multicell problem in " + << num_iters << " iterations." << std::endl; } double TransportModelCompressibleTwophase::fracFlow(double s, int cell) const { - double sat[2] = { s, 1.0 - s }; - double mob[2]; - props_.relperm(1, sat, &cell, mob, 0); - mob[0] /= visc_[2*cell + 0]; - mob[1] /= visc_[2*cell + 1]; - return mob[0]/(mob[0] + mob[1]); + double sat[2] = { s, 1.0 - s }; + double mob[2]; + props_.relperm(1, sat, &cell, mob, 0); + mob[0] /= visc_[2*cell + 0]; + mob[1] /= visc_[2*cell + 1]; + return mob[0]/(mob[0] + mob[1]); } @@ -307,19 +309,19 @@ namespace Opm // struct TransportModelCompressibleTwophase::GravityResidual { - int cell; + int cell; int nbcell[2]; - double s0; - double dtpv; // dt/pv(i) + double s0; + double dtpv; // dt/pv(i) double gf[2]; - const TransportModelCompressibleTwophase& tm; - explicit GravityResidual(const TransportModelCompressibleTwophase& tmodel, + const TransportModelCompressibleTwophase& tm; + explicit GravityResidual(const TransportModelCompressibleTwophase& tmodel, const std::vector& cells, const int pos, const double* gravflux) // Always oriented towards next in column. Size = colsize - 1. - : tm(tmodel) - { - cell = cells[pos]; + : tm(tmodel) + { + cell = cells[pos]; nbcell[0] = -1; gf[0] = 0.0; if (pos > 0) { @@ -332,40 +334,40 @@ namespace Opm nbcell[1] = cells[pos + 1]; gf[1] = gravflux[pos]; } - s0 = tm.saturation_[cell]; - dtpv = tm.dt_/tm.porevolume_[cell]; - - } - double operator()(double s) const - { - double res = s - s0; + s0 = tm.saturation_[cell]; + dtpv = tm.dt_/tm.porevolume0_[cell]; + + } + double operator()(double s) const + { + double res = s - s0; double mobcell[2]; tm.mobility(s, cell, mobcell); for (int nb = 0; nb < 2; ++nb) { - if (nbcell[nb] != -1) { + if (nbcell[nb] != -1) { double m[2]; if (gf[nb] < 0.0) { m[0] = mobcell[0]; m[1] = tm.mob_[2*nbcell[nb] + 1]; - } else { + } else { m[0] = tm.mob_[2*nbcell[nb]]; m[1] = mobcell[1]; - } - if (m[0] + m[1] > 0.0) { + } + if (m[0] + m[1] > 0.0) { res += -dtpv*gf[nb]*m[0]*m[1]/(m[0] + m[1]); } - } + } } return res; - } + } }; void TransportModelCompressibleTwophase::mobility(double s, int cell, double* mob) const { - double sat[2] = { s, 1.0 - s }; - props_.relperm(1, sat, &cell, mob, 0); - mob[0] /= visc_[0]; - mob[1] /= visc_[1]; + double sat[2] = { s, 1.0 - s }; + props_.relperm(1, sat, &cell, mob, 0); + mob[0] /= visc_[0]; + mob[1] /= visc_[1]; } @@ -407,7 +409,7 @@ namespace Opm saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used); } saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]); - mobility(saturation_[cell], cell, &mob_[2*cell]); + mobility(saturation_[cell], cell, &mob_[2*cell]); } @@ -418,17 +420,17 @@ namespace Opm const int nc = cells.size(); std::vector col_gravflux(nc - 1); for (int ci = 0; ci < nc - 1; ++ci) { - const int cell = cells[ci]; - const int next_cell = cells[ci + 1]; - for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) { - const int face = grid_.cell_faces[j]; - const int c1 = grid_.face_cells[2*face + 0]; + const int cell = cells[ci]; + const int next_cell = cells[ci + 1]; + for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) { + const int face = grid_.cell_faces[j]; + const int c1 = grid_.face_cells[2*face + 0]; const int c2 = grid_.face_cells[2*face + 1]; - if (c1 == next_cell || c2 == next_cell) { + if (c1 == next_cell || c2 == next_cell) { const double gf = gravflux_[face]; col_gravflux[ci] = (c1 == cell) ? gf : -gf; - } - } + } + } } // Store initial saturation s0 @@ -438,7 +440,7 @@ namespace Opm } // Solve single cell problems, repeating if necessary. - double max_s_change = 0.0; + double max_s_change = 0.0; int num_iters = 0; do { max_s_change = 0.0; @@ -454,12 +456,12 @@ namespace Opm std::fabs(saturation_[cells[ci2]] - old_s[1]))); } // std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change << std::endl; - } while (max_s_change > tol_ && ++num_iters < maxit_); + } while (max_s_change > tol_ && ++num_iters < maxit_); - if (max_s_change > tol_) { - THROW("In solveGravityColumn(), we did not converge after " - << num_iters << " iterations. Delta s = " << max_s_change); - } + if (max_s_change > tol_) { + THROW("In solveGravityColumn(), we did not converge after " + << num_iters << " iterations. Delta s = " << max_s_change); + } return num_iters + 1; } @@ -467,6 +469,7 @@ namespace Opm void TransportModelCompressibleTwophase::solveGravity(const std::vector >& columns, const double* pressure, + const double* porevolume0, const double* porevolume, const double dt, std::vector& saturation) @@ -489,6 +492,7 @@ namespace Opm } // Set up other variables. + porevolume0_ = porevolume0; porevolume_ = porevolume; dt_ = dt; toWaterSat(saturation, saturation_); diff --git a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp index d41deef6..fda83fbc 100644 --- a/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp +++ b/opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp @@ -47,6 +47,7 @@ namespace Opm /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] pressure Array of cell pressures /// \param[in] surfacevol0 Array of surface volumes at start of timestep + /// \param[in] porevolume0 Array of pore volumes at start of timestep. /// \param[in] porevolume Array of pore volumes. /// \param[in] source Transport source term. /// \param[in] dt Time step. @@ -54,6 +55,7 @@ namespace Opm void solve(const double* darcyflux, const double* pressure, const double* surfacevol0, + const double* porevolume0, const double* porevolume, const double* source, const double dt, @@ -71,6 +73,7 @@ namespace Opm /// \TODO: Implement this. void solveGravity(const std::vector >& columns, const double* pressure, + const double* porevolume0, const double* porevolume, const double dt, std::vector& saturation); @@ -96,6 +99,7 @@ namespace Opm const double* darcyflux_; // one flux per grid face const double* surfacevol0_; // one per phase per cell + const double* porevolume0_; // one volume per cell const double* porevolume_; // one volume per cell const double* source_; // one source per cell double dt_; diff --git a/tests/bo_resprop_test.cpp b/tests/bo_resprop_test.cpp index 316bd9f3..5bf3aea1 100644 --- a/tests/bo_resprop_test.cpp +++ b/tests/bo_resprop_test.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -39,14 +40,10 @@ int main(int argc, char** argv) // Parser. std::string ecl_file = param.get("filename"); Opm::EclipseGridParser deck(ecl_file); - Opm::EclipseGridInspector insp(deck); - std::tr1::array gs = insp.gridSize(); - int num_cells = gs[0]*gs[1]*gs[2]; - std::vector global_cell(num_cells); - for (int i = 0; i < num_cells; ++i) { - global_cell[i] = i; - } - Opm::BlackoilPropertiesFromDeck props(deck, global_cell); + UnstructuredGrid grid; + grid.number_of_cells = 1; + grid.global_cell = NULL; + Opm::BlackoilPropertiesFromDeck props(deck, grid); const int n = 1; double p[n] = { 150e5 };