Loading...
Searching...
No Matches
ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp
//===========================================================================
//
// File: ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp
//
// Created: Mon Oct 26 10:12:15 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B�rd Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
#define OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
namespace Opm
{
template <int dim>
template <class MatrixType>
void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::phaseMobility(int phase_index,
int cell_index,
double saturation,
MatrixType& phase_mob) const
{
const int region = Super::rock_.size() > 0 ? Super::cell_to_rock_[cell_index] : -1;
phaseMobilityByRock(phase_index, region, saturation, phase_mob);
}
template <int dim>
double ReservoirPropertyCapillaryAnisotropicRelperm<dim>::fractionalFlow(int cell_index, double saturation) const
{
// This method is a hack.
// Assumes that the relperm is diagonal.
Mobility m;
double ff_first = 0.0;
for (int direction = 0; direction < dim; ++direction) {
phaseMobility(0, cell_index, saturation, m.mob);
double l1 = m.mob(direction, direction);
phaseMobility(1, cell_index, saturation, m.mob);
double l2 = m.mob(direction, direction);
ff_first += l1/(l1 + l2);
}
ff_first /= double(dim);
return ff_first;
}
template <int dim>
template <class MatrixType>
void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::phaseMobilityByRock(int phase_index,
int rock_index,
double saturation,
MatrixType& phase_mob) const
{
assert ((0 <= phase_index) && (phase_index < Super::NumberOfPhases));
static_assert(Super::NumberOfPhases == 2, "");
double visc = phase_index == 0 ? Super::viscosity1_ : Super::viscosity2_;
if (rock_index != -1) {
assert (rock_index < int(Super::rock_.size()));
Super::rock_[rock_index].kr(phase_index, saturation, phase_mob);
std::transform(phase_mob.data(), phase_mob.data() + dim*dim,
phase_mob.data(), [visc](double param) { return param / visc; });
} else {
// HACK: With no rocks, we use quadratic relperm functions.
double kr = phase_index == 0 ? saturation*saturation : (1.0 - saturation)*(1.0 - saturation);
double mob = kr/visc;
zero(phase_mob);
for (int j = 0; j < dim; ++j) {
phase_mob(j,j) = mob;
}
}
}
template <int dim>
void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::cflFracFlows(int rock, int direction,
double s, double& ff_first, double& ff_gravity) const
{
// Assumes that the relperm is diagonal.
Mobility m;
phaseMobilityByRock(0, rock, s, m.mob);
double l1 = m.mob(direction, direction);
phaseMobilityByRock(1, rock, s, m.mob);
double l2 = m.mob(direction, direction);
ff_first = l1/(l1 + l2);
ff_gravity = l1*l2/(l1 + l2);
}
template <int dim>
std::array<double, 3> ReservoirPropertyCapillaryAnisotropicRelperm<dim>::computeSingleRockCflFactors(int rock) const
{
const int N = 257;
double delta = 1.0/double(N - 1);
double last_ff1, last_ffg;
double max_der1 = -1e100;
double max_derg = -1e100;
double min_ffg = -1e100;
for (int direction = 0; direction < dim; ++direction) {
cflFracFlows(rock, direction, 0.0, last_ff1, last_ffg);
min_ffg = std::min(min_ffg, last_ffg);
for (int i = 1; i < N; ++i) {
double s = double(i)*delta;
double ff1, ffg;
cflFracFlows(rock, direction, s, ff1, ffg);
double est_deriv_ff1 = std::fabs(ff1 - last_ff1)/delta;
double est_deriv_ffg = std::fabs(ffg - last_ffg)/delta;
max_der1 = std::max(max_der1, est_deriv_ff1);
max_derg = std::max(max_derg, est_deriv_ffg);
min_ffg = std::min(min_ffg, last_ffg);
last_ff1 = ff1;
last_ffg = ffg;
}
}
std::array<double, 3> retval = {{ 1.0/max_der1, 1.0/max_derg, min_ffg }};
return retval;
}
template <int dim>
void ReservoirPropertyCapillaryAnisotropicRelperm<dim>::computeCflFactors()
{
if (Super::rock_.empty()) {
std::array<double, 3> fac = computeSingleRockCflFactors(-1);
Super::cfl_factor_ = fac[0];
Super::cfl_factor_gravity_ = fac[1];
Super::cfl_factor_capillary_ = fac[2];
} else {
Super::cfl_factor_ = 1e100;
Super::cfl_factor_gravity_ = 1e100;
Super::cfl_factor_capillary_ = 1e100;
for (int r = 0; r < int(Super::rock_.size()); ++r) {
std::array<double, 3> fac = computeSingleRockCflFactors(r);
Super::cfl_factor_ = std::min(Super::cfl_factor_, fac[0]);
Super::cfl_factor_gravity_ = std::min(Super::cfl_factor_gravity_, fac[1]);
Super::cfl_factor_capillary_ = std::max(Super::cfl_factor_capillary_, fac[2]);
}
}
}
} // namespace Opm
#endif // OPENRS_RESERVOIRPROPERTYCAPILLARYANISOTROPICRELPERM_IMPL_HEADER
Opm::ReservoirPropertyCapillaryAnisotropicRelperm::fractionalFlow
double fractionalFlow(int cell_index, double saturation) const
Some approximation to a scalar fractional flow (of the first phase).
Definition ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:55
Opm::ReservoirPropertyCapillaryAnisotropicRelperm::Mobility
TensorMobility< dim > Mobility
The (tensorial) mobility type.
Definition ReservoirPropertyCapillaryAnisotropicRelperm.hpp:109
Opm::ReservoirPropertyCapillaryAnisotropicRelperm::phaseMobility
void phaseMobility(int phase_index, int cell_index, double saturation, MatrixType &phase_mob) const
Anisotropic phase mobility.
Definition ReservoirPropertyCapillaryAnisotropicRelperm_impl.hpp:44
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::zero
void zero(Matrix &A)
Zero-fill a.
Definition Matrix.hpp:602