From 4f8847b76c75d4bd5efb298a3befabfbead5ca56 Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Thu, 11 May 2023 16:19:48 +0200 Subject: [PATCH] KINSOL: do not copy vectors internally --- include/deal.II/sundials/kinsol.h | 14 --- source/sundials/kinsol.cc | 177 ++++++------------------------ 2 files changed, 34 insertions(+), 157 deletions(-) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 994f2b225a..56ce7b6994 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -739,20 +739,6 @@ namespace SUNDIALS SUNContext kinsol_ctx; # endif - /** - * KINSOL solution vector. - */ - N_Vector solution; - - /** - * KINSOL solution scale. - */ - N_Vector u_scale; - - /** - * KINSOL f scale. - */ - N_Vector f_scale; /** * Memory pool of vectors. diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index 1a89c03e6e..a2bd709493 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -36,7 +36,6 @@ # include # endif -# include # include // Make sure we #include the SUNDIALS config file... @@ -217,9 +216,6 @@ namespace SUNDIALS # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , kinsol_ctx(nullptr) # endif - , solution(nullptr) - , u_scale(nullptr) - , f_scale(nullptr) , pending_exception(nullptr) { set_functions_to_trigger_an_assert(); @@ -300,78 +296,32 @@ namespace SUNDIALS status = KINSetUserData(kinsol_mem, static_cast(this)); AssertKINSOL(status); - const auto system_size = initial_guess_and_solution.size(); - - // The solution is stored in solution. Here we take only a view of it. -# ifdef DEAL_II_WITH_MPI - if (is_serial_vector::value == false) - { - const IndexSet &is = - initial_guess_and_solution.locally_owned_elements(); - const auto local_system_size = is.n_elements(); - - solution = N_VNew_Parallel(mpi_communicator, - local_system_size, - system_size -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - - u_scale = N_VNew_Parallel(mpi_communicator, - local_system_size, - system_size -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - N_VConst_Parallel(1.e0, u_scale); - - f_scale = N_VNew_Parallel(mpi_communicator, - local_system_size, - system_size -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - N_VConst_Parallel(1.e0, f_scale); - } - else + // helper function to create N_Vectors compatible with different versions + // of SUNDIALS + const auto make_compatible_nvector_view = [this](auto &v) { + return internal::make_nvector_view(v +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + , + kinsol_ctx # endif + ); + }; + + + VectorType ones; + // Prepare constant vector for scaling f or u only if we need it + if (!get_function_scaling || !get_solution_scaling) { - Assert(is_serial_vector::value, ExcInternalError()); - solution = N_VNew_Serial(system_size -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - u_scale = N_VNew_Serial(system_size -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - N_VConst_Serial(1.e0, u_scale); - f_scale = N_VNew_Serial(system_size -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - N_VConst_Serial(1.e0, f_scale); + reinit_vector(ones); + ones = 1.0; } - if (get_solution_scaling) - internal::copy(u_scale, get_solution_scaling()); - - if (get_function_scaling) - internal::copy(f_scale, get_function_scaling()); + auto u_scale = make_compatible_nvector_view( + get_solution_scaling ? get_solution_scaling() : ones); + auto f_scale = make_compatible_nvector_view( + get_function_scaling ? get_function_scaling() : ones); - internal::copy(solution, initial_guess_and_solution); + auto solution = make_compatible_nvector_view(initial_guess_and_solution); // This must be called before KINSetMAA status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations); @@ -388,15 +338,9 @@ namespace SUNDIALS [](N_Vector yy, N_Vector FF, void *user_data) -> int { KINSOL &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src_yy(mem); - solver.reinit_vector(*src_yy); - - typename VectorMemory::Pointer dst_FF(mem); - solver.reinit_vector(*dst_FF); - internal::copy(*src_yy, yy); + auto src_yy = internal::unwrap_nvector_const(yy); + auto dst_FF = internal::unwrap_nvector(FF); Assert(solver.iteration_function, ExcInternalError()); @@ -406,8 +350,6 @@ namespace SUNDIALS *src_yy, *dst_FF); - internal::copy(FF, *dst_FF); - return err; }, solution); @@ -418,23 +360,15 @@ namespace SUNDIALS [](N_Vector yy, N_Vector FF, void *user_data) -> int { KINSOL &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src_yy(mem); - solver.reinit_vector(*src_yy); - - typename VectorMemory::Pointer dst_FF(mem); - solver.reinit_vector(*dst_FF); - internal::copy(*src_yy, yy); + auto src_yy = internal::unwrap_nvector_const(yy); + auto dst_FF = internal::unwrap_nvector(FF); Assert(solver.residual, ExcInternalError()); const int err = call_and_possibly_capture_exception( solver.residual, solver.pending_exception, *src_yy, *dst_FF); - internal::copy(FF, *dst_FF); - return err; }, solution); @@ -465,8 +399,8 @@ namespace SUNDIALS SUNLinearSolver LS = nullptr; if (solve_jacobian_system || - solve_with_jacobian) // user assigned a function object to the solver - // slot + solve_with_jacobian) // user assigned a function + // object to the solver slot { // Set the operations we care for in the sun_linear_solver object // and attach it to the KINSOL object. The functions that will get @@ -515,16 +449,8 @@ namespace SUNDIALS // signals to call. Let's first check the more modern one: if (solver.solve_with_jacobian) { - // Allocate temporary (deal.II-type) vectors into which to copy - // the N_vectors - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_b(mem); - typename VectorMemory::Pointer dst_x(mem); - - solver.reinit_vector(*src_b); - solver.reinit_vector(*dst_x); - - internal::copy(*src_b, b); + auto src_b = internal::unwrap_nvector_const(b); + auto dst_x = internal::unwrap_nvector(x); const int err = call_and_possibly_capture_exception(solver.solve_with_jacobian, @@ -533,8 +459,6 @@ namespace SUNDIALS *dst_x, tol); - internal::copy(x, *dst_x); - return err; } else @@ -544,18 +468,13 @@ namespace SUNDIALS // the old signal. Check this. Assert(solver.solve_jacobian_system, ExcInternalError()); - // Allocate temporary (deal.II-type) vectors into which to copy - // the N_vectors + // Allocate temporary (deal.II-type) dummy vectors GrowingVectorMemory mem; typename VectorMemory::Pointer src_ycur(mem); typename VectorMemory::Pointer src_fcur(mem); - typename VectorMemory::Pointer src_b(mem); - typename VectorMemory::Pointer dst_x(mem); - solver.reinit_vector(*src_b); - solver.reinit_vector(*dst_x); - - internal::copy(*src_b, b); + auto src_b = internal::unwrap_nvector_const(b); + auto dst_x = internal::unwrap_nvector(x); // Call the user-provided setup function with these arguments. // Note that Sundials 4.x and later no longer provide values for @@ -570,8 +489,6 @@ namespace SUNDIALS *src_b, *dst_x); - internal::copy(x, *dst_x); - return err; } }; @@ -631,16 +548,8 @@ namespace SUNDIALS const KINSOL &solver = *static_cast *>(user_data); - // Allocate temporary (deal.II-type) vectors into which to copy the - // N_vectors - GrowingVectorMemory mem; - typename VectorMemory::Pointer ycur(mem); - typename VectorMemory::Pointer fcur(mem); - solver.reinit_vector(*ycur); - solver.reinit_vector(*fcur); - - internal::copy(*ycur, u); - internal::copy(*fcur, f); + auto ycur = internal::unwrap_nvector_const(u); + auto fcur = internal::unwrap_nvector(f); // Call the user-provided setup function with these arguments: return call_and_possibly_capture_exception(solver.setup_jacobian, @@ -669,22 +578,6 @@ namespace SUNDIALS status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale); ScopeExit upon_exit([this, &J, &LS]() mutable { - // Free the vectors which are no longer used. -# ifdef DEAL_II_WITH_MPI - if (is_serial_vector::value == false) - { - N_VDestroy_Parallel(solution); - N_VDestroy_Parallel(u_scale); - N_VDestroy_Parallel(f_scale); - } - else -# endif - { - N_VDestroy_Serial(solution); - N_VDestroy_Serial(u_scale); - N_VDestroy_Serial(f_scale); - } - if (J != nullptr) SUNMatDestroy(J); if (LS != nullptr) @@ -714,8 +607,6 @@ namespace SUNDIALS } AssertKINSOL(status); - internal::copy(initial_guess_and_solution, solution); - long nniters; status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters); AssertKINSOL(status); -- 2.39.5