diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 313c51e58..740b85679 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -39,6 +39,7 @@ list (APPEND MAIN_SOURCE_FILES src/opm/common/utility/parameters/ParameterGroup.cpp src/opm/common/utility/parameters/ParameterTools.cpp src/opm/common/utility/numeric/calculateCellVol.cpp + src/opm/common/utility/numeric/RootFinders.cpp src/opm/common/utility/shmatch.cpp src/opm/common/utility/TimeService.cpp src/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.cpp diff --git a/opm/common/utility/numeric/RootFinders.hpp b/opm/common/utility/numeric/RootFinders.hpp index de597307f..3f2f17839 100644 --- a/opm/common/utility/numeric/RootFinders.hpp +++ b/opm/common/utility/numeric/RootFinders.hpp @@ -22,66 +22,42 @@ #define OPM_ROOTFINDERS_HEADER #include -#include -#include -#include -#include #include -#include +#include +#include namespace Opm { struct ThrowOnError { - static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) - { - std::ostringstream sstr; - sstr << "Error in parameters, zero not bracketed: [a, b] = [" - << x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1; - OpmLog::debug(sstr.str()); - OPM_THROW_NOLOG(std::runtime_error, sstr.str()); - return -1e100; // Never reached. - } - static double handleTooManyIterations(const double x0, const double x1, const int maxiter) - { - OPM_THROW(std::runtime_error, "Maximum number of iterations exceeded: " << maxiter << "\n" - << "Current interval is [" << std::min(x0, x1) << ", " - << std::max(x0, x1) << "] abs(x0-x1) " << std::abs(x0-x1)); - return -1e100; // Never reached. - } + static double handleBracketingFailure(const double x0, const double x1, + const double f0, const double f1); + + static double handleTooManyIterations(const double x0, + const double x1, const int maxiter); }; struct WarnAndContinueOnError { - static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) - { - OPM_REPORT; - std::cerr << "Error in parameters, zero not bracketed: [a, b] = [" - << x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1 - << ""; - return std::fabs(f0) < std::fabs(f1) ? x0 : x1; - } - static double handleTooManyIterations(const double x0, const double x1, const int maxiter) - { - OPM_REPORT; - std::cerr << "Maximum number of iterations exceeded: " << maxiter - << ", current interval is [" << std::min(x0, x1) << ", " - << std::max(x0, x1) << "] abs(x0-x1) " << std::abs(x0-x1); - return 0.5*(x0 + x1); - } + static double handleBracketingFailure(const double x0, const double x1, + const double f0, const double f1); + static double handleTooManyIterations(const double x0, + const double x1, const int maxiter); }; struct ContinueOnError { - static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1) + static double handleBracketingFailure(const double x0, const double x1, + const double f0, const double f1) { return std::fabs(f0) < std::fabs(f1) ? x0 : x1; } - static double handleTooManyIterations(const double x0, const double x1, const int /*maxiter*/) + static double handleTooManyIterations(const double x0, + const double x1, const int /*maxiter*/) { return 0.5*(x0 + x1); } @@ -427,7 +403,9 @@ namespace Opm cur_dx = -2.0*cur_dx; } if (i == max_iters) { - OPM_THROW(std::runtime_error, "Could not bracket zero in " << max_iters << "iterations."); + OPM_THROW(std::runtime_error, + "Could not bracket zero in " + + std::to_string(max_iters) + " iterations."); } if (cur_dx < 0.0) { a = x0 + cur_dx; @@ -438,13 +416,6 @@ namespace Opm } } - - - - } // namespace Opm - - - #endif // OPM_ROOTFINDERS_HEADER diff --git a/src/opm/common/utility/numeric/RootFinders.cpp b/src/opm/common/utility/numeric/RootFinders.cpp new file mode 100644 index 000000000..b644039c9 --- /dev/null +++ b/src/opm/common/utility/numeric/RootFinders.cpp @@ -0,0 +1,75 @@ +/* + Copyright 2010, 2019 SINTEF Digital + Copyright 2010, 2019 Equinor ASA + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include + +#include +#include + +#include + +#include + +namespace Opm +{ + +double ThrowOnError::handleBracketingFailure(const double x0, const double x1, + const double f0, const double f1) +{ + const auto str = fmt::format("Error in parameters, zero not bracketed: " + "[a, b] = [{}, {}] f(a) = {} f(b) = {}", + x0, x1, f0, f1); + OpmLog::debug(str); + OPM_THROW_NOLOG(std::runtime_error, str); + return -1e100; // Never reached. +} + +double ThrowOnError::handleTooManyIterations(const double x0, + const double x1, const int maxiter) +{ + OPM_THROW(std::runtime_error, + fmt::format("Maximum number of iterations exceeded: {}\n" + "Current interval is [{}, {}] abs(x0-91) {}", + maxiter, std::min(x0, x1), std::max(x0, x1), std::abs(x0-x1))); + return -1e100; // Never reached. +} + +double WarnAndContinueOnError::handleBracketingFailure(const double x0, const double x1, + const double f0, const double f1) +{ + OPM_REPORT; + OpmLog::warning(fmt::format("Error in parameters, zero not bracketed: " + "[a, b] = [{}, {}] f(a) = {} f(b) = {}", + x0, x1, f0, f1)); + return std::fabs(f0) < std::fabs(f1) ? x0 : x1; +} + +double WarnAndContinueOnError::handleTooManyIterations(const double x0, + const double x1, const int maxiter) +{ + OPM_REPORT; + OpmLog::warning(fmt::format("Maximum number of iterations exceeded: {}, " + "current interval is [{}, {}] abs(x0-x1) {}", + maxiter, std::abs(x0-x1))); + return 0.5*(x0 + x1); +} + +} // namespace Opm