Loading...
Searching...
No Matches
MimeticIPEvaluator.hpp
//===========================================================================
//
// File: MimeticIPEvaluator.hpp
//
// Created: Thu Jun 25 13:09:23 2009
//
// Author(s): B�rd Skaflestad <bard.skaflestad@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_MIMETICIPEVALUATOR_HEADER
#define OPENRS_MIMETICIPEVALUATOR_HEADER
#include <algorithm>
#include <vector>
#include <array>
#include <opm/common/ErrorMacros.hpp>
#include <opm/grid/utility/SparseTable.hpp>
#include <opm/porsol/common/fortran.hpp>
#include <opm/porsol/common/blas_lapack.hpp>
#include <opm/porsol/common/Matrix.hpp>
namespace Opm {
template <class GridInterface, class RockInterface>
class MimeticIPEvaluator
{
public:
static constexpr int dim = GridInterface::Dimension;
typedef typename GridInterface::CellIterator CellIter;
typedef typename CellIter::Scalar Scalar;
MimeticIPEvaluator()
: max_nf_(-1),
fa_ ( ),
t1_ ( ),
t2_ ( ),
Binv_ ( )
{}
explicit MimeticIPEvaluator(const int max_nf)
: max_nf_(max_nf ),
fa_ (max_nf * max_nf),
t1_ (max_nf * dim ),
t2_ (max_nf * dim ),
Binv_ ( ),
gflux_ ( )
{}
void init(const int max_nf)
{
max_nf_ = max_nf;
std::vector<double>(max_nf * max_nf).swap(fa_);
std::vector<double>(max_nf * dim ).swap(t1_);
std::vector<double>(max_nf * dim ).swap(t2_);
}
template<class Vector>
void reserveMatrices(const Vector& sz)
{
typedef typename Vector::value_type vt;
Vector sz2(sz.size());
std::transform(sz.begin(), sz.end(), sz2.begin(),
[](const vt& input) { return input*input; });
Binv_ .allocate(sz2.begin(), sz2.end());
gflux_.allocate(sz .begin(), sz .end());
}
void buildStaticContrib(const CellIter& c,
const RockInterface& r,
const typename CellIter::Vector& grav,
const int nf)
{
typedef typename CellIter::FaceIterator FI;
typedef typename CellIter::Vector CV;
typedef typename FI ::Vector FV;
const int ci = c->index();
static_assert (FV::dimension == int(dim), "");
assert (int(t1_.size()) >= nf * dim);
assert (int(t2_.size()) >= nf * dim);
assert (int(fa_.size()) >= nf * nf);
SharedFortranMatrix T1 (nf, dim, &t1_ [0]);
SharedFortranMatrix T2 (nf, dim, &t2_ [0]);
SharedFortranMatrix fa (nf, nf , &fa_ [0]);
SharedFortranMatrix Binv(nf, nf , &Binv_[ci][0]);
// Clear matrices of any residual data.
zero(Binv); zero(T1); zero(T2); zero(fa);
typename RockInterface::PermTensor K = r.permeability(ci);
const CV Kg = prod(K, grav);
// Setup: Binv <- I, T1 <- N, T2 <- C
// Complete: gflux <- N*K*g
const CV cc = c->centroid();
int i = 0;
for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
Binv(i,i) = Scalar(1.0);
fa(i,i) = f->area();
FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
FV fn = f->normal (); fn *= fa(i,i);
gflux_[ci][i] = fn * Kg;
for (int j = 0; j < dim; ++j) {
T1(i,j) = fn[j];
T2(i,j) = fc[j];
}
}
assert (i == nf);
// T2 <- orth(T2)
if (orthogonalizeColumns(T2) != 0) {
assert (false);
}
// Binv <- Binv - T2*T2' == I - Q*Q'
symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
// Binv <- diag(A) * Binv * diag(A)
symmetricUpdate(fa, Binv);
// T2 <- N*K
matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
// Binv <- (T2*N' + Binv) / vol(c)
// == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
//
// where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
//
Scalar t = Scalar(6.0) * trace(K) / dim;
matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
t / c->volume(), Binv );
}
template<class FluidInterface, class Sat>
void computeDynamicParams(const CellIter& c,
const FluidInterface& fl,
const std::vector<Sat>& s)
{
const int ci = c->index();
std::array<Scalar, FluidInterface::NumberOfPhases> mob ;
std::array<Scalar, FluidInterface::NumberOfPhases> rho ;
fl.phaseMobilities(ci, s[ci], mob);
fl.phaseDensities (ci, rho);
totmob_ = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
mob_dens_ = std::inner_product(rho.begin(), rho.end(), mob.begin(),
Scalar(0.0));
}
template<template<typename> class SP>
void getInverseMatrix(const CellIter& c,
FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
{
getInverseMatrix(c, totmob_, Binv);
}
template<template<typename> class SP>
void getInverseMatrix(const CellIter& c,
const Scalar totmob,
FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
{
const int ci = c->index();
std::transform(Binv_[ci].begin(), Binv_[ci].end(), Binv.data(),
[totmob](const Scalar& input) { return input*totmob; });
}
template<class PermTensor, template<typename> class SP>
void evaluate(const CellIter& c,
const PermTensor& K,
FullMatrix<Scalar,SP,FortranOrdering>& Binv)
{
typedef typename CellIter::FaceIterator FI;
typedef typename CellIter::Vector CV;
typedef typename FI ::Vector FV;
const int nf = Binv.numRows();
static_assert(FV::dimension == int(dim), "");
assert(Binv.numRows() <= max_nf_);
assert(Binv.numRows() == Binv.numCols());
assert(int(t1_.size()) >= nf * dim);
assert(int(t2_.size()) >= nf * dim);
assert(int(fa_.size()) >= nf * nf);
SharedFortranMatrix T1(nf, dim, &t1_[0]);
SharedFortranMatrix T2(nf, dim, &t2_[0]);
SharedFortranMatrix fa(nf, nf , &fa_[0]);
// Clear matrices of any residual data.
zero(Binv); zero(T1); zero(T2); zero(fa);
// Setup: Binv <- I, T1 <- N, T2 <- C
const CV cc = c->centroid();
int i = 0;
for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
Binv(i,i) = Scalar(1.0);
fa(i,i) = f->area();
FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
FV fn = f->normal (); fn *= fa(i,i);
for (int j = 0; j < dim; ++j) {
T1(i,j) = fn[j];
T2(i,j) = fc[j];
}
}
assert(i == nf);
// T2 <- orth(T2)
if (orthogonalizeColumns(T2) != 0) {
assert (false);
}
// Binv <- Binv - T2*T2' == I - Q*Q'
symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
// Binv <- diag(A) * Binv * diag(A)
symmetricUpdate(fa, Binv);
// T2 <- N*K
matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
// Binv <- (T2*N' + Binv) / vol(c)
// == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
//
// where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
//
Scalar t = Scalar(6.0) * trace(K) / dim;
matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
t / c->volume(), Binv );
}
template<class Vector>
void gravityTerm(const CellIter& c,
const typename CellIter::Vector& grav,
const Scalar omega,
Vector& gterm) const
{
typedef typename CellIter::FaceIterator FI;
typedef typename CellIter::Vector Point;
assert (gterm.size() <= max_nf_);
const Point cc = c->centroid();
int i = 0;
for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
Point fc = f->centroid();
fc -= cc;
gterm[i] = omega * (fc * grav);
}
}
template<class Vector>
void gravityTerm(const CellIter& c,
const typename CellIter::Vector& grav,
Vector& gterm) const
{
gravityTerm(c, grav, mob_dens_ / totmob_, gterm);
}
template<class FluidInterface, class Sat, class Vector>
void gravityTerm(const CellIter& c,
const FluidInterface& fl,
const std::vector<Sat>& s,
const typename CellIter::Vector& grav,
Vector& gterm) const
{
const int ci = c->index();
std::array<Scalar, FluidInterface::NumberOfPhases> mob;
std::array<Scalar, FluidInterface::NumberOfPhases> rho;
fl.phaseMobilities(ci, s[ci], mob);
fl.phaseDensities (ci, rho);
Scalar totmob = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
Scalar omega = std::inner_product(rho.begin(), rho.end(), mob.begin(),
Scalar(0.0)) / totmob;
gravityTerm(c, grav, omega, gterm);
}
template<class Vector>
void gravityFlux(const CellIter& c,
Vector& gflux) const
{
Scalar mob_dens = mob_dens_;
std::transform(gflux_[c->index()].begin(), gflux_[c->index()].end(),
gflux.begin(),
[mob_dens](const Scalar& input) { return input*mob_dens; });
}
private:
int max_nf_ ;
Scalar totmob_ ;
Scalar mob_dens_ ;
std::vector<Scalar> fa_, t1_, t2_;
Opm::SparseTable<Scalar> Binv_ ;
Opm::SparseTable<Scalar> gflux_ ;
};
} // namespace Opm
#endif // OPENRS_MIMETICIPEVALUATOR_HEADER
Opm::MimeticIPEvaluator::init
void init(const int max_nf)
Initialization routine.
Definition MimeticIPEvaluator.hpp:135
Opm::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator(const int max_nf)
Constructor.
Definition MimeticIPEvaluator.hpp:118
Opm::MimeticIPEvaluator::gravityTerm
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, const Scalar omega, Vector &gterm) const
Computes the mimetic discretization of the gravity term in Darcy's law.
Definition MimeticIPEvaluator.hpp:478
Opm::MimeticIPEvaluator::evaluate
void evaluate(const CellIter &c, const PermTensor &K, FullMatrix< Scalar, SP, FortranOrdering > &Binv)
Main evaluation routine.
Definition MimeticIPEvaluator.hpp:385
Opm::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator()
Default constructor.
Definition MimeticIPEvaluator.hpp:102
Opm::MimeticIPEvaluator::computeDynamicParams
void computeDynamicParams(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s)
Evaluate dynamic (saturation dependent) properties in single cell.
Definition MimeticIPEvaluator.hpp:298
Opm::MimeticIPEvaluator::CellIter
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition MimeticIPEvaluator.hpp:92
Opm::MimeticIPEvaluator::buildStaticContrib
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine.
Definition MimeticIPEvaluator.hpp:202
Opm::MimeticIPEvaluator::Scalar
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product.
Definition MimeticIPEvaluator.hpp:98
Opm::MimeticIPEvaluator::reserveMatrices
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells.
Definition MimeticIPEvaluator.hpp:158
Opm::MimeticIPEvaluator::gravityFlux
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition MimeticIPEvaluator.hpp:541
Opm::MimeticIPEvaluator::dim
static constexpr int dim
The number of space dimensions.
Definition MimeticIPEvaluator.hpp:89
Opm::MimeticIPEvaluator::getInverseMatrix
void getInverseMatrix(const CellIter &c, FullMatrix< Scalar, SP, FortranOrdering > &Binv) const
Retrieve the dynamic (mobility updated) inverse mimetic inner product matrix for specific cell.
Definition MimeticIPEvaluator.hpp:337
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::matMulAdd_NN
void matMulAdd_NN(const T &a1, const FullMatrix< T, SP1, FortranOrdering > &A, const FullMatrix< T, SP2, FortranOrdering > &B, const T &a2, FullMatrix< T, SP3, FortranOrdering > &C)
GEneral Matrix-Matrix product update of other matrix.
Definition Matrix.hpp:1136
Opm::orthogonalizeColumns
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space).
Definition Matrix.hpp:729
Opm::symmetricUpdate
void symmetricUpdate(const T &a1, const FullMatrix< T, StoragePolicy, FortranOrdering > &A, const T &a2, FullMatrix< T, StoragePolicy, FortranOrdering > &C)
Symmetric, rank update of symmetric matrix.
Definition Matrix.hpp:829
Opm::trace
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition Matrix.hpp:638
Opm::matMulAdd_NT
void matMulAdd_NT(const T &a1, const FullMatrix< T, SP1, FortranOrdering > &A, const FullMatrix< T, SP2, FortranOrdering > &B, const T &a2, FullMatrix< T, SP3, FortranOrdering > &C)
GEneral Matrix-Matrix product update of other matrix.
Definition Matrix.hpp:1194
Opm::zero
void zero(Matrix &A)
Zero-fill a.
Definition Matrix.hpp:602
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