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

Merge pull request #5103 from hnil/trueimpesanalytic

Browse Source 
add analytic true impes
This commit is contained in:
Atgeirr Flø Rasmussen 2024-01-18 22:51:57 +01:00 committed by GitHub
commit ae838d0b84
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
3 changed files with 138 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
opm/simulators/linalg/ISTLSolverEbos.cpp
View File

@ -81,49 +81,12 @@ void makeOverlapRowsInvalid(Matrix& matrix,
}
}
/// Return an appropriate weight function if a cpr preconditioner is asked for.
template<class Vector, class Matrix>
std::function<Vector()> getWeightsCalculator(const PropertyTree& prm,
const Matrix& matrix,
std::size_t pressureIndex,
std::function<Vector()> trueFunc)
{
std::function<Vector()> weightsCalculator;
using namespace std::string_literals;
auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
if (preconditionerType == "cpr" || preconditionerType == "cprt"
|| preconditionerType == "cprw" || preconditionerType == "cprwt") {
const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
if (weightsType == "quasiimpes") {
// weights will be created as default in the solver
// assignment p = pressureIndex prevent compiler warning about
// capturing variable with non-automatic storage duration
weightsCalculator = [matrix, transpose, pressureIndex]() {
return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
pressureIndex,
transpose);
};
} else if (weightsType == "trueimpes") {
weightsCalculator = trueFunc;
} else {
OPM_THROW(std::invalid_argument,
"Weights type " + weightsType +
"not implemented for cpr."
" Please use quasiimpes or trueimpes.");
}
}
return weightsCalculator;
}
template<class Matrix, class Vector, class Comm>
void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
bool parallel,
const PropertyTree& prm,
std::size_t pressureIndex,
std::function<Vector()> trueFunc,
std::function<Vector()> weightsCalculator,
const bool forceSerial,
[[maybe_unused]] Comm& comm)
@ -152,9 +115,6 @@ void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
}
}
std::function<Vector()> weightsCalculator =
getWeightsCalculator<Vector>(prm, matrix, pressureIndex, trueFunc);
if (parallel) {
#if HAVE_MPI
if (!wellOperator_) {

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

@ -102,7 +102,7 @@ struct FlexibleSolverInfo
bool parallel,
const PropertyTree& prm,
std::size_t pressureIndex,
std::function<Vector()> trueFunc,
std::function<Vector()> weightCalculator,
const bool forceSerial,
Comm& comm);
@ -460,22 +460,17 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
{
OPM_TIMEBLOCK(flexibleSolverPrepare);
if (shouldCreateSolver()) {
std::function<Vector()> trueFunc =
[this]
{
return this->getTrueImpesWeights(pressureIndex);
};
if (!useWellConn_) {
auto wellOp = std::make_unique<WellModelOperator>(simulator_.problem().wellModel());
flexibleSolver_[activeSolverNum_].wellOperator_ = std::move(wellOp);
}
std::function<Vector()> weightCalculator = this->getWeightsCalculator(prm_[activeSolverNum_], getMatrix(), pressureIndex);
OPM_TIMEBLOCK(flexibleSolverCreate);
flexibleSolver_[activeSolverNum_].create(getMatrix(),
isParallel(),
prm_[activeSolverNum_],
pressureIndex,
trueFunc,
weightCalculator,
forceSerial_,
*comm_);
}
@ -540,17 +535,62 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
// 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
std::function<Vector()> getWeightsCalculator(const PropertyTree& prm,
const Matrix& matrix,
std::size_t pressIndex) const
{
std::function<Vector()> weightsCalculator;
using namespace std::string_literals;
auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
if (preconditionerType == "cpr" || preconditionerType == "cprt"
|| preconditionerType == "cprw" || preconditionerType == "cprwt") {
const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
if (weightsType == "quasiimpes") {
// weights will be created as default in the solver
// assignment p = pressureIndex prevent compiler warning about
// capturing variable with non-automatic storage duration
weightsCalculator = [matrix, transpose, pressIndex]() {
return Amg::getQuasiImpesWeights<Matrix, Vector>(matrix,
pressIndex,
transpose);
};
} else if ( weightsType == "trueimpes" ) {
weightsCalculator =
[this, pressIndex]
{
OPM_TIMEBLOCK(getTrueImpesWeights);
Vector weights(rhs_->size());
ElementContext elemCtx(simulator_);
Amg::getTrueImpesWeights(pressureVarIndex, weights,
Amg::getTrueImpesWeights(pressIndex, weights,
simulator_.vanguard().gridView(),
elemCtx, simulator_.model(),
ThreadManager::threadId());
return weights;
};
} else if (weightsType == "trueimpesanalytic" ) {
weightsCalculator =
[this, pressIndex]
{
Vector weights(rhs_->size());
ElementContext elemCtx(simulator_);
Amg::getTrueImpesWeightsAnalytic(pressIndex, weights,
simulator_.vanguard().gridView(),
elemCtx, simulator_.model(),
ThreadManager::threadId());
return weights;
};
} else {
OPM_THROW(std::invalid_argument,
"Weights type " + weightsType +
"not implemented for cpr."
" Please use quasiimpes or trueimpes.");
}
}
return weightsCalculator;
}
Matrix& getMatrix()
{

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

@ -23,7 +23,7 @@
#include <dune/common/fvector.hh>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <algorithm>
#include <cmath>
@ -132,6 +132,83 @@ namespace Amg
}
OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
}
template <class Vector, class GridView, class ElementContext, class Model>
void getTrueImpesWeightsAnalytic(int /*pressureVarIndex*/,
Vector& weights,
const GridView& gridView,
ElementContext& elemCtx,
const Model& model,
std::size_t threadId)
{
// The sequential residual is a linear combination of the
// mass balance residuals, with coefficients equal to (for
// water, oil, gas):
// 1/bw,
// (1/bo - rs/bg)/(1-rs*rv)
// (1/bg - rv/bo)/(1-rs*rv)
// These coefficients must be applied for both the residual and
// Jacobian.
using FluidSystem = typename Model::FluidSystem;
using LhsEval = double;
using Indices = typename Model::Indices;
using PrimaryVariables = typename Model::PrimaryVariables;
using VectorBlockType = typename Vector::block_type;
using Evaluation =
typename std::decay_t<decltype(model.localLinearizer(threadId).localResidual().residual(0))>::block_type;
using Toolbox = MathToolbox<Evaluation>;
VectorBlockType rhs(0.0);
const auto& solution = model.solution(/*timeIdx*/ 0);
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (const auto& elem : elements(gridView)) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& index = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
VectorBlockType bweights;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
FluidSystem::solventComponentIndex(FluidSystem::waterPhaseIdx));
bweights[activeCompIdx]
= Toolbox::template decay<LhsEval>(1 / fs.invB(FluidSystem::waterPhaseIdx));
}
double denominator = 1.0;
double rs = Toolbox::template decay<double>(fs.Rs());
double rv = Toolbox::template decay<double>(fs.Rv());
const auto& priVars = solution[index];
if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
rs = 0.0;
}
if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
rv = 0.0;
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
&& FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
denominator = Toolbox::template decay<LhsEval>(1 - rs * rv);
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
FluidSystem::solventComponentIndex(FluidSystem::oilPhaseIdx));
bweights[activeCompIdx] = Toolbox::template decay<LhsEval>(
(1 / fs.invB(FluidSystem::oilPhaseIdx) - rs / fs.invB(FluidSystem::gasPhaseIdx))
/ denominator);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(
FluidSystem::solventComponentIndex(FluidSystem::gasPhaseIdx));
bweights[activeCompIdx] = Toolbox::template decay<LhsEval>(
(1 / fs.invB(FluidSystem::gasPhaseIdx) - rv / fs.invB(FluidSystem::oilPhaseIdx))
/ denominator);
}
weights[index] = bweights;
}
OPM_END_PARALLEL_TRY_CATCH("getTrueImpesAnalyticWeights() failed: ", elemCtx.simulator().vanguard().grid().comm());
}
} // namespace Amg
} // namespace Opm