21#ifdef DEAL_II_WITH_SUNDIALS
23# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
28# ifdef DEAL_II_WITH_TRILINOS
32# ifdef DEAL_II_WITH_PETSC
39# ifdef DEAL_II_SUNDIALS_WITH_IDAS
40# include <idas/idas_impl.h>
42# include <ida/ida_impl.h>
56 template <
typename VectorType>
58 t_dae_residual(realtype tt,
64 IDA<VectorType> &solver = *
static_cast<IDA<VectorType> *
>(user_data);
68 solver.reinit_vector(*src_yy);
71 solver.reinit_vector(*src_yp);
74 solver.reinit_vector(*residual);
79 int err = solver.residual(tt, *src_yy, *src_yp, *residual);
88 template <
typename VectorType>
90 t_dae_lsetup(IDAMem IDA_mem,
102 IDA<VectorType> &solver =
103 *
static_cast<IDA<VectorType> *
>(IDA_mem->ida_user_data);
107 solver.reinit_vector(*src_yy);
110 solver.reinit_vector(*src_yp);
115 int err = solver.setup_jacobian(IDA_mem->ida_tn,
124 template <
typename VectorType>
126 t_dae_solve(IDAMem IDA_mem,
137 IDA<VectorType> &solver =
138 *
static_cast<IDA<VectorType> *
>(IDA_mem->ida_user_data);
142 solver.reinit_vector(*src);
145 solver.reinit_vector(*dst);
149 int err = solver.solve_jacobian_system(*src, *dst);
157 template <
typename VectorType>
158 IDA<VectorType>::IDA(
const AdditionalData &data,
const MPI_Comm &mpi_comm)
169 set_functions_to_trigger_an_assert();
172 template <
typename VectorType>
173 IDA<VectorType>::~IDA()
177# ifdef DEAL_II_WITH_MPI
180 const int ierr = MPI_Comm_free(&communicator);
189 template <
typename VectorType>
191 IDA<VectorType>::solve_dae(VectorType &solution, VectorType &solution_dot)
193 unsigned int system_size = solution.size();
195 double t = data.initial_time;
196 double h = data.initial_step_size;
197 unsigned int step_number = 0;
205# ifdef DEAL_II_WITH_MPI
208 const IndexSet is = solution.locally_owned_elements();
209 const std::size_t local_system_size = is.
n_elements();
211 yy = N_VNew_Parallel(communicator, local_system_size, system_size);
213 yp = N_VNew_Parallel(communicator, local_system_size, system_size);
215 diff_id = N_VNew_Parallel(communicator, local_system_size, system_size);
218 N_VNew_Parallel(communicator, local_system_size, system_size);
225 "Trying to use a serial code with a parallel vector."));
226 yy = N_VNew_Serial(system_size);
227 yp = N_VNew_Serial(system_size);
228 diff_id = N_VNew_Serial(system_size);
229 abs_tolls = N_VNew_Serial(system_size);
231 reset(data.initial_time, data.initial_step_size, solution, solution_dot);
233 double next_time = data.initial_time;
235 output_step(0, solution, solution_dot, 0);
237 while (t < data.final_time)
239 next_time += data.output_period;
241 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
244 status = IDAGetLastStep(ida_mem, &h);
248 copy(solution_dot, yp);
250 while (solver_should_restart(t, solution, solution_dot))
251 reset(t, h, solution, solution_dot);
255 output_step(t, solution, solution_dot, step_number);
259# ifdef DEAL_II_WITH_MPI
262 N_VDestroy_Parallel(yy);
263 N_VDestroy_Parallel(yp);
264 N_VDestroy_Parallel(abs_tolls);
265 N_VDestroy_Parallel(diff_id);
270 N_VDestroy_Serial(yy);
271 N_VDestroy_Serial(yp);
272 N_VDestroy_Serial(abs_tolls);
273 N_VDestroy_Serial(diff_id);
279 template <
typename VectorType>
281 IDA<VectorType>::reset(
const double current_time,
282 const double current_time_step,
283 VectorType & solution,
284 VectorType & solution_dot)
286 unsigned int system_size;
287 bool first_step = (current_time == data.initial_time);
292 ida_mem = IDACreate();
298# ifdef DEAL_II_WITH_MPI
301 N_VDestroy_Parallel(yy);
302 N_VDestroy_Parallel(yp);
303 N_VDestroy_Parallel(abs_tolls);
304 N_VDestroy_Parallel(diff_id);
309 N_VDestroy_Serial(yy);
310 N_VDestroy_Serial(yp);
311 N_VDestroy_Serial(abs_tolls);
312 N_VDestroy_Serial(diff_id);
318 system_size = solution.size();
319# ifdef DEAL_II_WITH_MPI
322 const IndexSet is = solution.locally_owned_elements();
323 const std::size_t local_system_size = is.
n_elements();
325 yy = N_VNew_Parallel(communicator, local_system_size, system_size);
327 yp = N_VNew_Parallel(communicator, local_system_size, system_size);
329 diff_id = N_VNew_Parallel(communicator, local_system_size, system_size);
332 N_VNew_Parallel(communicator, local_system_size, system_size);
337 yy = N_VNew_Serial(system_size);
338 yp = N_VNew_Serial(system_size);
339 diff_id = N_VNew_Serial(system_size);
340 abs_tolls = N_VNew_Serial(system_size);
344 copy(yp, solution_dot);
346 status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
349 if (get_local_tolerances)
351 copy(abs_tolls, get_local_tolerances());
352 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tolls);
357 status = IDASStolerances(ida_mem,
358 data.relative_tolerance,
359 data.absolute_tolerance);
363 status = IDASetInitStep(ida_mem, current_time_step);
366 status = IDASetUserData(ida_mem,
this);
369 if (data.ic_type == AdditionalData::use_y_diff ||
370 data.reset_type == AdditionalData::use_y_diff ||
371 data.ignore_algebraic_terms_for_errors)
373 VectorType diff_comp_vector(solution);
374 diff_comp_vector = 0.0;
375 auto dc = differential_components();
376 for (
auto i = dc.begin(); i != dc.end(); ++i)
377 diff_comp_vector[*i] = 1.0;
379 copy(diff_id, diff_comp_vector);
380 status = IDASetId(ida_mem, diff_id);
384 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
388 status = IDASetStopTime(ida_mem, data.final_time);
391 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
395 auto IDA_mem =
static_cast<IDAMem
>(ida_mem);
397 IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
398 IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
399# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
400 IDA_mem->ida_setupNonNull =
true;
403 status = IDASetMaxOrd(ida_mem, data.maximum_order);
406 typename AdditionalData::InitialConditionCorrection type;
410 type = data.reset_type;
413 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
416 if (type == AdditionalData::use_y_dot)
420 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
423 status = IDAGetConsistentIC(ida_mem, yy, yp);
427 copy(solution_dot, yp);
429 else if (type == AdditionalData::use_y_diff)
432 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
435 status = IDAGetConsistentIC(ida_mem, yy, yp);
439 copy(solution_dot, yp);
443 template <
typename VectorType>
445 IDA<VectorType>::set_functions_to_trigger_an_assert()
447 reinit_vector = [](VectorType &) {
448 AssertThrow(
false, ExcFunctionNotProvided(
"reinit_vector"));
451 residual = [](
const double,
454 VectorType &) ->
int {
456 AssertThrow(
false, ExcFunctionNotProvided(
"residual"));
460 setup_jacobian = [](
const double,
463 const double) ->
int {
465 AssertThrow(
false, ExcFunctionNotProvided(
"setup_jacobian"));
469 solve_jacobian_system = [](
const VectorType &, VectorType &) ->
int {
471 AssertThrow(
false, ExcFunctionNotProvided(
"solve_jacobian_system"));
475 output_step = [](
const double,
478 const unsigned int) {
return; };
480 solver_should_restart =
481 [](
const double, VectorType &, VectorType &) ->
bool {
return false; };
483 differential_components = [&]() ->
IndexSet {
487 const unsigned int size = v->size();
492 template class IDA<Vector<double>>;
493 template class IDA<BlockVector<double>>;
495# ifdef DEAL_II_WITH_MPI
497# ifdef DEAL_II_WITH_TRILINOS
498 template class IDA<TrilinosWrappers::MPI::Vector>;
499 template class IDA<TrilinosWrappers::MPI::BlockVector>;
502# ifdef DEAL_II_WITH_PETSC
503# ifndef PETSC_USE_COMPLEX
504 template class IDA<PETScWrappers::MPI::Vector>;
505 template class IDA<PETScWrappers::MPI::BlockVector>;
size_type n_elements() const
IndexSet complete_index_set(const IndexSet::size_type N)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
void copy(const T *begin, const T *end, U *dest)