20#ifndef OPM_GET_QUASI_IMPES_WEIGHTS_HEADER_INCLUDED
21#define OPM_GET_QUASI_IMPES_WEIGHTS_HEADER_INCLUDED
23#include <dune/common/fvector.hh>
26#include <opm/material/common/MathToolbox.hpp>
35 template <
class DenseMatrix>
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);
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>
85 Vector weights(matrix.N());
90 template<
class Vector,
class Gr
idView,
class ElementContext,
class Model>
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;
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;
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;
136 template <
class Vector,
class Gr
idView,
class ElementContext,
class Model>
139 const GridView& gridView,
140 ElementContext& elemCtx,
142 std::size_t threadId)
152 using FluidSystem =
typename Model::FluidSystem;
153 using LhsEval = double;
154 using Indices =
typename Model::Indices;
156 using PrimaryVariables =
typename Model::PrimaryVariables;
157 using VectorBlockType =
typename Vector::block_type;
159 typename std::decay_t<
decltype(model.localLinearizer(threadId).localResidual().residual(0))>::block_type;
160 using Toolbox = MathToolbox<Evaluation>;
161 VectorBlockType rhs(0.0);
162 const auto& solution = model.solution( 0);
164 for (
const auto& elem : elements(gridView)) {
165 elemCtx.updatePrimaryStencil(elem);
166 elemCtx.updatePrimaryIntensiveQuantities(0);
167 const auto& index = elemCtx.globalSpaceIndex(0, 0);
168 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
169 const auto& fs = intQuants.fluidState();
170 VectorBlockType bweights;
172 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
173 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
174 FluidSystem::solventComponentIndex(FluidSystem::waterPhaseIdx));
175 bweights[activeCompIdx]
176 = Toolbox::template decay<LhsEval>(1 / fs.invB(FluidSystem::waterPhaseIdx));
179 double denominator = 1.0;
180 double rs = Toolbox::template decay<double>(fs.Rs());
181 double rv = Toolbox::template decay<double>(fs.Rv());
182 const auto& priVars = solution[index];
183 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
186 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
189 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
190 && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
191 denominator = Toolbox::template decay<LhsEval>(1 - rs * rv);
194 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
195 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
196 FluidSystem::solventComponentIndex(FluidSystem::oilPhaseIdx));
197 bweights[activeCompIdx] = Toolbox::template decay<LhsEval>(
198 (1 / fs.invB(FluidSystem::oilPhaseIdx) - rs / fs.invB(FluidSystem::gasPhaseIdx))
201 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
202 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
203 FluidSystem::solventComponentIndex(FluidSystem::gasPhaseIdx));
204 bweights[activeCompIdx] = Toolbox::template decay<LhsEval>(
205 (1 / fs.invB(FluidSystem::gasPhaseIdx) - rv / fs.invB(FluidSystem::oilPhaseIdx))
208 weights[index] = bweights;
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:192
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
void getTrueImpesWeightsAnalytic(int, Vector &weights, const GridView &gridView, ElementContext &elemCtx, const Model &model, std::size_t threadId)
Definition: getQuasiImpesWeights.hpp:137
void getTrueImpesWeights(int pressureVarIndex, Vector &weights, const GridView &gridView, ElementContext &elemCtx, const Model &model, std::size_t threadId)
Definition: getQuasiImpesWeights.hpp:91
void getQuasiImpesWeights(const Matrix &matrix, const int pressureVarIndex, const bool transpose, Vector &weights)
Definition: getQuasiImpesWeights.hpp:50
DenseMatrix transposeDenseMatrix(const DenseMatrix &M)
Definition: getQuasiImpesWeights.hpp:36
Definition: blackoilboundaryratevector.hh:37