00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
00021 #define OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
00022
00023
00024 #include <opm/core/wells.h>
00025 #include <opm/core/well_controls.h>
00026 #include <opm/common/ErrorMacros.hpp>
00027 #include <opm/autodiff/AutoDiffBlock.hpp>
00028
00029 #include <opm/autodiff/MultisegmentWells.hpp>
00030 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
00031 #include <vector>
00032 #include <cassert>
00033 #include <string>
00034 #include <utility>
00035 #include <map>
00036 #include <algorithm>
00037 #include <array>
00038
00039 namespace Opm
00040 {
00041
00043
00044
00045 class WellStateMultiSegment
00046 : public WellStateFullyImplicitBlackoil
00047 {
00048 public:
00049
00050 typedef WellStateFullyImplicitBlackoil Base;
00051
00052 typedef struct {
00053 int well_number;
00054 int start_segment;
00055 int number_of_segments;
00056 int start_perforation;
00057 int number_of_perforations;
00058 std::vector<int> start_perforation_segment;
00059 std::vector<int> number_of_perforations_segment;
00060 } SegmentedMapentryType;
00061
00062 typedef std::map<std::string, SegmentedMapentryType> SegmentedWellMapType;
00063
00067 template <class ReservoirState, class PrevWellState>
00068 void init(const MultisegmentWells& ms_wells, const ReservoirState& state, const PrevWellState& prevState, const Wells* legacy_wells_ptr)
00069 {
00070
00071 this->wells_.reset( clone_wells( legacy_wells_ptr ) );
00072
00073 const std::vector<WellMultiSegmentConstPtr>& wells = ms_wells.msWells();
00074 const int nw = wells.size();
00075 nseg_ = 0;
00076 nperf_ = 0;
00077 if (nw == 0) {
00078 perfPhaseRates().clear();
00079 perfPress().clear();
00080 segphaserates_.clear();
00081 segpress_.clear();
00082 return;
00083 }
00084
00085 const int np = wells[0]->numberOfPhases();
00086
00087 for (int iw = 0; iw < nw; ++iw) {
00088 nperf_ += wells[iw]->numberOfPerforations();
00089 nseg_ += wells[iw]->numberOfSegments();
00090 }
00091
00092 bhp().resize(nw);
00093 thp().resize(nw);
00094 top_segment_loc_.resize(nw);
00095 temperature().resize(nw, 273.15 + 20);
00096
00097 wellRates().resize(nw * np, 0.0);
00098
00099 currentControls().resize(nw);
00100 for(int iw = 0; iw < nw; ++iw) {
00101 currentControls()[iw] = well_controls_get_current(wells[iw]->wellControls());
00102 }
00103
00104 for (int iw = 0; iw < nw; ++iw) {
00105 assert((wells[iw]->wellType() == INJECTOR) || (wells[iw]->wellType() == PRODUCER));
00106 }
00107
00108 int start_segment = 0;
00109 int start_perforation = 0;
00110
00111 perfPhaseRates().resize(nperf_ * np, 0.0);
00112 perfPress().resize(nperf_, -1.0e100);
00113 segphaserates_.resize(nseg_ * np, 0.0);
00114 segpress_.resize(nseg_, -1.0e100);
00115
00116
00117 for (int w = 0; w < nw; ++w) {
00118 assert((wells[w]->wellType() == INJECTOR) || (wells[w]->wellType() == PRODUCER));
00119 const WellControls* ctrl = wells[w]->wellControls();
00120
00121 std::string well_name(wells[w]->name());
00122
00123 SegmentedMapentryType& wellMapEntry = segmentedWellMap_[well_name];
00124 wellMapEntry.well_number = w;
00125 wellMapEntry.start_segment = start_segment;
00126 wellMapEntry.number_of_segments = wells[w]->numberOfSegments();
00127 wellMapEntry.start_perforation = start_perforation;
00128 wellMapEntry.number_of_perforations = wells[w]->numberOfPerforations();
00129
00130 top_segment_loc_[w] = start_segment;
00131
00132 int start_perforation_segment = 0;
00133 wellMapEntry.start_perforation_segment.resize(wellMapEntry.number_of_segments);
00134 wellMapEntry.number_of_perforations_segment.resize(wellMapEntry.number_of_segments);
00135
00136 for (int i = 0; i < wellMapEntry.number_of_segments; ++i) {
00137 wellMapEntry.start_perforation_segment[i] = start_perforation_segment;
00138 wellMapEntry.number_of_perforations_segment[i] = wells[w]->segmentPerforations()[i].size();
00139 start_perforation_segment += wellMapEntry.number_of_perforations_segment[i];
00140 }
00141 assert(start_perforation_segment == wellMapEntry.number_of_perforations);
00142
00143
00144 if (well_controls_well_is_stopped(ctrl)) {
00145
00146
00147
00148 if (well_controls_get_current_type(ctrl) == BHP) {
00149 bhp()[w] = well_controls_get_current_target(ctrl);
00150 } else {
00151 const int first_cell = wells[0]->wellCells()[0];
00152 bhp()[w] = state.pressure()[first_cell];
00153 }
00154
00155
00156 if (well_controls_get_current_type(ctrl) == THP) {
00157 thp()[w] = well_controls_get_current_target( ctrl );
00158 } else {
00159 thp()[w] = bhp()[w];
00160 }
00161
00162
00163
00164 } else {
00165
00166
00167
00168
00169 if (well_controls_get_current_type(ctrl) == SURFACE_RATE) {
00170 const double rate_target = well_controls_get_current_target(ctrl);
00171 const double * distr = well_controls_get_current_distr( ctrl );
00172 for (int p = 0; p < np; ++p) {
00173 wellRates()[np * w + p] = rate_target * distr[p];
00174 }
00175 } else {
00176 const double small_rate = 1e-14;
00177 const double sign = (wells[w]->wellType() == INJECTOR) ? 1.0 : -1.0;
00178 for (int p = 0; p < np; ++p) {
00179 wellRates()[np * w + p] = small_rate * sign;
00180 }
00181 }
00182
00183
00184 if (well_controls_get_current_type(ctrl) == BHP) {
00185 bhp()[w] = well_controls_get_current_target(ctrl);
00186 } else {
00187 const int first_cell = wells[w]->wellCells()[0];
00188 const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99;
00189 bhp()[w] = safety_factor* state.pressure()[first_cell];
00190 }
00191
00192 if (well_controls_get_current_type(ctrl) == THP) {
00193 thp()[w] = well_controls_get_current_target(ctrl);
00194 } else {
00195 thp()[w] = bhp()[w];
00196 }
00197
00198
00199 int number_of_perforations = wellMapEntry.number_of_perforations;
00200 for (int i = 0; i < number_of_perforations; ++i) {
00201 for (int p = 0; p < np; ++p) {
00202 perfPhaseRates()[np * (i + start_perforation) + p] = wellRates()[np * w + p] / double(number_of_perforations);
00203 }
00204 if (wells[w]->isMultiSegmented()) {
00205 const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99;
00206 perfPress()[i + start_perforation] = safety_factor * state.pressure()[wells[w]->wellCells()[i]];
00207 } else {
00208 perfPress()[i + start_perforation] = state.pressure()[wells[w]->wellCells()[i]];
00209 }
00210 }
00211
00212
00213 int number_of_segments = wellMapEntry.number_of_segments;
00214
00215
00216
00217
00218 segpress_[start_segment] = bhp()[w];
00219 for (int i = 1; i < number_of_segments; ++i) {
00220 int first_perforation_segment = start_perforation + wellMapEntry.start_perforation_segment[i];
00221 segpress_[i + start_segment] = perfPress()[first_perforation_segment];
00222
00223 }
00224
00225 for (int p = 0; p < np; ++p) {
00226 Eigen::VectorXd v_perf_rates(number_of_perforations);
00227 for (int i = 0; i < number_of_perforations; ++i) {
00228 v_perf_rates[i] = perfPhaseRates()[np * (i + start_perforation) + p];
00229 }
00230
00231 Eigen::VectorXd v_segment_rates = wells[w]->wellOps().p2s_gather * v_perf_rates;
00232
00233 for (int i = 0; i < number_of_segments; ++i) {
00234 segphaserates_[np * (i + start_segment) + p] = v_segment_rates[i];
00235 }
00236 }
00237 }
00238
00239 start_segment += wellMapEntry.number_of_segments;
00240 start_perforation += wellMapEntry.number_of_perforations;
00241
00242 }
00243
00244
00245
00246
00247 currentControls().resize(nw);
00248 for (int w = 0; w < nw; ++w) {
00249 currentControls()[w] = well_controls_get_current(wells[w]->wellControls());
00250 }
00251
00252
00253
00254 if ( !(prevState.segmentedWellMap().empty()) )
00255 {
00256 typedef typename SegmentedWellMapType::const_iterator const_iterator;
00257 const_iterator end_old = prevState.segmentedWellMap().end();
00258
00259 for (int w = 0; w < nw; ++w) {
00260 std::string well_name(wells[w]->name());
00261 const_iterator it_old = prevState.segmentedWellMap().find(well_name);
00262 const_iterator it_this = segmentedWellMap().find(well_name);
00263
00264 assert(it_this != segmentedWellMap().end());
00265
00266 if (it_old != end_old) {
00267 const int oldIndex = (*it_old).second.well_number;
00268 const int newIndex = w;
00269
00270
00271 bhp()[newIndex] = prevState.bhp()[oldIndex];
00272
00273
00274 for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
00275 {
00276 wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
00277 }
00278
00279 const int num_seg_old_well = (*it_old).second.number_of_segments;
00280 const int num_perf_old_well = (*it_old).second.number_of_perforations;
00281
00282 const int num_seg_this_well = (*it_this).second.number_of_segments;
00283 const int num_perf_this_well = (*it_this).second.number_of_perforations;
00284
00285
00286
00287
00288 if ((num_perf_old_well == num_perf_this_well) && (num_seg_old_well == num_seg_this_well)) {
00289 const int old_start_perforation = (*it_old).second.start_perforation;
00290 const int old_start_segment = (*it_old).second.start_segment;
00291
00292 const int this_start_perforation = (*it_this).second.start_perforation;
00293 const int this_start_segment = (*it_this).second.start_segment;
00294
00295
00296 for (int i = 0; i < num_seg_this_well * np; ++i) {
00297 segphaserates_[this_start_segment * np + i] = prevState.segPhaseRates()[old_start_segment * np + i];
00298 }
00299
00300 for (int i = 0; i < num_perf_this_well * np; ++i) {
00301 perfPhaseRates()[this_start_perforation * np + i] = prevState.perfPhaseRates()[old_start_perforation * np + i];
00302 }
00303
00304
00305 for (int i = 0; i < num_perf_this_well; ++i) {
00306
00307 perfPress()[this_start_perforation + i] = prevState.perfPress()[old_start_perforation + i];
00308 }
00309
00310
00311 for (int i = 0; i < num_seg_this_well; ++i) {
00312
00313 segpress_[this_start_segment + i] = prevState.segPress()[old_start_segment + i];
00314 }
00315
00316
00317 const int old_control_index = prevState.currentControls()[ oldIndex ];
00318 if (old_control_index < well_controls_get_num(wells[w]->wellControls())) {
00319
00320
00321 currentControls()[ newIndex ] = old_control_index;
00322 }
00323 }
00324 }
00325 }
00326 }
00327 }
00328
00329
00330 std::vector<double>& segPress() { return segpress_; }
00331 const std::vector<double>& segPress() const { return segpress_; }
00332
00333 std::vector<double>& segPhaseRates() { return segphaserates_; }
00334 const std::vector<double>& segPhaseRates() const { return segphaserates_; }
00335
00336 const std::vector<int>& topSegmentLoc() const { return top_segment_loc_; };
00337
00338 const SegmentedWellMapType& segmentedWellMap() const { return segmentedWellMap_; }
00339 SegmentedWellMapType& segmentedWellMap() { return segmentedWellMap_; }
00340
00341 int numSegments() const { return nseg_; }
00342 int numPerforations() const { return nperf_; }
00343
00344 private:
00345
00346 std::vector<double> segpress_;
00347
00348 std::vector<double> segphaserates_;
00349
00350
00351
00352
00353 std::vector<int> top_segment_loc_;
00354
00355 SegmentedWellMapType segmentedWellMap_;
00356
00357 int nseg_;
00358 int nperf_;
00359 };
00360
00361 }
00362
00363
00364 #endif // OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDE