OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-28 18:21:00 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1346 from GitPaean/removing_more_multisegment_legacy

Browse Source 
removing more legacy Multisegment Wells related.
This commit is contained in:
Atgeirr Flø Rasmussen 2017-11-24 15:30:39 +01:00 committed by GitHub
commit 3fd27d9bd6
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 0 additions and 3135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists_files.cmake
View File

@ -46,8 +46,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/VFPProperties.cpp
opm/autodiff/VFPProdProperties.cpp
opm/autodiff/VFPInjProperties.cpp
opm/autodiff/WellMultiSegment.cpp
opm/autodiff/MultisegmentWells.cpp
opm/autodiff/MissingFeatures.cpp
opm/polymer/PolymerState.cpp
opm/polymer/PolymerBlackoilState.cpp
@ -92,7 +90,6 @@ list (APPEND TEST_SOURCE_FILES
tests/test_welldensitysegmented.cpp
tests/test_vfpproperties.cpp
tests/test_singlecellsolves.cpp
tests/test_multisegmentwells.cpp
tests/test_multiphaseupwind.cpp
tests/test_wellmodel.cpp
# tests/test_thresholdpressure.cpp
@ -213,10 +210,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/VFPHelpers.hpp
opm/autodiff/VFPProdProperties.hpp
opm/autodiff/VFPInjProperties.hpp
opm/autodiff/WellStateMultiSegment.hpp
opm/autodiff/WellMultiSegment.hpp
opm/autodiff/MultisegmentWells.hpp
opm/autodiff/MultisegmentWells_impl.hpp
opm/autodiff/WellHelpers.hpp
opm/autodiff/StandardWells.hpp
opm/autodiff/StandardWells_impl.hpp

442
opm/autodiff/MultisegmentWells.cpp
View File

@ -1,442 +0,0 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/autodiff/MultisegmentWells.hpp>
namespace Opm {
MultisegmentWells::
MultisegmentWellOps::MultisegmentWellOps(const std::vector<WellMultiSegmentConstPtr>& wells_ms)
{
// no multi-segment wells are involved by default.
has_multisegment_wells = false;
if (wells_ms.empty()) {
return;
}
// Count the total number of perforations and segments.
const int nw = wells_ms.size();
int total_nperf = 0;
int total_nseg = 0;
for (int w = 0; w < nw; ++w) {
total_nperf += wells_ms[w]->numberOfPerforations();
total_nseg += wells_ms[w]->numberOfSegments();
if (wells_ms[w]->isMultiSegmented()) {
has_multisegment_wells = true;
}
}
// Create well_cells and conn_trans_factors.
well_cells.reserve(total_nperf);
conn_trans_factors.resize(total_nperf);
int well_perf_start = 0;
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wells_ms[w];
well_cells.insert(well_cells.end(), well->wellCells().begin(), well->wellCells().end());
const std::vector<double>& perf_trans = well->wellIndex();
std::copy(perf_trans.begin(), perf_trans.end(), conn_trans_factors.data() + well_perf_start);
well_perf_start += well->numberOfPerforations();
}
assert(well_perf_start == total_nperf);
assert(int(well_cells.size()) == total_nperf);
// Create all the operator matrices,
// using the setFromTriplets() method.
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;
Vector topseg_zero = Vector::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);
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));
s2w_vector.push_back(Tri(w, seg_ind, 1.0));
if (seg == 0) {
topseg2w_vector.push_back(Tri(w, seg_ind, 1.0));
topseg_zero(seg_ind) = 0.0;
}
int seg_outlet = wells_ms[w]->outletSegment()[seg];
if (seg_outlet >= 0) {
const int outlet_ind = seg_start + seg_outlet;
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());
s2s_outlet.setFromTriplets(s2s_outlet_vector.begin(), s2s_outlet_vector.end());
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);
w2p = s2p * w2s;
p2w = s2w * p2s;
eliminate_topseg = AutoDiffMatrix(topseg_zero.matrix().asDiagonal());
}
MultisegmentWells::
MultisegmentWells(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const int time_step)
: wells_multisegment_( createMSWellVector(wells_arg, wells_ecl, time_step) )
, wops_ms_(wells_multisegment_)
, well_collection_(well_collection)
, well_perforation_efficiency_factors_(Vector::Ones(numWells()))
, num_phases_(wells_arg ? wells_arg->number_of_phases : 0)
, wells_(wells_arg)
, fluid_(nullptr)
, active_(nullptr)
, phase_condition_(nullptr)
, vfp_properties_(nullptr)
, well_segment_perforation_pressure_diffs_(ADB::null())
, well_segment_densities_(ADB::null())
, well_segment_pressures_delta_(ADB::null())
, segment_comp_surf_volume_initial_(num_phases_)
, segment_comp_surf_volume_current_(num_phases_, ADB::null())
, segment_mass_flow_rates_(ADB::null())
, segment_viscosities_(ADB::null())
{
const int nw = wells_multisegment_.size();
int nperf_total = 0;
int nseg_total = 0;
top_well_segments_.resize(nw);
for (int w = 0; w < nw; ++w) {
top_well_segments_[w] = nseg_total;
nperf_total += wells_multisegment_[w]->numberOfPerforations();
nseg_total += wells_multisegment_[w]->numberOfSegments();
}
nperf_total_ = nperf_total;
nseg_total_ = nseg_total;
}
std::vector<WellMultiSegmentConstPtr>
MultisegmentWells::createMSWellVector(const Wells* wells_arg,
const std::vector< const Well* >& wells_ecl,
const int time_step)
{
std::vector<WellMultiSegmentConstPtr> wells_multisegment;
if (wells_arg) {
// number of wells in wells_arg structure
const int nw = wells_arg->number_of_wells;
// number of wells in EclipseState
const int nw_ecl = wells_ecl.size();
wells_multisegment.reserve(nw);
for(int i = 0; i < nw_ecl; ++i) {
// not processing SHUT wells
if (wells_ecl[i]->getStatus(time_step) == WellCommon::SHUT) {
continue;
}
// checking if the well can be found in the wells
const std::string& well_name = wells_ecl[i]->name();
int index_well;
for (index_well = 0; index_well < nw; ++index_well) {
if (well_name == std::string(wells_arg->name[index_well])) {
break;
}
}
if (index_well != nw) { // found in the wells
wells_multisegment.push_back(std::make_shared<WellMultiSegment>(wells_ecl[i], time_step, wells_arg));
}
}
}
return wells_multisegment;
}
void
MultisegmentWells::init(const BlackoilPropsAdFromDeck* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg)
{
fluid_ = fluid_arg;
active_ = active_arg;
phase_condition_ = pc_arg;
vfp_properties_ = vfp_properties_arg;
gravity_ = gravity_arg;
perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);;
const int nw = wells_multisegment_.size();
// Calculating the depth difference between perforation and the cell center in the peforated cells.
std::vector<double> perf_depth_vec;
perf_depth_vec.reserve(nperf_total_);
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wells_multisegment_[w];
const std::vector<double>& perf_depth_well = well->perfDepth();
perf_depth_vec.insert(perf_depth_vec.end(), perf_depth_well.begin(), perf_depth_well.end());
}
assert(int(perf_depth_vec.size()) == nperf_total_);
const Vector perf_depth = Eigen::Map<Vector>(perf_depth_vec.data(), nperf_total_);
perf_cell_depth_diffs_ = perf_depth - perf_cell_depth_;
// Calculating the depth difference between segment nodes and perforations.
well_segment_perforation_depth_diffs_ = Vector::Constant(nperf_total_, -1e100);
int start_perforation = 0;
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = wells_multisegment_[w];
const int nseg = well->numberOfSegments();
const int nperf = well->numberOfPerforations();
const std::vector<std::vector<int>>& segment_perforations = well->segmentPerforations();
for (int s = 0; s < nseg; ++s) {
const int nperf_seg = segment_perforations[s].size();
const double segment_depth = well->segmentDepth()[s];
for (int perf = 0; perf < nperf_seg; ++perf) {
const int perf_number = segment_perforations[s][perf] + start_perforation;
well_segment_perforation_depth_diffs_[perf_number] = segment_depth - perf_depth[perf_number];
}
}
start_perforation += nperf;
}
assert(start_perforation == nperf_total_);
calculateEfficiencyFactors();
}
const std::vector<WellMultiSegmentConstPtr>&
MultisegmentWells::msWells() const
{
return wells_multisegment_;
}
const Wells&
MultisegmentWells::wells() const
{
assert(wells_ != nullptr);
return *(wells_);
}
const Wells*
MultisegmentWells::wellsPointer() const
{
return wells_;
}
const MultisegmentWells::MultisegmentWellOps&
MultisegmentWells::wellOps() const
{
return wops_ms_;
}
void
MultisegmentWells::
computeSegmentPressuresDelta(const double grav)
{
const int nw = msWells().size();
const int nseg_total = nseg_total_;
if ( !wellOps().has_multisegment_wells ) {
well_segment_pressures_delta_ = ADB::constant(Vector::Zero(nseg_total));
well_segment_perforation_pressure_diffs_ = wellOps().s2p * well_segment_pressures_delta_;
return;
}
// calculate the depth difference of the segments
// TODO: we need to store the following values somewhere to avoid recomputation.
Vector segment_depth_delta = Vector::Zero(nseg_total);
int start_segment = 0;
for (int w = 0; w < nw; ++w) {
WellMultiSegmentConstPtr well = msWells()[w];
const int nseg = well->numberOfSegments();
for (int s = 1; s < nseg; ++s) {
const int s_outlet = well->outletSegment()[s];
assert(s_outlet >= 0 && s_outlet < nseg);
segment_depth_delta[s + start_segment] = well->segmentDepth()[s_outlet] - well->segmentDepth()[s];
}
start_segment += nseg;
}
assert(start_segment == nseg_total);
const ADB grav_adb = ADB::constant(Vector::Constant(nseg_total, grav));
well_segment_pressures_delta_ = segment_depth_delta * grav_adb * well_segment_densities_;
ADB well_segment_perforation_densities = wellOps().s2p * well_segment_densities_;
well_segment_perforation_pressure_diffs_ = grav * well_segment_perforation_depth_diffs_ * well_segment_perforation_densities;
}
void
MultisegmentWells::
variableStateWellIndices(std::vector<int>& indices,
int& next) const
{
indices[Qs] = next++;
indices[Bhp] = next++;
}
std::vector<int>
MultisegmentWells::
variableWellStateIndices() const
{
// Black oil model standard is 5 equation.
// For the pure well solve, only the well equations are picked.
std::vector<int> indices(5, -1);
int next = 0;
variableStateWellIndices(indices, next);
assert(next == 2);
return indices;
}
WellCollection*
MultisegmentWells::
wellCollection() const {
return well_collection_;
}
void
MultisegmentWells::
calculateEfficiencyFactors()
{
if ( !localWellsActive() ) {
return;
}
// get efficiency factor for each well first
const int nw = wells_->number_of_wells;
Vector well_efficiency_factors = Vector::Ones(nw);
for (int w = 0; w < nw; ++w) {
const std::string well_name = wells_->name[w];
// get the well node in the well collection
WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
well_efficiency_factors(w) = well_node.getAccumulativeEfficiencyFactor();
}
// map them to the perforation.
well_perforation_efficiency_factors_ = wellOps().w2p * well_efficiency_factors.matrix();
}
const
MultisegmentWells::Vector&
MultisegmentWells::
wellPerfEfficiencyFactors() const
{
return well_perforation_efficiency_factors_;
}
} // end of namespace Opm

