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

Merge branch 'master' into reorder_tof

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2012-09-25 10:13:46 +02:00
parent d58314f624 4794341d06
commit 112aba50dd
15 changed files with 256 additions and 129 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@@ -29,6 +29,7 @@ ltmain.sh
m4/libtool.m4
m4/lt*.m4
missing
ar-lib
tutorials/tutorial[1-4]
# Ignoring executables
@@ -38,6 +39,7 @@ examples/spu_2p
examples/reorder-qfs
examples/refine_wells
examples/sim_2p_incomp_reorder
examples/sim_2p_comp_reorder
examples/sim_wateroil
examples/wells_example
tests/test_cfs_tpfa

3
Makefile.am
View File

@@ -306,3 +306,6 @@ opm/core/linalg/LinearSolverAGMG.hpp
lib_libopmcore_la_LDFLAGS += \
$(FCLIBS)
endif
pkgconfigdir = $(libdir)/pkgconfig
pkgconfig_DATA = lib/pkgconfig/opm-core.pc

51
README
View File

@@ -46,11 +46,17 @@ sudo apt-get install -y build-essential gfortran pkg-config libtool \
sudo apt-get install -y doxygen ghostscript texlive-latex-recommended pgf
# packages necessary for version control
sudo apt-get install -y git-core git-svn subversion
sudo apt-get install -y git-core
# libraries necessary for DUNE
# basic libraries necessary for both DUNE and OPM
sudo apt-get install -y libboost-all-dev libsuperlu3-dev libsuitesparse-dev
# add this repository for necessary backports (required for Ubuntu Precise)
sudo add-apt-repository -y ppa:opm/ppa
# parts of DUNE needed
sudo apt-get install -y libdune-common-dev libdune-istl-dev
# libraries necessary for OPM
sudo apt-get install -y libxml0-dev
@@ -59,33 +65,14 @@ DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS
-----------------------------------------
# libraries
sudo zypper install blas libblas3 lapack liblapack3 libboost libxml2 umfpack
sudo zypper in blas libblas3 lapack liblapack3 libboost libxml2 umfpack
# tools
sudo zypper install gcc automake autoconf git doxygen
sudo zypper in gcc automake autoconf git doxygen
RETRIEVING AND BUILDING DUNE PREREQUISITES
------------------------------------------
(only necessary if you want to use opm-core as a dune module)
# trust DUNE certificate (sic)
echo p | svn list https://svn.dune-project.org/svn/dune-common
# checkout DUNE libraries
for module in common istl geometry grid localfunctions; do
git svn clone -s \
https://svn.dune-project.org/svn/dune-$module/branches/release-2.2/ \
dune-$module
done
# building DUNE libraries
for module in common istl geometry grid localfunctions; do
env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=dune-$module \
--configure-opts="--enable-fieldvector-size-is-method" \
--make-opts="-j -l 0.8" autogen : configure : make
done
# DUNE libraries
sudo zypper ar http://download.opensuse.org/repositories/science/openSUSE_12.2/science.repo
sudo zypper in dune-common dune-istl
DOWNLOADING
@@ -100,20 +87,10 @@ If you want to contribute, fork OPM/opm-core on github.
BUILDING
--------
(standalone opm-core:)
cd ../opm-core
cd opm-core
autoreconf -i
./configure
make
sudo make install
(using opm-core as a dune module:)
# note: this is done from the parent directory of opm-core
env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=opm-core \
--configure-opts="" --make-opts="-j -l 0.8" autogen : configure : make
DOCUMENTATION

9
configure.ac
View File

