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

Merge from upstream.

Browse Source
This commit is contained in:
Bård Skaflestad
2012-03-08 18:47:37 +01:00
parent 47ae9a454d f919da439b
commit ac18454902
8 changed files with 212 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
examples/spu_2p.cpp
View File

@@ -96,15 +96,32 @@
class ReservoirState {
public:
ReservoirState(const UnstructuredGrid* g, const int num_phases = 2)
ReservoirState(const UnstructuredGrid* g, const double init_sat = 0.0)
: press_ (g->number_of_cells, 0.0),
fpress_(g->number_of_faces, 0.0),
flux_ (g->number_of_faces, 0.0),
sat_ (num_phases * g->number_of_cells, 0.0)
sat_ (2 * g->number_of_cells, 0.0)
{
for (int cell = 0; cell < g->number_of_cells; ++cell) {
sat_[num_phases*cell + num_phases - 1] = 1.0;
}
for (int cell = 0; cell < g->number_of_cells; ++cell) {
sat_[2*cell] = init_sat;
sat_[2*cell + 1] = 1.0 - init_sat;
}
}
void setToMinimumWaterSat(const Opm::IncompPropertiesInterface& props)
{
const int n = props.numCells();
std::vector<int> cells(n);
for (int i = 0; i < n; ++i) {
cells[i] = i;
}
std::vector<double> smin(2*n);
std::vector<double> smax(2*n);
props.satRange(n, &cells[0], &smin[0], &smax[0]);
for (int cell = 0; cell < n; ++cell) {
sat_[2*cell] = smin[2*cell];
sat_[2*cell + 1] = 1.0 - smin[2*cell];
}
}
int numPhases() const { return sat_.size()/press_.size(); }
@@ -378,7 +395,11 @@ main(int argc, char** argv)
int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> totmob;
std::vector<double> omega; // Will remain empty if no gravity.
ReservoirState state(grid->c_grid(), props->numPhases());
double init_sat = param.getDefault("init_sat", 0.0);
ReservoirState state(grid->c_grid(), init_sat);
if (!param.has("init_sat")) {
state.setToMinimumWaterSat(*props);
}
// We need a separate reorder_sat, because the reorder
// code expects a scalar sw, not both sw and so.
std::vector<double> reorder_sat(num_cells);

2
opm/core/GridManager.cpp
View File

@@ -50,7 +50,7 @@ namespace Opm
struct grdecl grdecl;
grdecl.zcorn = &zcorn[0];
grdecl.coord = &coord[0];
grdecl.actnum = &actnum[0];
grdecl.actnum = actnum.empty() ? 0 : &actnum[0];
grdecl.dims[0] = dims[0];
grdecl.dims[1] = dims[1];
grdecl.dims[2] = dims[2];

4
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
View File

@@ -61,7 +61,7 @@ namespace Opm
if (deck.hasField("PVTW")) {
const std::vector<double>& pvtw = deck.getPVTW().pvtw_[region_number];
if (pvtw[2] != 0.0 || pvtw[4] != 0.0) {
THROW("PvtPropertiesIncompFromDeck::init() -- must have no compressibility effects in PVTW.");
MESSAGE("Compressibility effects in PVTW are ignored.");
}
viscosity_[phase_usage.phase_pos[PhaseUsage::Aqua]] = pvtw[3];
} else {
@@ -74,7 +74,7 @@ namespace Opm
if (deck.hasField("PVCDO")) {
const std::vector<double>& pvcdo = deck.getPVCDO().pvcdo_[region_number];
if (pvcdo[2] != 0.0 || pvcdo[4] != 0.0) {
THROW("PvtPropertiesIncompFromDeck::init() -- must have no compressibility effects in PVCDO.");
MESSAGE("Compressibility effects in PVCDO are ignored.");
}
viscosity_[phase_usage.phase_pos[PhaseUsage::Liquid]] = pvcdo[3];
} else {

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

@@ -58,7 +58,7 @@ namespace Opm
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage_.phase_pos[Aqua]] = sw[0];
smax_[phase_usage_.phase_pos[Aqua]] = sw.back();
smax_[phase_usage_.phase_pos[Aqua]] = 1.0;
}
if (phase_usage_.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;

135
opm/core/pressure/FlowBCManager.cpp
View File

@@ -20,10 +20,19 @@
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/grid.h>
#include <vector>
namespace Opm
{
namespace
{
std::string sideString(FlowBCManager::Side s);
void findSideFaces(const UnstructuredGrid& grid,
const FlowBCManager::Side side,
std::vector<int>& faces);
} // anon namespace
/// Default constructor sets up empty boundary conditions.
@@ -74,6 +83,61 @@ namespace Opm
}
/// Add BC_PRESSURE boundary conditions to all faces on a given side.
/// The grid must have a logical cartesian structure, and grid
/// faces must be tagged (i.e. grid.cell_facetag must be
/// non-null). Only the set of faces adjacent to cells with
/// minimum/maximum I/J/K coordinate (depending on side) are
/// considered.
void FlowBCManager::pressureSide(const UnstructuredGrid& grid,
const Side side,
const double pressure)
{
std::vector<int> faces;
findSideFaces(grid, side, faces);
int ok = flow_conditions_append_multi(BC_PRESSURE, faces.size(), &faces[0], pressure, bc_);
if (!ok) {
THROW("Failed to append pressure boundary conditions for side " << sideString(side));
}
}
/// Add BC_FLUX_TOTVOL boundary conditions to all faces on a given side.
/// The grid must have a logical cartesian structure, and grid
/// faces must be tagged (i.e. grid.cell_facetag must be
/// non-null). Only the set of faces adjacent to cells with
/// minimum/maximum I/J/K coordinate (depending on side) are
/// considered.
/// The flux specified is taken to be the total flux through
/// the side, each individual face receiving a part of the
/// total flux in proportion to its area, so that all faces
/// will have identical normal velocities.
void FlowBCManager::fluxSide(const UnstructuredGrid& grid,
const Side side,
const double flux)
{
// Find side faces.
std::vector<int> faces;
findSideFaces(grid, side, faces);
// Compute total area of faces.
double tot_area = 0.0;
for (int fi = 0; fi < int(faces.size()); ++fi) {
tot_area += grid.face_areas[faces[fi]];
}
// Append flux conditions for all the faces individually.
for (int fi = 0; fi < int(faces.size()); ++fi) {
const double face_flux = flux * grid.face_areas[faces[fi]] / tot_area;
int ok = flow_conditions_append(BC_FLUX_TOTVOL, faces[fi], face_flux, bc_);
if (!ok) {
THROW("Failed to append flux boundary conditions for face " << faces[fi] << " on side " << sideString(side));
}
}
}
/// Access the managed boundary conditions.
/// The method is named similarly to c_str() in std::string,
/// to make it clear that we are returning a C-compatible struct.
@@ -84,4 +148,75 @@ namespace Opm
// ------ Utility functions ------
namespace
{
std::string sideString(FlowBCManager::Side s)
{
switch (s) {
case FlowBCManager::Xmin: return "Xmin";
case FlowBCManager::Xmax: return "Xmax";
case FlowBCManager::Ymin: return "Ymin";
case FlowBCManager::Ymax: return "Ymax";
case FlowBCManager::Zmin: return "Zmin";
case FlowBCManager::Zmax: return "Zmax";
default: THROW("Unknown side tag " << s);
}
}
void cartCoord(const int log_cart_coord,
const int* dims,
int* ijk)
{
ijk[2] = log_cart_coord / (dims[0]*dims[1]);
ijk[1] = (log_cart_coord - ijk[2]*(dims[0]*dims[1])) / dims[0];
ijk[0] = log_cart_coord % dims[0];
}
/// The grid must have a logical cartesian structure, and grid
/// faces must be tagged (i.e. grid.cell_facetag must be
/// non-null). Only the set of faces adjacent to cells with
/// minimum/maximum I/J/K coordinate (depending on side) are
/// considered.
void findSideFaces(const UnstructuredGrid& grid,
const FlowBCManager::Side side,
std::vector<int>& faces)
{
if (grid.cell_facetag == 0) {
THROW("Faces not tagged - cannot extract " << sideString(side) << " faces.");
}
// Get all boundary faces with the correct tag and with
// min/max i/j/k (depending on side).
const int correct_ijk = (side % 2) ? grid.cartdims[side/2] : 0;
for (int c = 0; c < grid.number_of_cells; ++c) {
int ijk[3] = { -1 };
cartCoord(grid.global_cell[c], grid.cartdims, ijk);
if (ijk[side/2] != correct_ijk) {
continue;
}
for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) {
if (grid.cell_facetag[hf] == side) {
// Tag is correct.
const int f = grid.cell_faces[hf];
if (grid.face_cells[2*f] == -1 || grid.face_cells[2*f + 1] == -1) {
// Face is on boundary.
faces.push_back(f);
} else {
THROW("Face not on boundary, even with correct tag and boundary cell. This should not occur.");
}
}
}
}
}
} // anon namespace
} // namespace Opm

