From: Stefano Zampini Date: Mon, 7 Aug 2023 09:50:52 +0000 (+0200) Subject: PETScWrappers:Align SNES with TS in comments and Jacobian handling X-Git-Tag: relicensing~610^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15849%2Fhead;p=dealii.git PETScWrappers:Align SNES with TS in comments and Jacobian handling --- diff --git a/include/deal.II/lac/petsc_snes.h b/include/deal.II/lac/petsc_snes.h index 7412b0ee02..02a64de044 100644 --- a/include/deal.II/lac/petsc_snes.h +++ b/include/deal.II/lac/petsc_snes.h @@ -30,6 +30,7 @@ # include # include + # if defined(DEAL_II_HAVE_CXX20) # include # endif @@ -207,19 +208,19 @@ namespace PETScWrappers * - deal.II style using NonlinearSolver::setup_jacobian and * NonlinearSolver::solve_with_jacobian. * The preconditioning matrix can be specified using - * NonlinearSolver::set_matrix. In case both approaches are implemented, the + * NonlinearSolver::set_matrix(). In case both approaches are implemented, the * deal.II style will be used. * - * NonlinearSolver::set_matrices must be used in case the user wants to + * NonlinearSolver::set_matrices() must be used in case the user wants to * provide the iteration matrix of the tangent system in the deal.II style * approach, thus replacing the matrix-free linearization. * * The correctness of the constructed Jacobians passed using - * NonlinearSolver::set_matrix can be checked using + * NonlinearSolver::set_matrix() can be checked using * @code * ./myApp -snes_test_jacobian * @endcode - * See NonlinearSolver::set_matrix and NonlinearSolver::set_matrices for + * See NonlinearSolver::set_matrix() and NonlinearSolver::set_matrices() for * additional details. * * The deal.II style approach still allows command line customization, like @@ -307,13 +308,14 @@ namespace PETScWrappers get_mpi_communicator() const; /** - * Reset the solver. It does not change the customization. + * Reset the solver, it does not change the customization. */ void reinit(); /** - * Reset solver. Change customization according to @p data. + * Reset solver. + * Change customization according to @p data. */ void reinit(const NonlinearSolverData &data); @@ -401,8 +403,8 @@ namespace PETScWrappers * 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)||. + * of each step. Input arguments are the current step number + * and the current value for $||F(x)||$. * * @note This variable represents a * @ref GlossUserProvidedCallBack "user provided callback". @@ -415,7 +417,7 @@ namespace PETScWrappers monitor; /** - * Callback to set up the Jacobian system. + * Callback for the set up of the Jacobian system. * * This callback gives full control to users to set up the tangent * operator $\dfrac{\partial F}{\partial x}$. diff --git a/include/deal.II/lac/petsc_snes.templates.h b/include/deal.II/lac/petsc_snes.templates.h index cbe8d09f26..a7a093a0e6 100644 --- a/include/deal.II/lac/petsc_snes.templates.h +++ b/include/deal.II/lac/petsc_snes.templates.h @@ -184,7 +184,7 @@ namespace PETScWrappers PMatrixType, Mat>)&&(concepts::is_dealii_petsc_matrix_type || std::constructible_from)) - MPI_Comm NonlinearSolver:: + inline MPI_Comm NonlinearSolver:: get_mpi_communicator() const { return PetscObjectComm(reinterpret_cast(snes)); @@ -416,12 +416,8 @@ namespace PETScWrappers PetscFunctionBeginUser; auto user = static_cast(ctx); - VectorType xdealii(x); - AMatrixType Adealii(A); - PMatrixType Pdealii(P); + VectorType xdealii(x); - user->A = &Adealii; - user->P = &Pdealii; const int lineno = __LINE__; const int err = call_and_possibly_capture_snes_exception( user->setup_jacobian, @@ -437,6 +433,13 @@ namespace PETScWrappers PETSC_ERR_LIB, PETSC_ERROR_INITIAL, "Failure in snes_jacobian_with_setup from dealii::PETScWrappers::NonlinearSolver"); + // The MatCopy calls below are 99% of the times dummy calls. + // They are only used in case we get different Mats then the one we passed + // to SNESSetJacobian. + if (user->P) + PetscCall(MatCopy(user->P->petsc_matrix(), P, SAME_NONZERO_PATTERN)); + if (user->A) + PetscCall(MatCopy(user->A->petsc_matrix(), A, SAME_NONZERO_PATTERN)); petsc_increment_state_counter(P); // Handle older versions of PETSc for which we cannot pass a MATSHELL @@ -474,12 +477,7 @@ namespace PETScWrappers VectorType xdealii(x); const int lineno = __LINE__; const int err = call_and_possibly_capture_snes_exception( - user->monitor, - user->pending_exception, - []() -> void {}, - xdealii, - it, - f); + user->monitor, user->pending_exception, {}, xdealii, it, f); if (err) return PetscError( PetscObjectComm((PetscObject)snes), @@ -540,10 +538,11 @@ namespace PETScWrappers P ? P->petsc_matrix() : nullptr, snes_jacobian_with_setup, this)); + + // Tell PETSc to set up a MFFD operator for the linear system matrix if (!A) set_use_matrix_free(snes, true, false); - // Do not waste memory by creating a dummy AIJ matrix inside PETSc. if (!P) { @@ -579,8 +578,8 @@ namespace PETScWrappers } // 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 + // preconditioner wrapping the user call. The default internal Krylov + // solver only applies the preconditioner. This choice // can be overridden 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 @@ -699,11 +698,13 @@ namespace PETScWrappers } } // namespace PETScWrappers + # undef AssertPETSc # if defined(undefPetscCall) # undef PetscCall # undef undefPetscCall # endif + DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_WITH_PETSC diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h index 050cb2e962..73d581a756 100644 --- a/include/deal.II/lac/petsc_ts.templates.h +++ b/include/deal.II/lac/petsc_ts.templates.h @@ -50,6 +50,7 @@ DEAL_II_NAMESPACE_OPEN CHKERRQ(ierr); \ } \ while (0) +# define undefPetscCall # endif namespace PETScWrappers @@ -796,7 +797,7 @@ namespace PETScWrappers user->need_dummy_assemble = false; // Handle the Jacobian-free case - // This call allow to resample the linearization point + // This call allows to resample the linearization point // of the MFFD tangent operator PetscBool flg; PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); @@ -1001,7 +1002,7 @@ namespace PETScWrappers // solver only applies the preconditioner. This choice // can be overridden 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 allow users + // 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. if (solve_with_jacobian)