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

Merge branch 'master' into ert

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-09-18 09:18:40 +02:00
commit 9153c7f87b
9 changed files with 120 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/sim_2p_comp_reorder.cpp
View File

@ -32,7 +32,6 @@
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/BlackoilPropertiesBasic.hpp>

2
examples/sim_wateroil.cpp
View File

@ -408,7 +408,7 @@ main(int argc, char** argv)
state.saturation(), state.surfacevol());
// Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
if (use_segregation_split) {
reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0],
reorder_model.solveGravity(columns,
stepsize, state.saturation(), state.surfacevol());
}
}

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

@ -36,6 +36,7 @@
#include <cmath>
#include <iostream>
#include <iomanip>
#include <numeric>
namespace Opm
{
@ -446,7 +447,7 @@ 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:
//
@ -458,18 +459,42 @@ namespace Opm
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];
}
}
}
}
}

7
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;
@ -245,6 +246,7 @@ namespace Opm
wells_(wells_manager.c_wells()),
src_(src),
bcs_(bcs),
gravity_(gravity),
psolver_(grid, props, rock_comp, linsolver,
param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
@ -445,14 +447,13 @@ namespace Opm
}
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &state.pressure()[0],
&porevol[0], &initial_porevol[0], &transport_src[0], stepsize,
&initial_porevol[0], &porevol[0], &transport_src[0], stepsize,
state.saturation(), state.surfacevol());
Opm::computeInjectedProduced(props_,
state.pressure(), state.surfacevol(), state.saturation(),
transport_src, stepsize, injected, produced);
if (gravity_ != 0 && use_segregation_split_) {
tsolver_.solveGravity(columns_, &state.pressure()[0], &initial_porevol[0],
stepsize, state.saturation(), state.surfacevol());
tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol());
}
}
transport_timer.stop();

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]);
}

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

@ -152,8 +152,8 @@ namespace Opm
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 ? B_cell*src_flux : 0.0;
outflux = !src_is_inflow ? B_cell*src_flux : 0.0;
influx = src_is_inflow ? src_flux : 0.0;
outflux = !src_is_inflow ? 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) {
@ -349,7 +349,7 @@ namespace Opm
gf[1] = gravflux[pos];
}
s0 = tm.saturation_[cell];
dtpv = tm.dt_/tm.porevolume0_[cell];
dtpv = tm.dt_/tm.porevolume_[cell];
}
double operator()(double s) const
@ -380,8 +380,8 @@ namespace Opm
{
double sat[2] = { s, 1.0 - s };
props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0];
mob[1] /= visc_[1];
mob[0] /= visc_[2*cell + 0];
mob[1] /= visc_[2*cell + 1];
}
@ -407,7 +407,7 @@ namespace Opm
{
// Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f)
// + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ]
// But b_w,i * rho_w,S = rho_w,i, which we conmpute with a call to props_.density().
// But b_w,i * rho_w,S = rho_w,i, which we compute with a call to props_.density().
// We assume that we already have stored T_ij in trans_.
// We also assume that the A_ matrices are updated from an earlier call to solve().
const int nc = grid_.number_of_cells;
@ -437,16 +437,15 @@ namespace Opm
void TransportModelCompressibleTwophase::solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux)
const int pos,
const double* gravflux)
{
const int cell = cells[pos];
GravityResidual res(*this, cells, pos, gravflux);
if (std::fabs(res(saturation_[cell])) > tol_) {
int iters_used;
saturation_[cell] = RootFinder::solve(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
saturation_[cell] = RootFinder::solve(res, saturation_[cell], 0.0, 1.0, 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]);
}
@ -506,8 +505,6 @@ namespace Opm
void TransportModelCompressibleTwophase::solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure,
const double* porevolume0,
const double dt,
std::vector<double>& saturation,
std::vector<double>& surfacevol)
@ -522,18 +519,22 @@ namespace Opm
cells[c] = c;
}
mob_.resize(2*nc);
std::vector<double> boths;
Opm::toBothSat(saturation, boths);
props_.relperm(cells.size(), &boths[0], &cells[0], &mob_[0], 0);
props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
for (int c = 0; c < nc; ++c) {
mob_[2*c + 0] /= visc_[2*c + 0];
mob_[2*c + 1] /= visc_[2*c + 1];
// props_.relperm(cells.size(), &saturation[0], &cells[0], &mob_[0], 0);
// props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
// for (int c = 0; c < nc; ++c) {
// mob_[2*c + 0] /= visc_[2*c + 0];
// mob_[2*c + 1] /= visc_[2*c + 1];
// }
const int np = props_.numPhases();
for (int cell = 0; cell < nc; ++cell) {
mobility(saturation_[cell], cell, &mob_[np*cell]);
}
// Set up other variables.
porevolume0_ = porevolume0;
dt_ = dt;
toWaterSat(saturation, saturation_);

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

@ -74,13 +74,10 @@ namespace Opm
/// vertical stack, that do not interact with other columns (for
/// gravity segregation.
/// \param[in] columns Vector of cell-columns.
/// \param[in] porevolume0 Array of pore volumes at start of timestep.
/// \param[in] dt Time step.
/// \param[in, out] saturation Phase saturations.
/// \param[in, out] surfacevol Surface volume densities for each phase.
void solveGravity(const std::vector<std::vector<int> >& columns,
const double* pressure,
const double* porevolume0,
const double dt,
std::vector<double>& saturation,
std::vector<double>& surfacevol);

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.

2
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;
}