341
opm/autodiff/MultisegmentWells.hpp
View File

@ -1,341 +0,0 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MULTISEGMENTWELLS_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELLS_HEADER_INCLUDED
#include <dune/common/parallel/mpihelper.hh>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <cassert>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/wells/WellCollection.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/WellHelpers.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/WellMultiSegment.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/simulators/WellSwitchingLogger.hpp>
namespace Opm {
/// Class for handling the multi-segment well model
class MultisegmentWells {
public:
// --------- Types ---------
using ADB = AutoDiffBlock<double>;
using Vector = ADB::V;
// Well operations and data needed.
struct MultisegmentWellOps {
explicit MultisegmentWellOps(const std::vector<WellMultiSegmentConstPtr>& wells_ms);
Eigen::SparseMatrix<double> w2p; // well -> perf (scatter)
Eigen::SparseMatrix<double> p2w; // perf -> well (gather)
Eigen::SparseMatrix<double> w2s; // well -> segment (scatter)
Eigen::SparseMatrix<double> s2w; // segment -> well (gather)
Eigen::SparseMatrix<double> s2p; // segment -> perf (scatter)
Eigen::SparseMatrix<double> p2s; // perf -> segment (gather)
Eigen::SparseMatrix<double> s2s_inlets; // segment -> its inlet segments
Eigen::SparseMatrix<double> s2s_outlet; // segment -> its outlet segment
Eigen::SparseMatrix<double> topseg2w; // top segment -> well
AutoDiffMatrix eliminate_topseg; // change the top segment related to be zero
std::vector<int> well_cells; // the set of perforated cells
Vector conn_trans_factors; // connection transmissibility factors
bool has_multisegment_wells; // flag indicating whether there is any muli-segment well
};
// copied from BlackoilModelBase
// should put to somewhere better
using DataBlock = Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor>;
using Communication =
Dune::CollectiveCommunication<typename Dune::MPIHelper
::MPICommunicator>;
// --------- Public methods ---------
// TODO: using a vector of WellMultiSegmentConstPtr for now
// TODO: it should use const Wells or something else later.
MultisegmentWells(const Wells* wells_arg,
WellCollection* well_collection,
const std::vector< const Well* >& wells_ecl,
const int time_step);
std::vector<WellMultiSegmentConstPtr> createMSWellVector(const Wells* wells_arg,
const std::vector< const Well* >& wells_ecl,
const int time_step);
void init(const BlackoilPropsAdFromDeck* fluid_arg,
const std::vector<bool>* active_arg,
const std::vector<PhasePresence>* pc_arg,
const VFPProperties* vfp_properties_arg,
const double gravity_arg,
const Vector& depth_arg);
const std::vector<WellMultiSegmentConstPtr>& msWells() const;
const MultisegmentWellOps& wellOps() const;
const Wells& wells() const;
const Wells* wellsPointer() const;
int numPhases() const { return num_phases_; };
int numWells() const { return wells_multisegment_.size(); }
int numSegment() const { return nseg_total_; };
int numPerf() const { return nperf_total_; };
bool wellsActive() const { return wells_active_; };
void setWellsActive(const bool wells_active) { wells_active_ = wells_active; };
bool localWellsActive() const { return ! wells_multisegment_.empty(); };
int numWellVars() const { return (num_phases_ + 1) * nseg_total_; };
template <class ReservoirResidualQuant, class SolutionState>
void extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
Vector& wellPerforationCellPressureDiffs() { return well_perforation_cell_pressure_diffs_; };
Vector& wellSegmentPerforationDepthDiffs() { return well_segment_perforation_depth_diffs_; };
const Vector& wellPerforationCellDensities() const { return well_perforation_cell_densities_; };
Vector& wellPerforationCellDensities() { return well_perforation_cell_densities_; };
const std::vector<Vector>& segmentCompSurfVolumeInitial() const { return segment_comp_surf_volume_initial_; };
std::vector<Vector>& segmentCompSurfVolumeInitial() { return segment_comp_surf_volume_initial_; };
const std::vector<ADB>& segmentCompSurfVolumeCurrent() const { return segment_comp_surf_volume_current_; };
const std::vector<int>& topWellSegments() const { return top_well_segments_; };
std::vector<int>& topWellSegments() { return top_well_segments_; };
Vector& segVDt() { return segvdt_; };
const Vector& wellPerforationDensities() const { return well_perforation_densities_; };
Vector& wellPerforationDensities() { return well_perforation_densities_; };
const Vector& wellPerforationPressureDiffs() const { return well_perforation_pressure_diffs_; };
Vector& wellPerforationPressureDiffs() { return well_perforation_pressure_diffs_; };
template <class WellState>
void
updateWellState(const Vector& dwells,
const double dpmaxrel,
WellState& well_state) const;
// TODO: some arguments can be removed later
// TODO: compi will be required in the multisegment wells
template <class SolutionState>
void
computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
std::vector<ADB>& cq_s) const;
template <class SolutionState, class WellState>
void updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const;
// Calculate the density of the mixture in the segments
// And the surface volume of the components in the segments by dt
template <class SolutionState>
void
computeSegmentFluidProperties(const SolutionState& state);
void
computeSegmentPressuresDelta(const double grav);
template <class SolutionState>
void
addWellFluxEq(const std::vector<ADB>& cq_s,
const SolutionState& state,
LinearisedBlackoilResidual& residual);
template <class SolutionState, class WellState>
void
addWellControlEq(const SolutionState& state,
const WellState& xw,
const Vector& aliveWells,
LinearisedBlackoilResidual& residual);
template <class WellState>
void
updateWellControls(WellState& xw) const;
// TODO: these code are same with the StandardWells
// to find a better solution later.
void
variableStateWellIndices(std::vector<int>& indices,
int& next) const;
template <class SolutionState>
void
variableStateExtractWellsVars(const std::vector<int>& indices,
std::vector<ADB>& vars,
SolutionState& state) const;
std::vector<int>
variableWellStateIndices() const;
template <class WellState>
void
variableWellStateInitials(const WellState& xw,
std::vector<Vector>& vars0) const;
template <class SolutionState, class WellState>
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw,
const std::vector<ADB>& kr_adb,
const std::vector<ADB>& fluid_density);
WellCollection* wellCollection() const;
void calculateEfficiencyFactors();
const Vector& wellPerfEfficiencyFactors() const;
// just return the passed well state
template<class WellState>
const WellState& wellState(const WellState& well_state) const { return well_state; }
protected:
// TODO: probably a wells_active_ will be required here.
bool wells_active_;
std::vector<WellMultiSegmentConstPtr> wells_multisegment_;
MultisegmentWellOps wops_ms_;
// It will probably need to be updated during running time.
WellCollection* well_collection_;
// The efficiency factor for each connection
// It is specified based on wells and groups
// By default, they should all be one.
Vector well_perforation_efficiency_factors_;
const int num_phases_;
int nseg_total_;
int nperf_total_;
// TODO: put the Wells here to simplify the interface
// TODO: at the moment, both wells_ and wells_multisegment_
// TODO: acutally contain all the wells
// TODO: they should be split eventually.
const Wells* wells_;
const BlackoilPropsAdFromDeck* fluid_;
const std::vector<bool>* active_;
const std::vector<PhasePresence>* phase_condition_;
const VFPProperties* vfp_properties_;
double gravity_;
// The depth of the all the cell centers
// It is different from the perforation depth in MultisegmentWells
Vector perf_cell_depth_;
// 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 only applies to the mutli-segmented wells.
Vector well_perforation_cell_pressure_diffs_;
// Pressure correction due to the depth differennce between segment depth and perforation depth.
ADB well_segment_perforation_pressure_diffs_;
// The depth difference between segment nodes and perforations
Vector well_segment_perforation_depth_diffs_;
// The depth difference between the perforations and the perforation cells.
Vector perf_cell_depth_diffs_;
// the average of the fluid densities in the grid block
// which is used to calculate the hydrostatic head correction due to the depth difference of the perforation
// and the cell center of the grid block
Vector well_perforation_cell_densities_;
// the density of the fluid mixture in the segments
// which is calculated in an implicit way
ADB well_segment_densities_;
// the hydrostatic pressure drop between segment nodes
// calculated with the above density of fluid mixtures
// for the top segment, they should always be zero for the moment.
ADB well_segment_pressures_delta_;
// the surface volume of components in the segments
// the initial value at the beginning of the time step
std::vector<Vector> segment_comp_surf_volume_initial_;
// the value within the current iteration.
std::vector<ADB> segment_comp_surf_volume_current_;
// the mass flow rate in the segments
ADB segment_mass_flow_rates_;
// 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?
ADB segment_viscosities_;
std::vector<int> top_well_segments_;
// segment volume by dt (time step)
// to handle the volume effects of the segment
Vector segvdt_;
// technically, they are only useful for standard wells
// since at the moment, we are handling both the standard
// wells and the multi-segment wells under the MultisegmentWells
// we need them to remove the dependency on StandardWells
Vector well_perforation_densities_;
Vector well_perforation_pressure_diffs_;
};
} // namespace Opm
#include "MultisegmentWells_impl.hpp"
#endif // OPM_MULTISEGMENTWELLS_HEADER_INCLUDED

