109 double t = data.initial_time;
110 double h = data.initial_step_size;
111 unsigned int step_number = 0;
116 reset(data.initial_time, data.initial_step_size, solution, solution_dot);
122 double next_time = data.initial_time;
124 output_step(0, solution, solution_dot, 0);
126 while (t < data.final_time)
128 next_time += data.output_period;
130 auto yy = internal::make_nvector_view(solution
136 auto yp = internal::make_nvector_view(solution_dot
147 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
148 if (pending_exception)
152 std::rethrow_exception(pending_exception);
156 pending_exception =
nullptr;
164 pending_exception =
nullptr;
170 status = IDAGetLastStep(ida_mem, &h);
173 while (solver_should_restart(t, solution, solution_dot))
174 reset(t, h, solution, solution_dot);
178 output_step(t, solution, solution_dot, step_number);
189 const double current_time_step,
190 VectorType &solution,
191 VectorType &solution_dot)
193 bool first_step = (current_time == data.initial_time);
198# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
199 status = SUNContext_Free(&ida_ctx);
204 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? SUN_COMM_NULL :
208# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
209 status = SUNContext_Free(&ida_ctx);
214 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
226# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
227 ida_mem = IDACreate();
229 ida_mem = IDACreate(ida_ctx);
232 auto yy = internal::make_nvector_view(solution
238 auto yp = internal::make_nvector_view(solution_dot
251 void *user_data) ->
int {
254 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
255 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
256 auto *residual = internal::unwrap_nvector<VectorType>(rr);
270 if (get_local_tolerances)
272 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
278 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
283 status = IDASStolerances(ida_mem,
284 data.relative_tolerance,
285 data.absolute_tolerance);
289 status = IDASetInitStep(ida_mem, current_time_step);
292 status = IDASetUserData(ida_mem,
this);
295 if (data.ic_type == AdditionalData::use_y_diff ||
296 data.reset_type == AdditionalData::use_y_diff ||
297 data.ignore_algebraic_terms_for_errors)
299 VectorType diff_comp_vector(solution);
300 diff_comp_vector = 0.0;
301 for (
const auto &component : differential_components())
302 diff_comp_vector[component] = 1.0;
305 const auto diff_id = internal::make_nvector_view(diff_comp_vector
311 status = IDASetId(ida_mem, diff_id);
315 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
319 status = IDASetStopTime(ida_mem, data.final_time);
322 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
326 SUNMatrix J =
nullptr;
333# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
334 LS = SUNLinSolNewEmpty();
336 LS = SUNLinSolNewEmpty(ida_ctx);
342 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
348 LS->content =
nullptr;
369 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
370 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
398# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
399 J = SUNMatNewEmpty();
401 J = SUNMatNewEmpty(ida_ctx);
405 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
406 return SUNMATRIX_CUSTOM;
409 J->ops->destroy = [](SUNMatrix A) {
412 A->content =
nullptr;
424 status = IDASetLinearSolver(ida_mem, LS, J);
427 status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
432 status = IDASetJacFn(
447 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
448 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
459 status = IDASetMaxOrd(ida_mem, data.maximum_order);
466 type = data.reset_type;
469 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
472 if (type == AdditionalData::use_y_dot)
476 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
479 status = IDAGetConsistentIC(ida_mem, yy, yp);
482 else if (type == AdditionalData::use_y_diff)
485 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
488 status = IDAGetConsistentIC(ida_mem, yy, yp);