]> https://gitweb.dealii.org/ - dealii.git/commitdiff
NonlinearSolverSelector class implementation
authorSean Ingimarson <singima@clemson.edu>
Wed, 5 Apr 2023 12:03:36 +0000 (08:03 -0400)
committerSean Ingimarson <singima@clemson.edu>
Sat, 13 May 2023 21:12:27 +0000 (17:12 -0400)
doc/news/changes/minor/20230501Ingimarson [new file with mode: 0644]
include/deal.II/numerics/nonlinear.h [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_01.cc [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_01.with_sundials=on.output [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_02.cc [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_02.with_trilinos=on.output [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_03.cc [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_03.with_sundials=on.mpirun=4.output [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_04.cc [new file with mode: 0644]
tests/numerics/nonlinear_solver_selector_04.with_trilinos=on.mpirun=4.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20230501Ingimarson b/doc/news/changes/minor/20230501Ingimarson
new file mode 100644 (file)
index 0000000..0da5689
--- /dev/null
@@ -0,0 +1,3 @@
+New: Added new nonlinear solver class NonlinearSolverSelector.
+<br>
+(Sean Ingimarson, 2023/05/01)
diff --git a/include/deal.II/numerics/nonlinear.h b/include/deal.II/numerics/nonlinear.h
new file mode 100644 (file)
index 0000000..11bb956
--- /dev/null
@@ -0,0 +1,558 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/exceptions.h>
+
+#include <deal.II/sundials/kinsol.h>
+
+#include <deal.II/trilinos/nox.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Selects a nonlinear solver by choosing between KINSOL or NOX.
+ * KINSOL and NOX are nonlinear solvers included in the SUNDIALS
+ * package and the Trilinos package, respectively. If no solver type is
+ * specified it will automaticlaly choose a solver based on what
+ * deal.II was configured with, KINSOL having priority.
+ *
+ * By calling the @p solve function of this @p NonlinearSolverSelector, it selects the
+ * @p solve function of that @p Solver that was specified in the constructor
+ * of this class, similar to the SolverSelector class.
+ *
+ * <h3>Usage</h3>
+ * An example of code one would run with this class is the
+ * following:
+ * @code
+ * // Generate a @p NonlinearSolverSelector that uses @p KINSOL
+ * NonlinearSolverSelector<Vector<double>>::AdditionalData additional_data;
+ * additional_data.solver_type =
+ *  NonlinearSolverSelector<Vector<double>>::AdditionalData::SolverType::kinsol;
+ *
+ * NonlinearSolverSelector<Vector<double>> nonlinear_solver(additional_data);
+ *
+ * // Functions that are required for solving a nonlinear problem,
+ * // which are utilized in both @p KINSOL and @p NOX.
+ * nonlinear_solver.reinit_vector = [&](Vector<double> &x) {...}
+ *
+ * nonlinear_solver.residual =
+ *           [&](const Vector<double> &evaluation_point,
+ *               Vector<double> &      residual) {...}
+ *
+ * nonlinear_solver.setup_jacobian =
+ *           [&](const Vector<double> &current_u,
+ *               const Vector<double> &current_f) {...}
+ *
+ * nonlinear_solver.solve_with_jacobian =
+ *           [&](const Vector<double> &rhs,
+ *               Vector<double> &      dst,
+ *               const double tolerance) {...}
+ *
+ * // Calling the @p solve function with an initial guess.
+ * nonlinear_solver.solve(current_solution);
+ * @endcode
+ */
+template <typename VectorType = Vector<double>>
+class NonlinearSolverSelector
+{
+public:
+  /**
+   * Additional data that will be sent through the NonlinearSolverSelector
+   * class, then into the specific nonlinear solver classes themselves.
+   * Many of these parameters can also be foud in the documentation for
+   * @p KINSOL and @p NOX.
+   */
+  class AdditionalData
+  {
+  public:
+    /**
+     * NonlinearSolverSelector solution strategy. The solver types included in
+     * this class both use a Newton-Krylov solver with a line search and even
+     * Picard iterations (KINSOL).
+     */
+    enum SolutionStrategy
+    {
+      /**
+       * Standard Newton iteration.
+       */
+      newton,
+      /**
+       * Newton iteration with linesearch.
+       */
+      linesearch,
+      /**
+       * Picard iteration.
+       */
+      picard,
+    };
+
+    enum SolverType
+    {
+      /**
+       * Default parameter, will use whatever solver is available
+       * with KINSOL as the priority.
+       */
+      automatic,
+      /**
+       * KINSOL nonlinear solver, part of the SUNDIALS package.
+       */
+      kinsol,
+      /**
+       * NOX nonlinear solver, part of the TRILINOS package.
+       */
+      nox,
+    };
+
+    /**
+     * Initialization parameters for NonlinearSolverSelector.
+     *
+     * @param solver_type Nonlinear solver type, can be 'auto', 'kinsol', or 'nox'.
+     * @param strategy Method of solving nonlinear problem, can be 'newton',
+     * 'linesearch', or 'picard'.
+     * @param maximum_non_linear_iterations Maximum number of nonlinear
+     * iterations. This parameter is shared between KINSOL and NOX.
+     * @param function_tolerance Function norm stopping tolerance.
+     * @param relative_tolerance Relative function norm stopping tolerance.
+     * @param step_tolerance Step tolerance for solution update.
+     * @param anderson_subspace_size Anderson acceleration subspace size
+     */
+    AdditionalData(const SolverType &      solver_type = automatic,
+                   const SolutionStrategy &strategy    = linesearch,
+                   const unsigned int      maximum_non_linear_iterations = 200,
+                   const double            function_tolerance            = 1e-8,
+                   const double            relative_tolerance            = 1e-5,
+                   const double            step_tolerance                = 0.0,
+                   const unsigned int      anderson_subspace_size        = 0);
+
+    /**
+     * The type of nonlinear solver to use. The default value is set to 'auto',
+     * which will use either KINSOL or NOX depending on which package is
+     * installed, KINSOL having priority.
+     */
+    SolverType solver_type;
+
+    /**
+     * The solution strategy to use. For this class, you can choose from
+     * SolutionStrategy::newton, SolutionStrategy::linesearch, or
+     * SolutionStrategy::picard. More details on this can be found on the
+     * @p KINSOL documentation.
+     */
+    SolutionStrategy strategy;
+
+    /**
+     * Maximum number of nonlinear iterations allowed.
+     */
+    unsigned int maximum_non_linear_iterations;
+
+    /**
+     * A scalar used as a stopping tolerance on the scaled
+     * maximum norm of the system function $F(u)$ or $G(u)$.
+     *
+     * If set to zero, default values provided by KINSOL will be used.
+     */
+    double function_tolerance;
+
+    /**
+     * A scalar used as a stopping tolerance on the minimum
+     * scaled step length.
+     *
+     * If set to zero, default values provided by KINSOL will be used.
+     */
+    double step_tolerance = 0.0;
+
+    /**
+     * Relative $l_2$ tolerance of the residual to be reached.
+     *
+     * @note Solver terminates successfully if either the absolute or
+     * the relative tolerance has been reached.
+     */
+    const double relative_tolerance;
+
+    /**
+     * The size of the subspace used with Anderson acceleration
+     * in conjunction with Picard or fixed-point iteration.
+     *
+     * If you set this to 0, no acceleration is used.
+     */
+    unsigned int anderson_subspace_size;
+  };
+
+  /**
+   * Constructor, filling in default values
+   */
+  NonlinearSolverSelector();
+
+  /**
+   * Constructor, selecting the solver and other parametersspecified in
+   * @p additional_data.
+   */
+  NonlinearSolverSelector(const AdditionalData &additional_data);
+
+  /**
+   * Constructor
+   *
+   * @param additional_data NonlinearSolverSelector configuration data
+   * @param mpi_communicator MPI communicator over which logging operations are
+   * computer.
+   */
+  NonlinearSolverSelector(const AdditionalData &additional_data,
+                          const MPI_Comm &      mpi_communicator);
+
+  /**
+   * Select a new nonlinear solver. All solver names used in this class are
+   * all lower case.
+   */
+  void
+  select(const typename AdditionalData::SolverType &type);
+
+/**
+ * Set the additional data. For more information see the @p Solver class.
+ */
+#ifdef DEAL_II_WITH_TRILINOS
+  void
+  set_data(
+    const typename TrilinosWrappers::NOXSolver<VectorType>::AdditionalData
+      &                                         additional_data,
+    const Teuchos::RCP<Teuchos::ParameterList> &parameters =
+      Teuchos::rcp(new Teuchos::ParameterList));
+#endif
+
+/**
+ * Set the additional data. For more information see the @p Solver class.
+ */
+#ifdef DEAL_II_WITH_SUNDIALS
+  void
+  set_data(const typename SUNDIALS::KINSOL<VectorType>::AdditionalData
+             &additional_data);
+#endif
+
+  /**
+   * Solve the nonlinear system. KINSOL uses the content of
+   * `initial_guess_and_solution` as an initial guess, and
+   * stores the final solution in the same vector.
+   *
+   * The functions herein are nearly identical in setup to what can be found
+   * in the KINSOL and NOX documentation.
+   */
+  void
+  solve(VectorType &initial_guess_and_solution);
+
+  /**
+   * A function object that users need to supply and that is intended to
+   * reinitize the given vector to its correct size, and block structure (if
+   * block vectors are used), along with any
+   * other properties necessary.
+   */
+  std::function<void(VectorType &)> reinit_vector;
+
+  /**
+   * A function object that users should supply and that is intended to
+   * compute the residual `dst = F(src)`.
+   *
+   * This function should return an int for either failure or success.
+   */
+  std::function<int(const VectorType &src, VectorType &dst)> residual;
+
+  /**
+   * A function object that users may supply and that is intended to
+   * prepare the linear solver for subsequent calls to
+   * solve_jacobian_system().
+   *
+   * The job of setup_jacobian() is to prepare the linear solver for
+   * subsequent calls to solve_with_jacobian(), in the solution of linear
+   * systems $Ax = b$. The exact nature of this system depends on the
+   * SolutionStrategy that has been selected.
+   *
+   * In the cases strategy = SolutionStrategy::newton or
+   * SolutionStrategy::linesearch, $A$ is the Jacobian $J = \partial
+   * F/\partial u$. If strategy = SolutionStrategy::picard, $A$ is the
+   * approximate Jacobian matrix $L$.
+   *
+   * The setup_jacobian() function may call a user-supplied function, or a
+   * function within the linear solver module, to compute Jacobian-related
+   * data that is required by the linear solver. It may also preprocess that
+   * data as needed for solve_with_jacobian(), which may involve calling a
+   * generic function (such as for LU factorization) or, more generally,
+   * build preconditioners from the assembled Jacobian. In any case, the
+   * data so generated may then be used whenever a linear system is solved.
+   *
+   * @param current_u Current value of $u$
+   */
+  std::function<int(const VectorType &current_u)> setup_jacobian;
+
+  /**
+   * A function object that users may supply and that is intended to solve
+   * a linear system with the Jacobian matrix.
+   *
+   * Specific details on this function can be found in the KINSOL
+   *
+   * Arguments to the function are:
+   *
+   * @param[in] rhs The system right hand side to solve for.
+   * @param[out] dst The solution of $J^{-1} * src$.
+   * @param[in] tolerance The tolerance with which to solve the linear system
+   *   of equations.
+   */
+  std::function<
+    int(const VectorType &rhs, VectorType &dst, const double tolerance)>
+    solve_with_jacobian;
+
+protected:
+  /**
+   * NonlinearSolverSelector configuration data.
+   */
+  AdditionalData additional_data;
+
+private:
+  /**
+   * The MPI communicator to be used by this solver, if any.
+   */
+  MPI_Comm mpi_communicator;
+
+/**
+ * KINSOL configuration data
+ */
+#ifdef DEAL_II_WITH_SUNDIALS
+  typename SUNDIALS::KINSOL<VectorType>::AdditionalData additional_data_kinsol;
+#endif
+
+/**
+ * NOX configuration data
+ */
+#ifdef DEAL_II_WITH_TRILINOS
+  typename TrilinosWrappers::NOXSolver<VectorType>::AdditionalData
+                                       additional_data_nox;
+  Teuchos::RCP<Teuchos::ParameterList> parameters_nox =
+    Teuchos::rcp(new Teuchos::ParameterList);
+#endif
+
+  /**
+   * Data transfer function
+   */
+  void
+  data_transfer(const AdditionalData &additional_data);
+};
+
+template <typename VectorType>
+void
+NonlinearSolverSelector<VectorType>::data_transfer(
+  const AdditionalData &additional_data)
+{
+#ifdef DEAL_II_WITH_SUNDIALS
+  // These if statements pass on the strategy to the other nonlinear solvers
+  if (additional_data.strategy ==
+      NonlinearSolverSelector<VectorType>::AdditionalData::linesearch)
+    additional_data_kinsol.strategy =
+      SUNDIALS::KINSOL<VectorType>::AdditionalData::linesearch;
+  else if (additional_data.strategy ==
+           NonlinearSolverSelector<VectorType>::AdditionalData::newton)
+    additional_data_kinsol.strategy =
+      SUNDIALS::KINSOL<VectorType>::AdditionalData::newton;
+  else if (additional_data.strategy ==
+           NonlinearSolverSelector<VectorType>::AdditionalData::picard)
+    additional_data_kinsol.strategy =
+      SUNDIALS::KINSOL<VectorType>::AdditionalData::picard;
+
+  // Setting data points in the KINSOL class from the NonlinearSolverSelector
+  // class
+  additional_data_kinsol.maximum_non_linear_iterations =
+    additional_data.maximum_non_linear_iterations;
+  additional_data_kinsol.function_tolerance =
+    additional_data.function_tolerance;
+  additional_data_kinsol.step_tolerance = additional_data.step_tolerance;
+  additional_data_kinsol.anderson_subspace_size =
+    additional_data.anderson_subspace_size;
+#endif
+
+// Do the same thing we did above but with NOX
+#ifdef DEAL_II_WITH_TRILINOS
+  // Some default settings for parameters.
+  parameters_nox->set("Nonlinear Solver", "Line Search Based");
+  Teuchos::ParameterList &Line_Search = parameters_nox->sublist("Line Search");
+  Line_Search.set("Method", "Full Step");
+
+  additional_data_nox.max_iter = additional_data.maximum_non_linear_iterations;
+  additional_data_nox.abs_tol  = additional_data.function_tolerance;
+  additional_data_nox.rel_tol  = additional_data.relative_tolerance;
+#endif
+}
+
+template <typename VectorType>
+NonlinearSolverSelector<VectorType>::NonlinearSolverSelector() = default;
+
+template <typename VectorType>
+NonlinearSolverSelector<VectorType>::NonlinearSolverSelector(
+  const AdditionalData &additional_data)
+  : additional_data(additional_data)
+{
+  data_transfer(additional_data);
+}
+
+template <typename VectorType>
+NonlinearSolverSelector<VectorType>::NonlinearSolverSelector(
+  const AdditionalData &additional_data,
+  const MPI_Comm &      mpi_communicator)
+  : additional_data(additional_data)
+  , mpi_communicator(mpi_communicator)
+{
+  data_transfer(additional_data);
+}
+
+
+template <typename VectorType>
+void
+NonlinearSolverSelector<VectorType>::select(
+  const typename AdditionalData::SolverType &type)
+{
+  additional_data.solver_type = type;
+}
+
+template <typename VectorType>
+NonlinearSolverSelector<VectorType>::AdditionalData::AdditionalData(
+  const SolverType &      solver_type,
+  const SolutionStrategy &strategy,
+  const unsigned int      maximum_non_linear_iterations,
+  const double            function_tolerance,
+  const double            relative_tolerance,
+  const double            step_tolerance,
+  const unsigned int      anderson_subspace_size)
+  : solver_type(solver_type)
+  , strategy(strategy)
+  , maximum_non_linear_iterations(maximum_non_linear_iterations)
+  , function_tolerance(function_tolerance)
+  , relative_tolerance(relative_tolerance)
+  , step_tolerance(step_tolerance)
+  , anderson_subspace_size(anderson_subspace_size)
+{}
+
+#ifdef DEAL_II_WITH_TRILINOS
+template <typename VectorType>
+void
+NonlinearSolverSelector<VectorType>::set_data(
+  const typename TrilinosWrappers::NOXSolver<VectorType>::AdditionalData
+    &                                         additional_data,
+  const Teuchos::RCP<Teuchos::ParameterList> &parameters)
+{
+  additional_data_nox = additional_data;
+  parameters_nox      = parameters;
+}
+#endif
+
+#ifdef DEAL_II_WITH_SUNDIALS
+template <typename VectorType>
+void
+NonlinearSolverSelector<VectorType>::set_data(
+  const typename SUNDIALS::KINSOL<VectorType>::AdditionalData &additional_data)
+{
+  additional_data_kinsol = additional_data;
+}
+#endif
+
+template <typename VectorType>
+void
+NonlinearSolverSelector<VectorType>::solve(
+  VectorType &initial_guess_and_solution)
+{
+  // The "auto" solver_type will default to kinsol, however if KINSOL is not
+  // available then we will use NOX.
+  if (additional_data.solver_type == AdditionalData::SolverType::automatic)
+    {
+#ifdef DEAL_II_WITH_TRILINOS
+      additional_data.solver_type = AdditionalData::SolverType::nox;
+#endif
+#ifdef DEAL_II_WITH_SUNDIALS
+      additional_data.solver_type = AdditionalData::SolverType::kinsol;
+#endif
+
+      // If "auto" is still the solver type we cannot solve the problem
+      if (additional_data.solver_type == AdditionalData::SolverType::automatic)
+        AssertThrow(false, ExcMessage("No valid solver type."));
+    }
+
+  if (additional_data.solver_type == AdditionalData::SolverType::kinsol)
+    {
+#ifdef DEAL_II_WITH_SUNDIALS
+      SUNDIALS::KINSOL<VectorType> nonlinear_solver(additional_data_kinsol,
+                                                    mpi_communicator);
+
+      // We set the KINSOL reinit vector equal to the same function
+      // defined for NonlinearSolverSelector.
+      nonlinear_solver.reinit_vector = reinit_vector;
+
+      nonlinear_solver.residual = residual;
+
+      // We cannot simply set these two functions equal to each other
+      // because they have a different number of inputs.
+      nonlinear_solver.setup_jacobian = [&](const VectorType &current_u,
+                                            const VectorType /*&current_f*/) {
+        return NonlinearSolverSelector<VectorType>::setup_jacobian(current_u);
+      };
+
+      nonlinear_solver.solve_with_jacobian = solve_with_jacobian;
+
+      nonlinear_solver.solve(initial_guess_and_solution);
+#else
+      AssertThrow(
+        false, ExcMessage("You do not have SUNDIALS configured with deal.II!"));
+#endif
+    }
+  else if (additional_data.solver_type == AdditionalData::SolverType::nox)
+    {
+#ifdef DEAL_II_WITH_TRILINOS
+      TrilinosWrappers::NOXSolver<VectorType> nonlinear_solver(
+        additional_data_nox, parameters_nox);
+
+      // Do the same thing for NOX that we did with KINSOL.
+      nonlinear_solver.residual = residual;
+
+      // setup_jacobian for NOX has the same number of arguments for the same
+      // function in NonlinearSolverSelector.
+      nonlinear_solver.setup_jacobian = setup_jacobian;
+
+      nonlinear_solver.solve_with_jacobian = solve_with_jacobian;
+
+      nonlinear_solver.solve(initial_guess_and_solution);
+#else
+      AssertThrow(
+        false, ExcMessage("You do not have Trilinos configured with deal.II"));
+#endif
+    }
+  else
+    {
+      std::string solver1;
+      std::string solver2;
+
+#ifdef DEAL_II_WITH_SUNDIALS
+      solver1 = "kinsol \n";
+#endif
+#ifdef DEAL_II_WITH_TRILINOS
+      solver2 = "nox \n";
+#endif
+
+      DeclException2(InvalidNonlinearSolver,
+                     std::string,
+                     std::string,
+                     "Invalid nonlinear solver specified, you may use:\n"
+                       << arg1 << arg2);
+
+      AssertThrow(false, InvalidNonlinearSolver(solver1, solver2));
+    }
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/numerics/nonlinear_solver_selector_01.cc b/tests/numerics/nonlinear_solver_selector_01.cc
new file mode 100644 (file)
index 0000000..266e6fe
--- /dev/null
@@ -0,0 +1,406 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests the class NonlinearSolverSelector using an example based on
+// the step-77 tutorial. The output will vary depending on what
+// packages are configured with deal.II.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/nonlinear.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+#ifdef DEAL_II_WITH_SUNDIALS
+#  include <deal.II/sundials/kinsol.h>
+#endif
+#ifdef DEAL_II_WITH_TRILINOS
+#  include <deal.II/trilinos/nox.h>
+#endif
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+namespace nonlinear_solver_selector_test
+{
+  using namespace dealii;
+
+  using NLSolve = NonlinearSolverSelector<Vector<double>>;
+
+#ifndef SOLVER
+#  define SOLVER NLSolve::AdditionalData::kinsol
+#endif
+
+  template <int dim>
+  class MinimalSurfaceProblem
+  {
+  public:
+    MinimalSurfaceProblem();
+    void
+    run();
+
+  private:
+    void
+    setup_system(const bool initial_step);
+    void
+    solve(const Vector<double> &rhs,
+          Vector<double> &      solution,
+          const double          tolerance);
+    void
+    refine_mesh();
+    void
+    output_results(const unsigned int refinement_cycle);
+    void
+    set_boundary_values();
+    void
+    compute_and_factorize_jacobian(const Vector<double> &evaluation_point);
+    void
+    compute_residual(const Vector<double> &evaluation_point,
+                     Vector<double> &      residual);
+
+    Triangulation<dim> triangulation;
+
+    DoFHandler<dim> dof_handler;
+    FE_Q<dim>       fe;
+
+    AffineConstraints<double> hanging_node_constraints;
+
+    SparsityPattern                      sparsity_pattern;
+    SparseMatrix<double>                 jacobian_matrix;
+    std::unique_ptr<SparseDirectUMFPACK> jacobian_matrix_factorization;
+
+    Vector<double> current_solution;
+  };
+
+
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+  };
+
+
+  template <int dim>
+  double
+  BoundaryValues<dim>::value(const Point<dim> &p,
+                             const unsigned int /*component*/) const
+  {
+    return std::sin(2 * numbers::PI * (p[0] + p[1]));
+  }
+
+
+  template <int dim>
+  MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
+    : dof_handler(triangulation)
+    , fe(1)
+  {}
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
+  {
+    if (initial_step)
+      {
+        dof_handler.distribute_dofs(fe);
+        current_solution.reinit(dof_handler.n_dofs());
+
+        hanging_node_constraints.clear();
+        DoFTools::make_hanging_node_constraints(dof_handler,
+                                                hanging_node_constraints);
+        hanging_node_constraints.close();
+      }
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler, dsp);
+
+    hanging_node_constraints.condense(dsp);
+
+    sparsity_pattern.copy_from(dsp);
+    jacobian_matrix.reinit(sparsity_pattern);
+    jacobian_matrix_factorization.reset();
+  }
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::compute_and_factorize_jacobian(
+    const Vector<double> &evaluation_point)
+  {
+    {
+      deallog << "  Computing Jacobian matrix" << std::endl;
+
+      const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+      jacobian_matrix = 0;
+
+      FEValues<dim> fe_values(fe,
+                              quadrature_formula,
+                              update_gradients | update_quadrature_points |
+                                update_JxW_values);
+
+      const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+      const unsigned int n_q_points    = quadrature_formula.size();
+
+      FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+
+      std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
+
+      std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          cell_matrix = 0;
+
+          fe_values.reinit(cell);
+
+          fe_values.get_function_gradients(evaluation_point,
+                                           evaluation_point_gradients);
+
+          for (unsigned int q = 0; q < n_q_points; ++q)
+            {
+              const double coeff =
+                1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
+                                      evaluation_point_gradients[q]);
+
+              for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                    cell_matrix(i, j) +=
+                      (((fe_values.shape_grad(i, q)    // ((\nabla \phi_i
+                         * coeff                       //   * a_n
+                         * fe_values.shape_grad(j, q)) //   * \nabla \phi_j)
+                        -                              //  -
+                        (fe_values.shape_grad(i, q)    //  (\nabla \phi_i
+                         * coeff * coeff * coeff       //   * a_n^3
+                         *
+                         (fe_values.shape_grad(j, q)       //   * (\nabla \phi_j
+                          * evaluation_point_gradients[q]) //      * \nabla u_n)
+                         * evaluation_point_gradients[q])) //   * \nabla u_n)))
+                       * fe_values.JxW(q));                // * dx
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+          hanging_node_constraints.distribute_local_to_global(cell_matrix,
+                                                              local_dof_indices,
+                                                              jacobian_matrix);
+        }
+
+      std::map<types::global_dof_index, double> boundary_values;
+      VectorTools::interpolate_boundary_values(dof_handler,
+                                               0,
+                                               Functions::ZeroFunction<dim>(),
+                                               boundary_values);
+      Vector<double> dummy_solution(dof_handler.n_dofs());
+      Vector<double> dummy_rhs(dof_handler.n_dofs());
+      MatrixTools::apply_boundary_values(boundary_values,
+                                         jacobian_matrix,
+                                         dummy_solution,
+                                         dummy_rhs);
+    }
+
+    {
+      deallog << "  Factorizing Jacobian matrix" << std::endl;
+
+      jacobian_matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
+      jacobian_matrix_factorization->factorize(jacobian_matrix);
+    }
+  }
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::compute_residual(
+    const Vector<double> &evaluation_point,
+    Vector<double> &      residual)
+  {
+    deallog << "  Computing residual vector..." << std::flush;
+
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+    FEValues<dim>     fe_values(fe,
+                            quadrature_formula,
+                            update_gradients | update_quadrature_points |
+                              update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    Vector<double>              cell_residual(dofs_per_cell);
+    std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        cell_residual = 0;
+        fe_values.reinit(cell);
+
+        fe_values.get_function_gradients(evaluation_point,
+                                         evaluation_point_gradients);
+
+
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            const double coeff =
+              1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
+                                    evaluation_point_gradients[q]);
+
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              cell_residual(i) +=
+                (fe_values.shape_grad(i, q)      // \nabla \phi_i
+                 * coeff                         // * a_n
+                 * evaluation_point_gradients[q] // * \nabla u_n
+                 * fe_values.JxW(q));            // * dx
+          }
+
+        cell->get_dof_indices(local_dof_indices);
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          residual(local_dof_indices[i]) += cell_residual(i);
+      }
+
+    hanging_node_constraints.condense(residual);
+
+    for (const types::global_dof_index i :
+         DoFTools::extract_boundary_dofs(dof_handler))
+      residual(i) = 0;
+
+    for (const types::global_dof_index i :
+         DoFTools::extract_hanging_node_dofs(dof_handler))
+      residual(i) = 0;
+
+    deallog << " norm=" << residual.l2_norm() << std::endl;
+  }
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
+                                    Vector<double> &      solution,
+                                    const double /*tolerance*/)
+  {
+    deallog << "  Solving linear system" << std::endl;
+
+    jacobian_matrix_factorization->vmult(solution, rhs);
+
+    hanging_node_constraints.distribute(solution);
+  }
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::set_boundary_values()
+  {
+    std::map<types::global_dof_index, double> boundary_values;
+    VectorTools::interpolate_boundary_values(dof_handler,
+                                             0,
+                                             BoundaryValues<dim>(),
+                                             boundary_values);
+    for (const auto &boundary_value : boundary_values)
+      current_solution(boundary_value.first) = boundary_value.second;
+
+    hanging_node_constraints.distribute(current_solution);
+  }
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::run()
+  {
+    GridGenerator::hyper_ball(triangulation);
+    triangulation.refine_global(2);
+
+    setup_system(/*initial_step=*/true);
+    set_boundary_values();
+
+    {
+      typename NLSolve::AdditionalData additional_data;
+      additional_data.solver_type = SOLVER;
+
+      NLSolve nonlinear_solver(additional_data);
+
+      nonlinear_solver.reinit_vector = [&](Vector<double> &x) {
+        x.reinit(dof_handler.n_dofs());
+      };
+
+      nonlinear_solver.residual = [&](const Vector<double> &evaluation_point,
+                                      Vector<double> &      residual) {
+        compute_residual(evaluation_point, residual);
+
+        return 0;
+      };
+
+      nonlinear_solver.setup_jacobian = [&](const Vector<double> &current_u) {
+        compute_and_factorize_jacobian(current_u);
+
+        return 0;
+      };
+
+      nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
+                                                 Vector<double> &      dst,
+                                                 const double tolerance) {
+        this->solve(rhs, dst, tolerance);
+
+        return 0;
+      };
+
+      nonlinear_solver.solve(current_solution);
+    }
+  }
+} // namespace nonlinear_solver_selector_test
+
+
+int
+main()
+{
+  initlog();
+
+  using namespace nonlinear_solver_selector_test;
+
+  MinimalSurfaceProblem<2> laplace_problem_2d;
+  laplace_problem_2d.run();
+
+  return 0;
+}
diff --git a/tests/numerics/nonlinear_solver_selector_01.with_sundials=on.output b/tests/numerics/nonlinear_solver_selector_01.with_sundials=on.output
new file mode 100644 (file)
index 0000000..c090d90
--- /dev/null
@@ -0,0 +1,71 @@
+
+DEAL::  Computing residual vector... norm=0.867975
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.867975
+DEAL::  Computing residual vector... norm=0.212073
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.212073
+DEAL::  Computing residual vector... norm=0.202631
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.202631
+DEAL::  Computing residual vector... norm=0.165773
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.165774
+DEAL::  Computing residual vector... norm=0.162594
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.162594
+DEAL::  Computing residual vector... norm=0.148175
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.148175
+DEAL::  Computing residual vector... norm=0.145391
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.145391
+DEAL::  Computing residual vector... norm=0.137551
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.137551
+DEAL::  Computing residual vector... norm=0.135366
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.135365
+DEAL::  Computing residual vector... norm=0.130367
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.130367
+DEAL::  Computing residual vector... norm=0.128704
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.128704
+DEAL::  Computing residual vector... norm=0.0302623
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.0302624
+DEAL::  Computing residual vector... norm=0.0126764
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.0126763
+DEAL::  Computing residual vector... norm=0.00488315
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.00488322
+DEAL::  Computing residual vector... norm=0.00195788
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.00195781
+DEAL::  Computing residual vector... norm=0.000773169
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.000773247
+DEAL::  Computing residual vector... norm=0.000307242
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.000307164
+DEAL::  Computing residual vector... norm=0.000121790
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.000121868
+DEAL::  Computing residual vector... norm=4.83248e-05
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=4.82467e-05
+DEAL::  Computing residual vector... norm=1.91672e-05
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=1.92453e-05
+DEAL::  Computing residual vector... norm=7.60355e-06
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=7.52545e-06
+DEAL::  Computing residual vector... norm=6.25728e-11
diff --git a/tests/numerics/nonlinear_solver_selector_02.cc b/tests/numerics/nonlinear_solver_selector_02.cc
new file mode 100644 (file)
index 0000000..78d0d9a
--- /dev/null
@@ -0,0 +1,24 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests the class NonlinearSolverSelector using an example based on
+// the test nonlinear_solver_selector_01. Here we use the nonlinear
+// solver NOX instead of KINSOL.
+
+#define SOLVER NonlinearSolverSelector<Vector<double>>::AdditionalData::nox
+
+#include "nonlinear_solver_selector_01.cc"
diff --git a/tests/numerics/nonlinear_solver_selector_02.with_trilinos=on.output b/tests/numerics/nonlinear_solver_selector_02.with_trilinos=on.output
new file mode 100644 (file)
index 0000000..64c2aa3
--- /dev/null
@@ -0,0 +1,18 @@
+
+DEAL::  Computing residual vector... norm=0.867975
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.212073
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.0189603
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=0.000314854
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::  Computing residual vector... norm=1.14048e-07
diff --git a/tests/numerics/nonlinear_solver_selector_03.cc b/tests/numerics/nonlinear_solver_selector_03.cc
new file mode 100644 (file)
index 0000000..036aa75
--- /dev/null
@@ -0,0 +1,491 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests the NonlinearSolverSelector class using an examplebased on the
+// step-77 tutorial program. This test checks the compatability of the
+// class with MPI.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+// Included from step-40
+#include <deal.II/lac/generic_linear_algebra.h>
+
+#define FORCE_USE_OF_TRILINOS
+
+namespace LA
+{
+#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
+  !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
+  using namespace dealii::LinearAlgebraPETSc;
+#  define USE_PETSC_LA
+#elif defined(DEAL_II_WITH_TRILINOS)
+  using namespace dealii::LinearAlgebraTrilinos;
+#else
+#  error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
+#endif
+} // namespace LA
+
+#include <deal.II/base/conditional_ostream.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/nonlinear.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+namespace MPI_nonlinear_solver_selector_test
+{
+  using namespace dealii;
+
+  using NLSolve = NonlinearSolverSelector<LA::MPI::Vector>;
+
+#ifndef SOLVER
+#  define SOLVER NLSolve::AdditionalData::kinsol
+#endif
+
+  template <int dim>
+  class MinimalSurfaceProblem
+  {
+  public:
+    MinimalSurfaceProblem();
+    void
+    run();
+
+  private:
+    void
+    setup_system(const bool initial_step);
+    void
+    solve(const LA::MPI::Vector &rhs,
+          LA::MPI::Vector &      solution,
+          const double           tolerance);
+    void
+    compute_and_factorize_jacobian(const LA::MPI::Vector &evaluation_point);
+    void
+    compute_residual(const LA::MPI::Vector &evaluation_point,
+                     LA::MPI::Vector &      residual);
+
+    MPI_Comm mpi_communicator;
+
+    parallel::distributed::Triangulation<dim> triangulation;
+
+    DoFHandler<dim> dof_handler;
+    FE_Q<dim>       fe;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    AffineConstraints<double> nonzero_constraints;
+    AffineConstraints<double> zero_constraints;
+
+    LA::MPI::SparseMatrix jacobian_matrix;
+
+    LA::MPI::Vector current_solution;
+  };
+
+
+  template <int dim>
+  MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
+    : mpi_communicator(MPI_COMM_WORLD)
+    , triangulation(mpi_communicator,
+                    typename Triangulation<dim>::MeshSmoothing(
+                      Triangulation<dim>::smoothing_on_refinement |
+                      Triangulation<dim>::smoothing_on_coarsening))
+    , dof_handler(triangulation)
+    , fe(1)
+  {}
+
+
+
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const override;
+  };
+
+
+  template <int dim>
+  double
+  BoundaryValues<dim>::value(const Point<dim> &p,
+                             const unsigned int /*component*/) const
+  {
+    return std::sin(2 * numbers::PI * (p[0] + p[1]));
+  };
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
+  {
+    if (initial_step)
+      {
+        dof_handler.distribute_dofs(fe);
+
+        locally_owned_dofs = dof_handler.locally_owned_dofs();
+        locally_relevant_dofs =
+          DoFTools::extract_locally_relevant_dofs(dof_handler);
+
+        current_solution.reinit(locally_owned_dofs, mpi_communicator);
+
+        {
+          nonzero_constraints.clear();
+          nonzero_constraints.reinit(locally_relevant_dofs);
+          DoFTools::make_hanging_node_constraints(dof_handler,
+                                                  nonzero_constraints);
+
+          nonzero_constraints.close();
+
+          nonzero_constraints.distribute(current_solution);
+
+          std::map<types::global_dof_index, double> boundary_values;
+          VectorTools::interpolate_boundary_values(dof_handler,
+                                                   0,
+                                                   BoundaryValues<dim>(),
+                                                   boundary_values);
+
+          for (const auto &boundary_value : boundary_values)
+            current_solution(boundary_value.first) = boundary_value.second;
+        }
+
+        {
+          zero_constraints.clear();
+          zero_constraints.reinit(locally_relevant_dofs);
+          DoFTools::make_hanging_node_constraints(dof_handler,
+                                                  zero_constraints);
+          VectorTools::interpolate_boundary_values(
+            dof_handler, 0, Functions::ZeroFunction<dim>(), zero_constraints);
+        }
+        zero_constraints.close();
+      }
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, zero_constraints, false);
+
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               dof_handler.locally_owned_dofs(),
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    jacobian_matrix.reinit(locally_owned_dofs,
+                           locally_owned_dofs,
+                           dsp,
+                           mpi_communicator);
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::compute_and_factorize_jacobian(
+    const LA::MPI::Vector &evaluation_point)
+  {
+    LA::MPI::Vector evaluation_point_1;
+    evaluation_point_1.reinit(locally_owned_dofs,
+                              locally_relevant_dofs,
+                              mpi_communicator);
+    evaluation_point_1 = evaluation_point;
+
+    {
+      deallog << "  Computing Jacobian matrix" << std::endl;
+
+      const QGauss<dim> quadrature_formula(fe.degree + 1);
+
+      jacobian_matrix = 0;
+
+      FEValues<dim> fe_values(fe,
+                              quadrature_formula,
+                              update_gradients | update_quadrature_points |
+                                update_JxW_values);
+
+      const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+      const unsigned int n_q_points    = quadrature_formula.size();
+
+      FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+
+      std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
+
+      std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+      for (const auto &cell : dof_handler.active_cell_iterators())
+        {
+          if (cell->is_locally_owned())
+            {
+              cell_matrix = 0.;
+
+              fe_values.reinit(cell);
+
+              fe_values.get_function_gradients(evaluation_point_1,
+                                               evaluation_point_gradients);
+
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                {
+                  const double coeff =
+                    1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
+                                          evaluation_point_gradients[q]);
+
+                  for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                    {
+                      for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                        {
+                          cell_matrix(i, j) +=
+                            (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
+                               * coeff                    //   * a_n
+                               *
+                               fe_values.shape_grad(j, q)) //   * \nabla \phi_j)
+                              -                            //  -
+                              (fe_values.shape_grad(i, q)  //  (\nabla \phi_i
+                               * coeff * coeff * coeff     //   * a_n^3
+                               *
+                               (fe_values.shape_grad(j, q) //   * (\nabla \phi_j
+                                *
+                                evaluation_point_gradients[q]) //      * \nabla
+                                                               //      u_n)
+                               * evaluation_point_gradients[q])) //   * \nabla
+                                                                 //   u_n)))
+                             * fe_values.JxW(q));                // * dx
+                        }
+                    }
+                }
+
+              cell->get_dof_indices(local_dof_indices);
+
+              zero_constraints.distribute_local_to_global(cell_matrix,
+                                                          local_dof_indices,
+                                                          jacobian_matrix);
+            }
+        }
+    }
+
+    jacobian_matrix.compress(VectorOperation::add);
+
+    deallog << "  Factorizing Jacobian matrix" << std::endl;
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::compute_residual(
+    const LA::MPI::Vector &evaluation_point,
+    LA::MPI::Vector &      residual)
+  {
+    deallog << "  Computing residual vector..." << std::flush;
+
+    LA::MPI::Vector evaluation_point_1;
+    evaluation_point_1.reinit(locally_owned_dofs,
+                              locally_relevant_dofs,
+                              mpi_communicator);
+    evaluation_point_1 = evaluation_point;
+
+    const QGauss<dim> quadrature_formula(fe.degree + 1);
+    FEValues<dim>     fe_values(fe,
+                            quadrature_formula,
+                            update_gradients | update_quadrature_points |
+                              update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    Vector<double>              cell_residual(dofs_per_cell);
+    std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        if (cell->is_locally_owned())
+          {
+            cell_residual = 0.;
+            fe_values.reinit(cell);
+
+            fe_values.get_function_gradients(evaluation_point_1,
+                                             evaluation_point_gradients);
+
+
+            for (unsigned int q = 0; q < n_q_points; ++q)
+              {
+                const double coeff =
+                  1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
+                                        evaluation_point_gradients[q]);
+
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  cell_residual(i) =
+                    (fe_values.shape_grad(i, q)      // \nabla \phi_i
+                     * coeff                         // * a_n
+                     * evaluation_point_gradients[q] // * \nabla u_n
+                     * fe_values.JxW(q));            // * dx
+              }
+
+            cell->get_dof_indices(local_dof_indices);
+
+            zero_constraints.distribute_local_to_global(cell_residual,
+                                                        local_dof_indices,
+                                                        residual);
+          }
+      }
+
+    zero_constraints.set_zero(residual);
+    residual.compress(VectorOperation::add);
+
+    deallog << " norm=" << residual.l2_norm() << std::endl;
+  }
+
+
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::solve(const LA::MPI::Vector &rhs,
+                                    LA::MPI::Vector &      solution,
+                                    const double /*tolerance*/)
+  {
+    deallog << "  Solving linear system" << std::endl;
+
+    SolverControl solver_control(dof_handler.n_dofs(), 1e-12);
+
+#ifdef USE_PETSC_LA
+    LA::SolverCG solver(solver_control, mpi_communicator);
+#else
+    LA::SolverCG solver(solver_control);
+#endif
+
+    LA::MPI::PreconditionAMG preconditioner;
+
+    LA::MPI::PreconditionAMG::AdditionalData data;
+
+#ifdef USE_PETSC_LA
+    data.symmetric_operator = true;
+#else
+    /* Trilinos defaults are good */
+#endif
+    preconditioner.initialize(jacobian_matrix, data);
+
+    solver.solve(jacobian_matrix, solution, rhs, preconditioner);
+
+    deallog << "   Solved in " << solver_control.last_step() << " iterations."
+            << std::endl;
+
+    zero_constraints.distribute(solution);
+  }
+
+  template <int dim>
+  void
+  MinimalSurfaceProblem<dim>::run()
+  {
+    deallog << "Running with "
+#ifdef USE_PETSC_LA
+            << "PETSc"
+#else
+            << "Trilinos"
+#endif
+            << " on " << Utilities::MPI::n_mpi_processes(mpi_communicator)
+            << " MPI rank(s)..." << std::endl;
+
+    GridGenerator::hyper_ball(triangulation);
+    triangulation.refine_global(4);
+
+    const bool initial_step = true;
+
+    setup_system(initial_step);
+
+    const double target_tolerance = 1e-3;
+    deallog << "  Target_tolerance: " << target_tolerance << std::endl
+            << std::endl;
+
+    typename NLSolve::AdditionalData additional_data;
+    additional_data.function_tolerance = target_tolerance;
+    additional_data.solver_type        = SOLVER;
+
+    NLSolve nonlinear_solver(additional_data, mpi_communicator);
+
+    nonlinear_solver.reinit_vector = [&](LA::MPI::Vector &x) {
+      x.reinit(locally_owned_dofs, mpi_communicator);
+    };
+
+    nonlinear_solver.residual = [&](const LA::MPI::Vector &evaluation_point,
+                                    LA::MPI::Vector &      residual) {
+      compute_residual(evaluation_point, residual);
+
+      return 0;
+    };
+
+    nonlinear_solver.setup_jacobian = [&](const LA::MPI::Vector &current_u) {
+      compute_and_factorize_jacobian(current_u);
+
+      return 0;
+    };
+
+    nonlinear_solver.solve_with_jacobian = [&](const LA::MPI::Vector &rhs,
+                                               LA::MPI::Vector &      dst,
+                                               const double tolerance) {
+      this->solve(rhs, dst, tolerance);
+
+      return 0;
+    };
+
+    nonlinear_solver.solve(current_solution);
+
+    deallog << std::endl;
+  }
+} // namespace MPI_nonlinear_solver_selector_test
+
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  using namespace MPI_nonlinear_solver_selector_test;
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  MinimalSurfaceProblem<2> laplace_problem_2d;
+  laplace_problem_2d.run();
+
+  return 0;
+}
diff --git a/tests/numerics/nonlinear_solver_selector_03.with_sundials=on.mpirun=4.output b/tests/numerics/nonlinear_solver_selector_03.with_sundials=on.mpirun=4.output
new file mode 100644 (file)
index 0000000..94ceaeb
--- /dev/null
@@ -0,0 +1,131 @@
+
+DEAL::Running with Trilinos on 4 MPI rank(s)...
+DEAL::  Target_tolerance: 0.00100000
+DEAL::
+DEAL::  Computing residual vector... norm=0.170956
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 9.85645e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.170956
+DEAL::  Computing residual vector... norm=0.129347
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 7.11784e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.129347
+DEAL::  Computing residual vector... norm=0.103614
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 5.28808e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.103614
+DEAL::  Computing residual vector... norm=0.0861914
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 3.98596e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0861914
+DEAL::  Computing residual vector... norm=0.0865470
+DEAL::  Computing residual vector... norm=0.0795265
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 3.14097e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0795265
+DEAL::  Computing residual vector... norm=0.0814816
+DEAL::  Computing residual vector... norm=0.0774703
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 2.72056e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0774703
+DEAL::  Computing residual vector... norm=0.0715560
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 2.76383e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0715560
+DEAL::  Computing residual vector... norm=0.0635000
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 4.09086e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0635000
+DEAL::  Computing residual vector... norm=0.0647100
+DEAL::  Computing residual vector... norm=0.0549834
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 3.01833e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0549834
+DEAL::  Computing residual vector... norm=0.0638162
+DEAL::  Computing residual vector... norm=0.0525536
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 3.51125e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0525536
+DEAL::  Computing residual vector... norm=0.0468639
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 3.11947e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0468639
+DEAL::  Computing residual vector... norm=0.0361431
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 2.33506e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0361430
+DEAL::  Computing residual vector... norm=0.0292451
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 1.77355e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0292451
+DEAL::  Computing residual vector... norm=0.0244098
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 1.36440e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0244098
+DEAL::  Computing residual vector... norm=0.0208326
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 1.06488e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0208326
+DEAL::  Computing residual vector... norm=0.0180857
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 8.45632e-14
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0180857
+DEAL::  Computing residual vector... norm=0.0159160
+DEAL::  Computing residual vector... norm=0.0139269
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 5.34166e-14
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0139268
+DEAL::  Computing residual vector... norm=0.0125307
+DEAL::  Computing residual vector... norm=0.0112272
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 3.69276e-14
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0112272
+DEAL::  Computing residual vector... norm=0.0102440
+DEAL::  Computing residual vector... norm=0.00931379
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 2.64848e-14
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.00931378
+DEAL::  Computing residual vector... norm=0.00857916
+DEAL::  Computing residual vector... norm=0.00787754
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 1.93828e-14
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.00787753
+DEAL::  Computing residual vector... norm=0.00730638
+DEAL::  Computing residual vector... norm=0.00675703
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 3.03169e-14
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.00675702
+DEAL::  Computing residual vector... norm=0.00509054
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 2.28021e-14
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.00509053
+DEAL::  Computing residual vector... norm=0.00400011
+DEAL::
diff --git a/tests/numerics/nonlinear_solver_selector_04.cc b/tests/numerics/nonlinear_solver_selector_04.cc
new file mode 100644 (file)
index 0000000..7173a5f
--- /dev/null
@@ -0,0 +1,24 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2007 - 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests the class NonlinearSolverSelector using an example based on
+// the test nonlinear_solver_selector_01. Here we use the nonlinear
+// solver NOX instead of KINSOL with MPI.
+
+#define SOLVER NonlinearSolverSelector<LA::MPI::Vector>::AdditionalData::nox
+
+#include "nonlinear_solver_selector_03.cc"
diff --git a/tests/numerics/nonlinear_solver_selector_04.with_trilinos=on.mpirun=4.output b/tests/numerics/nonlinear_solver_selector_04.with_trilinos=on.mpirun=4.output
new file mode 100644 (file)
index 0000000..e37a91e
--- /dev/null
@@ -0,0 +1,42 @@
+
+DEAL::Running with Trilinos on 4 MPI rank(s)...
+DEAL::  Target_tolerance: 0.00100000
+DEAL::
+DEAL::  Computing residual vector... norm=0.170956
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 9.85645e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.129347
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 6.57959e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.100710
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 4.51371e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0776513
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 2 value 2.44269e-16
+DEAL::   Solved in 2 iterations.
+DEAL::  Computing residual vector... norm=0.0593571
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 8.24123e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0459878
+DEAL::  Computing Jacobian matrix
+DEAL::  Factorizing Jacobian matrix
+DEAL::  Solving linear system
+DEAL::Convergence step 1 value 4.92582e-13
+DEAL::   Solved in 1 iterations.
+DEAL::  Computing residual vector... norm=0.0355434
+DEAL::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.