MimeticIPEvaluator.hpp
Go to the documentation of this file.
1 //===========================================================================
2 //
3 // File: MimeticIPEvaluator.hpp
4 //
5 // Created: Thu Jun 25 13:09:23 2009
6 //
7 // Author(s): B�rd Skaflestad <bard.skaflestad@sintef.no>
8 // Atgeirr F Rasmussen <atgeirr@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 2009, 2010 Statoil ASA.
19 
20  This file is part of The Open Reservoir Simulator Project (OpenRS).
21 
22  OpenRS is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation, either version 3 of the License, or
25  (at your option) any later version.
26 
27  OpenRS is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 #ifndef OPENRS_MIMETICIPEVALUATOR_HEADER
37 #define OPENRS_MIMETICIPEVALUATOR_HEADER
38 
39 #include <algorithm>
40 #include <vector>
41 
42 #include <boost/bind.hpp>
43 #include <array>
44 
45 #include <opm/common/ErrorMacros.hpp>
46 #include <opm/core/utility/SparseTable.hpp>
47 
48 #include <opm/porsol/common/fortran.hpp>
49 #include <opm/porsol/common/blas_lapack.hpp>
50 #include <opm/porsol/common/Matrix.hpp>
51 
52 namespace Opm {
84  template <class GridInterface, class RockInterface>
85  class MimeticIPEvaluator
86  {
87  public:
90  enum { dim = GridInterface::Dimension };
93  typedef typename GridInterface::CellIterator CellIter;
99  typedef typename CellIter::Scalar Scalar;
100 
101 
103  MimeticIPEvaluator()
104  : max_nf_(-1),
105  fa_ ( ),
106  t1_ ( ),
107  t2_ ( ),
108  Binv_ ( )
109  {}
110 
111 
119  MimeticIPEvaluator(const int max_nf)
120  : max_nf_(max_nf ),
121  fa_ (max_nf * max_nf),
122  t1_ (max_nf * dim ),
123  t2_ (max_nf * dim ),
124  Binv_ ( ),
125  gflux_ ( )
126  {}
127 
128 
136  void init(const int max_nf)
137  {
138  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_);
142  }
143 
144 
158  template<class Vector>
159  void reserveMatrices(const Vector& sz)
160  {
161  typedef typename Vector::value_type vt;
162 
163  Vector sz2(sz.size());
164 
165  std::transform(sz.begin(), sz.end(), sz2.begin(),
166  boost::bind(std::multiplies<vt>(), _1, _1));
167 
168  Binv_ .allocate(sz2.begin(), sz2.end());
169  gflux_.allocate(sz .begin(), sz .end());
170  }
171 
172 
201  void buildStaticContrib(const CellIter& c,
202  const RockInterface& r,
203  const typename CellIter::Vector& grav,
204  const int nf)
205  {
206  typedef typename CellIter::FaceIterator FI;
207  typedef typename CellIter::Vector CV;
208  typedef typename FI ::Vector FV;
209 
210  const int ci = c->index();
211 
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);
216 
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]);
221 
222  // Clear matrices of any residual data.
223  zero(Binv); zero(T1); zero(T2); zero(fa);
224 
225  typename RockInterface::PermTensor K = r.permeability(ci);
226  const CV Kg = prod(K, grav);
227 
228  // Setup: Binv <- I, T1 <- N, T2 <- C
229  // Complete: gflux <- N*K*g
230  const CV cc = c->centroid();
231  int i = 0;
232  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
233  Binv(i,i) = Scalar(1.0);
234  fa(i,i) = f->area();
235 
236  FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
237  FV fn = f->normal (); fn *= fa(i,i);
238 
239  gflux_[ci][i] = fn * Kg;
240 
241  for (int j = 0; j < dim; ++j) {
242  T1(i,j) = fn[j];
243  T2(i,j) = fc[j];
244  }
245  }
246  assert (i == nf);
247 
248  // T2 <- orth(T2)
249  if (orthogonalizeColumns(T2) != 0) {
250  assert (false);
251  }
252 
253  // Binv <- Binv - T2*T2' == I - Q*Q'
254  symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
255 
256  // Binv <- diag(A) * Binv * diag(A)
257  symmetricUpdate(fa, Binv);
258 
259  // T2 <- N*K
260  matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
261 
262  // Binv <- (T2*N' + Binv) / vol(c)
263  // == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
264  //
265  // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
266  //
267  Scalar t = Scalar(6.0) * trace(K) / dim;
268  matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
269  t / c->volume(), Binv );
270  }
271 
272 
296  template<class FluidInterface, class Sat>
297  void computeDynamicParams(const CellIter& c,
298  const FluidInterface& fl,
299  const std::vector<Sat>& s)
300  {
301  const int ci = c->index();
302 
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);
307 
308  totmob_ = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
309  mob_dens_ = std::inner_product(rho.begin(), rho.end(), mob.begin(),
310  Scalar(0.0));
311  }
312 
313 
335  template<template<typename> class SP>
336  void getInverseMatrix(const CellIter& c,
337  FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
338  {
339  getInverseMatrix(c, totmob_, Binv);
340  }
341 
342  template<template<typename> class SP>
343  void getInverseMatrix(const CellIter& c,
344  const Scalar totmob,
345  FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
346  {
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));
350  }
351 
383  template<class PermTensor, template<typename> class SP>
384  void evaluate(const CellIter& c,
385  const PermTensor& K,
386  FullMatrix<Scalar,SP,FortranOrdering>& Binv)
387  {
388  typedef typename CellIter::FaceIterator FI;
389  typedef typename CellIter::Vector CV;
390  typedef typename FI ::Vector FV;
391 
392  const int nf = Binv.numRows();
393 
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);
400 
401  SharedFortranMatrix T1(nf, dim, &t1_[0]);
402  SharedFortranMatrix T2(nf, dim, &t2_[0]);
403  SharedFortranMatrix fa(nf, nf , &fa_[0]);
404 
405  // Clear matrices of any residual data.
406  zero(Binv); zero(T1); zero(T2); zero(fa);
407 
408  // Setup: Binv <- I, T1 <- N, T2 <- C
409  const CV cc = c->centroid();
410  int i = 0;
411  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
412  Binv(i,i) = Scalar(1.0);
413  fa(i,i) = f->area();
414 
415  FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
416  FV fn = f->normal (); fn *= fa(i,i);
417 
418  for (int j = 0; j < dim; ++j) {
419  T1(i,j) = fn[j];
420  T2(i,j) = fc[j];
421  }
422  }
423  assert(i == nf);
424 
425  // T2 <- orth(T2)
426  if (orthogonalizeColumns(T2) != 0) {
427  assert (false);
428  }
429 
430  // Binv <- Binv - T2*T2' == I - Q*Q'
431  symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
432 
433  // Binv <- diag(A) * Binv * diag(A)
434  symmetricUpdate(fa, Binv);
435 
436  // T2 <- N*K
437  matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
438 
439  // Binv <- (T2*N' + Binv) / vol(c)
440  // == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
441  //
442  // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
443  //
444  Scalar t = Scalar(6.0) * trace(K) / dim;
445  matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
446  t / c->volume(), Binv );
447  }
448 
449 
476  template<class Vector>
477  void gravityTerm(const CellIter& c,
478  const typename CellIter::Vector& grav,
479  const Scalar omega,
480  Vector& gterm) const
481  {
482  typedef typename CellIter::FaceIterator FI;
483  typedef typename CellIter::Vector Point;
484 
485  assert (gterm.size() <= max_nf_);
486 
487  const Point cc = c->centroid();
488  int i = 0;
489  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
490  Point fc = f->centroid();
491  fc -= cc;
492  gterm[i] = omega * (fc * grav);
493  }
494  }
495 
496  template<class Vector>
497  void gravityTerm(const CellIter& c,
498  const typename CellIter::Vector& grav,
499  Vector& gterm) const
500  {
501  gravityTerm(c, grav, mob_dens_ / totmob_, gterm);
502  }
503 
504  template<class FluidInterface, class Sat, class Vector>
505  void gravityTerm(const CellIter& c,
506  const FluidInterface& fl,
507  const std::vector<Sat>& s,
508  const typename CellIter::Vector& grav,
509  Vector& gterm) const
510  {
511  const int ci = c->index();
512 
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);
517 
518  Scalar totmob = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
519  Scalar omega = std::inner_product(rho.begin(), rho.end(), mob.begin(),
520  Scalar(0.0)) / totmob;
521 
522  gravityTerm(c, grav, omega, gterm);
523  }
524 
525 
539  template<class Vector>
540  void gravityFlux(const CellIter& c,
541  Vector& gflux) const
542  {
543  std::transform(gflux_[c->index()].begin(), gflux_[c->index()].end(),
544  gflux.begin(),
545  boost::bind(std::multiplies<Scalar>(), _1, mob_dens_));
546  }
547 
548  private:
549  int max_nf_ ;
550  Scalar totmob_ ;
551  Scalar mob_dens_ ;
552  std::vector<Scalar> fa_, t1_, t2_;
553  Opm::SparseTable<Scalar> Binv_ ;
554  Opm::SparseTable<Scalar> gflux_ ;
555  };
556 } // namespace Opm
557 
558 #endif // OPENRS_MIMETICIPEVALUATOR_HEADER
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. Specifically, .
Definition: Matrix.hpp:830
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. Specificlly, .
Definition: Matrix.hpp:1195
Opm::MimeticIPEvaluator::evaluate
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
Opm
Definition: BlackoilFluid.hpp:31
Opm::MimeticIPEvaluator::gravityTerm
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, Vector &gterm) const
Definition: MimeticIPEvaluator.hpp:497
blas_lapack.hpp
Opm::MimeticIPEvaluator::gravityTerm
void gravityTerm(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s, const typename CellIter::Vector &grav, Vector &gterm) const
Definition: MimeticIPEvaluator.hpp:505
Opm::MimeticIPEvaluator
Definition: MimeticIPEvaluator.hpp:85
Opm::MimeticIPEvaluator::Scalar
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product. Assumed to be a floating ...
Definition: MimeticIPEvaluator.hpp:99
Opm::MimeticIPEvaluator::buildStaticContrib
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
Opm::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator()
Default constructor.
Definition: MimeticIPEvaluator.hpp:103
Opm::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPEvaluator.hpp:119
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:669
Opm::SharedFortranMatrix
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:590
Opm::MimeticIPEvaluator::dim
Definition: MimeticIPEvaluator.hpp:90
Opm::MimeticIPEvaluator::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: MimeticIPEvaluator.hpp:336
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. Specificlly, .
Definition: Matrix.hpp:1137
Opm::orthogonalizeColumns
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
Opm::MimeticIPEvaluator::gravityFlux
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPEvaluator.hpp:540
Matrix.hpp
Opm::trace
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639
Opm::MimeticIPEvaluator::getInverseMatrix
void getInverseMatrix(const CellIter &c, const Scalar totmob, FullMatrix< Scalar, SP, FortranOrdering > &Binv) const
Definition: MimeticIPEvaluator.hpp:343
Opm::MimeticIPEvaluator::init
void init(const int max_nf)
Initialization routine.
Definition: MimeticIPEvaluator.hpp:136
Opm::MimeticIPEvaluator::CellIter
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPEvaluator.hpp:93
Opm::MimeticIPEvaluator::reserveMatrices
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells...
Definition: MimeticIPEvaluator.hpp:159
Opm::MimeticIPEvaluator::gravityTerm
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, const Scalar omega, Vector &gterm) const
Computes the mimetic discretization of the gravity term in Darcy's law.
Definition: MimeticIPEvaluator.hpp:477
Opm::zero
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
Opm::MimeticIPEvaluator::computeDynamicParams
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