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 static constexpr int 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 explicit 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
202 void buildStaticContrib(const CellIter& c,
203 const RockInterface& r,
204 const typename CellIter::Vector& grav,
205 const int nf)
206 {
207 typedef typename CellIter::FaceIterator FI;
208 typedef typename CellIter::Vector CV;
209 typedef typename FI ::Vector FV;
210
211 const int ci = c->index();
212
213 static_assert (FV::dimension == int(dim), "");
214 assert (int(t1_.size()) >= nf * dim);
215 assert (int(t2_.size()) >= nf * dim);
216 assert (int(fa_.size()) >= nf * nf);
217
218 SharedFortranMatrix T1 (nf, dim, &t1_ [0]);
219 SharedFortranMatrix T2 (nf, dim, &t2_ [0]);
220 SharedFortranMatrix fa (nf, nf , &fa_ [0]);
221 SharedFortranMatrix Binv(nf, nf , &Binv_[ci][0]);
222
223 // Clear matrices of any residual data.
224 zero(Binv); zero(T1); zero(T2); zero(fa);
225
226 typename RockInterface::PermTensor K = r.permeability(ci);
227 const CV Kg = prod(K, grav);
228
229 // Setup: Binv <- I, T1 <- N, T2 <- C
230 // Complete: gflux <- N*K*g
231 const CV cc = c->centroid();
232 int i = 0;
233 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
234 Binv(i,i) = Scalar(1.0);
235 fa(i,i) = f->area();
236
237 FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
238 FV fn = f->normal (); fn *= fa(i,i);
239
240 gflux_[ci][i] = fn * Kg;
241
242 for (int j = 0; j < dim; ++j) {
243 T1(i,j) = fn[j];
244 T2(i,j) = fc[j];
245 }
246 }
247 assert (i == nf);
248
249 // T2 <- orth(T2)
250 if (orthogonalizeColumns(T2) != 0) {
251 assert (false);
252 }
253
254 // Binv <- Binv - T2*T2' == I - Q*Q'
255 symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
256
257 // Binv <- diag(A) * Binv * diag(A)
258 symmetricUpdate(fa, Binv);
259
260 // T2 <- N*K
261 matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
262
263 // Binv <- (T2*N' + Binv) / vol(c)
264 // == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
265 //
266 // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
267 //
268 Scalar t = Scalar(6.0) * trace(K) / dim;
269 matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
270 t / c->volume(), Binv );
271 }
272
273
297 template<class FluidInterface, class Sat>
298 void computeDynamicParams(const CellIter& c,
299 const FluidInterface& fl,
300 const std::vector<Sat>& s)
301 {
302 const int ci = c->index();
303
304 std::array<Scalar, FluidInterface::NumberOfPhases> mob ;
305 std::array<Scalar, FluidInterface::NumberOfPhases> rho ;
306 fl.phaseMobilities(ci, s[ci], mob);
307 fl.phaseDensities (ci, rho);
308
309 totmob_ = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
310 mob_dens_ = std::inner_product(rho.begin(), rho.end(), mob.begin(),
311 Scalar(0.0));
312 }
313
314
336 template<template<typename> class SP>
337 void getInverseMatrix(const CellIter& c,
338 FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
339 {
340 getInverseMatrix(c, totmob_, Binv);
341 }
342
343 template<template<typename> class SP>
344 void getInverseMatrix(const CellIter& c,
345 const Scalar totmob,
346 FullMatrix<Scalar,SP,FortranOrdering>& Binv) const
347 {
348 const int ci = c->index();
349 std::transform(Binv_[ci].begin(), Binv_[ci].end(), Binv.data(),
350 [totmob](const Scalar& input) { return input*totmob; });
351 }
352
384 template<class PermTensor, template<typename> class SP>
385 void evaluate(const CellIter& c,
386 const PermTensor& K,
387 FullMatrix<Scalar,SP,FortranOrdering>& Binv)
388 {
389 typedef typename CellIter::FaceIterator FI;
390 typedef typename CellIter::Vector CV;
391 typedef typename FI ::Vector FV;
392
393 const int nf = Binv.numRows();
394
395 static_assert(FV::dimension == int(dim), "");
396 assert(Binv.numRows() <= max_nf_);
397 assert(Binv.numRows() == Binv.numCols());
398 assert(int(t1_.size()) >= nf * dim);
399 assert(int(t2_.size()) >= nf * dim);
400 assert(int(fa_.size()) >= nf * nf);
401
402 SharedFortranMatrix T1(nf, dim, &t1_[0]);
403 SharedFortranMatrix T2(nf, dim, &t2_[0]);
404 SharedFortranMatrix fa(nf, nf , &fa_[0]);
405
406 // Clear matrices of any residual data.
407 zero(Binv); zero(T1); zero(T2); zero(fa);
408
409 // Setup: Binv <- I, T1 <- N, T2 <- C
410 const CV cc = c->centroid();
411 int i = 0;
412 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
413 Binv(i,i) = Scalar(1.0);
414 fa(i,i) = f->area();
415
416 FV fc = f->centroid(); fc -= cc; fc *= fa(i,i);
417 FV fn = f->normal (); fn *= fa(i,i);
418
419 for (int j = 0; j < dim; ++j) {
420 T1(i,j) = fn[j];
421 T2(i,j) = fc[j];
422 }
423 }
424 assert(i == nf);
425
426 // T2 <- orth(T2)
427 if (orthogonalizeColumns(T2) != 0) {
428 assert (false);
429 }
430
431 // Binv <- Binv - T2*T2' == I - Q*Q'
432 symmetricUpdate(Scalar(-1.0), T2, Scalar(1.0), Binv);
433
434 // Binv <- diag(A) * Binv * diag(A)
435 symmetricUpdate(fa, Binv);
436
437 // T2 <- N*K
438 matMulAdd_NN(Scalar(1.0), T1, K, Scalar(0.0), T2);
439
440 // Binv <- (T2*N' + Binv) / vol(c)
441 // == (N*K*N' + t*(diag(A) * (I - Q*Q') * diag(A))) / vol(c)
442 //
443 // where t = 6/d * TRACE(K) (== 2*TRACE(K) for 3D).
444 //
445 Scalar t = Scalar(6.0) * trace(K) / dim;
446 matMulAdd_NT(Scalar(1.0) / c->volume(), T2, T1,
447 t / c->volume(), Binv );
448 }
449
450
477 template<class Vector>
478 void gravityTerm(const CellIter& c,
479 const typename CellIter::Vector& grav,
480 const Scalar omega,
481 Vector& gterm) const
482 {
483 typedef typename CellIter::FaceIterator FI;
484 typedef typename CellIter::Vector Point;
485
486 assert (gterm.size() <= max_nf_);
487
488 const Point cc = c->centroid();
489 int i = 0;
490 for (FI f = c->facebegin(); f != c->faceend(); ++f, ++i) {
491 Point fc = f->centroid();
492 fc -= cc;
493 gterm[i] = omega * (fc * grav);
494 }
495 }
496
497 template<class Vector>
498 void gravityTerm(const CellIter& c,
499 const typename CellIter::Vector& grav,
500 Vector& gterm) const
501 {
502 gravityTerm(c, grav, mob_dens_ / totmob_, gterm);
503 }
504
505 template<class FluidInterface, class Sat, class Vector>
506 void gravityTerm(const CellIter& c,
507 const FluidInterface& fl,
508 const std::vector<Sat>& s,
509 const typename CellIter::Vector& grav,
510 Vector& gterm) const
511 {
512 const int ci = c->index();
513
514 std::array<Scalar, FluidInterface::NumberOfPhases> mob;
515 std::array<Scalar, FluidInterface::NumberOfPhases> rho;
516 fl.phaseMobilities(ci, s[ci], mob);
517 fl.phaseDensities (ci, rho);
518
519 Scalar totmob = std::accumulate (mob.begin(), mob.end(), Scalar(0.0));
520 Scalar omega = std::inner_product(rho.begin(), rho.end(), mob.begin(),
521 Scalar(0.0)) / totmob;
522
523 gravityTerm(c, grav, omega, gterm);
524 }
525
526
540 template<class Vector>
541 void gravityFlux(const CellIter& c,
542 Vector& gflux) const
543 {
544 Scalar mob_dens = mob_dens_;
545 std::transform(gflux_[c->index()].begin(), gflux_[c->index()].end(),
546 gflux.begin(),
547 [mob_dens](const Scalar& input) { return input*mob_dens; });
548 }
549
550 private:
551 int max_nf_ ;
552 Scalar totmob_ ;
553 Scalar mob_dens_ ;
554 std::vector<Scalar> fa_, t1_, t2_;
555 Opm::SparseTable<Scalar> Binv_ ;
556 Opm::SparseTable<Scalar> gflux_ ;
557 };
558} // namespace Opm
559
560#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:478
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:506
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:385
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:298
Opm::MimeticIPEvaluator::gravityTerm
void gravityTerm(const CellIter &c, const typename CellIter::Vector &grav, Vector &gterm) const
Definition: MimeticIPEvaluator.hpp:498
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:202
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:541
Opm::MimeticIPEvaluator::dim
static constexpr int dim
The number of space dimensions.
Definition: MimeticIPEvaluator.hpp:89
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:337
Opm::MimeticIPEvaluator::getInverseMatrix
void getInverseMatrix(const CellIter &c, const Scalar totmob, FullMatrix< Scalar, SP, FortranOrdering > &Binv) const
Definition: MimeticIPEvaluator.hpp:344
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