OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branches 'origin/master' and 'upstream/master'

Browse Source
This commit is contained in:
Markus Blatt
2013-01-18 13:48:55 +01:00
parent 87a5ba847e 08ecc6988c 862c489cc3
commit 0baa40fb88
12 changed files with 1518 additions and 444 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@@ -50,6 +50,9 @@ namespace Opm
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (deck.hasField("ENDSCALE") && threephase_model != "simple") {
THROW("Sorry, end point scaling currently available for the 'simple' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr

19
opm/core/fluid/SaturationPropsFromDeck.hpp
View File

@@ -107,9 +107,28 @@ namespace Opm
std::vector<SatFuncSet> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1
struct { // End point scaling parameters
std::vector<double> swl_;
std::vector<double> swcr_;
std::vector<double> swu_;
std::vector<double> sowcr_;
std::vector<double> krw_;
std::vector<double> krwr_;
std::vector<double> kro_;
std::vector<double> krorw_;
} eps_;
bool do_eps_; // ENDSCALE is active
bool do_3pt_; // SCALECRS: YES~true NO~false
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
void initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam);
void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const;
};

327
opm/core/fluid/SaturationPropsFromDeck_impl.hpp
View File

@@ -96,6 +96,35 @@ namespace Opm
for (int table = 0; table < num_tables; ++table) {
satfuncset_[table].init(deck, table, phase_usage_, samples);
}
// Saturation table scaling
do_eps_ = false;
do_3pt_ = false;
if (deck.hasField("ENDSCALE")) {
if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) {
THROW("Currently endpoint-scaling limited to oil-water systems without gas.");
}
if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) {
THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted.");
}
if (deck.getENDSCALE().revers_switch_ != std::string("REVERS")) {
THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'REVERS' accepted.");
}
if (deck.hasField("SCALECRS")) {
if (deck.getSCALECRS().scalecrs_ == std::string("YES")) {
do_3pt_ = true;
}
}
do_eps_ = true;
initEPS(deck, grid, std::string("SWCR"), eps_.swcr_);
initEPS(deck, grid, std::string("SWL"), eps_.swl_);
initEPS(deck, grid, std::string("SWU"), eps_.swu_);
initEPS(deck, grid, std::string("SOWCR"), eps_.sowcr_);
initEPS(deck, grid, std::string("KRW"), eps_.krw_);
initEPS(deck, grid, std::string("KRWR"), eps_.krwr_);
initEPS(deck, grid, std::string("KRO"), eps_.kro_);
initEPS(deck, grid, std::string("KRORW"), eps_.krorw_);
}
}
@@ -134,12 +163,20 @@ namespace Opm
if (dkrds) {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
if (do_eps_) {
relpermEPS(s + np*i, cells[i], kr + np*i, dkrds + np*np*i);
} else {
funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i);
}
}
} else {
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
if (do_eps_) {
relpermEPS(s + np*i, cells[i], kr + np*i);
} else {
funcForCell(cells[i]).evalKr(s + np*i, kr + np*i);
}
}
}
}
@@ -214,7 +251,293 @@ namespace Opm
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& scaleparam)
{
bool useKeyword = deck.hasField(keyword);
bool hasENPTVD = deck.hasField("ENPTVD");
bool hasENKRVD = deck.hasField("ENKRVD");
int itab = 0;
std::vector<std::vector<std::vector<double> > > table_dummy;
std::vector<std::vector<std::vector<double> > >& table = table_dummy;
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
if (keyword[0] == 'S' && (useKeyword || hasENPTVD)) {
if (keyword == std::string("SWL")) {
if (useKeyword || deck.getENPTVD().mask_[0]) {
itab = 1;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
}
} else if (keyword == std::string("SWCR")) {
if (useKeyword || deck.getENPTVD().mask_[1]) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).swcr_;
}
} else if (keyword == std::string("SWU")) {
if (useKeyword || deck.getENPTVD().mask_[2]) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
}
} else if (keyword == std::string("SOWCR")) {
if (useKeyword || deck.getENPTVD().mask_[3]) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sowcr_;
}
}else {
THROW(" -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
table = deck.getENPTVD().table_;
}
} else if (keyword[0] == 'K' && (useKeyword || hasENKRVD)) {
if (keyword == std::string("KRW")) {
if (useKeyword || deck.getENKRVD().mask_[0]) {
itab = 1;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
}
} else if (keyword == std::string("KRO")) {
if (useKeyword || deck.getENKRVD().mask_[1]) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
}
} else if (keyword == std::string("KRWR")) {
if (useKeyword || deck.getENKRVD().mask_[2]) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
}
} else if (keyword == std::string("KRORW")) {
if (useKeyword || deck.getENKRVD().mask_[3]) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
}
} else {
THROW(" -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
table = deck.getENKRVD().table_;
}
}
if (scaleparam.empty()) {
return;
} else if (useKeyword) {
// Keyword values from deck
std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
const int* gc = grid.global_cell;
const std::vector<double>& val = deck.getFloatingPointValue(keyword);
for (int c = 0; c < int(scaleparam.size()); ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c];
scaleparam[c] = val[deck_pos];
}
} else {
std::cout << "--- Scaling parameter '" << keyword << "' assigned via ";
if (keyword[0] == 'S')
deck.getENPTVD().write(std::cout);
else
deck.getENKRVD().write(std::cout);
const double* cc = grid.cell_centroids;
const int dim = grid.dimensions;
for (int cell = 0; cell < grid.number_of_cells; ++cell) {
int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
if (table[itab][jtab][0] != -1.0) {
std::vector<double>& depth = table[0][jtab];
std::vector<double>& val = table[itab][jtab];
double zc = cc[dim*cell+dim-1];
if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth, val, zc);
}
}
}
}
}
// Saturation scaling
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::relpermEPS(const double *s, const int cell, double *kr, double *dkrds) const
{
const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
double ss[PhaseUsage::MaxNumPhases];
if (do_3pt_) { // Three-point scaling
// Transforms for water saturation
if (eps_.swcr_.empty() && eps_.swu_.empty()) {
ss[wpos] = s[wpos];
} else {
double s_r = 1.0-funcForCell(cell).sowcr_;
double sr = eps_.sowcr_.empty() ? s_r : 1.0-eps_.sowcr_[cell];
if (s[wpos] <= sr) {
double sw_cr = funcForCell(cell).swcr_;
double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell];
ss[wpos] = (s[wpos] <= swcr) ? sw_cr : sw_cr+(s[wpos]-swcr)*(s_r-sw_cr)/(sr-swcr);
} else {
double sw_max = funcForCell(cell).smax_[wpos];
double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell];
ss[wpos] = (s[wpos] >= swmax) ? sw_max : s_r+(s[wpos]-sr)*(sw_max-s_r)/(swmax-sr);
}
}
// Transforms for oil saturation
if (eps_.sowcr_.empty() && eps_.swl_.empty()) {
ss[opos] = s[opos];
} else {
double s_r = 1.0-funcForCell(cell).swcr_;
double sr = eps_.swcr_.empty() ? s_r : 1.0-eps_.swcr_[cell];
if (s[opos] <= sr) {
double sow_cr = funcForCell(cell).sowcr_;
double sowcr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell];
ss[opos] = (s[opos] <= sowcr) ? sow_cr : sow_cr+(s[opos]-sowcr)*(s_r-sow_cr)/(sr-sowcr);
} else {
double sow_max = funcForCell(cell).smax_[opos];
double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]);
ss[opos] = (s[opos] >= sowmax) ? sow_max : s_r+(s[opos]-sr)*(sow_max-s_r)/(sowmax-sr);
}
}
} else { // Two-point scaling
// Transforms for water saturation
if (eps_.swcr_.empty() && eps_.swu_.empty()) {
ss[wpos] = s[wpos];
} else {
double sw_cr = funcForCell(cell).swcr_;
double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell];
if (s[wpos] <= swcr) {
ss[wpos] = sw_cr;
} else {
double sw_max = funcForCell(cell).smax_[wpos];
double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell];
ss[wpos] = (s[wpos] >= swmax) ? sw_max : sw_cr + (s[wpos]-swcr)*(sw_max-sw_cr)/(swmax-swcr);
}
}
// Transforms for oil saturation
if (eps_.sowcr_.empty() && eps_.swl_.empty()) {
ss[opos] = s[opos];
} else {
double sow_cr = funcForCell(cell).sowcr_;
double socr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell];
if (s[opos] <= socr) {
ss[opos] = sow_cr;
} else {
double sow_max = funcForCell(cell).smax_[opos];
double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]);
ss[opos] = (s[opos] >= sowmax) ? sow_max : sow_cr + (s[opos]-socr) *(sow_max-sow_cr)/(sowmax-socr);
}
}
}
// Evaluation of relperms
if (dkrds) {
THROW("Relperm derivatives not yet available in combination with end point scaling ...");
funcForCell(cell).evalKrDeriv(ss, kr, dkrds);
} else {
// Assume: sw_cr -> krw=0 sw_max -> krw=<max water relperm>
// sow_cr -> kro=0 sow_max -> kro=<max oil relperm>
funcForCell(cell).evalKr(ss, kr);
}
// Scaling of relperms values
// - Water
if (eps_.krw_.empty() && eps_.krwr_.empty()) { // No value scaling
} else if (eps_.krwr_.empty()) { // Two-point
kr[wpos] *= (eps_.krw_[cell]/funcForCell(cell).krwmax_);
} else {
double swcr = eps_.swcr_.empty() ? funcForCell(cell).swcr_ : eps_.swcr_[cell];
double swmax = eps_.swu_.empty() ? funcForCell(cell).smax_[wpos] : eps_.swu_[cell];
double sr;
if (do_3pt_) {
sr = eps_.sowcr_.empty() ? 1.0-funcForCell(cell).sowcr_ : 1.0-eps_.sowcr_[cell];
} else {
double sw_cr = funcForCell(cell).swcr_;
double sw_max = funcForCell(cell).smax_[wpos];
double s_r = 1.0-funcForCell(cell).sowcr_;
sr = swcr + (s_r-sw_cr)*(swmax-swcr)/(sw_max-sw_cr);
}
if (s[wpos] <= swcr) {
kr[wpos] = 0.0;
} else if (sr > swmax-1.0e-6) {
if (do_3pt_) { //Ignore krw and do two-point?
kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_;
} else if (!eps_.kro_.empty()){ //Ignore krwr and do two-point
kr[wpos] *= eps_.krw_[cell]/funcForCell(cell).krwmax_;
}
} else if (s[wpos] <= sr) {
kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_;
} else if (s[wpos] <= swmax) {
double krw_max = funcForCell(cell).krwmax_;
double krw = eps_.krw_.empty() ? krw_max : eps_.krw_[cell];
double krw_r = funcForCell(cell).krwr_;
double krwr = eps_.krwr_.empty() ? krw_r : eps_.krwr_[cell];
if (std::fabs(krw_max- krw_r) > 1.0e-6) {
kr[wpos] = krwr + (kr[wpos]-krw_r)*(krw-krwr)/(krw_max-krw_r);
} else {
kr[wpos] = krwr + (krw-krwr)*(s[wpos]-sr)/(swmax-sr);
}
} else {
kr[wpos] = eps_.krw_.empty() ? funcForCell(cell).krwmax_ : eps_.krw_[cell];
}
}
// - Oil
if (eps_.kro_.empty() && eps_.krorw_.empty()) { // No value scaling
} else if (eps_.krorw_.empty()) { // Two-point scaling
kr[opos] *= (eps_.kro_[cell]/funcForCell(cell).kromax_);
} else {
double sowcr = eps_.sowcr_.empty() ? funcForCell(cell).sowcr_ : eps_.sowcr_[cell];
double sowmax = eps_.swl_.empty() ? funcForCell(cell).smax_[opos] : 1.0-eps_.swl_[cell];
double sr;
if (do_3pt_) {
sr = eps_.swcr_.empty() ? 1.0-funcForCell(cell).swcr_ : 1.0-eps_.swcr_[cell];
} else {
double sow_cr = funcForCell(cell).sowcr_;
double sow_max = funcForCell(cell).smax_[opos];
double s_r = 1.0-funcForCell(cell).swcr_;
sr = sowcr + (s_r-sow_cr)*(sowmax-sowcr)/(sow_max-sow_cr);
}
if (s[opos] <= sowcr) {
kr[opos] = 0.0;
} else if (sr > sowmax-1.0e-6) {
if (do_3pt_) { //Ignore kro and do two-point?
kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_;
} else if (!eps_.kro_.empty()){ //Ignore krowr and do two-point
kr[opos] *= eps_.kro_[cell]/funcForCell(cell).kromax_;
}
} else if (s[opos] <= sr) {
kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_;
} else if (s[opos] <= sowmax) {
double kro_max = funcForCell(cell).kromax_;
double kro = eps_.kro_.empty() ? kro_max : eps_.kro_[cell];
double kro_rw = funcForCell(cell).krorw_;
double krorw = eps_.krorw_[cell];
if (std::fabs(kro_max- kro_rw) > 1.0e-6) {
kr[opos] = krorw + (kr[opos]- kro_rw)*(kro-krorw)/(kro_max- kro_rw);
} else {
kr[opos] = krorw + (kro-krorw)*(s[opos]-sr)/(sowmax-sr);
}
} else {
kr[opos] = eps_.kro_.empty() ? funcForCell(cell).kromax_ : eps_.kro_[cell];
}
}
}
} // namespace Opm

