36#ifndef OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
37#define OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
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>
72 typedef typename UpstreamSolver::FIt
FIt;
73 typedef typename UpstreamSolver::RP::PermTensor
PermTensor;
76 const UpstreamSolver&
s;
83 const std::vector<double>& sat,
85 const PressureSolution& psol,
86 std::vector<double>& res)
95 const double delta_rho =
s.preservoir_properties_->densityDifference();
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]) {
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]) {
137 const double loc_area = f->area();
139 const Vector loc_normal = f->normal();
183 typedef typename UpstreamSolver::RP::Mobility Mob;
187 =
arithAver(
s.preservoir_properties_->permeability(cell[0]),
188 s.preservoir_properties_->permeability(cell[1]));
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 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_) {
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]){
273 assert(cell[0] == cell[1]);
278 double rate =
s.pinjection_rates_->element(cell[0]);
282 rate *=
s.preservoir_properties_->fractionalFlow(cell[0], cell_sat[0]);
288 template <
typename Iter>
293 : iters_(iters), beg_(0), end_(iters_.size() - 1)
295 assert(iters_.size() >= 2);
301 int m = (r.beg_ + r.end_)/2;
313 return end_ - beg_ > 1;
324 const std::vector<Iter>& iters_;
329 template <
class Updater>
337 template <
class Range>
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>
379 preservoir_properties_ = &r;
387 template <
class GI,
class RP,
class BC>
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>
429 template <
class GI,
class RP,
class BC>
432 return *preservoir_properties_;
436 template <
class GI,
class RP,
class BC>
444 template <
class GI,
class RP,
class BC>
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>
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;
482 CellUpdater update_cell(*
this, saturation, gravity, pressure_sol, residual);
486 tbb::parallel_for(r, body);
495 template <
class GI,
class RP,
class BC>
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();
Definition: EulerUpstreamResidual.hpp:58
EulerUpstreamResidual()
Definition: EulerUpstreamResidual_impl.hpp:356
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary)
Definition: EulerUpstreamResidual_impl.hpp:376
void computeResidual(const std::vector< double > &saturation, const typename GridInterface::Vector &gravity, const FlowSolution &flow_sol, const Opm::SparseVector< double > &injection_rates, const bool method_viscous, const bool method_gravity, const bool method_capillary, std::vector< double > &sat_delta) const
const ReservoirProperties & reservoirProperties() const
Definition: EulerUpstreamResidual_impl.hpp:430
const GridInterface & grid() const
Definition: EulerUpstreamResidual_impl.hpp:422
const BoundaryConditions & boundaryConditions() const
Definition: EulerUpstreamResidual_impl.hpp:437
void computeCapPressures(const std::vector< double > &saturation) const
Definition: EulerUpstreamResidual_impl.hpp:445
ReservoirProperties RP
Definition: EulerUpstreamResidual.hpp:65
min[0]
Definition: elasticity_upscale_impl.hpp:146
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
FullMatrix< T, OwnData, OrderingPolicy > arithAver(const FullMatrix< T, StoragePolicy, OrderingPolicy > &m1, const FullMatrix< T, StoragePolicy, OrderingPolicy > &m2)
Definition: EulerUpstreamResidual_impl.hpp:61
Tresult arithmeticAverage(const T &t1, const T &t2)
Definition: Average.hpp:51
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
bool is_divisible() const
Definition: EulerUpstreamResidual_impl.hpp:311
IndirectRange(const std::vector< Iter > &iters)
Definition: EulerUpstreamResidual_impl.hpp:292
bool empty() const
Definition: EulerUpstreamResidual_impl.hpp:307
Iter end() const
Definition: EulerUpstreamResidual_impl.hpp:319
Iter Iterator
Definition: EulerUpstreamResidual_impl.hpp:291
Iter begin() const
Definition: EulerUpstreamResidual_impl.hpp:315
Definition: EulerUpstreamResidual_impl.hpp:70
UpstreamSolver::FIt FIt
Definition: EulerUpstreamResidual_impl.hpp:72
UpstreamSolver::RP::MutablePermTensor MutablePermTensor
Definition: EulerUpstreamResidual_impl.hpp:74
const UpstreamSolver & s
Definition: EulerUpstreamResidual_impl.hpp:76
UpdateForCell(const UpstreamSolver &solver, const std::vector< double > &sat, const Vector &grav, const PressureSolution &psol, std::vector< double > &res)
Definition: EulerUpstreamResidual_impl.hpp:82
const std::vector< double > & saturation
Definition: EulerUpstreamResidual_impl.hpp:77
const PressureSolution & pressure_sol
Definition: EulerUpstreamResidual_impl.hpp:79
std::vector< double > & residual
Definition: EulerUpstreamResidual_impl.hpp:80
void operator()(const CIt &c) const
Definition: EulerUpstreamResidual_impl.hpp:92
const Vector & gravity
Definition: EulerUpstreamResidual_impl.hpp:78
UpstreamSolver::RP::PermTensor PermTensor
Definition: EulerUpstreamResidual_impl.hpp:73
UpstreamSolver::Vector Vector
Definition: EulerUpstreamResidual_impl.hpp:71
Definition: EulerUpstreamResidual_impl.hpp:331
void operator()(const Range &r) const
Definition: EulerUpstreamResidual_impl.hpp:338
UpdateLoopBody(const Updater &upd)
Definition: EulerUpstreamResidual_impl.hpp:332
const Updater & updater
Definition: EulerUpstreamResidual_impl.hpp:336