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

a wIP version.

Browse Source
This commit is contained in:
Kai Bao 2015-08-28 15:57:27 +02:00
parent 62bdd301d3
commit 07709c52c6
6 changed files with 1647 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
opm/autodiff/BlackoilModel.hpp
View File

@ -27,6 +27,7 @@
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateMultiSegment.hpp>
namespace Opm {

119
opm/autodiff/BlackoilModelBase.hpp
View File

@ -34,6 +34,10 @@
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
// temporary usuage
#include <opm/autodiff/WellStateMultiSegment.hpp>
#include <opm/autodiff/WellMultiSegment.hpp>
#include <array>
struct Wells;
@ -74,7 +78,35 @@ namespace Opm {
};
struct MultiSegmentBlackoilSolutionState
{
typedef AutoDiffBlock<double> ADB;
explicit MultiSegmentBlackoilSolutionState(const int np)
: pressure ( ADB::null())
, temperature ( ADB::null())
, saturation (np, ADB::null())
, rs ( ADB::null())
, rv ( ADB::null())
, qs ( ADB::null())
, pseg ( ADB::null())
, canonical_phase_pressures(3, ADB::null())
{
}
ADB pressure;
ADB temperature;
std::vector<ADB> saturation;
ADB rs;
ADB rv;
// the flow rates for each segments
// the first one for each well is the flow rate
ADB qs;
// the pressure for the segments
// the first one for each well is the bhp
ADB pseg;
// Below are quantities stored in the state for optimization purposes.
std::vector<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
};
/// Class used for reporting the outcome of a nonlinearIteration() call.
struct IterationReport
@ -199,6 +231,10 @@ namespace Opm {
WellState& well_state,
const bool initial_assembly);
void assemble(const ReservoirState& reservoir_state,
WellStateMultiSegment& well_state,
const bool initial_assembly);
/// \brief Compute the residual norms of the mass balance for each phase,
/// the well flux, and the well equation.
/// \return a vector that contains for each phase the norm of the mass balance
@ -282,6 +318,9 @@ namespace Opm {
const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
const Wells* wells_;
// FOR TEMPORARY
// SHOUlD BE A REFERENCE
const std::vector<WellMultiSegment> wells_multi_segment_;
VFPProperties vfp_properties_;
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
@ -304,9 +343,48 @@ namespace Opm {
V isRs_;
V isRv_;
V isSg_;
// For the non-segmented well, it should be the density with AVG or SEG way.
// while usually SEG way
V well_perforation_densities_; //Density of each well perforation
// ADB version, when using AVG way, the calculation of the density and hydrostatic head
// is implicit
ADB well_perforation_densities_adb_;
// 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.
V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
// ADB version. Eventually, only ADB version will be kept.
ADB well_perforation_pressure_diffs_adb_;
// 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.
V well_perforation_pressure_cell_diffs_;
ADB well_perforation_pressure_cell_diffs_adb_;
// Pressure correction due to the depth differennce between segment depth and perforation depth.
// TODO: It should be able to be merge as a part of the perforation_pressure_diffs_.
ADB well_perforations_segment_pressure_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
V well_perforation_cell_densities_;
ADB well_perforation_cell_densities_adb_;
V well_perforatoin_cell_pressure_diffs_;
LinearisedBlackoilResidual residual_;
/// \brief Whether we print something to std::cout
@ -341,8 +419,10 @@ namespace Opm {
bool wellsActive() const { return wells_active_; }
// return true if wells are available on this process
bool localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; }
// return wells object
const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
const std::vector<WellMultiSegment>& wellsMultiSegment() const { return wells_multi_segment_; }
void
makeConstantState(SolutionState& state) const;
@ -351,9 +431,16 @@ namespace Opm {
variableState(const ReservoirState& x,
const WellState& xw) const;
SolutionState
variableState(const ReservoirState& x,
const WellStateMultiSegment& xw) const;
std::vector<V>
variableStateInitials(const ReservoirState& x,
const WellState& xw) const;
std::vector<V>
variableStateInitials(const ReservoirState& x,
const WellStateMultiSegment& xw) const;
void
variableReservoirStateInitials(const ReservoirState& x,
std::vector<V>& vars0) const;
@ -361,6 +448,13 @@ namespace Opm {
variableWellStateInitials(const WellState& xw,
std::vector<V>& vars0) const;
void variableWellStateInitials(const WellStateMultiSegment& xw,
std::vector<V>& vars0) const;
void
variableWellState(const WellStateMultiSegment& xw,
std::vector<V>& vars0) const;
std::vector<int>
variableStateIndices() const;
@ -384,6 +478,9 @@ namespace Opm {
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw);
void computeWellConnectionPressures(const SolutionState& state,
const WellStateMultiSegment& xw);
void
assembleMassBalanceEq(const SolutionState& state);
@ -401,14 +498,29 @@ namespace Opm {
std::vector<ADB>& cq_s);
void
computeWellFlux(const MultiSegmentBlackoilSolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s);
void
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw);
void
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const MultiSegmentBlackoilSolutionState& state,
WellStateMultiSegment& xw);
void
addWellFluxEq(const std::vector<ADB>& cq_s,
const SolutionState& state);
void
addWellFluxEq(const std::vector<ADB>& cq_s,
const MultiSegmentBlackoilSolutionState& state);
void
addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
const SolutionState& state,
@ -419,7 +531,14 @@ namespace Opm {
const WellState& xw,
const V& aliveWells);
void
addWellControlEq(const MultiSegmentBlackoilSolutionState& state,
const WellStateMultiSegment& xw,
const V& aliveWells);
void updateWellControls(WellState& xw) const;
void updateWellControls(WellStateMultiSegment& xw) const;
void updateWellState(const V& dwells,
WellState& well_state);

1074
opm/autodiff/BlackoilModelBase_impl.hpp
View File

File diff suppressed because it is too large Load Diff

82
opm/autodiff/WellMultiSegment.cpp
View File

@ -46,6 +46,7 @@ namespace Opm
// m_number_of_perforations_ from wells
// m_well_index_ from wells
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_);
@ -135,6 +136,13 @@ namespace Opm
}
assert(perf_count == m_number_of_perforations_);
// update m_inlet_segments_
for (size_t is = 0; is < m_number_of_segments_; ++is) {
const int index_outlet = m_outlet_segment_[is];
m_inlet_segments_[index_outlet].push_back(is);
}
// std::cin.ignore();
// how to build the relation between segments and completions
@ -207,6 +215,8 @@ namespace Opm
m_segment_perforations_[0][i] = i;
}
m_inlet_segments_.resize(1);
// std::cin.ignore();
// how to build the relation between segments and completions
@ -234,6 +244,68 @@ namespace Opm
// while for usual wells, it is typically SEG.
}
// update the wellOps (m_wops)
m_wops_.s2p = M(m_number_of_perforations_, m_number_of_segments_);
m_wops_.p2s = M(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_.size();
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 = M(m_number_of_segments_, m_number_of_segments_);
std::vector<Tri> s2s_gather;
s2s_gather.reserve(m_number_of_segments_ * m_number_of_segments_);
// a brutal way first
// will generate matrix with entries bigger than 1.0
// Then we need to normalize all the values.
for (int s = 0; s < (int)m_number_of_segments_; ++s) {
s2s_gather.push_back(Tri(s, s, 1.0));
int temp_s = s;
while (m_outlet_segment_[temp_s] >=0) {
s2s_gather.push_back(Tri(m_outlet_segment_[temp_s], temp_s, 1.0));
temp_s = m_outlet_segment_[temp_s];
}
}
M temp_s2s_gather(m_number_of_segments_, m_number_of_segments_);
M inverse_s2s_gather(m_number_of_segments_, m_number_of_segments_);
temp_s2s_gather.setFromTriplets(s2s_gather.begin(), s2s_gather.end());
inverse_s2s_gather = temp_s2s_gather.cwiseInverse();
m_wops_.s2s_gather = temp_s2s_gather.cwiseProduct(inverse_s2s_gather);
// p2w should be simple
m_wops_.p2s_gather = M(m_number_of_segments_, m_number_of_perforations_);
m_wops_.p2s_gather = m_wops_.s2s_gather * m_wops_.s2p;
// s2s_gather
}
const std::string& WellMultiSegment::name() const {
return m_well_name_;
}
const bool WellMultiSegment::isMultiSegmented() const {
return m_is_multi_segment_;
}
const enum WellType WellMultiSegment::wellType() const {
@ -248,7 +320,7 @@ namespace Opm
return m_number_of_perforations_;
}
const size_t WellMultiSegment::numberOfSegment() const {
const size_t WellMultiSegment::numberOfSegments() const {
return m_number_of_segments_;
}
@ -268,7 +340,7 @@ namespace Opm
return m_perf_depth_;
}
const std::vector<int>& WellMultiSegment::wellCell() const {
const std::vector<int>& WellMultiSegment::wellCells() const {
return m_well_cell_;
}
@ -296,7 +368,11 @@ namespace Opm
return m_segment_volume_;
}
const std::vector<std::vector<int>>& WellMultiSegment::segmentPerforatioins() const {
const std::vector<std::vector<int>>& WellMultiSegment::segmentPerforations() const {
return m_segment_perforations_;
}
const WellMultiSegment::WellOps& WellMultiSegment::wellOps() const {
return m_wops_;
}
}

