Merge pull request #4208 from akva2/well_convergence

Add class for checking well convergence
This commit is contained in:
Bård Skaflestad 2022-10-28 21:26:17 +02:00 committed by GitHub
commit d1d5ff2831
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
9 changed files with 253 additions and 242 deletions

View File

@ -102,6 +102,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/VFPHelpers.cpp
opm/simulators/wells/VFPProdProperties.cpp
opm/simulators/wells/VFPInjProperties.cpp
opm/simulators/wells/WellConvergence.cpp
opm/simulators/wells/WellGroupHelpers.cpp
opm/simulators/wells/WellHelpers.cpp
opm/simulators/wells/WellInterfaceEval.cpp
@ -380,6 +381,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/VFPProdProperties.hpp
opm/simulators/wells/VFPProperties.hpp
opm/simulators/wells/WellConnectionAuxiliaryModule.hpp
opm/simulators/wells/WellConvergence.hpp
opm/simulators/wells/WellGroupHelpers.hpp
opm/simulators/wells/WellHelpers.hpp
opm/simulators/wells/WellInterface.hpp

View File

@ -33,6 +33,7 @@
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
@ -142,89 +143,6 @@ initPrimaryVariablesEvaluation() const
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
const double tolerance_pressure_ms_wells,
const double tolerance_wells,
const double max_residual_allowed,
DeferredLogger& deferred_logger) const
{
double control_tolerance = 0.;
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
const int well_index = baseif_.indexOfWell();
const auto& ws = well_state.well(well_index);
if (baseif_.isInjector() )
{
auto current = ws.injection_cmode;
switch(current) {
case Well::InjectorCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = tolerance_pressure_ms_wells;
break;
case Well::InjectorCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = tolerance_pressure_ms_wells;
break;
case Well::InjectorCMode::RATE:
case Well::InjectorCMode::RESV:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerance_wells;
break;
case Well::InjectorCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerance_wells;
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger);
}
}
if (baseif_.isProducer() )
{
auto current = ws.production_cmode;
switch(current) {
case Well::ProducerCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = tolerance_pressure_ms_wells;
break;
case Well::ProducerCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = tolerance_pressure_ms_wells;
break;
case Well::ProducerCMode::ORAT:
case Well::ProducerCMode::WRAT:
case Well::ProducerCMode::GRAT:
case Well::ProducerCMode::LRAT:
case Well::ProducerCMode::RESV:
case Well::ProducerCMode::CRAT:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerance_wells;
break;
case Well::ProducerCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerance_wells;
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger);
}
}
const double well_control_residual = std::abs(resWell_[0][SPres]);
const int dummy_component = -1;
if (std::isnan(well_control_residual)) {
report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()});
} else if (well_control_residual > max_residual_allowed * 10.) {
report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, baseif_.name()});
} else if ( well_control_residual > control_tolerance) {
report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, baseif_.name()});
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
ConvergenceReport
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
@ -302,12 +220,16 @@ getWellConvergence(const WellState& well_state,
}
}
checkConvergenceControlEq(well_state,
report,
tolerance_pressure_ms_wells,
tolerance_wells,
max_residual_allowed,
deferred_logger);
WellConvergence(baseif_).
checkConvergenceControlEq(well_state,
{tolerance_pressure_ms_wells,
tolerance_pressure_ms_wells,
tolerance_wells,
tolerance_wells,
max_residual_allowed},
std::abs(resWell_[0][SPres]),
report,
deferred_logger);
return report;
}

View File

