36#ifndef OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
37#define OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
43 template <
class MatrixType>
47 MatrixType& phase_mob)
const
49 const int region = Super::rock_.size() > 0 ? Super::cell_to_rock_[cell_index] : -1;
50 phaseMobilityByRock(phase_index, region, saturation, phase_mob);
60 double ff_first = 0.0;
61 for (
int direction = 0; direction < dim; ++direction) {
62 phaseMobility(0, cell_index, saturation, m.
mob);
63 double l1 = m.
mob(direction, direction);
64 phaseMobility(1, cell_index, saturation, m.
mob);
65 double l2 = m.
mob(direction, direction);
66 ff_first += l1/(l1 + l2);
68 ff_first /= double(dim);
74 template <
class MatrixType>
78 MatrixType& phase_mob)
const
80 assert ((0 <= phase_index) && (phase_index < Super::NumberOfPhases));
81 static_assert(Super::NumberOfPhases == 2,
"");
83 double visc = phase_index == 0 ? Super::viscosity1_ : Super::viscosity2_;
84 if (rock_index != -1) {
85 assert (rock_index <
int(Super::rock_.size()));
86 Super::rock_[rock_index].kr(phase_index, saturation, phase_mob);
87 std::transform(phase_mob.data(), phase_mob.data() + dim*dim,
88 phase_mob.data(), [visc](
double param) { return param / visc; });
91 double kr = phase_index == 0 ? saturation*saturation : (1.0 - saturation)*(1.0 - saturation);
94 for (
int j = 0;
j < dim; ++
j) {
102 void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::cflFracFlows(
int rock,
int direction,
103 double s,
double& ff_first,
double& ff_gravity)
const
107 phaseMobilityByRock(0, rock, s, m.mob);
108 double l1 = m.mob(direction, direction);
109 phaseMobilityByRock(1, rock, s, m.mob);
110 double l2 = m.mob(direction, direction);
111 ff_first = l1/(l1 + l2);
112 ff_gravity = l1*l2/(l1 + l2);
119 std::array<double, 3> ReservoirPropertyCapillaryAnisotropicRelperm<dim>::computeSingleRockCflFactors(
int rock)
const
122 double delta = 1.0/double(N - 1);
123 double last_ff1, last_ffg;
124 double max_der1 = -1e100;
125 double max_derg = -1e100;
126 double min_ffg = -1e100;
127 for (
int direction = 0; direction < dim; ++direction) {
128 cflFracFlows(rock, direction, 0.0, last_ff1, last_ffg);
129 min_ffg =
std::min(min_ffg, last_ffg);
130 for (
int i = 1; i < N; ++i) {
131 double s = double(i)*delta;
133 cflFracFlows(rock, direction, s, ff1, ffg);
134 double est_deriv_ff1 = std::fabs(ff1 - last_ff1)/delta;
135 double est_deriv_ffg = std::fabs(ffg - last_ffg)/delta;
136 max_der1 = std::max(max_der1, est_deriv_ff1);
137 max_derg = std::max(max_derg, est_deriv_ffg);
138 min_ffg =
std::min(min_ffg, last_ffg);
143 std::array<double, 3> retval = {{ 1.0/max_der1, 1.0/max_derg, min_ffg }};
153 if (Super::rock_.empty()) {
154 std::array<double, 3> fac = computeSingleRockCflFactors(-1);
155 Super::cfl_factor_ = fac[0];
156 Super::cfl_factor_gravity_ = fac[1];
157 Super::cfl_factor_capillary_ = fac[2];
159 Super::cfl_factor_ = 1e100;
160 Super::cfl_factor_gravity_ = 1e100;
161 Super::cfl_factor_capillary_ = 1e100;
162 for (
int r = 0; r < int(Super::rock_.size()); ++r) {
163 std::array<double, 3> fac = computeSingleRockCflFactors(r);
164 Super::cfl_factor_ =
std::min(Super::cfl_factor_, fac[0]);
165 Super::cfl_factor_gravity_ =
std::min(Super::cfl_factor_gravity_, fac[1]);
166 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:55
void computeCflFactors()
Definition: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:151
void phaseMobility(int phase_index, int cell_index, double saturation, MatrixType &phase_mob) const
Anisotropic phase mobility.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:44
min[0]
Definition: elasticity_upscale_impl.hpp:146
int j
Definition: elasticity_upscale_impl.hpp:658
Definition: ImplicitAssembly.hpp:43
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:602
A wrapper for a tensor.
Definition: ReservoirPropertyCapillaryAnisotropicRelperm.hpp:50
FullMatrix< double, SharedData, COrdering > mob
Definition: ReservoirPropertyCapillaryAnisotropicRelperm.hpp:96