From: Daniel Arndt Date: Fri, 5 Nov 2021 15:29:41 +0000 (-0400) Subject: Revert "Moved to new interface in sundials." X-Git-Tag: v9.3.2~1^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1ee88b067aac0b5ee5be33c66c3724075bf29d8e;p=dealii.git Revert "Moved to new interface in sundials." This reverts commit b108fe75d06f08db71b33b031f2ea57daa7f4251. --- diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 5d006d969f..43d52e16af 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -883,11 +883,32 @@ namespace SUNDIALS */ void *ida_mem; + /** + * IDA solution vector. + */ + N_Vector yy; + + /** + * IDA solution derivative vector. + */ + N_Vector yp; + + /** + * IDA absolute tolerances vector. + */ + N_Vector abs_tolls; + + /** + * IDA differential components vector. + */ + N_Vector diff_id; + /** * Number of iteration required to solve the Jacobian system */ int n_iter; + /** * MPI communicator. SUNDIALS solver runs happily in * parallel. Note that if the library is compiled without MPI diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index 027fd017f1..0d2acde788 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -220,10 +220,6 @@ namespace SUNDIALS realtype weighted_rms_norm(N_Vector x, N_Vector w); - template - realtype - weighted_rms_norm_mask(N_Vector x, N_Vector w, N_Vector mask); - template realtype max_norm(N_Vector x); @@ -669,24 +665,6 @@ SUNDIALS::internal::NVectorOperations::weighted_rms_norm(N_Vector x, N_Vector w) -template -realtype -SUNDIALS::internal::NVectorOperations::weighted_rms_norm_mask(N_Vector x, - N_Vector w, - N_Vector mask) -{ - // TODO copy can be avoided by a custom kernel - VectorType tmp = *unwrap_nvector_const(x); - auto * w_dealii = unwrap_nvector_const(w); - auto * mask_dealii = unwrap_nvector_const(mask); - const auto n = tmp.size(); - tmp.scale(*w_dealii); - tmp.scale(*mask_dealii); - return tmp.l2_norm() / std::sqrt(n); -} - - - template realtype SUNDIALS::internal::NVectorOperations::max_norm(N_Vector x) @@ -962,11 +940,11 @@ SUNDIALS::internal::create_empty_nvector() v->ops->nvdotprod = NVectorOperations::dot_product; 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->nvwrmsnormmask = - NVectorOperations::weighted_rms_norm_mask; + ; // v->ops->nvcompare = undef; // v->ops->nvinvtest = undef; // v->ops->nvconstrmask = undef; diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index bd207108cf..6788c62fca 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -14,17 +14,12 @@ //----------------------------------------------------------- #include -#include - #include #ifdef DEAL_II_WITH_SUNDIALS # include # include - -# include -# include # ifdef DEAL_II_WITH_TRILINOS # include # include @@ -34,6 +29,7 @@ # include # endif +# include # include # include @@ -67,13 +63,24 @@ namespace SUNDIALS void * user_data) { IDA &solver = *static_cast *>(user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_yy(mem); + solver.reinit_vector(*src_yy); - auto *src_yy = internal::unwrap_nvector_const(yy); - auto *src_yp = internal::unwrap_nvector_const(yp); - auto *residual = internal::unwrap_nvector(rr); + typename VectorMemory::Pointer src_yp(mem); + solver.reinit_vector(*src_yp); + + typename VectorMemory::Pointer residual(mem); + solver.reinit_vector(*residual); + + copy(*src_yy, yy); + copy(*src_yp, yp); int err = solver.residual(tt, *src_yy, *src_yp, *residual); + copy(rr, *residual); + return err; } @@ -96,9 +103,16 @@ namespace SUNDIALS (void)resp; IDA &solver = *static_cast *>(IDA_mem->ida_user_data); + GrowingVectorMemory mem; - auto *src_yy = internal::unwrap_nvector_const(yy); - auto *src_yp = internal::unwrap_nvector_const(yp); + typename VectorMemory::Pointer src_yy(mem); + solver.reinit_vector(*src_yy); + + typename VectorMemory::Pointer src_yp(mem); + solver.reinit_vector(*src_yp); + + copy(*src_yy, yy); + copy(*src_yp, yp); int err = solver.setup_jacobian(IDA_mem->ida_tn, *src_yy, @@ -127,13 +141,16 @@ namespace SUNDIALS *static_cast *>(IDA_mem->ida_user_data); GrowingVectorMemory mem; + typename VectorMemory::Pointer src(mem); + solver.reinit_vector(*src); + typename VectorMemory::Pointer dst(mem); solver.reinit_vector(*dst); - auto *src = internal::unwrap_nvector(b); + copy(*src, b); int err = solver.solve_jacobian_system(*src, *dst); - *src = *dst; + copy(b, *dst); return err; } @@ -156,12 +173,21 @@ namespace SUNDIALS { Assert(user_data != nullptr, ExcInternalError()); IDA &solver = *static_cast *>(user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_yy(mem); + solver.reinit_vector(*src_yy); + + typename VectorMemory::Pointer src_yp(mem); + solver.reinit_vector(*src_yp); + + copy(*src_yy, yy); + copy(*src_yp, yp); - auto *src_yy = internal::unwrap_nvector_const(yy); - auto *src_yp = internal::unwrap_nvector_const(yp); int err = solver.setup_jacobian(tt, *src_yy, *src_yp, cj); + return err; } @@ -178,11 +204,21 @@ namespace SUNDIALS const IDA &solver = *static_cast *>(LS->content); - auto *src_b = internal::unwrap_nvector_const(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_jacobian_system(*src_b, *dst_x); + copy(x, *dst_x); + return err; } @@ -197,12 +233,20 @@ namespace SUNDIALS { IDA &solver = *static_cast *>(LS->content); - auto *src_b = internal::unwrap_nvector_const(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); int n_iter; const int err = solver.solve_with_jacobian(*src_b, *dst_x, n_iter, tol); solver.set_n_iter(n_iter > 0 ? n_iter : 1); + copy(x, *dst_x); return err; } @@ -215,6 +259,10 @@ namespace SUNDIALS IDA::IDA(const AdditionalData &data, const MPI_Comm mpi_comm) : data(data) , ida_mem(nullptr) + , yy(nullptr) + , yp(nullptr) + , abs_tolls(nullptr) + , diff_id(nullptr) , communicator(is_serial_vector::value ? MPI_COMM_SELF : Utilities::MPI::duplicate_communicator(mpi_comm)) @@ -243,6 +291,8 @@ namespace SUNDIALS unsigned int IDA::solve_dae(VectorType &solution, VectorType &solution_dot) { + unsigned int system_size = solution.size(); + double t = data.initial_time; double h = data.initial_step_size; unsigned int step_number = 0; @@ -251,11 +301,36 @@ namespace SUNDIALS int status; (void)status; - reset(data.initial_time, data.initial_step_size, solution, solution_dot); - // 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 = solution.locally_owned_elements(); + const std::size_t local_system_size = is.n_elements(); + + yy = N_VNew_Parallel(communicator, local_system_size, system_size); + + yp = N_VNew_Parallel(communicator, local_system_size, system_size); + + diff_id = N_VNew_Parallel(communicator, local_system_size, system_size); + + abs_tolls = + N_VNew_Parallel(communicator, local_system_size, system_size); + } + else +# endif + { + Assert(is_serial_vector::value, + ExcInternalError( + "Trying to use a serial code with a parallel vector.")); + yy = N_VNew_Serial(system_size); + yp = N_VNew_Serial(system_size); + diff_id = N_VNew_Serial(system_size); + abs_tolls = N_VNew_Serial(system_size); + } + reset(data.initial_time, data.initial_step_size, solution, solution_dot); double next_time = data.initial_time; @@ -265,15 +340,15 @@ namespace SUNDIALS { next_time += data.output_period; - auto yy = internal::make_nvector_view(solution); - auto yp = internal::make_nvector_view(solution_dot); - status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL); AssertIDA(status); status = IDAGetLastStep(ida_mem, &h); AssertIDA(status); + copy(solution, yy); + copy(solution_dot, yp); + while (solver_should_restart(t, solution, solution_dot)) reset(t, h, solution, solution_dot); @@ -282,6 +357,24 @@ namespace SUNDIALS output_step(t, solution, solution_dot, step_number); } + // Free the vectors which are no longer used. +# ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + { + N_VDestroy_Parallel(yy); + N_VDestroy_Parallel(yp); + N_VDestroy_Parallel(abs_tolls); + N_VDestroy_Parallel(diff_id); + } + else +# endif + { + N_VDestroy_Serial(yy); + N_VDestroy_Serial(yp); + N_VDestroy_Serial(abs_tolls); + N_VDestroy_Serial(diff_id); + } + return step_number; } @@ -294,26 +387,71 @@ namespace SUNDIALS VectorType & solution, VectorType & solution_dot) { - bool first_step = (current_time == data.initial_time); + unsigned int system_size; + bool first_step = (current_time == data.initial_time); if (ida_mem) IDAFree(&ida_mem); ida_mem = IDACreate(); + // Free the vectors which are no longer used. + if (yy) + { +# ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + { + N_VDestroy_Parallel(yy); + N_VDestroy_Parallel(yp); + N_VDestroy_Parallel(abs_tolls); + N_VDestroy_Parallel(diff_id); + } + else +# endif + { + N_VDestroy_Serial(yy); + N_VDestroy_Serial(yp); + N_VDestroy_Serial(abs_tolls); + N_VDestroy_Serial(diff_id); + } + } + int status; (void)status; + system_size = solution.size(); +# ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + { + const IndexSet is = solution.locally_owned_elements(); + const std::size_t local_system_size = is.n_elements(); + + yy = N_VNew_Parallel(communicator, local_system_size, system_size); - auto yy = internal::make_nvector_view(solution); - auto yp = internal::make_nvector_view(solution_dot); + yp = N_VNew_Parallel(communicator, local_system_size, system_size); + + diff_id = N_VNew_Parallel(communicator, local_system_size, system_size); + + abs_tolls = + N_VNew_Parallel(communicator, local_system_size, system_size); + } + else +# endif + { + yy = N_VNew_Serial(system_size); + yp = N_VNew_Serial(system_size); + diff_id = N_VNew_Serial(system_size); + abs_tolls = N_VNew_Serial(system_size); + } + + copy(yy, solution); + copy(yp, solution_dot); status = IDAInit(ida_mem, t_dae_residual, current_time, yy, yp); AssertIDA(status); if (get_local_tolerances) { - const auto abs_tols = - internal::make_nvector_view(get_local_tolerances()); - status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols); + copy(abs_tolls, get_local_tolerances()); + status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tolls); AssertIDA(status); } else @@ -339,10 +477,9 @@ namespace SUNDIALS auto dc = differential_components(); for (auto i = dc.begin(); i != dc.end(); ++i) 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); + copy(diff_id, diff_comp_vector); + status = IDASetId(ida_mem, diff_id); AssertIDA(status); } @@ -488,6 +625,9 @@ namespace SUNDIALS status = IDAGetConsistentIC(ida_mem, yy, yp); AssertIDA(status); + + copy(solution, yy); + copy(solution_dot, yp); } else if (type == AdditionalData::use_y_diff) { @@ -497,6 +637,9 @@ namespace SUNDIALS status = IDAGetConsistentIC(ida_mem, yy, yp); AssertIDA(status); + + copy(solution, yy); + copy(solution_dot, yp); } }