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> 48 #include <opm/porsol/common/fortran.hpp> 49 #include <opm/porsol/common/blas_lapack.hpp> 50 #include <opm/porsol/common/Matrix.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);
217 SharedFortranMatrix T1 (nf, dim, &t1_ [0]);
218 SharedFortranMatrix T2 (nf, dim, &t2_ [0]);
219 SharedFortranMatrix fa (nf, nf , &fa_ [0]);
220 SharedFortranMatrix Binv(nf, nf , &Binv_[ci][0]);
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);
401 SharedFortranMatrix T1(nf, dim, &t1_[0]);
402 SharedFortranMatrix T2(nf, dim, &t2_[0]);
403 SharedFortranMatrix fa(nf, nf , &fa_[0]);
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_ ;
558 #endif // OPENRS_MIMETICIPEVALUATOR_HEADER MimeticIPEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPEvaluator.hpp:119
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product.
Definition: MimeticIPEvaluator.hpp:99
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine.
Definition: MimeticIPEvaluator.hpp:201
Definition: MimeticIPEvaluator.hpp:85
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
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
MimeticIPEvaluator()
Default constructor.
Definition: MimeticIPEvaluator.hpp:103
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:1137
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
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 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:830
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPEvaluator.hpp:540
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space).
Definition: Matrix.hpp:730
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
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPEvaluator.hpp:93
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
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:1195
void init(const int max_nf)
Initialization routine.
Definition: MimeticIPEvaluator.hpp:136
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 evaluate(const CellIter &c, const PermTensor &K, FullMatrix< Scalar, SP, FortranOrdering > &Binv)
Main evaluation routine.
Definition: MimeticIPEvaluator.hpp:384