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
52namespace 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
blas_lapack.hpp
Opm::MimeticIPEvaluator
Definition: MimeticIPEvaluator.hpp:86
Opm::MimeticIPEvaluator::init
void init(const int max_nf)
Initialization routine.
Definition: MimeticIPEvaluator.hpp:136
Opm::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPEvaluator.hpp:119
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::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::dim
@ dim
Definition: MimeticIPEvaluator.hpp:90
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::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator()
Default constructor.
Definition: MimeticIPEvaluator.hpp:103
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
Opm::MimeticIPEvaluator::gravityTerm
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, Vector &gterm) const
Definition: MimeticIPEvaluator.hpp:497
Opm::MimeticIPEvaluator::CellIter
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPEvaluator.hpp:93
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::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::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::gravityFlux
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPEvaluator.hpp:540
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::MimeticIPEvaluator::getInverseMatrix
void getInverseMatrix(const CellIter &c, const Scalar totmob, FullMatrix< Scalar, SP, FortranOrdering > &Binv) const
Definition: MimeticIPEvaluator.hpp:343
Opm
Definition: BlackoilFluid.hpp:32
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::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::trace
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639
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::zero
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
Opm::SharedFortranMatrix
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:590
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