36#ifndef OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
37#define OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
47#include <opm/core/utility/Average.hpp>
54#include <tbb/parallel_for.h>
66 namespace EulerUpstreamResidualDetails
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)
73 return Opm::utils::arithmeticAverage<FullMatrix<T, StoragePolicy, OrderingPolicy>,
74 FullMatrix<T, OwnData, OrderingPolicy> >(m1, m2);
77 template <
class UpstreamSolver,
class PressureSolution>
80 typedef typename UpstreamSolver::Vector
Vector;
81 typedef typename UpstreamSolver::FIt
FIt;
82 typedef typename UpstreamSolver::RP::PermTensor
PermTensor;
85 const UpstreamSolver&
s;
92 const std::vector<double>& sat,
94 const PressureSolution& psol,
95 std::vector<double>& res)
104 const double delta_rho =
s.preservoir_properties_->densityDifference();
107 cell[0] = c->index();
111 for (
FIt f = c->facebegin(); f != c->faceend(); ++f) {
117 if (
s.pboundary_->satCond(*f).isPeriodic()) {
118 nbface =
s.bid_to_face_[
s.pboundary_->getPeriodicPartner(f->boundaryId())];
120 cell[1] = nbface->cellIndex();
121 assert(cell[0] != cell[1]);
125 if (cell[0] > cell[1]) {
131 assert(
s.pboundary_->satCond(*f).isDirichlet());
133 cell_sat[1] =
s.pboundary_->satCond(*f).saturation();
136 cell[1] = f->neighbourCellIndex();
137 assert(cell[0] != cell[1]);
138 if (cell[0] > cell[1]) {
146 const double loc_area = f->area();
148 const Vector loc_normal = f->normal();
192 typedef typename UpstreamSolver::RP::Mobility Mob;
193 using Opm::utils::arithmeticAverage;
196 =
arithAver(
s.preservoir_properties_->permeability(cell[0]),
197 s.preservoir_properties_->permeability(cell[1]));
200 grav_influence *= delta_rho;
203 const double G =
s.method_gravity_ ?
204 loc_area*inner(loc_normal, grav_influence)
206 const int triv_phase = G >= 0.0 ? 0 : 1;
207 const int ups_cell = loc_flux >= 0.0 ? 0 : 1;
210 s.preservoir_properties_->phaseMobility(triv_phase, cell[ups_cell],
211 cell_sat[ups_cell], m_ups[triv_phase].mob);
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));
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);
223 m_tot.setToSum(m_ups[0], m_ups[1]);
225 m_totinv.setToInverse(m_tot);
228 const double aver_sat
229 = Opm::utils::arithmeticAverage<double, double>(cell_sat[0], cell_sat[1]);
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);
237 m_aver[0].setToAverage(m1c0, m1c1);
238 m_aver[1].setToAverage(m2c0, m2c1);
240 m_aver_tot.setToSum(m_aver[0], m_aver[1]);
242 m_aver_totinv.setToInverse(m_aver_tot);
245 if (
s.method_viscous_) {
249 const double visc_change = inner(loc_normal, m_ups[0].multiply(m_totinv.multiply(v)));
256 if (
s.method_gravity_) {
257 if (cell[0] != cell[1]) {
259 const double grav_change = loc_area
260 *inner(loc_normal, m_ups[0].multiply(m_totinv.multiply(m_ups[1].multiply(grav_influence))));
268 if (
s.method_capillary_) {
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))));
284 if (cell[0] != cell[1]){
288 assert(cell[0] == cell[1]);
293 double rate =
s.pinjection_rates_->element(cell[0]);
297 rate *=
s.preservoir_properties_->fractionalFlow(cell[0], cell_sat[0]);
303 template <
typename Iter>
308 : iters_(iters), beg_(0), end_(iters_.size() - 1)
310 assert(iters_.size() >= 2);
316 int m = (r.beg_ + r.end_)/2;
328 return end_ - beg_ > 1;
339 const std::vector<Iter>& iters_;
344 template <
class Updater>
352 template <
class Range>
355 typename Range::Iterator c = r.begin();
356 typename Range::Iterator cend = r.end();
357 for (; c != cend; ++c) {
370 template <
class GI,
class RP,
class BC>
373 preservoir_properties_(0),
379 template <
class GI,
class RP,
class BC>
382 preservoir_properties_(&r),
390 template <
class GI,
class RP,
class BC>
394 preservoir_properties_ = &r;
402 template <
class GI,
class RP,
class BC>
407 for (
typename GI::CellIterator c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
408 for (
typename GI::CellIterator::FaceIterator f = c->facebegin(); f != c->faceend(); ++f) {
409 int bid = f->boundaryId();
410 maxbid = std::max(maxbid, bid);
413 bid_to_face_.clear();
414 bid_to_face_.resize(maxbid + 1);
415 for (
typename GI::CellIterator c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
416 for (
typename GI::CellIterator::FaceIterator f = c->facebegin(); f != c->faceend(); ++f) {
417 if (f->boundary() && pboundary_->satCond(*f).isPeriodic()) {
418 bid_to_face_[f->boundaryId()] = f;
424 const int num_cells_per_iter = std::min(50, pgrid_->numberOfCells());
426 for (
typename GI::CellIterator c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c, ++counter) {
427 if (counter % num_cells_per_iter == 0) {
428 cell_iters_.push_back(c);
431 cell_iters_.push_back(pgrid_->cellend());
436 template <
class GI,
class RP,
class BC>
444 template <
class GI,
class RP,
class BC>
447 return *preservoir_properties_;
451 template <
class GI,
class RP,
class BC>
459 template <
class GI,
class RP,
class BC>
462 int num_cells = saturation.size();
463 cap_pressures_.resize(num_cells);
464 for (
int cell = 0; cell < num_cells; ++cell) {
465 cap_pressures_[cell] = preservoir_properties_->capillaryPressure(cell, saturation[cell]);
472 template <
class GI,
class RP,
class BC>
473 template <
class PressureSolution>
476 const typename GI::Vector& gravity,
477 const PressureSolution& pressure_sol,
478 const Opm::SparseVector<double>& injection_rates,
479 const bool method_viscous,
480 const bool method_gravity,
481 const bool method_capillary,
482 std::vector<double>& residual)
const
486 residual.resize(saturation.size(), 0.0);
488 pinjection_rates_ = &injection_rates;
489 method_viscous_ = method_viscous;
490 method_gravity_ = method_gravity;
491 method_capillary_ = method_capillary;
497 CellUpdater update_cell(*
this, saturation, gravity, pressure_sol, residual);
501 tbb::parallel_for(r, body);
510 template <
class GI,
class RP,
class BC>
511 inline typename GI::Vector
515 typedef typename GI::CellIterator::FaceIterator Face;
516 typedef typename Face::Cell Cell;
517 typedef typename GI::Vector Vector;
521 if (f->boundary() && !pboundary_->satCond(*f).isPeriodic()) {
527 Cell nb = f->boundary() ? (f == nbf ? c : nbf->cell()) : f->neighbourCell();
532 Vector cell_c = c.centroid();
533 Vector nb_c = nb.centroid();
534 Vector f_c = f->centroid();
535 Vector nbf_c = nbf->centroid();
536 double d0 = (cell_c - f_c).two_norm();
537 double d1 = (nb_c - nbf_c).two_norm();
538 int cell = c.index();
539 int nbcell = nb.index();
540 double cp0 = cap_pressures_[cell];
541 double cp1 = cap_pressures_[nbcell];
542 double val = (cp1 - cp0)/(d0 + d1);
543 Vector res = nb_c - nbf_c + f_c - cell_c;
544 res /= res.two_norm();
Definition: EulerUpstreamResidual.hpp:60
EulerUpstreamResidual()
Definition: EulerUpstreamResidual_impl.hpp:371
void initObj(const GridInterface &grid, const ReservoirProperties &resprop, const BoundaryConditions &boundary)
Definition: EulerUpstreamResidual_impl.hpp:391
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:445
const GridInterface & grid() const
Definition: EulerUpstreamResidual_impl.hpp:437
const BoundaryConditions & boundaryConditions() const
Definition: EulerUpstreamResidual_impl.hpp:452
void computeCapPressures(const std::vector< double > &saturation) const
Definition: EulerUpstreamResidual_impl.hpp:460
ReservoirProperties RP
Definition: EulerUpstreamResidual.hpp:67
FullMatrix< T, OwnData, OrderingPolicy > arithAver(const FullMatrix< T, StoragePolicy, OrderingPolicy > &m1, const FullMatrix< T, StoragePolicy, OrderingPolicy > &m2)
Definition: EulerUpstreamResidual_impl.hpp:70
Definition: BlackoilFluid.hpp:32
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
Definition: EulerUpstreamResidual_impl.hpp:305
bool is_divisible() const
Definition: EulerUpstreamResidual_impl.hpp:326
IndirectRange(const std::vector< Iter > &iters)
Definition: EulerUpstreamResidual_impl.hpp:307
bool empty() const
Definition: EulerUpstreamResidual_impl.hpp:322
Iter end() const
Definition: EulerUpstreamResidual_impl.hpp:334
Iter Iterator
Definition: EulerUpstreamResidual_impl.hpp:306
Iter begin() const
Definition: EulerUpstreamResidual_impl.hpp:330
Definition: EulerUpstreamResidual_impl.hpp:79
UpstreamSolver::FIt FIt
Definition: EulerUpstreamResidual_impl.hpp:81
UpstreamSolver::RP::MutablePermTensor MutablePermTensor
Definition: EulerUpstreamResidual_impl.hpp:83
const UpstreamSolver & s
Definition: EulerUpstreamResidual_impl.hpp:85
UpdateForCell(const UpstreamSolver &solver, const std::vector< double > &sat, const Vector &grav, const PressureSolution &psol, std::vector< double > &res)
Definition: EulerUpstreamResidual_impl.hpp:91
const std::vector< double > & saturation
Definition: EulerUpstreamResidual_impl.hpp:86
const PressureSolution & pressure_sol
Definition: EulerUpstreamResidual_impl.hpp:88
std::vector< double > & residual
Definition: EulerUpstreamResidual_impl.hpp:89
void operator()(const CIt &c) const
Definition: EulerUpstreamResidual_impl.hpp:101
const Vector & gravity
Definition: EulerUpstreamResidual_impl.hpp:87
UpstreamSolver::RP::PermTensor PermTensor
Definition: EulerUpstreamResidual_impl.hpp:82
UpstreamSolver::Vector Vector
Definition: EulerUpstreamResidual_impl.hpp:80
Definition: EulerUpstreamResidual_impl.hpp:346
void operator()(const Range &r) const
Definition: EulerUpstreamResidual_impl.hpp:353
UpdateLoopBody(const Updater &upd)
Definition: EulerUpstreamResidual_impl.hpp:347
const Updater & updater
Definition: EulerUpstreamResidual_impl.hpp:351