From fe79a9fc07b0fc8324a8dc29b16f872099e9d320 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 15 Nov 2018 20:38:25 +0100 Subject: [PATCH 1/3] Refactor getConvergence() to use ConvergenceReport. Note: the communication and reduction for computing reservoir convergence is not done by gathering ConvergenceReports, but as before, using the convergenceReduction() method. --- opm/autodiff/BlackoilModelEbos.hpp | 110 +++++++++++++----------- opm/autodiff/BlackoilWellModel_impl.hpp | 8 -- 2 files changed, 61 insertions(+), 57 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 7c5e4b8ba..41c3095f9 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -780,64 +780,73 @@ 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& residual_norms) + ConvergenceReport getReservoirConvergence(const double dt, + const int iteration, + std::vector& B_avg, + std::vector& 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 CNV(numComp); + std::vector 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(); - - bool converged = converged_MB && converged_well; - - converged = converged && converged_CNV; + // Create convergence report. + ConvergenceReport report; + using CR = ConvergenceReport; + 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 int compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx); + 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}); + OpmLog::debug("NaN residual for " + compName + " equation."); + } else if (res[ii] > maxResidualAllowed()) { + report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); + OpmLog::debug("Too large residual for " + compName + " equation."); + } else if (res[ii] < 0.0) { + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + OpmLog::debug("Negative residual for " + compName + " 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 ); + std::vector< std::string > key( numEq ); for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { continue; @@ -883,29 +892,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 + bool getConvergence(const SimulatorTimerInterface& timer, + const int iteration, + std::vector& residual_norms) + { + // Get convergence reports for reservoir and wells. + std::vector 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.converged(); } diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 4269debbe..9c877c791 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -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; } From b42165b56080deead5f5c39388bb04a618222138 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 15 Nov 2018 21:21:25 +0100 Subject: [PATCH 2/3] Store convergence reports from all steps and iterations. --- opm/autodiff/BlackoilModelEbos.hpp | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 41c3095f9..3548bcad8 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -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()) { @@ -900,9 +907,9 @@ namespace Opm { /// \param[in] timer simulation timer /// \param[in] iteration current iteration number /// \param[out] residual_norms CNV residuals by phase - bool getConvergence(const SimulatorTimerInterface& timer, - const int iteration, - std::vector& residual_norms) + ConvergenceReport getConvergence(const SimulatorTimerInterface& timer, + const int iteration, + std::vector& residual_norms) { // Get convergence reports for reservoir and wells. std::vector B_avg(numEq, 0.0); @@ -917,7 +924,7 @@ namespace Opm { OPM_THROW(Opm::NumericalIssue, "Too large residual found!"); } - return report.converged(); + return report; } @@ -991,6 +998,8 @@ namespace Opm { std::unique_ptr matrix_for_preconditioner_; std::vector>> overlapRowAndColumns_; + + std::vector> convergence_reports_; public: /// return the StandardWells object BlackoilWellModel& From 83ce1eb9197e51cfbacdd94cf63b718b398bb110 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 19 Nov 2018 14:46:31 +0100 Subject: [PATCH 3/3] Ensure all residuals are checked for convergence. Looping over all components instead of phases, to handle polymer etc. correctly. Also slight refactoring of how component names for output are handled. --- opm/autodiff/BlackoilModelEbos.hpp | 75 ++++++++++++++++-------------- 1 file changed, 39 insertions(+), 36 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 3548bcad8..9d33891c2 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -816,16 +816,33 @@ namespace Opm { residual_norms.push_back(CNV[compIdx]); } + // Setup component names, only the first time the function is run. + static std::vector 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"; + } + } + // Create convergence report. ConvergenceReport report; using CR = ConvergenceReport; - 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 int compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx); + 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 }; @@ -833,13 +850,19 @@ namespace Opm { for (int ii : {0, 1}) { if (std::isnan(res[ii])) { report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx}); - OpmLog::debug("NaN residual for " + compName + " equation."); + if ( terminal_output_ ) { + OpmLog::debug("NaN residual for " + compNames[compIdx] + " equation."); + } } else if (res[ii] > maxResidualAllowed()) { report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); - OpmLog::debug("Too large residual for " + compName + " equation."); + 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}); - OpmLog::debug("Negative residual for " + compName + " equation."); + 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}); } @@ -852,35 +875,15 @@ namespace Opm { // Only rank 0 does print to std::cout if (iteration == 0) { std::string msg = "Iter"; - - std::vector< std::string > key( numEq ); - 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); }