From 8f29cae4cbe245020c0470565bb373a3853d8d46 Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 21 Jan 2015 15:02:11 +0100 Subject: [PATCH 1/3] FullyImplicitBlackoilSolver: throw Opm::NumericalProblem when one of the residuals is NaN. --- .../FullyImplicitBlackoilSolver_impl.hpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 244c04a88..6dc39933f 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1947,9 +1947,9 @@ namespace { RG_sum = RG.sum(); } - const double mass_balance_residual_water = fabs(BW_avg*RW_sum) * dt / pvSum; - const double mass_balance_residual_oil = fabs(BO_avg*RO_sum) * dt / pvSum; - const double mass_balance_residual_gas = fabs(BG_avg*RG_sum) * dt / pvSum; + const double mass_balance_residual_water = std::abs(BW_avg*RW_sum) * dt / pvSum; + const double mass_balance_residual_oil = std::abs(BO_avg*RO_sum) * dt / pvSum; + const double mass_balance_residual_gas = std::abs(BG_avg*RG_sum) * dt / pvSum; bool converged_MB = (mass_balance_residual_water < tol_mb) && (mass_balance_residual_oil < tol_mb) @@ -1963,6 +1963,19 @@ namespace { bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa); bool converged = converged_MB && converged_CNV && converged_Well; + // if one of the residuals is NaN, throw exception, so that the solver can be restarted + if( std::isnan(mass_balance_residual_water) || + std::isnan(mass_balance_residual_oil) || + std::isnan(mass_balance_residual_gas) || + std::isnan(CNVW) || + std::isnan(CNVO) || + std::isnan(CNVG) || + std::isnan(residualWellFlux) || + std::isnan(residualWell) ) + { + OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN"); + } + if (iteration == 0) { std::cout << "\nIter MB(OIL) MB(WATER) MB(GAS) CNVW CNVO CNVG WELL-FLOW WELL-CNTRL\n"; } From 894983fc7ba4bc3b5933a89401d7cd119847a48a Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 21 Jan 2015 15:43:18 +0100 Subject: [PATCH 2/3] FullyImplicitBlackoilSolver: add max_residual_allowed parameter to restart solver when residual is to large. --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 4 +++- .../FullyImplicitBlackoilSolver_impl.hpp | 21 +++++++++++-------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 3faab343e..0d4a39bc1 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -67,6 +67,7 @@ namespace Opm { double relax_max_; double relax_increment_; double relax_rel_tol_; + double max_residual_allowed_; int max_iter_; SolverParameter( const parameter::ParameterGroup& param ); @@ -363,7 +364,8 @@ namespace Opm { double relaxMax() const { return param_.relax_max_; }; double relaxIncrement() const { return param_.relax_increment_; }; double relaxRelTol() const { return param_.relax_rel_tol_; }; - double maxIter() const { return param_.max_iter_; } + double maxIter() const { return param_.max_iter_; } + double maxResidualAllowed() const { return param_.max_residual_allowed_; } }; } // namespace Opm diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 6dc39933f..a4968c7e9 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -41,6 +41,7 @@ #include #include #include +#include //#include // A debugging utility. @@ -169,6 +170,7 @@ namespace { relax_increment_ = 0.1; relax_rel_tol_ = 0.2; max_iter_ = 15; + max_residual_allowed_ = std::numeric_limits< double >::max(); } template @@ -192,6 +194,7 @@ namespace { dr_max_rel_ = param.getDefault("dr_max_rel", dr_max_rel_); relax_max_ = param.getDefault("relax_max", relax_max_); max_iter_ = param.getDefault("max_iter", max_iter_); + max_residual_allowed_ = param.getDefault("max_residual_allowed",max_residual_allowed_); std::string relaxation_type = param.getDefault("relax_type", std::string("dampen")); if (relaxation_type == "dampen") { @@ -1964,16 +1967,16 @@ namespace { bool converged = converged_MB && converged_CNV && converged_Well; // if one of the residuals is NaN, throw exception, so that the solver can be restarted - if( std::isnan(mass_balance_residual_water) || - std::isnan(mass_balance_residual_oil) || - std::isnan(mass_balance_residual_gas) || - std::isnan(CNVW) || - std::isnan(CNVO) || - std::isnan(CNVG) || - std::isnan(residualWellFlux) || - std::isnan(residualWell) ) + if( std::isnan(mass_balance_residual_water) || mass_balance_residual_water > maxResidualAllowed() || + std::isnan(mass_balance_residual_oil) || mass_balance_residual_oil > maxResidualAllowed() || + std::isnan(mass_balance_residual_gas) || mass_balance_residual_gas > maxResidualAllowed() || + std::isnan(CNVW) || CNVW > maxResidualAllowed() || + std::isnan(CNVO) || CNVO > maxResidualAllowed() || + std::isnan(CNVG) || CNVG > maxResidualAllowed() || + std::isnan(residualWellFlux) || residualWellFlux > maxResidualAllowed() || + std::isnan(residualWell) || residualWell > maxResidualAllowed() ) { - OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN"); + OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or to large!"); } if (iteration == 0) { From 3416cba92a481787df6693014764f768587ec886 Mon Sep 17 00:00:00 2001 From: Robert K Date: Mon, 26 Jan 2015 12:37:50 +0100 Subject: [PATCH 3/3] added white space after comma. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index a4968c7e9..4c8fed7b8 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -194,7 +194,7 @@ namespace { dr_max_rel_ = param.getDefault("dr_max_rel", dr_max_rel_); relax_max_ = param.getDefault("relax_max", relax_max_); max_iter_ = param.getDefault("max_iter", max_iter_); - max_residual_allowed_ = param.getDefault("max_residual_allowed",max_residual_allowed_); + max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_); std::string relaxation_type = param.getDefault("relax_type", std::string("dampen")); if (relaxation_type == "dampen") {