My Project
getQuasiImpesWeights.hpp
1 /*
2  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_GET_QUASI_IMPES_WEIGHTS_HEADER_INCLUDED
21 #define OPM_GET_QUASI_IMPES_WEIGHTS_HEADER_INCLUDED
22 
23 #include <dune/common/fvector.hh>
24 
25 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
26 
27 #include <algorithm>
28 #include <cmath>
29 
30 namespace Opm
31 {
32 
33 namespace Details
34 {
35  template <class DenseMatrix>
36  DenseMatrix transposeDenseMatrix(const DenseMatrix& M)
37  {
38  DenseMatrix tmp;
39  for (int i = 0; i < M.rows; ++i)
40  for (int j = 0; j < M.cols; ++j)
41  tmp[j][i] = M[i][j];
42 
43  return tmp;
44  }
45 } // namespace Details
46 
47 namespace Amg
48 {
49  template <class Matrix, class Vector>
50  void getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose, Vector& weights)
51  {
52  using VectorBlockType = typename Vector::block_type;
53  using MatrixBlockType = typename Matrix::block_type;
54  const Matrix& A = matrix;
55  VectorBlockType rhs(0.0);
56  rhs[pressureVarIndex] = 1.0;
57  const auto endi = A.end();
58  for (auto i = A.begin(); i != endi; ++i) {
59  const auto endj = (*i).end();
60  MatrixBlockType diag_block(0.0);
61  for (auto j = (*i).begin(); j != endj; ++j) {
62  if (i.index() == j.index()) {
63  diag_block = (*j);
64  break;
65  }
66  }
67  VectorBlockType bweights;
68  if (transpose) {
69  diag_block.solve(bweights, rhs);
70  } else {
71  auto diag_block_transpose = Details::transposeDenseMatrix(diag_block);
72  diag_block_transpose.solve(bweights, rhs);
73  }
74  double abs_max = *std::max_element(
75  bweights.begin(), bweights.end(), [](double a, double b) { return std::fabs(a) < std::fabs(b); });
76  bweights /= std::fabs(abs_max);
77  weights[i.index()] = bweights;
78  }
79  // return weights;
80  }
81 
82  template <class Matrix, class Vector>
83  Vector getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose)
84  {
85  Vector weights(matrix.N());
86  getQuasiImpesWeights(matrix, pressureVarIndex, transpose, weights);
87  return weights;
88  }
89 
90  template<class Vector, class GridView, class ElementContext, class Model>
91  void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const GridView& gridView,
92  ElementContext& elemCtx, const Model& model, std::size_t threadId)
93  {
94  using VectorBlockType = typename Vector::block_type;
95  using Matrix = typename std::decay_t<decltype(model.linearizer().jacobian())>;
96  using MatrixBlockType = typename Matrix::MatrixBlock;
97  constexpr int numEq = VectorBlockType::size();
98  using Evaluation = typename std::decay_t<decltype(model.localLinearizer(threadId).localResidual().residual(0))>
99  ::block_type;
100  VectorBlockType rhs(0.0);
101  rhs[pressureVarIndex] = 1.0;
102  int index = 0;
103  OPM_BEGIN_PARALLEL_TRY_CATCH();
104  for (const auto& elem : elements(gridView)) {
105  elemCtx.updatePrimaryStencil(elem);
106  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
107  Dune::FieldVector<Evaluation, numEq> storage;
108  model.localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0);
109  auto extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor();
110  auto scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor;
111  auto storage_scale = scvVolume / elemCtx.simulator().timeStepSize();
112  MatrixBlockType block;
113  double pressure_scale = 50e5;
114  for (int ii = 0; ii < numEq; ++ii) {
115  for (int jj = 0; jj < numEq; ++jj) {
116  block[ii][jj] = storage[ii].derivative(jj)/storage_scale;
117  if (jj == pressureVarIndex) {
118  block[ii][jj] *= pressure_scale;
119  }
120  }
121  }
122  VectorBlockType bweights;
123  MatrixBlockType block_transpose = Details::transposeDenseMatrix(block);
124  block_transpose.solve(bweights, rhs);
125  double abs_max = *std::max_element(
126  bweights.begin(), bweights.end(), [](double a, double b) { return std::fabs(a) < std::fabs(b); });
127  // probably a scaling which could give approximately total compressibility would be better
128  bweights /= std::fabs(abs_max); // given normal densities this scales weights to about 1.
129 
130  weights[index] = bweights;
131  ++index;
132  }
133  OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
134  }
135 } // namespace Amg
136 
137 } // namespace Opm
138 
139 #endif // OPM_GET_QUASI_IMPES_WEIGHTS_HEADER_INCLUDED
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27