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

Merge pull request #10 from atgeirr/master

Browse Source 
Improvements in grid creation and compressible solvers.
This commit is contained in:
Bård Skaflestad 2012-08-10 01:34:44 -07:00
commit 3d45ca3512
19 changed files with 472 additions and 367 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
examples/sim_2p_incomp_reorder.cpp
View File

@ -81,9 +81,7 @@ main(int argc, char** argv)
// Grid init // Grid init
grid.reset(new GridManager(*deck)); grid.reset(new GridManager(*deck));
// Rock and fluid init // Rock and fluid init
const int* gc = grid->c_grid()->global_cell; props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
props.reset(new IncompPropertiesFromDeck(*deck, global_cell));
// check_well_controls = param.getDefault("check_well_controls", false); // check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility. // Rock compressibility.

8
examples/sim_wateroil.cpp
View File

@ -178,9 +178,7 @@ main(int argc, char** argv)
// Grid init // Grid init
grid.reset(new Opm::GridManager(deck)); grid.reset(new Opm::GridManager(deck));
// Rock and fluid init // Rock and fluid init
const int* gc = grid->c_grid()->global_cell; props.reset(new Opm::BlackoilPropertiesFromDeck(deck, *grid->c_grid()));
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
props.reset(new Opm::BlackoilPropertiesFromDeck(deck, global_cell));
// Wells init. // Wells init.
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
check_well_controls = param.getDefault("check_well_controls", false); 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; std::cout << "Making " << num_transport_substeps << " transport substeps." << std::endl;
} }
for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) { 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], 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); // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
if (use_segregation_split) { if (use_segregation_split) {
THROW("Segregation not implemented yet."); THROW("Segregation not implemented yet.");

4
examples/spu_2p.cpp
View File

@ -309,9 +309,7 @@ main(int argc, char** argv)
// Grid init // Grid init
grid.reset(new Opm::GridManager(deck)); grid.reset(new Opm::GridManager(deck));
// Rock and fluid init // Rock and fluid init
const int* gc = grid->c_grid()->global_cell; props.reset(new Opm::IncompPropertiesFromDeck(deck, *grid->c_grid()));
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
// Wells init. // Wells init.
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
check_well_controls = param.getDefault("check_well_controls", false); check_well_controls = param.getDefault("check_well_controls", false);

4
examples/wells_example.cpp
View File

@ -38,10 +38,8 @@ int main(int argc, char** argv)
// Finally handle the wells // Finally handle the wells
WellsManager wells(parser, *grid.c_grid(), NULL); WellsManager wells(parser, *grid.c_grid(), NULL);
std::vector<int> 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<double>("gravity", 0.0)}; double gravity[3] = {0.0, 0.0, parameters.getDefault<double>("gravity", 0.0)};
IncompPropertiesFromDeck incomp_properties(parser, global_cells); IncompPropertiesFromDeck incomp_properties(parser, *grid.c_grid());
RockCompressibility rock_comp(parser); RockCompressibility rock_comp(parser);

166
opm/core/GridManager.cpp
View File

@ -22,6 +22,8 @@
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h> #include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/cornerpoint_grid.h> #include <opm/core/grid/cornerpoint_grid.h>
#include <algorithm>
#include <numeric>
@ -33,36 +35,25 @@ namespace Opm
/// Construct a 3d corner-point grid from a deck. /// Construct a 3d corner-point grid from a deck.
GridManager::GridManager(const Opm::EclipseGridParser& deck) GridManager::GridManager(const Opm::EclipseGridParser& deck)
{ {
// Extract data from deck. // We accept two different ways to specify the grid.
const std::vector<double>& zcorn = deck.getFloatingPointValue("ZCORN"); // 1. Corner point format.
const std::vector<double>& coord = deck.getFloatingPointValue("COORD"); // Requires ZCORN, COORDS, DIMENS or SPECGRID, optionally ACTNUM.
const int* actnum = 0; // For this format, we will verify that DXV, DYV, DZV,
if (deck.hasField("ACTNUM")) { // DEPTHZ and TOPS are not present.
actnum = &(deck.getIntegerValue("ACTNUM")[0]); // 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<int> 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. /// Construct a 2d cartesian grid with cells of unit size.
GridManager::GridManager(int nx, int ny) GridManager::GridManager(int nx, int ny)
{ {
ug_ = create_grid_cart2d(nx, ny); ug_ = create_grid_cart2d(nx, ny);
if (!ug_) { if (!ug_) {
THROW("Failed to construct grid."); THROW("Failed to construct grid.");
} }
} }
@ -83,10 +74,10 @@ namespace Opm
/// Construct a 3d cartesian grid with cells of unit size. /// Construct a 3d cartesian grid with cells of unit size.
GridManager::GridManager(int nx, int ny, int nz) GridManager::GridManager(int nx, int ny, int nz)
{ {
ug_ = create_grid_cart3d(nx, ny, nz); ug_ = create_grid_cart3d(nx, ny, nz);
if (!ug_) { if (!ug_) {
THROW("Failed to construct grid."); THROW("Failed to construct grid.");
} }
} }
@ -96,10 +87,10 @@ namespace Opm
GridManager::GridManager(int nx, int ny, int nz, GridManager::GridManager(int nx, int ny, int nz,
double dx, double dy, double dz) double dx, double dy, double dz)
{ {
ug_ = create_grid_hexa3d(nx, ny, nz, dx, dy, dz); ug_ = create_grid_hexa3d(nx, ny, nz, dx, dy, dz);
if (!ug_) { if (!ug_) {
THROW("Failed to construct grid."); THROW("Failed to construct grid.");
} }
} }
@ -108,7 +99,7 @@ namespace Opm
/// Destructor. /// Destructor.
GridManager::~GridManager() 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. /// to make it clear that we are returning a C-compatible struct.
const UnstructuredGrid* GridManager::c_grid() const 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<double>& zcorn = deck.getFloatingPointValue("ZCORN");
const std::vector<double>& coord = deck.getFloatingPointValue("COORD");
const int* actnum = 0;
if (deck.hasField("ACTNUM")) {
actnum = &(deck.getIntegerValue("ACTNUM")[0]);
}
std::vector<int> 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<double> coordsFromDeltas(const std::vector<double>& deltas)
{
std::vector<double> 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<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);
// 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()) {
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<double>& 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 } // namespace Opm

47
opm/core/GridManager.hpp
View File

@ -33,39 +33,46 @@ namespace Opm
/// encapsulates creation and destruction of the grid. /// encapsulates creation and destruction of the grid.
/// The following grid types can be constructed: /// The following grid types can be constructed:
/// - 3d corner-point grids (from deck input) /// - 3d corner-point grids (from deck input)
/// - 3d tensor grids (from deck input)
/// - 2d cartesian grids /// - 2d cartesian grids
/// - 3d cartesian grids /// - 3d cartesian grids
/// The resulting UnstructuredGrid is available through the c_grid() method. /// The resulting UnstructuredGrid is available through the c_grid() method.
class GridManager class GridManager
{ {
public: public:
/// Construct a 3d corner-point grid from a deck. /// Construct a 3d corner-point grid or tensor grid from a deck.
GridManager(const Opm::EclipseGridParser& deck); GridManager(const Opm::EclipseGridParser& deck);
/// Construct a 2d cartesian grid with cells of unit size. /// Construct a 2d cartesian grid with cells of unit size.
GridManager(int nx, int ny); GridManager(int nx, int ny);
/// Construct a 3d cartesian grid with cells of unit size. /// Construct a 3d cartesian grid with cells of unit size.
GridManager(int nx, int ny, int nz); GridManager(int nx, int ny, int nz);
/// Construct a 3d cartesian grid with cells of size [dx, dy, dz]. /// Construct a 3d cartesian grid with cells of size [dx, dy, dz].
GridManager(int nx, int ny, int nz, GridManager(int nx, int ny, int nz,
double dx, double dy, double dz); double dx, double dy, double dz);
/// Destructor. /// Destructor.
~GridManager(); ~GridManager();
/// Access the managed UnstructuredGrid. /// Access the managed UnstructuredGrid.
/// The method is named similarly to c_str() in std::string, /// The method is named similarly to c_str() in std::string,
/// to make it clear that we are returning a C-compatible struct. /// to make it clear that we are returning a C-compatible struct.
const UnstructuredGrid* c_grid() const; const UnstructuredGrid* c_grid() const;
private: private:
// Disable copying and assignment. // Disable copying and assignment.
GridManager(const GridManager& other); GridManager(const GridManager& other);
GridManager& operator=(const GridManager& other); GridManager& operator=(const GridManager& other);
// The managed UnstructuredGrid.
UnstructuredGrid* ug_; // 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 } // namespace Opm

8
opm/core/eclipse/EclipseGridParser.cpp
View File

@ -86,7 +86,7 @@ namespace EclipseKeywords
string("MULTPV"), string("PRESSURE"), string("SGAS"), string("MULTPV"), string("PRESSURE"), string("SGAS"),
string("SWAT"), string("SOIL"), string("RS"), string("SWAT"), string("SOIL"), string("RS"),
string("DXV"), string("DYV"), string("DZV"), string("DXV"), string("DYV"), string("DZV"),
string("DEPTHZ") string("DEPTHZ"), string("TOPS")
}; };
const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]);
@ -525,9 +525,9 @@ void EclipseGridParser::convertToSI()
// Find the right unit. // Find the right unit.
double unit = 1e100; double unit = 1e100;
bool do_convert = true; bool do_convert = true;
if (key == "COORD" || key == "ZCORN" || if (key == "COORD" || key == "ZCORN" ||
key == "DXV" || key == "DYV" || key == "DZV" || key == "DXV" || key == "DYV" || key == "DZV" ||
key == "DEPTHZ") { key == "DEPTHZ" || key == "TOPS") {
unit = units_.length; unit = units_.length;
} else if (key == "PERMX" || key == "PERMY" || key == "PERMZ" || } else if (key == "PERMX" || key == "PERMY" || key == "PERMZ" ||
key == "PERMXX" || key == "PERMYY" || key == "PERMZZ" || key == "PERMXX" || key == "PERMYY" || key == "PERMZZ" ||

6
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@ -23,11 +23,11 @@ namespace Opm
{ {
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const std::vector<int>& global_cell) const UnstructuredGrid& grid)
{ {
rock_.init(deck, global_cell); rock_.init(deck, grid);
pvt_.init(deck); pvt_.init(deck);
satprops_.init(deck, global_cell); satprops_.init(deck, grid);
if (pvt_.numPhases() != satprops_.numPhases()) { if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data (" THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");

11
opm/core/fluid/BlackoilPropertiesFromDeck.hpp
View File

@ -27,6 +27,8 @@
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
struct UnstructuredGrid;
namespace Opm namespace Opm
{ {
@ -35,12 +37,13 @@ namespace Opm
class BlackoilPropertiesFromDeck : public BlackoilPropertiesInterface class BlackoilPropertiesFromDeck : public BlackoilPropertiesInterface
{ {
public: public:
/// Construct from deck and cell mapping. /// Initialize from deck and grid.
/// \param deck eclipse input parser /// \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. /// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck, BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const std::vector<int>& global_cell); const UnstructuredGrid& grid);
/// Destructor. /// Destructor.
virtual ~BlackoilPropertiesFromDeck(); virtual ~BlackoilPropertiesFromDeck();

6
opm/core/fluid/IncompPropertiesFromDeck.cpp
View File

@ -27,11 +27,11 @@ namespace Opm
{ {
IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck, IncompPropertiesFromDeck::IncompPropertiesFromDeck(const EclipseGridParser& deck,
const std::vector<int>& global_cell) const UnstructuredGrid& grid)
{ {
rock_.init(deck, global_cell); rock_.init(deck, grid);
pvt_.init(deck); pvt_.init(deck);
satprops_.init(deck, global_cell); satprops_.init(deck, grid);
if (pvt_.numPhases() != satprops_.numPhases()) { if (pvt_.numPhases() != satprops_.numPhases()) {
THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data (" THROW("IncompPropertiesFromDeck::IncompPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ")."); << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");

11
opm/core/fluid/IncompPropertiesFromDeck.hpp
View File

@ -26,6 +26,8 @@
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
struct UnstructuredGrid;
namespace Opm namespace Opm
{ {
@ -43,12 +45,13 @@ namespace Opm
class IncompPropertiesFromDeck : public IncompPropertiesInterface class IncompPropertiesFromDeck : public IncompPropertiesInterface
{ {
public: public:
/// Construct from deck and cell mapping. /// Initialize from deck and grid.
/// \param deck eclipse input parser /// \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. /// to logical cartesian indices consistent with the deck.
IncompPropertiesFromDeck(const EclipseGridParser& deck, IncompPropertiesFromDeck(const EclipseGridParser& deck,
const std::vector<int>& global_cell); const UnstructuredGrid& grid);
/// Destructor. /// Destructor.
virtual ~IncompPropertiesFromDeck(); virtual ~IncompPropertiesFromDeck();

35
opm/core/fluid/RockFromDeck.cpp
View File

@ -19,7 +19,7 @@
#include <opm/core/fluid/RockFromDeck.hpp> #include <opm/core/fluid/RockFromDeck.hpp>
#include <opm/core/grid.h>
#include <tr1/array> #include <tr1/array>
namespace Opm namespace Opm
@ -53,28 +53,29 @@ namespace Opm
/// Initialize from deck and cell mapping. /// Initialize from deck and cell mapping.
/// \param deck Deck input parser /// \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. /// to logical cartesian indices consistent with the deck.
void RockFromDeck::init(const EclipseGridParser& deck, void RockFromDeck::init(const EclipseGridParser& deck,
const std::vector<int>& global_cell) const UnstructuredGrid& grid)
{ {
assignPorosity(deck, global_cell); assignPorosity(deck, grid);
permfield_valid_.assign(global_cell.size(), false); permfield_valid_.assign(grid.number_of_cells, false);
const double perm_threshold = 0.0; // Maybe turn into parameter? 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, void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
const std::vector<int>& 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")) { if (parser.hasField("PORO")) {
const std::vector<double>& poro = parser.getFloatingPointValue("PORO"); const std::vector<double>& poro = parser.getFloatingPointValue("PORO");
for (int c = 0; c < int(porosity_.size()); ++c) { 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, void RockFromDeck::assignPermeability(const EclipseGridParser& parser,
const std::vector<int>& global_cell, const UnstructuredGrid& grid,
double perm_threshold) double perm_threshold)
{ {
const int dim = 3; const int dim = 3;
const int num_global_cells = numGlobalCells(parser); const int num_global_cells = numGlobalCells(parser);
const int nc = grid.number_of_cells;
ASSERT (num_global_cells > 0); ASSERT (num_global_cells > 0);
permeability_.assign(dim * dim * global_cell.size(), 0.0); permeability_.assign(dim * dim * nc, 0.0);
std::vector<const std::vector<double>*> tensor; std::vector<const std::vector<double>*> tensor;
tensor.reserve(10); tensor.reserve(10);
@ -111,13 +114,13 @@ namespace Opm
// chosen) default value... // chosen) default value...
// //
if (tensor.size() > 1) { if (tensor.size() > 1) {
const int nc = global_cell.size(); const int* gc = grid.global_cell;
int off = 0; int off = 0;
for (int c = 0; c < nc; ++c, off += dim*dim) { for (int c = 0; c < nc; ++c, off += dim*dim) {
// SharedPermTensor K(dim, dim, &permeability_[off]); // SharedPermTensor K(dim, dim, &permeability_[off]);
int kix = 0; 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 i = 0; i < dim; ++i) {
for (int j = 0; j < dim; ++j, ++kix) { for (int j = 0; j < dim; ++j, ++kix) {

12
opm/core/fluid/RockFromDeck.hpp
View File

@ -24,6 +24,7 @@
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
#include <vector> #include <vector>
struct UnstructuredGrid;
namespace Opm namespace Opm
{ {
@ -34,12 +35,13 @@ namespace Opm
/// Default constructor. /// Default constructor.
RockFromDeck(); RockFromDeck();
/// Initialize from deck and cell mapping. /// Initialize from deck and grid.
/// \param deck Deck input parser /// \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. /// to logical cartesian indices consistent with the deck.
void init(const EclipseGridParser& deck, void init(const EclipseGridParser& deck,
const std::vector<int>& global_cell); const UnstructuredGrid& grid);
/// \return D, the number of spatial dimensions. Always 3 for deck input. /// \return D, the number of spatial dimensions. Always 3 for deck input.
int numDimensions() const int numDimensions() const
@ -69,9 +71,9 @@ namespace Opm
private: private:
void assignPorosity(const EclipseGridParser& parser, void assignPorosity(const EclipseGridParser& parser,
const std::vector<int>& global_cell); const UnstructuredGrid& grid);
void assignPermeability(const EclipseGridParser& parser, void assignPermeability(const EclipseGridParser& parser,
const std::vector<int>& global_cell, const UnstructuredGrid& grid,
const double perm_threshold); const double perm_threshold);
std::vector<double> porosity_; std::vector<double> porosity_;

9
opm/core/fluid/SaturationPropsFromDeck.cpp
View File

@ -18,6 +18,7 @@
*/ */
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp> #include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp> #include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
@ -33,7 +34,7 @@ namespace Opm
/// Initialize from deck. /// Initialize from deck.
void SaturationPropsFromDeck::init(const EclipseGridParser& deck, void SaturationPropsFromDeck::init(const EclipseGridParser& deck,
const std::vector<int>& global_cell) const UnstructuredGrid& grid)
{ {
phase_usage_ = phaseUsageFromDeck(deck); phase_usage_ = phaseUsageFromDeck(deck);
@ -49,10 +50,12 @@ namespace Opm
if (deck.hasField("SATNUM")) { if (deck.hasField("SATNUM")) {
const std::vector<int>& satnum = deck.getIntegerValue("SATNUM"); const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); 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); cell_to_func_.resize(num_cells);
const int* gc = grid.global_cell;
for (int cell = 0; cell < num_cells; ++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;
} }
} }

11
opm/core/fluid/SaturationPropsFromDeck.hpp
View File

@ -25,6 +25,8 @@
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp> #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <vector> #include <vector>
struct UnstructuredGrid;
namespace Opm namespace Opm
{ {
@ -34,10 +36,13 @@ namespace Opm
/// Default constructor. /// Default constructor.
SaturationPropsFromDeck(); SaturationPropsFromDeck();
/// Initialize from deck. /// Initialize from deck and grid.
/// global_cell maps from grid cells to their original logical Cartesian indices. /// \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, void init(const EclipseGridParser& deck,
const std::vector<int>& global_cell); const UnstructuredGrid& grid);
/// \return P, the number of phases. /// \return P, the number of phases.
int numPhases() const; int numPhases() const;

2
opm/core/pressure/CompressibleTpfa.cpp
View File

@ -341,7 +341,7 @@ namespace Opm
const double depth_diff = face_depth - grid_.cell_centroids[c[j]*dim + dim - 1]; 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]); props_.density(1, &cell_A_[np*np*c[j]], &gravcontrib[j][0]);
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
gravcontrib[j][p] *= depth_diff; gravcontrib[j][p] *= depth_diff*grav;
} }
} else { } else {
std::fill(gravcontrib[j].begin(), gravcontrib[j].end(), 0.0); std::fill(gravcontrib[j].begin(), gravcontrib[j].end(), 0.0);

478
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
View File

@ -39,71 +39,73 @@ namespace Opm
TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid, TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid,
const Opm::BlackoilPropertiesInterface& props, const Opm::BlackoilPropertiesInterface& props,
const double tol, const double tol,
const int maxit) const int maxit)
: grid_(grid), : grid_(grid),
props_(props), props_(props),
tol_(tol), tol_(tol),
maxit_(maxit), maxit_(maxit),
darcyflux_(0), darcyflux_(0),
source_(0), source_(0),
dt_(0.0), dt_(0.0),
saturation_(grid.number_of_cells, -1.0), saturation_(grid.number_of_cells, -1.0),
fractionalflow_(grid.number_of_cells, -1.0), fractionalflow_(grid.number_of_cells, -1.0),
mob_(2*grid.number_of_cells, -1.0), mob_(2*grid.number_of_cells, -1.0),
ia_upw_(grid.number_of_cells + 1, -1), ia_upw_(grid.number_of_cells + 1, -1),
ja_upw_(grid.number_of_faces, -1), ja_upw_(grid.number_of_faces, -1),
ia_downw_(grid.number_of_cells + 1, -1), ia_downw_(grid.number_of_cells + 1, -1),
ja_downw_(grid.number_of_faces, -1) ja_downw_(grid.number_of_faces, -1)
{ {
if (props.numPhases() != 2) { if (props.numPhases() != 2) {
THROW("Property object must have 2 phases"); THROW("Property object must have 2 phases");
} }
int np = props.numPhases(); int np = props.numPhases();
int num_cells = props.numCells(); int num_cells = props.numCells();
visc_.resize(np*num_cells); visc_.resize(np*num_cells);
A_.resize(np*np*num_cells); A_.resize(np*np*num_cells);
smin_.resize(np*num_cells); smin_.resize(np*num_cells);
smax_.resize(np*num_cells); smax_.resize(np*num_cells);
allcells_.resize(num_cells); allcells_.resize(num_cells);
for (int i = 0; i < num_cells; ++i) { for (int i = 0; i < num_cells; ++i) {
allcells_[i] = i; allcells_[i] = i;
} }
props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]); props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]);
} }
void TransportModelCompressibleTwophase::solve(const double* darcyflux, void TransportModelCompressibleTwophase::solve(const double* darcyflux,
const double* pressure, const double* pressure,
const double* surfacevol0, const double* surfacevol0,
const double* porevolume0,
const double* porevolume, const double* porevolume,
const double* source, const double* source,
const double dt, const double dt,
std::vector<double>& saturation) std::vector<double>& saturation)
{ {
darcyflux_ = darcyflux; darcyflux_ = darcyflux;
surfacevol0_ = surfacevol0; surfacevol0_ = surfacevol0;
porevolume0_ = porevolume0;
porevolume_ = porevolume; porevolume_ = porevolume;
source_ = source; source_ = source;
dt_ = dt; dt_ = dt;
toWaterSat(saturation, saturation_); toWaterSat(saturation, saturation_);
props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL); props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL); props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL);
std::vector<int> seq(grid_.number_of_cells); std::vector<int> seq(grid_.number_of_cells);
std::vector<int> comp(grid_.number_of_cells + 1); std::vector<int> comp(grid_.number_of_cells + 1);
int ncomp; int ncomp;
compute_sequence_graph(&grid_, darcyflux_, compute_sequence_graph(&grid_, darcyflux_,
&seq[0], &comp[0], &ncomp, &seq[0], &comp[0], &ncomp,
&ia_upw_[0], &ja_upw_[0]); &ia_upw_[0], &ja_upw_[0]);
const int nf = grid_.number_of_faces; const int nf = grid_.number_of_faces;
std::vector<double> neg_darcyflux(nf); std::vector<double> neg_darcyflux(nf);
std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate<double>()); std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate<double>());
compute_sequence_graph(&grid_, &neg_darcyflux[0], compute_sequence_graph(&grid_, &neg_darcyflux[0],
&seq[0], &comp[0], &ncomp, &seq[0], &comp[0], &ncomp,
&ia_downw_[0], &ja_downw_[0]); &ia_downw_[0], &ja_downw_[0]);
reorderAndTransport(grid_, darcyflux); reorderAndTransport(grid_, darcyflux);
toBothSat(saturation_, saturation); toBothSat(saturation_, saturation);
} }
@ -111,190 +113,190 @@ namespace Opm
// //
// [[ incompressible was: r(s) = s - s0 + dt/pv*( influx + outflux*f(s) ) ]] // [[ 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 // @@@ What about the source term
// //
// where influx is water influx, outflux is total outflux. // 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. // 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} - q = sum_{i->j} v_{ij} - q (as before) // 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. // Influxes are negative, outfluxes positive.
struct TransportModelCompressibleTwophase::Residual struct TransportModelCompressibleTwophase::Residual
{ {
int cell; int cell;
double s0; double B_cell;
double influx; // sum_j min(v_ij, 0)*f(s_j) + q_w // TODO: fix comment. double z0;
double outflux; // sum_j max(v_ij, 0) - q 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. // @@@ TODO: figure out change to rock-comp. terms with fluid compr.
// double comp_term; // q - sum_j v_ij double comp_term; // Now: used to be: q - sum_j v_ij
double dtpv; // dt/pv(i) double dtpv; // dt/pv(i)
const TransportModelCompressibleTwophase& tm; const TransportModelCompressibleTwophase& tm;
explicit Residual(const TransportModelCompressibleTwophase& tmodel, int cell_index) explicit Residual(const TransportModelCompressibleTwophase& tmodel, int cell_index)
: tm(tmodel) : tm(tmodel)
{ {
cell = cell_index; cell = cell_index;
s0 = tm.saturation_[cell];
const int np = tm.props_.numPhases(); 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]; double src_flux = -tm.source_[cell];
bool src_is_inflow = src_flux < 0.0; bool src_is_inflow = src_flux < 0.0;
influx = src_is_inflow ? src_flux : 0.0; influx = src_is_inflow ? B_cell*src_flux : 0.0;
outflux = !src_is_inflow ? src_flux : 0.0; outflux = !src_is_inflow ? B_cell*src_flux : 0.0;
// comp_term = tm.source_[cell]; // Note: this assumes that all source flux is water. comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell];
dtpv = tm.dt_/tm.porevolume_[cell]; dtpv = tm.dt_/tm.porevolume0_[cell];
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
const int f = tm.grid_.cell_faces[i]; const int f = tm.grid_.cell_faces[i];
double flux; double flux;
int other; int other;
// Compute cell flux // Compute cell flux
if (cell == tm.grid_.face_cells[2*f]) { if (cell == tm.grid_.face_cells[2*f]) {
flux = tm.darcyflux_[f]; flux = tm.darcyflux_[f];
other = tm.grid_.face_cells[2*f+1]; other = tm.grid_.face_cells[2*f+1];
} else { } else {
flux =-tm.darcyflux_[f]; flux =-tm.darcyflux_[f];
other = tm.grid_.face_cells[2*f]; other = tm.grid_.face_cells[2*f];
} }
// Add flux to influx or outflux, if interior. // Add flux to influx or outflux, if interior.
if (other != -1) { if (other != -1) {
if (flux < 0.0) { if (flux < 0.0) {
const double b_face = tm.A_[np*np*other + 0]; const double b_face = tm.A_[np*np*other + 0];
influx += B_cell*b_face*flux*tm.fractionalflow_[other]; influx += B_cell*b_face*flux*tm.fractionalflow_[other];
} else { } else {
outflux += flux; outflux += flux; // Because B_cell*b_face = 1 for outflow faces
} }
// comp_term -= flux; }
} }
} }
} double operator()(double s) const
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 + s*comp_term); return s - B_cell*z0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx) + s*comp_term;
return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx); }
}
}; };
void TransportModelCompressibleTwophase::solveSingleCell(const int cell) void TransportModelCompressibleTwophase::solveSingleCell(const int cell)
{ {
Residual res(*this, cell); Residual res(*this, cell);
int iters_used; int iters_used;
saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used); saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used);
fractionalflow_[cell] = fracFlow(saturation_[cell], cell); fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
} }
void TransportModelCompressibleTwophase::solveMultiCell(const int num_cells, const int* cells) void TransportModelCompressibleTwophase::solveMultiCell(const int num_cells, const int* cells)
{ {
// Experiment: when a cell changes more than the tolerance, // Experiment: when a cell changes more than the tolerance,
// mark all downwind cells as needing updates. After // mark all downwind cells as needing updates. After
// computing a single update in each cell, use marks // computing a single update in each cell, use marks
// to guide further updating. Clear mark in cell when // to guide further updating. Clear mark in cell when
// its solution gets updated. // its solution gets updated.
// Verdict: this is a good one! Approx. halved total time. // Verdict: this is a good one! Approx. halved total time.
std::vector<int> needs_update(num_cells, 1); std::vector<int> needs_update(num_cells, 1);
// This one also needs the mapping from all cells to // This one also needs the mapping from all cells to
// the strongly connected subset to filter out connections // the strongly connected subset to filter out connections
std::vector<int> pos(grid_.number_of_cells, -1); std::vector<int> pos(grid_.number_of_cells, -1);
for (int i = 0; i < num_cells; ++i) { for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i]; const int cell = cells[i];
pos[cell] = i; pos[cell] = i;
} }
// Note: partially copied from below. // Note: partially copied from below.
const double tol = 1e-9; const double tol = 1e-9;
const int max_iters = 300; const int max_iters = 300;
// Must store s0 before we start. // Must store s0 before we start.
std::vector<double> s0(num_cells); std::vector<double> s0(num_cells);
// Must set initial fractional flows before we start. // Must set initial fractional flows before we start.
// Also, we compute the # of upstream neighbours. // Also, we compute the # of upstream neighbours.
// std::vector<int> num_upstream(num_cells); // std::vector<int> num_upstream(num_cells);
for (int i = 0; i < num_cells; ++i) { for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i]; const int cell = cells[i];
fractionalflow_[cell] = fracFlow(saturation_[cell], cell); fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
s0[i] = saturation_[cell]; s0[i] = saturation_[cell];
// num_upstream[i] = ia_upw_[cell + 1] - ia_upw_[cell]; // num_upstream[i] = ia_upw_[cell + 1] - ia_upw_[cell];
} }
// Solve once in each cell. // Solve once in each cell.
// std::vector<int> fully_marked_stack; // std::vector<int> fully_marked_stack;
// fully_marked_stack.reserve(num_cells); // fully_marked_stack.reserve(num_cells);
int num_iters = 0; int num_iters = 0;
int update_count = 0; // Change name/meaning to cells_updated? int update_count = 0; // Change name/meaning to cells_updated?
do { do {
update_count = 0; // Must reset count for every iteration. update_count = 0; // Must reset count for every iteration.
for (int i = 0; i < num_cells; ++i) { for (int i = 0; i < num_cells; ++i) {
// while (!fully_marked_stack.empty()) { // while (!fully_marked_stack.empty()) {
// // std::cout << "# fully marked cells = " << fully_marked_stack.size() << std::endl; // // std::cout << "# fully marked cells = " << fully_marked_stack.size() << std::endl;
// const int fully_marked_ci = fully_marked_stack.back(); // const int fully_marked_ci = fully_marked_stack.back();
// fully_marked_stack.pop_back(); // fully_marked_stack.pop_back();
// ++update_count; // ++update_count;
// const int cell = cells[fully_marked_ci]; // const int cell = cells[fully_marked_ci];
// const double old_s = saturation_[cell]; // const double old_s = saturation_[cell];
// saturation_[cell] = s0[fully_marked_ci]; // saturation_[cell] = s0[fully_marked_ci];
// solveSingleCell(cell); // solveSingleCell(cell);
// const double s_change = std::fabs(saturation_[cell] - old_s); // const double s_change = std::fabs(saturation_[cell] - old_s);
// if (s_change > tol) { // if (s_change > tol) {
// // Mark downwind cells. // // Mark downwind cells.
// for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) { // for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
// const int downwind_cell = ja_downw_[j]; // const int downwind_cell = ja_downw_[j];
// int ci = pos[downwind_cell]; // int ci = pos[downwind_cell];
// ++needs_update[ci]; // ++needs_update[ci];
// if (needs_update[ci] == num_upstream[ci]) { // if (needs_update[ci] == num_upstream[ci]) {
// fully_marked_stack.push_back(ci); // fully_marked_stack.push_back(ci);
// } // }
// } // }
// } // }
// // Unmark this cell. // // Unmark this cell.
// needs_update[fully_marked_ci] = 0; // needs_update[fully_marked_ci] = 0;
// } // }
if (!needs_update[i]) { if (!needs_update[i]) {
continue; continue;
} }
++update_count; ++update_count;
const int cell = cells[i]; const int cell = cells[i];
const double old_s = saturation_[cell]; const double old_s = saturation_[cell];
saturation_[cell] = s0[i]; saturation_[cell] = s0[i];
solveSingleCell(cell); solveSingleCell(cell);
const double s_change = std::fabs(saturation_[cell] - old_s); const double s_change = std::fabs(saturation_[cell] - old_s);
if (s_change > tol) { if (s_change > tol) {
// Mark downwind cells. // Mark downwind cells.
for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) { for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
const int downwind_cell = ja_downw_[j]; const int downwind_cell = ja_downw_[j];
int ci = pos[downwind_cell]; int ci = pos[downwind_cell];
if (ci != -1) { if (ci != -1) {
needs_update[ci] = 1; needs_update[ci] = 1;
} }
// ++needs_update[ci]; // ++needs_update[ci];
// if (needs_update[ci] == num_upstream[ci]) { // if (needs_update[ci] == num_upstream[ci]) {
// fully_marked_stack.push_back(ci); // fully_marked_stack.push_back(ci);
// } // }
} }
} }
// Unmark this cell. // Unmark this cell.
needs_update[i] = 0; needs_update[i] = 0;
} }
// std::cout << "Iter = " << num_iters << " update_count = " << update_count // std::cout << "Iter = " << num_iters << " update_count = " << update_count
// << " # marked cells = " // << " # marked cells = "
// << std::accumulate(needs_update.begin(), needs_update.end(), 0) << std::endl; // << std::accumulate(needs_update.begin(), needs_update.end(), 0) << std::endl;
} while (update_count > 0 && ++num_iters < max_iters); } while (update_count > 0 && ++num_iters < max_iters);
// Done with iterations, check if we succeeded. // Done with iterations, check if we succeeded.
if (update_count > 0) { if (update_count > 0) {
THROW("In solveMultiCell(), we did not converge after " THROW("In solveMultiCell(), we did not converge after "
<< num_iters << " iterations. Remaining update count = " << update_count); << num_iters << " iterations. Remaining update count = " << update_count);
} }
std::cout << "Solved " << num_cells << " cell multicell problem in " std::cout << "Solved " << num_cells << " cell multicell problem in "
<< num_iters << " iterations." << std::endl; << num_iters << " iterations." << std::endl;
} }
double TransportModelCompressibleTwophase::fracFlow(double s, int cell) const double TransportModelCompressibleTwophase::fracFlow(double s, int cell) const
{ {
double sat[2] = { s, 1.0 - s }; double sat[2] = { s, 1.0 - s };
double mob[2]; double mob[2];
props_.relperm(1, sat, &cell, mob, 0); props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[2*cell + 0]; mob[0] /= visc_[2*cell + 0];
mob[1] /= visc_[2*cell + 1]; mob[1] /= visc_[2*cell + 1];
return mob[0]/(mob[0] + mob[1]); return mob[0]/(mob[0] + mob[1]);
} }
@ -307,19 +309,19 @@ namespace Opm
// //
struct TransportModelCompressibleTwophase::GravityResidual struct TransportModelCompressibleTwophase::GravityResidual
{ {
int cell; int cell;
int nbcell[2]; int nbcell[2];
double s0; double s0;
double dtpv; // dt/pv(i) double dtpv; // dt/pv(i)
double gf[2]; double gf[2];
const TransportModelCompressibleTwophase& tm; const TransportModelCompressibleTwophase& tm;
explicit GravityResidual(const TransportModelCompressibleTwophase& tmodel, explicit GravityResidual(const TransportModelCompressibleTwophase& tmodel,
const std::vector<int>& cells, const std::vector<int>& cells,
const int pos, const int pos,
const double* gravflux) // Always oriented towards next in column. Size = colsize - 1. const double* gravflux) // Always oriented towards next in column. Size = colsize - 1.
: tm(tmodel) : tm(tmodel)
{ {
cell = cells[pos]; cell = cells[pos];
nbcell[0] = -1; nbcell[0] = -1;
gf[0] = 0.0; gf[0] = 0.0;
if (pos > 0) { if (pos > 0) {
@ -332,40 +334,40 @@ namespace Opm
nbcell[1] = cells[pos + 1]; nbcell[1] = cells[pos + 1];
gf[1] = gravflux[pos]; gf[1] = gravflux[pos];
} }
s0 = tm.saturation_[cell]; s0 = tm.saturation_[cell];
dtpv = tm.dt_/tm.porevolume_[cell]; dtpv = tm.dt_/tm.porevolume0_[cell];
} }
double operator()(double s) const double operator()(double s) const
{ {
double res = s - s0; double res = s - s0;
double mobcell[2]; double mobcell[2];
tm.mobility(s, cell, mobcell); tm.mobility(s, cell, mobcell);
for (int nb = 0; nb < 2; ++nb) { for (int nb = 0; nb < 2; ++nb) {
if (nbcell[nb] != -1) { if (nbcell[nb] != -1) {
double m[2]; double m[2];
if (gf[nb] < 0.0) { if (gf[nb] < 0.0) {
m[0] = mobcell[0]; m[0] = mobcell[0];
m[1] = tm.mob_[2*nbcell[nb] + 1]; m[1] = tm.mob_[2*nbcell[nb] + 1];
} else { } else {
m[0] = tm.mob_[2*nbcell[nb]]; m[0] = tm.mob_[2*nbcell[nb]];
m[1] = mobcell[1]; 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]); res += -dtpv*gf[nb]*m[0]*m[1]/(m[0] + m[1]);
} }
} }
} }
return res; return res;
} }
}; };
void TransportModelCompressibleTwophase::mobility(double s, int cell, double* mob) const void TransportModelCompressibleTwophase::mobility(double s, int cell, double* mob) const
{ {
double sat[2] = { s, 1.0 - s }; double sat[2] = { s, 1.0 - s };
props_.relperm(1, sat, &cell, mob, 0); props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0]; mob[0] /= visc_[0];
mob[1] /= visc_[1]; 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] = 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]); 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(); const int nc = cells.size();
std::vector<double> col_gravflux(nc - 1); std::vector<double> col_gravflux(nc - 1);
for (int ci = 0; ci < nc - 1; ++ci) { for (int ci = 0; ci < nc - 1; ++ci) {
const int cell = cells[ci]; const int cell = cells[ci];
const int next_cell = cells[ci + 1]; const int next_cell = cells[ci + 1];
for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) { for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) {
const int face = grid_.cell_faces[j]; const int face = grid_.cell_faces[j];
const int c1 = grid_.face_cells[2*face + 0]; const int c1 = grid_.face_cells[2*face + 0];
const int c2 = grid_.face_cells[2*face + 1]; 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]; const double gf = gravflux_[face];
col_gravflux[ci] = (c1 == cell) ? gf : -gf; col_gravflux[ci] = (c1 == cell) ? gf : -gf;
} }
} }
} }
// Store initial saturation s0 // Store initial saturation s0
@ -438,7 +440,7 @@ namespace Opm
} }
// Solve single cell problems, repeating if necessary. // Solve single cell problems, repeating if necessary.
double max_s_change = 0.0; double max_s_change = 0.0;
int num_iters = 0; int num_iters = 0;
do { do {
max_s_change = 0.0; max_s_change = 0.0;
@ -454,12 +456,12 @@ namespace Opm
std::fabs(saturation_[cells[ci2]] - old_s[1]))); std::fabs(saturation_[cells[ci2]] - old_s[1])));
} }
// std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change << std::endl; // 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_) { if (max_s_change > tol_) {
THROW("In solveGravityColumn(), we did not converge after " THROW("In solveGravityColumn(), we did not converge after "
<< num_iters << " iterations. Delta s = " << max_s_change); << num_iters << " iterations. Delta s = " << max_s_change);
} }
return num_iters + 1; return num_iters + 1;
} }
@ -467,6 +469,7 @@ namespace Opm
void TransportModelCompressibleTwophase::solveGravity(const std::vector<std::vector<int> >& columns, void TransportModelCompressibleTwophase::solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure, const double* pressure,
const double* porevolume0,
const double* porevolume, const double* porevolume,
const double dt, const double dt,
std::vector<double>& saturation) std::vector<double>& saturation)
@ -489,6 +492,7 @@ namespace Opm
} }
// Set up other variables. // Set up other variables.
porevolume0_ = porevolume0;
porevolume_ = porevolume; porevolume_ = porevolume;
dt_ = dt; dt_ = dt;
toWaterSat(saturation, saturation_); toWaterSat(saturation, saturation_);

