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 <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
51namespace Opm {
83 template <class GridInterface, class RockInterface>
84 class MimeticIPEvaluator
85 {
86 public:
89 enum { dim = GridInterface::Dimension };
92 typedef typename GridInterface::CellIterator CellIter;
98 typedef typename CellIter::Scalar Scalar;
99
100
102 MimeticIPEvaluator()
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
200 void buildStaticContrib(const CellIter& c,
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>
296 void computeDynamicParams(const CellIter& c,
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
blas_lapack.hpp
Opm::MimeticIPEvaluator
Definition: MimeticIPEvaluator.hpp:85
Opm::MimeticIPEvaluator::init
void init(const int max_nf)
Initialization routine.
Definition: MimeticIPEvaluator.hpp:135
Opm::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator(const int max_nf)
Constructor.
Definition: MimeticIPEvaluator.hpp:118
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:476
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:504
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:383
Opm::MimeticIPEvaluator::MimeticIPEvaluator
MimeticIPEvaluator()
Default constructor.
Definition: MimeticIPEvaluator.hpp:102
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:296
Opm::MimeticIPEvaluator::gravityTerm
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, Vector &gterm) const
Definition: MimeticIPEvaluator.hpp:496
Opm::MimeticIPEvaluator::CellIter
GridInterface::CellIterator CellIter
The iterator type for iterating over grid cells.
Definition: MimeticIPEvaluator.hpp:92
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:200
Opm::MimeticIPEvaluator::dim
@ dim
Definition: MimeticIPEvaluator.hpp:89
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:98
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:158
Opm::MimeticIPEvaluator::gravityFlux
void gravityFlux(const CellIter &c, Vector &gflux) const
Compute gravity flux for all faces of single cell.
Definition: MimeticIPEvaluator.hpp:539
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:335
Opm::MimeticIPEvaluator::getInverseMatrix
void getInverseMatrix(const CellIter &c, const Scalar totmob, FullMatrix< Scalar, SP, FortranOrdering > &Binv) const
Definition: MimeticIPEvaluator.hpp:342
Opm::Elasticity::j
int j
Definition: elasticity_upscale_impl.hpp:658
Opm::Elasticity::Vector
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
Opm::Elasticity::rho
std::vector< double > rho
Definition: elasticity_upscale_impl.hpp:581
Opm
Definition: ImplicitAssembly.hpp:43
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:1136
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:729
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:829
Opm::trace
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:638
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:1194
Opm::zero
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:602
Opm::SharedFortranMatrix
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix
Definition: Matrix.hpp:589
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:668