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

Refactor getConvergence() to use ConvergenceReport.

Browse Source 
Note: the communication and reduction for computing reservoir
convergence is not done by gathering ConvergenceReports, but
as before, using the convergenceReduction() method.
This commit is contained in:
Atgeirr Flø Rasmussen 2018-11-15 20:38:25 +01:00
parent c006ea23f2
commit fe79a9fc07
2 changed files with 61 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

110
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -780,64 +780,73 @@ namespace Opm {
return pvSumLocal; return pvSumLocal;
} }
/// Compute convergence based on total mass balance (tol_mb) and maximum ConvergenceReport getReservoirConvergence(const double dt,
/// residual mass balance (tol_cnv). const int iteration,
/// \param[in] timer simulation timer std::vector<Scalar>& B_avg,
/// \param[in] dt timestep length std::vector<Scalar>& residual_norms)
/// \param[in] iteration current iteration number
bool getConvergence(const SimulatorTimerInterface& timer, const int iteration, std::vector<double>& residual_norms)
{ {
typedef std::vector< Scalar > Vector; typedef std::vector< Scalar > Vector;
const double dt = timer.currentStepLength(); const double tol_mb = param_.tolerance_mb_;
const double tol_mb = param_.tolerance_mb_; const double tol_cnv = (iteration < param_.max_strict_iter_) ? param_.tolerance_cnv_ : param_.tolerance_cnv_relaxed_;
const double tol_cnv = param_.tolerance_cnv_;
const double tol_cnv_relaxed = param_.tolerance_cnv_relaxed_;
const int numComp = numEq; const int numComp = numEq;
Vector R_sum(numComp, 0.0 ); Vector R_sum(numComp, 0.0 );
Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() ); Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
Vector B_avg(numComp, 0.0 );
const double pvSumLocal = localConvergenceData(R_sum, maxCoeff, B_avg); const double pvSumLocal = localConvergenceData(R_sum, maxCoeff, B_avg);
// compute global sum and max of quantities // compute global sum and max of quantities
const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal, const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
R_sum, maxCoeff, B_avg); R_sum, maxCoeff, B_avg);
Vector CNV(numComp);
Vector mass_balance_residual(numComp);
bool converged_MB = true;
bool converged_CNV = true;
// Finish computation // Finish computation
std::vector<Scalar> CNV(numComp);
std::vector<Scalar> mass_balance_residual(numComp);
for ( int compIdx = 0; compIdx < numComp; ++compIdx ) for ( int compIdx = 0; compIdx < numComp; ++compIdx )
{ {
CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx]; CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx];
mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum; 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]); residual_norms.push_back(CNV[compIdx]);
} }
const auto report_well = wellModel().getWellConvergence(B_avg); // Create convergence report.
const bool converged_well = report_well.converged(); ConvergenceReport report;
using CR = ConvergenceReport;
bool converged = converged_MB && converged_well; for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
converged = converged && converged_CNV; 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_ ) if ( terminal_output_ )
{ {
// Only rank 0 does print to std::cout // Only rank 0 does print to std::cout
if (iteration == 0) { if (iteration == 0) {
std::string msg = "Iter"; std::string msg = "Iter";
std::vector< std::string > key( numComp ); std::vector< std::string > key( numEq );
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) { if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue; continue;
@ -883,29 +892,32 @@ namespace Opm {
OpmLog::debug(ss.str()); OpmLog::debug(ss.str());
} }
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { return report;
if (!FluidSystem::phaseIsActive(phaseIdx)) }
continue;
const unsigned canonicalCompIdx = FluidSystem::solventComponentIndex(phaseIdx); /// Compute convergence based on total mass balance (tol_mb) and maximum
const std::string& compName = FluidSystem::componentName(canonicalCompIdx); /// residual mass balance (tol_cnv).
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(canonicalCompIdx); /// \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<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]) // Throw if any NaN or too large residual found.
|| std::isnan(CNV[compIdx])) { ConvergenceReport::Severity severity = report.severityOfWorstFailure();
OPM_THROW(Opm::NumericalIssue, "NaN residual for " << compName << " equation"); if (severity == ConvergenceReport::Severity::NotANumber) {
} OPM_THROW(Opm::NumericalIssue, "NaN residual found!");
if (mass_balance_residual[compIdx] > maxResidualAllowed() } else if (severity == ConvergenceReport::Severity::TooLarge) {
|| CNV[compIdx] > maxResidualAllowed()) { OPM_THROW(Opm::NumericalIssue, "Too large residual found!");
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");
}
} }
return converged; return report.converged();
} }

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; return report;
} }