236 const double current_time_step,
237 const VectorType &solution)
239 last_end_time = current_time;
243# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
244 status = SUNContext_Free(&arkode_ctx);
248 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
256 ARKStepFree(&arkode_mem);
262 auto initial_condition_nvector = internal::make_nvector_view(solution
269 Assert(explicit_function || implicit_function,
272 auto explicit_function_callback =
273 [](realtype tt, N_Vector yy, N_Vector yp,
void *user_data) ->
int {
278 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
279 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
290 auto implicit_function_callback =
291 [](realtype tt, N_Vector yy, N_Vector yp,
void *user_data) ->
int {
296 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
297 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
307 arkode_mem = ARKStepCreate(explicit_function ? explicit_function_callback :
309 implicit_function ? implicit_function_callback :
312 initial_condition_nvector
313# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
320 if (get_local_tolerances)
322 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
329 ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
334 status = ARKStepSStolerances(arkode_mem,
335 data.relative_tolerance,
336 data.absolute_tolerance);
341 status = ARKStepSetUserData(arkode_mem,
this);
344 setup_system_solver(solution);
346 setup_mass_solver(solution);
348 status = ARKStepSetInitStep(arkode_mem, current_time_step);
351 status = ARKStepSetStopTime(arkode_mem, data.final_time);
354 status = ARKStepSetOrder(arkode_mem, data.maximum_order);
358 custom_setup(arkode_mem);
368 if (!implicit_function)
375 if (jacobian_times_vector)
378 if (solve_linearized_system)
381 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
382 solve_linearized_system,
384# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
389 sun_linear_solver = *linear_solver;
395 auto y_template = internal::make_nvector_view(solution
401# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
403 SUNLinSol_SPGMR(y_template,
408 SUNLinSol_SPGMR(y_template,
414 status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver,
nullptr);
417 auto jacobian_times_vector_callback = [](N_Vector v,
428 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
429 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
430 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
432 auto *dst_Jv = internal::unwrap_nvector<VectorType>(Jv);
444 auto jacobian_times_vector_setup_callback =
445 [](realtype t, N_Vector y, N_Vector fy,
void *user_data) ->
int {
450 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
451 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
460 status = ARKStepSetJacTimes(arkode_mem,
461 jacobian_times_setup ?
462 jacobian_times_vector_setup_callback :
463 ARKLsJacTimesSetupFn(
nullptr),
464 jacobian_times_vector_callback);
466 if (jacobian_preconditioner_solve)
468 auto solve_with_jacobian_callback = [](realtype t,
476 void * user_data) ->
int {
481 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
482 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
483 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
485 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
500 auto jacobian_solver_setup_callback = [](realtype t,
504 booleantype *jcurPtr,
506 void *user_data) ->
int {
511 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
512 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
525 status = ARKStepSetPreconditioner(arkode_mem,
526 jacobian_preconditioner_setup ?
527 jacobian_solver_setup_callback :
528 ARKLsPrecSetupFn(
nullptr),
529 solve_with_jacobian_callback);
532 if (data.implicit_function_is_linear)
534 status = ARKStepSetLinear(
535 arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
541 auto y_template = internal::make_nvector_view(solution
548# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
549 SUNNonlinearSolver fixed_point_solver =
550 SUNNonlinSol_FixedPoint(y_template,
551 data.anderson_acceleration_subspace);
553 SUNNonlinearSolver fixed_point_solver =
554 SUNNonlinSol_FixedPoint(y_template,
555 data.anderson_acceleration_subspace,
559 status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
564 ARKStepSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
577 if (mass_times_vector)
583 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
587# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
592 sun_mass_linear_solver = *mass_solver;
596 auto y_template = internal::make_nvector_view(solution
602# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
603 sun_mass_linear_solver =
604 SUNLinSol_SPGMR(y_template,
608 sun_mass_linear_solver =
609 SUNLinSol_SPGMR(y_template,
615 booleantype mass_time_dependent =
616 data.mass_is_time_independent ? SUNFALSE : SUNTRUE;
617 status = ARKStepSetMassLinearSolver(arkode_mem,
618 sun_mass_linear_solver,
620 mass_time_dependent);
623 auto mass_matrix_times_vector_setup_callback =
624 [](realtype t,
void *mtimes_data) ->
int {
633 auto mass_matrix_times_vector_callback =
634 [](N_Vector v, N_Vector Mv, realtype t,
void *mtimes_data) ->
int {
639 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
640 auto *dst_Mv = internal::unwrap_nvector<VectorType>(Mv);
650 status = ARKStepSetMassTimes(arkode_mem,
652 mass_matrix_times_vector_setup_callback :
653 ARKLsMassTimesSetupFn(
nullptr),
654 mass_matrix_times_vector_callback,
658 if (mass_preconditioner_solve)
660 auto mass_matrix_solver_setup_callback =
661 [](realtype t,
void *user_data) ->
int {
670 auto solve_with_mass_matrix_callback = [](realtype t,
675 void *user_data) ->
int {
680 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
681 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
694 ARKStepSetMassPreconditioner(arkode_mem,
695 mass_preconditioner_setup ?
696 mass_matrix_solver_setup_callback :
697 ARKLsMassPrecSetupFn(
nullptr),
698 solve_with_mass_matrix_callback);
std::function< void(const double t, const VectorType &y, const VectorType &fy, const VectorType &r, VectorType &z, const double gamma, const double tol, const int lr) jacobian_preconditioner_solve)