From 83a437a5b8e36ed3f01d467b723b9d20d622d931 Mon Sep 17 00:00:00 2001 From: David Wells Date: Mon, 6 Jun 2022 14:21:56 -0400 Subject: [PATCH] Propagate SUNContext objects to linear solvers. --- include/deal.II/sundials/sunlinsol_wrapper.h | 51 +++++++- source/sundials/sunlinsol_wrapper.cc | 130 ++++++++++++++----- 2 files changed, 144 insertions(+), 37 deletions(-) diff --git a/include/deal.II/sundials/sunlinsol_wrapper.h b/include/deal.II/sundials/sunlinsol_wrapper.h index fb848dcd68..509596f629 100644 --- a/include/deal.II/sundials/sunlinsol_wrapper.h +++ b/include/deal.II/sundials/sunlinsol_wrapper.h @@ -51,6 +51,16 @@ namespace SUNDIALS void vmult(VectorType &dst, const VectorType &src) const; +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + /** + * Constructor. + * + * @param A_data Data required by @p a_times_fn + * @param a_times_fn A function pointer to the function that computes A*v + * @param linsol_ctx The context object used to set up the linear solver and all vectors + */ + SundialsOperator(void *A_data, ATimesFn a_times_fn, SUNContext linsol_ctx); +# else /** * Constructor. * @@ -58,6 +68,7 @@ namespace SUNDIALS * @param a_times_fn A function pointer to the function that computes A*v */ SundialsOperator(void *A_data, ATimesFn a_times_fn); +# endif private: /** @@ -70,6 +81,13 @@ namespace SUNDIALS * product. */ ATimesFn a_times_fn; + +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + /** + * Context object used for SUNDIALS logging. + */ + SUNContext linsol_ctx; +# endif }; @@ -92,6 +110,25 @@ namespace SUNDIALS void vmult(VectorType &dst, const VectorType &src) const; +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + /** + * Constructor. + * + * @param P_data Data required by @p p_solve_fn + * @param p_solve_fn A function pointer to the function that computes A*v + * @param linsol_ctx The context object used to set up the linear solver and all vectors + * @param tol Tolerance, that an iterative solver should use to judge + * convergence + * + * @note This function is only available with SUNDIALS 6.0.0 and later. + * 6.0.0. If you are using an earlier version of SUNDIALS then you need to + * use the other constructor. + */ + SundialsPreconditioner(void * P_data, + PSolveFn p_solve_fn, + SUNContext linsol_ctx, + double tol); +# else /** * Constructor. * @@ -99,8 +136,13 @@ namespace SUNDIALS * @param p_solve_fn A function pointer to the function that computes A*v * @param tol Tolerance, that an iterative solver should use to judge * convergence + * + * @note This function is only available with versions of SUNDIALS prior to + * 6.0.0. If you are using SUNDIALS 6 or later then you need to use the + * other constructor. */ SundialsPreconditioner(void *P_data, PSolveFn p_solve_fn, double tol); +# endif private: /** @@ -114,6 +156,13 @@ namespace SUNDIALS */ PSolveFn p_solve_fn; +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + /** + * Context object used for SUNDIALS logging. + */ + SUNContext linsol_ctx; +# endif + /** * Potential tolerance to use in the internal solve of the preconditioner * equation. @@ -169,7 +218,7 @@ namespace SUNDIALS explicit LinearSolverWrapper(LinearSolveFunction lsolve # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , - SUNContext &linsol_ctx + SUNContext linsol_ctx # endif ); diff --git a/source/sundials/sunlinsol_wrapper.cc b/source/sundials/sunlinsol_wrapper.cc index ae8897cc06..579ff6fd2f 100644 --- a/source/sundials/sunlinsol_wrapper.cc +++ b/source/sundials/sunlinsol_wrapper.cc @@ -60,10 +60,25 @@ namespace SUNDIALS template struct LinearSolverContent { + LinearSolverContent() + : a_times_fn(nullptr) + , preconditioner_setup(nullptr) + , preconditioner_solve(nullptr) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + , linsol_ctx(nullptr) +# endif + , P_data(nullptr) + , A_data(nullptr) + {} + ATimesFn a_times_fn; PSetupFn preconditioner_setup; PSolveFn preconditioner_solve; +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + SUNContext linsol_ctx; +# endif + LinearSolveFunction lsolve; void *P_data; @@ -109,10 +124,21 @@ namespace SUNDIALS auto *src_b = unwrap_nvector_const(b); auto *dst_x = unwrap_nvector(x); - SundialsOperator op(content->A_data, content->a_times_fn); + SundialsOperator op(content->A_data, + content->a_times_fn +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + , + content->linsol_ctx +# endif + ); SundialsPreconditioner preconditioner( - content->P_data, content->preconditioner_solve, tol); + content->P_data, + content->preconditioner_solve, +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + content->linsol_ctx, +# endif + tol); return content->lsolve(op, preconditioner, *dst_x, *src_b, tol); } @@ -173,17 +199,17 @@ namespace SUNDIALS template internal::LinearSolverWrapper::LinearSolverWrapper( LinearSolveFunction lsolve -# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , - SUNContext &linsol_ctx + SUNContext linsol_ctx # endif ) : content(std::make_unique>()) { -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) - sun_linear_solver = SUNLinSolNewEmpty(); -# else +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx); +# else + sun_linear_solver = SUNLinSolNewEmpty(); # endif sun_linear_solver->ops->gettype = arkode_linsol_get_type; @@ -194,7 +220,10 @@ namespace SUNDIALS sun_linear_solver->ops->setpreconditioner = arkode_linsol_set_preconditioner; - content->lsolve = lsolve; + content->lsolve = lsolve; +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + content->linsol_ctx = linsol_ctx; +# endif sun_linear_solver->content = content.get(); } @@ -219,8 +248,22 @@ namespace SUNDIALS +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + template + SundialsOperator::SundialsOperator(void * A_data, + ATimesFn a_times_fn, + SUNContext linsol_ctx) + : A_data(A_data) + , a_times_fn(a_times_fn) + , linsol_ctx(linsol_ctx) + + { + Assert(a_times_fn != nullptr, ExcInternalError()); + Assert(linsol_ctx != nullptr, ExcInternalError()); + } +# else template - SundialsOperator::SundialsOperator(void * A_data, + SundialsOperator::SundialsOperator(void *A_data, ATimesFn a_times_fn) : A_data(A_data) , a_times_fn(a_times_fn) @@ -228,6 +271,7 @@ namespace SUNDIALS { Assert(a_times_fn != nullptr, ExcInternalError()); } +# endif @@ -236,34 +280,48 @@ namespace SUNDIALS SundialsOperator::vmult(VectorType & dst, const VectorType &src) const { -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) - auto sun_dst = internal::make_nvector_view(dst); - auto sun_src = internal::make_nvector_view(src); - int status = a_times_fn(A_data, sun_src, sun_dst); + auto sun_dst = internal::make_nvector_view(dst +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + , + linsol_ctx +# endif + ); + auto sun_src = internal::make_nvector_view(src +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + , + linsol_ctx +# endif + ); + int status = a_times_fn(A_data, sun_src, sun_dst); (void)status; AssertSundialsSolver(status); -# else - // We don't currently know how to implement this: the - // internal::make_nvector_view() function called above requires a - // SUNContext object, but we don't actually have access to such an - // object here... - (void)dst; - (void)src; - Assert(false, ExcNotImplemented()); -# endif } +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + template + SundialsPreconditioner::SundialsPreconditioner( + void * P_data, + PSolveFn p_solve_fn, + SUNContext linsol_ctx, + double tol) + : P_data(P_data) + , p_solve_fn(p_solve_fn) + , linsol_ctx(linsol_ctx) + , tol(tol) + {} +# else template SundialsPreconditioner::SundialsPreconditioner( - void * P_data, + void *P_data, PSolveFn p_solve_fn, - double tol) + double tol) : P_data(P_data) , p_solve_fn(p_solve_fn) , tol(tol) {} +# endif @@ -279,24 +337,24 @@ namespace SUNDIALS return; } -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) - auto sun_dst = internal::make_nvector_view(dst); - auto sun_src = internal::make_nvector_view(src); + auto sun_dst = internal::make_nvector_view(dst +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + , + linsol_ctx +# endif + ); + auto sun_src = internal::make_nvector_view(src +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + , + linsol_ctx +# endif + ); // for custom preconditioners no distinction between left and right // preconditioning is made int status = p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/); (void)status; AssertSundialsSolver(status); -# else - // We don't currently know how to implement this: the - // internal::make_nvector_view() function called above requires a - // SUNContext object, but we don't actually have access to such an - // object here... - (void)dst; - (void)src; - Assert(false, ExcNotImplemented()); -# endif } template struct SundialsOperator>; -- 2.39.5