From 8c994ea8c0875c37d7f6a1801697922f26e0b991 Mon Sep 17 00:00:00 2001 From: Matthew Goodson Date: Fri, 19 Nov 2021 11:47:22 -0500 Subject: [PATCH 1/3] test: #313: Add a unit test for the Newton solver class. Add a unit test for the Newton solver class. Uses Rosenbrock function. See #313 --- tests/c++/test_newton_solver.cpp | 125 +++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 tests/c++/test_newton_solver.cpp diff --git a/tests/c++/test_newton_solver.cpp b/tests/c++/test_newton_solver.cpp new file mode 100644 index 00000000..c630c43a --- /dev/null +++ b/tests/c++/test_newton_solver.cpp @@ -0,0 +1,125 @@ +/* + * Copyright 2015-2020 von Karman Institute for Fluid Dynamics (VKI) + * + * This file is part of MUlticomponent Thermodynamic And Transport + * properties for IONized gases in C++ (Mutation++) software package. + * + * Mutation++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * Mutation++ 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with Mutation++. If not, see + * . + */ + +#include + +#include "mutation++.h" +#include "Configuration.h" +#include "TestMacros.h" +#include "NewtonSolver.h" +#include +#include + +//============================================================================== + +/// Helper class for Newton solver tests +class NewtonSolverTest : + public Mutation::Numerics::NewtonSolver< + Eigen::VectorXd, NewtonSolverTest> +{ +public: + NewtonSolverTest() + : mv_X(2), + mv_dX(2), + mv_f(2), + mv_f_old(2), + mv_df(2), + mv_trans(1, 2), + m_jac(2, 2), + m_inv_jac(2, 2), + m_tol(1.e-12) + { + // Setup NewtonSolver + setMaxIterations(10); + setWriteConvergenceHistory(true); + setEpsilon(m_tol); + } + + void updateFunction(Eigen::VectorXd& v_X) { + // Save the old f + mv_f_old = mv_f; + + // g(0) = 10.0*(x2 - x1^2) + // g(1) = 1.0 - x1 + mv_f(0) = 10.0*(v_X(1) - v_X(0)*v_X(0)); + mv_f(1) = 1.0 - v_X(0); + + // Compute and store difference from new to old + mv_df = mv_f - mv_f_old; + } + + void updateJacobian(Eigen::VectorXd& v_X) + { + m_jac(0, 0) = -20.0*v_X(0); + m_jac(0, 1) = 10.0; + m_jac(1, 0) = -1.0; + m_jac(1, 1) = 0.0; + } + + Eigen::VectorXd& systemSolution() + { + mv_dX = m_jac.fullPivLu().solve(mv_f); + return mv_dX; + } + + double norm() { + return mv_dX.lpNorm(); + } + +//============================================================================== + +private: + Eigen::VectorXd mv_X; + Eigen::VectorXd mv_dX; + Eigen::VectorXd mv_f; + Eigen::VectorXd mv_f_old; + Eigen::VectorXd mv_df; + Eigen::MatrixXd mv_trans; + Eigen::MatrixXd m_jac; + Eigen::MatrixXd m_inv_jac; + double m_tol; +}; + +//============================================================================== + +/** + * Checks that Newton's method solver works by solving the Rosenbrock function + */ +TEST_CASE("Newton solver returns the correct root", + "[numerics][NewtonSolver]" +) +{ + constexpr double tol = 1.0e-12; + NewtonSolverTest test; + Eigen::VectorXd x(2); + // Set the initial guess + x(0) = 10.0; + x(1) = 10.0; + // Solve + test.solve(x); + // Make sure the roots are 1.0 + INFO("x:\n" << x); + CHECK(x(0) == Catch::Detail::Approx(1.0).margin(tol)); + CHECK(x(1) == Catch::Detail::Approx(1.0).margin(tol)); +} + +//============================================================================== + From 2e0e620503b2deccb81fcb6a9c667dd8d1a0f89d Mon Sep 17 00:00:00 2001 From: Matthew Goodson Date: Thu, 9 Jul 2026 12:34:59 -0400 Subject: [PATCH 2/3] feat: #313: Return success or failure from NewtonSolver Update the return flag from the NewtonSolver's `solve` function to return a success/failure flag based on convergence, rather than returning a reference to `x`, which is redundant because `x` was already passed in by non-const reference. See #313. --- src/gsi/SurfaceBalanceSolverMass.cpp | 2 +- src/gsi/SurfaceBalanceSolverMassEnergy.cpp | 2 +- src/numerics/NewtonSolver.h | 17 ++++++++++------- tests/c++/test_newton_solver.cpp | 4 +++- 4 files changed, 15 insertions(+), 10 deletions(-) diff --git a/src/gsi/SurfaceBalanceSolverMass.cpp b/src/gsi/SurfaceBalanceSolverMass.cpp index d08681f6..371ee686 100644 --- a/src/gsi/SurfaceBalanceSolverMass.cpp +++ b/src/gsi/SurfaceBalanceSolverMass.cpp @@ -173,7 +173,7 @@ class SurfaceBalanceSolverMass : computeMoleFracfromPartialDens(mv_rhoi, mv_X); // Solving - mv_X = solve(mv_X); + const bool converged = solve(mv_X); applyTolerance(mv_X); computePartialDensfromMoleFrac(mv_X, mv_rhoi); diff --git a/src/gsi/SurfaceBalanceSolverMassEnergy.cpp b/src/gsi/SurfaceBalanceSolverMassEnergy.cpp index 1a7bd827..6648eb1c 100644 --- a/src/gsi/SurfaceBalanceSolverMassEnergy.cpp +++ b/src/gsi/SurfaceBalanceSolverMassEnergy.cpp @@ -215,7 +215,7 @@ class SurfaceBalanceSolverMassEnergy : applyTolerance(mv_X); // Solving - mv_X = solve(mv_X); + const bool converged = solve(mv_X); applyTolerance(mv_X); computePartialDensfromMoleFrac( diff --git a/src/numerics/NewtonSolver.h b/src/numerics/NewtonSolver.h index ef82a93b..1bd9defe 100644 --- a/src/numerics/NewtonSolver.h +++ b/src/numerics/NewtonSolver.h @@ -40,7 +40,7 @@ namespace Mutation { /** * Implements the skeleton framework for Newton's method for solving a nonlinear - * system of equations. Classes should extend this class and provide the + * system of equations. Classes should extend this class and provide the * methods for computing f(x), J(x) = df/dx, inv(J)*f, and norm(f). */ template @@ -61,7 +61,7 @@ class NewtonSolver /** * Uses Newton's method to compute the zero of f(x) given an initial guess. */ - T& solve(T& x); + bool solve(T& x); /** * Set the residual norm tolerance. @@ -116,7 +116,7 @@ NewtonSolver::NewtonSolver() //============================================================================== template -T& NewtonSolver::solve(T& x) +bool NewtonSolver::solve(T& x) { using std::cout; using std::endl; @@ -159,11 +159,14 @@ T& NewtonSolver::solve(T& x) cout << ", relative residual = " << resnorm << endl; } - if (resnorm > m_epsilon && m_conv_hist) { - cout << "Newton failed to converge after " << m_max_iter - << " iterations with a relative residual of " << resnorm << endl; + if (resnorm > m_epsilon) { + if (m_conv_hist) + cout << "Newton failed to converge after " << m_max_iter + << " iterations with a relative residual of " << resnorm << endl; + return false; + } else { + return true; } - return x; } //============================================================================== diff --git a/tests/c++/test_newton_solver.cpp b/tests/c++/test_newton_solver.cpp index c630c43a..e055786d 100644 --- a/tests/c++/test_newton_solver.cpp +++ b/tests/c++/test_newton_solver.cpp @@ -114,7 +114,9 @@ TEST_CASE("Newton solver returns the correct root", x(0) = 10.0; x(1) = 10.0; // Solve - test.solve(x); + bool converged = test.solve(x); + // Make sure it converged + CHECK(converged); // Make sure the roots are 1.0 INFO("x:\n" << x); CHECK(x(0) == Catch::Detail::Approx(1.0).margin(tol)); From 824e92180d731c746d432f1962bdf7475a6f7ba5 Mon Sep 17 00:00:00 2001 From: Matthew Goodson Date: Thu, 9 Jul 2026 12:36:21 -0400 Subject: [PATCH 3/3] feat: #313: Return success flag from GSI solver. Return a boolean indicating success/failure from the GSI `solveSurfaceBalance` routine. Previously was void. See #313. --- src/gsi/GasSurfaceInteraction.cpp | 4 ++-- src/gsi/GasSurfaceInteraction.h | 4 +++- src/gsi/Surface.h | 2 +- src/gsi/SurfaceBalanceSolverMass.cpp | 4 +++- src/gsi/SurfaceBalanceSolverMassEnergy.cpp | 4 +++- tests/c++/test_gsi_mass.cpp | 3 ++- tests/c++/test_gsi_mass_energy.cpp | 12 ++++++++---- 7 files changed, 22 insertions(+), 11 deletions(-) diff --git a/src/gsi/GasSurfaceInteraction.cpp b/src/gsi/GasSurfaceInteraction.cpp index 26685aac..2ec8334e 100644 --- a/src/gsi/GasSurfaceInteraction.cpp +++ b/src/gsi/GasSurfaceInteraction.cpp @@ -196,9 +196,9 @@ void GasSurfaceInteraction::setGasRadHeatFlux( //============================================================================== -void GasSurfaceInteraction::solveSurfaceBalance() +bool GasSurfaceInteraction::solveSurfaceBalance() { - mp_surf->solveSurfaceBalance(); + return mp_surf->solveSurfaceBalance(); } //============================================================================== diff --git a/src/gsi/GasSurfaceInteraction.h b/src/gsi/GasSurfaceInteraction.h index 14df0bcb..2888b565 100644 --- a/src/gsi/GasSurfaceInteraction.h +++ b/src/gsi/GasSurfaceInteraction.h @@ -163,8 +163,10 @@ class GasSurfaceInteraction * Function to be called in order to solve the mass and energy balances at * the surface according to the input model. The output state is stored in the * surface state and can be accessed by the getSurfaceState function. + + * @returns Boolean flag indicating whether solver converged */ - void solveSurfaceBalance(); + bool solveSurfaceBalance(); /** * Function which allows to change the number of iterations when solving diff --git a/src/gsi/Surface.h b/src/gsi/Surface.h index f9481cd1..2e16e0ce 100644 --- a/src/gsi/Surface.h +++ b/src/gsi/Surface.h @@ -165,7 +165,7 @@ class Surface * Virtual function to be called in order to solve the * surface balance. */ - virtual void solveSurfaceBalance() + virtual bool solveSurfaceBalance() { throw LogicError() << "solveSurfaceBalance can be called only when solving " diff --git a/src/gsi/SurfaceBalanceSolverMass.cpp b/src/gsi/SurfaceBalanceSolverMass.cpp index 371ee686..3de23c6b 100644 --- a/src/gsi/SurfaceBalanceSolverMass.cpp +++ b/src/gsi/SurfaceBalanceSolverMass.cpp @@ -158,7 +158,7 @@ class SurfaceBalanceSolverMass : //============================================================================= - void solveSurfaceBalance() + bool solveSurfaceBalance() { // errorUninitializedDiffusionModel errorSurfaceStateNotSet(); @@ -181,6 +181,8 @@ class SurfaceBalanceSolverMass : // Setting the state again m_surf_state.setSurfaceState( mv_rhoi.data(), mv_Tsurf.data(), set_state_with_rhoi_T); + + return converged; } //============================================================================== diff --git a/src/gsi/SurfaceBalanceSolverMassEnergy.cpp b/src/gsi/SurfaceBalanceSolverMassEnergy.cpp index 6648eb1c..21c09bca 100644 --- a/src/gsi/SurfaceBalanceSolverMassEnergy.cpp +++ b/src/gsi/SurfaceBalanceSolverMassEnergy.cpp @@ -195,7 +195,7 @@ class SurfaceBalanceSolverMassEnergy : //============================================================================= - void solveSurfaceBalance() + bool solveSurfaceBalance() { // errorUninitializedDiffusionModel errorSurfaceStateNotSet(); @@ -224,6 +224,8 @@ class SurfaceBalanceSolverMassEnergy : // Setting the state again m_surf_state.setSurfaceState( mv_rhoi.data(), mv_X.tail(m_nT).data(), set_state_with_rhoi_T); + + return converged; } //============================================================================== diff --git a/tests/c++/test_gsi_mass.cpp b/tests/c++/test_gsi_mass.cpp index 52b5c49e..88e91048 100644 --- a/tests/c++/test_gsi_mass.cpp +++ b/tests/c++/test_gsi_mass.cpp @@ -65,7 +65,8 @@ TEST_CASE("Solution of the MassBalanceSolver is converged.", mix.setSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); mix.setDiffusionModel(xi_e.data(), dx); - mix.solveSurfaceBalance(); + bool converged = mix.solveSurfaceBalance(); + CHECK(converged); mix.getSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); // Verifying the solution gives low residual in the balance equations diff --git a/tests/c++/test_gsi_mass_energy.cpp b/tests/c++/test_gsi_mass_energy.cpp index c0d9c849..2405c0f9 100644 --- a/tests/c++/test_gsi_mass_energy.cpp +++ b/tests/c++/test_gsi_mass_energy.cpp @@ -74,7 +74,8 @@ TEST_CASE("Solution of the MassEnergyBalanceSolver is converged.", "[gsi]") mix.setSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); // Solve balance and request solution - mix.solveSurfaceBalance(); + bool converged = mix.solveSurfaceBalance(); + CHECK(converged); mix.getSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); // Verifying the solution gives low residual in the balance equations @@ -156,7 +157,8 @@ TEST_CASE("Solution of the MassEnergyBalanceSolver is converged.", "[gsi]") mix.setSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); // Solve balance and request solution - mix.solveSurfaceBalance(); + bool converged = mix.solveSurfaceBalance(); + CHECK(converged); mix.getSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); // Verifying the solution gives low residual in the balance equations @@ -248,7 +250,8 @@ TEST_CASE("Solution of the MassEnergyBalanceSolver is converged.", "[gsi]") mix.setSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); // Solve balance and request solution - mix.solveSurfaceBalance(); + bool converged = mix.solveSurfaceBalance(); + CHECK(converged); mix.getSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); double rho = rhoi_s.sum(); @@ -348,7 +351,8 @@ TEST_CASE("Solution of the MassEnergyBalanceSolver is converged.", "[gsi]") mix.setSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); // Solve balance and request solution - mix.solveSurfaceBalance(); + bool converged = mix.solveSurfaceBalance(); + CHECK(converged); mix.getSurfaceState(rhoi_s.data(), T_s.data(), set_state_with_rhoi_T); double rho = rhoi_s.sum();