From: Daniel Arndt Date: Fri, 5 Nov 2021 15:29:45 +0000 (-0400) Subject: Revert "Remove copy from SUNDIALS::KINSOL" X-Git-Tag: v9.3.2~1^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=140856483c40ae91e63ba4a045ebb2305026542d;p=dealii.git Revert "Remove copy from SUNDIALS::KINSOL" This reverts commit 422d57b78a191bec3eccbb9258076c21222b3037. --- diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 939215b05f..8a33c0ea27 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -711,6 +711,21 @@ 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 0d2acde788..2b49e3efc9 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -168,14 +168,6 @@ 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); @@ -577,28 +569,6 @@ 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) @@ -941,10 +911,9 @@ 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 = NVectorOperations::weighted_l2_norm; - v->ops->nvl1norm = NVectorOperations::l1_norm; - ; + v->ops->nvmin = NVectorOperations::min_element; + // v->ops->nvwl2norm = undef; + // v->ops->nvl1norm = undef; // 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 cb12781b0e..ce3cb44077 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,9 +137,15 @@ namespace SUNDIALS { KINSOL &solver = *static_cast *>(user_data); + GrowingVectorMemory mem; - auto *src_yy = internal::unwrap_nvector_const(yy); - auto *dst_FF = internal::unwrap_nvector(FF); + 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); int err = 0; if (solver.residual) @@ -149,6 +155,8 @@ namespace SUNDIALS else Assert(false, ExcInternalError()); + copy(FF, *dst_FF); + return err; } @@ -161,11 +169,16 @@ namespace SUNDIALS { KINSOL &solver = *static_cast *>(kinsol_mem->kin_user_data); + GrowingVectorMemory mem; - auto *src_ycur = - internal::unwrap_nvector_const(kinsol_mem->kin_uu); - auto *src_fcur = - internal::unwrap_nvector_const(kinsol_mem->kin_fval); + 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); int err = solver.setup_jacobian(*src_ycur, *src_fcur); return err; @@ -183,15 +196,27 @@ namespace SUNDIALS { KINSOL &solver = *static_cast *>(kinsol_mem->kin_user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_ycur(mem); + solver.reinit_vector(*src_ycur); - 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); + 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); 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); @@ -218,8 +243,16 @@ namespace SUNDIALS const KINSOL &solver = *static_cast *>(user_data); - auto *ycur = internal::unwrap_nvector_const(u); - auto *fcur = internal::unwrap_nvector_const(f); + // 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); // Call the user-provided setup function with these arguments: solver.setup_jacobian(*ycur, *fcur); @@ -247,11 +280,21 @@ namespace SUNDIALS // signals to call. Let's first check the more modern one: if (solver.solve_with_jacobian) { - auto *src_b = internal::unwrap_nvector(b); - auto *dst_x = internal::unwrap_nvector(x); + // 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); const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol); + copy(x, *dst_x); + return err; } else @@ -266,9 +309,13 @@ 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); - auto *src_b = internal::unwrap_nvector_const(b); - auto *dst_x = internal::unwrap_nvector(x); + copy(*src_b, b); // Call the user-provided setup function with these arguments. Note // that Sundials 4.x and later no longer provide values for @@ -278,6 +325,8 @@ namespace SUNDIALS const int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x); + copy(x, *dst_x); + return err; } } @@ -292,6 +341,9 @@ 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)) @@ -325,29 +377,44 @@ namespace SUNDIALS { unsigned int system_size = initial_guess_and_solution.size(); - NVectorView u_scale, f_scale; + // 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(); - VectorType u_scale_temp, f_scale_temp; + solution = + N_VNew_Parallel(communicator, local_system_size, system_size); - if (get_solution_scaling) - u_scale = internal::make_nvector_view(get_solution_scaling()); + u_scale = N_VNew_Parallel(communicator, local_system_size, system_size); + N_VConst_Parallel(1.e0, u_scale); + + f_scale = N_VNew_Parallel(communicator, local_system_size, system_size); + N_VConst_Parallel(1.e0, f_scale); + } else +# endif { - reinit_vector(u_scale_temp); - u_scale_temp = 1.0; - u_scale = internal::make_nvector_view(u_scale_temp); + 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); } + if (get_solution_scaling) + copy(u_scale, get_solution_scaling()); + if (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(f_scale, get_function_scaling()); - auto solution = internal::make_nvector_view(initial_guess_and_solution); + copy(solution, initial_guess_and_solution); if (kinsol_mem) KINFree(&kinsol_mem); @@ -506,6 +573,24 @@ 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);