My Project
SteadyStateUpscaler_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_STEADYSTATEUPSCALER_IMPL_HEADER
37 #define OPM_STEADYSTATEUPSCALER_IMPL_HEADER
38 
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/parser/eclipse/Units/Units.hpp>
44 #include <algorithm>
45 #include <iostream>
46 
47 namespace Opm
48 {
49 
50 
51 
52  template <class Traits>
54  : Super(),
55  use_gravity_(false),
56  output_vtk_(false),
57  print_inoutflows_(false),
58  simulation_steps_(10),
59  stepsize_(0.1),
60  relperm_threshold_(1.0e-8),
61  maximum_mobility_contrast_(1.0e9),
62  sat_change_threshold_(0.0)
63  {
64  }
65 
66 
67 
68 
69  template <class Traits>
70  inline void SteadyStateUpscaler<Traits>::initImpl(const Opm::parameter::ParameterGroup& param)
71  {
72  Super::initImpl(param);
73  use_gravity_ = param.getDefault("use_gravity", use_gravity_);
74  output_vtk_ = param.getDefault("output_vtk", output_vtk_);
75  print_inoutflows_ = param.getDefault("print_inoutflows", print_inoutflows_);
76  simulation_steps_ = param.getDefault("simulation_steps", simulation_steps_);
77  stepsize_ = Opm::unit::convert::from(param.getDefault("stepsize", stepsize_),
78  Opm::unit::day);
79  relperm_threshold_ = param.getDefault("relperm_threshold", relperm_threshold_);
80  maximum_mobility_contrast_ = param.getDefault("maximum_mobility_contrast", maximum_mobility_contrast_);
81  sat_change_threshold_ = param.getDefault("sat_change_threshold", sat_change_threshold_);
82 
83  transport_solver_.init(param);
84  // Set viscosities and densities if given.
85  double v1_default = this->res_prop_.viscosityFirstPhase();
86  double v2_default = this->res_prop_.viscositySecondPhase();
87  this->res_prop_.setViscosities(param.getDefault("viscosity1", v1_default), param.getDefault("viscosity2", v2_default));
88  double d1_default = this->res_prop_.densityFirstPhase();
89  double d2_default = this->res_prop_.densitySecondPhase();
90  this->res_prop_.setDensities(param.getDefault("density1", d1_default), param.getDefault("density2", d2_default));
91  }
92 
93 
94 
95 
96  namespace {
97  double maxMobility(double m1, double m2)
98  {
99  return std::max(m1, m2);
100  }
101  // The matrix variant expects diagonal mobilities.
102  template <class SomeMatrixType>
103  double maxMobility(double m1, SomeMatrixType& m2)
104  {
105  double m = m1;
106  for (int i = 0; i < std::min(m2.numRows(), m2.numCols()); ++i) {
107  m = std::max(m, m2(i,i));
108  }
109  return m;
110  }
111  void thresholdMobility(double& m, double threshold)
112  {
113  m = std::max(m, threshold);
114  }
115  // The matrix variant expects diagonal mobilities.
116  template <class SomeMatrixType>
117  void thresholdMobility(SomeMatrixType& m, double threshold)
118  {
119  for (int i = 0; i < std::min(m.numRows(), m.numCols()); ++i) {
120  m(i, i) = std::max(m(i, i), threshold);
121  }
122  }
123  } // anon namespace
124 
125 
126 
127 
128  template <class Traits>
129  inline std::pair<typename SteadyStateUpscaler<Traits>::permtensor_t,
130  typename SteadyStateUpscaler<Traits>::permtensor_t>
132  upscaleSteadyState(const int flow_direction,
133  const std::vector<double>& initial_saturation,
134  const double boundary_saturation,
135  const double pressure_drop,
136  const permtensor_t& upscaled_perm)
137  {
138  static int count = 0;
139  ++count;
140  int num_cells = this->ginterf_.numberOfCells();
141  // No source or sink.
142  std::vector<double> src(num_cells, 0.0);
143  Opm::SparseVector<double> injection(num_cells);
144  // Gravity.
145  Dune::FieldVector<double, 3> gravity(0.0);
146  if (use_gravity_) {
147  gravity[2] = Opm::unit::gravity;
148  }
149  if (gravity.two_norm() > 0.0) {
150  OPM_MESSAGE("Warning: Gravity is experimental for flow solver.");
151  }
152 
153  // Set up initial saturation profile.
154  std::vector<double> saturation = initial_saturation;
155 
156  // Set up boundary conditions.
157  setupUpscalingConditions(this->ginterf_, this->bctype_, flow_direction,
158  pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
159 
160  // Set up solvers.
161  if (flow_direction == 0) {
162  this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
163  }
164  transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
165 
166  // Run pressure solver.
167  this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
168  this->residual_tolerance_, this->linsolver_verbosity_,
169  this->linsolver_type_, false,
170  this->linsolver_maxit_, this->linsolver_prolongate_factor_,
171  this->linsolver_smooth_steps_);
172  double max_mod = this->flow_solver_.postProcessFluxes();
173  std::cout << "Max mod = " << max_mod << std::endl;
174 
175  // Do a run till steady state. For now, we just do some pressure and transport steps...
176  std::vector<double> saturation_old = saturation;
177  for (int iter = 0; iter < simulation_steps_; ++iter) {
178  // Run transport solver.
179  transport_solver_.transportSolve(saturation, stepsize_, gravity, this->flow_solver_.getSolution(), injection);
180 
181  // Run pressure solver.
182  this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
183  this->residual_tolerance_, this->linsolver_verbosity_,
184  this->linsolver_type_, false,
185  this->linsolver_maxit_, this->linsolver_prolongate_factor_,
186  this->linsolver_smooth_steps_);
187  max_mod = this->flow_solver_.postProcessFluxes();
188  std::cout << "Max mod = " << max_mod << std::endl;
189 
190  // Print in-out flows if requested.
191  if (print_inoutflows_) {
192  std::pair<double, double> w_io, o_io;
193  computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
194  std::cout << "Pressure step " << iter
195  << "\nWater flow [in] " << w_io.first
196  << " [out] " << w_io.second
197  << "\nOil flow [in] " << o_io.first
198  << " [out] " << o_io.second
199  << std::endl;
200  }
201 
202  // Output.
203  if (output_vtk_) {
204  writeVtkOutput(this->ginterf_,
205  this->res_prop_,
206  this->flow_solver_.getSolution(),
207  saturation,
208  std::string("output-steadystate")
209  + '-' + std::to_string(count)
210  + '-' + std::to_string(flow_direction)
211  + '-' + std::to_string(iter));
212  }
213 
214  // Comparing old to new.
215  int num_cells = saturation.size();
216  double maxdiff = 0.0;
217  for (int i = 0; i < num_cells; ++i) {
218  maxdiff = std::max(maxdiff, std::fabs(saturation[i] - saturation_old[i]));
219  }
220 #ifdef VERBOSE
221  std::cout << "Maximum saturation change: " << maxdiff << std::endl;
222 #endif
223  if (maxdiff < sat_change_threshold_) {
224 #ifdef VERBOSE
225  std::cout << "Maximum saturation change is under steady state threshold." << std::endl;
226 #endif
227  break;
228  }
229 
230  // Copy to old.
231  saturation_old = saturation;
232  }
233 
234  // Compute phase mobilities.
235  // First: compute maximal mobilities.
236  typedef typename Super::ResProp::Mobility Mob;
237  Mob m;
238  double m1max = 0;
239  double m2max = 0;
240  for (int c = 0; c < num_cells; ++c) {
241  this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
242  m1max = maxMobility(m1max, m.mob);
243  this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
244  m2max = maxMobility(m2max, m.mob);
245  }
246  // Second: set thresholds.
247  const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
248  const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
249  const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
250  const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
251  const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
252  const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
253  // Third: extract and threshold.
254  std::vector<Mob> mob1(num_cells);
255  std::vector<Mob> mob2(num_cells);
256  for (int c = 0; c < num_cells; ++c) {
257  this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
258  thresholdMobility(mob1[c].mob, mob1_threshold);
259  this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
260  thresholdMobility(mob2[c].mob, mob2_threshold);
261  }
262 
263  // Compute upscaled relperm for each phase.
264  ReservoirPropertyFixedMobility<Mob> fluid_first(mob1);
265  permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
266  ReservoirPropertyFixedMobility<Mob> fluid_second(mob2);
267  permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
268 
269  // Set the steady state saturation fields for eventual outside access.
270  last_saturation_state_.swap(saturation);
271 
272  // Compute the (anisotropic) upscaled mobilities.
273  // eff_Kw := lambda_w*K
274  // => lambda_w = eff_Kw*inv(K);
275  permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
276  permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
277 
278  // Compute (anisotropic) upscaled relative permeabilities.
279  // lambda = k_r/mu
280  permtensor_t k_rw(lambda_w);
281  k_rw *= this->res_prop_.viscosityFirstPhase();
282  permtensor_t k_ro(lambda_o);
283  k_ro *= this->res_prop_.viscositySecondPhase();
284  return std::make_pair(k_rw, k_ro);
285  }
286 
287 
288 
289 
290  template <class Traits>
291  inline const std::vector<double>&
293  {
294  return last_saturation_state_;
295  }
296 
297 
298 
299 
300  template <class Traits>
302  {
303  typedef typename GridInterface::CellIterator CellIter;
304  double pore_vol = 0.0;
305  double sat_vol = 0.0;
306  for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
307  double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
308  pore_vol += cell_pore_vol;
309  sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
310  }
311  // Dividing by pore volume gives average saturations.
312  return sat_vol/pore_vol;
313  }
314 
315 
316  static inline double calc_frac_flow (
317  std::map< int, double >& frac_flow_by_bid,
318  const BCs& bcond,
319  const ResProp& res_prop,
320  GridInterface::CellIterator c,
321  CellIterator::FaceIterator f ) {
322 
323  const SatBC& sc = bcond.satCond( f );
324 
325  if( sc.isPeriodic() ) {
326  assert( sc.isDirichlet() );
327  return res_prop.fractionalFlow( c->index(), sc.saturation() );
328  }
329 
330  assert( sc.saturationDifference() == 0.0 );
331  const int partner_bid = bcond.getPeriodicPartner( f->boundaryId() );
332  const auto it = frac_flow_by_bid.find( partner_bid );
333 
334  // This is for the periodic case, so that we are sure all fractional flows have
335  // been set in frac_flow_by_bid.
336  if( it == frac_flow_by_bid.end() ) {
337  OPM_THROW(std::runtime_error, "Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
338  << " and partner bid = " << partner_bid);
339  }
340 
341  return it->second;
342  }
343 
344  template <class Traits>
345  template <class FlowSol>
346  void SteadyStateUpscaler<Traits>::computeInOutFlows(std::pair<double, double>& water_inout,
347  std::pair<double, double>& oil_inout,
348  const FlowSol& flow_solution,
349  const std::vector<double>& saturations) const
350  {
351 
352  struct sideflux {
353  double water;
354  double oil;
355 
356  sideflux& operator+=( const sideflux& other ) {
357  this->water += other.water;
358  this->oil += other.oil;
359  return *this;
360  }
361  };
362 
363  std::map< int, double > frac_flow_by_bid;
364  // Two passes: First pass, deal with outflow
365  sideflux outflow;
366  for( auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
367  const double frac_flow = this->res_prop_.fractional_flow( c->index(), saturations[ c->index() ] );
368 
369  for( auto f = c->facebegin(); f != c->faceend(); ++f ) {
370  if( !f->boundary() ) continue;
371 
372  const double flux = flow_solution.outflux( f );
373  if( flux < 0.0 ) continue;
374 
375  const SatBC& sc = this->bcond_.satCond( f );
376  if( sc.isPeriodic() ) frac_flow_by_bid[ f->boundaryId() ] = frac_flow;
377 
378  outflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
379  }
380  }
381 
382  // second pass, deal with inflow.
383  sideflux inflow;
384  for( auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
385  for( auto f = c->facebegin(); f != c->faceend(); ++f ) {
386  if( !f->boundary() ) continue;
387 
388  const double flux = flow_solution.outflux( f );
389  if( flux >= 0.0 ) continue;
390 
391  const double frac_flow = calc_frac_flow( frac_flow_by_bid, this->bcond_, this->res_prop_, c, f );
392 
393  inflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
394  }
395  }
396 
397  water_inout = std::make_pair( in.water, out.water );
398  oil_inout = std::make_pair( in.oil, out.oil );
399 
400  }
401 } // namespace Opm
402 
403 
404 #endif // OPM_STEADYSTATEUPSCALER_IMPL_HEADER
bool isPeriodic() const
Type query.
Definition: BoundaryConditions.hpp:97
bool isDirichlet() const
Type query.
Definition: BoundaryConditions.hpp:85
Definition: GridInterfaceEuler.hpp:367
Definition: ReservoirPropertyFixedMobility.hpp:46
A class for representing a saturation boundary condition.
Definition: BoundaryConditions.hpp:176
double saturation() const
Query a Dirichlet condition.
Definition: BoundaryConditions.hpp:197
double saturationDifference() const
Query a Periodic condition.
Definition: BoundaryConditions.hpp:204
virtual void initImpl(const Opm::parameter::ParameterGroup &param)
Override from superclass.
Definition: SteadyStateUpscaler_impl.hpp:70
const std::vector< double > & lastSaturationState() const
Accessor for the steady state saturation field.
Definition: SteadyStateUpscaler_impl.hpp:292
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)
Does a steady-state upscaling.
Definition: SteadyStateUpscaler_impl.hpp:132
double lastSaturationUpscaled() const
Computes the upscaled saturation corresponding to the saturation field returned by lastSaturationStat...
Definition: SteadyStateUpscaler_impl.hpp:301
SteadyStateUpscaler()
Default constructor.
Definition: SteadyStateUpscaler_impl.hpp:53
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