From cf96b10e73b13a2efbeebc09a3d2d8d33a3e86a5 Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Sun, 23 Jul 2023 12:07:37 +0200 Subject: [PATCH] PETScWrappers::TimeStepper support resizing while solving add hooks to completely support deal.II style improve documentation --- include/deal.II/lac/petsc_compatibility.h | 8 + include/deal.II/lac/petsc_ts.h | 125 +++++++- include/deal.II/lac/petsc_ts.templates.h | 345 +++++++++++++++++++--- source/lac/petsc_compatibility.cc | 7 + 4 files changed, 436 insertions(+), 49 deletions(-) diff --git a/include/deal.II/lac/petsc_compatibility.h b/include/deal.II/lac/petsc_compatibility.h index 449f68190f..9ed6b20533 100644 --- a/include/deal.II/lac/petsc_compatibility.h +++ b/include/deal.II/lac/petsc_compatibility.h @@ -175,6 +175,14 @@ namespace PETScWrappers + /** + * Reset DM (no public API). + */ + void + ts_reset_dm(TS ts); + + + /** * Set final time for ODE integration. */ diff --git a/include/deal.II/lac/petsc_ts.h b/include/deal.II/lac/petsc_ts.h index 92a088b309..72b6e818f7 100644 --- a/include/deal.II/lac/petsc_ts.h +++ b/include/deal.II/lac/petsc_ts.h @@ -263,18 +263,29 @@ namespace PETScWrappers * are ODE-solver specific. For details, consult the [PETSc * manual](https://petsc.org/release/manual/snes/#sec-nlmatrixfree). * - * In alternative, users can also provide the implementations of the - * *Jacobians*. This can be accomplished in two ways: + * Users can also provide the implementations of the *Jacobians*. This can be + * accomplished in two ways: * - PETSc style using TimeStepper::implicit_jacobian * and TimeStepper::explicit_jacobian. * - deal.II style using TimeStepper::setup_jacobian and - * TimeStepper::solve_with_jacobian + * TimeStepper::solve_with_jacobian. + * The preconditioning matrix can be specified using TimeStepper::set_matrix. + * In case both approaches are implemented, the deal.II style will be used. * - * In case both approaches are coded, the deal.II style - * will be used. + * TimeStepper::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 second approach is more in style with the deal.II philosophy - * and it also allows command line customization, like for example, + * The correctness of the constructed Jacobians passed using + * TimeStepper::set_matrix can be checked using + * @code + * ./myApp -snes_test_jacobian + * @endcode + * See TimeStepper::set_matrix and TimeStepper::set_matrices for + * additional details. + * + * The deal.II style approach still allows command line customization, like + * for example, * @code * ./myApp -snes_type newtontr -ksp_type cg * @endcode @@ -282,15 +293,14 @@ namespace PETScWrappers * a trust region solver and iterate on the tangent system with CG, * 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 - * faster turnaround for experimenting with linear solver preconditioning - * configurations via command line customizations, like for example, + * The PETSc style approach has instead the advantage that only the matrix + * assembly procedure has to be implemented, 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 - * See TimeStepper::set_matrix and TimeStepper::set_matrices for - * additional details. * * @dealiiConceptRequires{(concepts::is_dealii_petsc_vector_type * || std::constructible_from) && @@ -386,7 +396,8 @@ namespace PETScWrappers /** * Set both the linear system matrix and the preconditioning matrix - * that PETSc will use. + * that PETSc will use (can be the same matrix). In this case, the + * Jacobian-Free-Newton-Krylov approach will not be used. */ void set_matrices(AMatrixType &A, PMatrixType &P); @@ -403,6 +414,12 @@ namespace PETScWrappers real_type get_time_step(); + /** + * Return current step number. + */ + unsigned int + get_step_number(); + /** * Integrate the differential-algebraic equations starting from @p y. * @@ -438,6 +455,26 @@ namespace PETScWrappers unsigned int solve(VectorType &y, AMatrixType &A, PMatrixType &P); + /** + * Setup callbacks. + * + * This function is called inside TimeStepper::solve routines and does + * not need to be called by the user. It is used to reinitialize the + * solver if mesh adaption has been performed. + */ + void + setup_callbacks(); + + /** + * Setup algebraic constraints. + * + * This function is called inside TimeStepper::solve routines and does + * not need to be called by the user. It is used to reinitialize the + * solver if mesh adaption has been performed. + */ + void + setup_algebraic_constraints(const VectorType &y); + /** * Callback for the computation of the implicit residual $F(t, y, \dot y)$. * @@ -565,6 +602,61 @@ namespace PETScWrappers */ std::function algebraic_components; + /** + * Callback to distribute solution to hanging nodes. + * + * Implementation of this function is optional. + * It is called at the end of each successful stage. + * The same functionality can be equivalently implemented in + * TimeStepper::solve_with_jacobian. + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. + */ + std::function distribute; + + /** + * Callback to setup mesh adaption. + * + * Implementation of this function is optional. + * @p resize must be set to true if mesh adaption is to be performed, false + * otherwise. + * The @p y vector contains the current solution, @p t the current time, + * @ step the step number. + * Solution transfer and actual mesh adaption must be performed in a + * separate callback, TimeStepper::interpolate + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. + */ + std::function + prepare_for_coarsening_and_refinement; + + /** + * Callback to interpolate vectors and perform mesh adaption. + * + * Implementation of this function is mandatory if + * TimeStepper::prepare_for_coarsening_and_refinement is used. + * This function must perform mesh adaption and interpolate the discrete + * functions that are stored in @p all_in onto the refined and/or coarsenend grid. + * Output vectors must be created inside the callback. + * + * @note This variable represents a + * @ref GlossUserProvidedCallBack "user provided callback". + * See there for a description of how to deal with errors and other + * requirements and conventions. + */ + std::function &all_in, + std::vector & all_out)> + interpolate; + private: /** * The PETSc object. @@ -577,6 +669,11 @@ namespace PETScWrappers SmartPointer A; SmartPointer P; + /** + * Object to apply solve_with_jacobian. + */ + PreconditionShell solve_with_jacobian_pc; + /** * This flag is set when changing the customization and used within solve. */ diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h index 0d72db1db0..a403edb405 100644 --- a/include/deal.II/lac/petsc_ts.templates.h +++ b/include/deal.II/lac/petsc_ts.templates.h @@ -130,7 +130,8 @@ namespace PETScWrappers TimeStepper::TimeStepper( const TimeStepperData &data, const MPI_Comm mpi_comm) - : pending_exception(nullptr) + : solve_with_jacobian_pc(mpi_comm) + , pending_exception(nullptr) { AssertPETSc(TSCreate(mpi_comm, &ts)); AssertPETSc(TSSetApplicationContext(ts, this)); @@ -151,6 +152,12 @@ namespace PETScWrappers std::constructible_from)) TimeStepper::~TimeStepper() { + AssertPETSc(PetscObjectComposeFunction((PetscObject)ts, + "__dealii_ts_resize_setup__", + nullptr)); + AssertPETSc(PetscObjectComposeFunction((PetscObject)ts, + "__dealii_ts_resize_transfer__", + nullptr)); AssertPETSc(TSDestroy(&ts)); Assert(pending_exception == nullptr, ExcInternalError()); @@ -231,6 +238,24 @@ namespace PETScWrappers + template + DEAL_II_CXX20_REQUIRES( + (concepts::is_dealii_petsc_vector_type || + std::constructible_from< + VectorType, + Vec>)&&(concepts::is_dealii_petsc_matrix_type || + std::constructible_from< + PMatrixType, + Mat>)&&(concepts::is_dealii_petsc_matrix_type || + std::constructible_from)) + unsigned int TimeStepper:: + get_step_number() + { + return ts_get_step_number(ts); + } + + + template DEAL_II_CXX20_REQUIRES( (concepts::is_dealii_petsc_vector_type || @@ -401,12 +426,149 @@ namespace PETScWrappers PMatrixType, Mat>)&&(concepts::is_dealii_petsc_matrix_type || std::constructible_from)) - unsigned int TimeStepper::solve( - VectorType &y) + void TimeStepper::setup_callbacks() { - const auto ts_prestage = [](TS ts, PetscReal) -> PetscErrorCode { + const auto ts_resize_setup = + [](TS ts, PetscInt it, PetscReal t, Vec x, PetscBool *flg, void *ctx) + -> PetscErrorCode { + PetscFunctionBeginUser; + auto user = static_cast(ctx); + + VectorType xdealii(x); + bool resize = false; + + const int lineno = __LINE__; + const int err = call_and_possibly_capture_ts_exception( + user->prepare_for_coarsening_and_refinement, + user->pending_exception, + []() -> void {}, + t, + it, + xdealii, + resize); + *flg = resize ? PETSC_TRUE : PETSC_FALSE; + if (err) + return PetscError( + PetscObjectComm((PetscObject)ts), + lineno + 1, + "prepare_for_coarsening_and_refinement", + __FILE__, + PETSC_ERR_LIB, + PETSC_ERROR_INITIAL, + "Failure in ts_resize_setup from dealii::PETScWrappers::TimeStepper"); + + PetscFunctionReturn(PETSC_SUCCESS); + }; + + const auto ts_resize_transfer = [](TS ts, + PetscInt nv, + Vec vin[], + Vec vout[], + void * ctx) -> PetscErrorCode { + PetscFunctionBeginUser; + auto user = static_cast(ctx); + + std::vector all_in, all_out; + for (PetscInt i = 0; i < nv; i++) + { + all_in.push_back(VectorType(vin[i])); + } + + const int lineno = __LINE__; + const int err = call_and_possibly_capture_ts_exception( + user->interpolate, + user->pending_exception, + []() -> void {}, + all_in, + all_out); + if (err) + return PetscError( + PetscObjectComm((PetscObject)ts), + lineno + 1, + "interpolate", + __FILE__, + PETSC_ERR_LIB, + PETSC_ERROR_INITIAL, + "Failure in ts_resize_transfer from dealii::PETScWrappers::TimeStepper"); + + if (all_out.size() != all_in.size()) + SETERRQ(PetscObjectComm((PetscObject)ts), + PETSC_ERR_LIB, + "Broken all_out"); + + for (PetscInt i = 0; i < nv; i++) + { + vout[i] = all_out[i].petsc_vector(); + PetscCall(PetscObjectReference((PetscObject)vout[i])); + } + + try + { + user->setup_callbacks(); + } + catch (...) + { + SETERRQ(PetscObjectComm((PetscObject)ts), + PETSC_ERR_LIB, + "Exception in setup_callbacks"); + } + + try + { + user->setup_algebraic_constraints(all_out[0]); + } + catch (...) + { + SETERRQ(PetscObjectComm((PetscObject)ts), + PETSC_ERR_LIB, + "Exception in setup_algebraic_constraints"); + } + PetscFunctionReturn(PETSC_SUCCESS); + }; + + const auto ts_poststep_amr = [](TS ts) -> PetscErrorCode { + PetscFunctionBeginUser; void *ctx; + PetscCall(TSGetApplicationContext(ts, &ctx)); + auto user = static_cast(ctx); + + PetscErrorCode ( + *transfer_setup)(TS, PetscInt, PetscReal, Vec, PetscBool *, void *); + PetscErrorCode (*transfer)(TS, PetscInt, Vec[], Vec[], void *); + AssertPETSc(PetscObjectQueryFunction((PetscObject)ts, + "__dealii_ts_resize_setup__", + &transfer_setup)); + AssertPETSc(PetscObjectQueryFunction((PetscObject)ts, + "__dealii_ts_resize_transfer__", + &transfer)); + + PetscBool resize = PETSC_FALSE; + Vec x; + PetscCall(TSGetSolution(ts, &x)); + PetscCall(PetscObjectReference((PetscObject)x)); + PetscCall((*transfer_setup)( + ts, user->get_step_number(), user->get_time(), x, &resize, user)); + + if (resize) + { + Vec new_x; + + user->reinit(); + + // Reinitialize DM. Currently it is not possible to do with public API + ts_reset_dm(ts); + + PetscCall((*transfer)(ts, 1, &x, &new_x, user)); + PetscCall(TSSetSolution(ts, new_x)); + PetscCall(VecDestroy(&new_x)); + } + PetscCall(VecDestroy(&x)); + PetscFunctionReturn(PETSC_SUCCESS); + }; + + const auto ts_prestage = [](TS ts, PetscReal) -> PetscErrorCode { PetscFunctionBeginUser; + void *ctx; PetscCall(TSGetApplicationContext(ts, &ctx)); auto user = static_cast(ctx); user->error_in_function = false; @@ -418,17 +580,52 @@ namespace PETScWrappers PetscCall(TSGetSNES(ts, &snes)); snes_reset_domain_flags(snes); } - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); + }; + + const auto ts_poststage = + [](TS ts, PetscReal t, PetscInt i, Vec *xx) -> PetscErrorCode { + PetscFunctionBeginUser; + void *ctx; + PetscCall(TSGetApplicationContext(ts, &ctx)); + auto user = static_cast(ctx); + + VectorType xdealii(xx[i]); + + const int lineno = __LINE__; + const int err = call_and_possibly_capture_ts_exception( + user->distribute, + user->pending_exception, + [user, ts]() -> void { + user->error_in_function = true; + + SNES snes; + AssertPETSc(TSGetSNES(ts, &snes)); + AssertPETSc(SNESSetFunctionDomainError(snes)); + }, + t, + xdealii); + if (err) + return PetscError( + PetscObjectComm((PetscObject)ts), + lineno + 1, + "distribute", + __FILE__, + PETSC_ERR_LIB, + PETSC_ERROR_INITIAL, + "Failure in ts_poststage from dealii::PETScWrappers::TimeStepper"); + petsc_increment_state_counter(xx[i]); + PetscFunctionReturn(PETSC_SUCCESS); }; const auto ts_functiondomainerror = [](TS ts, PetscReal, Vec, PetscBool *accept) -> PetscErrorCode { - void *ctx; PetscFunctionBeginUser; + void *ctx; AssertPETSc(TSGetApplicationContext(ts, &ctx)); auto user = static_cast(ctx); *accept = user->error_in_function ? PETSC_FALSE : PETSC_TRUE; - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); }; const auto ts_ifunction = @@ -466,7 +663,7 @@ namespace PETScWrappers PETSC_ERROR_INITIAL, "Failure in ts_ifunction from dealii::PETScWrappers::TimeStepper"); petsc_increment_state_counter(f); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); }; const auto ts_ijacobian = [](TS ts, @@ -524,7 +721,7 @@ namespace PETScWrappers else petsc_increment_state_counter(A); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); }; const auto ts_ijacobian_with_setup = [](TS ts, @@ -598,7 +795,7 @@ namespace PETScWrappers else petsc_increment_state_counter(A); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); }; const auto ts_rhsfunction = @@ -636,7 +833,7 @@ namespace PETScWrappers PETSC_ERROR_INITIAL, "Failure in ts_rhsfunction from dealii::PETScWrappers::TimeStepper"); petsc_increment_state_counter(f); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); }; const auto ts_rhsjacobian = @@ -685,7 +882,7 @@ namespace PETScWrappers else petsc_increment_state_counter(A); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); }; const auto ts_monitor = @@ -712,27 +909,25 @@ namespace PETScWrappers PETSC_ERR_LIB, PETSC_ERROR_INITIAL, "Failure in ts_monitor from dealii::PETScWrappers::TimeStepper"); - PetscFunctionReturn(0); + PetscFunctionReturn(PETSC_SUCCESS); }; AssertThrow(explicit_function || implicit_function, StandardExceptions::ExcFunctionNotProvided( "explicit_function || implicit_function")); + // Handle recoverable errors. this->error_in_function = false; - - AssertPETSc(TSSetSolution(ts, y.petsc_vector())); - - // Handle recoverable errors AssertPETSc(TSSetPreStage(ts, ts_prestage)); AssertPETSc(TSSetFunctionDomainError(ts, ts_functiondomainerror)); + // Problem description. if (explicit_function) AssertPETSc(TSSetRHSFunction(ts, nullptr, ts_rhsfunction, this)); - if (implicit_function) AssertPETSc(TSSetIFunction(ts, nullptr, ts_ifunction, this)); + // Handle Jacobians. this->need_dummy_assemble = false; if (setup_jacobian) { @@ -801,12 +996,10 @@ namespace PETScWrappers // 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_with_jacobian) { - precond.vmult = [&](VectorBase & indst, - const VectorBase &insrc) -> void { + solve_with_jacobian_pc.vmult = [&](VectorBase & indst, + const VectorBase &insrc) -> void { VectorType dst(static_cast(indst)); const VectorType src(static_cast(insrc)); solve_with_jacobian(src, dst); @@ -818,25 +1011,78 @@ namespace PETScWrappers AssertPETSc(TSGetSNES(ts, &snes)); AssertPETSc(SNESGetKSP(snes, &ksp)); AssertPETSc(KSPSetType(ksp, KSPPREONLY)); - AssertPETSc(KSPSetPC(ksp, precond.get_pc())); + AssertPETSc(KSPSetPC(ksp, solve_with_jacobian_pc.get_pc())); } + if (distribute) + AssertPETSc(TSSetPostStage(ts, ts_poststage)); + // Attach user monitoring routine. if (monitor) AssertPETSc(TSMonitorSet(ts, ts_monitor, this, nullptr)); + // Handle AMR. + if (prepare_for_coarsening_and_refinement) + { + AssertThrow(interpolate, + StandardExceptions::ExcFunctionNotProvided("interpolate")); +# if DEAL_II_PETSC_VERSION_GTE(3, 20, 0) + (void)ts_poststep_amr; + AssertPETSc(TSSetResize(ts, ts_resize_setup, ts_resize_transfer, this)); +# else + AssertPETSc(PetscObjectComposeFunction( + (PetscObject)ts, + "__dealii_ts_resize_setup__", + static_cast( + ts_resize_setup))); + AssertPETSc(PetscObjectComposeFunction( + (PetscObject)ts, + "__dealii_ts_resize_transfer__", + static_cast( + ts_resize_transfer))); + AssertPETSc(TSSetPostStep(ts, ts_poststep_amr)); +# endif + } + // Allow command line customization. AssertPETSc(TSSetFromOptions(ts)); - // Handle algebraic components. + // By default PETSc does not check for Jacobian errors in optimized + // mode. Here we do it unconditionally. +# if DEAL_II_PETSC_VERSION_GTE(3, 11, 0) + { + SNES snes; + AssertPETSc(TSGetSNES(ts, &snes)); + AssertPETSc(SNESSetCheckJacobianDomainError(snes, PETSC_TRUE)); + } +# endif + } + + + + template + DEAL_II_CXX20_REQUIRES( + (concepts::is_dealii_petsc_vector_type || + std::constructible_from< + VectorType, + Vec>)&&(concepts::is_dealii_petsc_matrix_type || + std::constructible_from< + PMatrixType, + Mat>)&&(concepts::is_dealii_petsc_matrix_type || + std::constructible_from)) + void TimeStepper:: + setup_algebraic_constraints(const VectorType &y) + { if (algebraic_components && need_dae_tolerances) { PetscReal atol, rtol; AssertPETSc(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr)); Vec av, rv; - AssertPETSc(VecDuplicate(y.petsc_vector(), &av)); - AssertPETSc(VecDuplicate(y.petsc_vector(), &rv)); + Vec py = const_cast(y).petsc_vector(); + AssertPETSc(VecDuplicate(py, &av)); + AssertPETSc(VecDuplicate(py, &rv)); VectorBase avdealii(av); VectorBase rvdealii(rv); @@ -853,21 +1099,50 @@ namespace PETScWrappers AssertPETSc(VecDestroy(&av)); AssertPETSc(VecDestroy(&rv)); } + } + + + + template + DEAL_II_CXX20_REQUIRES( + (concepts::is_dealii_petsc_vector_type || + std::constructible_from< + VectorType, + Vec>)&&(concepts::is_dealii_petsc_matrix_type || + std::constructible_from< + PMatrixType, + Mat>)&&(concepts::is_dealii_petsc_matrix_type || + std::constructible_from)) + unsigned int TimeStepper::solve( + VectorType &y) + { + try + { + setup_callbacks(); + } + catch (...) + { + throw; + } + + // Set solution vector. + AssertPETSc(TSSetSolution(ts, y.petsc_vector())); // Setup internal data structures, including the internal nonlinear // solver SNES if using an implicit solver. AssertPETSc(TSSetUp(ts)); - // By default PETSc does not check for Jacobian errors in optimized - // mode. Here we do it unconditionally. -# if DEAL_II_PETSC_VERSION_GTE(3, 11, 0) - if (ts_has_snes(ts)) + // Handle algebraic components. This must be done when we know the layout + // of the solution vector. + try { - SNES snes; - AssertPETSc(TSGetSNES(ts, &snes)); - AssertPETSc(SNESSetCheckJacobianDomainError(snes, PETSC_TRUE)); + setup_algebraic_constraints(y); } -# endif + catch (...) + { + throw; + } + // Having set everything up, now do the actual work // and let PETSc do the time stepping. If there is // a pending exception, then one of the user callbacks @@ -877,7 +1152,7 @@ namespace PETScWrappers // a zero error code -- if so, just eat the exception and // continue on; otherwise, just rethrow the exception // and get outta here. - const int status = TSSolve(ts, nullptr); + const PetscErrorCode status = TSSolve(ts, nullptr); if (pending_exception) { try diff --git a/source/lac/petsc_compatibility.cc b/source/lac/petsc_compatibility.cc index b787052001..7a6a02893f 100644 --- a/source/lac/petsc_compatibility.cc +++ b/source/lac/petsc_compatibility.cc @@ -26,6 +26,7 @@ # include # include # include +# include // Shorthand notation for PETSc error codes. # define AssertPETSc(code) \ @@ -128,6 +129,12 @@ namespace PETScWrappers # endif } + void + ts_reset_dm(TS ts) + { + AssertPETSc(DMDestroy(&ts->dm)); + } + unsigned int ts_get_step_number(TS ts) { -- 2.39.5