OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into blattms-master-refactor-for-cpgrid-support

Browse Source 
Conflicts:
	examples/sim_fibo_ad.cpp
	opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

This brings the "CpGrid support" branch up to date with respect to
recent changes in opm-autodiff master.
This commit is contained in:
Bård Skaflestad
2014-04-08 16:31:17 +02:00
parent 8df314e90f 79afd9850c
commit a784d50f1e
14 changed files with 633 additions and 334 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
opm/autodiff/AutoDiffBlock.hpp
View File

@@ -194,6 +194,32 @@ namespace Opm
return *this;
}
/// Elementwise operator -=
AutoDiffBlock& operator-=(const AutoDiffBlock& rhs)
{
if (jac_.empty()) {
const int num_blocks = rhs.numBlocks();
jac_.resize(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jac_[block] = -rhs.jac_[block];
}
} else if (!rhs.jac_.empty()) {
assert (numBlocks() == rhs.numBlocks());
assert (value().size() == rhs.value().size());
const int num_blocks = numBlocks();
for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols());
jac_[block] -= rhs.jac_[block];
}
}
val_ -= rhs.val_;
return *this;
}
/// Elementwise operator +
AutoDiffBlock operator+(const AutoDiffBlock& rhs) const
{

34
opm/autodiff/AutoDiffHelpers.hpp
View File

@@ -288,7 +288,10 @@ spdiag(const AutoDiffBlock<double>::V& d)
public:
typedef AutoDiffBlock<Scalar> ADB;
Selector(const typename ADB::V& selection_basis)
enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero };
Selector(const typename ADB::V& selection_basis,
CriterionForLeftElement crit = GreaterEqualZero)
{
// Define selector structure.
const int n = selection_basis.size();
@@ -296,10 +299,33 @@ spdiag(const AutoDiffBlock<double>::V& d)
left_elems_.reserve(n);
right_elems_.reserve(n);
for (int i = 0; i < n; ++i) {
if (selection_basis[i] < 0.0) {
right_elems_.push_back(i);
} else {
bool chooseleft = false;
switch (crit) {
case GreaterEqualZero:
chooseleft = selection_basis[i] >= 0.0;
break;
case GreaterZero:
chooseleft = selection_basis[i] > 0.0;
break;
case Zero:
chooseleft = selection_basis[i] == 0.0;
break;
case NotEqualZero:
chooseleft = selection_basis[i] != 0.0;
break;
case LessZero:
chooseleft = selection_basis[i] < 0.0;
break;
case LessEqualZero:
chooseleft = selection_basis[i] <= 0.0;
break;
default:
OPM_THROW(std::logic_error, "No such criterion: " << crit);
}
if (chooseleft) {
left_elems_.push_back(i);
} else {
right_elems_.push_back(i);
}
}
}

21
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@@ -138,6 +138,7 @@ namespace Opm {
std::vector<ReservoirResidualQuant> rq_;
std::vector<PhasePresence> phaseCondition_;
V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
// The mass_balance vector has one element for each active phase,
// each of which has size equal to the number of cells.
@@ -161,10 +162,28 @@ namespace Opm {
computeAccum(const SolutionState& state,
const int aix );
void computeWellConnectionPressures(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw);
void
addOldWellEq(const SolutionState& state);
void
addWellControlEq(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw);
void
addWellEq(const SolutionState& state,
WellStateFullyImplicitBlackoil& xw);
void updateWellControls(const ADB& bhp,
const ADB& well_phase_flow_rate,
WellStateFullyImplicitBlackoil& xw) const;
void
assemble(const V& dtpv,
const BlackoilState& x,
const WellStateFullyImplicitBlackoil& xw);
WellStateFullyImplicitBlackoil& xw);
V solveJacobianSystem() const;

586
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@@ -19,12 +19,12 @@
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/core/grid.h>
@@ -32,6 +32,7 @@
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <opm/core/well_controls.h>
#include <cassert>
@@ -239,6 +240,7 @@ namespace {
{
const SolutionState state = constantState(x, xw);
computeAccum(state, 0);
computeWellConnectionPressures(state, xw);
}
const double atol = 1.0e-12;
@@ -345,91 +347,19 @@ namespace {
FullyImplicitBlackoilSolver<T>::constantState(const BlackoilState& x,
const WellStateFullyImplicitBlackoil& xw)
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
const int np = x.numPhases();
auto state = variableState(x, xw);
// The block pattern assumes the following primary variables:
// pressure
// water saturation (if water present)
// gas saturation, Rv (vapor oil/gas ratio) or Rs (solution gas/oil ratio) depending on hydrocarbon state
// Gas only (undersaturated gas): Rv
// Gas and oil: Sg
// Oil only (undersaturated oil): Rs
// well rates per active phase and well
// well bottom-hole pressure
// Note that oil is assumed to always be present, but is never
// a primary variable.
assert(active_[ Oil ]);
std::vector<int> bpat(np, nc);
bpat.push_back(xw.bhp().size() * np);
bpat.push_back(xw.bhp().size());
SolutionState state(np);
// Pressure.
assert (not x.pressure().empty());
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
state.pressure = ADB::constant(p, bpat);
// Saturation.
assert (not x.saturation().empty());
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
const Opm::PhaseUsage pu = fluid_.phaseUsage();
{
V so = V::Ones(nc, 1);
if (active_[ Water ]) {
const int pos = pu.phase_pos[ Water ];
const V sw = s.col(pos);
so -= sw;
state.saturation[pos] = ADB::constant(sw, bpat);
}
if (active_[ Gas ]) {
const int pos = pu.phase_pos[ Gas ];
const V sg = s.col(pos);
so -= sg;
state.saturation[pos] = ADB::constant(sg, bpat);
}
if (active_[ Oil ]) {
const int pos = pu.phase_pos[ Oil ];
state.saturation[pos] = ADB::constant(so, bpat);
}
}
// Solution Gas-oil ratio (rs).
if (active_[ Oil ] && active_[ Gas ]) {
const V rs = Eigen::Map<const V>(& x.gasoilratio()[0], x.gasoilratio().size());
state.rs = ADB::constant(rs, bpat);
} else {
const V Rs = V::Zero(nc, 1);
state.rs = ADB::constant(Rs, bpat);
}
// Vapor Oil-gas ratio (rv).
if (active_[ Oil ] && active_[ Gas ]) {
const V rv = Eigen::Map<const V>(& x.rv()[0], x.rv().size());
state.rv = ADB::constant(rv, bpat);
} else {
const V rv = V::Zero(nc, 1);
state.rv = ADB::constant(rv, bpat);
}
// Well rates.
assert (not xw.wellRates().empty());
// Need to reshuffle well rates, from ordered by wells, then phase,
// to ordered by phase, then wells.
const int nw = wells_.number_of_wells;
// The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
const V qs = Eigen::Map<const V>(wrates.data(), nw*np);
state.qs = ADB::constant(qs, bpat);
// Well bottom-hole pressure.
assert (not xw.bhp().empty());
const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
state.bhp = ADB::constant(bhp, bpat);
// HACK: throw away the derivatives. this may not be the most
// performant way to do things, but it will make the state
// automatically consistent with variableState() (and doing
// things automatically is all the rage in this module ;)
state.pressure = ADB::constant(state.pressure.value());
state.rs = ADB::constant(state.rs.value());
state.rv = ADB::constant(state.rv.value());
for (int phaseIdx= 0; phaseIdx < x.numPhases(); ++ phaseIdx)
state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value());
state.qs = ADB::constant(state.qs.value());
state.bhp = ADB::constant(state.bhp.value());
return state;
}
@@ -629,11 +559,81 @@ namespace {
template<class T>
void FullyImplicitBlackoilSolver<T>::computeWellConnectionPressures(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw)
{
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects.
const int nperf = wells_.well_connpos[wells_.number_of_wells];
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
// Compute b, rsmax, rvmax values for perforations.
const ADB perf_press = subset(state.pressure, well_cells);
std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid_.phaseUsage();
DataBlock b(nperf, pu.num_phases);
std::vector<double> rssat_perf(nperf, 0.0);
std::vector<double> rvsat_perf(nperf, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const ADB bw = fluid_.bWat(perf_press, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value();
const V rssat = fluidRsSat(perf_press.value(), well_cells);
rssat_perf.assign(rssat.data(), rssat.data() + nperf);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells);
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value();
const V rvsat = fluidRvSat(perf_press.value(), well_cells);
rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
}
// b is row major, so can just copy data.
std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
// Extract well connection depths.
const V depth = Eigen::Map<DataBlock>(grid_.cell_centroids, grid_.number_of_cells, grid_.dimensions).rightCols<1>();
const V pdepth = subset(depth, well_cells);
std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
// Surface density.
std::vector<double> surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
// Gravity
double grav = 0.0;
const double* g = geo_.gravity();
const int dim = grid_.dimensions;
if (g) {
// Guard against gravity in anything but last dimension.
for (int dd = 0; dd < dim - 1; ++dd) {
assert(g[dd] == 0.0);
}
grav = g[dim - 1];
}
// 2. Compute pressure deltas, and store the results.
std::vector<double> cdp = WellDensitySegmented
::computeConnectionPressureDelta(wells_, xw, fluid_.phaseUsage(),
b_perf, rssat_perf, rvsat_perf, perf_depth,
surf_dens, grav);
well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
}
template <class T>
void
FullyImplicitBlackoilSolver<T>::
assemble(const V& pvdt,
const BlackoilState& x ,
const WellStateFullyImplicitBlackoil& xw )
WellStateFullyImplicitBlackoil& xw )
{
using namespace Opm::AutoDiffGrid;
// Create the primary variables.
@@ -694,6 +694,367 @@ namespace {
}
addWellEq(state, xw); // Note: this can change xw (switching well controls).
addWellControlEq(state, xw);
}
void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state,
WellStateFullyImplicitBlackoil& xw)
{
const int nc = grid_.number_of_cells;
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
const int nperf = wells_.well_connpos[nw];
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
V Tw = Eigen::Map<const V>(wells_.WI, nperf);
const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
// pressure diffs computed already (once per step, not changing per iteration)
const V& cdp = well_perforation_pressure_diffs_;
// Extract variables for perforation cell pressures
// and corresponding perforation well pressures.
const ADB p_perfcell = subset(state.pressure, well_cells);
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcell - (wops_.w2p * state.bhp + cdp);
// current injecting connections
auto connInjInx = drawdown.value() < 0;
// injector == 1, producer == 0
V isInj = V::Zero(nw);
for (int w = 0; w < nw; ++w) {
if (wells_.type[w] == INJECTOR) {
isInj[w] = 1;
}
}
// // A cross-flow connection is defined as a connection which has opposite
// // flow-direction to the well total flow
// V isInjPerf = (wops_.w2p * isInj);
// auto crossFlowConns = (connInjInx != isInjPerf);
// bool allowCrossFlow = true;
// if (not allowCrossFlow) {
// auto closedConns = crossFlowConns;
// for (int c = 0; c < nperf; ++c) {
// if (closedConns[c]) {
// Tw[c] = 0;
// }
// }
// connInjInx = !closedConns;
// }
// TODO: not allow for crossflow
V isInjInx = V::Zero(nperf);
V isNotInjInx = V::Zero(nperf);
for (int c = 0; c < nperf; ++c){
if (connInjInx[c])
isInjInx[c] = 1;
else
isNotInjInx[c] = 1;
}
// HANDLE FLOW INTO WELLBORE
// compute phase volumerates standard conditions
std::vector<ADB> cq_ps(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const ADB& wellcell_mob = subset ( rq_[phase].mob, well_cells);
const ADB cq_p = -(isNotInjInx * Tw) * (wellcell_mob * drawdown);
cq_ps[phase] = subset(rq_[phase].b,well_cells) * cq_p;
}
if (active_[Oil] && active_[Gas]) {
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
ADB cq_psOil = cq_ps[oilpos];
ADB cq_psGas = cq_ps[gaspos];
cq_ps[gaspos] += subset(state.rs,well_cells) * cq_psOil;
cq_ps[oilpos] += subset(state.rv,well_cells) * cq_psGas;
}
// phase rates at std. condtions
std::vector<ADB> q_ps(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
q_ps[phase] = wops_.p2w * cq_ps[phase];
}
// total rates at std
ADB qt_s = ADB::constant(V::Zero(nw), state.bhp.blockPattern());
for (int phase = 0; phase < np; ++phase) {
qt_s += subset(state.qs, Span(nw, 1, phase*nw));
}
// compute avg. and total wellbore phase volumetric rates at std. conds
const DataBlock compi = Eigen::Map<const DataBlock>(wells_.comp_frac, nw, np);
std::vector<ADB> wbq(np, ADB::null());
ADB wbqt = ADB::constant(V::Zero(nw), state.pressure.blockPattern());
for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase];
wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase];
wbqt += wbq[phase];
}
// check for dead wells
V isNotDeadWells = V::Constant(nw,1.0);
for (int w = 0; w < nw; ++w) {
if (wbqt.value()[w] == 0) {
isNotDeadWells[w] = 0;
}
}
// compute wellbore mixture at std conds
Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
std::vector<ADB> mix_s(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase];
mix_s[phase] = notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
}
// HANDLE FLOW OUT FROM WELLBORE
// Total mobilities
ADB mt = subset(rq_[0].mob,well_cells);
for (int phase = 1; phase < np; ++phase) {
mt += subset(rq_[phase].mob,well_cells);
}
// injection connections total volumerates
ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown);
// compute volume ratio between connection at standard conditions
ADB volRat = ADB::constant(V::Zero(nperf), state.pressure.blockPattern());
std::vector<ADB> cmix_s(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
cmix_s[phase] = wops_.w2p * mix_s[phase];
}
ADB well_rv = subset(state.rv,well_cells);
ADB well_rs = subset(state.rs,well_cells);
ADB d = V::Constant(nperf,1.0) - well_rv * well_rs;
for (int phase = 0; phase < np; ++phase) {
ADB tmp = cmix_s[phase];
if (phase == Oil && active_[Gas]) {
const int gaspos = pu.phase_pos[Gas];
tmp = tmp - subset(state.rv,well_cells) * cmix_s[gaspos] / d;
}
if (phase == Gas && active_[Oil]) {
const int oilpos = pu.phase_pos[Oil];
tmp = tmp - subset(state.rs,well_cells) * cmix_s[oilpos] / d;
}
volRat += tmp / subset(rq_[phase].b,well_cells);
}
// injecting connections total volumerates at std cond
ADB cqt_is = cqt_i/volRat;
// connection phase volumerates at std cond
std::vector<ADB> cq_s(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
cq_s[phase] = cq_ps[phase] + (wops_.w2p * mix_s[phase])*cqt_is;
}
// Add well contributions to mass balance equations
for (int phase = 0; phase < np; ++phase) {
residual_.mass_balance[phase] -= superset(cq_s[phase],well_cells,nc);
}
// Add WELL EQUATIONS
ADB qs = state.qs;
for (int phase = 0; phase < np; ++phase) {
qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
}
residual_.well_flux_eq = qs;
// Updating the well controls is done from here, since we have
// access to the necessary bhp and rate data in this method.
updateWellControls(state.bhp, state.qs, xw);
}
namespace
{
double rateToCompare(const ADB& well_phase_flow_rate,
const int well,
const int num_phases,
const double* distr)
{
const int num_wells = well_phase_flow_rate.size() / num_phases;
double rate = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
// Important: well_phase_flow_rate is ordered with all rates for first
// phase coming first, then all for second phase etc.
rate += well_phase_flow_rate.value()[well + phase*num_wells] * distr[phase];
}
return rate;
}
bool constraintBroken(const ADB& bhp,
const ADB& well_phase_flow_rate,
const int well,
const int num_phases,
const WellType& well_type,
const WellControls* wc,
const int ctrl_index)
{
const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
const double target = well_controls_iget_target(wc, ctrl_index);
const double* distr = well_controls_iget_distr(wc, ctrl_index);
switch (well_type) {
case INJECTOR:
switch (ctrl_type) {
case BHP:
return bhp.value()[well] > target;
case SURFACE_RATE:
return rateToCompare(well_phase_flow_rate, well, num_phases, distr) > target;
case RESERVOIR_RATE:
// Intentional fall-through.
default:
OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls.");
}
break;
case PRODUCER:
switch (ctrl_type) {
case BHP:
return bhp.value()[well] < target;
case SURFACE_RATE:
// Note that the rates compared below are negative,
// so breaking the constraints means: too high flow rate
// (as for injection).
return rateToCompare(well_phase_flow_rate, well, num_phases, distr) < target;
case RESERVOIR_RATE:
// Intentional fall-through.
default:
OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls.");
}
break;
default:
OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
}
}
} // anonymous namespace
void FullyImplicitBlackoilSolver::updateWellControls(const ADB& bhp,
const ADB& well_phase_flow_rate,
WellStateFullyImplicitBlackoil& xw) const
{
std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = xw.currentControls()[w];
// Loop over all controls except the current one, and also
// skip any RESERVOIR_RATE controls, since we cannot
// handle those.
const int nwc = well_controls_get_num(wc);
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (ctrl_index == current) {
// This is the currently used control, so it is
// used as an equation. So this is not used as an
// inequality constraint, and therefore skipped.
continue;
}
if (well_controls_iget_type(wc, ctrl_index) == RESERVOIR_RATE) {
// We cannot handle this yet.
#ifdef OPM_VERBOSE
std::cout << "Warning: a RESERVOIR_RATE well control exists for well "
<< wells_.name[w] << ", but will never be checked." << std::endl;
#endif
continue;
}
if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells_.type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
}
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
std::cout << "Switching control mode for well " << wells_.name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
xw.currentControls()[w] = ctrl_index;
}
}
}
void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw)
{
// Handling BHP and SURFACE_RATE wells.
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
V bhp_targets(nw);
V rate_targets(nw);
M rate_distr(nw, np*nw);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = xw.currentControls()[w];
if (well_controls_iget_type(wc, current) == BHP) {
bhp_targets[w] = well_controls_iget_target(wc, current);
rate_targets[w] = -1e100;
} else if (well_controls_iget_type(wc, current) == SURFACE_RATE) {
bhp_targets[w] = -1e100;
rate_targets[w] = well_controls_iget_target(wc, current);
{
const double * distr = well_controls_iget_distr(wc, current);
for (int phase = 0; phase < np; ++phase) {
rate_distr.insert(w, phase*nw + w) = distr[phase];
}
}
} else {
OPM_THROW(std::runtime_error, "Can only handle BHP and SURFACE_RATE type controls.");
}
}
const ADB bhp_residual = state.bhp - bhp_targets;
const ADB rate_residual = rate_distr * state.qs - rate_targets;
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
// DUMP(residual_.well_eq);
}
void FullyImplicitBlackoilSolver::addOldWellEq(const SolutionState& state)
{
// -------- Well equation, and well contributions to the mass balance equations --------
// Contribution to mass balance will have to wait.
@@ -807,35 +1168,6 @@ namespace {
// Set the well flux equation
residual_.well_flux_eq = state.qs + well_rates_all;
// DUMP(residual_.well_flux_eq);
// Handling BHP and SURFACE_RATE wells.
V bhp_targets(nw);
V rate_targets(nw);
M rate_distr(nw, np*nw);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
if (well_controls_get_current_type(wc) == BHP) {
bhp_targets[w] = well_controls_get_current_target(wc);
rate_targets[w] = -1e100;
} else if (well_controls_get_current_type( wc ) == SURFACE_RATE) {
bhp_targets[w] = -1e100;
rate_targets[w] = well_controls_get_current_target(wc);
{
const double * distr = well_controls_get_current_distr( wc );
for (int phase = 0; phase < np; ++phase) {
rate_distr.insert(w, phase*nw + w) = distr[phase];
}
}
} else {
OPM_THROW(std::runtime_error, "Can only handle BHP and SURFACE_RATE type controls.");
}
}
const ADB bhp_residual = bhp - bhp_targets;
const ADB rate_residual = rate_distr * state.qs - rate_targets;
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
// DUMP(residual_.well_eq);
}
@@ -915,7 +1247,7 @@ namespace {
V isSg = V::Zero(nc,1);
bool disgas = false;
bool vapoil = false;
bool vapoil = false;
// this is a temporary hack to find if vapoil or disgas
// is a active component. Should be given directly from
@@ -1294,18 +1626,22 @@ namespace {
double
FullyImplicitBlackoilSolver<T>::residualNorm() const
{
double r = 0;
for (std::vector<ADB>::const_iterator
b = residual_.mass_balance.begin(),
e = residual_.mass_balance.end();
b != e; ++b)
double globalNorm = 0;
std::vector<ADB>::const_iterator quantityIt = residual_.mass_balance.begin();
const std::vector<ADB>::const_iterator endQuantityIt = residual_.mass_balance.end();
for (; quantityIt != endQuantityIt; ++quantityIt)
{
r = std::max(r, (*b).value().matrix().norm());
const double quantityResid = (*quantityIt).value().matrix().norm();
if (!std::isfinite(quantityResid)) {
OPM_THROW(Opm::NumericalProblem,
"Encountered a non-finite residual");
}
globalNorm = std::max(globalNorm, quantityResid);
}
r = std::max(r, residual_.well_flux_eq.value().matrix().norm());
r = std::max(r, residual_.well_eq.value().matrix().norm());
globalNorm = std::max(globalNorm, residual_.well_flux_eq.value().matrix().norm());
globalNorm = std::max(globalNorm, residual_.well_eq.value().matrix().norm());
return r;
return globalNorm;
}

4
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@@ -31,7 +31,6 @@ namespace Opm
{
namespace parameter { class ParameterGroup; }
class BlackoilPropsAdInterface;
class EclipseWriter;
class RockCompressibility;
class WellsManager;
class LinearSolverInterface;
@@ -75,8 +74,7 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity,
EclipseWriter &writer);
const double* gravity);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will

25
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@@ -33,7 +33,6 @@
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
@@ -67,8 +66,7 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity,
EclipseWriter &writer);
const double* gravity);
SimulatorReport run(SimulatorTimer& timer,
BlackoilState& state,
@@ -97,7 +95,6 @@ namespace Opm
FullyImplicitBlackoilSolver<Grid> solver_;
// Misc. data
std::vector<int> allcells_;
EclipseWriter &eclipseWriter_;
};
@@ -110,11 +107,10 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity,
EclipseWriter &eclipseWriter)
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity, eclipseWriter));
pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity));
}
@@ -196,8 +192,7 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
LinearSolverInterface& linsolver,
const double* gravity,
EclipseWriter &eclipseWriter)
const double* gravity)
: grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
@@ -205,9 +200,7 @@ namespace Opm
wells_(wells_manager.c_wells()),
gravity_(gravity),
geo_(grid_, props_, gravity_),
solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver),
eclipseWriter_(eclipseWriter)
solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver)
/* param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
@@ -350,12 +343,8 @@ namespace Opm
}
// advance to next timestep before reporting at this location
++timer;
// write an output file for later inspection
if (output_) {
eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState());
}
// ++timer; // Commented out since this has temporarily moved to the main() function.
break; // this is a temporary measure
}
total_timer.stop();

25
opm/autodiff/WellDensitySegmented.cpp
View File

@@ -19,22 +19,23 @@
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/core/wells.h>
#include <opm/core/simulator/WellState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <numeric>
#include <cmath>
std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
const WellState& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& z_perf,
const std::vector<double>& surf_dens,
const double gravity)
std::vector<double>
Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& z_perf,
const std::vector<double>& surf_dens,
const double gravity)
{
// Verify that we have consistent input.
const int np = wells.number_of_phases;
@@ -46,7 +47,7 @@ std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(co
if (surf_dens.size() != size_t(wells.number_of_phases)) {
OPM_THROW(std::logic_error, "Inconsistent input: surf_dens vs. phase_usage.");
}
if (nperf*np != int(wstate.perfRates().size())) {
if (nperf*np != int(wstate.perfPhaseRates().size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate.");
}
if (nperf*np != int(b_perf.size())) {
@@ -90,7 +91,7 @@ std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(co
q_out_perf[perf*np + phase] = q_out_perf[(perf+1)*np + phase];
}
// Subtract outflow through perforation.
q_out_perf[perf*np + phase] -= wstate.perfRates()[perf*np + phase];
q_out_perf[perf*np + phase] -= wstate.perfPhaseRates()[perf*np + phase];
}
}
}

4
opm/autodiff/WellDensitySegmented.hpp
View File

@@ -27,7 +27,7 @@ struct Wells;
namespace Opm
{
class WellState;
class WellStateFullyImplicitBlackoil;
struct PhaseUsage;
@@ -51,7 +51,7 @@ namespace Opm
/// \param[in] surf_dens surface densities for active components, size P
/// \param[in] gravity gravity acceleration constant
static std::vector<double> computeConnectionPressureDelta(const Wells& wells,
const WellState& wstate,
const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,

13
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@@ -70,6 +70,14 @@ namespace Opm
}
}
}
// Initialize current_controls_.
// The controls set in the Wells object are treated as defaults,
// and also used for initial values.
current_controls_.resize(nw);
for (int w = 0; w < nw; ++w) {
current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
}
}
/// One bhp pressure per well.
@@ -84,6 +92,10 @@ namespace Opm
std::vector<double>& perfPhaseRates() { return perfphaserates_; }
const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
/// One current control per well.
std::vector<int>& currentControls() { return current_controls_; }
const std::vector<int>& currentControls() const { return current_controls_; }
/// For interfacing with functions that take a WellState.
const WellState& basicWellState() const
{
@@ -93,6 +105,7 @@ namespace Opm
private:
WellState basic_well_state_;
std::vector<double> perfphaserates_;
std::vector<int> current_controls_;
};
} // namespace Opm