20#ifndef OPM_GET_QUASI_IMPES_WEIGHTS_HEADER_INCLUDED
21#define OPM_GET_QUASI_IMPES_WEIGHTS_HEADER_INCLUDED
23#include <dune/common/fvector.hh>
25#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
35 template <
class DenseMatrix>
36 DenseMatrix transposeDenseMatrix(
const DenseMatrix& M)
39 for (
int i = 0; i < M.rows; ++i)
40 for (
int j = 0; j < M.cols; ++j)
49 template <
class Matrix,
class Vector>
50 void getQuasiImpesWeights(
const Matrix& matrix,
const int pressureVarIndex,
const bool transpose, Vector& weights)
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()) {
67 VectorBlockType bweights;
69 diag_block.solve(bweights, rhs);
71 auto diag_block_transpose = Details::transposeDenseMatrix(diag_block);
72 diag_block_transpose.solve(bweights, rhs);
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;
82 template <
class Matrix,
class Vector>
83 Vector getQuasiImpesWeights(
const Matrix& matrix,
const int pressureVarIndex,
const bool transpose)
85 Vector weights(matrix.N());
86 getQuasiImpesWeights(matrix, pressureVarIndex, transpose, weights);
90 template<
class Vector,
class Gr
idView,
class ElementContext,
class Model>
91 void getTrueImpesWeights(
int pressureVarIndex, Vector& weights,
const GridView& gridView,
92 ElementContext& elemCtx,
const Model& model, std::size_t threadId)
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))>
100 VectorBlockType rhs(0.0);
101 rhs[pressureVarIndex] = 1.0;
103 OPM_BEGIN_PARALLEL_TRY_CATCH();
104 for (
const auto& elem : elements(gridView)) {
105 elemCtx.updatePrimaryStencil(elem);
106 elemCtx.updatePrimaryIntensiveQuantities(0);
107 Dune::FieldVector<Evaluation, numEq> storage;
108 model.localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,0, 0);
109 auto extrusionFactor = elemCtx.intensiveQuantities(0, 0).extrusionFactor();
110 auto scvVolume = elemCtx.stencil(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;
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); });
128 bweights /= std::fabs(abs_max);
130 weights[index] = bweights;
133 OPM_END_PARALLEL_TRY_CATCH(
"getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27