102 double t = data.initial_time;
103 double h = data.initial_step_size;
104 unsigned int step_number = 0;
109 reset(data.initial_time, data.initial_step_size, solution, solution_dot);
115 double next_time = data.initial_time;
117 output_step(0, solution, solution_dot, 0);
119 while (t < data.final_time)
121 next_time += data.output_period;
123 auto yy = internal::make_nvector_view(solution
129 auto yp = internal::make_nvector_view(solution_dot
140 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
141 if (pending_exception)
145 std::rethrow_exception(pending_exception);
149 pending_exception =
nullptr;
157 pending_exception =
nullptr;
163 status = IDAGetLastStep(ida_mem, &h);
166 while (solver_should_restart(t, solution, solution_dot))
167 reset(t, h, solution, solution_dot);
171 output_step(t, solution, solution_dot, step_number);
182 const double current_time_step,
183 VectorType & solution,
184 VectorType & solution_dot)
186 bool first_step = (current_time == data.initial_time);
191# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
192 status = SUNContext_Free(&ida_ctx);
197 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
209# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
210 ida_mem = IDACreate();
212 ida_mem = IDACreate(ida_ctx);
215 auto yy = internal::make_nvector_view(solution
221 auto yp = internal::make_nvector_view(solution_dot
230 [](realtype tt, N_Vector yy, N_Vector yp, N_Vector rr,
void *user_data)
234 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
235 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
236 auto *residual = internal::unwrap_nvector<VectorType>(rr);
250 if (get_local_tolerances)
252 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
258 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
263 status = IDASStolerances(ida_mem,
264 data.relative_tolerance,
265 data.absolute_tolerance);
269 status = IDASetInitStep(ida_mem, current_time_step);
272 status = IDASetUserData(ida_mem,
this);
275 if (data.ic_type == AdditionalData::use_y_diff ||
276 data.reset_type == AdditionalData::use_y_diff ||
277 data.ignore_algebraic_terms_for_errors)
279 VectorType diff_comp_vector(solution);
280 diff_comp_vector = 0.0;
281 for (
const auto &component : differential_components())
282 diff_comp_vector[component] = 1.0;
285 const auto diff_id = internal::make_nvector_view(diff_comp_vector
291 status = IDASetId(ida_mem, diff_id);
295 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
299 status = IDASetStopTime(ida_mem, data.final_time);
302 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
306 SUNMatrix J =
nullptr;
313# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
314 LS = SUNLinSolNewEmpty();
316 LS = SUNLinSolNewEmpty(ida_ctx);
322 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
328 LS->content =
nullptr;
340 AssertThrow(solve_jacobian_system || solve_with_jacobian,
342 "solve_jacobian_system or solve_with_jacobian"));
347 realtype tol) ->
int {
350 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
351 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
392# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
393 J = SUNMatNewEmpty();
395 J = SUNMatNewEmpty(ida_ctx);
399 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
400 return SUNMATRIX_CUSTOM;
403 J->ops->destroy = [](SUNMatrix A) {
406 A->content =
nullptr;
418 status = IDASetLinearSolver(ida_mem, LS, J);
421 status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
426 status = IDASetJacFn(
441 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
442 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
453 status = IDASetMaxOrd(ida_mem, data.maximum_order);
460 type = data.reset_type;
463 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
466 if (type == AdditionalData::use_y_dot)
470 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
473 status = IDAGetConsistentIC(ida_mem, yy, yp);
476 else if (type == AdditionalData::use_y_diff)
479 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
482 status = IDAGetConsistentIC(ida_mem, yy, yp);