29
opm/core/pressure/FlowBCManager.hpp
View File

@@ -22,6 +22,8 @@
#include <opm/core/pressure/flow_bc.h>
struct UnstructuredGrid;
namespace Opm
{
@@ -57,6 +59,33 @@ namespace Opm
const int face,
const double value);
/// Defines the canonical sides for logical cartesian grids.
enum Side { Xmin, Xmax, Ymin, Ymax, Zmin, Zmax };
/// Add BC_PRESSURE boundary conditions to all faces on a given side.
/// The grid must have a logical cartesian structure, and grid
/// faces must be tagged (i.e. grid.cell_facetag must be
/// non-null). Only the set of faces adjacent to cells with
/// minimum/maximum I/J/K coordinate (depending on side) are
/// considered.
void pressureSide(const UnstructuredGrid& grid,
const Side side,
const double pressure);
/// Add BC_FLUX_TOTVOL boundary conditions to all faces on a given side.
/// The grid must have a logical cartesian structure, and grid
/// faces must be tagged (i.e. grid.cell_facetag must be
/// non-null). Only the set of faces adjacent to cells with
/// minimum/maximum I/J/K coordinate (depending on side) are
/// considered.
/// The flux specified is taken to be the total flux through
/// the side, each individual face receiving a part of the
/// total flux in proportion to its area, so that all faces
/// will have identical normal velocities.
void fluxSide(const UnstructuredGrid& grid,
const Side side,
const double flux);
/// Access the managed boundary conditions.
/// The method is named similarly to c_str() in std::string,
/// to make it clear that we are returning a C-compatible struct.

