All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
SimulatorFullyImplicitBlackoilPolymer_impl.hpp
1 /*
2  Copyright 2013 SINTEF ICT, Applied Mathematics.
3  Copyright 2014 IRIS AS
4  Copyright 2014 STATOIL ASA.
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 namespace Opm
23 {
24  template <class GridT>
25  SimulatorFullyImplicitBlackoilPolymer<GridT>::
26  SimulatorFullyImplicitBlackoilPolymer(const ParameterGroup& param,
27  const GridT& grid,
28  DerivedGeology& geo,
29  BlackoilPropsAdFromDeck& props,
30  const PolymerPropsAd& polymer_props,
31  const RockCompressibility* rock_comp_props,
32  NewtonIterationBlackoilInterface& linsolver,
33  const double* gravity,
34  const bool has_disgas,
35  const bool has_vapoil,
36  const bool has_polymer,
37  const bool has_plyshlog,
38  const bool has_shrate,
39  std::shared_ptr<EclipseState> eclipse_state,
40  BlackoilOutputWriter& output_writer,
41  std::shared_ptr< Deck > deck,
42  const std::vector<double>& threshold_pressures_by_face)
43  : BaseType(param,
44  grid,
45  geo,
46  props,
47  rock_comp_props,
48  linsolver,
49  gravity,
50  has_disgas,
51  has_vapoil,
52  eclipse_state,
53  output_writer,
54  threshold_pressures_by_face,
55  // names of deactivated wells in parallel run
56  std::unordered_set<std::string>())
57  , polymer_props_(polymer_props)
58  , has_polymer_(has_polymer)
59  , has_plyshlog_(has_plyshlog)
60  , has_shrate_(has_shrate)
61  , deck_(deck)
62  {
63  }
64 
65  template <class GridT>
66  auto SimulatorFullyImplicitBlackoilPolymer<GridT>::
67  createSolver(const WellModel& well_model)
68  -> std::unique_ptr<Solver>
69  {
70  typedef typename Traits::Model Model;
71 
72 
73  auto model = std::unique_ptr<Model>(new Model(BaseType::model_param_,
74  BaseType::grid_,
75  BaseType::props_,
76  BaseType::geo_,
77  BaseType::rock_comp_props_,
78  polymer_props_,
79  well_model,
80  BaseType::solver_,
81  BaseType::eclipse_state_,
82  BaseType::has_disgas_,
83  BaseType::has_vapoil_,
84  has_polymer_,
85  has_plyshlog_,
86  has_shrate_,
87  wells_rep_radius_,
88  wells_perf_length_,
89  wells_bore_diameter_,
90  BaseType::terminal_output_));
91 
92  if (!BaseType::threshold_pressures_by_face_.empty()) {
93  model->setThresholdPressures(BaseType::threshold_pressures_by_face_);
94  }
95 
96  return std::unique_ptr<Solver>(new Solver(BaseType::solver_param_, std::move(model)));
97  }
98 
99 
100 
101 
102  template <class GridT>
103  void SimulatorFullyImplicitBlackoilPolymer<GridT>::
104  handleAdditionalWellInflow(SimulatorTimer& timer,
105  WellsManager& wells_manager,
106  typename BaseType::WellState& well_state,
107  const Wells* wells)
108  {
109  // compute polymer inflow
110  std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
111  if (deck_->hasKeyword("WPOLYMER")) {
112  if (wells_manager.c_wells() == 0) {
113  OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
114  }
115  polymer_inflow_ptr.reset(new PolymerInflowFromDeck(*BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
116  } else {
117  OPM_MESSAGE("Warning: simulating with no WPOLYMER in deck (no polymer will be injected).");
118  polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
119  1.0*Opm::unit::day,
120  0.0));
121  }
122  std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_));
123  polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
124  timer.simulationTimeElapsed() + timer.currentStepLength(),
125  polymer_inflow_c);
126  well_state.polymerInflow() = polymer_inflow_c;
127 
128  if (has_plyshlog_) {
129  computeRepRadiusPerfLength(*BaseType::eclipse_state_, timer.currentStepNum(), BaseType::grid_, wells_rep_radius_, wells_perf_length_, wells_bore_diameter_);
130  }
131  }
132 
133 
134  template <class GridT>
135  void SimulatorFullyImplicitBlackoilPolymer<GridT>::
136  setupCompressedToCartesian(const int* global_cell, int number_of_cells,
137  std::map<int,int>& cartesian_to_compressed )
138  {
139  if (global_cell) {
140  for (int i = 0; i < number_of_cells; ++i) {
141  cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
142  }
143  }
144  else {
145  for (int i = 0; i < number_of_cells; ++i) {
146  cartesian_to_compressed.insert(std::make_pair(i, i));
147  }
148  }
149 
150  }
151 
152 
153  template <class GridT>
154  void SimulatorFullyImplicitBlackoilPolymer<GridT>::
155  computeRepRadiusPerfLength(const Opm::EclipseState& eclipseState,
156  const size_t timeStep,
157  const GridT& grid,
158  std::vector<double>& wells_rep_radius,
159  std::vector<double>& wells_perf_length,
160  std::vector<double>& wells_bore_diameter)
161  {
162 
163  // TODO, the function does not work for parallel running
164  // to be fixed later.
165  int number_of_cells = Opm::UgGridHelpers::numCells(grid);
166  const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
167  const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
168  auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
169  auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
170 
171  if (eclipseState.getSchedule().numWells() == 0) {
172  OPM_MESSAGE("No wells specified in Schedule section, "
173  "initializing no wells");
174  return;
175  }
176 
177  const size_t n_perf = wells_rep_radius.size();
178 
179  wells_rep_radius.clear();
180  wells_perf_length.clear();
181  wells_bore_diameter.clear();
182 
183  wells_rep_radius.reserve(n_perf);
184  wells_perf_length.reserve(n_perf);
185  wells_bore_diameter.reserve(n_perf);
186 
187  std::map<int,int> cartesian_to_compressed;
188 
189  setupCompressedToCartesian(global_cell, number_of_cells,
190  cartesian_to_compressed);
191 
192  const auto& schedule = eclipseState.getSchedule();
193  auto wells = schedule.getWells(timeStep);
194 
195  int well_index = 0;
196 
197  for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
198  const auto* well = (*wellIter);
199 
200  if (well->getStatus(timeStep) == WellCommon::SHUT) {
201  continue;
202  }
203  { // COMPDAT handling
204  const auto& completionSet = well->getCompletions(timeStep);
205  for (size_t c=0; c<completionSet.size(); c++) {
206  const auto& completion = completionSet.get(c);
207  if (completion.getState() == WellCompletion::OPEN) {
208  int i = completion.getI();
209  int j = completion.getJ();
210  int k = completion.getK();
211 
212  const int* cpgdim = cart_dims;
213  int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
214  std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
215  if (cgit == cartesian_to_compressed.end()) {
216  OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
217  << k << " not found in grid (well = " << well->name() << ')');
218  }
219  int cell = cgit->second;
220 
221  {
222  double radius = 0.5*completion.getDiameter();
223  if (radius <= 0.0) {
224  radius = 0.5*unit::feet;
225  OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
226  }
227 
228  const std::array<double, 3> cubical =
229  WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
230 
231  WellCompletion::DirectionEnum direction = completion.getDirection();
232 
233  double re; // area equivalent radius of the grid block
234  double perf_length; // the length of the well perforation
235 
236  switch (direction) {
237  case Opm::WellCompletion::DirectionEnum::X:
238  re = std::sqrt(cubical[1] * cubical[2] / M_PI);
239  perf_length = cubical[0];
240  break;
241  case Opm::WellCompletion::DirectionEnum::Y:
242  re = std::sqrt(cubical[0] * cubical[2] / M_PI);
243  perf_length = cubical[1];
244  break;
245  case Opm::WellCompletion::DirectionEnum::Z:
246  re = std::sqrt(cubical[0] * cubical[1] / M_PI);
247  perf_length = cubical[2];
248  break;
249  default:
250  OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
251  }
252 
253  double repR = std::sqrt(re * radius);
254  wells_rep_radius.push_back(repR);
255  wells_perf_length.push_back(perf_length);
256  wells_bore_diameter.push_back(2. * radius);
257  }
258  } else {
259  if (completion.getState() != WellCompletion::SHUT) {
260  OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled");
261  }
262  }
263 
264  }
265  }
266  well_index++;
267  }
268  }
269 
270 } // namespace Opm