36 #ifndef OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER 37 #define OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER 43 #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;
118 fa_ (max_nf * max_nf),
138 std::vector<double>(max_nf * max_nf).swap(fa_);
139 std::vector<double>(max_nf * dim ).swap(t1_);
140 std::vector<double>(max_nf * dim ).swap(t2_);
157 template<
class Vector>
160 typedef typename Vector::value_type vt;
162 Vector sz2(sz.size());
164 std::transform(sz.begin(), sz.end(), sz2.begin(),
165 boost::bind(std::multiplies<vt>(), _1, _1));
167 second_term_.allocate(sz2.begin(), sz2.end());
169 std::transform(sz.begin(), sz.end(), sz2.begin(),
170 boost::bind(std::multiplies<vt>(), _1,
int(dim)));
172 n_.allocate(sz2.begin(), sz2.end());
174 std::fill(sz2.begin(), sz2.end(), vt(dim));
175 Kg_.allocate(sz2.begin(), sz2.end());
208 const RockInterface& r,
209 const typename CellIter::Vector& grav,
217 typedef typename CellIter::FaceIterator FI;
218 typedef typename CellIter::Vector CV;
219 typedef typename FI ::Vector FV;
225 const int ci = c->index();
227 static_assert (FV::dimension ==
int(dim),
"");
228 assert (
int(t1_.size()) >= nf * dim);
229 assert (
int(t2_.size()) >= nf * dim);
230 assert (
int(fa_.size()) >= nf * nf);
232 SharedFortranMatrix T2 (nf, dim, &t2_ [0]);
233 SharedFortranMatrix fa (nf, nf , &fa_ [0]);
234 SharedFortranMatrix second_term(nf, nf, &second_term_[ci][0]);
235 SharedFortranMatrix n(nf, dim, &n_[ci][0]);
241 const CV cc = c->centroid();
243 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
244 second_term(i,i) =
Scalar(1.0);
247 FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
248 FV fn = f->normal (); fn *= fa(i,i);
250 for (
int j = 0; j < dim; ++j) {
270 Scalar(0.0), &Kg_[ci][0]);
298 template<
class Flu
idInterface,
class Sat>
300 const FluidInterface& fl,
301 const std::vector<Sat>& s)
303 const int ci = c->index();
305 std::array<Scalar, dim * dim> lambda_t;
306 std::array<Scalar, dim * dim> pmob_data;
308 SharedFortranMatrix pmob(dim, dim, &pmob_data[0]);
309 SharedFortranMatrix Kg (dim, 1 , &Kg_[ci][0]);
311 std::array<Scalar, FluidInterface::NumberOfPhases> rho;
312 fl.phaseDensities(ci, rho);
314 std::fill(dyn_Kg_.begin(), dyn_Kg_.end(),
Scalar(0.0));
315 std::fill(lambda_t.begin(), lambda_t.end(), 0.0);
317 for (
int phase = 0; phase < FluidInterface::NumberOfPhases; ++phase) {
318 fl.phaseMobility(phase, ci, s[ci], pmob);
324 std::transform(lambda_t.begin(), lambda_t.end(), pmob_data.begin(),
326 std::plus<Scalar>());
330 SharedFortranMatrix lambdaT(dim, dim, lambda_t.data());
331 SharedFortranMatrix lambdaK(dim, dim, lambdaK_.data());
332 prod(lambdaT, prock_->permeability(ci), lambdaK);
362 template<
template<
typename>
class SP>
364 FullMatrix<Scalar,SP,FortranOrdering>& Binv)
const 371 int nf = Binv.numRows();
372 ImmutableFortranMatrix n(nf, dim, &n_[ci][0]);
373 ImmutableFortranMatrix t2(nf, nf, &second_term_[ci][0]);
375 ImmutableFortranMatrix lambdaK(dim, dim, lambdaK_.data());
376 SharedFortranMatrix T2(nf, dim, &t2_[0]);
388 t / c->volume(), Binv );
405 template<
class Vector>
409 const int ci = c->index();
410 const int nf = n_.rowSize(ci) / dim;
412 ImmutableFortranMatrix N(nf, dim, &n_[ci][0]);
421 mutable std::vector<Scalar> fa_, t1_, t2_;
422 Opm::SparseTable<Scalar> second_term_ ;
423 Opm::SparseTable<Scalar> n_ ;
424 Opm::SparseTable<Scalar> Kg_ ;
425 std::array<Scalar, dim> dyn_Kg_ ;
426 std::array<double, dim*dim> lambdaK_ ;
427 const RockInterface* prock_ ;
432 #endif // OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:93
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:406
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells...
Definition: MimeticIPAnisoRelpermEvaluator.hpp:158
MimeticIPAnisoRelpermEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:116
MimeticIPAnisoRelpermEvaluator()
Default constructor.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:103
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
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 init(const int max_nf)
Initialization routine.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:135
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
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
Definition: MimeticIPAnisoRelpermEvaluator.hpp:85
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:914
void computeDynamicParams(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s)
Evaluate dynamic (saturation dependent) properties in single cell.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:299
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:363
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space).
Definition: Matrix.hpp:730
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:207
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
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:99