21#ifdef DEAL_II_WITH_SUNDIALS
28# ifdef DEAL_II_WITH_TRILINOS
32# ifdef DEAL_II_WITH_PETSC
47 template <
typename VectorType>
54 template <
typename VectorType>
71# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
78# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
90 template <
typename VectorType>
94# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
95 const int status = SUNContext_Free(&
ida_ctx);
105 template <
typename VectorType>
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;
126 while (t <
data.final_time)
128 next_time +=
data.output_period;
147 status = IDASolve(
ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
170 status = IDAGetLastStep(
ida_mem, &h);
174 reset(t, h, solution, solution_dot);
178 output_step(t, solution, solution_dot, step_number);
181 status = IDAGetNumSteps(
ida_mem, &n_steps);
188 template <
typename VectorType>
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);
210# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
211 status = SUNContext_Free(&
ida_ctx);
228# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
253 void *user_data) ->
int {
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);
299 data.ignore_algebraic_terms_for_errors)
301 VectorType diff_comp_vector(solution);
302 diff_comp_vector = 0.0;
304 diff_comp_vector[component] = 1.0;
313 status = IDASetId(
ida_mem, diff_id);
317 status = IDASetSuppressAlg(
ida_mem,
data.ignore_algebraic_terms_for_errors);
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;
400# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
401 J = SUNMatNewEmpty();
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(
461 status = IDASetMaxOrd(
ida_mem,
data.maximum_order);
468 type =
data.reset_type;
471 IDASetMaxNumItersIC(
ida_mem,
data.maximum_non_linear_iterations_ic);
478 IDACalcIC(
ida_mem, IDA_Y_INIT, current_time + current_time_step);
481 status = IDAGetConsistentIC(
ida_mem, yy, yp);
487 IDACalcIC(
ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
490 status = IDAGetConsistentIC(
ida_mem, yy, yp);
495 template <
typename VectorType>
506 VectorType &) ->
int {
516 const unsigned int) {
return; };
519 [](
const double, VectorType &, VectorType &) ->
bool {
return false; };
525 return v->locally_owned_elements();
532# ifdef DEAL_II_WITH_MPI
534# ifdef DEAL_II_WITH_TRILINOS
539# ifdef DEAL_II_WITH_PETSC
540# ifndef PETSC_USE_COMPLEX
InitialConditionCorrection
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
MPI_Comm mpi_communicator
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
void set_functions_to_trigger_an_assert()
const AdditionalData data
void reset(const double t, const double h, VectorType &y, VectorType &yp)
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
std::function< void(const double t, const VectorType &sol, const VectorType &sol_dot, const unsigned int step_number)> output_step
std::function< IndexSet()> differential_components
std::function< VectorType &()> get_local_tolerances
std::function< bool(const double t, VectorType &sol, VectorType &sol_dot)> solver_should_restart
GrowingVectorMemory< VectorType > mem
std::function< void(VectorType &)> reinit_vector
IDA(const AdditionalData &data=AdditionalData())
std::exception_ptr pending_exception
NVectorView< VectorType > make_nvector_view(VectorType &vector, SUNContext nvector_context)
#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
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & RecoverableUserCallbackError()
#define AssertThrow(cond, exc)
int call_and_possibly_capture_exception(const F &f, std::exception_ptr &eptr, Args &&...args)
VectorType * unwrap_nvector(N_Vector v)
const VectorType * unwrap_nvector_const(N_Vector v)