My Project
MimeticIPEvaluator.hpp
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 <array>
43 
44 #include <opm/common/ErrorMacros.hpp>
45 #include <opm/grid/utility/SparseTable.hpp>
46 
47 #include <opm/porsol/common/fortran.hpp>
48 #include <opm/porsol/common/blas_lapack.hpp>
49 #include <opm/porsol/common/Matrix.hpp>
50 
51 namespace Opm {
83  template <class GridInterface, class RockInterface>
85  {
86  public:
89  enum { dim = GridInterface::Dimension };
92  typedef typename GridInterface::CellIterator CellIter;
98  typedef typename CellIter::Scalar Scalar;
99 
100 
103  : max_nf_(-1),
104  fa_ ( ),
105  t1_ ( ),
106  t2_ ( ),
107  Binv_ ( )
108  {}
109 
110 
118  MimeticIPEvaluator(const int max_nf)
119  : max_nf_(max_nf ),
120  fa_ (max_nf * max_nf),
121  t1_ (max_nf * dim ),
122  t2_ (max_nf * dim ),
123  Binv_ ( ),
124  gflux_ ( )
125  {}
126 
127 
135  void init(const int max_nf)
136  {
137  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_);
141  }
142 
143 
157  template<class Vector>
158  void reserveMatrices(const Vector& sz)
159  {
160  typedef typename Vector::value_type vt;
161 
162  Vector sz2(sz.size());
163 
164  std::transform(sz.begin(), sz.end(), sz2.begin(),
165  [](const vt& input) { return input*input; });
166 
167  Binv_ .allocate(sz2.begin(), sz2.end());
168  gflux_.allocate(sz .begin(), sz .end());
169  }
170 
171 
201  const RockInterface& r,
202  const typename CellIter::Vector& grav,
203  const int nf)
204  {
205  typedef typename CellIter::FaceIterator FI;
206  typedef typename CellIter::Vector CV;
207  typedef typename FI ::Vector FV;
208 
209  const int ci = c->index();
210 
211  static_assert (FV::dimension == int(dim), "");
212  assert (int(t1_.size()) >= nf * dim);
213  assert (int(t2_.size()) >= nf * dim);
214  assert (int(fa_.size()) >= nf * nf);
215 
216  SharedFortranMatrix T1 (nf, dim, &t1_ [0]);
217  SharedFortranMatrix T2 (nf, dim, &t2_ [0]);
218  SharedFortranMatrix fa (nf, nf , &fa_ [0]);
219  SharedFortranMatrix Binv(nf, nf , &Binv_[ci][0]);
220 
221  // Clear matrices of any residual data.
222  zero(Binv); zero(T1); zero(T2); zero(fa);
223 
224  typename RockInterface::PermTensor K = r.permeability(ci);
225  const CV Kg = prod(K, grav);
226 
227  // Setup: Binv <- I, T1 <- N, T2 <- C
228  // Complete: gflux <- N*K*g
229  const CV cc = c->centroid();
230  int i = 0;
231  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
232  Binv(i,i) = Scalar(1.0);
233  fa(i,i) = f->area();
234 
235  FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
236  FV fn = f->normal (); fn *= fa(i,i);
237 
238  gflux_[ci][i] = fn * Kg;
239 
240  for (int j = 0; j < dim; ++j) {
241  T1(i,j) = fn[j];
242  T2(i,j) = fc[j];
243  }
244  }
245  assert (i == nf);
246 
247  // T2 <- orth(T2)
248  if (orthogonalizeColumns(T2) != 0) {
249  assert (false);
250  }
251 
252  // Binv <- Binv - T2*T2' == I - Q*Q'
253  symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
254 
255  // Binv <- diag(A) * Binv * diag(A)
256  symmetricUpdate(fa, Binv);
257 
258  // T2 <- N*K
259  matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
260 
261  // Binv <- (T2*N' + Binv) / vol(c)
262  // == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
263  //
264  // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
265  //
266  Scalar t = Scalar(6.0) * trace(K) / dim;
267  matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
268  t / c->volume(), Binv );
269  }
270 
271 
295  template<class FluidInterface, class Sat>
297  const FluidInterface& fl,
298  const std::vector<Sat>& s)
299  {
300  const int ci = c->index();
301 
302  std::array<Scalar, FluidInterface::NumberOfPhases> mob ;
303  std::array<Scalar, FluidInterface::NumberOfPhases> rho ;
304  fl.phaseMobilities(ci, s[ci], mob);
305  fl.phaseDensities (ci, rho);
306 
307  totmob_ = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
308  mob_dens_ = std::inner_product(rho.begin(), rho.end(), mob.begin(),
309  Scalar(0.0));
310  }
311 
312 
334  template<template<typename> class SP>
335  void getInverseMatrix(const CellIter& c,
336  FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
337  {
338  getInverseMatrix(c, totmob_, Binv);
339  }
340 
341  template<template<typename> class SP>
342  void getInverseMatrix(const CellIter& c,
343  const Scalar totmob,
344  FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
345  {
346  const int ci = c->index();
347  std::transform(Binv_[ci].begin(), Binv_[ci].end(), Binv.data(),
348  [totmob](const Scalar& input) { return input*totmob; });
349  }
350 
382  template<class PermTensor, template<typename> class SP>
383  void evaluate(const CellIter& c,
384  const PermTensor& K,
385  FullMatrix<Scalar,SP,FortranOrdering>& Binv)
386  {
387  typedef typename CellIter::FaceIterator FI;
388  typedef typename CellIter::Vector CV;
389  typedef typename FI ::Vector FV;
390 
391  const int nf = Binv.numRows();
392 
393  static_assert(FV::dimension == int(dim), "");
394  assert(Binv.numRows() <= max_nf_);
395  assert(Binv.numRows() == Binv.numCols());
396  assert(int(t1_.size()) >= nf * dim);
397  assert(int(t2_.size()) >= nf * dim);
398  assert(int(fa_.size()) >= nf * nf);
399 
400  SharedFortranMatrix T1(nf, dim, &t1_[0]);
401  SharedFortranMatrix T2(nf, dim, &t2_[0]);
402  SharedFortranMatrix fa(nf, nf , &fa_[0]);
403 
404  // Clear matrices of any residual data.
405  zero(Binv); zero(T1); zero(T2); zero(fa);
406 
407  // Setup: Binv <- I, T1 <- N, T2 <- C
408  const CV cc = c->centroid();
409  int i = 0;
410  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
411  Binv(i,i) = Scalar(1.0);
412  fa(i,i) = f->area();
413 
414  FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
415  FV fn = f->normal (); fn *= fa(i,i);
416 
417  for (int j = 0; j < dim; ++j) {
418  T1(i,j) = fn[j];
419  T2(i,j) = fc[j];
420  }
421  }
422  assert(i == nf);
423 
424  // T2 <- orth(T2)
425  if (orthogonalizeColumns(T2) != 0) {
426  assert (false);
427  }
428 
429  // Binv <- Binv - T2*T2' == I - Q*Q'
430  symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
431 
432  // Binv <- diag(A) * Binv * diag(A)
433  symmetricUpdate(fa, Binv);
434 
435  // T2 <- N*K
436  matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
437 
438  // Binv <- (T2*N' + Binv) / vol(c)
439  // == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
440  //
441  // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
442  //
443  Scalar t = Scalar(6.0) * trace(K) / dim;
444  matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
445  t / c->volume(), Binv );
446  }
447 
448 
475  template<class Vector>
476  void gravityTerm(const CellIter& c,
477  const typename CellIter::Vector& grav,
478  const Scalar omega,
479  Vector& gterm) const
480  {
481  typedef typename CellIter::FaceIterator FI;
482  typedef typename CellIter::Vector Point;
483 
484  assert (gterm.size() <= max_nf_);
485 
486  const Point cc = c->centroid();
487  int i = 0;
488  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
489  Point fc = f->centroid();
490  fc -= cc;
491  gterm[i] = omega * (fc * grav);
492  }
493  }
494 
495  template<class Vector>
496  void gravityTerm(const CellIter& c,
497  const typename CellIter::Vector& grav,
498  Vector& gterm) const
499  {
500  gravityTerm(c, grav, mob_dens_ / totmob_, gterm);
501  }
502 
503  template<class FluidInterface, class Sat, class Vector>
504  void gravityTerm(const CellIter& c,
505  const FluidInterface& fl,
506  const std::vector<Sat>& s,
507  const typename CellIter::Vector& grav,
508  Vector& gterm) const
509  {
510  const int ci = c->index();
511 
512  std::array<Scalar, FluidInterface::NumberOfPhases> mob;
513  std::array<Scalar, FluidInterface::NumberOfPhases> rho;
514  fl.phaseMobilities(ci, s[ci], mob);
515  fl.phaseDensities (ci, rho);
516 
517  Scalar totmob = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
518  Scalar omega = std::inner_product(rho.begin(), rho.end(), mob.begin(),
519  Scalar(0.0)) / totmob;
520 
521  gravityTerm(c, grav, omega, gterm);
522  }
523 
524 
538  template<class Vector>
539  void gravityFlux(const CellIter& c,
540  Vector& gflux) const
541  {
542  Scalar mob_dens = mob_dens_;
543  std::transform(gflux_[c->index()].begin(), gflux_[c->index()].end(),
544  gflux.begin(),
545  [mob_dens](const Scalar& input) { return input*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
Definition: MimeticIPEvaluator.hpp:85
void init(const int max_nf)
Initialization routine.
Definition: MimeticIPEvaluator.hpp:135
MimeticIPEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPEvaluator.hpp:118
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:476
void evaluate(const CellIter &c, const PermTensor &K, FullMatrix< Scalar, SP, FortranOrdering > &Binv)
Main evaluation routine.
Definition: MimeticIPEvaluator.hpp:383
MimeticIPEvaluator()
Default constructor.
Definition: MimeticIPEvaluator.hpp:102
void computeDynamicParams(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s)
Evaluate dynamic (saturation dependent) properties in single cell.
Definition: MimeticIPEvaluator.hpp:296
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPEvaluator.hpp:92
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine.
Definition: MimeticIPEvaluator.hpp:200
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product.
Definition: MimeticIPEvaluator.hpp:98
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells.
Definition: MimeticIPEvaluator.hpp:158
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPEvaluator.hpp:539
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:335
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
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