36 #ifndef OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
37 #define OPENRS_EULERUPSTREAMRESIDUAL_IMPL_HEADER
47 #include <opm/core/utility/Average.hpp>
51 #include <opm/porsol/common/Matrix.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;
83 typedef typename UpstreamSolver::RP::MutablePermTensor MutablePermTensor;
85 const UpstreamSolver& s;
86 const std::vector<double>& saturation;
87 const Vector& gravity;
88 const PressureSolution& pressure_sol;
89 std::vector<double>& residual;
91 UpdateForCell(
const UpstreamSolver& solver,
92 const std::vector<double>& sat,
94 const PressureSolution& psol,
95 std::vector<double>& res)
96 : s(solver), saturation(sat), gravity(grav), pressure_sol(psol), residual(res)
101 void operator()(
const CIt& c)
const
104 const double delta_rho = s.preservoir_properties_->densityDifference();
107 cell[0] = c->index();
108 cell_sat[0] = saturation[cell[0]];
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]) {
129 cell_sat[1] = saturation[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]) {
142 cell_sat[1] = saturation[cell[1]];
146 const double loc_area = f->area();
147 const double loc_flux = pressure_sol.outflux(f);
148 const Vector loc_normal = f->normal();
192 typedef typename UpstreamSolver::RP::Mobility Mob;
193 using Opm::utils::arithmeticAverage;
195 const MutablePermTensor aver_perm
196 = arithAver(s.preservoir_properties_->permeability(cell[0]),
197 s.preservoir_properties_->permeability(cell[1]));
199 Vector grav_influence =
prod(aver_perm, gravity);
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_) {
247 Vector v(loc_normal);
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_) {
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))));
278 if (cell[0] != cell[1]){
279 residual[cell[0]] -= dS;
280 residual[cell[1]] += dS;
282 assert(cell[0] == cell[1]);
283 residual[cell[0]] -= dS;
287 double rate = s.pinjection_rates_->element(cell[0]);
291 rate *= s.preservoir_properties_->fractionalFlow(cell[0], cell_sat[0]);
293 residual[cell[0]] += rate;
297 template <
typename Iter>
300 typedef Iter Iterator;
302 : iters_(iters), beg_(0), end_(iters_.size() - 1)
304 assert(iters_.size() >= 2);
310 int m = (r.beg_ + r.end_)/2;
320 bool is_divisible()
const
322 return end_ - beg_ > 1;
333 const std::vector<Iter>& iters_;
338 template <
class Updater>
345 const Updater& updater;
346 template <
class Range>
347 void operator()(
const Range& r)
const
349 typename Range::Iterator c = r.begin();
350 typename Range::Iterator cend = r.end();
351 for (; c != cend; ++c) {
364 template <
class GI,
class RP,
class BC>
367 preservoir_properties_(0),
373 template <
class GI,
class RP,
class BC>
376 preservoir_properties_(&r),
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)
388 preservoir_properties_ = &r;
396 template <
class GI,
class RP,
class BC>
397 inline void EulerUpstreamResidual<GI, RP, BC>::initFinal()
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);
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;
418 const int num_cells_per_iter = std::min(50, pgrid_->numberOfCells());
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);
425 cell_iters_.push_back(pgrid_->cellend());
430 template <
class GI,
class RP,
class BC>
431 inline const GI& EulerUpstreamResidual<GI, RP, BC>::grid()
const
438 template <
class GI,
class RP,
class BC>
439 inline const RP& EulerUpstreamResidual<GI, RP, BC>::reservoirProperties()
const
441 return *preservoir_properties_;
445 template <
class GI,
class RP,
class BC>
446 inline const BC& EulerUpstreamResidual<GI, RP, BC>::boundaryConditions()
const
453 template <
class GI,
class RP,
class BC>
454 inline void EulerUpstreamResidual<GI, RP, BC>::computeCapPressures(
const std::vector<double>& saturation)
const
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]);
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
480 residual.resize(saturation.size(), 0.0);
482 pinjection_rates_ = &injection_rates;
483 method_viscous_ = method_viscous;
484 method_gravity_ = method_gravity;
485 method_capillary_ = method_capillary;
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_);
495 tbb::parallel_for(r, body);
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
509 typedef typename GI::CellIterator::FaceIterator Face;
510 typedef typename Face::Cell Cell;
511 typedef typename GI::Vector Vector;
515 if (f->boundary() && !pboundary_->satCond(*f).isPeriodic()) {
521 Cell nb = f->boundary() ? (f == nbf ? c : nbf->cell()) : f->neighbourCell();
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();
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