@ -144,13 +144,6 @@ protected:
WellState& well_state,
DeferredLogger& deferred_logger) const;
void checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
const double tolerance_pressure_ms_wells,
const double tolerance_wells,
const double max_residual_allowed,
DeferredLogger& deferred_logger) const;
/// check whether the well equations get converged for this well
ConvergenceReport getWellConvergence(const WellState& well_state,
const std::vector<double>& B_avg,

View File

@ -32,6 +32,7 @@
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
@ -47,7 +48,7 @@ namespace Opm
template<class FluidSystem, class Indices, class Scalar>
StandardWellEval<FluidSystem,Indices,Scalar>::
StandardWellEval(const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
: StandardWellGeneric<Scalar>(Bhp, baseif)
: StandardWellGeneric<Scalar>(baseif)
, baseif_(baseif)
, F0_(numWellConservationEq)
{
@ -831,7 +832,12 @@ getWellConvergence(const WellState& well_state,
}
}
this->checkConvergenceControlEq(well_state, report, deferred_logger, maxResidualAllowed);
WellConvergence(baseif_).
checkConvergenceControlEq(well_state,
{1.e3, 1.e4, 1.e-4, 1.e-6, maxResidualAllowed},
std::abs(this->resWell_[0][Bhp]),
report,
deferred_logger);
return report;
}

View File

@ -47,13 +47,11 @@ namespace Opm
template<class Scalar>
StandardWellGeneric<Scalar>::
StandardWellGeneric(int Bhp,
const WellInterfaceGeneric& baseif)
StandardWellGeneric(const WellInterfaceGeneric& baseif)
: baseif_(baseif)
, perf_densities_(baseif_.numPerfs())
, perf_pressure_diffs_(baseif_.numPerfs())
, parallelB_(duneB_, baseif_.parallelWellInfo())
, Bhp_(Bhp)
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise);
@ -393,131 +391,6 @@ computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>&
}
}
template<class Scalar>
void
StandardWellGeneric<Scalar>::
checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
DeferredLogger& deferred_logger,
const double max_residual_allowed) const
{
double control_tolerance = 0.;
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
const int well_index = baseif_.indexOfWell();
const auto& ws = well_state.well(well_index);
if (baseif_.wellIsStopped()) {
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-6; // use smaller tolerance for zero control?
}
else if (baseif_.isInjector() )
{
auto current = ws.injection_cmode;
switch(current) {
case Well::InjectorCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = 1.e4; // 0.1 bar
break;
case Well::InjectorCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = 1.e3; // 0.01 bar
break;
case Well::InjectorCMode::RATE:
case Well::InjectorCMode::RESV:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-4; //
break;
case Well::InjectorCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-6; //
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger);
}
}
else if (baseif_.isProducer() )
{
auto current = ws.production_cmode;
switch(current) {
case Well::ProducerCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = 1.e4; // 0.1 bar
break;
case Well::ProducerCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = 1.e3; // 0.01 bar
break;
case Well::ProducerCMode::ORAT:
case Well::ProducerCMode::WRAT:
case Well::ProducerCMode::GRAT:
case Well::ProducerCMode::LRAT:
case Well::ProducerCMode::RESV:
case Well::ProducerCMode::CRAT:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-4; // smaller tolerance for rate control
break;
case Well::ProducerCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-6; // smaller tolerance for rate control
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger);
}
}
const double well_control_residual = std::abs(this->resWell_[0][Bhp_]);
const int dummy_component = -1;
if (std::isnan(well_control_residual)) {
report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()});
} else if (well_control_residual > max_residual_allowed * 10.) {
report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, baseif_.name()});
} else if ( well_control_residual > control_tolerance) {
report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, baseif_.name()});
}
}
template<class Scalar>
void
StandardWellGeneric<Scalar>::
checkConvergencePolyMW(const std::vector<double>& res,
ConvergenceReport& report,
const double maxResidualAllowed) const
{
if (baseif_.isInjector()) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const int dummy_component = -1;
using CR = ConvergenceReport;
const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
const double wat_vel_residual = res[Bhp_ + 1 + perf];
if (std::isnan(wat_vel_residual)) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, baseif_.name()});
} else if (wat_vel_residual > maxResidualAllowed * 10.) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::TooLarge, dummy_component, baseif_.name()});
} else if (wat_vel_residual > wat_vel_tol) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::Normal, dummy_component, baseif_.name()});
}
}
// checking the convergence of the skin pressure
const double pskin_tol = 1000.; // 1000 pascal
const auto pskin_failure_type = CR::WellFailure::Type::Pressure;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
const double pskin_residual = res[Bhp_ + 1 + perf + baseif_.numPerfs()];
if (std::isnan(pskin_residual)) {
report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, baseif_.name()});
} else if (pskin_residual > maxResidualAllowed * 10.) {
report.setWellFailed({pskin_failure_type, CR::Severity::TooLarge, dummy_component, baseif_.name()});
} else if (pskin_residual > pskin_tol) {
report.setWellFailed({pskin_failure_type, CR::Severity::Normal, dummy_component, baseif_.name()});
}
}
}
}
template<class Scalar>
void
StandardWellGeneric<Scalar>::

View File

@ -69,8 +69,7 @@ public:
void getNumBlocks(unsigned int& _nnzs) const;
protected:
StandardWellGeneric(int Bhp,
const WellInterfaceGeneric& baseif);
StandardWellGeneric(const WellInterfaceGeneric& baseif);
// calculate a relaxation factor to avoid overshoot of total rates
static double relaxationFactorRate(const std::vector<double>& primary_variables,
@ -85,16 +84,6 @@ protected:
const double bhp,
DeferredLogger& deferred_logger) const;
// checking the convergence of the well control equations
void checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
DeferredLogger& deferred_logger,
const double max_residual_allowed) const;
void checkConvergencePolyMW(const std::vector<double>& res,
ConvergenceReport& report,
const double maxResidualAllowed) const;
void computeConnectionPressureDelta();
std::optional<double> computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
@ -132,9 +121,6 @@ protected:
mutable BVectorWell invDrw_;
double getRho() const { return perf_densities_[0]; }
private:
int Bhp_; // index of Bhp
};
}

View File

