From: Stefano Zampini Date: Sun, 16 Apr 2023 13:35:46 +0000 (+0300) Subject: PETScWrappers::TimeStepper fix backward compatibility to 3.7 X-Git-Tag: v9.5.0-rc1~299^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=41a0d407c8112bb0d8d9be24157fe1cd9334ad0b;p=dealii.git PETScWrappers::TimeStepper fix backward compatibility to 3.7 Align interface with SNES --- diff --git a/include/deal.II/lac/petsc_ts.h b/include/deal.II/lac/petsc_ts.h index a5f6153bae..3dc726407f 100644 --- a/include/deal.II/lac/petsc_ts.h +++ b/include/deal.II/lac/petsc_ts.h @@ -1,6 +1,6 @@ //----------------------------------------------------------- // -// Copyright (C) 2022 by the deal.II authors +// Copyright (C) 2023 by the deal.II authors // // This file is part of the deal.II library. // @@ -12,8 +12,6 @@ // the top level directory of deal.II. // //--------------------------------------------------------------- -// -// Author: Stefano Zampini, King Abdullah University of Science and Technology. #ifndef dealii_petsc_ts_h #define dealii_petsc_ts_h @@ -53,8 +51,8 @@ namespace PETScWrappers * * Running parameters: * - * @param opts_prefix The string indicating the options prefix for command line customization. - * @param tstype The string indicating the PETSc solver type. + * @param options_prefix The string indicating the options prefix for command line customization. + * @param ts_type The string indicating the PETSc solver type. * @param initial_time Initial simulation time. * @param final_time Final simulation time. * @param initial_step_size Initial step size. @@ -76,13 +74,15 @@ namespace PETScWrappers * Adaptive time stepping is disabled by default. * Negative values indicate using PETSc's default. * - * Note that all parameters values specified here can be overriden by + * @note All parameters values specified here can be overriden by * command line choices. + * + * @ingroup PETScWrappers */ TimeStepperData( // Running parameters - const std::string &opts_prefix = "", - const std::string &tstype = "", + const std::string &options_prefix = "", + const std::string &ts_type = "", const real_type initial_time = 0.0, const real_type final_time = 0.0, const real_type initial_step_size = 0.0, @@ -95,8 +95,8 @@ namespace PETScWrappers const real_type absolute_tolerance = -1.0, const real_type relative_tolerance = -1.0, const bool ignore_algebraic_lte = true) - : opts_prefix(opts_prefix) - , tstype(tstype) + : options_prefix(options_prefix) + , ts_type(ts_type) , initial_time(initial_time) , final_time(final_time) , initial_step_size(initial_step_size) @@ -110,68 +110,21 @@ namespace PETScWrappers , ignore_algebraic_lte(ignore_algebraic_lte) {} + /** + * Import parameter values. + */ void - add_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Running parameters"); - prm.add_parameter( - "options prefix", - opts_prefix, - "The string indicating the options prefix for command line customization."); - prm.add_parameter("solver type", - tstype, - "The string indicating the PETSc TS type."); - prm.add_parameter("initial time", - initial_time, - "The value for the initial time."); - prm.add_parameter("final time", - final_time, - "The value for the final time."); - prm.add_parameter("initial step size", - initial_step_size, - "The value for the initial time step."); - prm.add_parameter("maximum number of steps", - max_steps, - "Maximum number of time steps allowed."); - prm.add_parameter( - "match final time", - match_step, - "Whether or not to exactly stop at final time or step over it."); - prm.leave_subsection(); - - prm.enter_subsection("Error control"); - prm.add_parameter("adaptor type", - tsadapttype, - "The string for the TSAdapt type."); - prm.add_parameter("minimum step size", - minimum_step_size, - "Minimum time step size allowed."); - prm.add_parameter("maximum step size", - maximum_step_size, - "Maximum time step size allowed."); - prm.add_parameter("absolute error tolerance", - absolute_tolerance, - "Absolute error tolerance."); - prm.add_parameter("relative error tolerance", - relative_tolerance, - "Absolute error tolerance."); - prm.add_parameter( - "ignore algebraic lte", - ignore_algebraic_lte, - "Indicate whether or not to suppress algebraic variables " - "in the local truncation error test."); - prm.leave_subsection(); - } + add_parameters(ParameterHandler &prm); /** * Options database prefix. */ - std::string opts_prefix; + std::string options_prefix; /** * PETSc solver type. */ - std::string tstype; + std::string ts_type; /** * Initial time for the DAE. @@ -185,18 +138,20 @@ namespace PETScWrappers /** * Initial step size. - * Non-positive values are ignored. + * + * @note Non-positive values are ignored. */ real_type initial_step_size; /** * Maximum number of steps to be taken. - * Negative values are ignored. + * + * @note Negative values are ignored. */ int max_steps; /** - * Match final time requested? + * Flag to indicate to stop exactly at the requested final time. */ bool match_step; @@ -207,25 +162,29 @@ namespace PETScWrappers /** * Minimum allowed step size for adaptive time stepping. - * Non-positive values indicate to use PETSc's default. + * + * @note Non-positive values indicate to use PETSc's default. */ real_type minimum_step_size; /** * Maximum allowed step size for adaptive time stepping. - * Non-positive values indicate to use PETSc's default. + * + * @note Non-positive values indicate to use PETSc's default. */ real_type maximum_step_size; /** * Absolute error tolerance for adaptive time stepping. - * Negative values indicate to use PETSc's default. + * + * @note Negative values indicate to use PETSc's default. */ real_type absolute_tolerance; /** * Relative error tolerance for adaptive time stepping. - * Negative values indicate to use PETSc's default. + * + * @note Negative values indicate to use PETSc's default. */ real_type relative_tolerance; @@ -300,10 +259,10 @@ namespace PETScWrappers * * In alternative, users can also provide the implementations of the * *Jacobians*. This can be accomplished in two ways: - * - a-la-PETSc style using TimeStepper::implicit_jacobian + * - PETSc style using TimeStepper::implicit_jacobian * and TimeStepper::explicit_jacobian. - * - a-la-Deal.II style using TimeStepper::setup_jacobian and - * TimeStepper::solve_for_jacobian_system + * - deal.II style using TimeStepper::setup_jacobian and + * TimeStepper::solve_with_jacobian * * In case both approaches are coded, the deal.II style * will be used. @@ -315,7 +274,7 @@ namespace PETScWrappers * @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 TimeStepper::solve_for_jacobian_system as a preconditioner. + * still using TimeStepper::solve_with_jacobian as a preconditioner. * * The first approach has instead the advantage that only the matrix assembly * procedure has to be coded, thus allowing quicker implementations and @@ -324,6 +283,8 @@ namespace PETScWrappers * @code * ./myApp -ksp_type cg -pc_type gamg * @endcode + * + * @ingroup PETScWrappers */ template - solve_for_jacobian_system; + solve_with_jacobian; /** - * Return an index set containing the algebraic components. + * Callback to return an index set containing the algebraic components. * * Implementation of this function is optional. If your equation is also * algebraic (i.e., it contains algebraic constraints, or Lagrange @@ -533,55 +503,25 @@ namespace PETScWrappers */ std::function algebraic_components; - protected: + private: /** - * The PETSc object + * The PETSc object. */ TS ts; SmartPointer A; SmartPointer P; - bool need_dae_tolerances; - }; -# ifndef DOXYGEN - /* Inline functions */ - - template - inline TS - TimeStepper::petsc_ts() - { - return ts; - } - - template - inline MPI_Comm - TimeStepper::get_mpi_communicator() - const - { - return PetscObjectComm(reinterpret_cast(ts)); - } + /** + * This flag is set when changing the customization and used within solve. + */ + bool need_dae_tolerances; - template - inline typename TimeStepper::real_type - TimeStepper::get_time() - { - PetscReal t; - PetscErrorCode ierr = TSGetTime(ts, &t); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - return t; - } - - template - inline typename TimeStepper::real_type - TimeStepper::get_time_step() - { - PetscReal dt; - PetscErrorCode ierr = TSGetTimeStep(ts, &dt); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - return dt; - } -# endif + /** + * This flag is used to support versions of PETSc older than 3.13. + */ + bool need_dummy_assemble; + }; } // namespace PETScWrappers diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h index 93c0fb09c2..91d3268c64 100644 --- a/include/deal.II/lac/petsc_ts.templates.h +++ b/include/deal.II/lac/petsc_ts.templates.h @@ -1,6 +1,6 @@ //----------------------------------------------------------- // -// Copyright (C) 2022 by the deal.II authors +// Copyright (C) 2023 by the deal.II authors // // This file is part of the deal.II library. // @@ -12,8 +12,6 @@ // the top level directory of deal.II. // //--------------------------------------------------------------- -// -// Author: Stefano Zampini, King Abdullah University of Science and Technology. #ifndef dealii_petsc_ts_templates_h #define dealii_petsc_ts_templates_h @@ -32,7 +30,7 @@ DEAL_II_NAMESPACE_OPEN // Shorthand notation for PETSc error codes. -# define AssertTS(code) \ +# define AssertPETSc(code) \ do \ { \ PetscErrorCode ierr = (code); \ @@ -57,8 +55,8 @@ namespace PETScWrappers const TimeStepperData &data, const MPI_Comm & mpi_comm) { - AssertTS(TSCreate(mpi_comm, &ts)); - AssertTS(TSSetApplicationContext(ts, this)); + AssertPETSc(TSCreate(mpi_comm, &ts)); + AssertPETSc(TSSetApplicationContext(ts, this)); reinit(data); } @@ -67,7 +65,58 @@ namespace PETScWrappers template TimeStepper::~TimeStepper() { - AssertTS(TSDestroy(&ts)); + AssertPETSc(TSDestroy(&ts)); + } + + + + template + TimeStepper::operator TS() const + { + return ts; + } + + + + template + TS + TimeStepper::petsc_ts() + { + return ts; + } + + + + template + inline MPI_Comm + TimeStepper::get_mpi_communicator() + const + { + return PetscObjectComm(reinterpret_cast(ts)); + } + + + + template + typename TimeStepper::real_type + TimeStepper::get_time() + { + PetscReal t; + PetscErrorCode ierr = TSGetTime(ts, &t); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + return t; + } + + + + template + typename TimeStepper::real_type + TimeStepper::get_time_step() + { + PetscReal dt; + PetscErrorCode ierr = TSGetTimeStep(ts, &dt); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + return dt; } @@ -76,7 +125,7 @@ namespace PETScWrappers void TimeStepper::reinit() { - AssertTS(TSReset(ts)); + AssertPETSc(TSReset(ts)); } @@ -89,27 +138,28 @@ namespace PETScWrappers reinit(); // Solver type - if (data.tstype.size()) - AssertTS(TSSetType(ts, data.tstype.c_str())); + if (data.ts_type.size()) + AssertPETSc(TSSetType(ts, data.ts_type.c_str())); // Options prefix - if (data.opts_prefix.size()) - AssertTS(TSSetOptionsPrefix(ts, data.opts_prefix.c_str())); + if (data.options_prefix.size()) + AssertPETSc(TSSetOptionsPrefix(ts, data.options_prefix.c_str())); // Time and steps limits - AssertTS(TSSetTime(ts, data.initial_time)); + AssertPETSc(TSSetTime(ts, data.initial_time)); if (data.final_time > data.initial_time) - AssertTS(TSSetMaxTime(ts, data.final_time)); + ts_set_max_time(ts, data.final_time); if (data.initial_step_size > 0.0) - AssertTS(TSSetTimeStep(ts, data.initial_step_size)); + AssertPETSc(TSSetTimeStep(ts, data.initial_step_size)); if (data.max_steps >= 0) - AssertTS(TSSetMaxSteps(ts, data.max_steps)); + ts_set_max_steps(ts, data.max_steps); // Decide how to end the integration. Either stepover the final time or // match it. - AssertTS(TSSetExactFinalTime(ts, - data.match_step ? TS_EXACTFINALTIME_MATCHSTEP : - TS_EXACTFINALTIME_STEPOVER)); + AssertPETSc(TSSetExactFinalTime(ts, + data.match_step ? + TS_EXACTFINALTIME_MATCHSTEP : + TS_EXACTFINALTIME_STEPOVER)); // Adaptive tolerances const PetscReal atol = data.absolute_tolerance > 0.0 ? @@ -118,7 +168,7 @@ namespace PETScWrappers const PetscReal rtol = data.relative_tolerance > 0.0 ? data.relative_tolerance : static_cast(PETSC_DEFAULT); - AssertTS(TSSetTolerances(ts, atol, nullptr, rtol, nullptr)); + AssertPETSc(TSSetTolerances(ts, atol, nullptr, rtol, nullptr)); // At this point we do not know the problem size so we cannot // set variable tolerances for differential and algebratic equations @@ -127,13 +177,14 @@ namespace PETScWrappers // Adaptive time stepping TSAdapt tsadapt; - AssertTS(TSGetAdapt(ts, &tsadapt)); - AssertTS(TSAdaptSetType(tsadapt, data.tsadapttype.c_str())); + AssertPETSc(TSGetAdapt(ts, &tsadapt)); + AssertPETSc(TSAdaptSetType(tsadapt, data.tsadapttype.c_str())); // As of 3.19, PETSc does not propagate options prefixes to the // adaptors. - if (data.opts_prefix.size()) - AssertTS(TSAdaptSetOptionsPrefix(tsadapt, data.opts_prefix.c_str())); + if (data.options_prefix.size()) + AssertPETSc( + TSAdaptSetOptionsPrefix(tsadapt, data.options_prefix.c_str())); // Time step limits const PetscReal hmin = data.minimum_step_size > 0.0 ? @@ -142,15 +193,14 @@ namespace PETScWrappers const PetscReal hmax = data.maximum_step_size > 0.0 ? data.maximum_step_size : static_cast(PETSC_DEFAULT); - AssertTS(TSAdaptSetStepLimits(tsadapt, hmin, hmax)); + AssertPETSc(TSAdaptSetStepLimits(tsadapt, hmin, hmax)); } template void - TimeStepper::reinit_matrix( - PMatrixType &P) + TimeStepper::set_matrix(PMatrixType &P) { this->A = nullptr; this->P = &P; @@ -160,7 +210,7 @@ namespace PETScWrappers template void - TimeStepper::reinit_matrices( + TimeStepper::set_matrices( AMatrixType &A, PMatrixType &P) { @@ -174,11 +224,9 @@ namespace PETScWrappers unsigned int TimeStepper::solve(VectorType &y) { - auto ts_ifunction_ = - [](TS ts, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx) - -> PetscErrorCode { + const auto ts_ifunction = + [](TS, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx) -> PetscErrorCode { PetscFunctionBeginUser; - (void)ts; auto user = static_cast(ctx); VectorType xdealii(x); @@ -186,19 +234,14 @@ namespace PETScWrappers VectorType fdealii(f); AssertUser(user->implicit_function(t, xdealii, xdotdealii, fdealii), "implicit_function"); + petsc_increment_state_counter(f); PetscFunctionReturn(0); }; - auto ts_ijacobian_ = [](TS ts, - PetscReal t, - Vec x, - Vec xdot, - PetscReal s, - Mat A, - Mat P, - void * ctx) -> PetscErrorCode { + const auto ts_ijacobian = + [](TS, PetscReal t, Vec x, Vec xdot, PetscReal s, Mat A, Mat P, void *ctx) + -> PetscErrorCode { PetscFunctionBeginUser; - (void)ts; auto user = static_cast(ctx); VectorType xdealii(x); @@ -209,30 +252,27 @@ namespace PETScWrappers AssertUser( user->implicit_jacobian(t, xdealii, xdotdealii, s, Adealii, Pdealii), "implicit_jacobian"); + petsc_increment_state_counter(P); // Handle the Jacobian-free case // This call allow to resample the linearization point // of the MFFD tangent operator PetscBool flg; - AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); + AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); if (flg) { - AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); - AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); } + else + petsc_increment_state_counter(A); PetscFunctionReturn(0); }; - auto ts_ijacobian_with_setup_ = [](TS ts, - PetscReal t, - Vec x, - Vec xdot, - PetscReal s, - Mat A, - Mat P, - void * ctx) -> PetscErrorCode { + const auto ts_ijacobian_with_setup = + [](TS, PetscReal t, Vec x, Vec xdot, PetscReal s, Mat A, Mat P, void *ctx) + -> PetscErrorCode { PetscFunctionBeginUser; - (void)ts; auto user = static_cast(ctx); VectorType xdealii(x); @@ -244,24 +284,38 @@ namespace PETScWrappers user->P = &Pdealii; AssertUser(user->setup_jacobian(t, xdealii, xdotdealii, s), "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. + // What we need to do for older version of PETSc is instead to have + // a zero matrix with all diagonal entries present. + if (user->need_dummy_assemble) + { + AssertPETSc(MatZeroEntries(P)); + AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatShift(P, 0.0)); + } // Handle the Jacobian-free case // This call allow to resample the linearization point // of the MFFD tangent operator PetscBool flg; - AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); + AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); if (flg) { - AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); - AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); } + else + petsc_increment_state_counter(A); PetscFunctionReturn(0); }; - auto ts_rhsfunction_ = - [](TS ts, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode { + const auto ts_rhsfunction = + [](TS, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode { PetscFunctionBeginUser; - (void)ts; auto user = static_cast(ctx); VectorType xdealii(x); @@ -269,13 +323,13 @@ namespace PETScWrappers AssertUser(user->explicit_function(t, xdealii, fdealii), "explicit_function"); + petsc_increment_state_counter(f); PetscFunctionReturn(0); }; - auto ts_rhsjacobian_ = - [](TS ts, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode { + const auto ts_rhsjacobian = + [](TS, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode { PetscFunctionBeginUser; - (void)ts; auto user = static_cast(ctx); VectorType xdealii(x); @@ -284,24 +338,36 @@ namespace PETScWrappers AssertUser(user->explicit_jacobian(t, xdealii, Adealii, Pdealii), "explicit_jacobian"); + petsc_increment_state_counter(P); + + // Handle older versions of PETSc for which we cannot pass a MATSHELL + // matrix to DMSetMatType + if (user->need_dummy_assemble) + { + AssertPETSc(MatZeroEntries(P)); + AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatShift(P, 0.0)); + } // Handle the Jacobian-free case // This call allow to resample the linearization point // of the MFFD tangent operator PetscBool flg; - AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); + AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg)); if (flg) { - AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); - AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); + AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); } + else + petsc_increment_state_counter(A); PetscFunctionReturn(0); }; - auto ts_monitor_ = - [](TS ts, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode { + const auto ts_monitor = + [](TS, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode { PetscFunctionBeginUser; - (void)ts; auto user = static_cast(ctx); VectorType xdealii(x); @@ -313,56 +379,61 @@ namespace PETScWrappers StandardExceptions::ExcFunctionNotProvided( "explicit_function || implicit_function")); - AssertTS(TSSetSolution(ts, y.petsc_vector())); + AssertPETSc(TSSetSolution(ts, y.petsc_vector())); if (explicit_function) - AssertTS(TSSetRHSFunction(ts, nullptr, ts_rhsfunction_, this)); + AssertPETSc(TSSetRHSFunction(ts, nullptr, ts_rhsfunction, this)); if (implicit_function) - AssertTS(TSSetIFunction(ts, nullptr, ts_ifunction_, this)); + AssertPETSc(TSSetIFunction(ts, nullptr, ts_ifunction, this)); if (setup_jacobian) { - AssertTS(TSSetIJacobian(ts, - A ? A->petsc_matrix() : nullptr, - P ? P->petsc_matrix() : nullptr, - ts_ijacobian_with_setup_, - this)); + AssertPETSc(TSSetIJacobian(ts, + A ? A->petsc_matrix() : nullptr, + P ? P->petsc_matrix() : nullptr, + ts_ijacobian_with_setup, + this)); // Tell PETSc to setup a MFFD operator for the linear system matrix - SNES snes; - AssertTS(TSGetSNES(ts, &snes)); if (!A) - AssertTS(SNESSetUseMatrixFree(snes, PETSC_TRUE, PETSC_FALSE)); + set_use_matrix_free(ts, true, false); // Do not waste memory by creating a dummy AIJ matrix inside PETSc. + this->need_dummy_assemble = false; if (!P) { - DM dm; - AssertTS(SNESGetDM(snes, &dm)); - AssertTS(DMSetMatType(dm, MATSHELL)); +# if DEAL_II_PETSC_VERSION_GTE(3, 13, 0) + DM dm; + SNES snes; + AssertPETSc(TSGetSNES(ts, &snes)); + AssertPETSc(SNESGetDM(snes, &dm)); + AssertPETSc(DMSetMatType(dm, MATSHELL)); +# else + this->need_dummy_assemble = true; +# endif } } else { if (explicit_jacobian) { - AssertTS(TSSetRHSJacobian(ts, - A ? A->petsc_matrix() : - (P ? P->petsc_matrix() : nullptr), - P ? P->petsc_matrix() : nullptr, - ts_rhsjacobian_, - this)); + AssertPETSc(TSSetRHSJacobian(ts, + A ? A->petsc_matrix() : + (P ? P->petsc_matrix() : nullptr), + P ? P->petsc_matrix() : nullptr, + ts_rhsjacobian, + this)); } if (implicit_jacobian) { - AssertTS(TSSetIJacobian(ts, - A ? A->petsc_matrix() : - (P ? P->petsc_matrix() : nullptr), - P ? P->petsc_matrix() : nullptr, - ts_ijacobian_, - this)); + AssertPETSc(TSSetIJacobian(ts, + A ? A->petsc_matrix() : + (P ? P->petsc_matrix() : nullptr), + P ? P->petsc_matrix() : nullptr, + ts_ijacobian, + this)); } // The user did not set any Jacobian callback. PETSc default in this @@ -372,55 +443,53 @@ namespace PETScWrappers // be overriden from command line. if (!explicit_jacobian && !implicit_jacobian) { - SNES snes; - AssertTS(TSGetSNES(ts, &snes)); - AssertTS(SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE)); + set_use_matrix_free(ts, false, true); } } - // In case solve_for_jacobian_system is provided, create a shell + // In case solve_with_jacobian is provided, create a shell // preconditioner wrapping the user call. The default internal Krylov // solver only applies the preconditioner. 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_for_jacobian_system as a preconditioner allow users + // Using solve_with_jacobian as a preconditioner allow users // to provide approximate solvers and possibly iterate on a matrix-free // approximation of the tangent operator. PreconditionShell precond( PetscObjectComm(reinterpret_cast(ts))); - if (solve_for_jacobian_system) + 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_for_jacobian_system(src, dst); + return solve_with_jacobian(src, dst); }; // Default Krylov solver (preconditioner only) SNES snes; KSP ksp; - AssertTS(TSGetSNES(ts, &snes)); - AssertTS(SNESGetKSP(snes, &ksp)); - AssertTS(KSPSetType(ksp, KSPPREONLY)); - AssertTS(KSPSetPC(ksp, precond.get_pc())); + AssertPETSc(TSGetSNES(ts, &snes)); + AssertPETSc(SNESGetKSP(snes, &ksp)); + AssertPETSc(KSPSetType(ksp, KSPPREONLY)); + AssertPETSc(KSPSetPC(ksp, precond.get_pc())); } // Attach user monitoring routine. if (monitor) - AssertTS(TSMonitorSet(ts, ts_monitor_, this, nullptr)); + AssertPETSc(TSMonitorSet(ts, ts_monitor, this, nullptr)); // Allow command line customization. - AssertTS(TSSetFromOptions(ts)); + AssertPETSc(TSSetFromOptions(ts)); // Handle algebraic components. if (algebraic_components && need_dae_tolerances) { PetscReal atol, rtol; - AssertTS(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr)); + AssertPETSc(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr)); Vec av, rv; - AssertTS(VecDuplicate(y.petsc_vector(), &av)); - AssertTS(VecDuplicate(y.petsc_vector(), &rv)); + AssertPETSc(VecDuplicate(y.petsc_vector(), &av)); + AssertPETSc(VecDuplicate(y.petsc_vector(), &rv)); VectorBase avdealii(av); VectorBase rvdealii(rv); @@ -433,18 +502,27 @@ namespace PETScWrappers } avdealii.compress(VectorOperation::insert); rvdealii.compress(VectorOperation::insert); - AssertTS(TSSetTolerances(ts, atol, av, rtol, rv)); - AssertTS(VecDestroy(&av)); - AssertTS(VecDestroy(&rv)); + AssertPETSc(TSSetTolerances(ts, atol, av, rtol, rv)); + AssertPETSc(VecDestroy(&av)); + AssertPETSc(VecDestroy(&rv)); } // Having set everything up, now do the actual work // and let PETSc do the time stepping. - AssertTS(TSSolve(ts, nullptr)); + AssertPETSc(TSSolve(ts, nullptr)); + + // Get the number of steps taken. + auto nt = ts_get_step_number(ts); + + // Raise an exception if the solver has not converged + TSConvergedReason reason; + AssertPETSc(TSGetConvergedReason(ts, &reason)); + AssertThrow(reason > 0, + ExcMessage("TS solver did not converge after " + + std::to_string(nt) + " iterations with reason " + + TSConvergedReasons[reason])); - // Return the number of steps taken. - PetscInt nt; - AssertTS(TSGetStepNumber(ts, &nt)); + // Finally return return nt; } @@ -455,7 +533,7 @@ namespace PETScWrappers TimeStepper::solve(VectorType & y, PMatrixType &P) { - reinit_matrix(P); + set_matrix(P); return solve(y); } @@ -467,13 +545,13 @@ namespace PETScWrappers AMatrixType &A, PMatrixType &P) { - reinit_matrices(A, P); + set_matrices(A, P); return solve(y); } } // namespace PETScWrappers -# undef AssertTS +# undef AssertPETSc # undef AssertUser DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/petsc_ts.cc b/source/lac/petsc_ts.cc index d7018087d2..f54cdc0b20 100644 --- a/source/lac/petsc_ts.cc +++ b/source/lac/petsc_ts.cc @@ -24,6 +24,61 @@ DEAL_II_NAMESPACE_OPEN +namespace PETScWrappers +{ + void + TimeStepperData::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", + ts_type, + "The string indicating the PETSc TS type."); + prm.add_parameter("initial time", + initial_time, + "The value for the initial time."); + prm.add_parameter("final time", + final_time, + "The value for the final time."); + prm.add_parameter("initial step size", + initial_step_size, + "The value for the initial time step."); + prm.add_parameter("maximum number of steps", + max_steps, + "Maximum number of time steps allowed."); + prm.add_parameter( + "match final time", + match_step, + "Whether or not to exactly stop at final time or step over it."); + prm.leave_subsection(); + + prm.enter_subsection("Error control"); + prm.add_parameter("adaptor type", + tsadapttype, + "The string for the TSAdapt type."); + prm.add_parameter("minimum step size", + minimum_step_size, + "Minimum time step size allowed."); + prm.add_parameter("maximum step size", + maximum_step_size, + "Maximum time step size allowed."); + prm.add_parameter("absolute error tolerance", + absolute_tolerance, + "Absolute error tolerance."); + prm.add_parameter("relative error tolerance", + relative_tolerance, + "Absolute error tolerance."); + prm.add_parameter("ignore algebraic lte", + ignore_algebraic_lte, + "Indicate whether or not to suppress algebraic variables " + "in the local truncation error test."); + prm.leave_subsection(); + } + +} // namespace PETScWrappers template class PETScWrappers::TimeStepper<>; template class PETScWrappers::TimeStepper; diff --git a/tests/petsc/petsc_ts_00.cc b/tests/petsc/petsc_ts_00.cc index d95ae8aad6..98f197bc0d 100644 --- a/tests/petsc/petsc_ts_00.cc +++ b/tests/petsc/petsc_ts_00.cc @@ -1,6 +1,6 @@ //----------------------------------------------------------- // -// Copyright (C) 2022 by the deal.II authors +// Copyright (C) 2023 by the deal.II authors // // This file is part of the deal.II library. // @@ -12,8 +12,6 @@ // the top level directory of deal.II. // //----------------------------------------------------------- -// -// Author: Stefano Zampini, King Abdullah University of Science and Technology. #include @@ -136,8 +134,8 @@ public: // In the solve phase we se the stored shift to solve // for the implicit Jacobian system - time_stepper.solve_for_jacobian_system = - [&](const VectorType &src, VectorType &dst) -> int { + time_stepper.solve_with_jacobian = [&](const VectorType &src, + VectorType &dst) -> int { auto sf = 1. / (kappa * kappa + myshift * myshift); dst(0) = sf * (myshift * src(0) + src(1)); dst(1) = sf * (-kappa * kappa * src(0) + myshift * src(1)); @@ -207,7 +205,8 @@ public: y.compress(VectorOperation::insert); // Integrate the ODE. - time_stepper.solve(y); + auto nt = time_stepper.solve(y); + out << "# Number of steps taken: " << nt << std::endl; } private: diff --git a/tests/petsc/petsc_ts_00.output b/tests/petsc/petsc_ts_00.output index 3eba9270de..c6233771c3 100644 --- a/tests/petsc/petsc_ts_00.output +++ b/tests/petsc/petsc_ts_00.output @@ -116,6 +116,7 @@ end 5.875 -0.402097 (-0.396944) 0.914224 (0.917843) 5.9375 -0.344285 (-0.338842) 0.937523 (0.940843) 6 -0.285132 (-0.279415) 0.957169 (0.96017) +# Number of steps taken: 80 # Test implicit interface (J 0) 1 0.841471 (0.841471) 0.540302 (0.540302) 1.0625 0.873 (0.873575) 0.486385 (0.48669) @@ -198,6 +199,7 @@ end 5.875 -0.402097 (-0.396944) 0.914224 (0.917843) 5.9375 -0.344285 (-0.338842) 0.937523 (0.940843) 6 -0.285132 (-0.279415) 0.957169 (0.96017) +# Number of steps taken: 80 # Test explicit interface (J 1) 1 0.841471 (0.841471) 0.540302 (0.540302) 1.0625 0.873 (0.873575) 0.486385 (0.48669) @@ -280,6 +282,7 @@ end 5.875 -0.402097 (-0.396944) 0.914224 (0.917843) 5.9375 -0.344285 (-0.338842) 0.937523 (0.940843) 6 -0.285132 (-0.279415) 0.957169 (0.96017) +# Number of steps taken: 80 # Test implicit interface (J 1) 1 0.841471 (0.841471) 0.540302 (0.540302) 1.0625 0.873 (0.873575) 0.486385 (0.48669) @@ -362,6 +365,7 @@ end 5.875 -0.402097 (-0.396944) 0.914224 (0.917843) 5.9375 -0.344285 (-0.338842) 0.937523 (0.940843) 6 -0.285132 (-0.279415) 0.957169 (0.96017) +# Number of steps taken: 80 # Test user interface 1 0.841471 (0.841471) 0.540302 (0.540302) 1.0625 0.873 (0.873575) 0.486385 (0.48669) @@ -444,3 +448,4 @@ end 5.875 -0.402097 (-0.396944) 0.914224 (0.917843) 5.9375 -0.344285 (-0.338842) 0.937523 (0.940843) 6 -0.285132 (-0.279415) 0.957169 (0.96017) +# Number of steps taken: 80 diff --git a/tests/petsc/petsc_ts_01.cc b/tests/petsc/petsc_ts_01.cc index b2aee68d48..a20c22d47c 100644 --- a/tests/petsc/petsc_ts_01.cc +++ b/tests/petsc/petsc_ts_01.cc @@ -12,8 +12,6 @@ // the top level directory of deal.II. // //----------------------------------------------------------- -// -// Author: Stefano Zampini, King Abdullah University of Science and Technology. #include #include @@ -78,6 +76,7 @@ main(int argc, char **argv) auto ts = myode.petsc_ts(); auto t0 = myode.get_time(); auto dt = myode.get_time_step(); + AssertThrow(ts == static_cast(myode), ExcInternalError()); myode.solve(v, A); } catch (std::exception &exc) diff --git a/tests/petsc/petsc_ts_02.cc b/tests/petsc/petsc_ts_02.cc index 30ea46b1db..43315dc92f 100644 --- a/tests/petsc/petsc_ts_02.cc +++ b/tests/petsc/petsc_ts_02.cc @@ -12,8 +12,6 @@ // the top level directory of deal.II. // //----------------------------------------------------------- -// -// Author: Stefano Zampini, King Abdullah University of Science and Technology. //#include @@ -37,7 +35,7 @@ public: { // Customize solver: use BDF and basic adaptive time stepping PETScWrappers::TimeStepperData data; - data.tstype = "bdf"; + data.ts_type = "bdf"; data.final_time = 1.0; data.tsadapttype = "basic"; reinit(data);