1178
opm/autodiff/MultisegmentWells_impl.hpp
View File

File diff suppressed because it is too large Load Diff

402
opm/autodiff/WellMultiSegment.cpp
View File

@ -1,402 +0,0 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/WellMultiSegment.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/CompletionSet.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
namespace Opm
{
WellMultiSegment::WellMultiSegment(const Well* well, size_t time_step, const Wells* wells) {
m_well_name_ = well->name();
if (well->isMultiSegment(time_step)) {
initMultiSegmentWell(well, time_step, wells);
} else {
initNonMultiSegmentWell(well, time_step, wells);
}
updateWellOps();
}
void WellMultiSegment::initMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells) {
const auto& completion_set = well->getCompletions(time_step);
m_is_multi_segment_ = true;
const auto& segment_set = well->getSegmentSet(time_step);
m_number_of_segments_ = segment_set.numberSegment();
m_number_of_perforations_ = completion_set.size();
m_comp_pressure_drop_ = segment_set.compPressureDrop();
m_multiphase_model_ = segment_set.multiPhaseModel();
m_outlet_segment_.resize(m_number_of_segments_);
m_inlet_segments_.resize(m_number_of_segments_);
m_segment_length_.resize(m_number_of_segments_);
m_segment_depth_.resize(m_number_of_segments_);
m_segment_internal_diameter_.resize(m_number_of_segments_);
m_segment_roughness_.resize(m_number_of_segments_);
m_segment_cross_area_.resize(m_number_of_segments_, 0.);
m_segment_volume_.resize(m_number_of_segments_);
m_segment_perforations_.resize(m_number_of_segments_);
// we change the ID to location now for easier use later.
for (int i = 0; i < m_number_of_segments_; ++i) {
// The segment number for top segment is 0, the segment number of its outlet segment will be -1
m_outlet_segment_[i] = segment_set.segmentNumberToIndex(segment_set[i].outletSegment());
m_segment_length_[i] = segment_set[i].totalLength();
m_segment_depth_[i] = segment_set[i].depth();
m_segment_internal_diameter_[i] = segment_set[i].internalDiameter();
m_segment_roughness_[i] = segment_set[i].roughness();
m_segment_cross_area_[i] = segment_set[i].crossArea();
m_segment_volume_[i] = segment_set[i].volume();
}
// update the completion related information
// find the location of the well in wells
int index_well;
for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
if (m_well_name_ == std::string(wells->name[index_well])) {
break;
}
}
std::vector<int> temp_well_cell;
std::vector<double> temp_well_index;
if (index_well == wells->number_of_wells) {
throw std::runtime_error(" did not find the well " + m_well_name_ + "\n");
} else {
m_well_type_ = wells->type[index_well];
m_well_controls_ = wells->ctrls[index_well];
m_number_of_phases_ = wells->number_of_phases;
m_comp_frac_.resize(m_number_of_phases_);
std::copy(wells->comp_frac + index_well * m_number_of_phases_,
wells->comp_frac + (index_well + 1) * m_number_of_phases_, m_comp_frac_.begin());
int index_begin = wells->well_connpos[index_well];
int index_end = wells->well_connpos[index_well + 1];
for(int i = index_begin; i < index_end; ++i) {
// copy the WI and well_cell_ informatin to m_well_index_ and m_well_cell_
// maybe also the depth of the perforations.
temp_well_cell.push_back(wells->well_cells[i]);
temp_well_index.push_back(wells->WI[i]);
}
}
std::vector<double> temp_perf_depth;
temp_perf_depth.resize(m_number_of_perforations_);
for (int i = 0; i < (int)completion_set.size(); ++i) {
int i_segment = completion_set.get(i).getSegmentNumber();
// using the location of the segment in the array as the segment number/id.
// TODO: it can be helpful for output or postprocessing if we can keep the original number.
i_segment = segment_set.segmentNumberToIndex(i_segment);
m_segment_perforations_[i_segment].push_back(i);
temp_perf_depth[i] = completion_set.get(i).getCenterDepth();
}
// reordering the perforation related informations
// so that the perforations belong to the same segment will be continuous
m_well_cell_.resize(m_number_of_perforations_);
m_well_index_.resize(m_number_of_perforations_);
m_perf_depth_.resize(m_number_of_perforations_);
m_segment_cell_.resize(m_number_of_segments_, -1);
int perf_count = 0;
for (int is = 0; is < m_number_of_segments_; ++is) {
// TODO: the grid cell related to a segment should be calculated based on the location
// of the segment node.
// As the current temporary solution, the grid cell related to a segment determined by the
// first perforation cell related to the segment.
// when no perforation is related to the segment, use it outlet segment's cell.
const int nperf = m_segment_perforations_[is].size();
if (nperf > 0) {
const int first_perf_number = m_segment_perforations_[is][0];
m_segment_cell_[is] = temp_well_cell[first_perf_number];
for (int iperf = 0; iperf < nperf; ++iperf) {
const int perf_number = m_segment_perforations_[is][iperf];
m_well_cell_[perf_count] = temp_well_cell[perf_number];
m_well_index_[perf_count] = temp_well_index[perf_number];
m_perf_depth_[perf_count] = temp_perf_depth[perf_number];
m_segment_perforations_[is][iperf] = perf_count;
++perf_count;
}
} else {
// using the cell of its outlet segment
const int i_outlet_segment = m_outlet_segment_[is];
if (i_outlet_segment < 0) {
assert(is == 0); // it must be the top segment
OPM_THROW(std::logic_error, "Top segment is not related to any perforation, its related cell must be calculated based the location of its segment node, which is not implemented yet \n");
} else {
if (m_well_cell_[i_outlet_segment] < 0) {
OPM_THROW(std::logic_error, "The segment cell of its outlet segment is not determined yet, the current implementation does not support this \n");
} else {
m_segment_cell_[is] = m_segment_cell_[i_outlet_segment];
}
}
}
}
assert(perf_count == m_number_of_perforations_);
// update m_inlet_segments_
// top segment does not have a outlet segment
for (int is = 1; is < m_number_of_segments_; ++is) {
const int index_outlet = m_outlet_segment_[is];
m_inlet_segments_[index_outlet].push_back(is);
}
}
void WellMultiSegment::initNonMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells) {
const auto& completion_set = well->getCompletions(time_step);
m_is_multi_segment_ = false;
m_number_of_segments_ = 1;
m_comp_pressure_drop_ = WellSegment::H__;
m_multiphase_model_ = WellSegment::HO;
m_outlet_segment_.resize(m_number_of_segments_, -1);
m_segment_length_.resize(m_number_of_segments_, 0.);
m_segment_depth_.resize(m_number_of_segments_, 0.);
m_segment_internal_diameter_.resize(m_number_of_segments_, 0.);
m_segment_roughness_.resize(m_number_of_segments_, 0.);
m_segment_cross_area_.resize(m_number_of_segments_, 0.);
m_segment_volume_.resize(m_number_of_segments_, 0.);
m_segment_perforations_.resize(m_number_of_segments_);
// update the completion related information
int index_well;
for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
if (m_well_name_ == std::string(wells->name[index_well])) {
break;
}
}
if (index_well == wells->number_of_wells) {
throw std::runtime_error(" did not find the well " + m_well_name_ + "\n");
} else {
m_well_type_ = wells->type[index_well];
m_well_controls_ = wells->ctrls[index_well];
m_number_of_phases_ = wells->number_of_phases;
// set the segment depth to be the bhp reference depth
m_segment_depth_[0] = wells->depth_ref[index_well];
m_comp_frac_.resize(m_number_of_phases_);
std::copy(wells->comp_frac + index_well * m_number_of_phases_,
wells->comp_frac + (index_well + 1) * m_number_of_phases_, m_comp_frac_.begin());
int index_begin = wells->well_connpos[index_well];
int index_end = wells->well_connpos[index_well + 1];
m_number_of_perforations_ = index_end - index_begin;
for(int i = index_begin; i < index_end; ++i) {
m_well_cell_.push_back(wells->well_cells[i]);
m_well_index_.push_back(wells->WI[i]);
}
m_segment_cell_.resize(1, -1);
m_segment_cell_[0] = m_well_cell_[0];
}
// TODO: not sure if we need the perf_depth_.
m_perf_depth_.resize(m_number_of_perforations_, 0.);
m_segment_perforations_[0].resize(m_number_of_perforations_);
for (int i = 0; i < m_number_of_perforations_; ++i) {
m_segment_perforations_[0][i] = i;
m_perf_depth_[i] = completion_set.get(i).getCenterDepth();
}
m_inlet_segments_.resize(m_number_of_segments_);
}
void WellMultiSegment::updateWellOps() {
m_wops_.s2p = Matrix(m_number_of_perforations_, m_number_of_segments_);
m_wops_.p2s = Matrix(m_number_of_segments_, m_number_of_perforations_);
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> s2p;
std::vector<Tri> p2s;
s2p.reserve(m_number_of_perforations_);
p2s.reserve(m_number_of_perforations_);
for(int s = 0; s < (int)m_number_of_segments_; ++s) {
int temp_nperf = m_segment_perforations_[s].size();
// some segment may not have any perforation
assert(temp_nperf >= 0);
for (int perf = 0; perf < temp_nperf; ++perf) {
const int index_perf = m_segment_perforations_[s][perf];
s2p.push_back(Tri(index_perf, s, 1.0));
p2s.push_back(Tri(s, index_perf, 1.0));
}
}
m_wops_.s2p.setFromTriplets(s2p.begin(), s2p.end());
m_wops_.p2s.setFromTriplets(p2s.begin(), p2s.end());
m_wops_.s2s_gather = Matrix(m_number_of_segments_, m_number_of_segments_);
std::vector<Tri> s2s_gather;
s2s_gather.reserve(m_number_of_segments_ * m_number_of_segments_);
std::vector<Tri> s2s_inlets;
s2s_inlets.reserve(m_number_of_segments_);
std::vector<Tri> s2s_outlet;
s2s_outlet.reserve(m_number_of_segments_);
for (int s = 0; s < (int)m_number_of_segments_; ++s) {
s2s_gather.push_back(Tri(s, s, 1.0));
int s_outlet = m_outlet_segment_[s];
if (s_outlet >=0) {
s2s_inlets.push_back(Tri(s_outlet, s, 1.0));
s2s_outlet.push_back(Tri(s, s_outlet, 1.0));
}
int temp_s = s;
while (m_outlet_segment_[temp_s] >=0) {
s2s_gather.push_back(Tri(m_outlet_segment_[temp_s], s, 1.0));
temp_s = m_outlet_segment_[temp_s];
}
}
m_wops_.s2s_gather.setFromTriplets(s2s_gather.begin(), s2s_gather.end());
m_wops_.p2s_gather = Matrix(m_number_of_segments_, m_number_of_perforations_);
m_wops_.p2s_gather = m_wops_.s2s_gather * m_wops_.p2s;
m_wops_.s2s_inlets = Matrix(m_number_of_segments_, m_number_of_segments_);
m_wops_.s2s_inlets.setFromTriplets(s2s_inlets.begin(), s2s_inlets.end());
m_wops_.s2s_outlet = Matrix(m_number_of_segments_, m_number_of_segments_);
m_wops_.s2s_outlet.setFromTriplets(s2s_outlet.begin(), s2s_outlet.end());
m_wops_.p2s_average = Matrix(m_number_of_segments_, m_number_of_perforations_);
std::vector<Tri> p2s_average;
p2s_average.reserve(m_number_of_segments_);
for (int s = 0; s < (int)m_number_of_segments_; ++s) {
const int nperf = m_segment_perforations_[s].size();
if (nperf > 0) {
p2s_average.push_back(Tri(s, s, 1.0/nperf));
}
}
// constructing the diagonal matrix to do the averaging for p2s
Matrix temp_averaging_p2s = Matrix(m_number_of_segments_, m_number_of_segments_);
temp_averaging_p2s.setFromTriplets(p2s_average.begin(), p2s_average.end());
m_wops_.p2s_average = temp_averaging_p2s * m_wops_.p2s;
}
const std::string& WellMultiSegment::name() const {
return m_well_name_;
}
bool WellMultiSegment::isMultiSegmented() const {
return m_is_multi_segment_;
}
WellType WellMultiSegment::wellType() const {
return m_well_type_;
}
const WellControls* WellMultiSegment::wellControls() const {
return m_well_controls_;
}
int WellMultiSegment::numberOfPerforations() const {
return m_number_of_perforations_;
}
int WellMultiSegment::numberOfSegments() const {
return m_number_of_segments_;
}
std::string WellMultiSegment::compPressureDrop() const {
return WellSegment::CompPressureDropEnumToString(m_comp_pressure_drop_);
}
const std::vector<double>& WellMultiSegment::compFrac() const {
return m_comp_frac_;
}
int WellMultiSegment::numberOfPhases() const {
return m_number_of_phases_;
}
const std::vector<double>& WellMultiSegment::wellIndex() const {
return m_well_index_;
}
const std::vector<double>& WellMultiSegment::perfDepth() const {
return m_perf_depth_;
}
const std::vector<int>& WellMultiSegment::wellCells() const {
return m_well_cell_;
}
const std::vector<int>& WellMultiSegment::segmentCells() const {
return m_segment_cell_;
}
const std::vector<int>& WellMultiSegment::outletSegment() const {
return m_outlet_segment_;
}
const std::vector<std::vector<int>>& WellMultiSegment::inletSegments() const {
return m_inlet_segments_;
}
const std::vector<double>& WellMultiSegment::segmentLength() const {
return m_segment_length_;
}
const std::vector<double>& WellMultiSegment::segmentDepth() const {
return m_segment_depth_;
}
const std::vector<double>& WellMultiSegment::segmentDiameter() const {
return m_segment_internal_diameter_;
}
const std::vector<double>& WellMultiSegment::segmentCrossArea() const {
return m_segment_cross_area_;
}
const std::vector<double>& WellMultiSegment::segmentRoughness() const {
return m_segment_roughness_;
}
const std::vector<double>& WellMultiSegment::segmentVolume() const {
return m_segment_volume_;
}
const std::vector<std::vector<int>>& WellMultiSegment::segmentPerforations() const {
return m_segment_perforations_;
}
const WellMultiSegment::WellOps& WellMultiSegment::wellOps() const {
return m_wops_;
}
}

