From 741668572d067c7c9a1dbb5de224e64b776e820c Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Fri, 25 Nov 2022 19:45:47 +0300 Subject: [PATCH] PETScWrappers: Add support for nonlinear solver --- include/deal.II/lac/petsc_snes.h | 413 ++++++++++++++++++++ include/deal.II/lac/petsc_snes.templates.h | 425 +++++++++++++++++++++ source/lac/CMakeLists.txt | 1 + source/lac/petsc_snes.cc | 74 ++++ tests/petsc/petsc_snes_00.cc | 206 ++++++++++ tests/petsc/petsc_snes_00.output | 40 ++ tests/petsc/petsc_snes_00_in.prm | 8 + tests/petsc/petsc_snes_01.cc | 182 +++++++++ tests/petsc/petsc_snes_01.output | 9 + tests/petsc/petsc_snes_02.cc | 85 +++++ tests/petsc/petsc_snes_02.output | 11 + tests/petsc/petsc_snes_03.cc | 70 ++++ tests/petsc/petsc_snes_03.mpirun=3.output | 3 + tests/petsc/petsc_snes_03.output | 3 + tests/petsc/petsc_snes_04.cc | 65 ++++ tests/petsc/petsc_snes_04.mpirun=3.output | 12 + tests/petsc/petsc_snes_04.output | 12 + 17 files changed, 1619 insertions(+) create mode 100644 include/deal.II/lac/petsc_snes.h create mode 100644 include/deal.II/lac/petsc_snes.templates.h create mode 100644 source/lac/petsc_snes.cc create mode 100644 tests/petsc/petsc_snes_00.cc create mode 100644 tests/petsc/petsc_snes_00.output create mode 100644 tests/petsc/petsc_snes_00_in.prm create mode 100644 tests/petsc/petsc_snes_01.cc create mode 100644 tests/petsc/petsc_snes_01.output create mode 100644 tests/petsc/petsc_snes_02.cc create mode 100644 tests/petsc/petsc_snes_02.output create mode 100644 tests/petsc/petsc_snes_03.cc create mode 100644 tests/petsc/petsc_snes_03.mpirun=3.output create mode 100644 tests/petsc/petsc_snes_03.output create mode 100644 tests/petsc/petsc_snes_04.cc create mode 100644 tests/petsc/petsc_snes_04.mpirun=3.output create mode 100644 tests/petsc/petsc_snes_04.output diff --git a/include/deal.II/lac/petsc_snes.h b/include/deal.II/lac/petsc_snes.h new file mode 100644 index 0000000000..ebae3fa4da --- /dev/null +++ b/include/deal.II/lac/petsc_snes.h @@ -0,0 +1,413 @@ +//----------------------------------------------------------- +// +// 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. +// +//--------------------------------------------------------------- + +#ifndef dealii_petsc_snes_h +#define dealii_petsc_snes_h + +#include + +#ifdef DEAL_II_WITH_PETSC +# include +# include +# include + +# include +# include +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace PETScWrappers +{ + /** + * Additional parameters that can be passed to the NonlinearSolver class. + */ + class NonlinearSolverData + { + public: + /** + * Type that holds real-valued numbers. + * + * Used to represent norms. + */ + using real_type = PetscReal; + + /** + * Initialization parameters for NonlinearSolverData. + * + * Running parameters: + * + * @param options_prefix The string indicating the options prefix for command line customization. + * @param snes_type The string indicating the PETSc SNES solver type. + * @param snes_linesearch_type The string indicating the PETSc linesearch type. + * @param absolute_tolerance Absolute error tolerance. + * @param relative_tolerance Relative error tolerance. + * @param step_tolerance Step tolerance. + * @param maximum_non_linear_iterations Maximum number of iterations allowed. + * @param max_n_function_evaluations Maximum number of function evaluations allowed. + * + * @note All parameters values specified here can be overriden by + * command line choices. + * + * @ingroup PETScWrappers + */ + NonlinearSolverData( + // Running parameters + const std::string &options_prefix = "", + const std::string &snes_type = "", + const std::string &snes_linesearch_type = "", + const real_type absolute_tolerance = 0, + const real_type relative_tolerance = 0, + const real_type step_tolerance = 0, + const int maximum_non_linear_iterations = -1, + const int max_n_function_evaluations = -1) + : options_prefix(options_prefix) + , snes_type(snes_type) + , snes_linesearch_type(snes_linesearch_type) + , absolute_tolerance(absolute_tolerance) + , relative_tolerance(relative_tolerance) + , step_tolerance(step_tolerance) + , maximum_non_linear_iterations(maximum_non_linear_iterations) + , max_n_function_evaluations(max_n_function_evaluations) + {} + + /** + * Import parameter values. + */ + void + add_parameters(ParameterHandler &prm); + + /** + * Options database prefix. + */ + std::string options_prefix; + + /** + * PETSc solver type. + */ + std::string snes_type; + + /** + * Linesearch type. + */ + std::string snes_linesearch_type; + + /** + * Absolute error tolerance for function evaluation. + * + * @note Non-positive values indicate to use PETSc's default. + */ + real_type absolute_tolerance; + + /** + * Relative error tolerance for function evaluation. + * + * @note Non-positive values indicate to use PETSc's default. + */ + real_type relative_tolerance; + + /** + * Step tolerance for solution update. + * + * @note Non-positive values indicate to use PETSc's default. + */ + real_type step_tolerance; + + /** + * Maximum number of nonlinear iterations allowed. + * + * @note Negative values indicate to use PETSc's default. + */ + int maximum_non_linear_iterations; + + /** + * Maximum number of function evaluations allowed. + * + * @note Negative values indicate to use PETSc's default. + */ + int max_n_function_evaluations; + }; + + /** + * Interface to PETSc SNES solver for nonlinear equations. + * The SNES solver is described in the + * [PETSc manual](https://petsc.org/release/manual/snes/). + * + * This class solves the nonlinear system of algebraic + * equations $F(x) = 0$. + * + * The interface to PETSc is realized by means of std::function callbacks + * like in the TrilinosWrappers::NOXSolver and SUNDIALS::KINSOL classes. + * + * NonlinearSolver supports any vector and matrix type having constructors and + * methods: + * + * @code + * class VectorType : public Subscriptor + * ... + * explicit VectorType(Vec); + * ... + * Vec & petsc_vector(); + * ... + * @endcode + * + * @code + * class MatrixType : public Subscriptor + * ... + * explicit MatrixType(Mat); + * ... + * Mat & petsc_matrix(); + * ... + * @endcode + * + * In particular, the supported types are the ones that can *wrap* + * PETSc's native vector and matrix classes, that are able to modify + * them in place, and that can return PETSc native types when requested. + * + * To use the solvers the user needs to provide the implementation of $F$ via + * the NonlinearSolver::residual callback. + * + * The default linearization procedure of a solver instantiated with + * this class consists in using Jacobian-Free-Newton-Krylov; the action of + * tangent matrices inside a linear solver process are approximated via + * matrix-free finite-differencing of the nonlinear residual equations. + * For details, consult the [PETSc + * manual](https://petsc.org/release/manual/snes/#sec-nlmatrixfree). + * + * In alternative, users can also provide the implementation of the + * Jacobian. This can be accomplished in two ways: + * - PETSc style using NonlinearSolver::jacobian + * - deal.II style using NonlinearSolver::setup_jacobian and + * NonlinearSolver::solve_with_jacobian. + * + * In case both approaches are implemented, the deal.II style + * will be used. + * + * The second approach is more in style with the deal.II philosophy + * and it also allows command line customization, like for example, + * @code + * ./myApp -snes_type newtontr -ksp_type cg + * @endcode + * in case the user wants to change the default nonlinear solver to + * a trust region solver and iterate on the tangent system with CG, + * still using NonlinearSolver::solve_with_jacobian as a preconditioner. + * + * The first approach has instead the advantage that only the matrix assembly + * procedure has to be provided, thus allowing quicker implementations and + * faster turnaround for experimenting with linear solver preconditioning + * configurations via command line customizations, like for example, + * @code + * ./myApp -ksp_type cg -pc_type gamg + * @endcode + * + * In case the nonlinear equations are derived from energy minimization + * arguments, it may be beneficial to perform linesearch or test trust-region + * model reductions using the energy functional. In such cases, users can + * set an optional NonlinearSolver::energy callback. + * + * @ingroup PETScWrappers + */ + template + class NonlinearSolver + { + public: + /** + * Type that holds real-valued numbers. + * + * Used to represent norms. + */ + using real_type = PetscReal; + + /** + * Constructor. + */ + NonlinearSolver(const NonlinearSolverData &data = NonlinearSolverData(), + const MPI_Comm & mpi_comm = PETSC_COMM_WORLD); + + /** + * Destructor. + */ + ~NonlinearSolver(); + + /** + * Conversion operator to gain access to the underlying PETSc type. If you + * do this, you cut this class off some information it may need, so this + * conversion operator should only be used if you know what you do. + */ + operator SNES() const; + + /** + * Return the PETSc SNES object. + */ + SNES + petsc_snes(); + + /** + * Return the underlying MPI communicator. + */ + MPI_Comm + get_mpi_communicator() const; + + /** + * Reset the solver. It does not change the customization. + */ + void + reinit(); + + /** + * Reset solver. Change customization according to @p data. + */ + void + reinit(const NonlinearSolverData &data); + + /** + * Set the preconditioning matrix only. + * + * When used with NonlinearSolver::setup_jacobian and + * NonlinearSolver::solve_with_jacobian, PETSc will approximate the + * linear system matrix-vector product using an internal matrix-free + * representation. + * + * When used with NonlinearSolver::jacobian PETSc will use the same matrix + * for both preconditioning and matrix-vector products. + */ + void + set_matrix(PMatrixType &P); + + /** + * Set both the linear system matrix and the preconditioning matrix + * that PETSc will use. + */ + void + set_matrices(AMatrixType &A, PMatrixType &P); + + /** + * Solve the nonlinear system of equations $F(x) = 0$. + * + * This function returns the number of iterations. + * The vector @p x must contain the initial guess. + * Upon returning, the @p x vector contains the solution. + */ + unsigned int + solve(VectorType &x); + + /** + * Solve the nonlinear system of equations $F(x) = 0$. + * + * This function returns the number of iterations. + * The vector @p x must contain the initial guess. + * Upon returning, the @p x vector contains the solution. + * + * Here we also set the matrix to precondition the tangent system. + */ + unsigned int + solve(VectorType &x, PMatrixType &P); + + /** + * Solve the nonlinear system of equations $F(x) = 0$. + * + * This function returns the number of iterations. + * The vector @p x must contain the initial guess. + * Upon returning, the @p x vector contains the solution. + * + * Here we also set the matrices to describe and precondition + * the tangent system. + */ + unsigned int + solve(VectorType &x, AMatrixType &A, PMatrixType &P); + + /** + * Callback for the computation of the nonlinear residual $F(x)$. + */ + std::function residual; + + /** + * Callback for the computation of the Jacobian + * $\dfrac{\partial F}{\partial x}$. + */ + std::function + jacobian; + + /** + * Callback for monitoring the solution process. + * + * This function is called by NonlinearSolver at the beginning + * of each time step. Input arguments are the current step number + * and the current value for ||F(x)||. + */ + std::function + monitor; + + /** + * Callback to set up the Jacobian system. + * + * This callback gives full control to users to set up the tangent + * operator $\dfrac{\partial F}{\partial x}$. + * + * Solvers must be provided via NonlinearSolver::solve_with_jacobian. + */ + std::function setup_jacobian; + + /** + * Callback for the solution of the tangent system set up with + * NonlinearSolver::setup_jacobian. + * + * This is used as a preconditioner inside the Krylov process. + */ + std::function + solve_with_jacobian; + + /** + * Callback for the computation of the energy function. + * + * This is usually not needed, since by default SNES assumes that the + * objective function to be minimized is $\frac{1}{2} || F(x) ||^2 $. + * + * However, if the nonlinear equations are derived from energy arguments, it + * may be useful to use this callback to perform linesearch or to test for + * the reduction in a trust region step. + * + * The value of the energy function must be returned in @p energy_value. + */ + std::function energy; + + private: + /** + * The PETSc SNES object. + */ + SNES snes; + + SmartPointer A; + SmartPointer P; + + /** + * This flag is used to support versions of PETSc older than 3.13. + */ + bool need_dummy_assemble; + }; + +} // namespace PETScWrappers + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PETSC + +#endif diff --git a/include/deal.II/lac/petsc_snes.templates.h b/include/deal.II/lac/petsc_snes.templates.h new file mode 100644 index 0000000000..dff02f5623 --- /dev/null +++ b/include/deal.II/lac/petsc_snes.templates.h @@ -0,0 +1,425 @@ +//----------------------------------------------------------- +// +// 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. +// +//--------------------------------------------------------------- + +#ifndef dealii_petsc_snes_templates_h +#define dealii_petsc_snes_templates_h + +#include + +#ifdef DEAL_II_WITH_PETSC +# include + +# include +# include + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +// Shorthand notation for PETSc error codes. +# define AssertPETSc(code) \ + do \ + { \ + PetscErrorCode ierr = (code); \ + AssertThrow(ierr == 0, ExcPETScError(ierr)); \ + } \ + while (0) + +// Shorthand notation for User error codes. +# define AssertUser(code, name) \ + do \ + { \ + int ierr = (code); \ + AssertThrow(ierr == 0, \ + StandardExceptions::ExcFunctionNonzeroReturn(name, ierr)); \ + } \ + while (0) + +namespace PETScWrappers +{ + template + NonlinearSolver::NonlinearSolver( + const NonlinearSolverData &data, + const MPI_Comm & mpi_comm) + { + AssertPETSc(SNESCreate(mpi_comm, &snes)); + AssertPETSc(SNESSetApplicationContext(snes, this)); + reinit(data); + } + + + + template + NonlinearSolver::operator SNES() const + { + return snes; + } + + + + template + SNES + NonlinearSolver::petsc_snes() + { + return snes; + } + + + + template + MPI_Comm + NonlinearSolver::get_mpi_communicator() + const + { + return PetscObjectComm(reinterpret_cast(snes)); + } + + + + template + NonlinearSolver::~NonlinearSolver() + { + AssertPETSc(SNESDestroy(&snes)); + } + + + + template + void + NonlinearSolver::reinit() + { + AssertPETSc(SNESReset(snes)); + } + + + + template + void + NonlinearSolver::reinit( + const NonlinearSolverData &data) + { + reinit(); + + // Solver type + if (data.snes_type.size()) + AssertPETSc(SNESSetType(snes, data.snes_type.c_str())); + + // Linesearch type + if (data.snes_linesearch_type.size()) + { + SNESLineSearch linesearch; + AssertPETSc(SNESGetLineSearch(snes, &linesearch)); + AssertPETSc( + SNESLineSearchSetType(linesearch, data.snes_linesearch_type.c_str())); + } + + // Options prefix + if (data.options_prefix.size()) + AssertPETSc(SNESSetOptionsPrefix(snes, data.options_prefix.c_str())); + + // Solver tolerances + PetscReal atol = data.absolute_tolerance > 0.0 ? + data.absolute_tolerance : + static_cast(PETSC_DEFAULT); + PetscReal rtol = data.relative_tolerance > 0.0 ? + data.relative_tolerance : + static_cast(PETSC_DEFAULT); + PetscReal stol = data.step_tolerance > 0.0 ? + data.step_tolerance : + static_cast(PETSC_DEFAULT); + + // Maximum number of iterations and function evaluations. + PetscInt maxit = data.maximum_non_linear_iterations >= 0 ? + data.maximum_non_linear_iterations : + static_cast(PETSC_DEFAULT); + PetscInt maxfe = data.max_n_function_evaluations >= 0 ? + data.max_n_function_evaluations : + static_cast(PETSC_DEFAULT); + AssertPETSc(SNESSetTolerances(snes, atol, rtol, stol, maxit, maxfe)); + } + + + + template + void + NonlinearSolver::set_matrix( + PMatrixType &P) + { + this->A = nullptr; + this->P = &P; + } + + + + template + void + NonlinearSolver::set_matrices( + AMatrixType &A, + PMatrixType &P) + { + this->A = &A; + this->P = &P; + } + + + + template + unsigned int + NonlinearSolver::solve(VectorType &x) + { + const auto snes_function = + [](SNES, Vec x, Vec f, void *ctx) -> PetscErrorCode { + PetscFunctionBeginUser; + auto user = static_cast(ctx); + + VectorType xdealii(x); + VectorType fdealii(f); + AssertUser(user->residual(xdealii, fdealii), "residual"); + petsc_increment_state_counter(f); + PetscFunctionReturn(0); + }; + + const auto snes_jacobian = + [](SNES, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode { + PetscFunctionBeginUser; + auto user = static_cast(ctx); + + VectorType xdealii(x); + AMatrixType Adealii(A); + PMatrixType Pdealii(P); + AssertUser(user->jacobian(xdealii, Adealii, Pdealii), "jacobian"); + petsc_increment_state_counter(P); + + // Handle the Jacobian-free case + // This call allows to resample the linearization point + // of the MFFD tangent operator + PetscBool flg; + AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); + if (flg) + { + AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); + } + else + petsc_increment_state_counter(A); + PetscFunctionReturn(0); + }; + + const auto snes_jacobian_with_setup = + [](SNES, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode { + PetscFunctionBeginUser; + auto user = static_cast(ctx); + + VectorType xdealii(x); + AMatrixType Adealii(A); + PMatrixType Pdealii(P); + + user->A = &Adealii; + user->P = &Pdealii; + AssertUser(user->setup_jacobian(xdealii), "setup_jacobian"); + petsc_increment_state_counter(P); + + // Handle older versions of PETSc for which we cannot pass a MATSHELL + // matrix to DMSetMatType. This has been fixed from 3.13 on. + if (user->need_dummy_assemble) + { + AssertPETSc(MatZeroEntries(P)); + AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); + } + + // Handle the Jacobian-free case + // This call allows to resample the linearization point + // of the MFFD tangent operator + PetscBool flg; + AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); + if (flg) + { + AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); + } + else + petsc_increment_state_counter(A); + PetscFunctionReturn(0); + }; + + const auto snes_monitor = + [](SNES snes, PetscInt it, PetscReal f, void *ctx) -> PetscErrorCode { + PetscFunctionBeginUser; + auto user = static_cast(ctx); + + Vec x; + AssertPETSc(SNESGetSolution(snes, &x)); + VectorType xdealii(x); + AssertUser(user->monitor(xdealii, it, f), "monitor"); + PetscFunctionReturn(0); + }; + + const auto snes_objective = + [](SNES, Vec x, PetscReal *f, void *ctx) -> PetscErrorCode { + PetscFunctionBeginUser; + auto user = static_cast(ctx); + + real_type v; + VectorType xdealii(x); + AssertUser(user->energy(xdealii, v), "energy"); + *f = v; + PetscFunctionReturn(0); + }; + + AssertThrow(residual, + StandardExceptions::ExcFunctionNotProvided("residual")); + + AssertPETSc(SNESSetSolution(snes, x.petsc_vector())); + + AssertPETSc(SNESSetFunction(snes, nullptr, snes_function, this)); + + if (energy) + AssertPETSc(SNESSetObjective(snes, snes_objective, this)); + + if (setup_jacobian) + { + AssertPETSc(SNESSetJacobian(snes, + A ? A->petsc_matrix() : nullptr, + P ? P->petsc_matrix() : nullptr, + snes_jacobian_with_setup, + this)); + if (!A) + set_use_matrix_free(snes, true, false); + + + // Do not waste memory by creating a dummy AIJ matrix inside PETSc. + this->need_dummy_assemble = false; + if (!P) + { +# if DEAL_II_PETSC_VERSION_GTE(3, 13, 0) + DM dm; + AssertPETSc(SNESGetDM(snes, &dm)); + AssertPETSc(DMSetMatType(dm, MATSHELL)); +# else + this->need_dummy_assemble = true; +# endif + } + } + else + { + if (jacobian) + { + AssertPETSc(SNESSetJacobian(snes, + A ? A->petsc_matrix() : + (P ? P->petsc_matrix() : nullptr), + P ? P->petsc_matrix() : nullptr, + snes_jacobian, + this)); + } + else + // The user did not set any Jacobian callback. PETSc default in this + // case is to use FD and thus assemble a dense operator by finite + // differencing the residual callbacks. Here instead we decide to + // use a full matrix-free approach by default. This choice can always + // be overriden from command line. + { + set_use_matrix_free(snes, false, true); + } + } + + // In case solve_with_jacobian is provided, create a shell + // preconditioner wrapping the user call. The internal Krylov + // solver will apply the preconditioner only once. This choice + // can be overriden by command line and users can use any other + // Krylov method if their solve is not accurate enough. + // Using solve_with_jacobian as a preconditioner allows users + // to provide approximate solvers and possibly iterate on a matrix-free + // approximation of the tangent operator. + PreconditionShell precond( + PetscObjectComm(reinterpret_cast(snes))); + if (solve_with_jacobian) + { + precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int { + VectorType dst(static_cast(indst)); + const VectorType src(static_cast(insrc)); + return solve_with_jacobian(src, dst); + }; + + // Default Krylov solver (preconditioner only) + KSP ksp; + AssertPETSc(SNESGetKSP(snes, &ksp)); + AssertPETSc(KSPSetType(ksp, KSPPREONLY)); + AssertPETSc(KSPSetPC(ksp, precond.get_pc())); + } + + // Attach user monitoring routine + if (monitor) + AssertPETSc(SNESMonitorSet(snes, snes_monitor, this, nullptr)); + + // Allow command line customization. + AssertPETSc(SNESSetFromOptions(snes)); + + // Having set everything up, now do the actual work + // and let PETSc solve the system. + // Older versions of PETSc requires the solution vector specified even + // if we specified SNESSetSolution upfront + AssertPETSc(SNESSolve(snes, nullptr, x.petsc_vector())); + + // Get the number of steps taken. + PetscInt nt; + AssertPETSc(SNESGetIterationNumber(snes, &nt)); + + // Raise an exception if the solver has not converged + SNESConvergedReason reason; + AssertPETSc(SNESGetConvergedReason(snes, &reason)); + AssertThrow(reason > 0, + ExcMessage("SNES solver did not converge after " + + std::to_string(nt) + " iterations with reason " + + SNESConvergedReasons[reason])); + + // Finally return + return nt; + } + + + + template + unsigned int + NonlinearSolver::solve(VectorType & x, + PMatrixType &P) + { + set_matrix(P); + return solve(x); + } + + + + template + unsigned int + NonlinearSolver::solve(VectorType & x, + AMatrixType &A, + PMatrixType &P) + { + set_matrices(A, P); + return solve(x); + } + +} // namespace PETScWrappers + +# undef AssertPETSc +# undef AssertUser + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_PETSC + +#endif diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 7be2a079a7..157dc9d834 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -108,6 +108,7 @@ if(DEAL_II_WITH_PETSC) petsc_parallel_sparse_matrix.cc petsc_parallel_vector.cc petsc_precondition.cc + petsc_snes.cc petsc_solver.cc petsc_ts.cc petsc_sparse_matrix.cc diff --git a/source/lac/petsc_snes.cc b/source/lac/petsc_snes.cc new file mode 100644 index 0000000000..065af3abaf --- /dev/null +++ b/source/lac/petsc_snes.cc @@ -0,0 +1,74 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef DOXYGEN + +# include + +# ifdef DEAL_II_WITH_PETSC + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace PETScWrappers +{ + void + NonlinearSolverData::add_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Running parameters"); + prm.add_parameter( + "options prefix", + options_prefix, + "The string indicating the options prefix for command line customization."); + prm.add_parameter("solver type", + snes_type, + "The string indicating the PETSc SNES type."); + prm.add_parameter("linesearch type", + snes_linesearch_type, + "The string indicating the PETSc linesearch type."); + prm.add_parameter("absolute error tolerance", + absolute_tolerance, + "Absolute error tolerance."); + prm.add_parameter("relative error tolerance", + relative_tolerance, + "Relative error tolerance."); + prm.add_parameter("step tolerance", step_tolerance, "Step tolerance."); + prm.add_parameter("maximum iterations", + maximum_non_linear_iterations, + "Maximum number of iterations allowed."); + prm.add_parameter("maximum function evaluations", + max_n_function_evaluations, + "Maximum number of function evaluations allowed."); + prm.leave_subsection(); + } + +} // namespace PETScWrappers + +template class PETScWrappers::NonlinearSolver<>; +template class PETScWrappers::NonlinearSolver; +template class PETScWrappers::NonlinearSolver; +template class PETScWrappers::NonlinearSolver; +template class PETScWrappers::NonlinearSolver< + PETScWrappers::MPI::BlockVector, + PETScWrappers::MPI::BlockSparseMatrix>; + + +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_WITH_PETSC +#endif diff --git a/tests/petsc/petsc_snes_00.cc b/tests/petsc/petsc_snes_00.cc new file mode 100644 index 0000000000..ddc9c38752 --- /dev/null +++ b/tests/petsc/petsc_snes_00.cc @@ -0,0 +1,206 @@ +//----------------------------------------------------------- +// +// 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 + +#include + +#include "../tests.h" + +/** + * Solves the nonlinear system of equations + * + * (x - y^3 + 1)^3 - y^3 = 0 + * x + 2y - 3 = 0 + * + * using the PETScWrappers::NonlinearSolver class + * that interfaces PETSc SNES solver object. + * + * Here we test a pure SNES interface and a deal.II approach + * when users take control of solving the linear systems. + */ +using VectorType = PETScWrappers::MPI::Vector; +using MatrixType = PETScWrappers::MatrixBase; +using Solver = PETScWrappers::NonlinearSolver; +using real_type = Solver::real_type; + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + std::ofstream out("output"); + + PETScWrappers::NonlinearSolverData data; + ParameterHandler prm; + + data.add_parameters(prm); + out << "# Default Parameters" << std::endl; + prm.print_parameters(out, ParameterHandler::ShortText); + + std::ifstream ifile(SOURCE_DIR "/petsc_snes_00_in.prm"); + prm.parse_input(ifile); + + out << "# Testing Parameters" << std::endl; + prm.print_parameters(out, ParameterHandler::ShortText); + + { + out << "# Test PETSc interface (MFFD)" << std::endl; + + // Create and customize the solver + Solver solver(data); + + // Here we only pass the callback for the function evaluation + // The tangent system will be approximated via matrix-free finite + // differencing. + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + auto x = X[0]; + auto y = X[1]; + F(0) = std::pow(x - std::pow(y, 3) + 1, 3) - std::pow(y, 3); + F(1) = x + 2 * y - 3; + F.compress(VectorOperation::insert); + return 0; + }; + + // Test attaching a user-defined monitoring routine + solver.monitor = [&](const VectorType & X, + const unsigned int step, + const real_type fval) -> int { + out << "# " << step << ": " << fval << std::endl; + return 0; + }; + + // Create and initialize solution vector + VectorType x(MPI_COMM_SELF, 2, 2); + x = 0.0; + + // Solve the nonlinear system (defaults to Newton with cubic-backtracking) + auto nit = solver.solve(x); + + out << "# Solution " << x[0] << ", " << x[1] << std::endl; + out << "# Iterations " << nit << std::endl; + } + + { + out << "# Test PETSc interface" << std::endl; + + Solver solver(data); + + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + auto x = X[0]; + auto y = X[1]; + F(0) = std::pow(x - std::pow(y, 3) + 1, 3) - std::pow(y, 3); + F(1) = x + 2 * y - 3; + F.compress(VectorOperation::insert); + return 0; + }; + + // Here we use the Jacobian callback following the PETSc style, + // where matrices are used as output arguments. + // This test uses A == P, so we do not do anything on A + // The intermediate callbacks used to pass from PETSc to deal.II + // will handle the case of users requesting Jacobian-free Newton + // Krylov (i.e. using -snes_mf_operator) + solver.jacobian = + [&](const VectorType &X, MatrixType &A, MatrixType &P) -> int { + auto x = X[0]; + auto y = X[1]; + auto f0_x = 3 * std::pow(x - std::pow(y, 3) + 1, 2); + auto f0_y = -9 * std::pow(x - std::pow(y, 3) + 1, 2) * std::pow(y, 2) - + 3 * std::pow(y, 2); + P.set(0, 0, f0_x); + P.set(0, 1, f0_y); + P.set(1, 0, 1); + P.set(1, 1, 2); + P.compress(VectorOperation::insert); + return 0; + }; + + VectorType x(MPI_COMM_SELF, 2, 2); + x = 0.0; + + // Solve the nonlinear system without specifying matrices + // This will use the default matrix type in PETSc for SNES + // (MATDENSE), which is ok for our test. + // See petsc_snes_01.cc for a case where we pass the matrix + auto nit = solver.solve(x); + + out << "# Solution " << x[0] << ", " << x[1] << std::endl; + out << "# Iterations " << nit << std::endl; + } + + { + out << "# Test user interface" << std::endl; + + Solver solver(data); + + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + auto x = X[0]; + auto y = X[1]; + F(0) = std::pow(x - std::pow(y, 3) + 1, 3) - std::pow(y, 3); + F(1) = x + 2 * y - 3; + F.compress(VectorOperation::insert); + return 0; + }; + + // When users want to be in full control of the linear system + // solves, they need to use setup_jacobian and solve_with_jacobian + // For example, in this case we compute the inverse of the Jacobian + // during setup_jacobian and we use it in the solve phase + FullMatrix Jinv(2, 2); + + solver.setup_jacobian = [&](const VectorType &X) -> int { + auto x = X[0]; + auto y = X[1]; + auto f0_x = 3 * std::pow(x - std::pow(y, 3) + 1, 2); + auto f0_y = -9 * std::pow(x - std::pow(y, 3) + 1, 2) * std::pow(y, 2) - + 3 * std::pow(y, 2); + FullMatrix J(2, 2); + J(0, 0) = f0_x; + J(0, 1) = f0_y; + J(1, 0) = 1; + J(1, 1) = 2; + Jinv.invert(J); + return 0; + }; + + // Solve phase. By default, PETSc will use this callback as a preconditioner + // within a preconditioner only Krylov solve. Other Krylov + // solvers can still be used in a Jacobian-free way and selected at command + // line or within user code. + solver.solve_with_jacobian = [&](const VectorType &src, + VectorType & dst) -> int { + dst(0) = Jinv(0, 0) * src(0) + Jinv(0, 1) * src(1); + dst(1) = Jinv(1, 0) * src(0) + Jinv(1, 1) * src(1); + dst.compress(VectorOperation::insert); + return 0; + }; + + VectorType x(MPI_COMM_SELF, 2, 2); + x = 0.0; + + // We don't need to pass matrices, since in this case we don't need them + auto nit = solver.solve(x); + + out << "# Solution " << x[0] << ", " << x[1] << std::endl; + out << "# Iterations " << nit << std::endl; + } + + return 0; +} diff --git a/tests/petsc/petsc_snes_00.output b/tests/petsc/petsc_snes_00.output new file mode 100644 index 0000000000..dd9aea4bec --- /dev/null +++ b/tests/petsc/petsc_snes_00.output @@ -0,0 +1,40 @@ +# Default Parameters +subsection Running parameters + set absolute error tolerance = 0 + set linesearch type = + set maximum function evaluations = -1 + set maximum iterations = -1 + set options prefix = + set relative error tolerance = 0 + set solver type = + set step tolerance = 0 +end +# Testing Parameters +subsection Running parameters + set absolute error tolerance = 1.e-9 + set linesearch type = + set maximum function evaluations = 77 + set maximum iterations = 23 + set options prefix = + set relative error tolerance = 2.e-9 + set solver type = newtonls + set step tolerance = 3.e-9 +end +# Test PETSc interface (MFFD) +# 0: 3.16228 +# 1: 2.84158 +# 2: 2.5438 +# 3: 2.27647 +# 4: 1.36494 +# 5: 0.288774 +# 6: 0.0203366 +# 7: 7.92738e-05 +# 8: 1.22106e-09 +# Solution 1, 1 +# Iterations 8 +# Test PETSc interface +# Solution 1, 1 +# Iterations 8 +# Test user interface +# Solution 1, 1 +# Iterations 8 diff --git a/tests/petsc/petsc_snes_00_in.prm b/tests/petsc/petsc_snes_00_in.prm new file mode 100644 index 0000000000..100e7479d4 --- /dev/null +++ b/tests/petsc/petsc_snes_00_in.prm @@ -0,0 +1,8 @@ +subsection Running parameters + set solver type = newtonls + set absolute error tolerance = 1.e-9 + set relative error tolerance = 2.e-9 + set step tolerance = 3.e-9 + set maximum iterations = 23 + set maximum function evaluations = 77 +end diff --git a/tests/petsc/petsc_snes_01.cc b/tests/petsc/petsc_snes_01.cc new file mode 100644 index 0000000000..916d70fe1d --- /dev/null +++ b/tests/petsc/petsc_snes_01.cc @@ -0,0 +1,182 @@ +//----------------------------------------------------------- +// +// 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 + +#include + +#include "../tests.h" + +/** + * Solves the nonlinear system of equations + * + * (x - y^3 + 1)^3 - y^3 = 0 + * x + 2y - 3 = 0 + * + * using the PETScWrappers::NonlinearSolver class + * that interfaces PETSc SNES solver object. + * + * This code tests the block interface. + * See petsc_snes_00.cc for additional information on the callbacks + */ +using VectorType = PETScWrappers::MPI::BlockVector; +using MatrixType = PETScWrappers::MPI::BlockSparseMatrix; +using Solver = PETScWrappers::NonlinearSolver; + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + std::ofstream out("output"); + + { + out << "# Test PETSc interface (MFFD)" << std::endl; + + Solver solver; + + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + auto x = X.block(0)[0]; + auto y = X.block(1)[0]; + + F.block(0)[0] = std::pow(x - std::pow(y, 3) + 1, 3) - std::pow(y, 3); + F.block(1)[0] = x + 2 * y - 3; + F.compress(VectorOperation::insert); + return 0; + }; + + VectorType x(2, MPI_COMM_SELF, 1, 1); + x = 0.0; + + auto nit = solver.solve(x); + + out << "# Solution " << x[0] << ", " << x[1] << std::endl; + out << "# Iterations " << nit << std::endl; + } + + { + out << "# Test PETSc interface" << std::endl; + + Solver solver; + + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + auto x = X.block(0)[0]; + auto y = X.block(1)[0]; + + F.block(0)[0] = std::pow(x - std::pow(y, 3) + 1, 3) - std::pow(y, 3); + F.block(1)[0] = x + 2 * y - 3; + F.compress(VectorOperation::insert); + return 0; + }; + + solver.jacobian = + [&](const VectorType &X, MatrixType &A, MatrixType &P) -> int { + auto x = X.block(0)[0]; + auto y = X.block(1)[0]; + auto f0_x = 3 * std::pow(x - std::pow(y, 3) + 1, 2); + auto f0_y = -9 * std::pow(x - std::pow(y, 3) + 1, 2) * std::pow(y, 2) - + 3 * std::pow(y, 2); + P.block(0, 0).set(0, 0, f0_x); + P.block(0, 1).set(0, 0, f0_y); + P.block(1, 0).set(0, 0, 1); + P.block(1, 1).set(0, 0, 2); + P.compress(VectorOperation::insert); + return 0; + }; + + VectorType x(2, MPI_COMM_SELF, 1, 1); + x = 0.0; + + // We need to create a representative matrix upfront + MatrixType J; + J.reinit(2, 2); + DynamicSparsityPattern csp(1, 1); + csp.add(0, 0); + IndexSet indices(1); + indices.add_range(0, 1); + for (unsigned int row = 0; row < 2; ++row) + for (unsigned int col = 0; col < 2; ++col) + J.block(row, col).reinit(indices, indices, csp, MPI_COMM_SELF); + J.collect_sizes(); + + // Solve the nonlinear system and specify the matrix instance to + // use. The matrix will be returned back to the user for assembly. + // Note that this is currently done by wrapping PETSc's Mat objects + // within light-weight on-the-fly constructors. + // Any additional information stored in non standard members of J + // will be lost, unless it can be deducted from the Mat object. + auto nit = solver.solve(x, J); + + out << "# Solution " << x[0] << ", " << x[1] << std::endl; + out << "# Iterations " << nit << std::endl; + } + + { + out << "# Test user interface" << std::endl; + + Solver solver; + + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + auto x = X.block(0)[0]; + auto y = X.block(1)[0]; + + F.block(0)[0] = std::pow(x - std::pow(y, 3) + 1, 3) - std::pow(y, 3); + F.block(1)[0] = x + 2 * y - 3; + F.compress(VectorOperation::insert); + return 0; + }; + + FullMatrix Jinv(2, 2); + + solver.setup_jacobian = [&](const VectorType &X) -> int { + auto x = X.block(0)[0]; + auto y = X.block(1)[0]; + auto f0_x = 3 * std::pow(x - std::pow(y, 3) + 1, 2); + auto f0_y = -9 * std::pow(x - std::pow(y, 3) + 1, 2) * std::pow(y, 2) - + 3 * std::pow(y, 2); + FullMatrix J(2, 2); + J(0, 0) = f0_x; + J(0, 1) = f0_y; + J(1, 0) = 1; + J(1, 1) = 2; + Jinv.invert(J); + return 0; + }; + + solver.solve_with_jacobian = [&](const VectorType &src, + VectorType & dst) -> int { + dst.block(0)[0] = + Jinv(0, 0) * src.block(0)[0] + Jinv(0, 1) * src.block(1)[0]; + dst.block(1)[0] = + Jinv(1, 0) * src.block(0)[0] + Jinv(1, 1) * src.block(1)[0]; + dst.compress(VectorOperation::insert); + return 0; + }; + + VectorType x(2, MPI_COMM_SELF, 1, 1); + x = 0.0; + + auto nit = solver.solve(x); + + out << "# Solution " << x[0] << ", " << x[1] << std::endl; + out << "# Iterations " << nit << std::endl; + } + + return 0; +} diff --git a/tests/petsc/petsc_snes_01.output b/tests/petsc/petsc_snes_01.output new file mode 100644 index 0000000000..ac952147a7 --- /dev/null +++ b/tests/petsc/petsc_snes_01.output @@ -0,0 +1,9 @@ +# Test PETSc interface (MFFD) +# Solution 1, 1 +# Iterations 8 +# Test PETSc interface +# Solution 1, 1 +# Iterations 8 +# Test user interface +# Solution 1, 1 +# Iterations 8 diff --git a/tests/petsc/petsc_snes_02.cc b/tests/petsc/petsc_snes_02.cc new file mode 100644 index 0000000000..35fc9759d4 --- /dev/null +++ b/tests/petsc/petsc_snes_02.cc @@ -0,0 +1,85 @@ +//----------------------------------------------------------- +// +// 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 "../tests.h" + +/** + * Tests user defined Vector and Matrix types + * and exceptions handling for PETSCWrappers::NonlinearSolver. + */ + +class VectorType : public Subscriptor +{ +public: + explicit VectorType(Vec v) + : v(v) + {} + + Vec & + petsc_vector() + { + return v; + } + +private: + Vec v; +}; + +class MatrixType : public Subscriptor +{ +public: + explicit MatrixType(Mat A) + : A(A) + {} + + Mat & + petsc_matrix() + { + return A; + } + +private: + Mat A; +}; + +using NonlinearSolver = PETScWrappers::NonlinearSolver; + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + VectorType v(nullptr); + MatrixType A(nullptr); + + NonlinearSolver mysolver; + + try + { + auto snes = mysolver.petsc_snes(); + AssertThrow(snes == static_cast(mysolver), ExcInternalError()); + mysolver.solve(v, A); + } + catch (std::exception &exc) + { + deallog << exc.what() << std::endl; + } + return 0; +} diff --git a/tests/petsc/petsc_snes_02.output b/tests/petsc/petsc_snes_02.output new file mode 100644 index 0000000000..7d1d25f3b3 --- /dev/null +++ b/tests/petsc/petsc_snes_02.output @@ -0,0 +1,11 @@ + +DEAL:: +-------------------------------------------------------- +An error occurred in file in function + unsigned int dealii::PETScWrappers::NonlinearSolver::solve(VectorType&) [with VectorType = VectorType; PMatrixType = MatrixType; AMatrixType = MatrixType] +The violated condition was: + residual +Additional information: + Please provide an implementation for the function "residual" +-------------------------------------------------------- + diff --git a/tests/petsc/petsc_snes_03.cc b/tests/petsc/petsc_snes_03.cc new file mode 100644 index 0000000000..b7d48a2671 --- /dev/null +++ b/tests/petsc/petsc_snes_03.cc @@ -0,0 +1,70 @@ +//----------------------------------------------------------- +// +// 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 + +#include + +#include "../tests.h" + +/** + * Solves the minimization problem + * + * min_x ||x||^2 + * + * using the PETScWrappers::NonlinearSolver class + * that interfaces PETSc SNES solver object. + * + * This code tests the objective function (energy) callback. + */ +using VectorType = PETScWrappers::MPI::Vector; +using Solver = PETScWrappers::NonlinearSolver; +using real_type = Solver::real_type; + +int +main(int argc, char **argv) +{ + initlog(); + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + { + Solver solver; + + solver.energy = [&](const VectorType &X, real_type &v) -> int { + v = X.norm_sqr(); + return 0; + }; + + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + F.equ(2, X); + return 0; + }; + + auto commsize = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + auto commrank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + VectorType x(MPI_COMM_WORLD, 10, commrank == commsize - 1 ? 10 : 0); + x = 1.0; + + auto nit = solver.solve(x); + + deallog << "Iterations " << nit << std::endl; + x.print(deallog.get_file_stream()); + } + return 0; +} diff --git a/tests/petsc/petsc_snes_03.mpirun=3.output b/tests/petsc/petsc_snes_03.mpirun=3.output new file mode 100644 index 0000000000..9f16c5063e --- /dev/null +++ b/tests/petsc/petsc_snes_03.mpirun=3.output @@ -0,0 +1,3 @@ + +DEAL::Iterations 1 +[Proc2 0-9] -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 diff --git a/tests/petsc/petsc_snes_03.output b/tests/petsc/petsc_snes_03.output new file mode 100644 index 0000000000..3b1cde2b0b --- /dev/null +++ b/tests/petsc/petsc_snes_03.output @@ -0,0 +1,3 @@ + +DEAL::Iterations 1 +[Proc0 0-9] -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 -8.426e-09 diff --git a/tests/petsc/petsc_snes_04.cc b/tests/petsc/petsc_snes_04.cc new file mode 100644 index 0000000000..9f09ff4424 --- /dev/null +++ b/tests/petsc/petsc_snes_04.cc @@ -0,0 +1,65 @@ +//----------------------------------------------------------- +// +// 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 + +#include + +#include "../tests.h" + +/** + * Tests exception raise for non convergence. + */ +using VectorType = PETScWrappers::MPI::Vector; +using Solver = PETScWrappers::NonlinearSolver; +using real_type = Solver::real_type; + +int +main(int argc, char **argv) +{ + initlog(); + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + { + PETScWrappers::NonlinearSolverData data; + data.maximum_non_linear_iterations = 0; + + Solver solver(data); + + solver.residual = [&](const VectorType &X, VectorType &F) -> int { + F.equ(2, X); + return 0; + }; + + auto commsize = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + auto commrank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + VectorType x(MPI_COMM_WORLD, 10, commrank == commsize - 1 ? 10 : 0); + x = 1.0; + + try + { + auto nit = solver.solve(x); + } + catch (std::exception &exc) + { + deallog << exc.what() << std::endl; + } + } + return 0; +} diff --git a/tests/petsc/petsc_snes_04.mpirun=3.output b/tests/petsc/petsc_snes_04.mpirun=3.output new file mode 100644 index 0000000000..f17430fb0a --- /dev/null +++ b/tests/petsc/petsc_snes_04.mpirun=3.output @@ -0,0 +1,12 @@ + +DEAL:: +-------------------------------------------------------- +An error occurred in file in function + unsigned int dealii::PETScWrappers::NonlinearSolver::solve(VectorType&) [with VectorType = dealii::PETScWrappers::MPI::Vector; PMatrixType = dealii::PETScWrappers::MatrixBase; AMatrixType = dealii::PETScWrappers::MatrixBase] +The violated condition was: + reason > 0 +Additional information: + SNES solver did not converge after 0 iterations with reason + DIVERGED_MAX_IT +-------------------------------------------------------- + diff --git a/tests/petsc/petsc_snes_04.output b/tests/petsc/petsc_snes_04.output new file mode 100644 index 0000000000..f17430fb0a --- /dev/null +++ b/tests/petsc/petsc_snes_04.output @@ -0,0 +1,12 @@ + +DEAL:: +-------------------------------------------------------- +An error occurred in file in function + unsigned int dealii::PETScWrappers::NonlinearSolver::solve(VectorType&) [with VectorType = dealii::PETScWrappers::MPI::Vector; PMatrixType = dealii::PETScWrappers::MatrixBase; AMatrixType = dealii::PETScWrappers::MatrixBase] +The violated condition was: + reason > 0 +Additional information: + SNES solver did not converge after 0 iterations with reason + DIVERGED_MAX_IT +-------------------------------------------------------- + -- 2.39.5