114 double t =
data.initial_time;
115 double h =
data.initial_step_size;
116 unsigned int step_number = 0;
121 reset(
data.initial_time,
data.initial_step_size, solution, solution_dot);
127 double next_time =
data.initial_time;
129 output_step(0, solution, solution_dot, 0);
131 while (t <
data.final_time)
133 next_time +=
data.output_period;
135 auto yy = internal::make_nvector_view(solution
141 auto yp = internal::make_nvector_view(solution_dot
152 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
153 if (pending_exception)
157 std::rethrow_exception(pending_exception);
161 pending_exception =
nullptr;
169 pending_exception =
nullptr;
175 status = IDAGetLastStep(ida_mem, &h);
178 while (solver_should_restart(t, solution, solution_dot))
179 reset(t, h, solution, solution_dot);
183 output_step(t, solution, solution_dot, step_number);
186 status = IDAGetNumSteps(ida_mem, &n_steps);
196 const double current_time_step,
197 VectorType &solution,
198 VectorType &solution_dot)
200 bool first_step = (current_time ==
data.initial_time);
205# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
206 status = SUNContext_Free(&ida_ctx);
211 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? SUN_COMM_NULL :
215# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
216 status = SUNContext_Free(&ida_ctx);
221 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
233# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
234 ida_mem = IDACreate();
236 ida_mem = IDACreate(ida_ctx);
239 auto yy = internal::make_nvector_view(solution
245 auto yp = internal::make_nvector_view(solution_dot
258 void *user_data) ->
int {
261 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
262 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
263 auto *residual = internal::unwrap_nvector<VectorType>(rr);
277 if (get_local_tolerances)
279 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
285 status = IDASVtolerances(ida_mem,
data.relative_tolerance, abs_tols);
290 status = IDASStolerances(ida_mem,
291 data.relative_tolerance,
292 data.absolute_tolerance);
296 status = IDASetInitStep(ida_mem, current_time_step);
299 status = IDASetUserData(ida_mem,
this);
302 if (
data.ic_type == AdditionalData::use_y_diff ||
303 data.reset_type == AdditionalData::use_y_diff ||
304 data.ignore_algebraic_terms_for_errors)
306 VectorType diff_comp_vector(solution);
307 diff_comp_vector = 0.0;
308 for (
const auto &component : differential_components())
309 diff_comp_vector[component] = 1.0;
312 const auto diff_id = internal::make_nvector_view(diff_comp_vector
318 status = IDASetId(ida_mem, diff_id);
322 status = IDASetSuppressAlg(ida_mem,
data.ignore_algebraic_terms_for_errors);
326 status = IDASetStopTime(ida_mem,
data.final_time);
329 status = IDASetMaxNonlinIters(ida_mem,
data.maximum_non_linear_iterations);
333 SUNMatrix J =
nullptr;
340# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
341 LS = SUNLinSolNewEmpty();
343 LS = SUNLinSolNewEmpty(ida_ctx);
349 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
355 LS->content =
nullptr;
376 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
377 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
405# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
406 J = SUNMatNewEmpty();
408 J = SUNMatNewEmpty(ida_ctx);
412 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
413 return SUNMATRIX_CUSTOM;
416 J->ops->destroy = [](SUNMatrix A) {
419 A->content =
nullptr;
431 status = IDASetLinearSolver(ida_mem, LS, J);
434 status = IDASetLSNormFactor(ida_mem,
data.ls_norm_factor);
439 status = IDASetJacFn(
454 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
455 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
466 status = IDASetMaxOrd(ida_mem,
data.maximum_order);
473 type =
data.reset_type;
476 IDASetMaxNumItersIC(ida_mem,
data.maximum_non_linear_iterations_ic);
479 if (type == AdditionalData::use_y_dot)
483 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
486 status = IDAGetConsistentIC(ida_mem, yy, yp);
489 else if (type == AdditionalData::use_y_diff)
492 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
495 status = IDAGetConsistentIC(ida_mem, yy, yp);