From 435cae94b4ec8b9630a2264a9340f6f7d6550288 Mon Sep 17 00:00:00 2001 From: Sean Ingimarson Date: Wed, 5 Apr 2023 08:03:36 -0400 Subject: [PATCH] NonlinearSolverSelector class implementation --- doc/news/changes/minor/20230501Ingimarson | 3 + include/deal.II/numerics/nonlinear.h | 558 ++++++++++++++++++ .../numerics/nonlinear_solver_selector_01.cc | 406 +++++++++++++ ...solver_selector_01.with_sundials=on.output | 71 +++ .../numerics/nonlinear_solver_selector_02.cc | 24 + ...solver_selector_02.with_trilinos=on.output | 18 + .../numerics/nonlinear_solver_selector_03.cc | 491 +++++++++++++++ ...lector_03.with_sundials=on.mpirun=4.output | 131 ++++ .../numerics/nonlinear_solver_selector_04.cc | 24 + ...lector_04.with_trilinos=on.mpirun=4.output | 42 ++ 10 files changed, 1768 insertions(+) create mode 100644 doc/news/changes/minor/20230501Ingimarson create mode 100644 include/deal.II/numerics/nonlinear.h create mode 100644 tests/numerics/nonlinear_solver_selector_01.cc create mode 100644 tests/numerics/nonlinear_solver_selector_01.with_sundials=on.output create mode 100644 tests/numerics/nonlinear_solver_selector_02.cc create mode 100644 tests/numerics/nonlinear_solver_selector_02.with_trilinos=on.output create mode 100644 tests/numerics/nonlinear_solver_selector_03.cc create mode 100644 tests/numerics/nonlinear_solver_selector_03.with_sundials=on.mpirun=4.output create mode 100644 tests/numerics/nonlinear_solver_selector_04.cc create mode 100644 tests/numerics/nonlinear_solver_selector_04.with_trilinos=on.mpirun=4.output diff --git a/doc/news/changes/minor/20230501Ingimarson b/doc/news/changes/minor/20230501Ingimarson new file mode 100644 index 0000000000..0da568961e --- /dev/null +++ b/doc/news/changes/minor/20230501Ingimarson @@ -0,0 +1,3 @@ +New: Added new nonlinear solver class NonlinearSolverSelector. +
+(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 index 0000000000..11bb956a8c --- /dev/null +++ b/include/deal.II/numerics/nonlinear.h @@ -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 + +#include + +#include + +#include + + +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. + * + *

Usage

+ * An example of code one would run with this class is the + * following: + * @code + * // Generate a @p NonlinearSolverSelector that uses @p KINSOL + * NonlinearSolverSelector>::AdditionalData additional_data; + * additional_data.solver_type = + * NonlinearSolverSelector>::AdditionalData::SolverType::kinsol; + * + * NonlinearSolverSelector> 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 &x) {...} + * + * nonlinear_solver.residual = + * [&](const Vector &evaluation_point, + * Vector & residual) {...} + * + * nonlinear_solver.setup_jacobian = + * [&](const Vector ¤t_u, + * const Vector ¤t_f) {...} + * + * nonlinear_solver.solve_with_jacobian = + * [&](const Vector &rhs, + * Vector & dst, + * const double tolerance) {...} + * + * // Calling the @p solve function with an initial guess. + * nonlinear_solver.solve(current_solution); + * @endcode + */ +template > +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::AdditionalData + & additional_data, + const Teuchos::RCP ¶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::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 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 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 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::AdditionalData additional_data_kinsol; +#endif + +/** + * NOX configuration data + */ +#ifdef DEAL_II_WITH_TRILINOS + typename TrilinosWrappers::NOXSolver::AdditionalData + additional_data_nox; + Teuchos::RCP parameters_nox = + Teuchos::rcp(new Teuchos::ParameterList); +#endif + + /** + * Data transfer function + */ + void + data_transfer(const AdditionalData &additional_data); +}; + +template +void +NonlinearSolverSelector::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::AdditionalData::linesearch) + additional_data_kinsol.strategy = + SUNDIALS::KINSOL::AdditionalData::linesearch; + else if (additional_data.strategy == + NonlinearSolverSelector::AdditionalData::newton) + additional_data_kinsol.strategy = + SUNDIALS::KINSOL::AdditionalData::newton; + else if (additional_data.strategy == + NonlinearSolverSelector::AdditionalData::picard) + additional_data_kinsol.strategy = + SUNDIALS::KINSOL::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 +NonlinearSolverSelector::NonlinearSolverSelector() = default; + +template +NonlinearSolverSelector::NonlinearSolverSelector( + const AdditionalData &additional_data) + : additional_data(additional_data) +{ + data_transfer(additional_data); +} + +template +NonlinearSolverSelector::NonlinearSolverSelector( + const AdditionalData &additional_data, + const MPI_Comm & mpi_communicator) + : additional_data(additional_data) + , mpi_communicator(mpi_communicator) +{ + data_transfer(additional_data); +} + + +template +void +NonlinearSolverSelector::select( + const typename AdditionalData::SolverType &type) +{ + additional_data.solver_type = type; +} + +template +NonlinearSolverSelector::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 +void +NonlinearSolverSelector::set_data( + const typename TrilinosWrappers::NOXSolver::AdditionalData + & additional_data, + const Teuchos::RCP ¶meters) +{ + additional_data_nox = additional_data; + parameters_nox = parameters; +} +#endif + +#ifdef DEAL_II_WITH_SUNDIALS +template +void +NonlinearSolverSelector::set_data( + const typename SUNDIALS::KINSOL::AdditionalData &additional_data) +{ + additional_data_kinsol = additional_data; +} +#endif + +template +void +NonlinearSolverSelector::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 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::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 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 index 0000000000..266e6fe7fb --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_01.cc @@ -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 +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#ifdef DEAL_II_WITH_SUNDIALS +# include +#endif +#ifdef DEAL_II_WITH_TRILINOS +# include +#endif + +#include +#include + +#include "../tests.h" + + +namespace nonlinear_solver_selector_test +{ + using namespace dealii; + + using NLSolve = NonlinearSolverSelector>; + +#ifndef SOLVER +# define SOLVER NLSolve::AdditionalData::kinsol +#endif + + template + class MinimalSurfaceProblem + { + public: + MinimalSurfaceProblem(); + void + run(); + + private: + void + setup_system(const bool initial_step); + void + solve(const Vector &rhs, + Vector & 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 &evaluation_point); + void + compute_residual(const Vector &evaluation_point, + Vector & residual); + + Triangulation triangulation; + + DoFHandler dof_handler; + FE_Q fe; + + AffineConstraints hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix jacobian_matrix; + std::unique_ptr jacobian_matrix_factorization; + + Vector current_solution; + }; + + + template + class BoundaryValues : public Function + { + public: + virtual double + value(const Point &p, const unsigned int component = 0) const override; + }; + + + template + double + BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const + { + return std::sin(2 * numbers::PI * (p[0] + p[1])); + } + + + template + MinimalSurfaceProblem::MinimalSurfaceProblem() + : dof_handler(triangulation) + , fe(1) + {} + + + template + void + MinimalSurfaceProblem::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 + void + MinimalSurfaceProblem::compute_and_factorize_jacobian( + const Vector &evaluation_point) + { + { + deallog << " Computing Jacobian matrix" << std::endl; + + const QGauss quadrature_formula(fe.degree + 1); + + jacobian_matrix = 0; + + FEValues 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 cell_matrix(dofs_per_cell, dofs_per_cell); + + std::vector> evaluation_point_gradients(n_q_points); + + std::vector 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 boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + boundary_values); + Vector dummy_solution(dof_handler.n_dofs()); + Vector 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(); + jacobian_matrix_factorization->factorize(jacobian_matrix); + } + } + + + template + void + MinimalSurfaceProblem::compute_residual( + const Vector &evaluation_point, + Vector & residual) + { + deallog << " Computing residual vector..." << std::flush; + + const QGauss quadrature_formula(fe.degree + 1); + FEValues 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 cell_residual(dofs_per_cell); + std::vector> evaluation_point_gradients(n_q_points); + + std::vector 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 + void + MinimalSurfaceProblem::solve(const Vector &rhs, + Vector & solution, + const double /*tolerance*/) + { + deallog << " Solving linear system" << std::endl; + + jacobian_matrix_factorization->vmult(solution, rhs); + + hanging_node_constraints.distribute(solution); + } + + + template + void + MinimalSurfaceProblem::set_boundary_values() + { + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + BoundaryValues(), + boundary_values); + for (const auto &boundary_value : boundary_values) + current_solution(boundary_value.first) = boundary_value.second; + + hanging_node_constraints.distribute(current_solution); + } + + + template + void + MinimalSurfaceProblem::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 &x) { + x.reinit(dof_handler.n_dofs()); + }; + + nonlinear_solver.residual = [&](const Vector &evaluation_point, + Vector & residual) { + compute_residual(evaluation_point, residual); + + return 0; + }; + + nonlinear_solver.setup_jacobian = [&](const Vector ¤t_u) { + compute_and_factorize_jacobian(current_u); + + return 0; + }; + + nonlinear_solver.solve_with_jacobian = [&](const Vector &rhs, + Vector & 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 index 0000000000..c090d903f8 --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_01.with_sundials=on.output @@ -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 index 0000000000..78d0d9a401 --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_02.cc @@ -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>::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 index 0000000000..64c2aa3c0c --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_02.with_trilinos=on.output @@ -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 index 0000000000..036aa753ab --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_03.cc @@ -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 +#include +#include + +// Included from step-40 +#include + +#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 + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + +namespace MPI_nonlinear_solver_selector_test +{ + using namespace dealii; + + using NLSolve = NonlinearSolverSelector; + +#ifndef SOLVER +# define SOLVER NLSolve::AdditionalData::kinsol +#endif + + template + 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 triangulation; + + DoFHandler dof_handler; + FE_Q fe; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints nonzero_constraints; + AffineConstraints zero_constraints; + + LA::MPI::SparseMatrix jacobian_matrix; + + LA::MPI::Vector current_solution; + }; + + + template + MinimalSurfaceProblem::MinimalSurfaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , dof_handler(triangulation) + , fe(1) + {} + + + + template + class BoundaryValues : public Function + { + public: + virtual double + value(const Point &p, const unsigned int component = 0) const override; + }; + + + template + double + BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const + { + return std::sin(2 * numbers::PI * (p[0] + p[1])); + }; + + + template + void + MinimalSurfaceProblem::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 boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + BoundaryValues(), + 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(), 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 + void + MinimalSurfaceProblem::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 quadrature_formula(fe.degree + 1); + + jacobian_matrix = 0; + + FEValues 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 cell_matrix(dofs_per_cell, dofs_per_cell); + + std::vector> evaluation_point_gradients(n_q_points); + + std::vector 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 + void + MinimalSurfaceProblem::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 quadrature_formula(fe.degree + 1); + FEValues 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 cell_residual(dofs_per_cell); + std::vector> evaluation_point_gradients(n_q_points); + + std::vector 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 + void + MinimalSurfaceProblem::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 + void + MinimalSurfaceProblem::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; +} 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 index 0000000000..94ceaebddf --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_03.with_sundials=on.mpirun=4.output @@ -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 index 0000000000..7173a5f45e --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_04.cc @@ -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::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 index 0000000000..e37a91eac9 --- /dev/null +++ b/tests/numerics/nonlinear_solver_selector_04.with_trilinos=on.mpirun=4.output @@ -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:: -- 2.39.5