310
opm/core/transport/reorder/DGBasis.cpp Normal file
View File

@@ -0,0 +1,310 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/DGBasis.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <numeric>
namespace Opm
{
/// Virtual destructor.
DGBasisInterface::~DGBasisInterface()
{
}
// ---------------- Methods for class DGBasisBoundedTotalDegree ----------------
/// Constructor.
/// \param[in] grid grid on which basis is used (cell-wise)
/// \param[in] degree polynomial degree of basis
DGBasisBoundedTotalDegree::DGBasisBoundedTotalDegree(const UnstructuredGrid& grid,
const int degree_arg)
: grid_(grid),
degree_(degree_arg)
{
if (grid_.dimensions > 3) {
THROW("Grid dimension must be 1, 2 or 3.");
}
if (degree_ > 1 || degree_ < 0) {
THROW("Degree must be 0 or 1.");
}
}
/// Destructor.
DGBasisBoundedTotalDegree::~DGBasisBoundedTotalDegree()
{
}
/// The number of basis functions per cell.
int DGBasisBoundedTotalDegree::numBasisFunc() const
{
switch (dimensions()) {
case 1:
return degree_ + 1;
case 2:
return (degree_ + 2)*(degree_ + 1)/2;
case 3:
return (degree_ + 3)*(degree_ + 2)*(degree_ + 1)/6;
default:
THROW("Dimensions must be 1, 2 or 3.");
}
}
/// The number of space dimensions.
int DGBasisBoundedTotalDegree::dimensions() const
{
return grid_.dimensions;
}
/// The polynomial degree of the basis functions.
int DGBasisBoundedTotalDegree::degree() const
{
return degree_;
}
/// Evaluate all basis functions associated with cell at x,
/// writing to f_x. The array f_x must have size equal to
/// numBasisFunc().
void DGBasisBoundedTotalDegree::eval(const int cell,
const double* x,
double* f_x) const
{
const int dim = dimensions();
const double* cc = grid_.cell_centroids + dim*cell;
// Note intentional fallthrough in this switch statement!
switch (degree_) {
case 1:
for (int ix = 0; ix < dim; ++ix) {
f_x[1 + ix] = x[ix] - cc[ix];
}
case 0:
f_x[0] = 1;
break;
default:
THROW("Maximum degree is 1 for now.");
}
}
/// Evaluate gradients of all basis functions associated with
/// cell at x, writing to grad_f_x. The array grad_f_x must
/// have size numBasisFunc() * dimensions(). The dimensions()
/// components of the first basis function gradient come
/// before the components of the second etc.
void DGBasisBoundedTotalDegree::evalGrad(const int /*cell*/,
const double* /*x*/,
double* grad_f_x) const
{
const int dim = dimensions();
const int num_basis = numBasisFunc();
std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
if (degree_ == 1) {
for (int ix = 0; ix < dim; ++ix) {
grad_f_x[dim*(ix + 1) + ix] = 1.0;
}
}
}
/// Modify basis coefficients to add to the function value.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to (f + increment) by modifying the c_i. This is done without
/// modifying its gradient.
/// \param[in] increment Add this value to the function.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
void DGBasisBoundedTotalDegree::addConstant(const double increment,
double* coefficients) const
{
coefficients[0] += increment;
}
/// Modify basis coefficients to change the function's slope.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to a function g with the property that grad g = factor * grad f
/// by modifying the c_i. This is done without modifying the average,
/// i.e. the integrals of g and f over the cell are the same.
/// \param[in] factor Multiply gradient by this factor.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
void DGBasisBoundedTotalDegree::multiplyGradient(const double factor,
double* coefficients) const
{
const int nb = numBasisFunc();
for (int ix = 1; ix < nb; ++ix) {
coefficients[ix] *= factor;
}
}
// ---------------- Methods for class DGBasisMultilin ----------------
/// Constructor.
/// \param[in] grid grid on which basis is used (cell-wise)
/// \param[in] degree polynomial degree of basis
DGBasisMultilin::DGBasisMultilin(const UnstructuredGrid& grid,
const int degree_arg)
: grid_(grid),
degree_(degree_arg)
{
if (grid_.dimensions > 3) {
THROW("Grid dimension must be 1, 2 or 3.");
}
if (degree_ > 1 || degree_ < 0) {
THROW("Degree must be 0 or 1.");
}
}
/// Destructor.
DGBasisMultilin::~DGBasisMultilin()
{
}
/// The number of basis functions per cell.
int DGBasisMultilin::numBasisFunc() const
{
switch (dimensions()) {
case 1:
return degree_ + 1;
case 2:
return (degree_ + 1)*(degree_ + 1);
case 3:
return (degree_ + 1)*(degree_ + 1)*(degree_ + 1);
default:
THROW("Dimensions must be 1, 2 or 3.");
}
}
/// The number of space dimensions.
int DGBasisMultilin::dimensions() const
{
return grid_.dimensions;
}
/// The polynomial degree of the basis functions.
int DGBasisMultilin::degree() const
{
return degree_;
}
/// Evaluate all basis functions associated with cell at x,
/// writing to f_x. The array f_x must have size equal to
/// numBasisFunc().
void DGBasisMultilin::eval(const int cell,
const double* x,
double* f_x) const
{
const int dim = dimensions();
const int num_basis = numBasisFunc();
const double* cc = grid_.cell_centroids + dim*cell;
switch (degree_) {
case 0:
f_x[0] = 1;
break;
case 1:
std::fill(f_x, f_x + num_basis, 1.0);
for (int dd = 0; dd < dim; ++dd) {
const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] };
const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example.
for (int ix = 0; ix < num_basis; ++ix) {
f_x[ix] *= f[(ix/divi) % 2];
}
}
break;
default:
THROW("Maximum degree is 1 for now.");
}
}
/// Evaluate gradients of all basis functions associated with
/// cell at x, writing to grad_f_x. The array grad_f_x must
/// have size numBasisFunc() * dimensions(). The dimensions()
/// components of the first basis function gradient come
/// before the components of the second etc.
void DGBasisMultilin::evalGrad(const int cell,
const double* x,
double* grad_f_x) const
{
const int dim = dimensions();
const int num_basis = numBasisFunc();
const double* cc = grid_.cell_centroids + dim*cell;
switch (degree_) {
case 0:
std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
break;
case 1:
std::fill(grad_f_x, grad_f_x + num_basis*dim, 1.0);
for (int dd = 0; dd < dim; ++dd) {
const double f[2] = { 0.5 - x[dd] + cc[dd], 0.5 + x[dd] - cc[dd] };
const double fder[2] = { -1.0, 1.0 };
const int divi = 1 << (dim - dd - 1); // { 4, 2, 1 } for 3d, for example.
for (int ix = 0; ix < num_basis; ++ix) {
const int ind = (ix/divi) % 2;
for (int dder = 0; dder < dim; ++dder) {
grad_f_x[ix*dim + dder] *= (dder == dd ? fder[ind] : f[ind]);
}
}
}
break;
default:
THROW("Maximum degree is 1 for now.");
}
}
/// Modify basis coefficients to add to the function value.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to (f + increment) by modifying the c_i. This is done without
/// modifying its gradient.
/// \param[in] increment Add this value to the function.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
void DGBasisMultilin::addConstant(const double increment,
double* coefficients) const
{
const int nb = numBasisFunc();
const double term = increment/double(nb);
for (int ix = 0; ix < nb; ++ix) {
coefficients[ix] += term;
}
}
/// Modify basis coefficients to change the function's slope.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to a function g with the property that grad g = factor * grad f
/// by modifying the c_i. This is done without modifying the average,
/// i.e. the integrals of g and f over the cell are the same.
/// \param[in] factor Multiply gradient by this factor.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
void DGBasisMultilin::multiplyGradient(const double factor,
double* coefficients) const
{
const int nb = numBasisFunc();
const double average = std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb);
for (int ix = 0; ix < nb; ++ix) {
coefficients[ix] = factor*(coefficients[ix] - average) + average;
}
}
} // namespace Opm

