OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Do not change true impes weights.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2022-08-09 09:56:58 +02:00
parent 727bf8d01f
commit 90e8a9af69
2 changed files with 11 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -92,7 +92,6 @@ namespace Opm
using AbstractPreconditionerType = Dune::PreconditionerWithUpdate<Vector, Vector>;
using WellModelOperator = WellModelAsLinearOperator<WellModel, Vector, Vector>;
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
constexpr static std::size_t pressureIndex = GetPropType<TypeTag, Properties::Indices>::pressureSwitchIdx;
#if HAVE_CUDA || HAVE_OPENCL || HAVE_FPGA || HAVE_AMGCL
@ -506,8 +505,7 @@ namespace Opm
// assignment p = pressureIndex prevent compiler warning about
// capturing variable with non-automatic storage duration
weightsCalculator = [this, p = pressureIndex]() {
ElementContext elemCtx(this->simulator_);
return this->getTrueImpesWeights(p, elemCtx);
return this->getTrueImpesWeights(p);
};
} else {
OPM_THROW(std::invalid_argument,
@ -522,24 +520,16 @@ namespace Opm
// Weights to make approximate pressure equations.
// Calculated from the storage terms (only) of the
// conservation equations, ignoring all other terms.
template<class ElemCtx>
Vector getTrueImpesWeights(int pressureVarIndex,ElemCtx& elemCtx) const
Vector getTrueImpesWeights(int pressureVarIndex) const
{
Vector weights(rhs_->size());
Amg::getTrueImpesWeights<Evaluation>(pressureVarIndex, weights, simulator_.vanguard().gridView(),
elemCtx, simulator_.model(),
ThreadManager::threadId());
ElementContext elemCtx(simulator_);
Amg::getTrueImpesWeights(pressureVarIndex, weights, simulator_.vanguard().gridView(),
elemCtx, simulator_.model(),
ThreadManager::threadId());
return weights;
}
#if 0
Vector getTrueImpesWeights(int pressureVarIndex,SmallElementContext<TypeTag>& /*elemCtx*/) const
{
Vector weights(rhs_->size());
Amg::getTrueImpesWeights<Evaluation>(pressureVarIndex, weights, simulator_.model());
return weights;
}
#endif
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1.0).

53
opm/simulators/linalg/getQuasiImpesWeights.hpp
View File

@ -87,7 +87,7 @@ namespace Amg
return weights;
}
template<class Evaluation, class Vector, class GridView, class ElementContext, class Model>
template<class Vector, class GridView, class ElementContext, class Model>
void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const GridView& gridView,
ElementContext& elemCtx, const Model& model, std::size_t threadId)
{
@ -95,8 +95,8 @@ namespace Amg
using Matrix = typename std::decay_t<decltype(model.linearizer().jacobian())>;
using MatrixBlockType = typename Matrix::MatrixBlock;
constexpr int numEq = VectorBlockType::size();
// using Evaluation = typename std::decay_t<decltype(model.localLinearizer(threadId).localResidual().residual(0))>
// ::block_type;
using Evaluation = typename std::decay_t<decltype(model.localLinearizer(threadId).localResidual().residual(0))>
::block_type;
VectorBlockType rhs(0.0);
rhs[pressureVarIndex] = 1.0;
int index = 0;
@ -108,8 +108,8 @@ namespace Amg
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
Dune::FieldVector<Evaluation, numEq> storage;
model.localLinearizer(threadId).localResidual().computeStorage(storage,elemCtx,/*spaceIdx=*/0, /*timeIdx=*/0);
//auto extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor();
auto scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume();// * extrusionFactor;
auto extrusionFactor = elemCtx.intensiveQuantities(0, /*timeIdx=*/0).extrusionFactor();
auto scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(0).volume() * extrusionFactor;
auto storage_scale = scvVolume / elemCtx.simulator().timeStepSize();
MatrixBlockType block;
double pressure_scale = 50e5;
@ -130,49 +130,6 @@ namespace Amg
}
OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
}
template<class Evaluation, class Vector, class Model>
void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const Model& model)
{
using VectorBlockType = typename Vector::block_type;
using Matrix = typename std::decay_t<decltype(model.linearizer().jacobian())>;
using MatrixBlockType = typename Matrix::MatrixBlock;
constexpr int numEq = VectorBlockType::size();
unsigned numCells = model.numTotalDof();
VectorBlockType rhs(0.0);
rhs[pressureVarIndex] = 1.0;
//NB !!OPM_BEGIN_PARALLEL_TRY_CATCH();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for(unsigned globI = 0; globI < numCells; globI++){
Dune::FieldVector<Evaluation, numEq> storage;
const auto* intQuantsInP = model.cachedIntensiveQuantities(globI, /*timeIdx*/0);
assert(intQuantsInP);
const auto& intQuantsIn = *intQuantsInP;
Model::LocalResidual::computeStorage(storage,intQuantsIn, 0);
double scvVolume = model.dofTotalVolume(globI);
double dt = 3600*24;
auto storage_scale = scvVolume / dt;
MatrixBlockType block;
double pressure_scale = 50e5;
for (int ii = 0; ii < numEq; ++ii) {
for (int jj = 0; jj < numEq; ++jj) {
block[ii][jj] = storage[ii].derivative(jj)/storage_scale;
if (jj == pressureVarIndex) {
block[ii][jj] *= pressure_scale;
}
}
}
VectorBlockType bweights;
MatrixBlockType block_transpose = Details::transposeDenseMatrix(block);
block_transpose.solve(bweights, rhs);
bweights /= 1000.0; // given normal densities this scales weights to about 1.
weights[globI] = bweights;
}
//NB!! OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
}
} // namespace Amg
} // namespace Opm