41
opm/autodiff/WellMultiSegment.hpp
View File

@ -21,6 +21,11 @@
#define OPM_WELLMULTISEGMENT_HEADER_INCLUDED
#include <opm/core/utility/platform_dependent/disable_warnings.h>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
@ -43,10 +48,15 @@ namespace Opm
{
public:
typedef Eigen::Array<double, Eigen::Dynamic, 1> V;
typedef Eigen::SparseMatrix<double> M;
WellMultiSegment(WellConstPtr well, size_t time_step, const Wells* wells);
const std::string& name() const;
const bool isMultiSegmented() const;
const size_t numberOfPerforations() const;
const size_t numberOfSegment() const;
const size_t numberOfSegments() const;
const struct WellControls* wellControls() const;
const std::vector<double>& compFrac() const;
@ -56,14 +66,31 @@ namespace Opm
const enum WellType wellType() const;
const std::vector<double>& wellIndex() const;
const std::vector<double>& perfDepth() const;
const std::vector<int>& wellCell() const;
const std::vector<int>& wellCells() const;
const std::vector<int>& outletSegment() const;
const std::vector<std::vector<int>>& inletSegments() const;
const std::vector<double>& segmentLength() const;
const std::vector<double>& segmentDepth() const;
const std::vector<double>& segmentCrossArea() const;
const std::vector<double>& segmentRoughness() const;
const std::vector<double>& segmentVolume() const;
const std::vector<std::vector<int>>& segmentPerforatioins() const;
const std::vector<std::vector<int>>& segmentPerforations() const;
struct WellOps {
M s2p; // segment -> perf (scatter)
M p2s; // perf -> segment (gather)
// M w2p; // well -> perf (scatter)
// M p2w; // perf - > well (gather)
// but since only one well, so it is just an arrary
// not needed now.
// M w2s; // well -> segment (scatter)
M 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
M p2s_gather; // perforation -> segment (in an accumative way)
};
const WellOps& wellOps() const;
private:
// for the moment, we use the information from wells.
@ -106,8 +133,12 @@ namespace Opm
// 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_;
std::map<int, int> m_number_to_location_;
// 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
@ -128,6 +159,8 @@ namespace Opm
// 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> WellMutliSegmentPtr;

352
opm/autodiff/WellStateMultiSegment.hpp
View File

@ -25,6 +25,7 @@
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/WellMultiSegment.hpp>
#include <vector>
#include <cassert>
@ -38,24 +39,35 @@ namespace Opm
{
/// The state of a set of multi-sgemnet wells
/// since we are avoiding to use the old wells structure
/// it makes it might be a good idea not to relate this State to the WellState
class WellStateMultiSegment
: public WellState
// : public WellState
{
public:
typedef WellMultiSegment::V V;
// typedef std::array< int, 3 > mapentry_t;
// typedef std::map< std::string, mapentry_t > WellMapType;
// this map needs to change a little bit?
/* struct mapentry {
typedef struct {
int well_number;
int start_segment;
int number_of_segments;
std::vector<int> number_of_performations;
} */
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
} MapentryType;
typedef std::map<std::string, MapentryType> WellMapType;
// MAYNOT NEED THIS
/// 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
/// the PrevState here must be the same with State
template <class State, class PrevState>
void init(const std::vector<WellMultiSegment>& wells, const State& state, const PrevState& prevState)
{
@ -70,13 +82,21 @@ namespace Opm
int nseg = 0; // the nubmer of the segments
for (int iw = 0; iw < nw; ++iw) {
nperf += wells[i].numberOfPerforations();
nseg += wells[i].numberOfSegment();
nperf += wells[iw].numberOfPerforations();
nseg += wells[iw].numberOfSegments();
}
bhp_.resize(nw);
thp_.resize(nw);
wellrates_.resize(nw * np, 0.0);
// deciding to add the following variables temporarily
// TODO: making it better later
np_ = np;
nseg_ = nseg;
nperf_ = nperf_;
nwells_ = nw;
well_rates_.resize(nw * np, 0.0);
current_controls_.resize(nw);
for(int iw = 0; iw < nw; ++iw) {
@ -85,28 +105,336 @@ namespace Opm
for (int iw = 0; iw < nw; ++iw) {
assert((wells[i].wellType() == INJECTOR) || (wells[i].wellType() == PRODUCER));
const WellControls* ctrl = wells[iw]->wellControls();
const WellControls* ctrl = wells[iw].wellControls();
}
int start_segment = 0;
int start_perforation = 0;
// Map is used to map the value from the previous state to the current state as the initial values
// TODO: handle this later.
// Trying to figure out the work flow first.
perf_phaserates_.clear();
perf_phaserates_.resize(nperf * np, 0.0);
perf_pressures_.clear();
perf_pressures_.resize(nperf, -1.0e100);
seg_phaserates_.clear();
seg_phaserates_.resize(nseg * np, 0.0);
seg_pressures_.clear();
seg_pressures_.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
MapentryType& wellMapEntry = wellMap_[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();
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) {
well_rates_[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) {
well_rates_[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[0].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) {
perf_phaserates_[np * (i + start_perforation) + p] = well_rates_[np * w + p] / double(number_of_perforations);
}
perf_pressures_[i + start_perforation] = state.pressure()[wells[w].wellCells()[i + start_perforation]];
}
// 5. Segment rates and pressures
int number_of_segments = wellMapEntry.number_of_segments;
// the seg_pressure is the same with the first pref_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
for (int i = 0; i < number_of_segments; ++i) {
/* for (int p = 0; p < np; ++p) {
seg_phaserates_[np * (i + start_segment) + p] = 0.;
} */
int first_perforation_segment = start_perforation + wellMapEntry.start_perforation_segment[i];
seg_pressures_[i + start_segment] = perf_pressures_[first_perforation_segment];
// the segmnent pressure of the top segment should be the bhp
}
for (int p = 0; p < np; ++p) {
// std::vector<double> v_phase_rates(number_of_perforations);
// V v_perf_rates = V::Zero(number_of_perforations);
Eigen::VectorXd v_perf_rates(number_of_perforations);
for (int i = 0; i < number_of_perforations; ++i) {
v_perf_rates[i] = perf_phaserates_[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) {
seg_phaserates_[np * (i + start_segment) + p] = v_segment_rates[i];
}
}
// initialize the segmnet rates.
// it works in the analog way with the usual wells.
// How to initialize the perforation rates and the segment rates.?
// Perforation pressures can be set to the pressure of the corresponding grid cells?
// deviding the well rates by the number of the perforations
// then calculating the segment rate based on the rates for perforations and
// make sure the flow rate for the top segment is consistent with the well flow rates
// for pressure it is not that trival
}
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.
current_controls_.resize(nw);
for (int w = 0; w < nw; ++w) {
current_controls_[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.wellMap().empty()) )
{
typedef typename WellMapType::const_iterator const_iterator;
const_iterator end_old = prevState.WellMap().end;
const_iterator end_this = wellMap().end;
for (int w = 0; w < nw; ++w) {
std::string well_name(wells[w].name());
const_iterator it_old = prevState.wellMap().find(well_name);
const_iterator it_this = wellMap().find(well_name);
assert(it_this != end_this); // 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_perf_old_well = (*it_old).second.number_of_perforations;
const int num_seg_old_well = (*it_old).second.number_of_segments;
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_seg_old_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) {
seg_phaserates_[this_start_segment * np + i] = prevState.segPhaseRates()[old_start_segment * np + i];
}
for (int i = 0; i < num_perf_this_well * np; ++i) {
perf_phaserates_[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
perf_pressures_[this_start_perforation + i] = prevState.perfPressures()[old_start_perforation + i];
}
// seg_pressures_
for (int i = 0; i < num_seg_this_well; ++i) {
// p
seg_pressures_[this_start_segment + i] = prevState.segPressures()[old_start_segment + i];
}
// current controls
}
// else {
// deviding the well rates by the number of the perforations
// then calculating the segment rate based on the rates for perforations and
// make sure the flow rate for the top segment is consistent with the well flow rates
// for pressure it is not that trival
// }
/* 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
} MapentryType; */
// peforation rates
// segment rates
// It really depends on if the structures on the segments and perforations are changed.
// TODO: if it will be reasonable to assume that if the numbers of segments and perforations are same, then
// the structures of the wells are not changed.
// Obviously it is not true.
// for the perforation rates, it is Okay to calculate by deviding the well rates by the perforation numbers.
// Then the segment rates are calculated based on the perforation rates and the well rates.
// The segment rates of top segments should be the same with the well rates.
}
}
}
// TODO: maybe we should store the values of np, nw, nseg, nperf for the states for later use.
// assert(start_perforation == total_perforation);
// assert(start_segment == total_segment);
/* if (well_controls_well_is_stopped(ctrl)) {
// shut well: all the rates are zero.
} else {
// Initialize the phase rates for each perforation by deviding the well rates by the number of perofrations
// Then using the perf rates to initialize the rates for the segments
} */
}
std::vector<double>& segPhaseRates() { return seg_phaserates_; }
const std::vector<double>& segPhaseRates() const { return seg_phaserates_; }
std::vector<double>& segPressures() { return seg_pressures_; }
const std::vector<double>& segPressures() const { return seg_pressures_; }
std::vector<double>& perfPressures() { return perf_pressures_; }
const std::vector<double>& perfPressures() const { return perf_pressures_; }
std::vector<double>& perfPhaseRates() { return perf_phaserates_; }
const std::vector<double>& perfPhaseRates() const { return perf_phaserates_; }
std::vector<double>& bhp() { return bhp_; }
const std::vector<double>& bhp() const { return bhp_; }
std::vector<double>& thp() { return thp_; }
const std::vector<double>& thp() const { return thp_; }
std::vector<double>& wellRates() { return well_rates_; }
const std::vector<double>& wellRates() const { return well_rates_; }
std::vector<int>& currentControls() { return current_controls_; }
const std::vector<int>& currentControls() const { return current_controls_; }
// wellrate should be the out segment rates for the top segments
const WellMapType& wellMap() const { return wellMap_; }
WellMapType& wellMap() { return wellMap_; }
const int numberOfPhases() const { return np_; }
const int numberOfSegments() const { return nseg_; }
const int numberOfPerforations() const { return nperf_; }
const int numberOfWells() const { return nwells_; }
private:
std::vector<double> bhp_;
std::vector<double> thp_;
std::vector<double> well_rates_;
// pressure for the segment nodes
std::vector<double> seg_pressure_;
std::vector<double> seg_pressures_;
// phase rates for the segments
std::vector<double> seg_phaserates_;
// phase rates for the completions
std::vector<double> perf_phaserates_;
// pressure for the perforatins
std::vector<double> perf_pressures_;
// TODO: MIGHT NOT USE THE FOLLOWING VARIABLES AT THE
// fractions for each segments (W, O, G)
std::vector<double> seg_phasefrac_;
// total flow rates for each segments, G_T
std::vector<double> seg_totalrate_;
std::vector<int> current_controls_;
WellMapType wellMap_;
int nseg_;
int np_;
int nperf_;
int nwells_;
};
} // namespace Opm