From b124fb2887c1643b8673759e0d0b4a71479b1ba2 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 22 Feb 2022 20:59:07 -0700 Subject: [PATCH] Support SUNDIALS 6.0 and later. --- include/deal.II/sundials/arkode.h | 9 +- include/deal.II/sundials/ida.h | 9 +- include/deal.II/sundials/kinsol.h | 12 ++ include/deal.II/sundials/n_vector.h | 18 ++- include/deal.II/sundials/n_vector.templates.h | 66 +++++++-- include/deal.II/sundials/sunlinsol_wrapper.h | 9 +- source/sundials/arkode.cc | 136 +++++++++++++++--- source/sundials/ida.cc | 97 ++++++++++--- source/sundials/kinsol.cc | 86 +++++++++-- source/sundials/n_vector.inst.in | 56 ++++++-- source/sundials/sunlinsol_wrapper.cc | 34 ++++- 11 files changed, 463 insertions(+), 69 deletions(-) diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index ad2ca7d863..eb31f2e22f 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -1267,12 +1267,19 @@ namespace SUNDIALS */ void *arkode_mem; +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + /** + * A context object associated with the ARKode solver. + */ + SUNContext arkode_ctx; +# endif + /** * MPI communicator. SUNDIALS solver runs happily in * parallel. Note that if the library is compiled without MPI * support, MPI_Comm is aliased as int. */ - MPI_Comm communicator; + MPI_Comm mpi_communicator; /** * The final time in the last call to solve_ode(). diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 16f6f5070c..bc50c9b7e4 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -882,12 +882,19 @@ namespace SUNDIALS */ void *ida_mem; +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + /** + * A context object associated with the IDA solver. + */ + SUNContext ida_ctx; +# endif + /** * MPI communicator. SUNDIALS solver runs happily in * parallel. Note that if the library is compiled without MPI * support, MPI_Comm is aliased as int. */ - MPI_Comm communicator; + MPI_Comm mpi_communicator; /** * Memory pool of vectors. diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 89c5a146c5..df2f6b16c7 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -684,11 +684,23 @@ namespace SUNDIALS */ AdditionalData data; + /** + * The MPI communicator to be used by this solver, if any. + */ + MPI_Comm mpi_communicator; + /** * KINSOL memory object. */ void *kinsol_mem; +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + /** + * A context object associated with the KINSOL solver. + */ + SUNContext kinsol_ctx; +# endif + /** * Memory pool of vectors. */ diff --git a/include/deal.II/sundials/n_vector.h b/include/deal.II/sundials/n_vector.h index df2c27f129..6f959444e6 100644 --- a/include/deal.II/sundials/n_vector.h +++ b/include/deal.II/sundials/n_vector.h @@ -56,13 +56,22 @@ namespace SUNDIALS * @tparam VectorType Type of the viewed vector. This parameter can be * deduced automatically and will respect a potential const-qualifier. * @param vector The vector to view. +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + * @param nvector_context A SUNDIALS context object to be used for the + * operations we want to do on this vector. +# endif * @return A NVectorView of the @p vector. * * @related NVectorView */ template NVectorView - make_nvector_view(VectorType &vector); + make_nvector_view(VectorType &vector +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext nvector_context +# endif + ); /** * Retrieve the underlying vector attached to N_Vector @p v. This call will @@ -129,7 +138,12 @@ namespace SUNDIALS /** * Constructor. Create view of @p vector. */ - NVectorView(VectorType &vector); + NVectorView(VectorType &vector +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext nvector_context +# endif + ); /** * Move assignment. diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index 002eb5eb66..ee5220d023 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -125,7 +125,12 @@ namespace SUNDIALS */ template N_Vector - create_nvector(NVectorContent *content); + create_nvector(NVectorContent *content +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext nvector_context +# endif + ); /** * Helper to create an empty vector with all operations but no content. @@ -133,7 +138,11 @@ namespace SUNDIALS */ template N_Vector - create_empty_nvector(); + create_empty_nvector( +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + SUNContext nvector_context +# endif + ); /** * Collection of all operations specified by SUNDIALS N_Vector @@ -356,9 +365,19 @@ namespace SUNDIALS template NVectorView - make_nvector_view(VectorType &vector) + make_nvector_view(VectorType &vector +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext nvector_context +# endif + ) { - return NVectorView(vector); + return NVectorView(vector +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + nvector_context +# endif + ); } @@ -390,11 +409,21 @@ namespace SUNDIALS template - NVectorView::NVectorView(VectorType &vector) + NVectorView::NVectorView(VectorType &vector +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext nvector_context +# endif + ) : vector_ptr( create_nvector( new NVectorContent::type>( - &vector)), + &vector) +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + nvector_context +# endif + ), [](N_Vector v) { N_VDestroy(v); }) {} @@ -421,10 +450,19 @@ namespace SUNDIALS template N_Vector - create_nvector(NVectorContent *content) + create_nvector(NVectorContent *content +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext nvector_context +# endif + ) { // Create an N_Vector with operators attached and empty content +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) N_Vector v = create_empty_nvector(); +# else + N_Vector v = create_empty_nvector(nvector_context); +# endif Assert(v != nullptr, ExcInternalError()); v->content = content; @@ -447,7 +485,11 @@ namespace SUNDIALS clone_empty(N_Vector w) { Assert(w != nullptr, ExcInternalError()); +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) N_Vector v = N_VNewEmpty(); +# else + N_Vector v = N_VNewEmpty(w->sunctx); +# endif Assert(v != nullptr, ExcInternalError()); int status = N_VCopyOps(w, v); @@ -940,9 +982,17 @@ namespace SUNDIALS template N_Vector - create_empty_nvector() + create_empty_nvector( +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + SUNContext nvector_context +# endif + ) { +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) N_Vector v = N_VNewEmpty(); +# else + N_Vector v = N_VNewEmpty(nvector_context); +# endif Assert(v != nullptr, ExcInternalError()); /* constructors, destructors, and utility operations */ diff --git a/include/deal.II/sundials/sunlinsol_wrapper.h b/include/deal.II/sundials/sunlinsol_wrapper.h index 4b94e1dc9b..7ae539e82c 100644 --- a/include/deal.II/sundials/sunlinsol_wrapper.h +++ b/include/deal.II/sundials/sunlinsol_wrapper.h @@ -156,7 +156,7 @@ namespace SUNDIALS namespace internal { - /*! + /** * Attach wrapper functions to SUNDIALS' linear solver interface. We pretend * that the user-supplied linear solver is matrix-free, even though it can * be matrix-based. This way SUNDIALS does not need to understand our matrix @@ -166,7 +166,12 @@ namespace SUNDIALS class LinearSolverWrapper { public: - explicit LinearSolverWrapper(LinearSolveFunction lsolve); + explicit LinearSolverWrapper(LinearSolveFunction lsolve +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext &linsol_ctx +# endif + ); ~LinearSolverWrapper(); diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index c68e57abd7..83cd451abf 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -368,9 +368,9 @@ namespace SUNDIALS const MPI_Comm & mpi_comm) : data(data) , arkode_mem(nullptr) - , communicator(is_serial_vector::value ? - MPI_COMM_SELF : - Utilities::MPI::duplicate_communicator(mpi_comm)) + , mpi_communicator(is_serial_vector::value ? + MPI_COMM_SELF : + Utilities::MPI::duplicate_communicator(mpi_comm)) , last_end_time(data.initial_time) { set_functions_to_trigger_an_assert(); @@ -382,15 +382,23 @@ namespace SUNDIALS ARKode::~ARKode() { if (arkode_mem) + { # if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - ARKodeFree(&arkode_mem); + ARKodeFree(&arkode_mem); # else - ARKStepFree(&arkode_mem); + ARKStepFree(&arkode_mem); +# endif + +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + const int status = SUNContext_Free(&arkode_ctx); + (void)status; + AssertARKode(status); # endif + } # ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) - Utilities::MPI::free_communicator(communicator); + Utilities::MPI::free_communicator(mpi_communicator); # endif } @@ -443,7 +451,12 @@ namespace SUNDIALS time.get_step_number()); } - auto solution_nvector = internal::make_nvector_view(solution); + auto solution_nvector = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); while (!time.is_at_end()) { @@ -510,7 +523,12 @@ namespace SUNDIALS // just a view on the memory in solution, all write operations on yy by // ARKODE will automatically be mirrored to solution - auto initial_condition_nvector = internal::make_nvector_view(solution); + auto initial_condition_nvector = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); status = ARKodeInit( arkode_mem, @@ -522,8 +540,12 @@ namespace SUNDIALS if (get_local_tolerances) { - const auto abs_tols = - internal::make_nvector_view(get_local_tolerances()); + const auto abs_tols = internal::make_nvector_view(get_local_tolerances() +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); status = ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols); AssertARKode(status); @@ -605,30 +627,58 @@ namespace SUNDIALS { last_end_time = current_time; if (arkode_mem) - ARKStepFree(&arkode_mem); + { + ARKStepFree(&arkode_mem); + +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + const int status = SUNContext_Free(&arkode_ctx); + (void)status; + AssertARKode(status); +# endif + } int status; (void)status; // just a view on the memory in solution, all write operations on yy by // ARKODE will automatically be mirrored to solution - auto initial_condition_nvector = internal::make_nvector_view(solution); + auto initial_condition_nvector = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); Assert(explicit_function || implicit_function, ExcFunctionNotProvided("explicit_function || implicit_function")); +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) arkode_mem = ARKStepCreate( explicit_function ? &explicit_function_callback : nullptr, implicit_function ? &implicit_function_callback : nullptr, current_time, initial_condition_nvector); +# else + status = SUNContext_Create(&mpi_communicator, &arkode_ctx); + AssertARKode(status); + arkode_mem = ARKStepCreate( + explicit_function ? &explicit_function_callback : nullptr, + implicit_function ? &implicit_function_callback : nullptr, + current_time, + initial_condition_nvector, + arkode_ctx); +# endif Assert(arkode_mem != nullptr, ExcInternalError()); if (get_local_tolerances) { - const auto abs_tols = - internal::make_nvector_view(get_local_tolerances()); + const auto abs_tols = internal::make_nvector_view(get_local_tolerances() +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); status = ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols); AssertARKode(status); @@ -683,18 +733,36 @@ namespace SUNDIALS { linear_solver = std::make_unique>( - solve_linearized_system); + solve_linearized_system +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); sun_linear_solver = *linear_solver; } else { // use default solver from SUNDIALS // TODO give user options - auto y_template = internal::make_nvector_view(solution); + auto y_template = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) sun_linear_solver = SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/); +# else + sun_linear_solver = + SUNLinSol_SPGMR(y_template, + PREC_NONE, + 0 /*krylov subvectors, 0 uses default*/, + arkode_ctx); +# endif } status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr); AssertARKode(status); @@ -724,11 +792,23 @@ namespace SUNDIALS } else { - auto y_template = internal::make_nvector_view(solution); + auto y_template = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) SUNNonlinearSolver fixed_point_solver = SUNNonlinSol_FixedPoint(y_template, data.anderson_acceleration_subspace); +# else + SUNNonlinearSolver fixed_point_solver = + SUNNonlinSol_FixedPoint(y_template, + data.anderson_acceleration_subspace, + arkode_ctx); +# endif status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver); AssertARKode(status); @@ -755,16 +835,34 @@ namespace SUNDIALS { mass_solver = std::make_unique>( - solve_mass); + solve_mass +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); sun_mass_linear_solver = *mass_solver; } else { - auto y_template = internal::make_nvector_view(solution); + auto y_template = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + arkode_ctx +# endif + ); +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) sun_mass_linear_solver = SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/); +# else + sun_mass_linear_solver = + SUNLinSol_SPGMR(y_template, + PREC_NONE, + 0 /*krylov subvectors, 0 uses default*/, + arkode_ctx); +# endif } booleantype mass_time_dependent = data.mass_is_time_independent ? SUNFALSE : SUNTRUE; diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index fe212f4e3e..8c886bca8d 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -197,21 +197,32 @@ namespace SUNDIALS IDA::IDA(const AdditionalData &data, const MPI_Comm &mpi_comm) : data(data) , ida_mem(nullptr) - , communicator(is_serial_vector::value ? - MPI_COMM_SELF : - Utilities::MPI::duplicate_communicator(mpi_comm)) + , mpi_communicator(is_serial_vector::value ? + MPI_COMM_SELF : + Utilities::MPI::duplicate_communicator(mpi_comm)) { set_functions_to_trigger_an_assert(); } + + template IDA::~IDA() { if (ida_mem) - IDAFree(&ida_mem); + { + IDAFree(&ida_mem); + +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + const int status = SUNContext_Free(&ida_ctx); + (void)status; + AssertIDA(status); +# endif + } + # ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) - Utilities::MPI::free_communicator(communicator); + Utilities::MPI::free_communicator(mpi_communicator); # endif } @@ -242,8 +253,18 @@ namespace SUNDIALS { next_time += data.output_period; - auto yy = internal::make_nvector_view(solution); - auto yp = internal::make_nvector_view(solution_dot); + auto yy = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + ida_ctx +# endif + ); + auto yp = internal::make_nvector_view(solution_dot +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + ida_ctx +# endif + ); status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL); AssertIDA(status); @@ -274,23 +295,53 @@ namespace SUNDIALS bool first_step = (current_time == data.initial_time); if (ida_mem) - IDAFree(&ida_mem); + { + IDAFree(&ida_mem); - ida_mem = IDACreate(); +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + const int status = SUNContext_Free(&ida_ctx); + (void)status; + AssertIDA(status); +# endif + } int status; (void)status; - auto yy = internal::make_nvector_view(solution); - auto yp = internal::make_nvector_view(solution_dot); + +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + ida_mem = IDACreate(); +# else + status = SUNContext_Create(&mpi_communicator, &ida_ctx); + AssertIDA(status); + + ida_mem = IDACreate(ida_ctx); +# endif + + auto yy = internal::make_nvector_view(solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + ida_ctx +# endif + ); + auto yp = internal::make_nvector_view(solution_dot +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + ida_ctx +# endif + ); status = IDAInit(ida_mem, residual_callback, current_time, yy, yp); AssertIDA(status); if (get_local_tolerances) { - const auto abs_tols = - internal::make_nvector_view(get_local_tolerances()); + const auto abs_tols = internal::make_nvector_view(get_local_tolerances() +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + ida_ctx +# endif + ); status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols); AssertIDA(status); } @@ -319,8 +370,13 @@ namespace SUNDIALS diff_comp_vector[*i] = 1.0; diff_comp_vector.compress(VectorOperation::insert); - const auto diff_id = internal::make_nvector_view(diff_comp_vector); - status = IDASetId(ida_mem, diff_id); + const auto diff_id = internal::make_nvector_view(diff_comp_vector +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + ida_ctx +# endif + ); + status = IDASetId(ida_mem, diff_id); AssertIDA(status); } @@ -352,7 +408,12 @@ namespace SUNDIALS // called do not actually receive the IDAMEM object, just the LS // object, so we have to store a pointer to the current // object in the LS object - LS = SUNLinSolNewEmpty(); +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + LS = SUNLinSolNewEmpty(); +# else + LS = SUNLinSolNewEmpty(ida_ctx); +# endif + LS->content = this; LS->ops->gettype = [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type { @@ -398,7 +459,11 @@ namespace SUNDIALS // if we don't set it, it won't call the functions that set up // the matrix object (i.e., the argument to the 'IDASetJacFn' // function below). +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) J = SUNMatNewEmpty(); +# else + J = SUNMatNewEmpty(ida_ctx); +# endif J->content = this; J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID { diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index d26f1bfe12..36cbc51c4b 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -303,8 +303,12 @@ namespace SUNDIALS template - KINSOL::KINSOL(const AdditionalData &data, const MPI_Comm &) + KINSOL::KINSOL(const AdditionalData &data, + const MPI_Comm & mpi_comm) : data(data) + , mpi_communicator(is_serial_vector::value ? + MPI_COMM_SELF : + Utilities::MPI::duplicate_communicator(mpi_comm)) , kinsol_mem(nullptr) { set_functions_to_trigger_an_assert(); @@ -315,8 +319,23 @@ namespace SUNDIALS template KINSOL::~KINSOL() { + int status = 0; + (void)status; + if (kinsol_mem) - KINFree(&kinsol_mem); + { + KINFree(&kinsol_mem); + +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + status = SUNContext_Free(&kinsol_ctx); + AssertKINSOL(status); +# endif + } + +# ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + Utilities::MPI::free_communicator(mpi_communicator); +# endif } @@ -330,21 +349,41 @@ namespace SUNDIALS VectorType u_scale_temp, f_scale_temp; if (get_solution_scaling) - u_scale = internal::make_nvector_view(get_solution_scaling()); + u_scale = internal::make_nvector_view(get_solution_scaling() +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + kinsol_ctx +# endif + ); else { reinit_vector(u_scale_temp); u_scale_temp = 1.0; - u_scale = internal::make_nvector_view(u_scale_temp); + u_scale = internal::make_nvector_view(u_scale_temp +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + kinsol_ctx +# endif + ); } if (get_function_scaling) - f_scale = internal::make_nvector_view(get_function_scaling()); + f_scale = internal::make_nvector_view(get_function_scaling() +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + kinsol_ctx +# endif + ); else { reinit_vector(f_scale_temp); f_scale_temp = 1.0; - f_scale = internal::make_nvector_view(f_scale_temp); + f_scale = internal::make_nvector_view(f_scale_temp +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + kinsol_ctx +# endif + ); } // Make sure we have what we need @@ -361,15 +400,34 @@ namespace SUNDIALS "solve_jacobian_system || solve_with_jacobian")); } - auto solution = internal::make_nvector_view(initial_guess_and_solution); + auto solution = internal::make_nvector_view(initial_guess_and_solution +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + kinsol_ctx +# endif + ); + + int status = 0; + (void)status; if (kinsol_mem) - KINFree(&kinsol_mem); + { + KINFree(&kinsol_mem); +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + status = SUNContext_Free(&kinsol_ctx); + AssertKINSOL(status); +# endif + } + +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) kinsol_mem = KINCreate(); +# else + status = SUNContext_Create(&mpi_communicator, &kinsol_ctx); + AssertKINSOL(status); - int status = 0; - (void)status; + kinsol_mem = KINCreate(kinsol_ctx); +# endif status = KINSetUserData(kinsol_mem, static_cast(this)); AssertKINSOL(status); @@ -443,7 +501,11 @@ namespace SUNDIALS // called do not actually receive the KINSOL object, just the LS // object, so we have to store a pointer to the current // object in the LS object +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) LS = SUNLinSolNewEmpty(); +# else + LS = SUNLinSolNewEmpty(kinsol_ctx); +# endif LS->content = this; LS->ops->gettype = @@ -473,7 +535,11 @@ namespace SUNDIALS // if we don't set it, it won't call the functions that set up // the matrix object (i.e., the argument to the 'KINSetJacFn' // function below). +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) J = SUNMatNewEmpty(); +# else + J = SUNMatNewEmpty(kinsol_ctx); +# endif J->content = this; J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID { diff --git a/source/sundials/n_vector.inst.in b/source/sundials/n_vector.inst.in index 899195b795..53757bffaa 100644 --- a/source/sundials/n_vector.inst.in +++ b/source/sundials/n_vector.inst.in @@ -17,9 +17,19 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) { template SUNDIALS::internal::NVectorView> - SUNDIALS::internal::make_nvector_view<>(V &); + SUNDIALS::internal::make_nvector_view<>(V & +#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +#endif + ); template SUNDIALS::internal::NVectorView> - SUNDIALS::internal::make_nvector_view<>(const V &); + SUNDIALS::internal::make_nvector_view<>(const V & +#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +#endif + ); template V * SUNDIALS::internal::unwrap_nvector>(N_Vector); template const V *SUNDIALS::internal::unwrap_nvector_const>( N_Vector); @@ -32,11 +42,21 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) { template SUNDIALS::internal::NVectorView> - SUNDIALS::internal::make_nvector_view<>(LinearAlgebra::distributed::V &); + SUNDIALS::internal::make_nvector_view<>(LinearAlgebra::distributed::V & +#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +#endif + ); template SUNDIALS::internal::NVectorView< const LinearAlgebra::distributed::V> SUNDIALS::internal::make_nvector_view<>( - const LinearAlgebra::distributed::V &); + const LinearAlgebra::distributed::V & +#if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +#endif + ); template LinearAlgebra::distributed::V *SUNDIALS::internal::unwrap_nvector>( N_Vector); @@ -54,9 +74,19 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) for (V : DEAL_II_VEC_TEMPLATES) { template SUNDIALS::internal::NVectorView - SUNDIALS::internal::make_nvector_view<>(TrilinosWrappers::MPI::V &); + SUNDIALS::internal::make_nvector_view<>(TrilinosWrappers::MPI::V & +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +# endif + ); template SUNDIALS::internal::NVectorView - SUNDIALS::internal::make_nvector_view<>(const TrilinosWrappers::MPI::V &); + SUNDIALS::internal::make_nvector_view<>(const TrilinosWrappers::MPI::V & +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +# endif + ); template TrilinosWrappers::MPI::V *SUNDIALS::internal::unwrap_nvector(N_Vector); template const TrilinosWrappers::MPI::V @@ -74,9 +104,19 @@ for (V : DEAL_II_VEC_TEMPLATES) for (V : DEAL_II_VEC_TEMPLATES) { template SUNDIALS::internal::NVectorView - SUNDIALS::internal::make_nvector_view<>(PETScWrappers::MPI::V &); + SUNDIALS::internal::make_nvector_view<>(PETScWrappers::MPI::V & +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +# endif + ); template SUNDIALS::internal::NVectorView - SUNDIALS::internal::make_nvector_view<>(const PETScWrappers::MPI::V &); + SUNDIALS::internal::make_nvector_view<>(const PETScWrappers::MPI::V & +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext +# endif + ); template PETScWrappers::MPI::V *SUNDIALS::internal::unwrap_nvector(N_Vector); template const PETScWrappers::MPI::V * diff --git a/source/sundials/sunlinsol_wrapper.cc b/source/sundials/sunlinsol_wrapper.cc index 78bc3db2fd..506be94b7e 100644 --- a/source/sundials/sunlinsol_wrapper.cc +++ b/source/sundials/sunlinsol_wrapper.cc @@ -176,10 +176,20 @@ namespace SUNDIALS template internal::LinearSolverWrapper::LinearSolverWrapper( - LinearSolveFunction lsolve) + LinearSolveFunction lsolve +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + SUNContext &linsol_ctx +# endif + ) : content(std::make_unique>()) { - sun_linear_solver = SUNLinSolNewEmpty(); +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + sun_linear_solver = SUNLinSolNewEmpty(); +# else + sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx); +# endif + sun_linear_solver->ops->gettype = arkode_linsol_get_type; sun_linear_solver->ops->solve = arkode_linsol_solve; sun_linear_solver->ops->setup = arkode_linsol_setup; @@ -230,11 +240,21 @@ 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); (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 } @@ -263,6 +283,7 @@ 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); // for custom preconditioners no distinction between left and right @@ -271,6 +292,15 @@ namespace SUNDIALS 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