removing the use of wops_ from the BlackoilModelBase in ms well

Also some cleaning up of comments in BlackoilMultiSegmentModel
This commit is contained in:
Kai Bao 2015-11-20 13:14:21 +01:00
parent 61cccfebf8
commit f90b96abf6
2 changed files with 47 additions and 94 deletions

View File

@ -64,18 +64,19 @@ namespace Opm {
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
/// \param[in] param parameters
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] vfp_properties Vertical flow performance tables
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
/// \param[in] param parameters
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] vfp_properties Vertical flow performance tables
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
/// \param[in] wells_multisegment a vector of multisegment wells
BlackoilMultiSegmentModel(const typename Base::ModelParameters& param,
const Grid& grid ,
const BlackoilPropsAdInterface& fluid,
@ -111,16 +112,10 @@ namespace Opm {
using Base::materialName;
protected:
/*
// --------- Types and enums ---------
// using Base::DataBlock;
// using Base::ReservoirResidualQuant;
*/
// --------- Data members ---------
// For the non-segmented well, it should be the density with AVG or SEG way.
// while usually SEG way
using Base::wops_; // Only for well_cells, the w2p and p2w members may have wrong perf order.
// For non-segmented wells, it should be the density calculated with AVG or SEG way.
// while usually SEG way by default.
using Base::well_perforation_densities_; //Density of each well perforation
using Base::pvdt_;
using Base::geo_;
@ -141,37 +136,21 @@ namespace Opm {
using Base::param_;
using Base::linsolver_;
// Diff to the pressure of the related segment.
// When the well is a usual well, the bhp will be the pressure of the top segment
// For mutlti-segmented wells, only AVG is allowed.
// For non-segmented wells, typically SEG is used. AVG way might not have been
// implemented yet.
// Diff to bhp for each well perforation. only for usual wells.
// For segmented wells, they are zeros.
using Base::well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
// ADB version of the densities, when using AVG way, the calculation of the density and hydrostatic head
// is implicit
// ADB well_perforation_densities_adb_; // TODO: NOT NEEDED for explicit SEG way
// ADB version. Eventually, only ADB version will be kept.
// ADB well_perforation_pressure_diffs_adb_; // TODO: NOT NEEDED for explicit SEG way
using Base::well_perforation_pressure_diffs_;
// Pressure correction due to the different depth of the perforation
// and the cell center of the grid block
// For the non-segmented wells, since the perforation are forced to be
// at the center of the grid cell, it should be ZERO.
// It should only apply to the mutli-segmented wells.
// It only applies to the mutli-segmented wells.
V well_perforation_cell_pressure_diffs_;
// Pressure correction due to the depth differennce between segment depth and perforation depth.
// TODO: It will only be able to be merge as a part of the perforation_pressure_diffs_
ADB well_segment_perforation_pressure_diffs_;
// The depth difference between segment nodes and perforations
// TODO: it should be a member in a global wells class later
V well_segment_perforation_depth_diffs_;
// the average of the fluid densities in the grid block
@ -179,7 +158,6 @@ namespace Opm {
// and the cell center of the grid block
V well_perforation_cell_densities_;
// the density of the fluid mixture in the segments
// which is calculated in an implicit way
ADB well_segment_densities_;
@ -190,9 +168,10 @@ namespace Opm {
ADB well_segment_pressures_delta_;
// the surface volume of components in the segments
// initial one at the beginning of the time step
// the initial value at the beginning of the time step
std::vector<V> segment_comp_surf_volume_initial_;
// the current one at the current iteration.
// the value within the current iteration.
std::vector<ADB> segment_comp_surf_volume_current_;
// the mass flow rate in the segments
@ -200,7 +179,7 @@ namespace Opm {
// the viscosity of the fluid mixture in the segments
// TODO: it is only used to calculate the Reynolds number as we know
// maybe it is not better just to store the Reynolds number here?
// maybe it is not better just to store the Reynolds number?
ADB segment_viscosities_;
const std::vector<WellMultiSegmentConstPtr> wells_multisegment_;

View File

@ -98,22 +98,6 @@ namespace Opm {
, wells_multisegment_(wells_multisegment)
, wops_ms_(wells_multisegment)
{
// Modify the wops_.well_cell member, since the
// order of perforations is different for multi-segment
// wells, since the segments may have been reordered.
// It should already have the correct capacity, though.
// Note that we need the const_cast since wops_ is declared
// as a 'const WellOps'.
std::vector<int>& well_cells = const_cast<typename Base::WellOps&>(wops_).well_cells;
const size_t old_sz = well_cells.size();
well_cells.clear();
const int nw = wells_multisegment.size();
for (int i = 0; i < nw; ++i) {
const std::vector<int>& temp_well_cells = wellsMultiSegment()[i]->wellCells();
well_cells.insert(well_cells.end(), temp_well_cells.begin(), temp_well_cells.end());
}
assert(old_sz == well_cells.size());
static_cast<void>(old_sz); // Avoid warning in release mode.
}
@ -152,58 +136,36 @@ namespace Opm {
assert(well_perf_start == total_nperf);
assert(int(well_cells.size()) == total_nperf);
// Create the segment <-> perforation operator matrices,
// Create all the operator matrices,
// using the setFromTriplets() method.
s2p = Eigen::SparseMatrix<double>(total_nperf, total_nseg);
p2s = Eigen::SparseMatrix<double>(total_nseg, total_nperf);
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> scatter, gather;
scatter.reserve(total_nperf);
gather .reserve(total_nperf);
int seg_start = 0;
int perf_start = 0;
for (int w = 0; w < nw; ++w) {
const int ns = wells_ms[w]->numberOfSegments();
const int np = wells_ms[w]->numberOfPerforations();
for (int seg = 0; seg < ns; ++seg) {
const auto& segPerf = wells_ms[w]->segmentPerforations()[seg];
// the number of perforations related to this segment
const int npseg = segPerf.size();
for (int perf = 0; perf < npseg; ++perf) {
const int perf_ind = perf_start + segPerf[perf];
const int seg_ind = seg_start + seg;
scatter.push_back(Tri(perf_ind, seg_ind, 1.0));
gather .push_back(Tri(seg_ind, perf_ind, 1.0));
}
}
perf_start += np;
seg_start += ns;
}
s2p.setFromTriplets(scatter.begin(), scatter.end());
p2s.setFromTriplets(gather .begin(), gather .end());
// Crate the segment -> its inlet segments and segment-> its outlet segment operator matrices,
// It can be done togehter with other operator generation process.
// Just do seperately for the moment for parallel development.
s2s_inlets = Eigen::SparseMatrix<double>(total_nseg, total_nseg);
s2s_outlet = Eigen::SparseMatrix<double>(total_nseg, total_nseg);
s2w = Eigen::SparseMatrix<double>(nw, total_nseg);
w2s = Eigen::SparseMatrix<double>(total_nseg, nw);
topseg2w = Eigen::SparseMatrix<double>(nw, total_nseg);
s2p = Eigen::SparseMatrix<double>(total_nperf, total_nseg);
p2s = Eigen::SparseMatrix<double>(total_nseg, total_nperf);
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> s2s_inlets_vector;
std::vector<Tri> s2s_outlet_vector;
std::vector<Tri> s2w_vector;
std::vector<Tri> w2s_vector;
std::vector<Tri> topseg2w_vector;
std::vector<Tri> s2p_vector;
std::vector<Tri> p2s_vector;
V topseg_zero = V::Ones(total_nseg);
s2s_inlets_vector.reserve(total_nseg);
s2s_outlet_vector.reserve(total_nseg);
s2w_vector.reserve(total_nseg);
w2s_vector.reserve(total_nseg);
seg_start = 0;
topseg2w_vector.reserve(nw);
s2p_vector.reserve(total_nperf);
p2s_vector.reserve(total_nperf);
int seg_start = 0;
int perf_start = 0;
for (int w = 0; w < nw; ++w) {
const int ns = wells_ms[w]->numberOfSegments();
const int np = wells_ms[w]->numberOfPerforations();
for (int seg = 0; seg < ns; ++seg) {
const int seg_ind = seg_start + seg;
w2s_vector.push_back(Tri(seg_ind, w, 1.0));
@ -218,8 +180,18 @@ namespace Opm {
s2s_inlets_vector.push_back(Tri(outlet_ind, seg_ind, 1.0));
s2s_outlet_vector.push_back(Tri(seg_ind, outlet_ind, 1.0));
}
const auto& seg_perf = wells_ms[w]->segmentPerforations()[seg];
// the number of perforations related to this segment
const int npseg = seg_perf.size();
for (int perf = 0; perf < npseg; ++perf) {
const int perf_ind = perf_start + seg_perf[perf];
s2p_vector.push_back(Tri(perf_ind, seg_ind, 1.0));
p2s_vector.push_back(Tri(seg_ind, perf_ind, 1.0));
}
}
seg_start += ns;
perf_start += np;
}
s2s_inlets.setFromTriplets(s2s_inlets_vector.begin(), s2s_inlets_vector.end());
@ -227,6 +199,8 @@ namespace Opm {
w2s.setFromTriplets(w2s_vector.begin(), w2s_vector.end());
s2w.setFromTriplets(s2w_vector.begin(), s2w_vector.end());
topseg2w.setFromTriplets(topseg2w_vector.begin(), topseg2w_vector.end());
s2p.setFromTriplets(s2p_vector.begin(), s2p_vector.end());
p2s.setFromTriplets(p2s_vector.begin(), p2s_vector.end());
w2p = Eigen::SparseMatrix<double>(total_nperf, nw);
p2w = Eigen::SparseMatrix<double>(nw, total_nperf);
@ -808,7 +782,7 @@ namespace Opm {
is_multisegment_well[w] = double(wellsMultiSegment()[w]->isMultiSegmented());
}
// Take one flag per well and expand to one flag per perforation.
V is_multisegment_perf = wops_.w2p * is_multisegment_well.matrix();
V is_multisegment_perf = wops_ms_.w2p * is_multisegment_well.matrix();
Selector<double> msperf_selector(is_multisegment_perf, Selector<double>::NotEqualZero);
// Compute drawdown.
@ -937,7 +911,7 @@ namespace Opm {
for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase];
const ADB phase_fraction = wops_ms_.topseg2w * (wbq[phase] / wbqt);
cmix_s[phase] = wops_.w2p * aliveWells_selector.select(phase_fraction, ADB::constant(compi.col(pos)));
cmix_s[phase] = wops_ms_.w2p * aliveWells_selector.select(phase_fraction, ADB::constant(compi.col(pos)));
}
// compute volume ration between connection at standard conditions