36#ifndef OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER
37#define OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER
43#include <opm/common/ErrorMacros.hpp>
44#include <opm/grid/utility/SparseTable.hpp>
46#include <opm/porsol/common/fortran.hpp>
47#include <opm/porsol/common/blas_lapack.hpp>
48#include <opm/porsol/common/Matrix.hpp>
82 template<
class Gr
idInterface,
class RockInterface>
88 enum { dim = GridInterface::Dimension };
91 typedef typename GridInterface::CellIterator
CellIter;
97 typedef typename CellIter::Scalar
Scalar;
116 fa_ (max_nf * max_nf),
136 std::vector<double>(max_nf * max_nf).swap(fa_);
137 std::vector<double>(max_nf * dim ).swap(t1_);
138 std::vector<double>(max_nf * dim ).swap(t2_);
155 template<
class Vector>
158 typedef typename Vector::value_type vt;
160 Vector sz2(sz.size());
162 std::transform(sz.begin(), sz.end(), sz2.begin(),
163 [](
const vt& input) { return input*input; });
165 second_term_.allocate(sz2.begin(), sz2.end());
168 std::transform(sz.begin(), sz.end(), sz2.begin(),
169 [idim](
const vt& input) { return input*idim; });
171 n_.allocate(sz2.begin(), sz2.end());
173 std::fill(sz2.begin(), sz2.end(), vt(dim));
174 Kg_.allocate(sz2.begin(), sz2.end());
207 const RockInterface& r,
208 const typename CellIter::Vector& grav,
216 typedef typename CellIter::FaceIterator FI;
217 typedef typename CellIter::Vector CV;
218 typedef typename FI ::Vector FV;
224 const int ci = c->index();
226 static_assert (FV::dimension == int(dim),
"");
227 assert (
int(t1_.size()) >= nf * dim);
228 assert (
int(t2_.size()) >= nf * dim);
229 assert (
int(fa_.size()) >= nf * nf);
231 SharedFortranMatrix T2 (nf, dim, &t2_ [0]);
232 SharedFortranMatrix fa (nf, nf , &fa_ [0]);
233 SharedFortranMatrix second_term(nf, nf, &second_term_[ci][0]);
234 SharedFortranMatrix n(nf, dim, &n_[ci][0]);
240 const CV cc = c->centroid();
242 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
243 second_term(i,i) =
Scalar(1.0);
246 FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
247 FV fn = f->normal (); fn *= fa(i,i);
249 for (
int j = 0; j < dim; ++j) {
269 Scalar(0.0), &Kg_[ci][0]);
297 template<
class Flu
idInterface,
class Sat>
299 const FluidInterface& fl,
300 const std::vector<Sat>& s)
302 const int ci = c->index();
304 std::array<Scalar, dim * dim> lambda_t;
305 std::array<Scalar, dim * dim> pmob_data;
307 SharedFortranMatrix pmob(dim, dim, &pmob_data[0]);
308 SharedFortranMatrix Kg (dim, 1 , &Kg_[ci][0]);
310 std::array<Scalar, FluidInterface::NumberOfPhases> rho;
311 fl.phaseDensities(ci, rho);
313 std::fill(dyn_Kg_.begin(), dyn_Kg_.end(),
Scalar(0.0));
314 std::fill(lambda_t.begin(), lambda_t.end(), 0.0);
316 for (
int phase = 0; phase < FluidInterface::NumberOfPhases; ++phase) {
317 fl.phaseMobility(phase, ci, s[ci], pmob);
323 std::transform(lambda_t.begin(), lambda_t.end(), pmob_data.begin(),
325 std::plus<Scalar>());
329 SharedFortranMatrix lambdaT(dim, dim, lambda_t.data());
330 SharedFortranMatrix lambdaK(dim, dim, lambdaK_.data());
331 prod(lambdaT, prock_->permeability(ci), lambdaK);
361 template<
template<
typename>
class SP>
363 FullMatrix<Scalar,SP,FortranOrdering>& Binv)
const
370 int nf = Binv.numRows();
371 ImmutableFortranMatrix n(nf, dim, &n_[ci][0]);
372 ImmutableFortranMatrix t2(nf, nf, &second_term_[ci][0]);
374 ImmutableFortranMatrix lambdaK(dim, dim, lambdaK_.data());
375 SharedFortranMatrix T2(nf, dim, &t2_[0]);
387 t / c->volume(), Binv );
404 template<
class Vector>
408 const int ci = c->index();
409 const int nf = n_.rowSize(ci) / dim;
411 ImmutableFortranMatrix N(nf, dim, &n_[ci][0]);
420 mutable std::vector<Scalar> fa_, t1_, t2_;
421 Opm::SparseTable<Scalar> second_term_ ;
422 Opm::SparseTable<Scalar> n_ ;
423 Opm::SparseTable<Scalar> Kg_ ;
424 std::array<Scalar, dim> dyn_Kg_ ;
425 std::array<double, dim*dim> lambdaK_ ;
426 const RockInterface* prock_ ;
Definition MimeticIPAnisoRelpermEvaluator.hpp:84
void init(const int max_nf)
Initialization routine.
Definition MimeticIPAnisoRelpermEvaluator.hpp:133
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product.
Definition MimeticIPAnisoRelpermEvaluator.hpp:97
MimeticIPAnisoRelpermEvaluator(const int max_nf)
Constructor.
Definition MimeticIPAnisoRelpermEvaluator.hpp:114
MimeticIPAnisoRelpermEvaluator()
Default constructor.
Definition MimeticIPAnisoRelpermEvaluator.hpp:101
void computeDynamicParams(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s)
Evaluate dynamic (saturation dependent) properties in single cell.
Definition MimeticIPAnisoRelpermEvaluator.hpp:298
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition MimeticIPAnisoRelpermEvaluator.hpp:91
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine.
Definition MimeticIPAnisoRelpermEvaluator.hpp:206
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition MimeticIPAnisoRelpermEvaluator.hpp:405
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells.
Definition MimeticIPAnisoRelpermEvaluator.hpp:156
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:362
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
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
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space).
Definition Matrix.hpp:729
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
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition Matrix.hpp:638
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
void zero(Matrix &A)
Zero-fill a.
Definition Matrix.hpp:602
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
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