From: Sebastian Proell Date: Sat, 30 Jan 2021 11:23:31 +0000 (+0100) Subject: SUNDIALS: Use N_Vector module in ARKode wrapper X-Git-Tag: v9.3.0-rc1~532^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11655%2Fhead;p=dealii.git SUNDIALS: Use N_Vector module in ARKode wrapper --- diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index 5f565e3640..bd5bda0dcf 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -41,6 +41,8 @@ # ifdef DEAL_II_WITH_MPI # include # endif +# include + # include # include @@ -610,26 +612,7 @@ namespace SUNDIALS * @param[in,out] y The new initial solution */ void - reset(const double t, const double h, const VectorType &y); - - /*! - * Create a new SUNDIALS vector from a given template. - * - * @note Vectors created this way should be freed with free_vector(). - * - * @param template_vector The vector to use as a template for the layout of - * a new vector. - */ - N_Vector - create_vector(const VectorType &template_vector) const; - - /** - * Free a SUNDIALS vector created with create_vector(). - * - * @param vector The vector to free - */ - void - free_vector(N_Vector vector) const; + reset(const double t, const double h, VectorType &y); /** * Provides user access to the internally used ARKODE memory. @@ -1299,7 +1282,7 @@ namespace SUNDIALS /** * Set up the (non)linear solver and preconditioners in the ARKODE memory * object based on the user-specified functions. - * @param solution The solution vector whihc is used as a template to create + * @param solution The solution vector which is used as a template to create * new vectors. */ void @@ -1335,12 +1318,12 @@ namespace SUNDIALS /** * ARKode solution vector. */ - N_Vector yy; + internal::NVectorView yy; /** * ARKode absolute tolerances vector. */ - N_Vector abs_tolls; + internal::NVectorView abs_tolls; /** * MPI communicator. SUNDIALS solver runs happily in diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index 0a52dfeced..7cf6d94a4d 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -33,7 +33,7 @@ # include # endif -# include +# include # if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) # include @@ -72,21 +72,11 @@ namespace SUNDIALS Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src_yy(mem); - solver.reinit_vector(*src_yy); - - typename VectorMemory::Pointer dst_yp(mem); - solver.reinit_vector(*dst_yp); - - copy(*src_yy, yy); - int err = solver.explicit_function(tt, *src_yy, *dst_yp); + auto *src_yy = internal::unwrap_nvector_const(yy); + auto *dst_yp = internal::unwrap_nvector(yp); - copy(yp, *dst_yp); - - return err; + return solver.explicit_function(tt, *src_yy, *dst_yp); } @@ -101,21 +91,11 @@ namespace SUNDIALS Assert(user_data != nullptr, ExcInternalError()); ARKode &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 *dst_yp = internal::unwrap_nvector(yp); - typename VectorMemory::Pointer dst_yp(mem); - solver.reinit_vector(*dst_yp); - - copy(*src_yy, yy); - - int err = solver.implicit_function(tt, *src_yy, *dst_yp); - - copy(yp, *dst_yp); - - return err; + return solver.implicit_function(tt, *src_yy, *dst_yp); } @@ -135,17 +115,9 @@ namespace SUNDIALS Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(arkode_mem->ark_user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src_ypred(mem); - solver.reinit_vector(*src_ypred); - - typename VectorMemory::Pointer src_fpred(mem); - solver.reinit_vector(*src_fpred); - - copy(*src_ypred, ypred); - copy(*src_fpred, fpred); + auto *src_ypred = internal::unwrap_nvector_const(ypred); + auto *src_fpred = internal::unwrap_nvector_const(fpred); // avoid reinterpret_cast bool jcurPtr_tmp = false; int err = solver.setup_jacobian(convfail, @@ -178,31 +150,21 @@ namespace SUNDIALS Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(arkode_mem->ark_user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src(mem); - solver.reinit_vector(*src); - typename VectorMemory::Pointer src_ycur(mem); - solver.reinit_vector(*src_ycur); + auto *dst = internal::unwrap_nvector(b); + auto *src_ycur = internal::unwrap_nvector_const(ycur); + auto *src_fcur = internal::unwrap_nvector_const(fcur); - typename VectorMemory::Pointer src_fcur(mem); - solver.reinit_vector(*src_fcur); - - typename VectorMemory::Pointer dst(mem); - solver.reinit_vector(*dst); - - copy(*src, b); - copy(*src_ycur, ycur); - copy(*src_fcur, fcur); + // make a temporary copy to work on in the user call + VectorType src = *dst; int err = solver.solve_jacobian_system(arkode_mem->ark_tn, arkode_mem->ark_gamma, *src_ycur, *src_fcur, - *src, + src, *dst); - copy(b, *dst); + return err; } @@ -236,18 +198,13 @@ namespace SUNDIALS Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(arkode_mem->ark_user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer src(mem); - solver.reinit_vector(*src); - typename VectorMemory::Pointer dst(mem); - solver.reinit_vector(*dst); + auto *dst = internal::unwrap_nvector(b); - copy(*src, b); + // make a temporary copy to work on in the user call + VectorType src = *dst; - int err = solver.solve_mass_system(*src, *dst); - copy(b, *dst); + int err = solver.solve_mass_system(src, *dst); return err; } @@ -268,22 +225,13 @@ namespace SUNDIALS *static_cast *>(user_data); GrowingVectorMemory mem; - typename VectorMemory::Pointer src_v(mem); - solver.reinit_vector(*src_v); - typename VectorMemory::Pointer dst_Jv(mem); - solver.reinit_vector(*dst_Jv); - typename VectorMemory::Pointer src_y(mem); - solver.reinit_vector(*src_y); - typename VectorMemory::Pointer src_fy(mem); - solver.reinit_vector(*src_fy); - copy(*src_v, v); - copy(*src_y, y); - copy(*src_fy, fy); - - int err = - solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy); - copy(Jv, *dst_Jv); - return err; + auto *src_v = internal::unwrap_nvector_const(v); + auto *src_y = internal::unwrap_nvector_const(y); + auto *src_fy = internal::unwrap_nvector_const(fy); + + auto *dst_Jv = internal::unwrap_nvector(Jv); + + return solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy); } @@ -298,15 +246,11 @@ namespace SUNDIALS Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_y(mem); - solver.reinit_vector(*src_y); - typename VectorMemory::Pointer src_fy(mem); - solver.reinit_vector(*src_fy); - copy(*src_y, y); - copy(*src_fy, fy); - int err = solver.jacobian_times_setup(t, *src_y, *src_fy); - return err; + + auto *src_y = internal::unwrap_nvector_const(y); + auto *src_fy = internal::unwrap_nvector_const(fy); + + return solver.jacobian_times_setup(t, *src_y, *src_fy); } @@ -327,23 +271,14 @@ namespace SUNDIALS ARKode &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_y(mem); - solver.reinit_vector(*src_y); - typename VectorMemory::Pointer src_fy(mem); - solver.reinit_vector(*src_fy); - typename VectorMemory::Pointer src_r(mem); - solver.reinit_vector(*src_r); - typename VectorMemory::Pointer dst_z(mem); - solver.reinit_vector(*dst_z); - copy(*src_y, y); - copy(*src_fy, fy); - copy(*src_r, r); - - int status = solver.jacobian_preconditioner_solve( + auto *src_y = internal::unwrap_nvector_const(y); + auto *src_fy = internal::unwrap_nvector_const(fy); + auto *src_r = internal::unwrap_nvector_const(r); + + auto *dst_z = internal::unwrap_nvector(z); + + return solver.jacobian_preconditioner_solve( t, *src_y, *src_fy, *src_r, *dst_z, gamma, delta, lr); - copy(z, *dst_z); - return status; } @@ -361,13 +296,10 @@ namespace SUNDIALS Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_y(mem); - solver.reinit_vector(*src_y); - typename VectorMemory::Pointer src_fy(mem); - solver.reinit_vector(*src_fy); - copy(*src_y, y); - copy(*src_fy, fy); + + auto *src_y = internal::unwrap_nvector_const(y); + auto *src_fy = internal::unwrap_nvector_const(fy); + return solver.jacobian_preconditioner_setup( t, *src_y, *src_fy, jok, *jcurPtr, gamma); } @@ -384,17 +316,11 @@ namespace SUNDIALS Assert(mtimes_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(mtimes_data); - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_v(mem); - solver.reinit_vector(*src_v); - typename VectorMemory::Pointer dst_Mv(mem); - solver.reinit_vector(*dst_Mv); - copy(*src_v, v); - int err = solver.mass_times_vector(t, *src_v, *dst_Mv); - copy(Mv, *dst_Mv); + auto *src_v = internal::unwrap_nvector_const(v); + auto *dst_Mv = internal::unwrap_nvector(Mv); - return err; + return solver.mass_times_vector(t, *src_v, *dst_Mv); } @@ -425,17 +351,10 @@ namespace SUNDIALS ARKode &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; - typename VectorMemory::Pointer src_r(mem); - solver.reinit_vector(*src_r); - typename VectorMemory::Pointer dst_z(mem); - solver.reinit_vector(*dst_z); - copy(*src_r, r); - - int status = - solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr); - copy(z, *dst_z); - return status; + auto *src_r = internal::unwrap_nvector_const(r); + auto *dst_z = internal::unwrap_nvector(z); + + return solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr); } @@ -506,23 +425,17 @@ namespace SUNDIALS auto content = access_content(LS); ARKode &solver = *(content->solver); - GrowingVectorMemory mem; - typename VectorMemory::Pointer rhs(mem); - solver.reinit_vector(*rhs); - typename VectorMemory::Pointer dst(mem); - solver.reinit_vector(*dst); + auto *src_b = internal::unwrap_nvector_const(b); + auto *dst_x = internal::unwrap_nvector(x); - copy(*rhs, b); - - SundialsOperator op(solver, + SundialsOperator op(solver, content->A_data, content->a_times_fn); + SundialsPreconditioner preconditioner( solver, content->P_data, content->preconditioner_solve, tol); - int err = content->lsolve(op, preconditioner, *dst, *rhs, tol); - copy(x, *dst); - return err; + return content->lsolve(op, preconditioner, *dst_x, *src_b, tol); } @@ -630,8 +543,6 @@ namespace SUNDIALS const MPI_Comm & mpi_comm) : data(data) , arkode_mem(nullptr) - , yy(nullptr) - , abs_tolls(nullptr) , communicator(is_serial_vector::value ? MPI_COMM_SELF : Utilities::MPI::duplicate_communicator(mpi_comm)) @@ -665,30 +576,6 @@ namespace SUNDIALS unsigned int ARKode::solve_ode(VectorType &solution) { - const unsigned int system_size = 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 = solution.locally_owned_elements(); - const std::size_t local_system_size = is.n_elements(); - - yy = 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); - abs_tolls = N_VNew_Serial(system_size); - } - DiscreteTime time(data.initial_time, data.final_time, data.initial_step_size); @@ -711,8 +598,6 @@ namespace SUNDIALS (void)status; AssertARKode(status); - copy(solution, yy); - time.set_next_step_size(actual_next_time - time.get_current_time()); time.advance_time(); @@ -727,70 +612,31 @@ namespace SUNDIALS time.get_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(abs_tolls); - } - else -# endif - { - N_VDestroy_Serial(yy); - N_VDestroy_Serial(abs_tolls); - } - return time.get_step_number(); } # if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) template void - ARKode::reset(const double current_time, - const double current_time_step, - const VectorType &solution) + ARKode::reset(const double current_time, + const double current_time_step, + VectorType & solution) { - unsigned int system_size; - if (arkode_mem) ARKodeFree(&arkode_mem); arkode_mem = ARKodeCreate(); - // Free the vectors which are no longer used. - if (yy) - { - free_vector(yy); - free_vector(abs_tolls); - } - 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); - - abs_tolls = - N_VNew_Parallel(communicator, local_system_size, system_size); - } - else -# endif - { - yy = N_VNew_Serial(system_size); - abs_tolls = N_VNew_Serial(system_size); - } - - copy(yy, solution); Assert(explicit_function || implicit_function, ExcFunctionNotProvided("explicit_function || implicit_function")); + // just a view on the memory in solution, all write operations on yy by + // ARKODE will automatically be mirrored to solution + yy = internal::make_nvector_view(solution); + status = ARKodeInit( arkode_mem, explicit_function ? &t_arkode_explicit_function : nullptr, @@ -801,7 +647,7 @@ namespace SUNDIALS if (get_local_tolerances) { - copy(abs_tolls, get_local_tolerances()); + abs_tolls = make_nvector_view(get_local_tolerances()); status = ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls); AssertARKode(status); @@ -885,23 +731,17 @@ namespace SUNDIALS void ARKode::reset(const double current_time, const double current_time_step, - const VectorType &solution) + VectorType &solution) { if (arkode_mem) ARKStepFree(&arkode_mem); - // Free the vectors which are no longer used. - if (yy) - { - free_vector(yy); - free_vector(abs_tolls); - } - int status; (void)status; - yy = create_vector(solution); - abs_tolls = create_vector(solution); - copy(yy, solution); + + // just a view on the memory in solution, all write operations on yy by + // ARKODE will automatically be mirrored to solution + yy = internal::make_nvector_view(solution); Assert(explicit_function || implicit_function, ExcFunctionNotProvided("explicit_function || implicit_function")); @@ -912,9 +752,11 @@ namespace SUNDIALS current_time, yy); + Assert(arkode_mem != nullptr, ExcInternalError()); + if (get_local_tolerances) { - copy(abs_tolls, get_local_tolerances()); + abs_tolls = internal::make_nvector_view(get_local_tolerances()); status = ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls); AssertARKode(status); @@ -1009,13 +851,12 @@ namespace SUNDIALS } else { - N_Vector y_template = create_vector(solution); + auto y_template = internal::make_nvector_view(solution); SUNNonlinearSolver fixed_point_solver = SUNNonlinSol_FixedPoint(y_template, data.anderson_acceleration_subspace); - free_vector(y_template); status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver); AssertARKode(status); } @@ -1098,43 +939,6 @@ namespace SUNDIALS - template - N_Vector - ARKode::create_vector(const VectorType &template_vector) const - { - N_Vector new_vector; - auto system_size = template_vector.size(); - -# ifdef DEAL_II_WITH_MPI - if (is_serial_vector::value == false) - { - new_vector = - N_VNew_Parallel(communicator, - template_vector.locally_owned_elements().n_elements(), - system_size); - } - else -# endif - { - new_vector = N_VNew_Serial(system_size); - } - return new_vector; - } - - template - void - ARKode::free_vector(N_Vector vector) const - { -# ifdef DEAL_II_WITH_MPI - if (is_serial_vector::value == false) - N_VDestroy_Parallel(vector); - else -# endif - N_VDestroy_Serial(vector); - } - - - template void * ARKode::get_arkode_memory() const @@ -1163,17 +967,11 @@ namespace SUNDIALS SundialsOperator::vmult(VectorType & dst, const VectorType &src) const { - N_Vector sun_dst = solver.create_vector(dst); - N_Vector sun_src = solver.create_vector(src); - copy(sun_src, src); - int status = a_times_fn(A_data, sun_src, sun_dst); + auto sun_dst = internal::make_nvector_view(dst); + auto sun_src = internal::make_nvector_view(src); + int status = a_times_fn(A_data, sun_src, sun_dst); (void)status; AssertARKode(status); - - copy(dst, sun_dst); - - solver.free_vector(sun_dst); - solver.free_vector(sun_src); } @@ -1204,20 +1002,14 @@ namespace SUNDIALS return; } - N_Vector sun_dst = solver.create_vector(dst); - N_Vector sun_src = solver.create_vector(src); - copy(sun_src, src); + auto sun_dst = internal::make_nvector_view(dst); + auto sun_src = internal::make_nvector_view(src); // for custom preconditioners no distinction between left and right // preconditioning is made int status = p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/); (void)status; AssertARKode(status); - - copy(dst, sun_dst); - - solver.free_vector(sun_dst); - solver.free_vector(sun_src); } # endif diff --git a/tests/sundials/arkode_03.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_03.with_sundials.geq.5.4.0.output index 9e2ed6e8a6..503f3a999d 100644 --- a/tests/sundials/arkode_03.with_sundials.geq.5.4.0.output +++ b/tests/sundials/arkode_03.with_sundials.geq.5.4.0.output @@ -1,12 +1,12 @@ -DEAL::0 3.9 1.1 2.8 -DEAL::1 2.075274578 1.059601971 2.499947283 -DEAL::2 1.100460018 1.723610111 2.499972388 -DEAL::3 0.7896202372 2.326668043 2.499979776 -DEAL::4 0.8117726433 2.725583965 2.499987084 -DEAL::5 1.187002916 2.595541878 2.499973736 -DEAL::6 1.887311324 1.494098057 2.499952698 -DEAL::7 1.323271341 1.619016372 2.499965777 -DEAL::8 0.9217874641 2.130219386 2.499970637 -DEAL::9 0.8495425771 2.539475069 2.499991366 -DEAL::10 1.064969831 2.595955794 2.499973383 \ No newline at end of file +DEAL::0.00000 3.900000000 1.100000000 2.800000000 +DEAL::1.000000000 2.075274578 1.059601971 2.499947291 +DEAL::2.000000000 1.100460018 1.723610111 2.499972248 +DEAL::3.000000000 0.7896202372 2.326668043 2.499980243 +DEAL::4.000000000 0.8117726433 2.725583965 2.499982784 +DEAL::5.000000000 1.187002916 2.595541878 2.499973234 +DEAL::6.000000000 1.887311324 1.494098057 2.499934582 +DEAL::7.000000000 1.323271340 1.619016372 2.499963670 +DEAL::8.000000000 0.9217874641 2.130219386 2.499972986 +DEAL::9.000000000 0.8495425770 2.539475069 2.499983024 +DEAL::10.00000000 1.064969831 2.595955794 2.499973386