00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef OPM_WELLHELPERS_HEADER_INCLUDED
00023 #define OPM_WELLHELPERS_HEADER_INCLUDED
00024
00025 #include <opm/core/wells.h>
00026
00027
00028 #include <vector>
00029
00030 namespace Opm {
00031
00032
00033
00034
00035 namespace wellhelpers
00036 {
00037
00038
00039 inline
00040 double rateToCompare(const std::vector<double>& well_phase_flow_rate,
00041 const int well,
00042 const int num_phases,
00043 const double* distr)
00044 {
00045 double rate = 0.0;
00046 for (int phase = 0; phase < num_phases; ++phase) {
00047
00048
00049 rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase];
00050 }
00051 return rate;
00052 }
00053
00054 inline
00055 bool constraintBroken(const std::vector<double>& bhp,
00056 const std::vector<double>& thp,
00057 const std::vector<double>& well_phase_flow_rate,
00058 const int well,
00059 const int num_phases,
00060 const WellType& well_type,
00061 const WellControls* wc,
00062 const int ctrl_index)
00063 {
00064 const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
00065 const double target = well_controls_iget_target(wc, ctrl_index);
00066 const double* distr = well_controls_iget_distr(wc, ctrl_index);
00067
00068 bool broken = false;
00069
00070 switch (well_type) {
00071 case INJECTOR:
00072 {
00073 switch (ctrl_type) {
00074 case BHP:
00075 broken = bhp[well] > target;
00076 break;
00077
00078 case THP:
00079 broken = thp[well] > target;
00080 break;
00081
00082 case RESERVOIR_RATE:
00083 case SURFACE_RATE:
00084 broken = rateToCompare(well_phase_flow_rate,
00085 well, num_phases, distr) > target;
00086 break;
00087 }
00088 }
00089 break;
00090
00091 case PRODUCER:
00092 {
00093 switch (ctrl_type) {
00094 case BHP:
00095 broken = bhp[well] < target;
00096 break;
00097
00098 case THP:
00099 broken = thp[well] < target;
00100 break;
00101
00102 case RESERVOIR_RATE:
00103 case SURFACE_RATE:
00104
00105
00106
00107 broken = rateToCompare(well_phase_flow_rate,
00108 well, num_phases, distr) < target;
00109 break;
00110 }
00111 }
00112 break;
00113
00114 default:
00115 OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
00116 }
00117
00118 return broken;
00119 }
00120
00121
00122
00123
00132 inline
00133 double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
00134 const double& rho, const double gravity) {
00135 if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
00136 {
00137
00138
00139
00140
00141 return 0;
00142 }
00143 const double well_ref_depth = wells.depth_ref[w];
00144 const double dh = vfp_ref_depth - well_ref_depth;
00145 const double dp = rho*gravity*dh;
00146
00147 return dp;
00148 }
00149
00150 inline
00151 double computeHydrostaticCorrection(const double well_ref_depth, const double vfp_ref_depth,
00152 const double rho, const double gravity) {
00153 const double dh = vfp_ref_depth - well_ref_depth;
00154 const double dp = rho * gravity * dh;
00155
00156 return dp;
00157 }
00158
00159 template <class Vector>
00160 inline
00161 Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,
00162 const Vector& well_perforation_densities, const double gravity) {
00163 const int nw = wells.number_of_wells;
00164 Vector retval = Vector::Zero(nw);
00165
00166 #pragma omp parallel for schedule(static)
00167 for (int i=0; i<nw; ++i) {
00168 const int perf = wells.well_connpos[i];
00169 retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
00170 }
00171
00172 return retval;
00173 }
00174
00175 inline
00176 std::vector<double> computeHydrostaticCorrection(const Wells& wells, const std::vector<double>& vfp_ref_depth,
00177 const std::vector<double>& well_perforation_densities, const double gravity) {
00178 const int nw = wells.number_of_wells;
00179 std::vector<double> retval(nw,0.0);
00180
00181 #pragma omp parallel for schedule(static)
00182 for (int i=0; i<nw; ++i) {
00183 const int perf = wells.well_connpos[i];
00184 retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
00185 }
00186
00187 return retval;
00188 }
00189
00190 }
00191
00192 }
00193
00194 #endif