101 double t = data.initial_time;
102 double h = data.initial_step_size;
103 unsigned int step_number = 0;
108 reset(data.initial_time, data.initial_step_size, solution, solution_dot);
114 double next_time = data.initial_time;
116 output_step(0, solution, solution_dot, 0);
118 while (t < data.final_time)
120 next_time += data.output_period;
122 auto yy = internal::make_nvector_view(solution
128 auto yp = internal::make_nvector_view(solution_dot
139 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
140 if (pending_exception)
144 std::rethrow_exception(pending_exception);
148 pending_exception =
nullptr;
156 pending_exception =
nullptr;
162 status = IDAGetLastStep(ida_mem, &h);
165 while (solver_should_restart(t, solution, solution_dot))
166 reset(t, h, solution, solution_dot);
170 output_step(t, solution, solution_dot, step_number);
181 const double current_time_step,
182 VectorType &solution,
183 VectorType &solution_dot)
185 bool first_step = (current_time == data.initial_time);
190# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
191 status = SUNContext_Free(&ida_ctx);
196 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
208# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
209 ida_mem = IDACreate();
211 ida_mem = IDACreate(ida_ctx);
214 auto yy = internal::make_nvector_view(solution
220 auto yp = internal::make_nvector_view(solution_dot
233 void *user_data) ->
int {
236 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
237 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
238 auto *residual = internal::unwrap_nvector<VectorType>(rr);
252 if (get_local_tolerances)
254 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
260 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
265 status = IDASStolerances(ida_mem,
266 data.relative_tolerance,
267 data.absolute_tolerance);
271 status = IDASetInitStep(ida_mem, current_time_step);
274 status = IDASetUserData(ida_mem,
this);
277 if (data.ic_type == AdditionalData::use_y_diff ||
278 data.reset_type == AdditionalData::use_y_diff ||
279 data.ignore_algebraic_terms_for_errors)
281 VectorType diff_comp_vector(solution);
282 diff_comp_vector = 0.0;
283 for (
const auto &component : differential_components())
284 diff_comp_vector[component] = 1.0;
287 const auto diff_id = internal::make_nvector_view(diff_comp_vector
293 status = IDASetId(ida_mem, diff_id);
297 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
301 status = IDASetStopTime(ida_mem, data.final_time);
304 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
308 SUNMatrix J =
nullptr;
315# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
316 LS = SUNLinSolNewEmpty();
318 LS = SUNLinSolNewEmpty(ida_ctx);
324 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
330 LS->content =
nullptr;
351 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
352 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
380# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
381 J = SUNMatNewEmpty();
383 J = SUNMatNewEmpty(ida_ctx);
387 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
388 return SUNMATRIX_CUSTOM;
391 J->ops->destroy = [](SUNMatrix A) {
394 A->content =
nullptr;
406 status = IDASetLinearSolver(ida_mem, LS, J);
409 status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
414 status = IDASetJacFn(
429 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
430 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
441 status = IDASetMaxOrd(ida_mem, data.maximum_order);
448 type = data.reset_type;
451 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
454 if (type == AdditionalData::use_y_dot)
458 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
461 status = IDAGetConsistentIC(ida_mem, yy, yp);
464 else if (type == AdditionalData::use_y_diff)
467 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
470 status = IDAGetConsistentIC(ida_mem, yy, yp);