From 5c4c71cbc9d2e413256c996f98c01a4b99914609 Mon Sep 17 00:00:00 2001 From: David Wells Date: Fri, 3 Jun 2022 11:03:39 -0400 Subject: [PATCH] Revert SUNDIALS copy to fix KINSOL, part 1 --- include/deal.II/sundials/kinsol.h | 15 +++ source/sundials/kinsol.cc | 205 +++++++++++++++++++++--------- 2 files changed, 161 insertions(+), 59 deletions(-) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 5b2c50cd19..9bcdfda83c 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -691,6 +691,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; + # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) /** * A context object associated with the KINSOL solver. diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index 66b40dfc84..921755853f 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -34,6 +34,7 @@ # include # endif +# include # include // Make sure we #include the SUNDIALS config file... @@ -134,9 +135,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); + + internal::copy(*src_yy, yy); int err = 0; if (solver.residual) @@ -144,6 +151,8 @@ namespace SUNDIALS else Assert(false, ExcInternalError()); + internal::copy(FF, *dst_FF); + return err; } @@ -155,9 +164,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); + + internal::copy(*src_yy, yy); int err = 0; if (solver.iteration_function) @@ -165,6 +180,8 @@ namespace SUNDIALS else Assert(false, ExcInternalError()); + internal::copy(FF, *dst_FF); + return err; } @@ -185,8 +202,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); + + internal::copy(*ycur, u); + internal::copy(*fcur, f); // Call the user-provided setup function with these arguments: solver.setup_jacobian(*ycur, *fcur); @@ -214,11 +239,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); + + internal::copy(*src_b, b); const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol); + internal::copy(x, *dst_x); + return err; } else @@ -233,9 +268,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); + internal::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 @@ -245,6 +284,8 @@ namespace SUNDIALS const int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x); + internal::copy(x, *dst_x); + return err; } } @@ -260,6 +301,9 @@ namespace SUNDIALS MPI_COMM_SELF : Utilities::MPI::duplicate_communicator(mpi_comm)) , kinsol_mem(nullptr) + , solution(nullptr) + , u_scale(nullptr) + , f_scale(nullptr) { set_functions_to_trigger_an_assert(); @@ -299,48 +343,6 @@ namespace SUNDIALS unsigned int KINSOL::solve(VectorType &initial_guess_and_solution) { - internal::NVectorView u_scale, f_scale; - - VectorType u_scale_temp, f_scale_temp; - - if (get_solution_scaling) - u_scale = internal::make_nvector_view(get_solution_scaling() -# if DEAL_II_SUNDIALS_VERSION_GTE(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 -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - } - - if (get_function_scaling) - f_scale = internal::make_nvector_view(get_function_scaling() -# if DEAL_II_SUNDIALS_VERSION_GTE(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 -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - } - // Make sure we have what we need if (data.strategy == AdditionalData::fixed_point) { @@ -355,13 +357,7 @@ namespace SUNDIALS "solve_jacobian_system || solve_with_jacobian")); } - auto solution = internal::make_nvector_view(initial_guess_and_solution -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - kinsol_ctx -# endif - ); - + // Create a new solver object: int status = 0; (void)status; @@ -387,6 +383,79 @@ 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 +# endif + { + 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); + } + + if (get_solution_scaling) + internal::copy(u_scale, get_solution_scaling()); + + if (get_function_scaling) + internal::copy(f_scale, get_function_scaling()); + + internal::copy(solution, initial_guess_and_solution); + // This must be called before KINSetMAA status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations); AssertKINSOL(status); @@ -512,6 +581,24 @@ namespace SUNDIALS status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale); AssertKINSOL(status); + internal::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