227
opm/autodiff/WellMultiSegment.hpp
View File

@ -1,227 +0,0 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLMULTISEGMENT_HEADER_INCLUDED
#define OPM_WELLMULTISEGMENT_HEADER_INCLUDED
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/MSW/SegmentSet.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <vector>
#include <cassert>
#include <string>
#include <utility>
#include <map>
#include <algorithm>
#include <array>
namespace Opm
{
class WellMultiSegment
{
public:
typedef Eigen::SparseMatrix<double> Matrix;
/// Constructor of WellMultiSegment
/// \param[in] well information from EclipseState
/// \param[in] current time step
/// \param[in[ pointer to Wells structure, to be removed eventually
WellMultiSegment(const Well* well, size_t time_step, const Wells* wells);
/// Well name.
const std::string& name() const;
/// Flag indicating if the well is a multi-segment well.
bool isMultiSegmented() const;
/// Number of the perforations.
int numberOfPerforations() const;
/// Number of the segments.
int numberOfSegments() const;
/// Components of the pressure drop invloved.
/// HFA Hydrostatic + friction + acceleration
/// HF- Hydrostatic + friction
/// H-- Hydrostatic only.
std::string compPressureDrop() const;
/// Well control.
const WellControls* wellControls() const;
/// Component fractions for each well.
const std::vector<double>& compFrac() const;
/// Number of phases.
int numberOfPhases() const;
/// Well type.
WellType wellType() const;
/// Well productivity index.
const std::vector<double>& wellIndex() const;
/// Depth of the perforations.
const std::vector<double>& perfDepth() const;
/// Indices of the grid blocks that perforations are completed in.
const std::vector<int>& wellCells() const;
/// Indices of the gird blocks that segments locate at.
const std::vector<int>& segmentCells() const;
/// Outlet segments, a segment (except top segment) can only have one outlet segment.
/// For top segment, its outlet segments is -1 always, which means no outlet segment for top segment.
const std::vector<int>& outletSegment() const;
/// Inlet segments. a segment can have more than one inlet segments.
const std::vector<std::vector<int>>& inletSegments() const;
/// The length of the segment nodes down the wellbore from the reference point.
const std::vector<double>& segmentLength() const;
/// The depth of the segment nodes.
const std::vector<double>& segmentDepth() const;
/// Tubing internal diameter.
const std::vector<double>& segmentDiameter() const;
/// Cross-sectional area.
const std::vector<double>& segmentCrossArea() const;
/// Effective absolute roughness of the tube.
const std::vector<double>& segmentRoughness() const;
/// Volume of segment.
const std::vector<double>& segmentVolume() const;
/// Perforations related to each segment.
const std::vector<std::vector<int>>& segmentPerforations() const;
/// Struct for the well operator matrices.
/// All the operator matrics only apply to the one specifi well.
struct WellOps {
Matrix s2p; // segment -> perf (scatter)
Matrix p2s; // perf -> segment (gather)
Matrix p2s_average; // perf -> segment (avarage)
Matrix s2s_gather; // segment -> segment (in an accumlative way)
// means the outlet segments will gather all the contribution
// from all the inlet segments in a recurisive way
Matrix p2s_gather; // perforation -> segment (in an accumative way)
Matrix s2s_inlets; // segment -> its inlet segments
Matrix s2s_outlet; // segment -> its outlet segment
};
/// Well operator matrics
const WellOps& wellOps() const;
private:
// for the moment, we use the information from wells.
// TODO: remove the dependency on wells from opm-core.
void initMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells);
void initNonMultiSegmentWell(const Well* well, size_t time_step, const Wells* wells);
void updateWellOps();
private:
// name of the well
std::string m_well_name_;
// flag to indicate if this well is a
// multi-segmented well
bool m_is_multi_segment_;
// well type
// INJECTOR or PRODUCER
enum WellType m_well_type_;
// number of phases
int m_number_of_phases_;
// component fractions for each well
std::vector<double> m_comp_frac_;
// controls for this well
// using pointer for temporary
// changing it when figuring out how to using it
struct WellControls *m_well_controls_;
// components of the pressure drop to be included
WellSegment::CompPressureDropEnum m_comp_pressure_drop_;
// multi-phase flow model
WellSegment::MultiPhaseModelEnum m_multiphase_model_;
// number of perforation for this well
int m_number_of_perforations_;
// number of segments for this well
int m_number_of_segments_;
// well index for each completion
std::vector<double> m_well_index_;
// depth for each completion // form the keyword COMPSEGS?
std::vector<double> m_perf_depth_;
// well cell for each completion
std::vector<int> m_well_cell_;
// cell for each segment
std::vector<int> m_segment_cell_;
// how to organize the segment structure here?
// indicate the outlet segment for each segment
// maybe here we can use the location in the vector
// at the moment, we still use the ID number
// then a mapping from the ID number to the actual location will be required
// The ID is already changed to the location now.
std::vector<int> m_outlet_segment_;
// for convinience, we store the inlet segments for each segment
std::vector<std::vector<int>> m_inlet_segments_;
// this one is not necessary any more, since the segment number is its location.
// std::map<int, int> m_number_to_location_;
// has not decided to use the absolute length from the well head
// or the length of this single segment
// using the absolute length first
std::vector<double> m_segment_length_;
// the depth of the segmnet node
std::vector<double> m_segment_depth_;
// the internal diameter of the segment
std::vector<double> m_segment_internal_diameter_;
// the roughness of the segment
std::vector<double> m_segment_roughness_;
// the cross area of the segment
std::vector<double> m_segment_cross_area_;
// the volume of the segment
std::vector<double> m_segment_volume_;
// the completions that is related to each segment
// the completions's ids are their location in the vector m_well_index_
// m_well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
std::vector<std::vector<int>> m_segment_perforations_;
WellOps m_wops_;
};
typedef std::shared_ptr<WellMultiSegment> WellMultiSegmentPtr;
typedef std::shared_ptr<const WellMultiSegment> WellMultiSegmentConstPtr;
} // namespace Opm
#endif // OPM_WELLMULTISEGMENT_HEADER_INCLUDE

