--- /dev/null
+New: Added new nonlinear solver class NonlinearSolverSelector.
+<br>
+(Sean Ingimarson, 2023/05/01)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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> ¤t_u,
+ * const Vector<double> ¤t_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> ¶meters =
+ 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 ¤t_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> ¶meters)
+{
+ 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 ¤t_u,
+ const VectorType /*¤t_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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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> ¤t_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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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"
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ¤t_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;
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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"
--- /dev/null
+
+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::
projects / dealii.git / commitdiff