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