My Project
MimeticIPAnisoRelpermEvaluator.hpp
1 //===========================================================================
2 //
3 // File: MimeticIPAnisoRelpermEvaluator.hpp
4 //
5 // Created: Mon Oct 19 10:22:22 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // B�rd Skaflestad <bard.skaflestad@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_MIMETICIPANISORELPERMEVALUATOR_HEADER
37 #define OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER
38 
39 
40 #include <algorithm>
41 #include <vector>
42 
43 #include <opm/common/ErrorMacros.hpp>
44 #include <opm/grid/utility/SparseTable.hpp>
45 
46 #include <opm/porsol/common/fortran.hpp>
47 #include <opm/porsol/common/blas_lapack.hpp>
48 #include <opm/porsol/common/Matrix.hpp>
49 
50 namespace Opm {
82  template<class GridInterface, class RockInterface>
84  {
85  public:
88  enum { dim = GridInterface::Dimension };
91  typedef typename GridInterface::CellIterator CellIter;
97  typedef typename CellIter::Scalar Scalar;
98 
99 
102  : max_nf_(-1),
103  prock_(0)
104  {}
105 
106 
115  : max_nf_ (max_nf ),
116  fa_ (max_nf * max_nf),
117  t1_ (max_nf * dim ),
118  t2_ (max_nf * dim ),
119  second_term_ ( ),
120  n_ ( ),
121  Kg_ ( ),
122  prock_ ( 0 )
123  {}
124 
125 
133  void init(const int max_nf)
134  {
135  max_nf_ = max_nf;
136  std::vector<double>(max_nf * max_nf).swap(fa_);
137  std::vector<double>(max_nf * dim ).swap(t1_);
138  std::vector<double>(max_nf * dim ).swap(t2_);
139  }
140 
141 
155  template<class Vector>
156  void reserveMatrices(const Vector& sz)
157  {
158  typedef typename Vector::value_type vt;
159 
160  Vector sz2(sz.size());
161 
162  std::transform(sz.begin(), sz.end(), sz2.begin(),
163  [](const vt& input) { return input*input; });
164 
165  second_term_.allocate(sz2.begin(), sz2.end());
166 
167  int idim = int(dim);
168  std::transform(sz.begin(), sz.end(), sz2.begin(),
169  [idim](const vt& input) { return input*idim; });
170 
171  n_.allocate(sz2.begin(), sz2.end());
172 
173  std::fill(sz2.begin(), sz2.end(), vt(dim));
174  Kg_.allocate(sz2.begin(), sz2.end());
175  }
176 
177 
207  const RockInterface& r,
208  const typename CellIter::Vector& grav,
209  const int nf)
210  {
211  // Binv = (N*lambda*K*N' + t*diag(A)*(I - Q*Q')*diag(A))/vol
212  // ^ ^^^^^^^^^^^^^^^^^^^^^^^^^^
213  // precompute: n_ precompute: second_term_
214  // t = 6/dim * trace(lambda*K)
215 
216  typedef typename CellIter::FaceIterator FI;
217  typedef typename CellIter::Vector CV;
218  typedef typename FI ::Vector FV;
219 
220  // Now we need to remember the rocks, since we will need
221  // the permeability for dynamic assembly.
222  prock_ = &r;
223 
224  const int ci = c->index();
225 
226  static_assert (FV::dimension == int(dim), "");
227  assert (int(t1_.size()) >= nf * dim);
228  assert (int(t2_.size()) >= nf * dim);
229  assert (int(fa_.size()) >= nf * nf);
230 
231  SharedFortranMatrix T2 (nf, dim, &t2_ [0]);
232  SharedFortranMatrix fa (nf, nf , &fa_ [0]);
233  SharedFortranMatrix second_term(nf, nf, &second_term_[ci][0]);
234  SharedFortranMatrix n(nf, dim, &n_[ci][0]);
235 
236  // Clear matrices of any residual data.
237  zero(second_term); zero(n); zero(T2); zero(fa);
238 
239  // Setup: second_term <- I, n <- N, T2 <- C
240  const CV cc = c->centroid();
241  int i = 0;
242  for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
243  second_term(i,i) = Scalar(1.0);
244  fa(i,i) = f->area();
245 
246  FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
247  FV fn = f->normal (); fn *= fa(i,i);
248 
249  for (int j = 0; j < dim; ++j) {
250  n (i,j) = fn[j];
251  T2(i,j) = fc[j];
252  }
253  }
254  assert (i == nf);
255 
256  // T2 <- orth(T2)
257  if (orthogonalizeColumns(T2) != 0) {
258  assert (false);
259  }
260 
261  // second_term <- second_term - T2*T2' == I - Q*Q'
262  symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), second_term);
263 
264  // second_term <- diag(A) * second_term * diag(A)
265  symmetricUpdate(fa, second_term);
266 
267  // Gravity term: Kg_ = K * grav
268  vecMulAdd_N(Scalar(1.0), r.permeability(ci), &grav[0],
269  Scalar(0.0), &Kg_[ci][0]);
270  }
271 
272 
297  template<class FluidInterface, class Sat>
299  const FluidInterface& fl,
300  const std::vector<Sat>& s)
301  {
302  const int ci = c->index();
303 
304  std::array<Scalar, dim * dim> lambda_t;
305  std::array<Scalar, dim * dim> pmob_data;
306 
307  SharedFortranMatrix pmob(dim, dim, &pmob_data[0]);
308  SharedFortranMatrix Kg (dim, 1 , &Kg_[ci][0]);
309 
310  std::array<Scalar, FluidInterface::NumberOfPhases> rho;
311  fl.phaseDensities(ci, rho);
312 
313  std::fill(dyn_Kg_.begin(), dyn_Kg_.end(), Scalar(0.0));
314  std::fill(lambda_t.begin(), lambda_t.end(), 0.0);
315 
316  for (int phase = 0; phase < FluidInterface::NumberOfPhases; ++phase) {
317  fl.phaseMobility(phase, ci, s[ci], pmob);
318 
319  // dyn_Kg_ += (\rho_phase \lambda_phase) Kg
320  vecMulAdd_N(rho[phase], pmob, Kg.data(), Scalar(1.0), dyn_Kg_.data());
321 
322  // \lambda_t += \lambda_phase
323  std::transform(lambda_t.begin(), lambda_t.end(), pmob_data.begin(),
324  lambda_t.begin(),
325  std::plus<Scalar>());
326  }
327 
328  // lambdaK_ = (\sum_i \lambda_i) K
329  SharedFortranMatrix lambdaT(dim, dim, lambda_t.data());
330  SharedFortranMatrix lambdaK(dim, dim, lambdaK_.data());
331  prod(lambdaT, prock_->permeability(ci), lambdaK);
332  }
333 
334 
361  template<template<typename> class SP>
362  void getInverseMatrix(const CellIter& c,
363  FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
364  {
365  // Binv = (N*lambda*K*N' + t*diag(A)*(I - Q*Q')*diag(A))/vol
366  // ^ ^^^^^^^^^^^^^^^^^^^^^^^^^^
367  // precomputed: n_ precomputed: second_term_
368  // t = 6/dim * trace(lambda*K)
369  int ci = c->index();
370  int nf = Binv.numRows();
371  ImmutableFortranMatrix n(nf, dim, &n_[ci][0]);
372  ImmutableFortranMatrix t2(nf, nf, &second_term_[ci][0]);
373  Binv = t2;
374  ImmutableFortranMatrix lambdaK(dim, dim, lambdaK_.data());
375  SharedFortranMatrix T2(nf, dim, &t2_[0]);
376 
377  // T2 <- N*lambda*K
378  matMulAdd_NN(Scalar(1.0), n, lambdaK, Scalar(0.0), T2);
379 
380  // Binv <- (T2*N' + t*Binv) / vol(c)
381  // == (N*lambda*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
382  //
383  // where t = 6/d * TRACE(lambda*K) (== 2*TRACE(lambda*K) for 3D).
384  //
385  Scalar t = Scalar(6.0) * trace(lambdaK) / dim;
386  matMulAdd_NT(Scalar(1.0) / c->volume(), T2, n,
387  t / c->volume(), Binv );
388  }
389 
390 
404  template<class Vector>
405  void gravityFlux(const CellIter& c,
406  Vector& gflux) const
407  {
408  const int ci = c->index();
409  const int nf = n_.rowSize(ci) / dim;
410 
411  ImmutableFortranMatrix N(nf, dim, &n_[ci][0]);
412 
413  // gflux = N (\sum_i \rho_i \lambda_i) Kg
414  vecMulAdd_N(Scalar(1.0), N, &dyn_Kg_[0],
415  Scalar(0.0), &gflux[0]);
416  }
417 
418  private:
419  int max_nf_ ;
420  mutable std::vector<Scalar> fa_, t1_, t2_;
421  Opm::SparseTable<Scalar> second_term_ ;
422  Opm::SparseTable<Scalar> n_ ;
423  Opm::SparseTable<Scalar> Kg_ ;
424  std::array<Scalar, dim> dyn_Kg_ ;
425  std::array<double, dim*dim> lambdaK_ ;
426  const RockInterface* prock_ ;
427  };
428 } // namespace Opm
429 
430 
431 #endif // OPENRS_MIMETICIPANISORELPERMEVALUATOR_HEADER
Definition: MimeticIPAnisoRelpermEvaluator.hpp:84
void init(const int max_nf)
Initialization routine.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:133
CellIter::Scalar Scalar
The element type of the matrix representation of the mimetic inner product.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:97
MimeticIPAnisoRelpermEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:114
MimeticIPAnisoRelpermEvaluator()
Default constructor.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:101
void computeDynamicParams(const CellIter &c, const FluidInterface &fl, const std::vector< Sat > &s)
Evaluate dynamic (saturation dependent) properties in single cell.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:298
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:91
void buildStaticContrib(const CellIter &c, const RockInterface &r, const typename CellIter::Vector &grav, const int nf)
Main evaluation routine.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:206
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:405
void reserveMatrices(const Vector &sz)
Reserve internal space for storing values of (static) IP contributions for given set of cells.
Definition: MimeticIPAnisoRelpermEvaluator.hpp:156
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:362
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
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:913
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