adding computePropertiesForWellConnectionPressures to StandardWells

tested with flow. Not sure how other flow siblings work.
This commit is contained in:
Kai Bao 2016-04-07 16:17:25 +02:00
parent 3bcfc905bd
commit 266f6f2df5
3 changed files with 109 additions and 1 deletions

View File

@ -849,7 +849,7 @@ namespace detail {
std::vector<double> rsmax_perf;
std::vector<double> rvmax_perf;
std::vector<double> surf_dens_perf;
asImpl().computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// 2. Compute densities
std::vector<double> cd =

View File

@ -31,6 +31,8 @@
#include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
namespace Opm {
@ -39,6 +41,13 @@ namespace Opm {
typedef AutoDiffBlock<double> ADB;
typedef ADB::V Vector;
// copied from BlackoilModelBase
// should put to somewhere better
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
/// Class for handling the standard well model.
template <class SolutionState, class WellState>
class StandardWells {
@ -72,6 +81,17 @@ namespace Opm {
Vector& wellPerforationPressureDiffs();
const Vector& wellPerforationPressureDiffs() const;
void computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& pc,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf);
protected:
bool wells_active_;
const Wells* wells_;

View File

@ -23,6 +23,7 @@
namespace Opm
{
@ -185,4 +186,91 @@ namespace Opm
return well_perforation_pressure_diffs_;
}
template<class SolutionState, class WellState>
void
StandardWells<SolutionState, WellState>::
computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const BlackoilPropsAdInterface& fluid,
const std::vector<bool>& active,
const std::vector<PhasePresence>& pc,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf)
{
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
// Compute the average pressure in each well block
const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf);
Vector avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
const std::vector<int>& well_cells = wellOps().well_cells;
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
const ADB avg_press_ad = ADB::constant(avg_press);
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);
if (pu.phase_used[BlackoilPhases::Aqua]) {
const Vector bw = fluid.bWat(avg_press_ad, perf_temp, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
}
assert(active[Oil]);
const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells);
const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
// const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
} else {
rsmax_perf.assign(nperf, 0.0);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
const Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
// const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
} else {
rvmax_perf.assign(nperf, 0.0);
}
// b is row major, so can just copy data.
b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
// Surface density.
// The compute density segment wants the surface densities as
// an np * number of wells cells array
Vector rho = superset(fluid.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
rho += superset(fluid.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
}
surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases);
}
}