Loading...
Searching...
No Matches
MimeticIPAnisoRelpermEvaluator.hpp
//===========================================================================
//
// File: MimeticIPAnisoRelpermEvaluator.hpp
//
// Created: Mon Oct 19 10:22:22 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B�rd Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER
#define OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER
#include <algorithm>
#include <vector>
#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 MimeticIPAnisoRelpermEvaluator
{
public:
static constexpr int dim = GridInterface::Dimension;
typedef typename GridInterface::CellIterator CellIter;
typedef typename CellIter::Scalar Scalar;
MimeticIPAnisoRelpermEvaluator()
: max_nf_(-1),
prock_(0)
{}
explicit MimeticIPAnisoRelpermEvaluator(const int max_nf)
: max_nf_ (max_nf ),
fa_ (max_nf * max_nf),
t1_ (max_nf * dim ),
t2_ (max_nf * dim ),
second_term_ ( ),
n_ ( ),
Kg_ ( ),
prock_ ( 0 )
{}
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; });
second_term_.allocate(sz2.begin(), sz2.end());
int idim = int(dim);
std::transform(sz.begin(), sz.end(), sz2.begin(),
[idim](const vt& input) { return input*idim; });
n_.allocate(sz2.begin(), sz2.end());
std::fill(sz2.begin(), sz2.end(), vt(dim));
Kg_.allocate(sz2.begin(), sz2.end());
}
void buildStaticContrib(const CellIter& c,
const RockInterface& r,
const typename CellIter::Vector& grav,
const int nf)
{
// Binv = (N*lambda*K*N' + t*diag(A)*(I - Q*Q')*diag(A))/vol
// ^ ^^^^^^^^^^^^^^^^^^^^^^^^^^
// precompute: n_ precompute: second_term_
// t = 6/dim * trace(lambda*K)
typedef typename CellIter::FaceIterator FI;
typedef typename CellIter::Vector CV;
typedef typename FI ::Vector FV;
// Now we need to remember the rocks, since we will need
// the permeability for dynamic assembly.
prock_ = &r;
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 T2 (nf, dim, &t2_ [0]);
SharedFortranMatrix fa (nf, nf , &fa_ [0]);
SharedFortranMatrix second_term(nf, nf, &second_term_[ci][0]);
SharedFortranMatrix n(nf, dim, &n_[ci][0]);
// Clear matrices of any residual data.
zero(second_term); zero(n); zero(T2); zero(fa);
// Setup: second_term <- I, n <- N, T2 <- C
const CV cc = c->centroid();
int i = 0;
for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
second_term(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) {
n (i,j) = fn[j];
T2(i,j) = fc[j];
}
}
assert (i == nf);
// T2 <- orth(T2)
if (orthogonalizeColumns(T2) != 0) {
assert (false);
}
// second_term <- second_term - T2*T2' == I - Q*Q'
symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), second_term);
// second_term <- diag(A) * second_term * diag(A)
symmetricUpdate(fa, second_term);
// Gravity term: Kg_ = K * grav
vecMulAdd_N(Scalar(1.0), r.permeability(ci), &grav[0],
Scalar(0.0), &Kg_[ci][0]);
}
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, dim * dim> lambda_t;
std::array<Scalar, dim * dim> pmob_data;
SharedFortranMatrix pmob(dim, dim, &pmob_data[0]);
SharedFortranMatrix Kg (dim, 1 , &Kg_[ci][0]);
std::array<Scalar, FluidInterface::NumberOfPhases> rho;
fl.phaseDensities(ci, rho);
std::fill(dyn_Kg_.begin(), dyn_Kg_.end(), Scalar(0.0));
std::fill(lambda_t.begin(), lambda_t.end(), 0.0);
for (int phase = 0; phase < FluidInterface::NumberOfPhases; ++phase) {
fl.phaseMobility(phase, ci, s[ci], pmob);
// dyn_Kg_ += (\rho_phase \lambda_phase) Kg
vecMulAdd_N(rho[phase], pmob, Kg.data(), Scalar(1.0), dyn_Kg_.data());
// \lambda_t += \lambda_phase
std::transform(lambda_t.begin(), lambda_t.end(), pmob_data.begin(),
lambda_t.begin(),
std::plus<Scalar>());
}
// lambdaK_ = (\sum_i \lambda_i) K
SharedFortranMatrix lambdaT(dim, dim, lambda_t.data());
SharedFortranMatrix lambdaK(dim, dim, lambdaK_.data());
prod(lambdaT, prock_->permeability(ci), lambdaK);
}
template<template<typename> class SP>
void getInverseMatrix(const CellIter& c,
FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
{
// Binv = (N*lambda*K*N' + t*diag(A)*(I - Q*Q')*diag(A))/vol
// ^ ^^^^^^^^^^^^^^^^^^^^^^^^^^
// precomputed: n_ precomputed: second_term_
// t = 6/dim * trace(lambda*K)
int ci = c->index();
int nf = Binv.numRows();
ImmutableFortranMatrix n(nf, dim, &n_[ci][0]);
ImmutableFortranMatrix t2(nf, nf, &second_term_[ci][0]);
Binv = t2;
ImmutableFortranMatrix lambdaK(dim, dim, lambdaK_.data());
SharedFortranMatrix T2(nf, dim, &t2_[0]);
// T2 <- N*lambda*K
matMulAdd_NN(Scalar(1.0), n, lambdaK, Scalar(0.0), T2);
// Binv <- (T2*N' + t*Binv) / vol(c)
// == (N*lambda*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
//
// where t = 6/d * TRACE(lambda*K) (== 2*TRACE(lambda*K) for 3D).
//
Scalar t = Scalar(6.0) * trace(lambdaK) / dim;
matMulAdd_NT(Scalar(1.0) / c->volume(), T2, n,
t / c->volume(), Binv );
}
template<class Vector>
void gravityFlux(const CellIter& c,
Vector& gflux) const
{
const int ci = c->index();
const int nf = n_.rowSize(ci) / dim;
ImmutableFortranMatrix N(nf, dim, &n_[ci][0]);
// gflux = N (\sum_i \rho_i \lambda_i) Kg
vecMulAdd_N(Scalar(1.0), N, &dyn_Kg_[0],
Scalar(0.0), &gflux[0]);
}
private:
int max_nf_ ;
mutable std::vector<Scalar> fa_, t1_, t2_;
Opm::SparseTable<Scalar> second_term_ ;
Opm::SparseTable<Scalar> n_ ;
Opm::SparseTable<Scalar> Kg_ ;
std::array<Scalar, dim> dyn_Kg_ ;
std::array<double, dim*dim> lambdaK_ ;
const RockInterface* prock_ ;
};
} // namespace Opm
#endif // OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER
Opm::MimeticIPAnisoRelpermEvaluator::init
void init(const int max_nf)
Initialization routine.
Definition MimeticIPAnisoRelpermEvaluator.hpp:133
Opm::MimeticIPAnisoRelpermEvaluator::Scalar
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product.
Definition MimeticIPAnisoRelpermEvaluator.hpp:97
Opm::MimeticIPAnisoRelpermEvaluator::MimeticIPAnisoRelpermEvaluator
MimeticIPAnisoRelpermEvaluator(const int max_nf)
Constructor.
Definition MimeticIPAnisoRelpermEvaluator.hpp:114
Opm::MimeticIPAnisoRelpermEvaluator::dim
static constexpr int dim
The number of space dimensions.
Definition MimeticIPAnisoRelpermEvaluator.hpp:88
Opm::MimeticIPAnisoRelpermEvaluator::MimeticIPAnisoRelpermEvaluator
MimeticIPAnisoRelpermEvaluator()
Default constructor.
Definition MimeticIPAnisoRelpermEvaluator.hpp:101
Opm::MimeticIPAnisoRelpermEvaluator::computeDynamicParams
void computeDynamicParams(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s)
Evaluate dynamic (saturation dependent) properties in single cell.
Definition MimeticIPAnisoRelpermEvaluator.hpp:300
Opm::MimeticIPAnisoRelpermEvaluator::CellIter
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition MimeticIPAnisoRelpermEvaluator.hpp:91
Opm::MimeticIPAnisoRelpermEvaluator::buildStaticContrib
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine.
Definition MimeticIPAnisoRelpermEvaluator.hpp:208
Opm::MimeticIPAnisoRelpermEvaluator::gravityFlux
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition MimeticIPAnisoRelpermEvaluator.hpp:407
Opm::MimeticIPAnisoRelpermEvaluator::reserveMatrices
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells.
Definition MimeticIPAnisoRelpermEvaluator.hpp:156
Opm::MimeticIPAnisoRelpermEvaluator::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 MimeticIPAnisoRelpermEvaluator.hpp:364
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::vecMulAdd_N
void vecMulAdd_N(const T &a1, const FullMatrix< T, SP, FortranOrdering > &A, const std::vector< T > &x, const T &a2, std::vector< T > &y)
GEneral Matrix-Vector product (GAXPY operation).
Definition Matrix.hpp:913
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