36#ifndef OPENRS_MIMETICIPEVALUATOR_HEADER
37#define OPENRS_MIMETICIPEVALUATOR_HEADER
42#include <boost/bind.hpp>
45#include <opm/common/ErrorMacros.hpp>
46#include <opm/core/utility/SparseTable.hpp>
84 template <
class Gr
idInterface,
class RockInterface>
90 enum {
dim = GridInterface::Dimension };
93 typedef typename GridInterface::CellIterator
CellIter;
99 typedef typename CellIter::Scalar
Scalar;
121 fa_ (max_nf * max_nf),
139 std::vector<double>(max_nf * max_nf).swap(fa_);
140 std::vector<double>(max_nf *
dim ).swap(t1_);
141 std::vector<double>(max_nf *
dim ).swap(t2_);
158 template<
class Vector>
161 typedef typename Vector::value_type vt;
163 Vector sz2(sz.size());
165 std::transform(sz.begin(), sz.end(), sz2.begin(),
166 boost::bind(std::multiplies<vt>(), _1, _1));
168 Binv_ .allocate(sz2.begin(), sz2.end());
169 gflux_.allocate(sz .begin(), sz .end());
202 const RockInterface& r,
203 const typename CellIter::Vector& grav,
206 typedef typename CellIter::FaceIterator FI;
207 typedef typename CellIter::Vector CV;
208 typedef typename FI ::Vector FV;
210 const int ci = c->index();
212 static_assert (FV::dimension == int(
dim),
"");
213 assert (
int(t1_.size()) >= nf *
dim);
214 assert (
int(t2_.size()) >= nf *
dim);
215 assert (
int(fa_.size()) >= nf * nf);
225 typename RockInterface::PermTensor K = r.permeability(ci);
226 const CV Kg =
prod(K, grav);
230 const CV cc = c->centroid();
232 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
236 FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
237 FV fn = f->normal (); fn *= fa(i,i);
239 gflux_[ci][i] = fn * Kg;
241 for (
int j = 0; j <
dim; ++j) {
269 t / c->volume(), Binv );
296 template<
class Flu
idInterface,
class Sat>
298 const FluidInterface& fl,
299 const std::vector<Sat>& s)
301 const int ci = c->index();
303 std::array<Scalar, FluidInterface::NumberOfPhases> mob ;
304 std::array<Scalar, FluidInterface::NumberOfPhases> rho ;
305 fl.phaseMobilities(ci, s[ci], mob);
306 fl.phaseDensities (ci, rho);
308 totmob_ = std::accumulate (mob.begin(), mob.end(),
Scalar(0.0));
309 mob_dens_ = std::inner_product(rho.begin(), rho.end(), mob.begin(),
335 template<
template<
typename>
class SP>
337 FullMatrix<Scalar,SP,FortranOrdering>& Binv)
const
342 template<
template<
typename>
class SP>
345 FullMatrix<Scalar,SP,FortranOrdering>& Binv)
const
347 const int ci = c->index();
348 std::transform(Binv_[ci].begin(), Binv_[ci].end(), Binv.data(),
349 boost::bind(std::multiplies<Scalar>(), _1, totmob));
383 template<
class PermTensor,
template<
typename>
class SP>
386 FullMatrix<Scalar,SP,FortranOrdering>& Binv)
388 typedef typename CellIter::FaceIterator FI;
389 typedef typename CellIter::Vector CV;
390 typedef typename FI ::Vector FV;
392 const int nf = Binv.numRows();
394 static_assert(FV::dimension == int(
dim),
"");
395 assert(Binv.numRows() <= max_nf_);
396 assert(Binv.numRows() == Binv.numCols());
397 assert(
int(t1_.size()) >= nf *
dim);
398 assert(
int(t2_.size()) >= nf *
dim);
399 assert(
int(fa_.size()) >= nf * nf);
409 const CV cc = c->centroid();
411 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
415 FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
416 FV fn = f->normal (); fn *= fa(i,i);
418 for (
int j = 0; j <
dim; ++j) {
446 t / c->volume(), Binv );
476 template<
class Vector>
478 const typename CellIter::Vector& grav,
482 typedef typename CellIter::FaceIterator FI;
483 typedef typename CellIter::Vector Point;
485 assert (gterm.size() <= max_nf_);
487 const Point cc = c->centroid();
489 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
490 Point fc = f->centroid();
492 gterm[i] = omega * (fc * grav);
496 template<
class Vector>
498 const typename CellIter::Vector& grav,
504 template<
class Flu
idInterface,
class Sat,
class Vector>
506 const FluidInterface& fl,
507 const std::vector<Sat>& s,
508 const typename CellIter::Vector& grav,
511 const int ci = c->index();
513 std::array<Scalar, FluidInterface::NumberOfPhases> mob;
514 std::array<Scalar, FluidInterface::NumberOfPhases> rho;
515 fl.phaseMobilities(ci, s[ci], mob);
516 fl.phaseDensities (ci, rho);
518 Scalar totmob = std::accumulate (mob.begin(), mob.end(),
Scalar(0.0));
519 Scalar omega = std::inner_product(rho.begin(), rho.end(), mob.begin(),
539 template<
class Vector>
543 std::transform(gflux_[c->index()].begin(), gflux_[c->index()].end(),
545 boost::bind(std::multiplies<Scalar>(), _1, mob_dens_));
552 std::vector<Scalar> fa_, t1_, t2_;
553 Opm::SparseTable<Scalar> Binv_ ;
554 Opm::SparseTable<Scalar> gflux_ ;
Definition: MimeticIPEvaluator.hpp:86
void init(const int max_nf)
Initialization routine.
Definition: MimeticIPEvaluator.hpp:136
MimeticIPEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPEvaluator.hpp:119
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, const Scalar omega, Vector >erm) const
Computes the mimetic discretization of the gravity term in Darcy's law.
Definition: MimeticIPEvaluator.hpp:477
void gravityTerm(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s, const typename CellIter::Vector &grav, Vector >erm) const
Definition: MimeticIPEvaluator.hpp:505
@ dim
Definition: MimeticIPEvaluator.hpp:90
void evaluate(const CellIter &c, const PermTensor &K, FullMatrix< Scalar, SP, FortranOrdering > &Binv)
Main evaluation routine. Computes the inverse of the matrix representation of the mimetic inner produ...
Definition: MimeticIPEvaluator.hpp:384
MimeticIPEvaluator()
Default constructor.
Definition: MimeticIPEvaluator.hpp:103
void computeDynamicParams(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s)
Evaluate dynamic (saturation dependent) properties in single cell.
Definition: MimeticIPEvaluator.hpp:297
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, Vector >erm) const
Definition: MimeticIPEvaluator.hpp:497
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPEvaluator.hpp:93
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine. Computes the inverse of the matrix representation of the mimetic inner produ...
Definition: MimeticIPEvaluator.hpp:201
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product. Assumed to be a floating ...
Definition: MimeticIPEvaluator.hpp:99
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells.
Definition: MimeticIPEvaluator.hpp:159
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPEvaluator.hpp:540
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:336
void getInverseMatrix(const CellIter &c, const Scalar totmob, FullMatrix< Scalar, SP, FortranOrdering > &Binv) const
Definition: MimeticIPEvaluator.hpp:343
Definition: BlackoilFluid.hpp:32
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. Specificlly, .
Definition: Matrix.hpp:1137
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space). Based on a QR factorization of the...
Definition: Matrix.hpp:730
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. Specifically, .
Definition: Matrix.hpp:830
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639
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. Specificlly, .
Definition: Matrix.hpp:1195
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:590
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:669