From 422d57b78a191bec3eccbb9258076c21222b3037 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 16 May 2021 11:14:16 +0200 Subject: [PATCH] Remove copy from SUNDIALS::KINSOL --- include/deal.II/sundials/kinsol.h | 15 -- include/deal.II/sundials/n_vector.templates.h | 37 ++++- source/sundials/kinsol.cc | 153 ++++-------------- 3 files changed, 68 insertions(+), 137 deletions(-) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 8a33c0ea27..939215b05f 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -711,21 +711,6 @@ namespace SUNDIALS */ void *kinsol_mem; - /** - * KINSOL solution vector. - */ - N_Vector solution; - - /** - * KINSOL solution scale. - */ - N_Vector u_scale; - - /** - * KINSOL f scale. - */ - N_Vector f_scale; - /** * MPI communicator. SUNDIALS solver runs happily in parallel. */ diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index a6d6562ff6..49ea702018 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -168,6 +168,14 @@ namespace SUNDIALS realtype dot_product(N_Vector x, N_Vector y); + template + realtype + weighted_l2_norm(N_Vector x, N_Vector y); + + template + realtype + l1_norm(N_Vector x); + template void elementwise_product(N_Vector x, N_Vector y, N_Vector z); @@ -569,6 +577,28 @@ SUNDIALS::internal::NVectorOperations::dot_product(N_Vector x, N_Vector y) +template +realtype +SUNDIALS::internal::NVectorOperations::weighted_l2_norm(N_Vector x, N_Vector w) +{ + // TODO copy can be avoided by a custom kernel + VectorType tmp = *unwrap_nvector_const(x); + auto * w_dealii = unwrap_nvector_const(w); + tmp.scale(*w_dealii); + return tmp.l2_norm(); +} + + + +template +realtype +SUNDIALS::internal::NVectorOperations::l1_norm(N_Vector x) +{ + return unwrap_nvector_const(x)->l1_norm(); +} + + + template void SUNDIALS::internal::NVectorOperations::set_constant(realtype c, N_Vector v) @@ -911,9 +941,10 @@ SUNDIALS::internal::create_empty_nvector() v->ops->nvmaxnorm = NVectorOperations::max_norm; v->ops->nvwrmsnorm = NVectorOperations::weighted_rms_norm; // v->ops->nvwrmsnormmask = undef; - v->ops->nvmin = NVectorOperations::min_element; - // v->ops->nvwl2norm = undef; - // v->ops->nvl1norm = undef; + v->ops->nvmin = NVectorOperations::min_element; + v->ops->nvwl2norm = NVectorOperations::weighted_l2_norm; + v->ops->nvl1norm = NVectorOperations::l1_norm; + ; // v->ops->nvcompare = undef; // v->ops->nvinvtest = undef; // v->ops->nvconstrmask = undef; diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index ce3cb44077..cb12781b0e 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -34,7 +34,7 @@ # include # endif -# include +# include // Make sure we #include the SUNDIALS config file... # include @@ -137,15 +137,9 @@ namespace SUNDIALS { 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); - - copy(*src_yy, yy); + auto *src_yy = internal::unwrap_nvector_const(yy); + auto *dst_FF = internal::unwrap_nvector(FF); int err = 0; if (solver.residual) @@ -155,8 +149,6 @@ namespace SUNDIALS else Assert(false, ExcInternalError()); - copy(FF, *dst_FF); - return err; } @@ -169,16 +161,11 @@ namespace SUNDIALS { KINSOL &solver = *static_cast *>(kinsol_mem->kin_user_data); - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_ycur(mem); - solver.reinit_vector(*src_ycur); - - typename VectorMemory::Pointer src_fcur(mem); - solver.reinit_vector(*src_fcur); - - copy(*src_ycur, kinsol_mem->kin_uu); - copy(*src_fcur, kinsol_mem->kin_fval); + auto *src_ycur = + internal::unwrap_nvector_const(kinsol_mem->kin_uu); + auto *src_fcur = + internal::unwrap_nvector_const(kinsol_mem->kin_fval); int err = solver.setup_jacobian(*src_ycur, *src_fcur); return err; @@ -196,27 +183,15 @@ namespace SUNDIALS { KINSOL &solver = *static_cast *>(kinsol_mem->kin_user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src_ycur(mem); - solver.reinit_vector(*src_ycur); - typename VectorMemory::Pointer src_fcur(mem); - solver.reinit_vector(*src_fcur); - - copy(*src_ycur, kinsol_mem->kin_uu); - copy(*src_fcur, kinsol_mem->kin_fval); - - typename VectorMemory::Pointer src(mem); - solver.reinit_vector(*src); - - typename VectorMemory::Pointer dst(mem); - solver.reinit_vector(*dst); - - copy(*src, b); + auto *src_ycur = + internal::unwrap_nvector_const(kinsol_mem->kin_uu); + auto *src_fcur = + internal::unwrap_nvector_const(kinsol_mem->kin_fval); + auto *src = internal::unwrap_nvector_const(b); + auto *dst = internal::unwrap_nvector(x); int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src, *dst); - copy(x, *dst); *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale); N_VProd(b, kinsol_mem->kin_fscale, b); @@ -243,16 +218,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); - - copy(*ycur, u); - copy(*fcur, f); + auto *ycur = internal::unwrap_nvector_const(u); + auto *fcur = internal::unwrap_nvector_const(f); // Call the user-provided setup function with these arguments: solver.setup_jacobian(*ycur, *fcur); @@ -280,21 +247,11 @@ 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); - - copy(*src_b, b); + auto *src_b = internal::unwrap_nvector(b); + auto *dst_x = internal::unwrap_nvector(x); const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol); - copy(x, *dst_x); - return err; } else @@ -309,13 +266,9 @@ namespace SUNDIALS 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); - 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 @@ -325,8 +278,6 @@ namespace SUNDIALS const int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x); - copy(x, *dst_x); - return err; } } @@ -341,9 +292,6 @@ namespace SUNDIALS const MPI_Comm & mpi_comm) : data(data) , kinsol_mem(nullptr) - , solution(nullptr) - , u_scale(nullptr) - , f_scale(nullptr) , communicator(is_serial_vector::value ? MPI_COMM_SELF : Utilities::MPI::duplicate_communicator(mpi_comm)) @@ -377,44 +325,29 @@ namespace SUNDIALS { unsigned int 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 unsigned int local_system_size = is.n_elements(); - - solution = - N_VNew_Parallel(communicator, local_system_size, system_size); + NVectorView u_scale, f_scale; - u_scale = N_VNew_Parallel(communicator, local_system_size, system_size); - N_VConst_Parallel(1.e0, u_scale); + VectorType u_scale_temp, f_scale_temp; - f_scale = N_VNew_Parallel(communicator, local_system_size, system_size); - N_VConst_Parallel(1.e0, f_scale); - } + if (get_solution_scaling) + u_scale = internal::make_nvector_view(get_solution_scaling()); else -# endif { - Assert(is_serial_vector::value, - ExcInternalError( - "Trying to use a serial code with a parallel vector.")); - solution = N_VNew_Serial(system_size); - u_scale = N_VNew_Serial(system_size); - N_VConst_Serial(1.e0, u_scale); - f_scale = N_VNew_Serial(system_size); - N_VConst_Serial(1.e0, f_scale); + reinit_vector(u_scale_temp); + u_scale_temp = 1.0; + u_scale = internal::make_nvector_view(u_scale_temp); } - if (get_solution_scaling) - copy(u_scale, get_solution_scaling()); - if (get_function_scaling) - copy(f_scale, get_function_scaling()); + f_scale = internal::make_nvector_view(get_function_scaling()); + else + { + reinit_vector(f_scale_temp); + f_scale_temp = 1.0; + f_scale = internal::make_nvector_view(f_scale_temp); + } - copy(solution, initial_guess_and_solution); + auto solution = internal::make_nvector_view(initial_guess_and_solution); if (kinsol_mem) KINFree(&kinsol_mem); @@ -573,24 +506,6 @@ namespace SUNDIALS status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale); AssertKINSOL(status); - copy(initial_guess_and_solution, solution); - - // 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); - } - long nniters; status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters); AssertKINSOL(status); -- 2.39.5