SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
1 /*
2  Copyright 2013 SINTEF ICT, Applied Mathematics.
3  Copyright 2014 IRIS AS
4  Copyright 2015 Andreas Lauser
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 
23 namespace Opm
24 {
25 
26  template <class GridT>
27  auto SimulatorFullyImplicitBlackoilMultiSegment<GridT>::
28  createSolver(const WellModel& well_model)
29  -> std::unique_ptr<Solver>
30  {
31  typedef typename Traits::Model Model;
32 
33  auto model = std::unique_ptr<Model>(new Model(model_param_,
34  grid_,
35  props_,
36  geo_,
37  rock_comp_props_,
38  well_model,
39  solver_,
40  eclipse_state_,
41  has_disgas_,
42  has_vapoil_,
43  terminal_output_));
44 
45  if (!Base::threshold_pressures_by_face_.empty()) {
46  model->setThresholdPressures(Base::threshold_pressures_by_face_);
47  }
48 
49  return std::unique_ptr<ThisType::Solver>(new Solver(Base::solver_param_, std::move(model)));
50 
51  }
52 
53  template <class GridT>
54  SimulatorReport SimulatorFullyImplicitBlackoilMultiSegment<GridT>::run(SimulatorTimer& timer,
55  ReservoirState& state)
56  {
57  WellState prev_well_state;
58 
59  // Create timers and file for writing timing info.
60  Opm::time::StopWatch solver_timer;
61  double stime = 0.0;
62  Opm::time::StopWatch step_timer;
63  Opm::time::StopWatch total_timer;
64  total_timer.start();
65  std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
66  std::ofstream tstep_os(tstep_filename.c_str());
67 
68  // adaptive time stepping
69  const auto& events = eclipse_state_->getSchedule().getEvents();
70  std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
71  if( param_.getDefault("timestep.adaptive", true ) )
72  {
73  adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
74  }
75 
76  std::string restorefilename = param_.getDefault("restorefile", std::string("") );
77  if( ! restorefilename.empty() )
78  {
79  // -1 means that we'll take the last report step that was written
80  const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
81  output_writer_.restore( timer,
82  state,
83  prev_well_state,
84  restorefilename,
85  desiredRestoreStep );
86  }
87 
88  unsigned int totalNonlinearIterations = 0;
89  unsigned int totalLinearIterations = 0;
90  DynamicListEconLimited dynamic_list_econ_limited;
91 
92  bool ooip_computed = false;
93  std::vector<int> fipnum_global = eclipse_state_->get3DProperties().getIntGridProperty("FIPNUM").getData();
94  //Get compressed cell fipnum.
95  std::vector<int> fipnum(AutoDiffGrid::numCells(grid_));
96  if (fipnum_global.empty()) {
97  std::fill(fipnum.begin(), fipnum.end(), 0);
98  } else {
99  for (size_t c = 0; c < fipnum.size(); ++c) {
100  fipnum[c] = fipnum_global[AutoDiffGrid::globalCell(grid_)[c]];
101  }
102  }
103  std::vector<std::vector<double> > OOIP;
104 
105  // Main simulation loop.
106  while (!timer.done()) {
107  // Report timestep.
108  step_timer.start();
109  if ( terminal_output_ )
110  {
111  timer.report(std::cout);
112  }
113 
114  // Create wells and well state.
115  WellsManager wells_manager(*eclipse_state_,
116  timer.currentStepNum(),
117  Opm::UgGridHelpers::numCells(grid_),
118  Opm::UgGridHelpers::globalCell(grid_),
119  Opm::UgGridHelpers::cartDims(grid_),
120  Opm::UgGridHelpers::dimensions(grid_),
121  Opm::UgGridHelpers::cell2Faces(grid_),
122  Opm::UgGridHelpers::beginFaceCentroids(grid_),
123  dynamic_list_econ_limited,
124  is_parallel_run_,
125  // We need to pass the optionaly arguments
126  // as we get the following error otherwise
127  // with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
128  // converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructor
129  Base::defunct_well_names_);
130  const Wells* wells = wells_manager.c_wells();
131  WellState well_state;
132  // well_state.init(wells, state, prev_well_state);
133 
134  const auto wells_ecl = eclipse_state_->getSchedule().getWells(timer.currentStepNum());
135  const int current_time_step = timer.currentStepNum();
136 
137  const WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, current_time_step);
138 
139  well_state.init(well_model, state, prev_well_state, wells);
140 
141  // give the polymer and surfactant simulators the chance to do their stuff
142  Base::asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
143 
144  // write the inital state at the report stage
145  if (timer.initialStep()) {
146  // No per cell data is written for initial step, but will be
147  // for subsequent steps, when we have started simulating
148  output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state, {}, {} );
149  }
150 
151  // Max oil saturation (for VPPARS), hysteresis update.
152  props_.updateSatOilMax(state.saturation());
153  props_.updateSatHyst(state.saturation(), allcells_);
154 
155  // Compute reservoir volumes for RESV controls.
156  Base::asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state);
157 
158  // Run a multiple steps of the solver depending on the time step control.
159  solver_timer.start();
160 
161  auto solver = createSolver(well_model);
162 
163  // Compute orignal FIP;
164  if (!ooip_computed) {
165  OOIP = solver->computeFluidInPlace(state, fipnum);
166  Base::FIPUnitConvert(eclipse_state_->getUnits(), OOIP);
167  ooip_computed = true;
168  }
169 
170  // If sub stepping is enabled allow the solver to sub cycle
171  // in case the report steps are too large for the solver to converge
172  //
173  // \Note: The report steps are met in any case
174  // \Note: The sub stepping will require a copy of the state variables
175  if( adaptiveTimeStepping ) {
176  bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
177  events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
178  events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
179  events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
180  adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_);
181  }
182  else {
183  // solve for complete report step
184  solver->step(timer, state, well_state);
185  }
186 
187  // take time that was used to solve system for this reportStep
188  solver_timer.stop();
189 
190  // accumulate the number of nonlinear and linear Iterations
191  totalNonlinearIterations += solver->nonlinearIterations();
192  totalLinearIterations += solver->linearIterations();
193 
194  // Report timing.
195  const double st = solver_timer.secsSinceStart();
196 
197  // Compute current FIP.
198  std::vector<std::vector<double> > COIP;
199  COIP = solver->computeFluidInPlace(state, fipnum);
200  std::vector<double> OOIP_totals = Base::FIPTotals(OOIP, state);
201  std::vector<double> COIP_totals = Base::FIPTotals(COIP, state);
202 
203  //Convert to correct units
204  Base::FIPUnitConvert(eclipse_state_->getUnits(), COIP);
205  Base::FIPUnitConvert(eclipse_state_->getUnits(), OOIP_totals);
206  Base::FIPUnitConvert(eclipse_state_->getUnits(), COIP_totals);
207 
208  if ( terminal_output_ )
209  {
210  Base::outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
211  for (size_t reg = 0; reg < OOIP.size(); ++reg) {
212  Base::outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
213  }
214  }
215 
216  if ( terminal_output_ )
217  {
218  std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
219  }
220 
221  stime += st;
222  if ( output_writer_.output() ) {
223  SimulatorReport step_report;
224  step_report.solver_time = st;
225  step_report.total_time = step_timer.secsSinceStart();
226  step_report.reportParam(tstep_os);
227  }
228 
229  // Increment timer, remember well state.
230  ++timer;
231 
232 
233  // write simulation state at the report stage
234  const auto& physicalModel = solver->model();
235  output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
236 
237  prev_well_state = well_state;
238  }
239 
240  // Stop timer and create timing report
241  total_timer.stop();
242  SimulatorReport report;
243  report.total_time = total_timer.secsSinceStart();
244  report.solver_time = stime;
245  report.total_newton_iterations = totalNonlinearIterations;
246  report.total_linear_iterations = totalLinearIterations;
247  return report;
248  }
249 
250 } // namespace Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22