--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+#ifndef dealii_petsc_snes_h
+#define dealii_petsc_snes_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PETSC
+# include <deal.II/base/mpi.h>
+# include <deal.II/base/parameter_handler.h>
+# include <deal.II/base/smartpointer.h>
+
+# include <deal.II/lac/petsc_matrix_base.h>
+# include <deal.II/lac/petsc_precondition.h>
+# include <deal.II/lac/petsc_vector_base.h>
+
+# include <petscsnes.h>
+
+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 <typename VectorType = PETScWrappers::VectorBase,
+ typename PMatrixType = PETScWrappers::MatrixBase,
+ typename AMatrixType = PMatrixType>
+ 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<int(const VectorType &x, VectorType &res)> residual;
+
+ /**
+ * Callback for the computation of the Jacobian
+ * $\dfrac{\partial F}{\partial x}$.
+ */
+ std::function<int(const VectorType &x, AMatrixType &A, PMatrixType &P)>
+ 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<int(const VectorType & x,
+ const unsigned int step_number,
+ const real_type f_norm)>
+ 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<int(const VectorType &x)> 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<int(const VectorType &src, VectorType &dst)>
+ 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<int(const VectorType &x, real_type &energy_value)> energy;
+
+ private:
+ /**
+ * The PETSc SNES object.
+ */
+ SNES snes;
+
+ SmartPointer<AMatrixType, NonlinearSolver> A;
+ SmartPointer<PMatrixType, NonlinearSolver> 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
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+
+#ifndef dealii_petsc_snes_templates_h
+#define dealii_petsc_snes_templates_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PETSC
+# include <deal.II/base/exceptions.h>
+
+# include <deal.II/lac/petsc_precondition.h>
+# include <deal.II/lac/petsc_snes.h>
+
+# include <petscdm.h>
+# include <petscsnes.h>
+
+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 <typename VectorType, typename PMatrixType, typename AMatrixType>
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::NonlinearSolver(
+ const NonlinearSolverData &data,
+ const MPI_Comm & mpi_comm)
+ {
+ AssertPETSc(SNESCreate(mpi_comm, &snes));
+ AssertPETSc(SNESSetApplicationContext(snes, this));
+ reinit(data);
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::operator SNES() const
+ {
+ return snes;
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ SNES
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::petsc_snes()
+ {
+ return snes;
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ MPI_Comm
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::get_mpi_communicator()
+ const
+ {
+ return PetscObjectComm(reinterpret_cast<PetscObject>(snes));
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::~NonlinearSolver()
+ {
+ AssertPETSc(SNESDestroy(&snes));
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::reinit()
+ {
+ AssertPETSc(SNESReset(snes));
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::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<PetscReal>(PETSC_DEFAULT);
+ PetscReal rtol = data.relative_tolerance > 0.0 ?
+ data.relative_tolerance :
+ static_cast<PetscReal>(PETSC_DEFAULT);
+ PetscReal stol = data.step_tolerance > 0.0 ?
+ data.step_tolerance :
+ static_cast<PetscReal>(PETSC_DEFAULT);
+
+ // Maximum number of iterations and function evaluations.
+ PetscInt maxit = data.maximum_non_linear_iterations >= 0 ?
+ data.maximum_non_linear_iterations :
+ static_cast<PetscInt>(PETSC_DEFAULT);
+ PetscInt maxfe = data.max_n_function_evaluations >= 0 ?
+ data.max_n_function_evaluations :
+ static_cast<PetscInt>(PETSC_DEFAULT);
+ AssertPETSc(SNESSetTolerances(snes, atol, rtol, stol, maxit, maxfe));
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::set_matrix(
+ PMatrixType &P)
+ {
+ this->A = nullptr;
+ this->P = &P;
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ void
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::set_matrices(
+ AMatrixType &A,
+ PMatrixType &P)
+ {
+ this->A = &A;
+ this->P = &P;
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ unsigned int
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType &x)
+ {
+ const auto snes_function =
+ [](SNES, Vec x, Vec f, void *ctx) -> PetscErrorCode {
+ PetscFunctionBeginUser;
+ auto user = static_cast<NonlinearSolver *>(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<NonlinearSolver *>(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<NonlinearSolver *>(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<NonlinearSolver *>(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<NonlinearSolver *>(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<PetscObject>(snes)));
+ if (solve_with_jacobian)
+ {
+ precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int {
+ VectorType dst(static_cast<const Vec &>(indst));
+ const VectorType src(static_cast<const Vec &>(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 <typename VectorType, typename PMatrixType, typename AMatrixType>
+ unsigned int
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType & x,
+ PMatrixType &P)
+ {
+ set_matrix(P);
+ return solve(x);
+ }
+
+
+
+ template <typename VectorType, typename PMatrixType, typename AMatrixType>
+ unsigned int
+ NonlinearSolver<VectorType, PMatrixType, AMatrixType>::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
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef DOXYGEN
+
+# include <deal.II/base/config.h>
+
+# ifdef DEAL_II_WITH_PETSC
+
+# include <deal.II/lac/petsc_block_sparse_matrix.h>
+# include <deal.II/lac/petsc_snes.templates.h>
+
+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<PETScWrappers::MPI::Vector>;
+template class PETScWrappers::NonlinearSolver<PETScWrappers::MPI::BlockVector>;
+template class PETScWrappers::NonlinearSolver<PETScWrappers::MPI::Vector,
+ PETScWrappers::MPI::SparseMatrix>;
+template class PETScWrappers::NonlinearSolver<
+ PETScWrappers::MPI::BlockVector,
+ PETScWrappers::MPI::BlockSparseMatrix>;
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+# endif // DEAL_II_WITH_PETSC
+#endif
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_matrix_base.h>
+#include <deal.II/lac/petsc_snes.h>
+#include <deal.II/lac/petsc_vector.h>
+
+#include <cmath>
+
+#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<VectorType, MatrixType>;
+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<double> 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<double> 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;
+}
--- /dev/null
+# 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
--- /dev/null
+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
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_block_sparse_matrix.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_snes.h>
+
+#include <cmath>
+
+#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<VectorType, MatrixType>;
+
+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<double> 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<double> 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;
+}
--- /dev/null
+# 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
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/lac/exceptions.h>
+#include <deal.II/lac/petsc_snes.templates.h>
+
+#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<VectorType, MatrixType>;
+
+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<SNES>(mysolver), ExcInternalError());
+ mysolver.solve(v, A);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << exc.what() << std::endl;
+ }
+ return 0;
+}
--- /dev/null
+
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_snes.templates.h> in function
+ unsigned int dealii::PETScWrappers::NonlinearSolver<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = VectorType; PMatrixType = MatrixType; AMatrixType = MatrixType]
+The violated condition was:
+ residual
+Additional information:
+ Please provide an implementation for the function "residual"
+--------------------------------------------------------
+
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_block_sparse_matrix.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_snes.h>
+
+#include <cmath>
+
+#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<VectorType>;
+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;
+}
--- /dev/null
+
+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
--- /dev/null
+
+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
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_block_sparse_matrix.h>
+#include <deal.II/lac/petsc_block_vector.h>
+#include <deal.II/lac/petsc_snes.h>
+
+#include <cmath>
+
+#include "../tests.h"
+
+/**
+ * Tests exception raise for non convergence.
+ */
+using VectorType = PETScWrappers::MPI::Vector;
+using Solver = PETScWrappers::NonlinearSolver<VectorType>;
+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;
+}
--- /dev/null
+
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_snes.templates.h> in function
+ unsigned int dealii::PETScWrappers::NonlinearSolver<VectorType, PMatrixType, AMatrixType>::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
+--------------------------------------------------------
+
--- /dev/null
+
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_snes.templates.h> in function
+ unsigned int dealii::PETScWrappers::NonlinearSolver<VectorType, PMatrixType, AMatrixType>::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
+--------------------------------------------------------
+
projects / dealii.git / commitdiff