fixed true impes for new code

This commit is contained in:
hnil 2022-06-22 18:27:17 +02:00 committed by Atgeirr Flø Rasmussen
parent 4f6755025c
commit f751279841
2 changed files with 18 additions and 8 deletions

View File

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

View File

@ -87,7 +87,7 @@ namespace Amg
return weights;
}
template<class Vector, class GridView, class ElementContext, class Model>
template<class Evaluation, 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;
@ -131,7 +131,7 @@ namespace Amg
OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
}
template<class Vector, class Model, class Evaluation>
template<class Evaluation, class Vector, class Model>
void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const Model& model)
{
using VectorBlockType = typename Vector::block_type;
@ -150,7 +150,7 @@ namespace Amg
const auto* intQuantsInP = model.cachedIntensiveQuantities(globI, /*timeIdx*/0);
assert(intQuantsInP);
const auto& intQuantsIn = *intQuantsInP;
//NB !!!!! LocalResidual::computeStorage(storage,intQuantsIn, 0);
Model::LocalResidual::computeStorage(storage,intQuantsIn, 0);
double scvVolume = model.dofTotalVolume(globI);
double dt = 3600*24;
auto storage_scale = scvVolume / dt;