From 355be6c26c56289a220def2708c5ae180a884979 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 22 Aug 2017 14:49:30 +0200 Subject: [PATCH] collecting the NaN and too large well residuals make sure all the processes will throw if there is any of the processes found abnormal residual values. --- opm/autodiff/StandardWell.hpp | 8 +++-- opm/autodiff/StandardWell_impl.hpp | 33 +++++++++++------ opm/autodiff/StandardWellsDense.hpp | 2 ++ opm/autodiff/StandardWellsDense_impl.hpp | 39 +++++++++++++++++--- opm/autodiff/WellInterface.hpp | 45 ++++++++++++++++++++---- 5 files changed, 102 insertions(+), 25 deletions(-) diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 42ba78262..d1c4c0a67 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -58,6 +58,8 @@ namespace Opm }; using typename Base::Scalar; + using typename Base::ConvergenceReport; + using Base::numEq; using Base::has_solvent; @@ -123,9 +125,9 @@ namespace Opm wellhelpers::WellSwitchingLogger& logger) const; /// check whether the well equations get converged for this well - virtual bool getWellConvergence(Simulator& ebosSimulator, - const std::vector& B_avg, - const ModelParameters& param) const; + virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, + const std::vector& B_avg, + const ModelParameters& param) const; /// computing the accumulation term for later use in well mass equations virtual void computeAccumWell(); diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index c86e60b94..70f8782fa 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -1457,7 +1457,7 @@ namespace Opm template - bool + typename StandardWell::ConvergenceReport StandardWell:: getWellConvergence(Simulator& ebosSimulator, const std::vector& B_avg, @@ -1485,30 +1485,43 @@ namespace Opm } Vector well_flux_residual(numComp); - bool converged_Well = true; // Finish computation for ( int compIdx = 0; compIdx < numComp; ++compIdx ) { well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx]; - converged_Well = converged_Well && (well_flux_residual[compIdx] < tol_wells); } - // if one of the residuals is NaN, throw exception, so that the solver can be restarted + ConvergenceReport report; + // checking if any NaN or too large residuals found // TODO: not understand why phase here while component in other places. for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); if (std::isnan(well_flux_residual[phaseIdx])) { - OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << phaseName << " for well " << name()); - } - - if (well_flux_residual[phaseIdx] > maxResidualAllowed) { - OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << phaseName << " for well " << name()); + report.nan_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName}; + report.nan_residual_wells.push_back(problem_well); + } else { + if (well_flux_residual[phaseIdx] > maxResidualAllowed) { + report.too_large_residual_found = true; + const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName}; + report.too_large_residual_wells.push_back(problem_well); + } } } - return converged_Well; + if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found + // check convergence + for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + { + report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells); + } + } else { // abnormal values found and no need to check the convergence + report.converged = false; + } + + return report; } diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 2441d1ea2..48114afa5 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -168,6 +168,8 @@ namespace Opm { // later, might make share_ptr const later. std::vector well_container_; + using ConvergenceReport = typename WellInterface::ConvergenceReport; + // create the well container static std::vector createWellContainer(const Wells* wells, const std::vector& wells_ecl, diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index f343693db..ba841d2a5 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -479,22 +479,51 @@ namespace Opm { getWellConvergence(Simulator& ebosSimulator, const std::vector& B_avg) const { - bool converged_well = true; + ConvergenceReport report; - // currently, if there is any well not converged, we consider the well eqautions do not get converged for (const auto& well : well_container_) { - if ( !well->getWellConvergence(ebosSimulator, B_avg, param_) ) { - converged_well = false; + report += well->getWellConvergence(ebosSimulator, B_avg, param_); + } + + // checking NaN residuals + { + bool nan_residual_found = report.nan_residual_found; + const auto& grid = ebosSimulator.gridManager().grid(); + int value = nan_residual_found ? 1 : 0; + + nan_residual_found = grid.comm().max(value); + + if (nan_residual_found) { + for (const auto& well : report.nan_residual_wells) { + OpmLog::debug("NaN residual found with phase " + well.phase_name + " for well " + well.well_name); + } + OPM_THROW(Opm::NumericalProblem, "NaN residual found!"); } } + // checking too large residuals + { + bool too_large_residual_found = report.too_large_residual_found; + const auto& grid = ebosSimulator.gridManager().grid(); + int value = too_large_residual_found ? 1 : 0; + + too_large_residual_found = grid.comm().max(value); + if (too_large_residual_found) { + for (const auto& well : report.too_large_residual_wells) { + OpmLog::debug("Too large residual found with phase " + well.phase_name + " fow well " + well.well_name); + } + OPM_THROW(Opm::NumericalProblem, "Too large residual found!"); + } + } + + // checking convergence + bool converged_well = report.converged; { const auto& grid = ebosSimulator.gridManager().grid(); int value = converged_well ? 1 : 0; converged_well = grid.comm().min(value); } - // TODO: to think about the output here. return converged_well; } diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index aa698a582..6d6c0de8b 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -109,9 +109,40 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const = 0; - virtual bool getWellConvergence(Simulator& ebosSimulator, - const std::vector& B_avg, - const ModelParameters& param) const = 0; + /// a struct to collect information about the convergence checking + struct ConvergenceReport { + struct ProblemWell { + std::string well_name; + std::string phase_name; + }; + bool converged = true; + bool nan_residual_found = false; + std::vector nan_residual_wells; + // We consider Inf is large residual here + bool too_large_residual_found = false; + std::vector too_large_residual_wells; + + ConvergenceReport& operator+=(const ConvergenceReport& rhs) { + converged = converged && rhs.converged; + nan_residual_found = nan_residual_found || rhs.nan_residual_found; + if (rhs.nan_residual_found) { + for (const ProblemWell& well : rhs.nan_residual_wells) { + nan_residual_wells.push_back(well); + } + } + too_large_residual_found = too_large_residual_found || rhs.too_large_residual_found; + if (rhs.too_large_residual_found) { + for (const ProblemWell& well : rhs.too_large_residual_wells) { + too_large_residual_wells.push_back(well); + } + } + return *this; + } + }; + + virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, + const std::vector& B_avg, + const ModelParameters& param) const = 0; virtual void solveEqAndUpdateWellState(const ModelParameters& param, WellState& well_state) = 0; @@ -128,15 +159,15 @@ namespace Opm void computeRepRadiusPerfLength(const Grid& grid, const std::map& cartesian_to_compressed); - // using the solution x to recover the solution xw for wells and applying - // xw to update Well State + /// using the solution x to recover the solution xw for wells and applying + /// xw to update Well State virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param, WellState& well_state) const = 0; - // Ax = Ax - C D^-1 B x + /// Ax = Ax - C D^-1 B x virtual void apply(const BVector& x, BVector& Ax) const = 0; - // r = r - C D^-1 Rw + /// r = r - C D^-1 Rw virtual void apply(BVector& r) const = 0; virtual void computeWellPotentials(const Simulator& ebosSimulator,