Closes #313. Return success flag from GSI solver

src/gsi/GasSurfaceInteraction.cpp

@@ -196,9 +196,9 @@ void GasSurfaceInteraction::setGasRadHeatFlux(
 
 //==============================================================================
 
-void GasSurfaceInteraction::solveSurfaceBalance()
+bool GasSurfaceInteraction::solveSurfaceBalance()
 {
-    mp_surf->solveSurfaceBalance();
+    return mp_surf->solveSurfaceBalance();
 }
 
 //==============================================================================

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

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 "

src/gsi/SurfaceBalanceSolverMass.cpp

@@ -158,7 +158,7 @@ class SurfaceBalanceSolverMass :
 
 //=============================================================================
 
-    void solveSurfaceBalance()
+    bool solveSurfaceBalance()
     {
         // errorUninitializedDiffusionModel
         errorSurfaceStateNotSet();
@@ -173,14 +173,16 @@ 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);
 
         // Setting the state again
         m_surf_state.setSurfaceState(
             mv_rhoi.data(), mv_Tsurf.data(), set_state_with_rhoi_T);
+
+        return converged;
     }
 
 //==============================================================================

src/gsi/SurfaceBalanceSolverMassEnergy.cpp

@@ -195,7 +195,7 @@ class SurfaceBalanceSolverMassEnergy :
 
 //=============================================================================
 
-    void solveSurfaceBalance()
+    bool solveSurfaceBalance()
     {
         // errorUninitializedDiffusionModel
         errorSurfaceStateNotSet();
@@ -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(
@@ -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;
     }
 
 //==============================================================================

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 <typename T, typename Solver>
@@ -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<T, Solver>::NewtonSolver()
 //==============================================================================
 
 template <typename T, typename Solver>
-T& NewtonSolver<T, Solver>::solve(T& x)
+bool NewtonSolver<T, Solver>::solve(T& x)
 {
     using std::cout;
     using std::endl;
@@ -159,11 +159,14 @@ T& NewtonSolver<T, Solver>::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;
 }
 
 //==============================================================================

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

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();
 

tests/c++/test_newton_solver.cpp

@@ -0,0 +1,127 @@
+/*
+ * 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
+ * <http://www.gnu.org/licenses/>.
+ */
+
+#include <iostream>
+
+#include "mutation++.h"
+#include "Configuration.h"
+#include "TestMacros.h"
+#include "NewtonSolver.h"
+#include <catch.hpp>
+#include <Eigen/Dense>
+
+//==============================================================================
+
+/// 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<Eigen::Infinity>();
+    }
+
+//==============================================================================
+
+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
+    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));
+    CHECK(x(1) == Catch::Detail::Approx(1.0).margin(tol));
+}
+
+//==============================================================================
+