36#ifndef OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
37#define OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
39#include <boost/lambda/lambda.hpp>
45 template <
class MatrixType>
49 MatrixType& phase_mob)
const
51 const int region = Super::rock_.size() > 0 ? Super::cell_to_rock_[cell_index] : -1;
52 phaseMobilityByRock(phase_index, region, saturation, phase_mob);
62 double ff_first = 0.0;
63 for (
int direction = 0; direction < dim; ++direction) {
64 phaseMobility(0, cell_index, saturation, m.
mob);
65 double l1 = m.
mob(direction, direction);
66 phaseMobility(1, cell_index, saturation, m.
mob);
67 double l2 = m.
mob(direction, direction);
68 ff_first += l1/(l1 + l2);
70 ff_first /= double(dim);
76 template <
class MatrixType>
80 MatrixType& phase_mob)
const
82 assert ((0 <= phase_index) && (phase_index < Super::NumberOfPhases));
83 static_assert(Super::NumberOfPhases == 2,
"");
85 double visc = phase_index == 0 ? Super::viscosity1_ : Super::viscosity2_;
86 if (rock_index != -1) {
87 assert (rock_index <
int(Super::rock_.size()));
88 Super::rock_[rock_index].kr(phase_index, saturation, phase_mob);
90 std::transform(phase_mob.data(), phase_mob.data() + dim*dim,
91 phase_mob.data(), boost::lambda::_1/visc);
94 double kr = phase_index == 0 ? saturation*saturation : (1.0 - saturation)*(1.0 - saturation);
97 for (
int j = 0; j < dim; ++j) {
105 void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::cflFracFlows(
int rock,
int direction,
106 double s,
double& ff_first,
double& ff_gravity)
const
110 phaseMobilityByRock(0, rock, s, m.mob);
111 double l1 = m.mob(direction, direction);
112 phaseMobilityByRock(1, rock, s, m.mob);
113 double l2 = m.mob(direction, direction);
114 ff_first = l1/(l1 + l2);
115 ff_gravity = l1*l2/(l1 + l2);
122 std::array<double, 3> ReservoirPropertyCapillaryAnisotropicRelperm<dim>::computeSingleRockCflFactors(
int rock)
const
125 double delta = 1.0/double(N - 1);
126 double last_ff1, last_ffg;
127 double max_der1 = -1e100;
128 double max_derg = -1e100;
129 double min_ffg = -1e100;
130 for (
int direction = 0; direction < dim; ++direction) {
131 cflFracFlows(rock, direction, 0.0, last_ff1, last_ffg);
132 min_ffg = std::min(min_ffg, last_ffg);
133 for (
int i = 1; i < N; ++i) {
134 double s = double(i)*delta;
136 cflFracFlows(rock, direction, s, ff1, ffg);
137 double est_deriv_ff1 = std::fabs(ff1 - last_ff1)/delta;
138 double est_deriv_ffg = std::fabs(ffg - last_ffg)/delta;
139 max_der1 = std::max(max_der1, est_deriv_ff1);
140 max_derg = std::max(max_derg, est_deriv_ffg);
141 min_ffg = std::min(min_ffg, last_ffg);
146 std::array<double, 3> retval = {{ 1.0/max_der1, 1.0/max_derg, min_ffg }};
156 if (Super::rock_.empty()) {
157 std::array<double, 3> fac = computeSingleRockCflFactors(-1);
158 Super::cfl_factor_ = fac[0];
159 Super::cfl_factor_gravity_ = fac[1];
160 Super::cfl_factor_capillary_ = fac[2];
162 Super::cfl_factor_ = 1e100;
163 Super::cfl_factor_gravity_ = 1e100;
164 Super::cfl_factor_capillary_ = 1e100;
165 for (
int r = 0; r < int(Super::rock_.size()); ++r) {
166 std::array<double, 3> fac = computeSingleRockCflFactors(r);
167 Super::cfl_factor_ = std::min(Super::cfl_factor_, fac[0]);
168 Super::cfl_factor_gravity_ = std::min(Super::cfl_factor_gravity_, fac[1]);
169 Super::cfl_factor_capillary_ = std::max(Super::cfl_factor_capillary_, fac[2]);
A property class for incompressible two-phase flow.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm.hpp:106
double fractionalFlow(int cell_index, double saturation) const
Some approximation to a scalar fractional flow (of the first phase).
Definition: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:57
void computeCflFactors()
Definition: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:154
void phaseMobility(int phase_index, int cell_index, double saturation, MatrixType &phase_mob) const
Anisotropic phase mobility.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:46
Definition: BlackoilFluid.hpp:32
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
A wrapper for a tensor.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm.hpp:50
FullMatrix< double, SharedData, COrdering > mob
Definition: ReservoirPropertyCapillaryAnisotropicRelperm.hpp:96