4
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
View File

@ -47,6 +47,7 @@ namespace Opm
/// \param[in] darcyflux Array of signed face fluxes. /// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] pressure Array of cell pressures /// \param[in] pressure Array of cell pressures
/// \param[in] surfacevol0 Array of surface volumes at start of timestep /// \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] porevolume Array of pore volumes.
/// \param[in] source Transport source term. /// \param[in] source Transport source term.
/// \param[in] dt Time step. /// \param[in] dt Time step.
@ -54,6 +55,7 @@ namespace Opm
void solve(const double* darcyflux, void solve(const double* darcyflux,
const double* pressure, const double* pressure,
const double* surfacevol0, const double* surfacevol0,
const double* porevolume0,
const double* porevolume, const double* porevolume,
const double* source, const double* source,
const double dt, const double dt,
@ -71,6 +73,7 @@ namespace Opm
/// \TODO: Implement this. /// \TODO: Implement this.
void solveGravity(const std::vector<std::vector<int> >& columns, void solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure, const double* pressure,
const double* porevolume0,
const double* porevolume, const double* porevolume,
const double dt, const double dt,
std::vector<double>& saturation); std::vector<double>& saturation);
@ -96,6 +99,7 @@ namespace Opm
const double* darcyflux_; // one flux per grid face const double* darcyflux_; // one flux per grid face
const double* surfacevol0_; // one per phase per cell const double* surfacevol0_; // one per phase per cell
const double* porevolume0_; // one volume per cell
const double* porevolume_; // one volume per cell const double* porevolume_; // one volume per cell
const double* source_; // one source per cell const double* source_; // one source per cell
double dt_; double dt_;

13
tests/bo_resprop_test.cpp
View File

@ -24,6 +24,7 @@
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/eclipse/EclipseGridInspector.hpp> #include <opm/core/eclipse/EclipseGridInspector.hpp>
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp> #include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/grid.h>
#include <algorithm> #include <algorithm>
#include <iostream> #include <iostream>
@ -39,14 +40,10 @@ int main(int argc, char** argv)
// Parser. // Parser.
std::string ecl_file = param.get<std::string>("filename"); std::string ecl_file = param.get<std::string>("filename");
Opm::EclipseGridParser deck(ecl_file); Opm::EclipseGridParser deck(ecl_file);
Opm::EclipseGridInspector insp(deck); UnstructuredGrid grid;
std::tr1::array<int, 3> gs = insp.gridSize(); grid.number_of_cells = 1;
int num_cells = gs[0]*gs[1]*gs[2]; grid.global_cell = NULL;
std::vector<int> global_cell(num_cells); Opm::BlackoilPropertiesFromDeck props(deck, grid);
for (int i = 0; i < num_cells; ++i) {
global_cell[i] = i;
}
Opm::BlackoilPropertiesFromDeck props(deck, global_cell);
const int n = 1; const int n = 1;
double p[n] = { 150e5 }; double p[n] = { 150e5 };