232
opm/core/transport/reorder/DGBasis.hpp Normal file
View File

@@ -0,0 +1,232 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_DGBASIS_HEADER_INCLUDED
#define OPM_DGBASIS_HEADER_INCLUDED
struct UnstructuredGrid;
namespace Opm
{
/// Base class for Discontinuous Galerkin bases, intended for time-of-flight computations.
class DGBasisInterface
{
public:
/// Virtual destructor.
virtual ~DGBasisInterface();
/// The number of basis functions per cell.
virtual int numBasisFunc() const = 0;
/// The number of space dimensions.
virtual int dimensions() const = 0;
/// The polynomial degree of the basis functions.
virtual int degree() const = 0;
/// Evaluate all basis functions associated with cell at x,
/// writing to f_x. The array f_x must have size equal to
/// numBasisFunc().
virtual void eval(const int cell,
const double* x,
double* f_x) const = 0;
/// Evaluate gradients of all basis functions associated with
/// cell at x, writing to grad_f_x. The array grad_f_x must
/// have size numBasisFunc() * dimensions(). The dimensions()
/// components of the first basis function gradient come
/// before the components of the second etc.
virtual void evalGrad(const int cell,
const double* x,
double* grad_f_x) const = 0;
/// Modify basis coefficients to add to the function value.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to (f + increment) by modifying the c_i. This is done without
/// modifying its gradient.
/// \param[in] increment Add this value to the function.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
virtual void addConstant(const double increment,
double* coefficients) const = 0;
/// Modify basis coefficients to change the function's slope.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to a function g with the property that grad g = factor * grad f
/// by modifying the c_i. This is done without modifying the average,
/// i.e. the integrals of g and f over the cell are the same.
/// \param[in] factor Multiply gradient by this factor.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
virtual void multiplyGradient(const double factor,
double* coefficients) const = 0;
};
/// A class providing discontinuous Galerkin basis functions
/// of bounded total degree.
///
/// The basis functions are the following for each cell (example for 3d):
/// Degree 0: 1.
/// Degree 1: 1, x - xc, y - yc, z - zc
/// where (xc, yc, zc) are the coordinates of the cell centroid.
/// Further degrees await development.
class DGBasisBoundedTotalDegree : public DGBasisInterface
{
public:
/// Constructor.
/// \param[in] grid grid on which basis is used (cell-wise)
/// \param[in] degree polynomial degree of basis
DGBasisBoundedTotalDegree(const UnstructuredGrid& grid, const int degree);
/// Destructor.
virtual ~DGBasisBoundedTotalDegree();
/// The number of basis functions per cell.
virtual int numBasisFunc() const;
/// The number of space dimensions.
virtual int dimensions() const;
/// The polynomial degree of the basis functions.
virtual int degree() const;
/// Evaluate all basis functions associated with cell at x,
/// writing to f_x. The array f_x must have size equal to
/// numBasisFunc().
virtual void eval(const int cell,
const double* x,
double* f_x) const;
/// Evaluate gradients of all basis functions associated with
/// cell at x, writing to grad_f_x. The array grad_f_x must
/// have size numBasisFunc() * dimensions(). The dimensions()
/// components of the first basis function gradient come
/// before the components of the second etc.
virtual void evalGrad(const int cell,
const double* x,
double* grad_f_x) const;
/// Modify basis coefficients to add to the function value.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to (f + increment) by modifying the c_i. This is done without
/// modifying its gradient.
/// \param[in] increment Add this value to the function.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
virtual void addConstant(const double increment,
double* coefficients) const;
/// Modify basis coefficients to change the function's slope.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to a function g with the property that grad g = factor * grad f
/// by modifying the c_i. This is done without modifying the average,
/// i.e. the integrals of g and f over the cell are the same.
/// \param[in] factor Multiply gradient by this factor.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
virtual void multiplyGradient(const double factor,
double* coefficients) const;
private:
const UnstructuredGrid& grid_;
const int degree_;
};
/// A class providing discontinuous Galerkin basis functions of
/// multi-degree 1 (bilinear or trilinear functions).
///
/// The basis functions for a cell are the following
/// Degree 0: 1.
/// (for 2 dims:)
/// (Bi)degree 1: (x-)(y-), (x-)(y+), (x+)(y-), (x+)(y+)
/// where (x-) = (1/2 - x + xc), (x+) = (1/2 + x - xc)
/// and xc is the x-coordinate of the cell centroid.
/// Similar for (y-), (y+).
class DGBasisMultilin : public DGBasisInterface
{
public:
/// Constructor.
/// \param[in] grid grid on which basis is used (cell-wise)
/// \param[in] degree polynomial degree of basis (in each coordinate)
DGBasisMultilin(const UnstructuredGrid& grid, const int degree);
/// Destructor.
virtual ~DGBasisMultilin();
/// The number of basis functions per cell.
virtual int numBasisFunc() const;
/// The number of space dimensions.
virtual int dimensions() const;
/// The polynomial degree of the basis functions.
virtual int degree() const;
/// Evaluate all basis functions associated with cell at x,
/// writing to f_x. The array f_x must have size equal to
/// numBasisFunc().
virtual void eval(const int cell,
const double* x,
double* f_x) const;
/// Evaluate gradients of all basis functions associated with
/// cell at x, writing to grad_f_x. The array grad_f_x must
/// have size numBasisFunc() * dimensions(). The dimensions()
/// components of the first basis function gradient come
/// before the components of the second etc.
virtual void evalGrad(const int cell,
const double* x,
double* grad_f_x) const;
/// Modify basis coefficients to add to the function value.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to (f + increment) by modifying the c_i. This is done without
/// modifying its gradient.
/// \param[in] increment Add this value to the function.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
virtual void addConstant(const double increment,
double* coefficients) const;
/// Modify basis coefficients to change the function's slope.
/// A function f = sum_i c_i b_i is assumed, and we change
/// it to a function g with the property that grad g = factor * grad f
/// by modifying the c_i. This is done without modifying the average,
/// i.e. the integrals of g and f over the cell are the same.
/// \param[in] factor Multiply gradient by this factor.
/// \param[out] coefficients Coefficients {c_i} for a single cell.
virtual void multiplyGradient(const double factor,
double* coefficients) const;
private:
const UnstructuredGrid& grid_;
const int degree_;
};
} // namespace Opm
#endif // OPM_DGBASIS_HEADER_INCLUDED

