36#ifndef OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
37#define OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
41#include <opm/core/utility/Average.hpp>
42#include <opm/porsol/common/Matrix.hpp>
45#include <tbb/parallel_for.h>
57 namespace EulerUpstreamResidualDetails
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)
64 return Opm::utils::arithmeticAverage<FullMatrix<T, StoragePolicy, OrderingPolicy>,
65 FullMatrix<T, OwnData, OrderingPolicy> >(m1, m2);
68 template <
class UpstreamSolver,
class PressureSolution>
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;
76 const UpstreamSolver& s;
77 const std::vector<double>& saturation;
78 const Vector& gravity;
79 const PressureSolution& pressure_sol;
80 std::vector<double>& residual;
83 const std::vector<double>& sat,
85 const PressureSolution& psol,
86 std::vector<double>& res)
87 : s(solver), saturation(sat), gravity(grav), pressure_sol(psol), residual(res)
92 void operator()(
const CIt& c)
const
95 const double delta_rho = s.preservoir_properties_->densityDifference();
99 cell_sat[0] = saturation[cell[0]];
102 for (FIt f = c->facebegin(); f != c->faceend(); ++f) {
108 if (s.pboundary_->satCond(*f).isPeriodic()) {
109 nbface = s.bid_to_face_[s.pboundary_->getPeriodicPartner(f->boundaryId())];
111 cell[1] = nbface->cellIndex();
112 assert(cell[0] != cell[1]);
116 if (cell[0] > cell[1]) {
120 cell_sat[1] = saturation[cell[1]];
122 assert(s.pboundary_->satCond(*f).isDirichlet());
124 cell_sat[1] = s.pboundary_->satCond(*f).saturation();
127 cell[1] = f->neighbourCellIndex();
128 assert(cell[0] != cell[1]);
129 if (cell[0] > cell[1]) {
133 cell_sat[1] = saturation[cell[1]];
137 const double loc_area = f->area();
138 const double loc_flux = pressure_sol.outflux(f);
139 const Vector loc_normal = f->normal();
183 typedef typename UpstreamSolver::RP::Mobility Mob;
184 using Opm::utils::arithmeticAverage;
186 const MutablePermTensor aver_perm
187 = arithAver(s.preservoir_properties_->permeability(cell[0]),
188 s.preservoir_properties_->permeability(cell[1]));
190 Vector grav_influence =
prod(aver_perm, gravity);
191 grav_influence *= delta_rho;
194 const double G = s.method_gravity_ ?
195 loc_area*inner(loc_normal, grav_influence)
197 const int triv_phase = G >= 0.0 ? 0 : 1;
198 const int ups_cell = loc_flux >= 0.0 ? 0 : 1;
201 s.preservoir_properties_->phaseMobility(triv_phase, cell[ups_cell],
202 cell_sat[ups_cell], m_ups[triv_phase].mob);
204 const 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));
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);
214 m_tot.setToSum(m_ups[0], m_ups[1]);
216 m_totinv.setToInverse(m_tot);
219 const double aver_sat
220 = Opm::utils::arithmeticAverage<double, double>(cell_sat[0], cell_sat[1]);
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);
228 m_aver[0].setToAverage(m1c0, m1c1);
229 m_aver[1].setToAverage(m2c0, m2c1);
231 m_aver_tot.setToSum(m_aver[0], m_aver[1]);
233 m_aver_totinv.setToInverse(m_aver_tot);
236 if (s.method_viscous_) {
238 Vector v(loc_normal);
240 const double visc_change = inner(loc_normal, m_ups[0].multiply(m_totinv.multiply(v)));
247 if (s.method_gravity_) {
248 if (cell[0] != cell[1]) {
250 const double grav_change = loc_area
251 *inner(loc_normal, m_ups[0].multiply(m_totinv.multiply(m_ups[1].multiply(grav_influence))));
259 if (s.method_capillary_) {
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))));
269 if (cell[0] != cell[1]){
270 residual[cell[0]] -= dS;
271 residual[cell[1]] += dS;
273 assert(cell[0] == cell[1]);
274 residual[cell[0]] -= dS;
278 double rate = s.pinjection_rates_->element(cell[0]);
282 rate *= s.preservoir_properties_->fractionalFlow(cell[0], cell_sat[0]);
284 residual[cell[0]] += rate;
288 template <
typename Iter>
291 typedef Iter Iterator;
293 : iters_(iters), beg_(0), end_(iters_.size() - 1)
295 assert(iters_.size() >= 2);
301 int m = (r.beg_ + r.end_)/2;
311 bool is_divisible()
const
313 return end_ - beg_ > 1;
324 const std::vector<Iter>& iters_;
329 template <
class Updater>
336 const Updater& updater;
337 template <
class Range>
338 void operator()(
const Range& r)
const
340 typename Range::Iterator c = r.begin();
341 typename Range::Iterator cend = r.end();
342 for (; c != cend; ++c) {
355 template <
class GI,
class RP,
class BC>
358 preservoir_properties_(0),
364 template <
class GI,
class RP,
class BC>
367 preservoir_properties_(&r),
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)
379 preservoir_properties_ = &r;
387 template <
class GI,
class RP,
class BC>
388 inline void EulerUpstreamResidual<GI, RP, BC>::initFinal()
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);
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;
409 const int num_cells_per_iter = std::min(50, pgrid_->numberOfCells());
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);
416 cell_iters_.push_back(pgrid_->cellend());
421 template <
class GI,
class RP,
class BC>
422 inline const GI& EulerUpstreamResidual<GI, RP, BC>::grid()
const
429 template <
class GI,
class RP,
class BC>
430 inline const RP& EulerUpstreamResidual<GI, RP, BC>::reservoirProperties()
const
432 return *preservoir_properties_;
436 template <
class GI,
class RP,
class BC>
437 inline const BC& EulerUpstreamResidual<GI, RP, BC>::boundaryConditions()
const
444 template <
class GI,
class RP,
class BC>
445 inline void EulerUpstreamResidual<GI, RP, BC>::computeCapPressures(
const std::vector<double>& saturation)
const
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]);
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
471 residual.resize(saturation.size(), 0.0);
473 pinjection_rates_ = &injection_rates;
474 method_viscous_ = method_viscous;
475 method_gravity_ = method_gravity;
476 method_capillary_ = method_capillary;
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_);
486 tbb::parallel_for(r, body);
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
502 if (f->boundary() && !pboundary_->satCond(*f).isPeriodic()) {
508 auto nb = f->boundary() ? (f == nbf ? c : nbf->cell()) : f->neighbourCell();
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();
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