Loading...
Searching...
No Matches
CflCalculator.hpp
//===========================================================================
//
// File: CflCalculator.hpp
//
// Created: Wed Jul 1 09:35:59 2009
//
// Author(s): Halvor M Nilsen <hnil@sintef.no>
// Atgeirr F Rasmussen <atgeirr@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_CFLCALCULATOR_HEADER
#define OPENRS_CFLCALCULATOR_HEADER
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/Average.hpp>
#include <opm/porsol/common/Matrix.hpp>
#include <opm/porsol/common/MatrixInverse.hpp>
namespace Opm {
namespace cfl_calculator {
template <class Grid, class ReservoirProperties, class PressureSolution>
double findCFLtimeVelocity(const Grid& grid,
const ReservoirProperties& resprop,
const PressureSolution& pressure_sol)
{
double dt = 1e100;
typename Grid::CellIterator c = grid.cellbegin();
for (; c != grid.cellend(); ++c) {
double flux_p = 0.0;
double flux_n = 0.0;
typename Grid::CellIterator::FaceIterator f = c->facebegin();
for (; f != c->faceend(); ++f) {
const double loc_flux = pressure_sol.outflux(f);
if (loc_flux > 0) {
flux_p += loc_flux;
} else {
flux_n -= loc_flux;
}
}
double flux = std::max(flux_n, flux_p);
double loc_dt = (resprop.cflFactor()*c->volume()*resprop.porosity(c->index()))/flux;
if (loc_dt == 0.0) {
OPM_THROW(std::runtime_error, "Cfl computation gave dt = 0.0");
}
if (loc_dt < dt){
dt = loc_dt;
}
}
return dt;
}
template <class Grid, class ReservoirProperties>
double findCFLtimeGravity(const Grid& grid,
const ReservoirProperties& resprop,
const typename Grid::Vector& gravity)
{
typedef typename ReservoirProperties::PermTensor PermTensor;
typedef typename ReservoirProperties::MutablePermTensor MutablePermTensor;
const int dimension = Grid::Vector::dimension;
double dt = 1e100;
typename Grid::CellIterator c = grid.cellbegin();
for (; c != grid.cellend(); ++c) {
double flux = 0.0;
typename Grid::CellIterator::FaceIterator f = c->facebegin();
for (; f != c->faceend(); ++f) {
// UGLY WARNING
MutablePermTensor loc_perm_aver;
const double* permdata = 0;
if (!f->boundary()) {
PermTensor K0 = resprop.permeability(f->cellIndex());
PermTensor K1 = resprop.permeability(f->neighbourCellIndex());
loc_perm_aver = Opm::utils::arithmeticAverage<PermTensor, MutablePermTensor>(K0, K1);
permdata = loc_perm_aver.data();
} else {
permdata = resprop.permeability(f->cellIndex()).data();
}
PermTensor loc_perm(dimension, dimension, permdata);
typename Grid::Vector loc_halfface_normal = f->normal();
double loc_gravity_flux = 0.0;
for (int k = 0; k < dimension; ++k) {
for (int q = 0; q < dimension; ++q) {
loc_gravity_flux += loc_halfface_normal[q]*(loc_perm(q,k)*gravity[k]*resprop.densityDifference());
}
}
loc_gravity_flux *= f->area();
if (loc_gravity_flux > 0) {
flux += loc_gravity_flux;
}
}
double loc_dt = (resprop.cflFactorGravity()*c->volume()*resprop.porosity(c->index()))/flux;
if (loc_dt < dt){
dt = loc_dt;
}
}
return dt;
}
template <class Grid, class ReservoirProperties>
double findCFLtimeCapillary(const Grid& grid,
const ReservoirProperties& resprop)
{
typedef typename ReservoirProperties::PermTensor PermTensor;
typedef typename ReservoirProperties::MutablePermTensor MutablePermTensor;
const int dimension = Grid::Vector::dimension;
double dt = 1e100;
typename Grid::CellIterator c = grid.cellbegin();
for (; c != grid.cellend(); ++c) {
typename Grid::CellIterator::FaceIterator f = c->facebegin();
for (; f != c->faceend(); ++f) {
// UGLY WARNING
MutablePermTensor loc_perm_aver;
const double* permdata = 0;
if (!f->boundary()) {
PermTensor K0 = resprop.permeability(f->cellIndex());
PermTensor K1 = resprop.permeability(f->neighbourCellIndex());
loc_perm_aver = Opm::utils::arithmeticAverage<PermTensor, MutablePermTensor>(K0, K1);
permdata = loc_perm_aver.data();
} else {
permdata = resprop.permeability(f->cellIndex()).data();
}
// PermTensor loc_perm(dimension, dimension, permdata);
MutablePermTensor loc_perm(dimension, dimension, permdata);
MutablePermTensor loc_perm_inv = inverse3x3(loc_perm);
typename Grid::Vector loc_centroid = f->centroid();
loc_centroid -= c->centroid();
double spatial_contrib = loc_centroid*prod(loc_perm_inv, loc_centroid);
double loc_dt = spatial_contrib/resprop.cflFactorCapillary();
dt = std::min(dt, loc_dt);
}
}
return dt;
}
} // namespace cfl_calculator
} // namespace Opm
#endif // OPENRS_CFLCALCULATOR_HEADER
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::prod
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:668