@@ -2,7 +2,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.59])
AC_INIT([OPM Core Library], [0.1], [atgeirr@sintef.no],dnl
AC_INIT([OPM Core Library], [0.1], [atgeirr@sintef.no],
[opmcore], [https://public.ict.sintef.no/opm/hg/opmcore])
AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects])
@@ -20,8 +20,11 @@ AC_CONFIG_HEADERS([config.h])
AC_PROG_CC
AM_PROG_CC_C_O
dnl Initialize libtool; the funny indentation here is to
dnl satisfy libtoolize' check for the presence of this macro
m4_ifdef([LT_INIT],
[LT_INIT[]dnl
[
LT_INIT[]dnl
LT_LANG([C++])dnl
LT_LANG([Fortran 77])dnl
LT_LANG([Fortran])dnl
@@ -51,6 +54,8 @@ AC_CONFIG_FILES([
tests/Makefile
examples/Makefile
tutorials/Makefile
opm-core.pc
lib/pkgconfig/opm-core.pc
])
AC_OUTPUT

11
lib/pkgconfig/opm-core.pc.in Normal file
View File

@@ -0,0 +1,11 @@
prefix=@prefix@
exec_prefix=@exec_prefix@
libdir=@libdir@
includedir=@includedir@
Name: @PACKAGE_NAME@
Description: @PACKAGE_STRING@
Version: @PACKAGE_VERSION@
URL: @PACKAGE_URL@
Libs: -L${libdir} -l@PACKAGE@
Cflags: -I${includedir}

18
opm-core.pc.in Normal file
View File

@@ -0,0 +1,18 @@
# This is the configuration for local builds. Use this by putting the
# compilation output path (the directory in which you ran ./configure)
# into the environment variable PKG_CONFIG_PATH. This will enable you
# to use pkg-config in your code while making changes to opm-core.
# This is NOT the file that is installed in the system directories when
# you do `make install`. That is the one in lib/pkgconfig. However, if
# you make changes here, you should consider that one as well.
libdir=@abs_top_builddir@/lib/.libs
includedir=@abs_top_srcdir@
Name: @PACKAGE_NAME@
Description: @PACKAGE_STRING@
Version: @PACKAGE_VERSION@
URL: @PACKAGE_URL@
Libs: -L${libdir} -l@PACKAGE@
Cflags: -I${includedir}

4
opm/core/linalg/LinearSolverIstl.cpp
View File

@@ -26,6 +26,10 @@
#include <opm/core/utility/have_boost_redef.hpp>
// Silence compatibility warning from DUNE headers since we don't use
// the deprecated member anyway (in this compilation unit)
#define DUNE_COMMON_FIELDVECTOR_SIZE_IS_METHOD 1
// TODO: clean up includes.
#include <dune/common/deprecated.hh>
#include <dune/istl/bvector.hh>

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

@@ -36,6 +36,7 @@
#include <cmath>
#include <iostream>
#include <iomanip>
#include <numeric>
namespace Opm
{
@@ -122,7 +123,7 @@ namespace Opm
WellState& well_state)
{
const int nc = grid_.number_of_cells;
const int nw = wells_->number_of_wells;
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
// Set up dynamic data.
computePerSolveDynamicData(dt, state, well_state);
@@ -446,30 +447,54 @@ namespace Opm
/// Compute per-iteration dynamic properties for wells.
void CompressibleTpfa::computeWellDynamicData(const double /*dt*/,
const BlackoilState& /*state*/,
const WellState& /*well_state*/)
const WellState& well_state)
{
// These are the variables that get computed by this function:
//
// std::vector<double> wellperf_A_;
// std::vector<double> wellperf_phasemob_;
const int np = props_.numPhases();
const int nw = wells_->number_of_wells;
const int nperf = wells_->well_connpos[nw];
const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0;
wellperf_A_.resize(nperf*np*np);
wellperf_phasemob_.resize(nperf*np);
// The A matrix is set equal to the perforation grid cells'
// matrix, for both injectors and producers.
// matrix for producers, computed from bhp and injection
// component fractions from
// The mobilities are set equal to the perforation grid cells'
// mobilities, for both injectors and producers.
// mobilities for producers.
std::vector<double> mu(np);
for (int w = 0; w < nw; ++w) {
bool producer = (wells_->type[w] == PRODUCER);
const double* comp_frac = &wells_->comp_frac[np*w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
const int c = wells_->well_cells[j];
const double* cA = &cell_A_[np*np*c];
double* wpA = &wellperf_A_[np*np*j];
std::copy(cA, cA + np*np, wpA);
const double* cM = &cell_phasemob_[np*c];
double* wpM = &wellperf_phasemob_[np*j];
std::copy(cM, cM + np, wpM);
if (producer) {
const double* cA = &cell_A_[np*np*c];
std::copy(cA, cA + np*np, wpA);
const double* cM = &cell_phasemob_[np*c];
std::copy(cM, cM + np, wpM);
} else {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
// Hack warning: comp_frac is used as a component
// surface-volume variable in calls to matrix() and
// viscosity(), but as a saturation in the call to
// relperm(). This is probably ok as long as injectors
// only inject pure fluids.
props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL);
props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL);
ASSERT(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6);
props_.relperm (1, comp_frac, &c, wpM , NULL);
for (int phase = 0; phase < np; ++phase) {
wpM[phase] /= mu[phase];
}
}
}
}
}
@@ -487,9 +512,9 @@ namespace Opm
const double* z = &state.surfacevol()[0];
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = &wellperf_gpot_[0];
completion_data.A = &wellperf_A_[0];
completion_data.phasemob = &wellperf_phasemob_[0];
completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0;
completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
cfs_tpfa_res_wells wells_tmp;
wells_tmp.W = const_cast<Wells*>(wells_);
wells_tmp.data = &completion_data;
@@ -574,9 +599,9 @@ namespace Opm
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
CompletionData completion_data;
completion_data.gpot = const_cast<double*>(&wellperf_gpot_[0]);
completion_data.A = const_cast<double*>(&wellperf_A_[0]);
completion_data.phasemob = const_cast<double*>(&wellperf_phasemob_[0]);
completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast<double*>(&wellperf_gpot_[0]) : 0;
completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
cfs_tpfa_res_wells wells_tmp;
wells_tmp.W = const_cast<Wells*>(wells_);
wells_tmp.data = &completion_data;
@@ -584,6 +609,9 @@ namespace Opm
forces.wells = &wells_tmp;
forces.src = NULL;
double* wpress = ! well_state.bhp ().empty() ? & well_state.bhp ()[0] : 0;
double* wflux = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0;
cfs_tpfa_res_flux(gg,
&forces,
props_.numPhases(),
@@ -592,9 +620,9 @@ namespace Opm
&face_phasemob_[0],
&face_gravcap_[0],
&state.pressure()[0],
&well_state.bhp()[0],
wpress,
&state.faceflux()[0],
&well_state.perfRates()[0]);
wflux);
cfs_tpfa_res_fpress(gg,
props_.numPhases(),
&htrans_[0],

84
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@@ -142,7 +142,7 @@ impl_allocate(struct UnstructuredGrid *G ,
nnu = G->number_of_cells;
nwperf = 0;
if (wells != NULL) { assert (wells->W != NULL);
if ((wells != NULL) && (wells->W != NULL)) {
nnu += wells->W->number_of_wells;
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
}
@@ -185,13 +185,15 @@ construct_matrix(struct UnstructuredGrid *G ,
/* ---------------------------------------------------------------------- */
{
int f, c1, c2, w, i, nc, nnu;
int wells_present;
size_t nnz;
struct CSRMatrix *A;
nc = nnu = G->number_of_cells;
if (wells != NULL) {
assert (wells->W != NULL);
wells_present = (wells != NULL) && (wells->W != NULL);
if (wells_present) {
nnu += wells->W->number_of_wells;
}
@@ -214,7 +216,7 @@ construct_matrix(struct UnstructuredGrid *G ,
}
}
if (wells != NULL) {
if (wells_present) {
/* Well <-> cell connections */
struct Wells *W = wells->W;
@@ -252,7 +254,7 @@ construct_matrix(struct UnstructuredGrid *G ,
}
}
if (wells != NULL) {
if (wells_present) {
/* Fill well <-> cell connections */
struct Wells *W = wells->W;
@@ -741,6 +743,29 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt,
}
/* ---------------------------------------------------------------------- */
static void
welleq_coeff_shut(int np, struct cfs_tpfa_res_data *h,
double *res, double *w2c, double *w2w)
/* ---------------------------------------------------------------------- */
{
int p;
double fwi;
const double *dpflux_w;
/* Sum reservoir phase flux derivatives set by
* compute_darcyflux_and_deriv(). */
dpflux_w = h->pimpl->flux_work + (1 * np);
for (p = 0, fwi = 0.0; p < np; p++) {
fwi += dpflux_w[ p ];
}
*res = 0.0;
*w2c = 0.0;
*w2w = fwi;
}
/* ---------------------------------------------------------------------- */
static void
welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h,
@@ -815,23 +840,25 @@ assemble_completion_to_well(int w, int c, int nc, int np,
ctrl = W->ctrls[ w ];
if (ctrl->current < 0) {
assert (0); /* Shut wells currently not supported */
/* Interpreting a negative current control index to mean a shut well */
welleq_coeff_shut(np, h, &res, &w2c, &w2w);
}
else {
switch (ctrl->type[ ctrl->current ]) {
case BHP :
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
h, &res, &w2c, &w2w);
break;
switch (ctrl->type[ ctrl->current ]) {
case BHP :
welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
h, &res, &w2c, &w2w);
break;
case RESERVOIR_RATE:
assert (W->number_of_phases == np);
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
break;
case RESERVOIR_RATE:
assert (W->number_of_phases == np);
welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w);
break;
case SURFACE_RATE:
assert (0); /* Surface rate currently not supported */
break;
case SURFACE_RATE:
assert (0); /* Surface rate currently not supported */
break;
}
}
/* Assemble completion contributions */
@@ -854,7 +881,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
struct cfs_tpfa_res_data *h )
{
int w, i, c, np, np2, nc;
int is_neumann;
int is_neumann, is_open;
double pw, dp;
double *WI, *gpot, *pmobp;
const double *Ac, *dAc;
@@ -876,6 +903,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
for (w = i = 0; w < W->number_of_wells; w++) {
pw = wpress[ w ];
is_open = W->ctrls[w]->current >= 0;
for (; i < W->well_connpos[w + 1];
i++, gpot += np, pmobp += np) {
@@ -888,7 +916,9 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
init_completion_contrib(i, np, Ac, dAc, h->pimpl);
assemble_completion_to_cell(c, nc + w, np, dt, h);
if (is_open) {
assemble_completion_to_cell(c, nc + w, np, dt, h);
}
/* Prepare for RESV controls */
compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
@@ -1127,8 +1157,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid *G ,
nf = G->number_of_faces;
nwperf = 0;
if (wells != NULL) {
assert (wells->W != NULL);
if ((wells != NULL) && (wells->W != NULL)) {
nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
}
@@ -1194,7 +1223,9 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid *G ,
assemble_cell_contrib(G, c, h);
}
if ((forces != NULL) && (forces->wells != NULL)) {
if ((forces != NULL) &&
(forces->wells != NULL) &&
(forces->wells->W != NULL)) {
compute_well_compflux_and_deriv(forces->wells, cq->nphases,
cpress, wpress, h->pimpl);
@@ -1297,8 +1328,9 @@ cfs_tpfa_res_flux(struct UnstructuredGrid *G ,
{
compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux);
if ((forces != NULL) && (forces->wells != NULL) &&
(wpress != NULL) && (wflux != NULL)) {
if ((forces != NULL) &&
(forces->wells != NULL) &&
(forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) {
compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux);
}

1
opm/core/simulator/SimulatorCompressibleTwophase.cpp
View File

@@ -175,6 +175,7 @@ namespace Opm
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["surfvolume"] = &state.surfacevol();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;

13
opm/core/simulator/WellState.hpp
View File

@@ -37,10 +37,17 @@ namespace Opm
if (wells) {
const int nw = wells->number_of_wells;
bhp_.resize(nw);
// Initialize bhp to be pressure in first perforation cell.
// Initialize bhp to be target pressure
// if bhp-controlled well, otherwise set
// to pressure in first perforation cell.
for (int w = 0; w < nw; ++w) {
const int cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = state.pressure()[cell];
const WellControls* ctrl = wells->ctrls[w];
if (ctrl->type[ctrl->current] == BHP) {
bhp_[w] = ctrl->target[ctrl->current];
} else {
const int cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = state.pressure()[cell];
}
}
perfrates_.resize(wells->well_connpos[nw]);
}

1
opm/core/utility/Units.hpp
View File

@@ -101,6 +101,7 @@ namespace Opm
/// @{
const double gallon = 231 * cubic(inch);
const double stb = 42 * gallon;
const double liter = 1 * cubic(prefix::deci*meter);
/// @}
/// \name Mass

66
opm/core/utility/initState_impl.hpp
View File

@@ -29,6 +29,7 @@
#include <opm/core/utility/Units.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <cmath>
namespace Opm
@@ -512,15 +513,23 @@ namespace Opm
State& state)
{
const int num_phases = props.numPhases();
const PhaseUsage pu = phaseUsageFromDeck(deck);
if (num_phases != pu.num_phases) {
THROW("initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
state.init(grid, num_phases);
if (deck.hasField("EQUIL")) {
if (num_phases != 2) {
THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
THROW("initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
}
// Set saturations depending on oil-water contact.
const EQUIL& equil= deck.getEQUIL();
if (equil.equil.size() != 1) {
THROW("No region support yet.");
THROW("initStateFromDeck(): No region support yet.");
}
const double woc = equil.equil[0].water_oil_contact_depth_;
initWaterOilContact(grid, props, woc, WaterBelow, state);
@@ -528,37 +537,64 @@ namespace Opm
const double datum_z = equil.equil[0].datum_depth_;
const double datum_p = equil.equil[0].datum_depth_pressure_;
initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
} else if (deck.hasField("SWAT") && deck.hasField("PRESSURE")) {
// Set saturations from SWAT, pressure from PRESSURE.
} else if (deck.hasField("PRESSURE")) {
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
std::vector<double>& s = state.saturation();
std::vector<double>& p = state.pressure();
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
const int num_cells = grid.number_of_cells;
if (num_phases == 2) {
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c] = sw_deck[c_deck];
s[2*c + 1] = 1.0 - s[2*c];
p[c] = p_deck[c_deck];
if (!pu.phase_used[BlackoilPhases::Aqua]) {
// oil-gas: we require SGAS
if (!deck.hasField("SGAS")) {
THROW("initStateFromDeck(): missing SGAS keyword in 2-phase init");
}
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c + gpos] = sg_deck[c_deck];
s[2*c + opos] = 1.0 - sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
// water-oil or water-gas: we require SWAT
if (!deck.hasField("SWAT")) {
THROW("initStateFromDeck(): missing SWAT keyword in 2-phase init");
}
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int nwpos = (wpos + 1) % 2;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c + wpos] = sw_deck[c_deck];
s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
p[c] = p_deck[c_deck];
}
}
} else if (num_phases == 3) {
if (!deck.hasField("SGAS")) {
THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found).");
const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS");
if (!has_swat_sgas) {
THROW("initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
}
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[3*c] = sw_deck[c_deck];
s[3*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
s[3*c + 2] = sg_deck[c_deck];
s[3*c + wpos] = sw_deck[c_deck];
s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
s[3*c + gpos] = sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
}
} else {
THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE.");
THROW("initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
}
// Finally, init face pressures.

5
opm/core/utility/miscUtilities.cpp
View File

@@ -374,7 +374,10 @@ namespace Opm
if (perf_rate > 0.0) {
// perf_rate is a total inflow rate, we want a water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well with index " << w << " ignored." << std::endl;
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);

53
opm/core/wells/WellsManager.cpp
View File

@@ -544,7 +544,7 @@ namespace Opm
cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
} else if (wci_line.injector_type_[0] == 'G') {
if (!pu.phase_used[BlackoilPhases::Vapour]) {
THROW("Water phase not used, yet found water-injecting well.");
THROW("Gas phase not used, yet found gas-injecting well.");
}
cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
}
@@ -716,32 +716,31 @@ namespace Opm
#endif
if (deck.hasField("WELOPEN")) {
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
int index = it->second;
if (line.openshutflag_ == "SHUT") {
// We currently don't care if the well is open or not.
/// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }?
//ASSERT(w_->ctrls[index]->current < 0);
} else if (line.openshutflag_ == "OPEN") {
//ASSERT(w_->ctrls[index]->current >= 0);
} else {
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
// We revert back to it's original control.
// Note that this is OK as ~~ = id.
w_->ctrls[index]->current = ~w_->ctrls[index]->current;
}
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
const int index = it->second;
if (line.openshutflag_ == "SHUT") {
int& cur_ctrl = w_->ctrls[index]->current;
if (cur_ctrl >= 0) {
cur_ctrl = ~cur_ctrl;
}
ASSERT(w_->ctrls[index]->current < 0);
} else if (line.openshutflag_ == "OPEN") {
int& cur_ctrl = w_->ctrls[index]->current;
if (cur_ctrl < 0) {
cur_ctrl = ~cur_ctrl;
}
ASSERT(w_->ctrls[index]->current >= 0);
} else {
THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
}
}
// Build the well_collection_ well group hierarchy.