OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-11 08:51:55 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1639 from atgeirr/store-convrep

Browse Source 
Store ConvergenceReport objects
This commit is contained in:
Atgeirr Flø Rasmussen 2018-11-19 15:13:14 +01:00 committed by GitHub
commit 4337483ac6
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 98 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

172
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -173,6 +173,7 @@ namespace Opm {
{
OPM_THROW(std::logic_error,"solver down cast to ISTLSolver failed");
}
convergence_reports_.reserve(300); // Often insufficient, but avoids frequent moves.
}
bool isParallel() const
@ -242,6 +243,8 @@ namespace Opm {
residual_norms_history_.clear();
current_relaxation_ = 1.0;
dx_old_ = 0.0;
convergence_reports_.push_back({});
convergence_reports_.back().reserve(11);
}
report.total_linearizations = 1;
@ -261,7 +264,11 @@ namespace Opm {
perfTimer.reset();
perfTimer.start();
// the step is not considered converged until at least minIter iterations is done
report.converged = getConvergence(timer, iteration,residual_norms) && iteration > nonlinear_solver.minIter();
{
auto convrep = getConvergence(timer, iteration,residual_norms);
report.converged = convrep.converged() && iteration > nonlinear_solver.minIter();;
convergence_reports_.back().push_back(std::move(convrep));
}
// checking whether the group targets are converged
if (wellModel().wellCollection().groupControlActive()) {
@ -780,91 +787,103 @@ namespace Opm {
return pvSumLocal;
}
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
/// \param[in] timer simulation timer
/// \param[in] dt timestep length
/// \param[in] iteration current iteration number
bool getConvergence(const SimulatorTimerInterface& timer, const int iteration, std::vector<double>& residual_norms)
ConvergenceReport getReservoirConvergence(const double dt,
const int iteration,
std::vector<Scalar>& B_avg,
std::vector<Scalar>& residual_norms)
{
typedef std::vector< Scalar > Vector;
const double dt = timer.currentStepLength();
const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = param_.tolerance_cnv_;
const double tol_cnv_relaxed = param_.tolerance_cnv_relaxed_;
const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = (iteration < param_.max_strict_iter_) ? param_.tolerance_cnv_ : param_.tolerance_cnv_relaxed_;
const int numComp = numEq;
Vector R_sum(numComp, 0.0 );
Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
Vector B_avg(numComp, 0.0 );
const double pvSumLocal = localConvergenceData(R_sum, maxCoeff, B_avg);
// compute global sum and max of quantities
const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
R_sum, maxCoeff, B_avg);
Vector CNV(numComp);
Vector mass_balance_residual(numComp);
bool converged_MB = true;
bool converged_CNV = true;
// Finish computation
std::vector<Scalar> CNV(numComp);
std::vector<Scalar> mass_balance_residual(numComp);
for ( int compIdx = 0; compIdx < numComp; ++compIdx )
{
CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[compIdx] < tol_mb);
if (iteration < param_.max_strict_iter_)
converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv);
else
converged_CNV = converged_CNV && (CNV[compIdx] < tol_cnv_relaxed);
residual_norms.push_back(CNV[compIdx]);
}
const auto report_well = wellModel().getWellConvergence(B_avg);
const bool converged_well = report_well.converged();
// Setup component names, only the first time the function is run.
static std::vector<std::string> compNames;
if (compNames.empty()) {
compNames.resize(numComp);
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
compNames[compIdx] = FluidSystem::componentName(canonicalCompIdx);
}
if (has_solvent_) {
compNames[solventSaturationIdx] = "Solvent";
}
if (has_polymer_) {
compNames[polymerConcentrationIdx] = "Polymer";
}
if (has_energy_) {
compNames[temperatureIdx] = "Energy";
}
}
bool converged = converged_MB && converged_well;
converged = converged && converged_CNV;
// Create convergence report.
ConvergenceReport report;
using CR = ConvergenceReport;
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance,
CR::ReservoirFailure::Type::Cnv };
double tol[2] = { tol_mb, tol_cnv };
for (int ii : {0, 1}) {
if (std::isnan(res[ii])) {
report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx});
if ( terminal_output_ ) {
OpmLog::debug("NaN residual for " + compNames[compIdx] + " equation.");
}
} else if (res[ii] > maxResidualAllowed()) {
report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx});
if ( terminal_output_ ) {
OpmLog::debug("Too large residual for " + compNames[compIdx] + " equation.");
}
} else if (res[ii] < 0.0) {
report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
if ( terminal_output_ ) {
OpmLog::debug("Negative residual for " + compNames[compIdx] + " equation.");
}
} else if (res[ii] > tol[ii]) {
report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
}
}
}
// Output of residuals.
if ( terminal_output_ )
{
// Only rank 0 does print to std::cout
if (iteration == 0) {
std::string msg = "Iter";
std::vector< std::string > key( numComp );
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
const std::string& compName = FluidSystem::componentName(canonicalCompIdx);
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
key[ compIdx ] = std::toupper( compName.front() );
}
if (has_solvent_) {
key[ solventSaturationIdx ] = "S";
}
if (has_polymer_) {
key[ polymerConcentrationIdx ] = "P";
}
if (has_energy_) {
key[ temperatureIdx ] = "E";
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " MB(" + key[ compIdx ] + ") ";
msg += " MB(";
msg += compNames[compIdx][0];
msg += ") ";
}
for (int compIdx = 0; compIdx < numComp; ++compIdx) {
msg += " CNV(" + key[ compIdx ] + ") ";
msg += " CNV(";
msg += compNames[compIdx][0];
msg += ") ";
}
OpmLog::debug(msg);
}
@ -883,29 +902,32 @@ namespace Opm {
OpmLog::debug(ss.str());
}
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
return report;
}
const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
const std::string& compName = FluidSystem::componentName(canonicalCompIdx);
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx);
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
/// \param[in] timer simulation timer
/// \param[in] iteration current iteration number
/// \param[out] residual_norms CNV residuals by phase
ConvergenceReport getConvergence(const SimulatorTimerInterface& timer,
const int iteration,
std::vector<double>& residual_norms)
{
// Get convergence reports for reservoir and wells.
std::vector<Scalar> B_avg(numEq, 0.0);
auto report = getReservoirConvergence(timer.currentStepLength(), iteration, B_avg, residual_norms);
report += wellModel().getWellConvergence(B_avg);
if (std::isnan(mass_balance_residual[compIdx])
|| std::isnan(CNV[compIdx])) {
OPM_THROW(Opm::NumericalIssue, "NaN residual for " << compName << " equation");
}
if (mass_balance_residual[compIdx] > maxResidualAllowed()
|| CNV[compIdx] > maxResidualAllowed()) {
OPM_THROW(Opm::NumericalIssue, "Too large residual for " << compName << " equation");
}
if (mass_balance_residual[compIdx] < 0
|| CNV[compIdx] < 0) {
OPM_THROW(Opm::NumericalIssue, "Negative residual for " << compName << " equation");
}
// Throw if any NaN or too large residual found.
ConvergenceReport::Severity severity = report.severityOfWorstFailure();
if (severity == ConvergenceReport::Severity::NotANumber) {
OPM_THROW(Opm::NumericalIssue, "NaN residual found!");
} else if (severity == ConvergenceReport::Severity::TooLarge) {
OPM_THROW(Opm::NumericalIssue, "Too large residual found!");
}
return converged;
return report;
}
@ -979,6 +1001,8 @@ namespace Opm {
std::unique_ptr<Mat> matrix_for_preconditioner_;
std::vector<std::pair<int,std::vector<int>>> overlapRowAndColumns_;
std::vector<std::vector<ConvergenceReport>> convergence_reports_;
public:
/// return the StandardWells object
BlackoilWellModel<TypeTag>&

8
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@ -816,14 +816,6 @@ namespace Opm {
}
}
// Throw if any NaN or too large residual found.
ConvergenceReport::Severity severity = report.severityOfWorstFailure();
if (severity == ConvergenceReport::Severity::NotANumber) {
OPM_THROW(Opm::NumericalIssue, "NaN residual found!");
} else if (severity == ConvergenceReport::Severity::TooLarge) {
OPM_THROW(Opm::NumericalIssue, "Too large residual found!");
}
return report;
}