364
opm/autodiff/WellStateMultiSegment.hpp
View File

@ -1,364 +0,0 @@
/*
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
#define OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
// #include <opm/autodiff/WellMultiSegment.hpp>
#include <opm/autodiff/MultisegmentWells.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <vector>
#include <cassert>
#include <string>
#include <utility>
#include <map>
#include <algorithm>
#include <array>
namespace Opm
{
/// The state of a set of multi-sgemnet wells
// Since we are avoiding to use the old Wells structure,
// it might be a good idea not to relate this State to the WellState much.
class WellStateMultiSegment
: public WellStateFullyImplicitBlackoil
{
public:
typedef WellStateFullyImplicitBlackoil Base;
typedef struct {
int well_number;
int start_segment;
int number_of_segments;
int start_perforation;
int number_of_perforations;
std::vector<int> start_perforation_segment; // the starting position of perforation inside the segment
std::vector<int> number_of_perforations_segment; // the numbers for perforations for the segments
} SegmentedMapentryType;
typedef std::map<std::string, SegmentedMapentryType> SegmentedWellMapType;
/// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls
template <class ReservoirState, class PrevWellState>
void init(const MultisegmentWells& ms_wells, const ReservoirState& state, const PrevWellState& prevState, const Wells* legacy_wells_ptr)
{
// Used by output facilities.
this->wells_.reset( clone_wells( legacy_wells_ptr ) );
const std::vector<WellMultiSegmentConstPtr>& wells = ms_wells.msWells();
const int nw = wells.size();
nseg_ = 0;
nperf_ = 0;
if (nw == 0) {
perfPhaseRates().clear();
perfPress().clear();
segphaserates_.clear();
segpress_.clear();
return;
}
const int np = wells[0]->numberOfPhases(); // number of the phases
for (int iw = 0; iw < nw; ++iw) {
nperf_ += wells[iw]->numberOfPerforations();
nseg_ += wells[iw]->numberOfSegments();
}
bhp().resize(nw);
thp().resize(nw);
top_segment_loc_.resize(nw);
temperature().resize(nw, 273.15 + 20); // standard temperature for now
wellRates().resize(nw * np, 0.0);
currentControls().resize(nw);
for(int iw = 0; iw < nw; ++iw) {
currentControls()[iw] = well_controls_get_current(wells[iw]->wellControls());
}
for (int iw = 0; iw < nw; ++iw) {
assert((wells[iw]->wellType() == INJECTOR) || (wells[iw]->wellType() == PRODUCER));
}
int start_segment = 0;
int start_perforation = 0;
perfPhaseRates().resize(nperf_ * np, 0.0);
perfPress().resize(nperf_, -1.0e100);
segphaserates_.resize(nseg_ * np, 0.0);
segpress_.resize(nseg_, -1.0e100);
for (int w = 0; w < nw; ++w) {
assert((wells[w]->wellType() == INJECTOR) || (wells[w]->wellType() == PRODUCER));
const WellControls* ctrl = wells[w]->wellControls();
std::string well_name(wells[w]->name());
// Initialize the wellMap_ here
SegmentedMapentryType& wellMapEntry = segmentedWellMap_[well_name];
wellMapEntry.well_number = w;
wellMapEntry.start_segment = start_segment;
wellMapEntry.number_of_segments = wells[w]->numberOfSegments();
wellMapEntry.start_perforation = start_perforation;
wellMapEntry.number_of_perforations = wells[w]->numberOfPerforations();
top_segment_loc_[w] = start_segment;
int start_perforation_segment = 0;
wellMapEntry.start_perforation_segment.resize(wellMapEntry.number_of_segments);
wellMapEntry.number_of_perforations_segment.resize(wellMapEntry.number_of_segments);
for (int i = 0; i < wellMapEntry.number_of_segments; ++i) {
wellMapEntry.start_perforation_segment[i] = start_perforation_segment;
wellMapEntry.number_of_perforations_segment[i] = wells[w]->segmentPerforations()[i].size();
start_perforation_segment += wellMapEntry.number_of_perforations_segment[i];
}
assert(start_perforation_segment == wellMapEntry.number_of_perforations);
if (well_controls_well_is_stopped(ctrl)) {
// 1. WellRates: 0
// 2. Bhp: assign bhp equal to bhp control, if applicable, otherwise
// assign equal to first perforation cell pressure.
if (well_controls_get_current_type(ctrl) == BHP) {
bhp()[w] = well_controls_get_current_target(ctrl);
} else {
const int first_cell = wells[0]->wellCells()[0];
bhp()[w] = state.pressure()[first_cell];
}
// 3. Thp: assign thp equal to thp control, if applicable,
// otherwise assign equal to bhp value.
if (well_controls_get_current_type(ctrl) == THP) {
thp()[w] = well_controls_get_current_target( ctrl );
} else {
thp()[w] = bhp()[w];
}
// 4. Perforation pressures and phase rates
// 5. Segment pressures and phase rates
} else {
// Open Wells
// 1. Rates: initialize well rates to match controls if type is SURFACE_RATE. Otherwise, we
// cannot set the correct value here, so we aasign a small rate with the correct sign so that any
// logic depending on that sign will work as expected.
if (well_controls_get_current_type(ctrl) == SURFACE_RATE) {
const double rate_target = well_controls_get_current_target(ctrl);
const double * distr = well_controls_get_current_distr( ctrl );
for (int p = 0; p < np; ++p) {
wellRates()[np * w + p] = rate_target * distr[p];
}
} else {
const double small_rate = 1e-14;
const double sign = (wells[w]->wellType() == INJECTOR) ? 1.0 : -1.0;
for (int p = 0; p < np; ++p) {
wellRates()[np * w + p] = small_rate * sign;
}
}
// 2. Bhp:
if (well_controls_get_current_type(ctrl) == BHP) {
bhp()[w] = well_controls_get_current_target(ctrl);
} else {
const int first_cell = wells[w]->wellCells()[0];
const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99;
bhp()[w] = safety_factor* state.pressure()[first_cell];
}
// 3. Thp:
if (well_controls_get_current_type(ctrl) == THP) {
thp()[w] = well_controls_get_current_target(ctrl);
} else {
thp()[w] = bhp()[w];
}
// 4. Perf rates and pressures
int number_of_perforations = wellMapEntry.number_of_perforations;
for (int i = 0; i < number_of_perforations; ++i) {
for (int p = 0; p < np; ++p) {
perfPhaseRates()[np * (i + start_perforation) + p] = wellRates()[np * w + p] / double(number_of_perforations);
}
if (wells[w]->isMultiSegmented()) {
const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99;
perfPress()[i + start_perforation] = safety_factor * state.pressure()[wells[w]->wellCells()[i]];
} else {
perfPress()[i + start_perforation] = state.pressure()[wells[w]->wellCells()[i]];
}
}
// 5. Segment rates and pressures
int number_of_segments = wellMapEntry.number_of_segments;
// the seg_pressure is the same with the first perf_pressure. For the top segment, it is the same with bhp,
// when under bhp control.
// the seg_rates will related to the sum of the perforation rates, and also trying to keep consistent with the
// well rates. Most importantly, the segment rates of the top segment is the same with the well rates.
segpress_[start_segment] = bhp()[w];
for (int i = 1; i < number_of_segments; ++i) {
int first_perforation_segment = start_perforation + wellMapEntry.start_perforation_segment[i];
segpress_[i + start_segment] = perfPress()[first_perforation_segment];
// the segmnent pressure of the top segment should be the bhp
}
for (int p = 0; p < np; ++p) {
Eigen::VectorXd v_perf_rates(number_of_perforations);
for (int i = 0; i < number_of_perforations; ++i) {
v_perf_rates[i] = perfPhaseRates()[np * (i + start_perforation) + p];
}
Eigen::VectorXd v_segment_rates = wells[w]->wellOps().p2s_gather * v_perf_rates;
for (int i = 0; i < number_of_segments; ++i) {
segphaserates_[np * (i + start_segment) + p] = v_segment_rates[i];
}
}
}
start_segment += wellMapEntry.number_of_segments;
start_perforation += wellMapEntry.number_of_perforations;
}
// Initialize current_controls_.
// The controls set in the Wells object are treated as defaults,
// and also used for initial values.
currentControls().resize(nw);
for (int w = 0; w < nw; ++w) {
currentControls()[w] = well_controls_get_current(wells[w]->wellControls());
}
// initialize wells that have been there before
// order can change so the mapping is based on the well names
if ( !(prevState.segmentedWellMap().empty()) )
{
typedef typename SegmentedWellMapType::const_iterator const_iterator;
const_iterator end_old = prevState.segmentedWellMap().end();
for (int w = 0; w < nw; ++w) {
std::string well_name(wells[w]->name());
const_iterator it_old = prevState.segmentedWellMap().find(well_name);
const_iterator it_this = segmentedWellMap().find(well_name);
assert(it_this != segmentedWellMap().end()); // the current well must be present in the current well map
if (it_old != end_old) {
const int oldIndex = (*it_old).second.well_number;
const int newIndex = w;
// bhp
bhp()[newIndex] = prevState.bhp()[oldIndex];
// well rates
for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
{
wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
}
const int num_seg_old_well = (*it_old).second.number_of_segments;
const int num_perf_old_well = (*it_old).second.number_of_perforations;
const int num_seg_this_well = (*it_this).second.number_of_segments;
const int num_perf_this_well = (*it_this).second.number_of_perforations;
// determing if the structure of the wells has been changed by comparing the number of segments and perforations
// may not be very safe.
// The strategy HAS to be changed later with experiments and analysis
if ((num_perf_old_well == num_perf_this_well) && (num_seg_old_well == num_seg_this_well)) {
const int old_start_perforation = (*it_old).second.start_perforation;
const int old_start_segment = (*it_old).second.start_segment;
const int this_start_perforation = (*it_this).second.start_perforation;
const int this_start_segment = (*it_this).second.start_segment;
// this is not good when the the well rates changed dramatically
for (int i = 0; i < num_seg_this_well * np; ++i) {
segphaserates_[this_start_segment * np + i] = prevState.segPhaseRates()[old_start_segment * np + i];
}
for (int i = 0; i < num_perf_this_well * np; ++i) {
perfPhaseRates()[this_start_perforation * np + i] = prevState.perfPhaseRates()[old_start_perforation * np + i];
}
// perf_pressures_
for (int i = 0; i < num_perf_this_well; ++i) {
// p
perfPress()[this_start_perforation + i] = prevState.perfPress()[old_start_perforation + i];
}
// segpress_
for (int i = 0; i < num_seg_this_well; ++i) {
// p
segpress_[this_start_segment + i] = prevState.segPress()[old_start_segment + i];
}
// current controls
const int old_control_index = prevState.currentControls()[ oldIndex ];
if (old_control_index < well_controls_get_num(wells[w]->wellControls())) {
// If the set of controls have changed, this may not be identical
// to the last control, but it must be a valid control.
currentControls()[ newIndex ] = old_control_index;
}
}
}
}
}
}
std::vector<double>& segPress() { return segpress_; }
const std::vector<double>& segPress() const { return segpress_; }
std::vector<double>& segPhaseRates() { return segphaserates_; }
const std::vector<double>& segPhaseRates() const { return segphaserates_; }
const std::vector<int>& topSegmentLoc() const { return top_segment_loc_; };
const SegmentedWellMapType& segmentedWellMap() const { return segmentedWellMap_; }
SegmentedWellMapType& segmentedWellMap() { return segmentedWellMap_; }
int numSegments() const { return nseg_; }
int numPerforations() const { return nperf_; }
private:
// pressure for the segment nodes
std::vector<double> segpress_;
// phase rates for the segments
std::vector<double> segphaserates_;
// the location of the top segments within the whole segment list
// it is better in the Wells class if we have a class instead of
// using a vector for all the wells
std::vector<int> top_segment_loc_;
SegmentedWellMapType segmentedWellMap_;
int nseg_;
int nperf_;
};
} // namespace Opm
#endif // OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDE

174
tests/test_multisegmentwells.cpp
View File

@ -1,174 +0,0 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE MultisegmentWellsTest
#define BOOST_TEST_NO_MAIN
#include <vector>
#include <unordered_set>
#include <memory>
#include <array>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <boost/filesystem.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/DynamicListEconLimited.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/MultisegmentWells.hpp>
struct SetupMSW {
using Grid = UnstructuredGrid;
using GridInit = Opm::GridInit<Grid>;
using FluidProps = Opm::BlackoilPropsAdFromDeck;
using MaterialLawManager = FluidProps::MaterialLawManager;
SetupMSW()
{
Opm::ParseContext parse_context;
Opm::Parser parser;
auto deck = parser.parseFile("msw.data", parse_context);
Opm::EclipseState ecl_state(deck , parse_context);
Opm::Schedule schedule(deck, ecl_state.getInputGrid(), ecl_state.get3DProperties(), Opm::Phases(true, true, true), parse_context );
// Create grid.
const std::vector<double>& porv =
ecl_state.get3DProperties().getDoubleGridProperty("PORV").getData();
std::unique_ptr<GridInit> grid_init(new GridInit(ecl_state, porv));
const Grid& grid = grid_init->grid();
const size_t current_timestep = 0;
// dummy_dynamic_list_econ_lmited
const Opm::DynamicListEconLimited dummy_dynamic_list;
// Create wells.
Opm::WellsManager wells_manager(ecl_state,
schedule,
current_timestep,
Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
Opm::UgGridHelpers::dimensions(grid),
Opm::UgGridHelpers::cell2Faces(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
dummy_dynamic_list,
false,
// We need to pass the optionaly arguments
// as we get the following error otherwise
// with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
// converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructor
std::unordered_set<std::string>());
const Wells* wells = wells_manager.c_wells();
const auto& wells_ecl = schedule.getWells(current_timestep);
ms_wells.reset(new Opm::MultisegmentWells(wells, &(wells_manager.wellCollection()), wells_ecl, current_timestep));
};
std::shared_ptr<const Opm::MultisegmentWells> ms_wells;
};
// number of wells for this case
const int nw = 2;
// number of segments for this case
const int nseg = 7;
// number of perforations for this case
const int nperf = 8;
BOOST_AUTO_TEST_CASE(testOperators)
{
SetupMSW msw_setup;
const std::shared_ptr<const Opm::MultisegmentWells>& ms_wells = msw_setup.ms_wells;
const Opm::MultisegmentWells::MultisegmentWellOps& wops_ms = ms_wells->wellOps();
BOOST_CHECK_EQUAL(true, wops_ms.has_multisegment_wells);
BOOST_CHECK_EQUAL(nperf, wops_ms.well_cells.size());
BOOST_CHECK_EQUAL(nperf, wops_ms.conn_trans_factors.size());
BOOST_CHECK_EQUAL(nseg, wops_ms.s2s_outlet.rows());
BOOST_CHECK_EQUAL(nseg, wops_ms.s2s_inlets.cols());
BOOST_CHECK_EQUAL(nw, wops_ms.s2w.rows());
BOOST_CHECK_EQUAL(nseg, wops_ms.s2w.cols());
BOOST_CHECK_EQUAL(nseg, wops_ms.topseg2w.cols());
BOOST_CHECK_EQUAL(nperf, wops_ms.s2p.rows());
BOOST_CHECK_EQUAL(nseg, wops_ms.p2s.rows());
BOOST_CHECK_EQUAL(nperf, wops_ms.w2p.rows());
BOOST_CHECK_EQUAL(nw, wops_ms.w2p.cols());
BOOST_CHECK_EQUAL(nperf, wops_ms.p2w.cols());
}
BOOST_AUTO_TEST_CASE(testStructure)
{
SetupMSW msw_setup;
const std::shared_ptr<const Opm::MultisegmentWells>& ms_wells = msw_setup.ms_wells;
BOOST_CHECK_EQUAL(3, ms_wells->numPhases());
BOOST_CHECK_EQUAL(nseg, ms_wells->numSegment());
BOOST_CHECK_EQUAL(nperf, ms_wells->numPerf());
BOOST_CHECK_EQUAL(nw, ms_wells->topWellSegments().size());
BOOST_CHECK_EQUAL(nw, ms_wells->msWells().size());
BOOST_CHECK_EQUAL(0, ms_wells->topWellSegments()[0]);
BOOST_CHECK_EQUAL(1, ms_wells->topWellSegments()[1]);
}
bool
init_unit_test_func()
{
return true;
}
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
boost::unit_test::unit_test_main(&init_unit_test_func,
argc, argv);
}