All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
EulerUpstreamResidual_impl.hpp
1 //===========================================================================
2 //
3 // File: EulerUpstreamResidual_impl.hpp
4 //
5 // Created: Thu May 6 11:22:04 2010
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // Jostein R Natvig <jostein.r.natvig@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 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 OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
37 #define OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
38 
39 
40 
41 
42 //#include <cassert>
43 //#include <cmath>
44 //#include <algorithm>
45 
46 //#include <opm/common/ErrorMacros.hpp>
47 #include <opm/core/utility/Average.hpp>
48 //#include <opm/parser/eclipse/Units/Units.hpp>
49 //#include <dune/grid/common/Volumes.hpp>
50 
51 #include <opm/porsol/common/Matrix.hpp>
52 
53 #ifdef USE_TBB
54 #include <tbb/parallel_for.h>
55 #endif
56 
57 #include <iostream>
58 
59 
60 namespace Opm
61 {
62 
63 
64 
65 
66  namespace EulerUpstreamResidualDetails
67  {
68  template <typename T, template <typename> class StoragePolicy, class OrderingPolicy>
69  FullMatrix<T, OwnData, OrderingPolicy>
70  arithAver(const FullMatrix<T, StoragePolicy, OrderingPolicy>& m1,
71  const FullMatrix<T, StoragePolicy, OrderingPolicy>& m2)
72  {
73  return Opm::utils::arithmeticAverage<FullMatrix<T, StoragePolicy, OrderingPolicy>,
74  FullMatrix<T, OwnData, OrderingPolicy> >(m1, m2);
75  }
76 
77  template <class UpstreamSolver, class PressureSolution>
78  struct UpdateForCell
79  {
80  typedef typename UpstreamSolver::Vector Vector;
81  typedef typename UpstreamSolver::FIt FIt;
82  typedef typename UpstreamSolver::RP::PermTensor PermTensor;
83  typedef typename UpstreamSolver::RP::MutablePermTensor MutablePermTensor;
84 
85  const UpstreamSolver& s;
86  const std::vector<double>& saturation;
87  const Vector& gravity;
88  const PressureSolution& pressure_sol;
89  std::vector<double>& residual;
90 
91  UpdateForCell(const UpstreamSolver& solver,
92  const std::vector<double>& sat,
93  const Vector& grav,
94  const PressureSolution& psol,
95  std::vector<double>& res)
96  : s(solver), saturation(sat), gravity(grav), pressure_sol(psol), residual(res)
97  {
98  }
99 
100  template <class CIt>
101  void operator()(const CIt& c) const
102  {
103  // This is constant for the whole run.
104  const double delta_rho = s.preservoir_properties_->densityDifference();
105  int cell[2];
106  double cell_sat[2];
107  cell[0] = c->index();
108  cell_sat[0] = saturation[cell[0]];
109 
110  // Loop over all cell faces.
111  for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
112  // Neighbour face, will be changed if on a periodic boundary.
113  FIt nbface = f;
114  double dS = 0.0;
115  // Compute cell[1], cell_sat[1]
116  if (f->boundary()) {
117  if (s.pboundary_->satCond(*f).isPeriodic()) {
118  nbface = s.bid_to_face_[s.pboundary_->getPeriodicPartner(f->boundaryId())];
119  assert(nbface != f);
120  cell[1] = nbface->cellIndex();
121  assert(cell[0] != cell[1]);
122  // Periodic faces will be visited twice, but only once
123  // should they contribute. We make sure that we skip the
124  // periodic faces half the time.
125  if (cell[0] > cell[1]) {
126  // We skip this face.
127  continue;
128  }
129  cell_sat[1] = saturation[cell[1]];
130  } else {
131  assert(s.pboundary_->satCond(*f).isDirichlet());
132  cell[1] = cell[0];
133  cell_sat[1] = s.pboundary_->satCond(*f).saturation();
134  }
135  } else {
136  cell[1] = f->neighbourCellIndex();
137  assert(cell[0] != cell[1]);
138  if (cell[0] > cell[1]) {
139  // We skip this face.
140  continue;
141  }
142  cell_sat[1] = saturation[cell[1]];
143  }
144 
145  // Get some local properties.
146  const double loc_area = f->area();
147  const double loc_flux = pressure_sol.outflux(f);
148  const Vector loc_normal = f->normal();
149 
150  // We will now try to establish the upstream directions for each
151  // phase. They may be the same, or different (due to gravity).
152  // Recall the equation for v_w (water phase velocity):
153  // v_w = lambda_w * (lambda_o + lambda_w)^{-1}
154  // * (v + lambda_o * K * grad p_{cow} + lambda_o * K * (rho_w - rho_o) * g)
155  // ^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
156  // viscous term capillary term gravity term
157  //
158  // For the purpose of upstream weighting, we only consider the viscous and gravity terms.
159  // The question is, in which direction does v_w and v_o point? That is, what is the sign
160  // of v_w*loc_normal and v_o*loc_normal?
161  //
162  // For the case when the mobilities are scalar, the following analysis applies:
163  // The viscous contribution to v_w is loc_area*loc_normal*f_w*v == f_w*loc_flux.
164  // Then the phase fluxes become
165  // flux_w = f_w*(loc_flux + loc_area*loc_normal*lambda_o*K*(rho_w - rho_o)*g)
166  // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
167  // := lambda_o*G (only scalar case)
168  // flux_o = f_o*(loc_flux - lambda_w*G)
169  // In the above, we must decide where to evaluate K, and for this purpose (deciding
170  // upstream directions) we use a K averaged between the two cells.
171  // Since all mobilities and fractional flow functions are positive, the sign
172  // of one of these cases is trivial. If G >= 0, flux_w is in the same direction as
173  // loc_flux, if G <= 0, flux_o is in the same direction as loc_flux.
174  // The phase k for which flux_k and loc_flux are of same sign, is called the trivial
175  // phase in the code below.
176  //
177  // Assuming for the moment that G >=0, we know the direction of the water flux
178  // (same as loc_flux) and evaluate lambda_w in the upstream cell. Then we may use
179  // that lambda_w to evaluate flux_o using the above formula. Knowing flux_o, we know
180  // the direction of the oil flux, and can evaluate lambda_o in the corresponding
181  // upstream cell. Finally, we can use the equation for flux_w to compute that flux.
182  // The opposite case is similar.
183  //
184  // What about tensorial mobilities? In the following code, we make the assumption
185  // that the directions of vectors are not so changed by the multiplication with
186  // mobility tensors that upstream directions change. In other words, we let all
187  // the upstream logic stand as it is. This assumption may need to be revisited.
188  // A worse problem is that
189  // 1) we do not have v, just loc_area*loc_normal*v,
190  // 2) we cannot define G, since the lambdas do not commute with the dot product.
191 
192  typedef typename UpstreamSolver::RP::Mobility Mob;
193  using Opm::utils::arithmeticAverage;
194  // Doing arithmetic averages. Should we consider harmonic or geometric instead?
195  const MutablePermTensor aver_perm
196  = arithAver(s.preservoir_properties_->permeability(cell[0]),
197  s.preservoir_properties_->permeability(cell[1]));
198  // Computing the raw gravity influence vector = (rho_w - rho_o)Kg
199  Vector grav_influence = prod(aver_perm, gravity);
200  grav_influence *= delta_rho;
201  // Computing G. Note that we do not multiply with the mobility,
202  // so this G is wrong in case of anisotropic relperm.
203  const double G = s.method_gravity_ ?
204  loc_area*inner(loc_normal, grav_influence)
205  : 0.0;
206  const int triv_phase = G >= 0.0 ? 0 : 1;
207  const int ups_cell = loc_flux >= 0.0 ? 0 : 1;
208  // Compute mobility of the trivial phase.
209  Mob m_ups[2];
210  s.preservoir_properties_->phaseMobility(triv_phase, cell[ups_cell],
211  cell_sat[ups_cell], m_ups[triv_phase].mob);
212  // Compute gravity flow of the nontrivial phase.
213  double sign_G[2] = { -1.0, 1.0 };
214  double grav_flux_nontriv = sign_G[triv_phase]*loc_area
215  *inner(loc_normal, m_ups[triv_phase].multiply(grav_influence));
216  // Find flow direction of nontrivial phase.
217  const int ups_cell_nontriv = (loc_flux + grav_flux_nontriv >= 0.0) ? 0 : 1;
218  const int nontriv_phase = (triv_phase + 1) % 2;
219  s.preservoir_properties_->phaseMobility(nontriv_phase, cell[ups_cell_nontriv],
220  cell_sat[ups_cell_nontriv], m_ups[nontriv_phase].mob);
221  // Now we have the upstream phase mobilities in m_ups[].
222  Mob m_tot;
223  m_tot.setToSum(m_ups[0], m_ups[1]);
224  Mob m_totinv;
225  m_totinv.setToInverse(m_tot);
226 
227 
228  const double aver_sat
229  = Opm::utils::arithmeticAverage<double, double>(cell_sat[0], cell_sat[1]);
230 
231  Mob m1c0, m1c1, m2c0, m2c1;
232  s.preservoir_properties_->phaseMobility(0, cell[0], aver_sat, m1c0.mob);
233  s.preservoir_properties_->phaseMobility(0, cell[1], aver_sat, m1c1.mob);
234  s.preservoir_properties_->phaseMobility(1, cell[0], aver_sat, m2c0.mob);
235  s.preservoir_properties_->phaseMobility(1, cell[1], aver_sat, m2c1.mob);
236  Mob m_aver[2];
237  m_aver[0].setToAverage(m1c0, m1c1);
238  m_aver[1].setToAverage(m2c0, m2c1);
239  Mob m_aver_tot;
240  m_aver_tot.setToSum(m_aver[0], m_aver[1]);
241  Mob m_aver_totinv;
242  m_aver_totinv.setToInverse(m_aver_tot);
243 
244  // Viscous (pressure driven) term.
245  if (s.method_viscous_) {
246  // v is not correct for anisotropic relperm.
247  Vector v(loc_normal);
248  v *= loc_flux;
249  const double visc_change = inner(loc_normal, m_ups[0].multiply(m_totinv.multiply(v)));
250  // const double visc_change = (m_ups[0].mob/(m_ups[1].mob + m_ups[0].mob))*loc_flux;
251  // std::cout << "New: " << visc_change_2 << " old: " << visc_change << '\n';
252  dS += visc_change;
253  }
254 
255  // Gravity term.
256  if (s.method_gravity_) {
257  if (cell[0] != cell[1]) {
258  // We only add gravity flux on internal or periodic faces.
259  const double grav_change = loc_area
260  *inner(loc_normal, m_ups[0].multiply(m_totinv.multiply(m_ups[1].multiply(grav_influence))));
261  // const double grav_change = (lambda_one*lambda_two/(lambda_two+lambda_one))*G;
262  // const double grav_change = (lambda_one*lambda_two/(lambda_two+lambda_one))*loc_gravity_flux;
263  dS += grav_change;
264  }
265  }
266 
267  // Capillary term.
268  if (s.method_capillary_) {
269  // J(s_w) = \frac{p_c(s_w)\sqrt{k/\phi}}{\sigma \cos\theta}
270  // p_c = \frac{J \sigma \cos\theta}{\sqrt{k/\phi}}
271  Vector cap_influence = prod(aver_perm, s.estimateCapPressureGradient(f, nbface));
272  const double cap_change = loc_area
273  *inner(loc_normal, m_aver[0].multiply(m_aver_totinv.multiply(m_aver[1].multiply(cap_influence))));
274  dS += cap_change;
275  }
276 
277  // Modify saturation.
278  if (cell[0] != cell[1]){
279  residual[cell[0]] -= dS;
280  residual[cell[1]] += dS;
281  } else {
282  assert(cell[0] == cell[1]);
283  residual[cell[0]] -= dS;
284  }
285  }
286  // Source term.
287  double rate = s.pinjection_rates_->element(cell[0]);
288  if (rate < 0.0) {
289  // For anisotropic relperm, fractionalFlow does not really make sense
290  // as a scalar
291  rate *= s.preservoir_properties_->fractionalFlow(cell[0], cell_sat[0]);
292  }
293  residual[cell[0]] += rate;
294  }
295  };
296 
297  template <typename Iter>
299  {
300  typedef Iter Iterator;
301  IndirectRange(const std::vector<Iter>& iters)
302  : iters_(iters), beg_(0), end_(iters_.size() - 1)
303  {
304  assert(iters_.size() >= 2);
305  }
306 #ifdef USE_TBB
307  IndirectRange(IndirectRange& r, tbb::split)
308  : iters_(r.iters_)
309  {
310  int m = (r.beg_ + r.end_)/2;
311  beg_ = m;
312  end_ = r.end_;
313  r.end_ = m;
314  }
315 #endif
316  bool empty() const
317  {
318  return beg_ == end_;
319  }
320  bool is_divisible() const
321  {
322  return end_ - beg_ > 1;
323  }
324  Iter begin() const
325  {
326  return iters_[beg_];
327  }
328  Iter end() const
329  {
330  return iters_[end_];
331  }
332  private:
333  const std::vector<Iter>& iters_;
334  int beg_;
335  int end_;
336  };
337 
338  template <class Updater>
340  {
341  UpdateLoopBody(const Updater& upd)
342  : updater(upd)
343  {
344  }
345  const Updater& updater;
346  template <class Range>
347  void operator()(const Range& r) const
348  {
349  typename Range::Iterator c = r.begin();
350  typename Range::Iterator cend = r.end();
351  for (; c != cend; ++c) {
352  updater(c);
353  }
354  }
355  };
356 
357  } // namespace EulerUpstreamResidualDetails
358 
359 
360  // --------- Member functions -----------
361 
362 
363 
364  template <class GI, class RP, class BC>
366  : pgrid_(0),
367  preservoir_properties_(0),
368  pboundary_(0)
369  {
370  }
371 
372 
373  template <class GI, class RP, class BC>
374  inline EulerUpstreamResidual<GI, RP, BC>::EulerUpstreamResidual(const GI& g, const RP& r, const BC& b)
375  : pgrid_(&g),
376  preservoir_properties_(&r),
377  pboundary_(&b)
378  {
379  initFinal();
380  }
381 
382 
383 
384  template <class GI, class RP, class BC>
385  inline void EulerUpstreamResidual<GI, RP, BC>::initObj(const GI& g, const RP& r, const BC& b)
386  {
387  pgrid_ = &g;
388  preservoir_properties_ = &r;
389  pboundary_ = &b;
390  initFinal();
391  }
392 
393 
394 
395 
396  template <class GI, class RP, class BC>
397  inline void EulerUpstreamResidual<GI, RP, BC>::initFinal()
398  {
399  // Build bid_to_face_ mapping for handling periodic conditions.
400  int maxbid = 0;
401  for (typename GI::CellIterator c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
402  for (typename GI::CellIterator::FaceIterator f = c->facebegin(); f != c->faceend(); ++f) {
403  int bid = f->boundaryId();
404  maxbid = std::max(maxbid, bid);
405  }
406  }
407  bid_to_face_.clear();
408  bid_to_face_.resize(maxbid + 1);
409  for (typename GI::CellIterator c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
410  for (typename GI::CellIterator::FaceIterator f = c->facebegin(); f != c->faceend(); ++f) {
411  if (f->boundary() && pboundary_->satCond(*f).isPeriodic()) {
412  bid_to_face_[f->boundaryId()] = f;
413  }
414  }
415  }
416 
417  // Build cell_iters_.
418  const int num_cells_per_iter = std::min(50, pgrid_->numberOfCells());
419  int counter = 0;
420  for (typename GI::CellIterator c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++counter) {
421  if (counter % num_cells_per_iter == 0) {
422  cell_iters_.push_back(c);
423  }
424  }
425  cell_iters_.push_back(pgrid_->cellend());
426  }
427 
428 
429 
430  template <class GI, class RP, class BC>
431  inline const GI& EulerUpstreamResidual<GI, RP, BC>::grid() const
432  {
433  return *pgrid_;
434  }
435 
436 
437 
438  template <class GI, class RP, class BC>
439  inline const RP& EulerUpstreamResidual<GI, RP, BC>::reservoirProperties() const
440  {
441  return *preservoir_properties_;
442  }
443 
444 
445  template <class GI, class RP, class BC>
446  inline const BC& EulerUpstreamResidual<GI, RP, BC>::boundaryConditions() const
447  {
448  return *pboundary_;
449  }
450 
451 
452 
453  template <class GI, class RP, class BC>
454  inline void EulerUpstreamResidual<GI, RP, BC>::computeCapPressures(const std::vector<double>& saturation) const
455  {
456  int num_cells = saturation.size();
457  cap_pressures_.resize(num_cells);
458  for (int cell = 0; cell < num_cells; ++cell) {
459  cap_pressures_[cell] = preservoir_properties_->capillaryPressure(cell, saturation[cell]);
460  }
461  }
462 
463 
464 
465 
466  template <class GI, class RP, class BC>
467  template <class PressureSolution>
468  inline void EulerUpstreamResidual<GI, RP, BC>::
469  computeResidual(const std::vector<double>& saturation,
470  const typename GI::Vector& gravity,
471  const PressureSolution& pressure_sol,
472  const Opm::SparseVector<double>& injection_rates,
473  const bool method_viscous,
474  const bool method_gravity,
475  const bool method_capillary,
476  std::vector<double>& residual) const
477  {
478  // Make sure sat_change is zero, and has the right size.
479  residual.clear();
480  residual.resize(saturation.size(), 0.0);
481 
482  pinjection_rates_ = &injection_rates;
483  method_viscous_ = method_viscous;
484  method_gravity_ = method_gravity;
485  method_capillary_ = method_capillary;
486 
487  // For every face, we will modify residual for adjacent cells.
488  // We loop over every cell and intersection, and modify only if
489  // this cell has lower index than the neighbour, or we are on the boundary.
490  typedef EulerUpstreamResidualDetails::UpdateForCell<EulerUpstreamResidual<GI,RP,BC>, PressureSolution> CellUpdater;
491  CellUpdater update_cell(*this, saturation, gravity, pressure_sol, residual);
492  EulerUpstreamResidualDetails::UpdateLoopBody<CellUpdater> body(update_cell);
493  EulerUpstreamResidualDetails::IndirectRange<CIt> r(cell_iters_);
494 #ifdef USE_TBB
495  tbb::parallel_for(r, body);
496 #else
497  body(r);
498 #endif
499  }
500 
501 
502 
503 
504  template <class GI, class RP, class BC>
505  inline typename GI::Vector
506  EulerUpstreamResidual<GI, RP, BC>::
507  estimateCapPressureGradient(const FIt& f, const FIt& nbf) const
508  {
509  typedef typename GI::CellIterator::FaceIterator Face;
510  typedef typename Face::Cell Cell;
511  typedef typename GI::Vector Vector;
512 
513  // At nonperiodic boundaries, we return a zero gradient.
514  // That is (sort of) a trivial Neumann (noflow) condition for the capillary pressure.
515  if (f->boundary() && !pboundary_->satCond(*f).isPeriodic()) {
516  return Vector(0.0);
517  }
518  // Find neighbouring cell and face: nbc and nbf.
519  // If we are not on a periodic boundary, nbf is of course equal to f.
520  Cell c = f->cell();
521  Cell nb = f->boundary() ? (f == nbf ? c : nbf->cell()) : f->neighbourCell();
522 
523  // Estimate the gradient like a finite difference between
524  // cell centers, except that in order to handle periodic
525  // conditions we pass through the face centroid(s).
526  Vector cell_c = c.centroid();
527  Vector nb_c = nb.centroid();
528  Vector f_c = f->centroid();
529  Vector nbf_c = nbf->centroid();
530  double d0 = (cell_c - f_c).two_norm();
531  double d1 = (nb_c - nbf_c).two_norm();
532  int cell = c.index();
533  int nbcell = nb.index();
534  double cp0 = cap_pressures_[cell];
535  double cp1 = cap_pressures_[nbcell];
536  double val = (cp1 - cp0)/(d0 + d1);
537  Vector res = nb_c - nbf_c + f_c - cell_c;
538  res /= res.two_norm();
539  res *= val;
540  return res;
541  }
542 
543 
544 } // namespace Opm
545 
546 
547 #endif // OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
Definition: EulerUpstreamResidual_impl.hpp:339
Definition: EulerUpstreamResidual_impl.hpp:298
EulerUpstreamResidual()
Definition: EulerUpstreamResidual_impl.hpp:365
Dune::FieldVector< typename Matrix::value_type, rows > prod(const Matrix &A, const Dune::FieldVector< typename Matrix::value_type, rows > &x)
Matrix applied to a vector.
Definition: Matrix.hpp:669