All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
WellHelpers.hpp
1 /*
2  Copyright 2016 SINTEF ICT, Applied Mathematics.
3  Copyright 2016 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #ifndef OPM_WELLHELPERS_HEADER_INCLUDED
23 #define OPM_WELLHELPERS_HEADER_INCLUDED
24 
25 #include <opm/core/wells.h>
26 // #include <opm/autodiff/AutoDiffHelpers.hpp>
27 
28 #include <vector>
29 
30 namespace Opm {
31 
32 
33 
34 
35  namespace wellhelpers
36  {
37 
38 
39  inline
40  double rateToCompare(const std::vector<double>& well_phase_flow_rate,
41  const int well,
42  const int num_phases,
43  const double* distr)
44  {
45  double rate = 0.0;
46  for (int phase = 0; phase < num_phases; ++phase) {
47  // Important: well_phase_flow_rate is ordered with all phase rates for first
48  // well first, then all phase rates for second well etc.
49  rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase];
50  }
51  return rate;
52  }
53 
54  inline
55  bool constraintBroken(const std::vector<double>& bhp,
56  const std::vector<double>& thp,
57  const std::vector<double>& well_phase_flow_rate,
58  const int well,
59  const int num_phases,
60  const WellType& well_type,
61  const WellControls* wc,
62  const int ctrl_index)
63  {
64  const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
65  const double target = well_controls_iget_target(wc, ctrl_index);
66  const double* distr = well_controls_iget_distr(wc, ctrl_index);
67 
68  bool broken = false;
69 
70  switch (well_type) {
71  case INJECTOR:
72  {
73  switch (ctrl_type) {
74  case BHP:
75  broken = bhp[well] > target;
76  break;
77 
78  case THP:
79  broken = thp[well] > target;
80  break;
81 
82  case RESERVOIR_RATE: // Intentional fall-through
83  case SURFACE_RATE:
84  broken = rateToCompare(well_phase_flow_rate,
85  well, num_phases, distr) > target;
86  break;
87  }
88  }
89  break;
90 
91  case PRODUCER:
92  {
93  switch (ctrl_type) {
94  case BHP:
95  broken = bhp[well] < target;
96  break;
97 
98  case THP:
99  broken = thp[well] < target;
100  break;
101 
102  case RESERVOIR_RATE: // Intentional fall-through
103  case SURFACE_RATE:
104  // Note that the rates compared below are negative,
105  // so breaking the constraints means: too high flow rate
106  // (as for injection).
107  broken = rateToCompare(well_phase_flow_rate,
108  well, num_phases, distr) < target;
109  break;
110  }
111  }
112  break;
113 
114  default:
115  OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
116  }
117 
118  return broken;
119  }
120 
121 
122  // --------- Types ---------
123 
132  inline
133  double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
134  const double& rho, const double gravity) {
135  if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
136  {
137  // This is a well with no perforations.
138  // If this is the last well we would subscript over the
139  // bounds below.
140  // we assume well_perforation_densities to be 0
141  return 0;
142  }
143  const double well_ref_depth = wells.depth_ref[w];
144  const double dh = vfp_ref_depth - well_ref_depth;
145  const double dp = rho*gravity*dh;
146 
147  return dp;
148  }
149 
150  inline
151  double computeHydrostaticCorrection(const double well_ref_depth, const double vfp_ref_depth,
152  const double rho, const double gravity) {
153  const double dh = vfp_ref_depth - well_ref_depth;
154  const double dp = rho * gravity * dh;
155 
156  return dp;
157  }
158 
159  template <class Vector>
160  inline
161  Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,
162  const Vector& well_perforation_densities, const double gravity) {
163  const int nw = wells.number_of_wells;
164  Vector retval = Vector::Zero(nw);
165 
166 #pragma omp parallel for schedule(static)
167  for (int i=0; i<nw; ++i) {
168  const int perf = wells.well_connpos[i];
169  retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
170  }
171 
172  return retval;
173  }
174 
175  inline
176  std::vector<double> computeHydrostaticCorrection(const Wells& wells, const std::vector<double>& vfp_ref_depth,
177  const std::vector<double>& well_perforation_densities, const double gravity) {
178  const int nw = wells.number_of_wells;
179  std::vector<double> retval(nw,0.0);
180 
181 #pragma omp parallel for schedule(static)
182  for (int i=0; i<nw; ++i) {
183  const int perf = wells.well_connpos[i];
184  retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
185  }
186 
187  return retval;
188  }
189 
190  } // namespace wellhelpers
191 
192 }
193 
194 #endif