27
opm/core/transport/reorder/TransportModelInterface.cpp
View File

@@ -29,15 +29,18 @@
void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux)
{
// Compute reordered sequence of single-cell problems
std::vector<int> sequence(grid.number_of_cells);
std::vector<int> components(grid.number_of_cells + 1);
sequence_.resize(grid.number_of_cells);
components_.resize(grid.number_of_cells + 1);
int ncomponents;
time::StopWatch clock;
clock.start();
compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
compute_sequence(&grid, darcyflux, &sequence_[0], &components_[0], &ncomponents);
clock.stop();
std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl;
// Make vector's size match actual used data.
components_.resize(ncomponents + 1);
// Invoke appropriate solve method for each interdependent component.
for (int comp = 0; comp < ncomponents; ++comp) {
#if 0
@@ -50,11 +53,23 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g
}
#endif
#endif
const int comp_size = components[comp + 1] - components[comp];
const int comp_size = components_[comp + 1] - components_[comp];
if (comp_size == 1) {
solveSingleCell(sequence[components[comp]]);
solveSingleCell(sequence_[components_[comp]]);
} else {
solveMultiCell(comp_size, &sequence[components[comp]]);
solveMultiCell(comp_size, &sequence_[components_[comp]]);
}
}
}
const std::vector<int>& Opm::TransportModelInterface::sequence() const
{
return sequence_;
}
const std::vector<int>& Opm::TransportModelInterface::components() const
{
return components_;
}

10
opm/core/transport/reorder/TransportModelInterface.hpp
View File

@@ -20,6 +20,8 @@
#ifndef OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
#define OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
#include <vector>
struct UnstructuredGrid;
namespace Opm
@@ -31,7 +33,8 @@ namespace Opm
/// method that will have an interface geared to the model's
/// needs. (The solve() method is therefore not virtual in this
/// class.) The reorderAndTransport() method is provided as an aid
/// to implementing solve() in subclasses.
/// to implementing solve() in subclasses, together with the
/// sequence() and components() methods for accessing the ordering.
class TransportModelInterface
{
public:
@@ -41,6 +44,11 @@ namespace Opm
virtual void solveMultiCell(const int num_cells, const int* cells) = 0;
protected:
void reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux);
const std::vector<int>& sequence() const;
const std::vector<int>& components() const;
private:
std::vector<int> sequence_;
std::vector<int> components_;
};

110
opm/core/transport/reorder/TransportModelTracerTof.cpp
View File

