Add another test.

diff --git a/tests/trilinos/nox_solver_03.cc b/tests/trilinos/nox_solver_03.cc
new file mode 100644 (file)
index 0000000..a9b381e
--- /dev/null
+++ b/tests/trilinos/nox_solver_03.cc
@@ -0,0 +1,144 @@
+// ---------------------------------------------------------------------
+//
+// 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. This test solves the same problem as
+// sundials/kinsol_05_new_interface, which describes the problem as follows:
+
+// Similar to the _02 test, but we're now actually providing a solver with
+// the Jacobian matrix. For the current case,
+//
+//   F(u) = [  cos(u1 + u2) - 1 + 2*u1  ]
+//          [  sin(u1 - u2)     + 2*u2   ]
+//
+// the Jacobian is the 2x2 matrix
+//
+//   J(u) = [ -sin(u1 + u2) + 2      -sin(u1 + u2)]
+//          [  cos(u1 - u2)          -cos(u1 - u2) + 2]
+//
+// The addition of the +2u_i to the function F does not move the solution
+// (it is still u=0) but it makes sure that the Jacobian at the solution
+// remains non-singular
+
+
+#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 = 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
+  solver.residual = [](const VectorType &u, VectorType &F) -> int {
+    deallog << "Evaluating the solution at u=(" << u[0] << ',' << u[1] << ')'
+            << std::endl;
+
+    F(0) = std::cos(u[0] + u[1]) - 1 + 2 * u[0];
+    F(1) = std::sin(u[0] - u[1]) + 2 * u[1];
+    return 0;
+  };
+
+  FullMatrix<double> J(2, 2);
+  FullMatrix<double> J_inverse(2, 2);
+
+  solver.setup_jacobian = [&J, &J_inverse](const VectorType &u) -> int {
+    // We don't do any kind of set-up in this program, but we can at least
+    // say that we're here
+    deallog << "Setting up Jacobian system at u=(" << u[0] << ',' << u[1] << ')'
+            << std::endl;
+
+    J(0, 0) = -std::sin(u[0] + u[1]) + 2;
+    J(0, 1) = -std::sin(u[0] + u[1]);
+    J(1, 0) = std::cos(u[0] - u[1]);
+    J(1, 1) = -std::cos(u[0] - u[1]) + 2;
+
+    J_inverse.invert(J);
+
+    return 0;
+  };
+
+  solver.apply_jacobian = [&](const VectorType &src, VectorType &dst) {
+    J.vmult(dst, src);
+    return -0;
+  };
+
+  solver.solve_with_jacobian = [&J_inverse](const VectorType &rhs,
+                                            VectorType &      dst,
+                                            const double /*tolerance*/) -> int {
+    deallog << "Solving Jacobian system with rhs=(" << rhs[0] << ',' << rhs[1]
+            << ')' << std::endl;
+
+    J_inverse.vmult(dst, rhs);
+
+    return 0;
+  };
+
+  // initial guess
+  const unsigned int N = 2;
+  VectorType         solution(N);
+  solution[0] = 0.5;
+  solution[1] = 1.234;
+
+  // solve with the given initial guess
+  auto niter = solver.solve(solution);
+  solution.print(deallog.get_file_stream());
+  deallog << "Converged in " << niter << " iterations." << std::endl;
+}

diff --git a/tests/trilinos/nox_solver_03.with_trilinos_with_nox=true.output b/tests/trilinos/nox_solver_03.with_trilinos_with_nox=true.output
new file mode 100644 (file)
index 0000000..6a11606
--- /dev/null
+++ b/tests/trilinos/nox_solver_03.with_trilinos_with_nox=true.output
@@ -0,0 +1,16 @@
+
+DEAL::Evaluating the solution at u=(0.500000,1.23400)
+DEAL::Setting up Jacobian system at u=(0.500000,1.23400)
+DEAL::Solving Jacobian system with rhs=(-0.162480,1.79816)
+DEAL::Evaluating the solution at u=(-0.282294,0.265967)
+DEAL::Setting up Jacobian system at u=(-0.282294,0.265967)
+DEAL::Solving Jacobian system with rhs=(-0.564722,0.0107297)
+DEAL::Evaluating the solution at u=(-0.000445202,0.0468183)
+DEAL::Setting up Jacobian system at u=(-0.000445202,0.0468183)
+DEAL::Solving Jacobian system with rhs=(-0.00196544,0.0463907)
+DEAL::Evaluating the solution at u=(-0.000536539,0.000570488)
+DEAL::Setting up Jacobian system at u=(-0.000536539,0.000570488)
+DEAL::Solving Jacobian system with rhs=(-0.00107308,3.39491e-05)
+DEAL::Evaluating the solution at u=(-2.88123e-10,7.40347e-10)
+-2.881e-10 7.403e-10 
+DEAL::Converged in 4 iterations.