22#ifdef DEAL_II_WITH_SUNDIALS
29# ifdef DEAL_II_WITH_TRILINOS
33# ifdef DEAL_II_WITH_PETSC
49 template <
typename VectorType>
51 residual_callback(realtype tt,
57 IDA<VectorType> &solver = *
static_cast<IDA<VectorType> *
>(user_data);
59 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
60 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
61 auto *residual = internal::unwrap_nvector<VectorType>(rr);
63 int err = solver.residual(tt, *src_yy, *src_yp, *residual);
70 template <
typename VectorType>
72 setup_jacobian_callback(realtype tt,
84 IDA<VectorType> &solver = *
static_cast<IDA<VectorType> *
>(user_data);
86 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
87 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
89 int err = solver.setup_jacobian(tt, *src_yy, *src_yp, cj);
96 template <
typename VectorType>
104 IDA<VectorType> &solver = *
static_cast<IDA<VectorType> *
>(LS->content);
106 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
107 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
109 if (solver.solve_with_jacobian)
110 err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
111 else if (solver.solve_jacobian_system)
112 err = solver.solve_jacobian_system(*src_b, *dst_x);
123 template <
typename VectorType>
125 :
IDA(data, MPI_COMM_SELF)
130 template <
typename VectorType>
137 , mpi_communicator(mpi_comm)
144# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
157 template <
typename VectorType>
161# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
162 const int status = SUNContext_Free(&ida_ctx);
170 template <
typename VectorType>
174 double t = data.initial_time;
175 double h = data.initial_step_size;
176 unsigned int step_number = 0;
181 reset(data.initial_time, data.initial_step_size, solution, solution_dot);
187 double next_time = data.initial_time;
189 output_step(0, solution, solution_dot, 0);
191 while (t < data.final_time)
193 next_time += data.output_period;
195 auto yy = internal::make_nvector_view(solution
201 auto yp = internal::make_nvector_view(solution_dot
208 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
211 status = IDAGetLastStep(ida_mem, &h);
214 while (solver_should_restart(t, solution, solution_dot))
215 reset(t, h, solution, solution_dot);
219 output_step(t, solution, solution_dot, step_number);
227 template <
typename VectorType>
230 const double current_time_step,
231 VectorType & solution,
232 VectorType & solution_dot)
234 bool first_step = (current_time == data.initial_time);
239# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
240 status = SUNContext_Free(&ida_ctx);
245 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
257# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
258 ida_mem = IDACreate();
260 ida_mem = IDACreate(ida_ctx);
263 auto yy = internal::make_nvector_view(solution
269 auto yp = internal::make_nvector_view(solution_dot
277 IDAInit(ida_mem, residual_callback<VectorType>, current_time, yy, yp);
279 if (get_local_tolerances)
281 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
287 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
292 status = IDASStolerances(ida_mem,
293 data.relative_tolerance,
294 data.absolute_tolerance);
298 status = IDASetInitStep(ida_mem, current_time_step);
301 status = IDASetUserData(ida_mem,
this);
304 if (data.ic_type == AdditionalData::use_y_diff ||
305 data.reset_type == AdditionalData::use_y_diff ||
306 data.ignore_algebraic_terms_for_errors)
308 VectorType diff_comp_vector(solution);
309 diff_comp_vector = 0.0;
310 auto dc = differential_components();
311 for (
auto i = dc.begin(); i != dc.end(); ++i)
312 diff_comp_vector[*i] = 1.0;
315 const auto diff_id = internal::make_nvector_view(diff_comp_vector
321 status = IDASetId(ida_mem, diff_id);
325 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
329 status = IDASetStopTime(ida_mem, data.final_time);
332 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
336 SUNMatrix J =
nullptr;
343# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
344 LS = SUNLinSolNewEmpty();
346 LS = SUNLinSolNewEmpty(ida_ctx);
352 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
358 LS->content =
nullptr;
370 AssertThrow(solve_jacobian_system || solve_with_jacobian,
371 ExcFunctionNotProvided(
372 "solve_jacobian_system or solve_with_jacobian"));
373 LS->ops->solve = solve_with_jacobian_callback<VectorType>;
394# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
395 J = SUNMatNewEmpty();
397 J = SUNMatNewEmpty(ida_ctx);
401 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
402 return SUNMATRIX_CUSTOM;
405 J->ops->destroy = [](SUNMatrix A) {
408 A->content =
nullptr;
420 status = IDASetLinearSolver(ida_mem, LS, J);
423 status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
428 status = IDASetJacFn(ida_mem, &setup_jacobian_callback<VectorType>);
430 status = IDASetMaxOrd(ida_mem, data.maximum_order);
437 type = data.reset_type;
440 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
443 if (type == AdditionalData::use_y_dot)
447 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
450 status = IDAGetConsistentIC(ida_mem, yy, yp);
453 else if (type == AdditionalData::use_y_diff)
456 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
459 status = IDAGetConsistentIC(ida_mem, yy, yp);
464 template <
typename VectorType>
468 reinit_vector = [](VectorType &) {
469 AssertThrow(
false, ExcFunctionNotProvided(
"reinit_vector"));
472 residual = [](
const double,
475 VectorType &) ->
int {
477 AssertThrow(
false, ExcFunctionNotProvided(
"residual"));
482 output_step = [](
const double,
485 const unsigned int) {
return; };
487 solver_should_restart =
488 [](
const double, VectorType &, VectorType &) ->
bool {
return false; };
490 differential_components = [&]() ->
IndexSet {
494 return v->locally_owned_elements();
501# ifdef DEAL_II_WITH_MPI
503# ifdef DEAL_II_WITH_TRILINOS
508# ifdef DEAL_II_WITH_PETSC
509# ifndef PETSC_USE_COMPLEX
InitialConditionCorrection
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
MPI_Comm mpi_communicator
void set_functions_to_trigger_an_assert()
void reset(const double t, const double h, VectorType &y, VectorType &yp)
IDA(const AdditionalData &data=AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
#define DEAL_II_SUNDIALS_VERSION_LT(major, minor, patch)
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)