10
opm/core/pressure/flow_bc.c
View File

@@ -47,7 +47,7 @@ alloc_size(size_t n, size_t c)
/* ---------------------------------------------------------------------- */
/* Put structure in a well-defined, initial state */
/* ---------------------------------------------------------------------- */
static void
static int
initialise_structure(struct FlowBoundaryConditions *fbc)
/* ---------------------------------------------------------------------- */
{
@@ -55,10 +55,12 @@ initialise_structure(struct FlowBoundaryConditions *fbc)
fbc->cond_cpty = 0;
fbc->face_cpty = 0;
fbc->cond_pos = NULL;
fbc->cond_pos = malloc(1 * sizeof *fbc->cond_pos);
fbc->type = NULL;
fbc->value = NULL;
fbc->face = NULL;
return fbc->cond_pos != NULL;
}
@@ -130,9 +132,9 @@ flow_conditions_construct(size_t nbc)
fbc = malloc(1 * sizeof *fbc);
if (fbc != NULL) {
initialise_structure(fbc);
ok = initialise_structure(fbc);
ok = expand_tables(nbc, nbc, fbc);
ok = ok && expand_tables(nbc, nbc, fbc);
if (! ok) {
flow_conditions_destroy(fbc);

16
opm/core/utility/cpgpreprocess/preprocess.c
View File

@@ -434,12 +434,18 @@ void process_grdecl(const struct grdecl *in,
/* Permute actnum */
iptr = actnum;
for (j=0; j<ny; ++j){
for (i=0; i<nx; ++i){
for (k=0; k<nz; ++k){
*iptr++ = in->actnum[i+nx*(j+ny*k)];
if (in->actnum != NULL) {
for (j=0; j<ny; ++j){
for (i=0; i<nx; ++i){
for (k=0; k<nz; ++k){
*iptr++ = in->actnum[i+nx*(j+ny*k)];
}
}
}
} else {
for (i=0; i<nx*ny*nz; ++i){
*iptr++ = 1;
}
}
g.actnum = actnum;
@@ -460,7 +466,7 @@ void process_grdecl(const struct grdecl *in,
int c1 = i/2 + nx*(j/2 + ny*(k/2));
int c2 = i/2 + nx*(j/2 + ny*((k+1)/2));
if (in->actnum[c1] && in->actnum[c2] && (z2 < z1)){
if (((in->actnum == NULL) || (in->actnum[c1] && in->actnum[c2])) && (z2 < z1)){
fprintf(stderr, "\nZCORN should be strictly nondecreasing along pillars!\n");
error = 1;