36#ifndef OPENRS_CFLCALCULATOR_HEADER
37#define OPENRS_CFLCALCULATOR_HEADER
39#include <opm/common/ErrorMacros.hpp>
40#include <opm/core/utility/Average.hpp>
47 namespace cfl_calculator {
54 template <
class Gr
id,
class ReservoirProperties,
class PressureSolution>
56 const ReservoirProperties& resprop,
57 const PressureSolution& pressure_sol)
60 typename Grid::CellIterator c = grid.cellbegin();
61 for (; c != grid.cellend(); ++c) {
64 typename Grid::CellIterator::FaceIterator f = c->facebegin();
65 for (; f != c->faceend(); ++f) {
66 const double loc_flux = pressure_sol.outflux(f);
73 double flux = std::max(flux_n, flux_p);
74 double loc_dt = (resprop.cflFactor()*c->volume()*resprop.porosity(c->index()))/flux;
76 OPM_THROW(std::runtime_error,
"Cfl computation gave dt = 0.0");
90 template <
class Gr
id,
class ReservoirProperties>
92 const ReservoirProperties& resprop,
93 const typename Grid::Vector& gravity)
95 typedef typename ReservoirProperties::PermTensor PermTensor;
96 typedef typename ReservoirProperties::MutablePermTensor MutablePermTensor;
97 const int dimension = Grid::Vector::dimension;
99 typename Grid::CellIterator c = grid.cellbegin();
100 for (; c != grid.cellend(); ++c) {
102 typename Grid::CellIterator::FaceIterator f = c->facebegin();
103 for (; f != c->faceend(); ++f) {
105 MutablePermTensor loc_perm_aver;
106 const double* permdata = 0;
107 if (!f->boundary()) {
108 PermTensor K0 = resprop.permeability(f->cellIndex());
109 PermTensor K1 = resprop.permeability(f->neighbourCellIndex());
110 loc_perm_aver = Opm::utils::arithmeticAverage<PermTensor, MutablePermTensor>(K0, K1);
111 permdata = loc_perm_aver.data();
113 permdata = resprop.permeability(f->cellIndex()).data();
115 PermTensor loc_perm(dimension, dimension, permdata);
116 typename Grid::Vector loc_halfface_normal = f->normal();
117 double loc_gravity_flux = 0.0;
118 for (
int k = 0; k < dimension; ++k) {
119 for (
int q = 0; q < dimension; ++q) {
120 loc_gravity_flux += loc_halfface_normal[q]*(loc_perm(q,k)*gravity[k]*resprop.densityDifference());
123 loc_gravity_flux *= f->area();
124 if (loc_gravity_flux > 0) {
125 flux += loc_gravity_flux;
128 double loc_dt = (resprop.cflFactorGravity()*c->volume()*resprop.porosity(c->index()))/flux;
142 template <
class Gr
id,
class ReservoirProperties>
144 const ReservoirProperties& resprop)
146 typedef typename ReservoirProperties::PermTensor PermTensor;
147 typedef typename ReservoirProperties::MutablePermTensor MutablePermTensor;
148 const int dimension = Grid::Vector::dimension;
150 typename Grid::CellIterator c = grid.cellbegin();
151 for (; c != grid.cellend(); ++c) {
152 typename Grid::CellIterator::FaceIterator f = c->facebegin();
153 for (; f != c->faceend(); ++f) {
155 MutablePermTensor loc_perm_aver;
156 const double* permdata = 0;
157 if (!f->boundary()) {
158 PermTensor K0 = resprop.permeability(f->cellIndex());
159 PermTensor K1 = resprop.permeability(f->neighbourCellIndex());
160 loc_perm_aver = Opm::utils::arithmeticAverage<PermTensor, MutablePermTensor>(K0, K1);
161 permdata = loc_perm_aver.data();
163 permdata = resprop.permeability(f->cellIndex()).data();
166 MutablePermTensor loc_perm(dimension, dimension, permdata);
167 MutablePermTensor loc_perm_inv =
inverse3x3(loc_perm);
168 typename Grid::Vector loc_centroid = f->centroid();
169 loc_centroid -= c->centroid();
170 double spatial_contrib = loc_centroid*
prod(loc_perm_inv, loc_centroid);
171 double loc_dt = spatial_contrib/resprop.cflFactorCapillary();
172 dt = std::min(dt, loc_dt);
double findCFLtimeVelocity(const Grid &grid, const ReservoirProperties &resprop, const PressureSolution &pressure_sol)
Definition: CflCalculator.hpp:55
double findCFLtimeGravity(const Grid &grid, const ReservoirProperties &resprop, const typename Grid::Vector &gravity)
Definition: CflCalculator.hpp:91
double findCFLtimeCapillary(const Grid &grid, const ReservoirProperties &resprop)
Definition: CflCalculator.hpp:143
Definition: BlackoilFluid.hpp:32
M inverse3x3(const M &m)
Definition: MatrixInverse.hpp:86
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