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