@@ -32,7 +32,14 @@ namespace Opm
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid,
const bool use_multidim_upwind)
: grid_(grid), use_multidim_upwind_(use_multidim_upwind)
: grid_(grid),
darcyflux_(0),
porevolume_(0),
source_(0),
tof_(0),
tracer_(0),
num_tracers_(0),
use_multidim_upwind_(use_multidim_upwind)
{
}
@@ -68,6 +75,62 @@ namespace Opm
face_tof_.resize(grid_.number_of_faces);
std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
}
num_tracers_ = 0;
reorderAndTransport(grid_, darcyflux);
}
/// Solve for time-of-flight and a number of tracers.
/// One tracer will be used for each inflow flux specified in
/// the source parameter.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values (1 per cell).
/// \param[out] tracer Array of tracer values (N per cell, where N is
/// the number of cells c for which source[c] > 0.0).
void TransportModelTracerTof::solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof,
std::vector<double>& tracer)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
}
#endif
tof.resize(grid_.number_of_cells);
std::fill(tof.begin(), tof.end(), 0.0);
tof_ = &tof[0];
// Find the tracer heads (injectors).
std::vector<int> tracerheads;
for (int c = 0; c < grid_.number_of_cells; ++c) {
if (source[c] > 0.0) {
tracerheads.push_back(c);
}
}
num_tracers_ = tracerheads.size();
tracer.resize(grid_.number_of_cells*num_tracers_);
std::fill(tracer.begin(), tracer.end(), 0.0);
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer[tracerheads[tr]*num_tracers_ + tr] = 1.0;
}
tracer_ = &tracer[0];
if (use_multidim_upwind_) {
face_tof_.resize(grid_.number_of_faces);
std::fill(face_tof_.begin(), face_tof_.end(), 0.0);
THROW("Multidimensional upwind not yet implemented for tracer.");
}
reorderAndTransport(grid_, darcyflux);
}
@@ -100,18 +163,32 @@ namespace Opm
flux =-darcyflux_[f];
other = grid_.face_cells[2*f];
}
// Add flux to upwind_term or downwind_flux, if interior.
if (other != -1) {
if (flux < 0.0) {
// Add flux to upwind_term or downwind_flux
if (flux < 0.0) {
// Using tof == 0 on inflow, so we only add a
// nonzero contribution if we are on an internal
// face.
if (other != -1) {
upwind_term += flux*tof_[other];
} else {
downwind_flux += flux;
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer_[num_tracers_*cell + tr] += flux*tracer_[num_tracers_*other + tr];
}
}
} else {
downwind_flux += flux;
}
}
// Compute tof.
tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
// Compute tracers (if any).
// Do not change tracer solution in source cells.
if (source_[cell] <= 0.0) {
for (int tr = 0; tr < num_tracers_; ++tr) {
tracer_[num_tracers_*cell + tr] *= -1.0/downwind_flux;
}
}
}
@@ -131,25 +208,20 @@ namespace Opm
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
int f = grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == grid_.face_cells[2*f]) {
flux = darcyflux_[f];
other = grid_.face_cells[2*f+1];
} else {
flux =-darcyflux_[f];
other = grid_.face_cells[2*f];
}
// Add flux to upwind_term or downwind_flux, if interior.
if (other != -1) {
if (flux < 0.0) {
upwind_term += flux*face_tof_[f];
} else {
double fterm, cterm_factor;
multidimUpwindTerms(f, cell, fterm, cterm_factor);
downwind_term_face += fterm*flux;
downwind_term_cell_factor += cterm_factor*flux;
}
// Add flux to upwind_term or downwind_term_[face|cell_factor].
if (flux < 0.0) {
upwind_term += flux*face_tof_[f];
} else {
double fterm, cterm_factor;
multidimUpwindTerms(f, cell, fterm, cterm_factor);
downwind_term_face += fterm*flux;
downwind_term_cell_factor += cterm_factor*flux;
}
}

19
opm/core/transport/reorder/TransportModelTracerTof.hpp
View File

@@ -59,6 +59,23 @@ namespace Opm
const double* source,
std::vector<double>& tof);
/// Solve for time-of-flight and a number of tracers.
/// One tracer will be used for each inflow flux specified in
/// the source parameter.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values (1 per cell).
/// \param[out] tracer Array of tracer values (N per cell, where N is
/// the number of cells c for which source[c] > 0.0).
void solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof,
std::vector<double>& tracer);
private:
virtual void solveSingleCell(const int cell);
void solveSingleCellMultidimUpwind(const int cell);
@@ -73,6 +90,8 @@ namespace Opm
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
double* tof_;
double* tracer_;
int num_tracers_;
bool use_multidim_upwind_;
std::vector<double> face_tof_; // For multidim upwind face tofs.
mutable std::vector<int> adj_faces_; // For multidim upwind logic.

827
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
View File

