diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index cac48e5cc..49e059aeb 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -495,8 +495,73 @@ namespace Opm const std::vector& B_avg, const ModelParameters& param) const { - // TODO: it will be very similar + // assert((int(B_avg.size()) == numComponents()) || has_polymer); + assert( (int(B_avg.size()) == numComponents()) ); + + // checking if any residual is NaN or too large. The two large one is only handled for the well flux + std::vector> residual(numberOfSegments(), std::vector(numWellEq, 0.0)); + for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]); + } + } + + std::vector maximum_residual(numComponents(), 0.0); + ConvergenceReport report; + // TODO: the following is a little complicated, maybe can be simplified in some way? + for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + if (eq_idx < numComponents()) { // phase or component mass equations + const double flux_residual = B_avg[eq_idx] * residual[seg][eq_idx]; + // TODO: the report can not handle the segment number yet. + if (std::isnan(flux_residual)) { + report.nan_residual_found = true; + const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); + const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; + report.nan_residual_wells.push_back(problem_well); + } else if (flux_residual > param.max_residual_allowed_) { + report.too_large_residual_found = true; + const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx)); + const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name}; + report.nan_residual_wells.push_back(problem_well); + } else { // it is a normal residual + if (flux_residual > maximum_residual[eq_idx]) { + maximum_residual[eq_idx] = flux_residual; + } + } + } else { // pressure equation + const double pressure_residal = residual[seg][eq_idx]; + const std::string eq_name("Pressure"); + if (std::isnan(pressure_residal)) { + report.nan_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name}; + report.nan_residual_wells.push_back(problem_well); + } else if (std::isinf(pressure_residal)) { + report.too_large_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), eq_name}; + report.nan_residual_wells.push_back(problem_well); + } else { // it is a normal residual + if (pressure_residal > maximum_residual[eq_idx]) { + maximum_residual[eq_idx] = pressure_residal; + } + } + } + } + } + + if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found + // check convergence + for ( int comp_idx = 0; comp_idx < numComponents(); ++comp_idx) + { + report.converged = report.converged && (maximum_residual[comp_idx] < param.tolerance_wells_); + } + + report.converged = report.converged && (maximum_residual[SPres] < param.tolerance_well_control_); + } else { // abnormal values found and no need to check the convergence + report.converged = false; + } + return report; } diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index a6e1a254c..8c2c06110 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -408,7 +408,6 @@ namespace Opm wellhelpers::WellSwitchingLogger& logger) const { const int np = number_of_phases_; - const int nw = well_state.bhp().size(); const int w = index_of_well_; const int old_control_index = well_state.currentControls()[w];