@ -23,6 +23,7 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <algorithm>
#include <functional>
@ -2501,7 +2502,8 @@ namespace Opm
// checking the convergence of the extra equations related to polymer injectivity
if constexpr (Base::has_polymermw) {
this->checkConvergencePolyMW(res, report, this->param_.max_residual_allowed_);
WellConvergence(*this).
checkConvergencePolyMW(res, Bhp, this->param_.max_residual_allowed_, report);
}
}

View File

@ -0,0 +1,158 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 2017 IRIS AS.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/wells/WellConvergence.hpp>
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <cmath>
#include <stdexcept>
namespace Opm
{
void WellConvergence::
checkConvergenceControlEq(const WellState& well_state,
const Tolerances& tolerances,
const double well_control_residual,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const
{
double control_tolerance = 0.;
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
const int well_index = well_.indexOfWell();
const auto& ws = well_state.well(well_index);
if (well_.wellIsStopped()) {
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerances.rates; // use smaller tolerance for zero control?
}
else if (well_.isInjector() )
{
auto current = ws.injection_cmode;
switch(current) {
case Well::InjectorCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = tolerances.thp;
break;
case Well::InjectorCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = tolerances.bhp;
break;
case Well::InjectorCMode::RATE:
case Well::InjectorCMode::RESV:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerances.rates;
break;
case Well::InjectorCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerances.grup;
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << well_.name(), deferred_logger);
}
}
else if (well_.isProducer() )
{
auto current = ws.production_cmode;
switch(current) {
case Well::ProducerCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = tolerances.thp;
break;
case Well::ProducerCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = tolerances.bhp;
break;
case Well::ProducerCMode::ORAT:
case Well::ProducerCMode::WRAT:
case Well::ProducerCMode::GRAT:
case Well::ProducerCMode::LRAT:
case Well::ProducerCMode::RESV:
case Well::ProducerCMode::CRAT:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerances.rates;
break;
case Well::ProducerCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = tolerances.grup;
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << well_.name(), deferred_logger);
}
}
const int dummy_component = -1;
if (std::isnan(well_control_residual)) {
report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, well_.name()});
} else if (well_control_residual > tolerances.max_residual_allowed * 10.) {
report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, well_.name()});
} else if (well_control_residual > control_tolerance) {
report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, well_.name()});
}
}
void
WellConvergence::
checkConvergencePolyMW(const std::vector<double>& res,
const int Bhp,
const double maxResidualAllowed,
ConvergenceReport& report) const
{
if (well_.isInjector()) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const int dummy_component = -1;
using CR = ConvergenceReport;
const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
for (int perf = 0; perf < well_.numPerfs(); ++perf) {
const double wat_vel_residual = res[Bhp + 1 + perf];
if (std::isnan(wat_vel_residual)) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, well_.name()});
} else if (wat_vel_residual > maxResidualAllowed * 10.) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::TooLarge, dummy_component, well_.name()});
} else if (wat_vel_residual > wat_vel_tol) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::Normal, dummy_component, well_.name()});
}
}
// checking the convergence of the skin pressure
const double pskin_tol = 1000.; // 1000 pascal
const auto pskin_failure_type = CR::WellFailure::Type::Pressure;
for (int perf = 0; perf < well_.numPerfs(); ++perf) {
const double pskin_residual = res[Bhp + 1 + perf + well_.numPerfs()];
if (std::isnan(pskin_residual)) {
report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, well_.name()});
} else if (pskin_residual > maxResidualAllowed * 10.) {
report.setWellFailed({pskin_failure_type, CR::Severity::TooLarge, dummy_component, well_.name()});
} else if (pskin_residual > pskin_tol) {
report.setWellFailed({pskin_failure_type, CR::Severity::Normal, dummy_component, well_.name()});
}
}
}
}
}

View File

@ -0,0 +1,69 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 2017 IRIS AS.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELL_CONVERGENCE_HEADER_INCLUDED
#define OPM_WELL_CONVERGENCE_HEADER_INCLUDED
#include <vector>
namespace Opm
{
class ConvergenceReport;
class DeferredLogger;
class WellInterfaceGeneric;
class WellState;
class WellConvergence
{
public:
WellConvergence(const WellInterfaceGeneric& well)
: well_(well)
{}
struct Tolerances {
double bhp; //!< Tolerance for bhp controlled well
double thp; //!< Tolerance for thp controlled well
double rates; //!< Tolerance for a rate controlled well
double grup; //!< Tolerance for a grup controlled well
double max_residual_allowed; //!< Max residual allowd
};
// checking the convergence of the well control equations
void checkConvergenceControlEq(const WellState& well_state,
const Tolerances& tolerances,
const double well_control_residual,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const;
void checkConvergencePolyMW(const std::vector<double>& res,
const int Bhp,
const double maxResidualAllowed,
ConvergenceReport& report) const;
private:
const WellInterfaceGeneric& well_;
};
}
#endif // OPM_WELL_CONVERGENCE_HEADER_INCLUDED