From fdb9fa9a864738249f3604e0bf549a4bb4df8728 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Thu, 4 May 2023 12:58:59 -0600
Subject: [PATCH] Add a test.

---
 tests/trilinos/nox_solver_02.cc               | 110 ++++++++++++++++++
 ...lver_02.with_trilinos_with_nox=true.output |  24 ++++
 2 files changed, 134 insertions(+)
 create mode 100644 tests/trilinos/nox_solver_02.cc
 create mode 100644 tests/trilinos/nox_solver_02.with_trilinos_with_nox=true.output

diff --git a/tests/trilinos/nox_solver_02.cc b/tests/trilinos/nox_solver_02.cc
new file mode 100644
index 0000000000..c8365508ae
--- /dev/null
+++ b/tests/trilinos/nox_solver_02.cc
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, 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 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Check TrilinosWrappers::NOXSolver by solving f(x) = atan(x)-0.5 = 0, starting
+// at x=10. This corresponds also to what the sundials/kinsol_06 test checks.
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/trilinos/nox.h>
+
+// as reference solution
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#include "../tests.h"
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  using Number     = double;
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  // set up solver control
+  const unsigned int n_max_iterations  = 100;
+  const double       abs_tolerance     = 1e-9;
+  const double       rel_tolerance     = 1e-5;
+  const double       lin_rel_tolerance = 1e-3;
+
+  TrilinosWrappers::NOXSolver<VectorType>::AdditionalData additional_data(
+    n_max_iterations, abs_tolerance, rel_tolerance);
+
+  // set up parameters
+  Teuchos::RCP<Teuchos::ParameterList> non_linear_parameters =
+    Teuchos::rcp(new Teuchos::ParameterList);
+
+  non_linear_parameters->set("Nonlinear Solver", "Line Search Based");
+  non_linear_parameters->sublist("Printing").set("Output Information", 15);
+  non_linear_parameters->sublist("Direction").set("Method", "Newton");
+  non_linear_parameters->sublist("Direction")
+    .sublist("Newton")
+    .sublist("Linear Solver")
+    .set("Tolerance", lin_rel_tolerance);
+  non_linear_parameters->sublist("Line Search").set("Method", "Polynomial");
+
+
+  // set up solver
+  TrilinosWrappers::NOXSolver<VectorType> solver(additional_data,
+                                                 non_linear_parameters);
+
+  // ... helper functions
+  double J = 0.0;
+
+  solver.residual = [](const auto &src, auto &dst) {
+    // compute residual
+    deallog << "Evaluating residual at u=" << src[0] << std::endl;
+    dst[0] = std::atan(src[0]) - 0.5;
+    return 0;
+  };
+
+  solver.setup_jacobian = [&](const auto &src) {
+    // compute Jacobian
+    deallog << "Setting up Jacobian at u=" << src[0] << std::endl;
+    J = 1. / (1 + src[0] * src[0]);
+    return 0;
+  };
+
+  solver.apply_jacobian = [&](const auto &src, auto &dst) {
+    // solve with Jacobian
+    dst[0] = src[0] * J;
+    return 0;
+  };
+
+  solver.solve_with_jacobian = [&](const auto &src, auto &dst, const auto) {
+    // solve with Jacobian
+    dst[0] = src[0] / J;
+    return 0;
+  };
+
+  // initial guess
+  VectorType solution(1);
+  solution[0] = 10.0;
+
+  // solve with the given initial guess
+  solver.solve(solution);
+
+  deallog << "The solution is: " << solution[0] << std::endl;
+}
diff --git a/tests/trilinos/nox_solver_02.with_trilinos_with_nox=true.output b/tests/trilinos/nox_solver_02.with_trilinos_with_nox=true.output
new file mode 100644
index 0000000000..ea24c462c4
--- /dev/null
+++ b/tests/trilinos/nox_solver_02.with_trilinos_with_nox=true.output
@@ -0,0 +1,24 @@
+
+DEAL::Evaluating residual at u=10.0000
+DEAL::Setting up Jacobian at u=10.0000
+DEAL::Evaluating residual at u=-88.0839
+DEAL::Evaluating residual at u=-7.84235
+DEAL::Evaluating residual at u=8.21577
+DEAL::Setting up Jacobian at u=8.21577
+DEAL::Evaluating residual at u=-56.8358
+DEAL::Evaluating residual at u=-3.24854
+DEAL::Evaluating residual at u=7.06933
+DEAL::Setting up Jacobian at u=7.06933
+DEAL::Evaluating residual at u=-40.3518
+DEAL::Evaluating residual at u=-1.05454
+DEAL::Evaluating residual at u=6.16035
+DEAL::Setting up Jacobian at u=6.16035
+DEAL::Evaluating residual at u=-29.2791
+DEAL::Evaluating residual at u=0.264030
+DEAL::Setting up Jacobian at u=0.264030
+DEAL::Evaluating residual at u=0.522751
+DEAL::Setting up Jacobian at u=0.522751
+DEAL::Evaluating residual at u=0.546074
+DEAL::Setting up Jacobian at u=0.546074
+DEAL::Evaluating residual at u=0.546302
+DEAL::The solution is: 0.546302
-- 
2.39.5