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);
181 status = IDAGetNumSteps(ida_mem, &n_steps);
191 const double current_time_step,
192 VectorType &solution,
193 VectorType &solution_dot)
195 bool first_step = (current_time == data.initial_time);
200# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
201 status = SUNContext_Free(&ida_ctx);
206 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? SUN_COMM_NULL :
210# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
211 status = SUNContext_Free(&ida_ctx);
216 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
228# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
229 ida_mem = IDACreate();
231 ida_mem = IDACreate(ida_ctx);
234 auto yy = internal::make_nvector_view(solution
240 auto yp = internal::make_nvector_view(solution_dot
253 void *user_data) ->
int {
256 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
257 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
258 auto *residual = internal::unwrap_nvector<VectorType>(rr);
272 if (get_local_tolerances)
274 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
280 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
285 status = IDASStolerances(ida_mem,
286 data.relative_tolerance,
287 data.absolute_tolerance);
291 status = IDASetInitStep(ida_mem, current_time_step);
294 status = IDASetUserData(ida_mem,
this);
297 if (data.ic_type == AdditionalData::use_y_diff ||
298 data.reset_type == AdditionalData::use_y_diff ||
299 data.ignore_algebraic_terms_for_errors)
301 VectorType diff_comp_vector(solution);
302 diff_comp_vector = 0.0;
303 for (
const auto &component : differential_components())
304 diff_comp_vector[component] = 1.0;
307 const auto diff_id = internal::make_nvector_view(diff_comp_vector
313 status = IDASetId(ida_mem, diff_id);
317 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
321 status = IDASetStopTime(ida_mem, data.final_time);
324 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
328 SUNMatrix J =
nullptr;
335# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
336 LS = SUNLinSolNewEmpty();
338 LS = SUNLinSolNewEmpty(ida_ctx);
344 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
350 LS->content =
nullptr;
371 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
372 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
400# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
401 J = SUNMatNewEmpty();
403 J = SUNMatNewEmpty(ida_ctx);
407 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
408 return SUNMATRIX_CUSTOM;
411 J->ops->destroy = [](SUNMatrix A) {
414 A->content =
nullptr;
426 status = IDASetLinearSolver(ida_mem, LS, J);
429 status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
434 status = IDASetJacFn(
449 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
450 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
461 status = IDASetMaxOrd(ida_mem, data.maximum_order);
468 type = data.reset_type;
471 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
474 if (type == AdditionalData::use_y_dot)
478 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
481 status = IDAGetConsistentIC(ida_mem, yy, yp);
484 else if (type == AdditionalData::use_y_diff)
487 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
490 status = IDAGetConsistentIC(ida_mem, yy, yp);