@@ -17,10 +17,14 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/grid/CellQuadrature.hpp>
#include <opm/core/grid/FaceQuadrature.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/transport/reorder/DGBasis.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/VelocityInterpolation.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/linalg/blas_lapack.h>
#include <algorithm>
#include <cmath>
@@ -30,399 +34,80 @@ namespace Opm
{
// --------------- Helpers for TransportModelTracerTofDiscGal ---------------
/// A class providing discontinuous Galerkin basis functions.
struct DGBasis
{
static int numBasisFunc(const int dimensions,
const int degree)
{
switch (dimensions) {
case 1:
return degree + 1;
case 2:
return (degree + 2)*(degree + 1)/2;
case 3:
return (degree + 3)*(degree + 2)*(degree + 1)/6;
default:
THROW("Dimensions must be 1, 2 or 3.");
}
}
/// Evaluate all nonzero basis functions at x,
/// writing to f_x. The array f_x must have
/// size numBasisFunc(grid.dimensions, degree).
///
/// The basis functions are the following
/// Degree 0: 1.
/// Degree 1: x - xc, y - yc, z - zc etc.
/// Further degrees await development.
static void eval(const UnstructuredGrid& grid,
const int cell,
const int degree,
const double* x,
double* f_x)
{
const int dim = grid.dimensions;
const double* cc = grid.cell_centroids + dim*cell;
// Note intentional fallthrough in this switch statement!
switch (degree) {
case 1:
for (int ix = 0; ix < dim; ++ix) {
f_x[1 + ix] = x[ix] - cc[ix];
}
case 0:
f_x[0] = 1;
break;
default:
THROW("Maximum degree is 1 for now.");
}
}
/// Evaluate gradients of all nonzero basis functions at x,
/// writing to grad_f_x. The array grad_f_x must have size
/// numBasisFunc(grid.dimensions, degree) * grid.dimensions.
/// The <grid.dimensions> components of the first basis function
/// gradient come before the components of the second etc.
static void evalGrad(const UnstructuredGrid& grid,
const int /*cell*/,
const int degree,
const double* /*x*/,
double* grad_f_x)
{
const int dim = grid.dimensions;
const int num_basis = numBasisFunc(dim, degree);
std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
if (degree > 1) {
THROW("Maximum degree is 1 for now.");
} else if (degree == 1) {
for (int ix = 0; ix < dim; ++ix) {
grad_f_x[dim*(ix + 1) + ix] = 1.0;
}
}
}
};
static void cross(const double* a, const double* b, double* res)
{
res[0] = a[1]*b[2] - a[2]*b[1];
res[1] = a[2]*b[0] - a[0]*b[2];
res[2] = a[0]*b[1] - a[1]*b[0];
}
static double triangleArea3d(const double* p0,
const double* p1,
const double* p2)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double cr[3];
cross(a, b, cr);
return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]);
}
/// Calculates the determinant of a 3 x 3 matrix, represented as
/// three three-dimensional arrays.
static double determinantOf(const double* a0,
const double* a1,
const double* a2)
{
return
a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) -
a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) +
a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]);
}
/// Computes the volume of a tetrahedron consisting of 4 vertices
/// with 3-dimensional coordinates
static double tetVolume(const double* p0,
const double* p1,
const double* p2,
const double* p3)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
return std::fabs(determinantOf(a, b, c) / 6.0);
}
/// A class providing numerical quadrature for cells.
/// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by cell volume,
/// so weights always sum to cell volume.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = cell volume, x_0 = cell centroid
/// Degree 2 method:
/// Based on subdivision of each cell face into triangles
/// with the face centroid as a common vertex, and then
/// subdividing the cell into tetrahedra with the cell
/// centroid as a common vertex. Then apply the tetrahedron
/// rule with the following 4 nodes (uniform weights):
/// a = 0.138196601125010515179541316563436
/// x_i has all barycentric coordinates = a, except for
/// the i'th coordinate which is = 1 - 3a.
/// This rule is from http://nines.cs.kuleuven.be/ecf,
/// it is the second degree 2 4-point rule for tets,
/// referenced to Stroud(1971).
/// The tetrahedra are numbered T_{i,j}, and are given by the
/// cell centroid C, the face centroid FC_i, and two nodes
/// of face i: FN_{i,j}, FN_{i,j+1}.
class CellQuadrature
{
public:
CellQuadrature(const UnstructuredGrid& grid,
const int cell,
const int degree)
: grid_(grid), cell_(cell), degree_(degree)
{
if (degree > 2) {
THROW("CellQuadrature exact for polynomial degrees > 1 not implemented.");
}
if (degree == 2) {
// Prepare subdivision.
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
int sumnodes = 0;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
}
return 4*sumnodes;
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
if (degree_ < 2) {
std::copy(cc, cc + dim, coord);
return;
}
// Degree 2 case.
int tetindex = index / 4;
const int subindex = index % 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
const double a = 0.138196601125010515179541316563436;
// Barycentric coordinates of our point in the tet.
double baryc[4] = { a, a, a, a };
baryc[subindex] = 1.0 - 3.0*a;
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd];
}
return;
}
THROW("Should never reach this point.");
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.cell_volumes[cell_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
int tetindex = index / 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
return 0.25*tetVolume(cc, fc, n0c, n1c);
}
THROW("Should never reach this point.");
}
private:
const UnstructuredGrid& grid_;
const int cell_;
const int degree_;
};
/// A class providing numerical quadrature for faces.
/// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by face area,
/// so weights always sum to face area.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = face area, x_0 = face centroid
/// Degree 2 method:
/// Based on subdivision of the face into triangles,
/// with the centroid as a common vertex, and the triangle
/// edge midpoint rule.
/// Triangle i consists of the centroid C, nodes N_i and N_{i+1}.
/// Its area is A_i.
/// n = 2 * nn (nn = num nodes in face)
/// For i = 0..(nn-1):
/// w_i = 1/3 A_i.
/// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i
/// x_i = (N_i + N_{i+1})/2
/// x_{nn+i} = (C + N_i)/2
/// All N and A indices are interpreted cyclic, modulus nn.
class FaceQuadrature
{
public:
FaceQuadrature(const UnstructuredGrid& grid,
const int face,
const int degree)
: grid_(grid), face_(face), degree_(degree)
{
if (grid_.dimensions != 3) {
THROW("FaceQuadrature only implemented for 3D case.");
}
if (degree_ > 2) {
THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented.");
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
if (degree_ < 2) {
std::copy(fc, fc + dim, coord);
return;
}
// Degree 2 case.
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]);
}
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node = fnodes[index - nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]);
}
}
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.face_areas[face_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
return area / 3.0;
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node0 = fnodes[(index - 1) % nn];
const int node1 = fnodes[index - nn];
const int node2 = fnodes[(index + 1) % nn];
const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc);
return (area0 + area1) / 3.0;
}
}
private:
const UnstructuredGrid& grid_;
const int face_;
const int degree_;
};
// --------------- Methods of TransportModelTracerTofDiscGal ---------------
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
/// \param[in] param Parameters for the solver.
/// The following parameters are accepted (defaults):
/// use_cvi (false) Use ECVI velocity interpolation.
/// use_limiter (false) Use a slope limiter. If true, the next three parameters are used.
/// limiter_relative_flux_threshold (1e-3) Ignore upstream fluxes below this threshold, relative to total cell flux.
/// limiter_method ("MinUpwindFace") Limiter method used. Accepted methods are:
/// MinUpwindFace Limit cell tof to >= inflow face tofs.
/// limiter_usage ("DuringComputations") Usage pattern for limiter. Accepted choices are:
/// DuringComputations Apply limiter to cells as they are computed,
/// so downstream cells' solutions may be affected
/// by limiting in upstream cells.
/// AsPostProcess Apply in dependency order, but only after
/// computing (unlimited) solution.
/// AsSimultaneousPostProcess Apply to each cell independently, using un-
/// limited solution in neighbouring cells.
TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi)
const parameter::ParameterGroup& param)
: grid_(grid),
use_cvi_(false),
use_limiter_(false),
limiter_relative_flux_threshold_(1e-3),
limiter_method_(MinUpwindAverage),
limiter_usage_(DuringComputations),
coord_(grid.dimensions),
velocity_(grid.dimensions)
{
if (use_cvi) {
velocity_interpolation_.reset(new VelocityInterpolationECVI(grid));
const int dg_degree = param.getDefault("dg_degree", 0);
const bool use_tensorial_basis = param.getDefault("use_tensorial_basis", false);
if (use_tensorial_basis) {
basis_func_.reset(new DGBasisMultilin(grid_, dg_degree));
} else {
velocity_interpolation_.reset(new VelocityInterpolationConstant(grid));
basis_func_.reset(new DGBasisBoundedTotalDegree(grid_, dg_degree));
}
use_cvi_ = param.getDefault("use_cvi", use_cvi_);
use_limiter_ = param.getDefault("use_limiter", use_limiter_);
if (use_limiter_) {
limiter_relative_flux_threshold_ = param.getDefault("limiter_relative_flux_threshold",
limiter_relative_flux_threshold_);
const std::string limiter_method_str = param.getDefault<std::string>("limiter_method", "MinUpwindAverage");
if (limiter_method_str == "MinUpwindFace") {
limiter_method_ = MinUpwindFace;
} else if (limiter_method_str == "MinUpwindAverage") {
limiter_method_ = MinUpwindAverage;
} else {
THROW("Unknown limiter method: " << limiter_method_str);
}
const std::string limiter_usage_str = param.getDefault<std::string>("limiter_usage", "DuringComputations");
if (limiter_usage_str == "DuringComputations") {
limiter_usage_ = DuringComputations;
} else if (limiter_usage_str == "AsPostProcess") {
limiter_usage_ = AsPostProcess;
} else if (limiter_usage_str == "AsSimultaneousPostProcess") {
limiter_usage_ = AsSimultaneousPostProcess;
} else {
THROW("Unknown limiter usage spec: " << limiter_usage_str);
}
}
// A note about the use_cvi_ member variable:
// In principle, we should not need it, since the choice of velocity
// interpolation is made below, but we may need to use higher order
// quadrature to exploit CVI, so we store the choice.
// An alternative would be to add a virtual method isConstant() to
// the VelocityInterpolationInterface.
if (use_cvi_) {
velocity_interpolation_.reset(new VelocityInterpolationECVI(grid_));
} else {
velocity_interpolation_.reset(new VelocityInterpolationConstant(grid_));
}
}
@@ -445,7 +130,6 @@ namespace Opm
void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff)
{
darcyflux_ = darcyflux;
@@ -455,11 +139,11 @@ namespace Opm
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
// THROW("Sources do not sum to zero: " << cum_src);
MESSAGE("Warning: sources do not sum to zero: " << cum_src);
}
#endif
degree_ = degree;
const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_);
const int num_basis = basis_func_->numBasisFunc();
tof_coeff.resize(num_basis*grid_.number_of_cells);
std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0);
tof_coeff_ = &tof_coeff[0];
@@ -471,6 +155,19 @@ namespace Opm
grad_basis_.resize(num_basis*grid_.dimensions);
velocity_interpolation_->setupFluxes(darcyflux);
reorderAndTransport(grid_, darcyflux);
switch (limiter_usage_) {
case AsPostProcess:
applyLimiterAsPostProcess();
break;
case AsSimultaneousPostProcess:
applyLimiterAsSimultaneousPostProcess();
break;
case DuringComputations:
// Do nothing.
break;
default:
THROW("Unknown limiter usage choice: " << limiter_usage_);
}
}
@@ -489,14 +186,25 @@ namespace Opm
// equal to Res = Jac*c - rhs, and we compute rhs directly.
const int dim = grid_.dimensions;
const int num_basis = DGBasis::numBasisFunc(dim, degree_);
const int num_basis = basis_func_->numBasisFunc();
std::fill(rhs_.begin(), rhs_.end(), 0.0);
std::fill(jac_.begin(), jac_.end(), 0.0);
// Compute cell residual contribution.
// Note: Assumes that \int_K b_j = 0 for all j > 0
rhs_[0] += porevolume_[cell];
{
const int deg_needed = basis_func_->degree();
CellQuadrature quad(grid_, cell, deg_needed);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// Integral of: b_i \phi
quad.quadPtCoord(quad_pt, &coord_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
rhs_[j] += w * basis_[j] * porevolume_[cell] / grid_.cell_volumes[cell];
}
}
}
// Compute upstream residual contribution.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
@@ -510,23 +218,28 @@ namespace Opm
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
continue;
}
if (flux >= 0.0) {
// This is an outflow boundary.
continue;
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds
// (where u_h^{ext} is the upstream unknown (tof)).
// Quadrature degree set to 2*D, since u_h^{ext} varies
// with degree D, and b_j too. We assume that the normal
// velocity is constant (this assumption may have to go
// for higher order than DG1).
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, degree_);
const int deg_needed = 2*basis_func_->degree();
FaceQuadrature quad(grid_, face, deg_needed);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
basis_func_->eval(upstream_cell, &coord_[0], &basis_nb_[0]);
const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
const double w = quad.quadPtWeight(quad_pt);
@@ -539,12 +252,22 @@ namespace Opm
// Compute cell jacobian contribution. We use Fortran ordering
// for jac_, i.e. rows cycling fastest.
{
CellQuadrature quad(grid_, cell, 2*degree_ - 1);
// Even with ECVI velocity interpolation, degree of precision 1
// is sufficient for optimal convergence order for DG1 when we
// use linear (total degree 1) basis functions.
// With bi(tri)-linear basis functions, it still seems sufficient
// for convergence order 2, but the solution looks much better and
// has significantly lower error with degree of precision 2.
// For now, we err on the side of caution, and use 2*degree, even
// though this is wasteful for the pure linear basis functions.
// const int deg_needed = 2*basis_func_->degree() - 1;
const int deg_needed = 2*basis_func_->degree();
CellQuadrature quad(grid_, cell, deg_needed);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// b_i (v \cdot \grad b_j)
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
basis_func_->evalGrad(cell, &coord_[0], &grad_basis_[0]);
velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
@@ -573,11 +296,11 @@ namespace Opm
// Do quadrature over the face to compute
// \int_{\partial K} b_i (v(x) \cdot n) b_j ds
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, 2*degree_);
FaceQuadrature quad(grid_, face, 2*basis_func_->degree());
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// u^ext flux B (B = {b_j})
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
@@ -598,10 +321,10 @@ namespace Opm
const double flux_density = flux / grid_.cell_volumes[cell];
// Do quadrature over the cell to compute
// \int_{K} b_i flux b_j dx
CellQuadrature quad(grid_, cell, 2*degree_);
CellQuadrature quad(grid_, cell, 2*basis_func_->degree());
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
basis_func_->eval(cell, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
@@ -619,6 +342,7 @@ namespace Opm
MAT_SIZE_T ldb = num_basis;
MAT_SIZE_T info = 0;
orig_jac_ = jac_;
orig_rhs_ = rhs_;
dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
if (info != 0) {
// Print the local matrix and rhs.
@@ -632,12 +356,41 @@ namespace Opm
}
std::cerr << "and b = \n";
for (int row = 0; row < n; ++row) {
std::cerr << " " << rhs_[row] << '\n';
std::cerr << " " << orig_rhs_[row] << '\n';
}
THROW("Lapack error: " << info << " encountered in cell " << cell);
}
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
// Apply limiter.
if (basis_func_->degree() > 0 && use_limiter_ && limiter_usage_ == DuringComputations) {
#ifdef EXTRA_VERBOSE
std::cout << "Cell: " << cell << " ";
std::cout << "v = ";
for (int dd = 0; dd < dim; ++dd) {
std::cout << velocity_[dd] << ' ';
}
std::cout << " grad tau = ";
for (int dd = 0; dd < dim; ++dd) {
std::cout << tof_coeff_[num_basis*cell + dd + 1] << ' ';
}
const double prod = std::inner_product(velocity_.begin(), velocity_.end(),
tof_coeff_ + num_basis*cell + 1, 0.0);
const double vv = std::inner_product(velocity_.begin(), velocity_.end(),
velocity_.begin(), 0.0);
const double gg = std::inner_product(tof_coeff_ + num_basis*cell + 1,
tof_coeff_ + num_basis*cell + num_basis,
tof_coeff_ + num_basis*cell + 1, 0.0);
std::cout << " prod = " << std::inner_product(velocity_.begin(), velocity_.end(),
tof_coeff_ + num_basis*cell + 1, 0.0);
std::cout << " normalized = " << prod/std::sqrt(vv*gg);
std::cout << " angle = " << std::acos(prod/std::sqrt(vv*gg))*360.0/(2.0*M_PI);
std::cout << std::endl;
#endif
applyLimiter(cell, tof_coeff_);
}
}
@@ -654,4 +407,272 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyLimiter(const int cell, double* tof)
{
switch (limiter_method_) {
case MinUpwindFace:
applyMinUpwindFaceLimiter(cell, tof);
break;
case MinUpwindAverage:
applyMinUpwindAverageLimiter(cell, tof);
break;
default:
THROW("Limiter type not implemented: " << limiter_method_);
}
}
void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof)
{
if (basis_func_->degree() != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
}
// Limiter principles:
// 1. Let M be the minimum TOF value on the upstream faces,
// evaluated in the upstream cells. Then the value at all
// points in this cell shall be at least M.
// Upstream faces whose flux does not exceed the relative
// flux threshold are not considered for this minimum.
// 2. The TOF shall not be below zero in any point.
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
double downstream_flux = 0.0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux < 0.0) {
upstream_flux += flux;
} else {
downstream_flux += flux;
}
}
// In the presence of sources, significant fluxes may be missing from the computed fluxes,
// setting the total flux to the (positive) maximum avoids this: since source is either
// inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
const double total_flux = std::max(-upstream_flux, downstream_flux);
// Find minimum tof on upstream faces and for this cell.
const int dim = grid_.dimensions;
const int num_basis = basis_func_->numBasisFunc();
double min_upstream_tof = 1e100;
double min_here_tof = 1e100;
int num_upstream_faces = 0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
int upstream_cell = -1;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
upstream_cell = grid_.face_cells[2*face+1];
} else {
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_);
if (upstream) {
++num_upstream_faces;
}
bool interior = (upstream_cell >= 0);
// Evaluate the solution in all corners.
for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
basis_func_->eval(cell, nc, &basis_[0]);
const double tof_here = std::inner_product(basis_.begin(), basis_.end(),
tof_coeff_ + num_basis*cell, 0.0);
min_here_tof = std::min(min_here_tof, tof_here);
if (upstream) {
if (interior) {
basis_func_->eval(upstream_cell, nc, &basis_nb_[0]);
const double tof_upstream
= std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
min_upstream_tof = std::min(min_upstream_tof, tof_upstream);
} else {
// Allow tof down to 0 on inflow boundaries.
min_upstream_tof = std::min(min_upstream_tof, 0.0);
}
}
}
}
// Compute slope multiplier (limiter).
if (num_upstream_faces == 0) {
min_upstream_tof = 0.0;
min_here_tof = 0.0;
}
if (min_upstream_tof < 0.0) {
min_upstream_tof = 0.0;
}
const double tof_c = tof_coeff_[num_basis*cell];
double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
if (tof_c < min_upstream_tof) {
// Handle by setting a flat solution.
std::cout << "Trouble in cell " << cell << std::endl;
limiter = 0.0;
basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell);
}
ASSERT(limiter >= 0.0);
// Actually do the limiting (if applicable).
if (limiter < 1.0) {
// std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
basis_func_->multiplyGradient(limiter, tof + num_basis*cell);
} else {
// std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
}
}
void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof)
{
if (basis_func_->degree() != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
}
// Limiter principles:
// 1. Let M be the average TOF value of the upstream cells.
/// Then the value at all points in this cell shall be at least M.
// Upstream faces whose flux does not exceed the relative
// flux threshold are not considered for this minimum.
// 2. The TOF shall not be below zero in any point.
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
double downstream_flux = 0.0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux < 0.0) {
upstream_flux += flux;
} else {
downstream_flux += flux;
}
}
// In the presence of sources, significant fluxes may be missing from the computed fluxes,
// setting the total flux to the (positive) maximum avoids this: since source is either
// inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
const double total_flux = std::max(-upstream_flux, downstream_flux);
// Find minimum tof on upstream faces and for this cell.
const int dim = grid_.dimensions;
const int num_basis = basis_func_->numBasisFunc();
double min_upstream_tof = 1e100;
double min_here_tof = 1e100;
int num_upstream_faces = 0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
int upstream_cell = -1;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
upstream_cell = grid_.face_cells[2*face+1];
} else {
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_);
if (upstream) {
++num_upstream_faces;
}
bool interior = (upstream_cell >= 0);
// Evaluate the solution in all corners.
for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
basis_func_->eval(cell, nc, &basis_[0]);
const double tof_here = std::inner_product(basis_.begin(), basis_.end(),
tof_coeff_ + num_basis*cell, 0.0);
min_here_tof = std::min(min_here_tof, tof_here);
if (upstream) {
if (interior) {
min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]);
} else {
// Allow tof down to 0 on inflow boundaries.
min_upstream_tof = std::min(min_upstream_tof, 0.0);
}
}
}
}
// Compute slope multiplier (limiter).
if (num_upstream_faces == 0) {
min_upstream_tof = 0.0;
min_here_tof = 0.0;
}
if (min_upstream_tof < 0.0) {
min_upstream_tof = 0.0;
}
const double tof_c = tof_coeff_[num_basis*cell];
double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
if (tof_c < min_upstream_tof) {
// Handle by setting a flat solution.
std::cout << "Trouble in cell " << cell << std::endl;
limiter = 0.0;
basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell);
}
ASSERT(limiter >= 0.0);
// Actually do the limiting (if applicable).
if (limiter < 1.0) {
// std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
basis_func_->multiplyGradient(limiter, tof + num_basis*cell);
} else {
// std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
}
}
void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess()
{
// Apply the limiter sequentially to all cells.
// This means that a cell's limiting behaviour may be affected by
// any limiting applied to its upstream cells.
const std::vector<int>& seq = TransportModelInterface::sequence();
const int nc = seq.size();
ASSERT(nc == grid_.number_of_cells);
for (int i = 0; i < nc; ++i) {
const int cell = seq[i];
applyLimiter(cell, tof_coeff_);
}
}
void TransportModelTracerTofDiscGal::applyLimiterAsSimultaneousPostProcess()
{
// Apply the limiter simultaneously to all cells.
// This means that each cell is limited independently from all other cells,
// we write the resulting dofs to a new array instead of writing to tof_coeff_.
// Afterwards we copy the results back to tof_coeff_.
const int num_basis = basis_func_->numBasisFunc();
std::vector<double> tof_coeffs_new(tof_coeff_, tof_coeff_ + num_basis*grid_.number_of_cells);
for (int c = 0; c < grid_.number_of_cells; ++c) {
applyLimiter(c, &tof_coeffs_new[0]);
}
std::copy(tof_coeffs_new.begin(), tof_coeffs_new.end(), tof_coeff_);
}
} // namespace Opm

