All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
SteadyStateUpscalerImplicit_impl.hpp
1 //===========================================================================
2 //
3 // File: SteadyStateUpscaler_impl.hpp
4 //
5 // Created: Fri Aug 28 14:07:51 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // Brd Skaflestad <bard.skaflestad@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 2009, 2010 Statoil ASA.
19 
20  This file is part of The Open Reservoir Simulator Project (OpenRS).
21 
22  OpenRS is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation, either version 3 of the License, or
25  (at your option) any later version.
26 
27  OpenRS is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 #ifndef OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
37 #define OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
38 
39 
40 #include <boost/lexical_cast.hpp>
41 #include <opm/porsol/common/MatrixInverse.hpp>
42 #include <opm/porsol/common/SimulatorUtilities.hpp>
43 #include <opm/porsol/common/ReservoirPropertyFixedMobility.hpp>
44 #include <opm/porsol/euler/MatchSaturatedVolumeFunctor.hpp>
45 #include <opm/parser/eclipse/Units/Units.hpp>
46 
47 #include <opm/upscaling/writeECLData.hpp>
48 #include <opm/core/utility/miscUtilities.hpp>
49 #include <algorithm>
50 #include <iostream>
51 
52 namespace Opm
53 {
54 
55 
56 
57  template <class Traits>
59  : Super(),
60  use_gravity_(false),
61  output_vtk_(false),
62  output_ecl_(false),
63  print_inoutflows_(false),
64  simulation_steps_(10),
65  init_stepsize_(10),
66  relperm_threshold_(1.0e-8),
67  maximum_mobility_contrast_(1.0e9),
68  sat_change_year_(1.0e-5),
69  max_it_(100),
70  max_stepsize_(1e4),
71  dt_sat_tol_(1e-2),
72  use_maxdiff_(true)
73  {
74  }
75 
76 
77 
78 
79  template <class Traits>
80  inline void SteadyStateUpscalerImplicit<Traits>::initImpl(const Opm::ParameterGroup& param)
81  {
82  Super::initImpl(param);
83  use_gravity_ = param.getDefault("use_gravity", use_gravity_);
84  output_vtk_ = param.getDefault("output_vtk", output_vtk_);
85  output_ecl_ = param.getDefault("output_ecl", output_ecl_);
86  if (output_ecl_) {
87  grid_adapter_.init(Super::grid());
88  }
89  print_inoutflows_ = param.getDefault("print_inoutflows", print_inoutflows_);
90  simulation_steps_ = param.getDefault("simulation_steps", simulation_steps_);
91  init_stepsize_ = Opm::unit::convert::from(param.getDefault("init_stepsize", init_stepsize_),
92  Opm::unit::day);
93  relperm_threshold_ = param.getDefault("relperm_threshold", relperm_threshold_);
94  maximum_mobility_contrast_ = param.getDefault("maximum_mobility_contrast", maximum_mobility_contrast_);
95  sat_change_year_ = param.getDefault("sat_change_year", sat_change_year_);
96  dt_sat_tol_ = param.getDefault("dt_sat_tol", dt_sat_tol_);
97  max_it_ = param.getDefault("max_it", max_it_);
98  max_stepsize_ = Opm::unit::convert::from(param.getDefault("max_stepsize", max_stepsize_),Opm::unit::year);
99  use_maxdiff_ = param.getDefault("use_maxdiff", use_maxdiff_);
100  transport_solver_.init(param);
101  // Set viscosities and densities if given.
102  double v1_default = this->res_prop_.viscosityFirstPhase();
103  double v2_default = this->res_prop_.viscositySecondPhase();
104  this->res_prop_.setViscosities(param.getDefault("viscosity1", v1_default), param.getDefault("viscosity2", v2_default));
105  double d1_default = this->res_prop_.densityFirstPhase();
106  double d2_default = this->res_prop_.densitySecondPhase();
107  this->res_prop_.setDensities(param.getDefault("density1", d1_default), param.getDefault("density2", d2_default));
108  }
109 
110 
111 
112 
113  namespace {
114  double maxMobility(double m1, double m2)
115  {
116  return std::max(m1, m2);
117  }
118  // The matrix variant expects diagonal mobilities.
119  template <class SomeMatrixType>
120  double maxMobility(double m1, SomeMatrixType& m2)
121  {
122  double m = m1;
123  for (int i = 0; i < std::min(m2.numRows(), m2.numCols()); ++i) {
124  m = std::max(m, m2(i,i));
125  }
126  return m;
127  }
128  void thresholdMobility(double& m, double threshold)
129  {
130  m = std::max(m, threshold);
131  }
132  // The matrix variant expects diagonal mobilities.
133  template <class SomeMatrixType>
134  void thresholdMobility(SomeMatrixType& m, double threshold)
135  {
136  for (int i = 0; i < std::min(m.numRows(), m.numCols()); ++i) {
137  m(i, i) = std::max(m(i, i), threshold);
138  }
139  }
140 
141  // Workaround: duplicate of Opm::destripe in opm-simulators.
142  inline std::vector< double > destripe( const std::vector< double >& v,
143  size_t stride,
144  size_t offset )
145  {
146  std::vector< double > dst( v.size() / stride );
147  size_t di = 0;
148  for( size_t i = offset; i < v.size(); i += stride ) {
149  dst[ di++ ] = v[ i ];
150  }
151  return dst;
152  }
153 
154  } // anon namespace
155 
156 
157 
158 
159  template <class Traits>
160  inline std::pair<typename SteadyStateUpscalerImplicit<Traits>::permtensor_t,
161  typename SteadyStateUpscalerImplicit<Traits>::permtensor_t>
163  upscaleSteadyState(const int flow_direction,
164  const std::vector<double>& initial_saturation,
165  const double boundary_saturation,
166  const double pressure_drop,
167  const permtensor_t& upscaled_perm,
168  bool& success)
169  {
170  static int count = 0;
171  ++count;
172  int num_cells = this->ginterf_.numberOfCells();
173  // No source or sink.
174  std::vector<double> src(num_cells, 0.0);
175  Opm::SparseVector<double> injection(num_cells);
176  // Gravity.
177  Dune::FieldVector<double, 3> gravity(0.0);
178  if (use_gravity_) {
179  gravity[2] = Opm::unit::gravity;
180  }
181 
182  if (gravity.two_norm() > 0.0) {
183  OPM_MESSAGE("Warning: Gravity is experimental for flow solver.");
184  }
185 
186  // Put pore volume in vector.
187  std::vector<double> pore_vol;
188  pore_vol.reserve(num_cells);
189  double tot_pore_vol = 0.0;
190  typedef typename GridInterface::CellIterator CellIter;
191  for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
192  double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
193  pore_vol.push_back(cell_pore_vol);
194  tot_pore_vol += cell_pore_vol;
195  }
196 
197  // Set up initial saturation profile.
198  std::vector<double> saturation = initial_saturation;
199 
200  // Set up boundary conditions.
201  setupUpscalingConditions(this->ginterf_, this->bctype_, flow_direction,
202  pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
203 
204  // Set up solvers.
205  if (flow_direction == 0) {
206  this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
207  }
208  transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
209 
210  // Run pressure solver.
211  this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
212  this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_);
213  double max_mod = this->flow_solver_.postProcessFluxes();
214  std::cout << "Max mod = " << max_mod << std::endl;
215 
216  // Do a run till steady state.
217  std::vector<double> saturation_old = saturation;
218  bool stationary = false;
219  int it_count = 0;
220  double stepsize = init_stepsize_;
221  double ecl_time = 0.0;
222  std::vector<double> ecl_sat;
223  std::vector<double> ecl_press;
224  std::vector<double> init_saturation(saturation);
225  while ((!stationary) && (it_count < max_it_)) { // && transport_cost < max_transport_cost_)
226  // Run transport solver.
227  std::cout << "Running transport step " << it_count << " with stepsize "
228  << stepsize/Opm::unit::year << " years." << std::endl;
229  bool converged = transport_solver_.transportSolve(saturation, stepsize, gravity,
230  this->flow_solver_.getSolution(), injection);
231  // Run pressure solver.
232  if (converged) {
233  init_saturation = saturation;
234  // this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
235  // this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_);
236  // max_mod = this->flow_solver_.postProcessFluxes();
237  // std::cout << "Max mod of fluxes= " << max_mod << std::endl;
238  this->flow_solver_.postProcessFluxes();
239  // Print in-out flows if requested.
240  if (print_inoutflows_) {
241  std::pair<double, double> w_io, o_io;
242  computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
243  std::cout << "Pressure step " << it_count
244  << "\nWater flow [in] " << w_io.first
245  << " [out] " << w_io.second
246  << "\nOil flow [in] " << o_io.first
247  << " [out] " << o_io.second
248  << std::endl;
249  }
250  // Output.
251  if (output_vtk_) {
252  writeVtkOutput(this->ginterf_,
253  this->res_prop_,
254  this->flow_solver_.getSolution(),
255  saturation,
256  std::string("output-steadystate")
257  + '-' + boost::lexical_cast<std::string>(count)
258  + '-' + boost::lexical_cast<std::string>(flow_direction)
259  + '-' + boost::lexical_cast<std::string>(it_count));
260  }
261  if (output_ecl_) {
262  const char* fd = "xyz";
263  std::string basename = std::string("ecldump-steadystate")
264  + '-' + boost::lexical_cast<std::string>(boundary_saturation)
265  + '-' + boost::lexical_cast<std::string>(fd[flow_direction])
266  + '-' + boost::lexical_cast<std::string>(pressure_drop);
267  Opm::toBothSat(saturation, ecl_sat);
268  getCellPressure(ecl_press, this->ginterf_, this->flow_solver_.getSolution());
269  data::Solution solution;
270  solution.insert( "PRESSURE" , UnitSystem::measure::pressure , ecl_press , data::TargetType::RESTART_SOLUTION );
271  solution.insert( "SWAT" , UnitSystem::measure::identity , destripe( ecl_sat, 2, 0) , data::TargetType::RESTART_SOLUTION );
272  ecl_time += stepsize;
273  boost::posix_time::ptime ecl_startdate( boost::gregorian::date(2012, 1, 1) );
274  boost::posix_time::ptime ecl_curdate = ecl_startdate + boost::posix_time::seconds(int(ecl_time));
275  boost::posix_time::ptime epoch( boost::gregorian::date( 1970, 1, 1 ) );
276  auto ecl_posix_time = ( ecl_curdate - epoch ).total_seconds();
277  const auto* cgrid = grid_adapter_.c_grid();
278  Opm::writeECLData(cgrid->cartdims[ 0 ],
279  cgrid->cartdims[ 1 ],
280  cgrid->cartdims[ 2 ],
281  cgrid->number_of_cells,
282  solution, it_count,
283  ecl_time, ecl_posix_time,
284  "./", basename);
285  }
286  // Comparing old to new.
287  double maxdiff = 0.0;
288  double euclidean_diff = 0.0;
289  for (int i = 0; i < num_cells; ++i) {
290  const double sat_diff_cell = saturation[i] - saturation_old[i];
291  maxdiff = std::max(maxdiff, std::fabs(sat_diff_cell));
292  euclidean_diff += sat_diff_cell * sat_diff_cell * pore_vol[i];
293  }
294  euclidean_diff = std::sqrt(euclidean_diff / tot_pore_vol);
295  double ds_year;
296  if (use_maxdiff_) {
297  ds_year = maxdiff*Opm::unit::year/stepsize;
298  std::cout << "Maximum saturation change/year: " << ds_year << std::endl;
299  if (maxdiff < dt_sat_tol_) {
300  stepsize=std::min(max_stepsize_,2*stepsize);
301  }
302  }
303  else {
304  ds_year = euclidean_diff*Opm::unit::year/stepsize;
305  std::cout << "Euclidean saturation change/year: " << ds_year << std::endl;
306  if (euclidean_diff < dt_sat_tol_) {
307  stepsize=std::min(max_stepsize_,2*stepsize);
308  }
309  }
310  if (ds_year < sat_change_year_) {
311  stationary = true;
312  }
313  } else {
314  std::cerr << "Cutting time step\n";
315  init_saturation = saturation_old;
316  stepsize=stepsize/2.0;
317  }
318  ++it_count;
319  // Copy to old.
320  saturation_old = saturation;
321  }
322  success = stationary;
323 
324  // Compute phase mobilities.
325  // First: compute maximal mobilities.
326  typedef typename Super::ResProp::Mobility Mob;
327  Mob m;
328  double m1max = 0;
329  double m2max = 0;
330  for (int c = 0; c < num_cells; ++c) {
331  this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
332  m1max = maxMobility(m1max, m.mob);
333  this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
334  m2max = maxMobility(m2max, m.mob);
335  }
336  // Second: set thresholds.
337  const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
338  const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
339  const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
340  const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
341  const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
342  const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
343  // Third: extract and threshold.
344  std::vector<Mob> mob1(num_cells);
345  std::vector<Mob> mob2(num_cells);
346  for (int c = 0; c < num_cells; ++c) {
347  this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
348  thresholdMobility(mob1[c].mob, mob1_threshold);
349  this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
350  thresholdMobility(mob2[c].mob, mob2_threshold);
351  }
352 
353  // Compute upscaled relperm for each phase.
354  ReservoirPropertyFixedMobility<Mob> fluid_first(mob1);
355  permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
356  ReservoirPropertyFixedMobility<Mob> fluid_second(mob2);
357  permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
358 
359  // Set the steady state saturation fields for eventual outside access.
360  last_saturation_state_.swap(saturation);
361 
362  // Compute the (anisotropic) upscaled mobilities.
363  // eff_Kw := lambda_w*K
364  // => lambda_w = eff_Kw*inv(K);
365  permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
366  permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
367 
368  // Compute (anisotropic) upscaled relative permeabilities.
369  // lambda = k_r/mu
370  permtensor_t k_rw(lambda_w);
371  k_rw *= this->res_prop_.viscosityFirstPhase();
372  permtensor_t k_ro(lambda_o);
373  k_ro *= this->res_prop_.viscositySecondPhase();
374  return std::make_pair(k_rw, k_ro);
375  }
376 
377 
378 
379 
380  template <class Traits>
381  inline const std::vector<double>&
383  {
384  return last_saturation_state_;
385  }
386 
387 
388 
389 
390  template <class Traits>
392  {
393  typedef typename GridInterface::CellIterator CellIter;
394  double pore_vol = 0.0;
395  double sat_vol = 0.0;
396  for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
397  double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
398  pore_vol += cell_pore_vol;
399  sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
400  }
401  // Dividing by pore volume gives average saturations.
402  return sat_vol/pore_vol;
403  }
404 
405 
406  template <class Traits>
407  void SteadyStateUpscalerImplicit<Traits>::initSatLimits(std::vector<double>& s) const
408  {
409  for (int c = 0; c < int (s.size()); c++ ) {
410  double s_min = this->res_prop_.s_min(c);
411  double s_max = this->res_prop_.s_max(c);
412  s[c] = std::max(s_min+1e-4, s[c]);
413  s[c] = std::min(s_max-1e-4, s[c]);
414  }
415  }
416 
417  template <class Traits>
418  void SteadyStateUpscalerImplicit<Traits>::setToCapillaryLimit(double average_s, std::vector<double>& s) const
419  {
420  int num_cells = this->ginterf_.numberOfCells();
421  std::vector<double> s_orig(num_cells, average_s);
422 
423  initSatLimits(s_orig);
424  std::vector<double> cap_press(num_cells, 0.0);
425  typedef typename UpscalerBase<Traits>::ResProp ResProp;
426  MatchSaturatedVolumeFunctor<GridInterface, ResProp> func(this->ginterf_, this->res_prop_, s_orig, cap_press);
427  double cap_press_range = 1e2;
428  double mod_low = 1e100;
429  double mod_high = -1e100;
430  Opm::bracketZero(func, 0.0, cap_press_range, mod_low, mod_high);
431  const int max_iter = 40;
432  const double nonlinear_tolerance = 1e-12;
433  int iterations_used = -1;
434  typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
435  double mod_correct = RootFinder::solve(func, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
436  std::cout << "Moved capillary pressure solution by " << mod_correct << " after "
437  << iterations_used << " iterations." << std::endl;
438  s = func.lastSaturations();
439  }
440 
441 
442  template <class Traits>
443  template <class FlowSol>
444  void SteadyStateUpscalerImplicit<Traits>::computeInOutFlows(std::pair<double, double>& water_inout,
445  std::pair<double, double>& oil_inout,
446  const FlowSol& flow_solution,
447  const std::vector<double>& saturations) const
448  {
449  typedef typename GridInterface::CellIterator CellIter;
450  typedef typename CellIter::FaceIterator FaceIter;
451 
452  double side1_flux = 0.0;
453  double side2_flux = 0.0;
454  double side1_flux_oil = 0.0;
455  double side2_flux_oil = 0.0;
456  std::map<int, double> frac_flow_by_bid;
457  int num_cells = this->ginterf_.numberOfCells();
458  std::vector<double> cell_inflows_w(num_cells, 0.0);
459  std::vector<double> cell_outflows_w(num_cells, 0.0);
460 
461  // Two passes: First pass, deal with outflow, second pass, deal with inflow.
462  // This is for the periodic case, so that we are sure all fractional flows have
463  // been set in frac_flow_by_bid.
464  for (int pass = 0; pass < 2; ++pass) {
465  for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
466  for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
467  if (f->boundary()) {
468  double flux = flow_solution.outflux(f);
469  const SatBC& sc = this->bcond_.satCond(f);
470  if (flux < 0.0 && pass == 1) {
471  // This is an inflow face.
472  double frac_flow = 0.0;
473  if (sc.isPeriodic()) {
474  assert(sc.saturationDifference() == 0.0);
475  int partner_bid = this->bcond_.getPeriodicPartner(f->boundaryId());
476  std::map<int, double>::const_iterator it = frac_flow_by_bid.find(partner_bid);
477  if (it == frac_flow_by_bid.end()) {
478  OPM_THROW(std::runtime_error, "Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
479  << " and partner bid = " << partner_bid);
480  }
481  frac_flow = it->second;
482  } else {
483  assert(sc.isDirichlet());
484  frac_flow = this->res_prop_.fractionalFlow(c->index(), sc.saturation());
485  }
486  cell_inflows_w[c->index()] += flux*frac_flow;
487  side1_flux += flux*frac_flow;
488  side1_flux_oil += flux*(1.0 - frac_flow);
489  } else if (flux >= 0.0 && pass == 0) {
490  // This is an outflow face.
491  double frac_flow = this->res_prop_.fractionalFlow(c->index(), saturations[c->index()]);
492  if (sc.isPeriodic()) {
493  frac_flow_by_bid[f->boundaryId()] = frac_flow;
494  // std::cout << "Inserted bid " << f->boundaryId() << std::endl;
495  }
496  cell_outflows_w[c->index()] += flux*frac_flow;
497  side2_flux += flux*frac_flow;
498  side2_flux_oil += flux*(1.0 - frac_flow);
499  }
500  }
501  }
502  }
503  }
504  water_inout = std::make_pair(side1_flux, side2_flux);
505  oil_inout = std::make_pair(side1_flux_oil, side2_flux_oil);
506  }
507 
508 
509 
510 
511 } // namespace Opm
512 
513 
514 #endif // OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
const std::vector< double > & lastSaturationState() const
Accessor for the steady state saturation field.
Definition: SteadyStateUpscalerImplicit_impl.hpp:382
void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition: setupBoundaryConditions.hpp:99
Definition: MatchSaturatedVolumeFunctor.hpp:68
Definition: GridInterfaceEuler.hpp:376
A base class for upscaling.
Definition: UpscalerBase.hpp:55
void initSatLimits(std::vector< double > &s) const
Ensure saturations are not outside table.
Definition: SteadyStateUpscalerImplicit_impl.hpp:407
void getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition: SimulatorUtilities.hpp:184
SteadyStateUpscalerImplicit()
Default constructor.
Definition: SteadyStateUpscalerImplicit_impl.hpp:58
double lastSaturationUpscaled() const
Computes the upscaled saturation corresponding to the saturation field returned by lastSaturationStat...
Definition: SteadyStateUpscalerImplicit_impl.hpp:391
std::pair< permtensor_t, permtensor_t > upscaleSteadyState(const int flow_direction, const std::vector< double > &initial_saturation, const double boundary_saturation, const double pressure_drop, const permtensor_t &upscaled_perm, bool &success)
Does a steady-state upscaling.
Definition: SteadyStateUpscalerImplicit_impl.hpp:163
virtual void initImpl(const Opm::ParameterGroup &param)
Override from superclass.
Definition: SteadyStateUpscalerImplicit_impl.hpp:80
A class for doing steady state upscaling.
Definition: SteadyStateUpscalerImplicit.hpp:52
Definition: ReservoirPropertyFixedMobility.hpp:46