48
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
View File

@@ -33,6 +33,8 @@ namespace Opm
class IncompPropertiesInterface;
class VelocityInterpolationInterface;
class DGBasisInterface;
namespace parameter { class ParameterGroup; }
/// Implements a discontinuous Galerkin solver for
/// (single-phase) time-of-flight using reordering.
@@ -48,10 +50,26 @@ namespace Opm
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_cvi If true, use corner point velocity interpolation.
/// Otherwise, use the basic constant interpolation.
/// \param[in] param Parameters for the solver.
/// The following parameters are accepted (defaults):
/// dg_degree (0) Polynomial degree of basis functions.
/// use_tensorial_basis (false) Use tensor-product basis, interpreting dg_degree as
/// bi/tri-degree not total degree.
/// use_cvi (false) Use ECVI velocity interpolation.
/// use_limiter (false) Use a slope limiter. If true, the next three parameters are used.
/// limiter_relative_flux_threshold (1e-3) Ignore upstream fluxes below this threshold, relative to total cell flux.
/// limiter_method ("MinUpwindFace") Limiter method used. Accepted methods are:
/// MinUpwindFace Limit cell tof to >= inflow face tofs.
/// limiter_usage ("DuringComputations") Usage pattern for limiter. Accepted choices are:
/// DuringComputations Apply limiter to cells as they are computed,
/// so downstream cells' solutions may be affected
/// by limiting in upstream cells.
/// AsPostProcess Apply in dependency order, but only after
/// computing (unlimited) solution.
/// AsSimultaneousPostProcess Apply to each cell independently, using un-
/// limited solution in neighbouring cells.
TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const bool use_cvi);
const parameter::ParameterGroup& param);
/// Solve for time-of-flight.
@@ -60,7 +78,6 @@ namespace Opm
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] degree Polynomial degree of DG basis functions used.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
@@ -70,7 +87,6 @@ namespace Opm
void solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff);
private:
@@ -82,15 +98,24 @@ namespace Opm
TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&);
TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&);
// Data members
const UnstructuredGrid& grid_;
boost::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
bool use_cvi_;
bool use_limiter_;
double limiter_relative_flux_threshold_;
enum LimiterMethod { MinUpwindFace, MinUpwindAverage };
LimiterMethod limiter_method_;
enum LimiterUsage { DuringComputations, AsPostProcess, AsSimultaneousPostProcess };
LimiterUsage limiter_usage_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
int degree_;
boost::shared_ptr<DGBasisInterface> basis_func_;
double* tof_coeff_;
std::vector<double> rhs_; // single-cell right-hand-side
std::vector<double> jac_; // single-cell jacobian
std::vector<double> orig_rhs_; // single-cell right-hand-side (copy)
std::vector<double> orig_jac_; // single-cell jacobian (copy)
// Below: storage for quantities needed by solveSingleCell().
std::vector<double> coord_;
@@ -98,6 +123,17 @@ namespace Opm
std::vector<double> basis_nb_;
std::vector<double> grad_basis_;
std::vector<double> velocity_;
// Private methods
// Apply some limiter, writing to array tof
// (will read data from tof_coeff_, it is ok to call
// with tof_coeff as tof argument.
void applyLimiter(const int cell, double* tof);
void applyMinUpwindFaceLimiter(const int cell, double* tof);
void applyMinUpwindAverageLimiter(const int cell, double* tof);
void applyLimiterAsPostProcess();
void applyLimiterAsSimultaneousPostProcess();
};
} // namespace Opm