Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
ida.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2015 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#include <deal.II/base/config.h>
17
19
21
22#ifdef DEAL_II_WITH_SUNDIALS
24
26
29# ifdef DEAL_II_WITH_TRILINOS
32# endif
33# ifdef DEAL_II_WITH_PETSC
36# endif
37
40
41# include <iomanip>
42# include <iostream>
43
45
46namespace SUNDIALS
47{
48 template <typename VectorType>
50 : IDA(data, MPI_COMM_SELF)
51 {}
52
53
54
55 template <typename VectorType>
56 IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm mpi_comm)
57 : data(data)
58 , ida_mem(nullptr)
60 , ida_ctx(nullptr)
61# endif
62 , mpi_communicator(mpi_comm)
63 , pending_exception(nullptr)
64 {
65 // SUNDIALS will always duplicate communicators if we provide them. This
66 // can cause problems if SUNDIALS is configured with MPI and we pass along
67 // MPI_COMM_SELF in a serial application as MPI won't be
68 // initialized. Hence, work around that by just not providing a
69 // communicator in that case.
70# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
71 const int status =
72 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
74 &ida_ctx);
75 (void)status;
76 AssertIDA(status);
77# endif
79 }
80
81
82
83 template <typename VectorType>
85 {
86 IDAFree(&ida_mem);
87# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
88 const int status = SUNContext_Free(&ida_ctx);
89 (void)status;
90 AssertIDA(status);
91# endif
92
93 Assert(pending_exception == nullptr, ExcInternalError());
94 }
95
96
97
98 template <typename VectorType>
99 unsigned int
100 IDA<VectorType>::solve_dae(VectorType &solution, VectorType &solution_dot)
101 {
102 double t = data.initial_time;
103 double h = data.initial_step_size;
104 unsigned int step_number = 0;
105
106 int status;
107 (void)status;
108
109 reset(data.initial_time, data.initial_step_size, solution, solution_dot);
110
111 // The solution is stored in
112 // solution. Here we take only a
113 // view of it.
114
115 double next_time = data.initial_time;
116
117 output_step(0, solution, solution_dot, 0);
118
119 while (t < data.final_time)
120 {
121 next_time += data.output_period;
122
123 auto yy = internal::make_nvector_view(solution
124# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
125 ,
126 ida_ctx
127# endif
128 );
129 auto yp = internal::make_nvector_view(solution_dot
130# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
131 ,
132 ida_ctx
133# endif
134 );
135
136 // Execute time steps. If we ended up with a pending exception,
137 // see if it was recoverable; if it was, and if IDA recovered,
138 // just continue on. If IDA did not recover, rethrow the exception.
139 // Do the same if the exception was not recoverable.
140 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
141 if (pending_exception)
142 {
143 try
144 {
145 std::rethrow_exception(pending_exception);
146 }
147 catch (const RecoverableUserCallbackError &exc)
148 {
149 pending_exception = nullptr;
150 if (status == 0)
151 /* just eat the exception and continue */;
152 else
153 throw;
154 }
155 catch (...)
156 {
157 pending_exception = nullptr;
158 throw;
159 }
160 }
161 AssertIDA(status);
162
163 status = IDAGetLastStep(ida_mem, &h);
164 AssertIDA(status);
165
166 while (solver_should_restart(t, solution, solution_dot))
167 reset(t, h, solution, solution_dot);
168
169 step_number++;
170
171 output_step(t, solution, solution_dot, step_number);
172 }
173
174 return step_number;
175 }
176
177
178
179 template <typename VectorType>
180 void
181 IDA<VectorType>::reset(const double current_time,
182 const double current_time_step,
183 VectorType & solution,
184 VectorType & solution_dot)
185 {
186 bool first_step = (current_time == data.initial_time);
187
188 int status;
189 (void)status;
190
191# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
192 status = SUNContext_Free(&ida_ctx);
193 AssertIDA(status);
194
195 // Same comment applies as in class constructor:
196 status =
197 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
198 &mpi_communicator,
199 &ida_ctx);
200 AssertIDA(status);
201# endif
202
203 if (ida_mem)
204 {
205 IDAFree(&ida_mem);
206 // Initialization is version-dependent: do that in a moment
207 }
208
209# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
210 ida_mem = IDACreate();
211# else
212 ida_mem = IDACreate(ida_ctx);
213# endif
214
215 auto yy = internal::make_nvector_view(solution
216# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
217 ,
218 ida_ctx
219# endif
220 );
221 auto yp = internal::make_nvector_view(solution_dot
222# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
223 ,
224 ida_ctx
225# endif
226 );
227
228 status = IDAInit(
229 ida_mem,
230 [](realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
231 -> int {
232 IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
233
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);
237
239 solver.residual,
240 solver.pending_exception,
241 tt,
242 *src_yy,
243 *src_yp,
244 *residual);
245 },
246 current_time,
247 yy,
248 yp);
249 AssertIDA(status);
250 if (get_local_tolerances)
251 {
252 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
253# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
254 ,
255 ida_ctx
256# endif
257 );
258 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
259 AssertIDA(status);
260 }
261 else
262 {
263 status = IDASStolerances(ida_mem,
264 data.relative_tolerance,
265 data.absolute_tolerance);
266 AssertIDA(status);
267 }
268
269 status = IDASetInitStep(ida_mem, current_time_step);
270 AssertIDA(status);
271
272 status = IDASetUserData(ida_mem, this);
273 AssertIDA(status);
274
275 if (data.ic_type == AdditionalData::use_y_diff ||
276 data.reset_type == AdditionalData::use_y_diff ||
277 data.ignore_algebraic_terms_for_errors)
278 {
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;
283 diff_comp_vector.compress(VectorOperation::insert);
284
285 const auto diff_id = internal::make_nvector_view(diff_comp_vector
286# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
287 ,
288 ida_ctx
289# endif
290 );
291 status = IDASetId(ida_mem, diff_id);
292 AssertIDA(status);
293 }
294
295 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
296 AssertIDA(status);
297
298 // status = IDASetMaxNumSteps(ida_mem, max_steps);
299 status = IDASetStopTime(ida_mem, data.final_time);
300 AssertIDA(status);
301
302 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
303 AssertIDA(status);
304
305 // Initialize solver
306 SUNMatrix J = nullptr;
307 SUNLinearSolver LS = nullptr;
308
309 // and attach it to the SUNLinSol object. The functions that will get
310 // called do not actually receive the IDAMEM object, just the LS
311 // object, so we have to store a pointer to the current
312 // object in the LS object
313# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
314 LS = SUNLinSolNewEmpty();
315# else
316 LS = SUNLinSolNewEmpty(ida_ctx);
317# endif
318
319 LS->content = this;
320
321 LS->ops->gettype = [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
322 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
323 };
324
325 LS->ops->free = [](SUNLinearSolver LS) -> int {
326 if (LS->content)
327 {
328 LS->content = nullptr;
329 }
330 if (LS->ops)
331 {
332 free(LS->ops);
333 LS->ops = nullptr;
334 }
335 free(LS);
336 LS = nullptr;
337 return 0;
338 };
339
340 AssertThrow(solve_jacobian_system || solve_with_jacobian,
342 "solve_jacobian_system or solve_with_jacobian"));
343 LS->ops->solve = [](SUNLinearSolver LS,
344 SUNMatrix /*ignored*/,
345 N_Vector x,
346 N_Vector b,
347 realtype tol) -> int {
348 IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
349
350 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
351 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
352 if (solver.solve_with_jacobian)
354 solver.solve_with_jacobian,
355 solver.pending_exception,
356 *src_b,
357 *dst_x,
358 tol);
359 else if (solver.solve_jacobian_system)
362 solver.pending_exception,
363 *src_b,
364 *dst_x);
365 else
366 {
367 // We have already checked this outside, so we should never get here.
368 Assert(false, ExcInternalError());
369 return -1;
370 }
371 };
372
373 // When we set an iterative solver IDA requires that resid is provided. From
374 // SUNDIALS docs If an iterative method computes the preconditioned initial
375 // residual and returns with a successful solve without performing any
376 // iterations (i.e., either the initial guess or the preconditioner is
377 // sufficiently accurate), then this optional routine may be called by the
378 // SUNDIALS package. This routine should return the N_Vector containing the
379 // preconditioned initial residual vector.
380 LS->ops->resid = [](SUNLinearSolver /*ignored*/) -> N_Vector {
381 return nullptr;
382 };
383 // When we set an iterative solver IDA requires that last number of
384 // iteration is provided. Since we can't know what kind of solver the user
385 // has provided we set 1. This is clearly suboptimal.
386 LS->ops->numiters = [](SUNLinearSolver /*ignored*/) -> int { return 1; };
387 // Even though we don't use it, IDA still wants us to set some
388 // kind of matrix object for the nonlinear solver. This is because
389 // if we don't set it, it won't call the functions that set up
390 // the matrix object (i.e., the argument to the 'IDASetJacFn'
391 // function below).
392# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
393 J = SUNMatNewEmpty();
394# else
395 J = SUNMatNewEmpty(ida_ctx);
396# endif
397 J->content = this;
398
399 J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
400 return SUNMATRIX_CUSTOM;
401 };
402
403 J->ops->destroy = [](SUNMatrix A) {
404 if (A->content)
405 {
406 A->content = nullptr;
407 }
408 if (A->ops)
409 {
410 free(A->ops);
411 A->ops = nullptr;
412 }
413 free(A);
414 A = nullptr;
415 };
416
417 // Now set the linear system and Jacobian objects in the solver:
418 status = IDASetLinearSolver(ida_mem, LS, J);
419 AssertIDA(status);
420
421 status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
422 AssertIDA(status);
423 // Finally tell IDA about
424 // it as well. The manual says that this must happen *after*
425 // calling IDASetLinearSolver
426 status = IDASetJacFn(
427 ida_mem,
428 [](realtype tt,
429 realtype cj,
430 N_Vector yy,
431 N_Vector yp,
432 N_Vector /* residual */,
433 SUNMatrix /* ignored */,
434 void *user_data,
435 N_Vector /* tmp1 */,
436 N_Vector /* tmp2 */,
437 N_Vector /* tmp3 */) -> int {
438 Assert(user_data != nullptr, ExcInternalError());
439 IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
440
441 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
442 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
443
445 solver.setup_jacobian,
446 solver.pending_exception,
447 tt,
448 *src_yy,
449 *src_yp,
450 cj);
451 });
452 AssertIDA(status);
453 status = IDASetMaxOrd(ida_mem, data.maximum_order);
454 AssertIDA(status);
455
457 if (first_step)
458 type = data.ic_type;
459 else
460 type = data.reset_type;
461
462 status =
463 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
464 AssertIDA(status);
465
466 if (type == AdditionalData::use_y_dot)
467 {
468 // (re)initialization of the vectors
469 status =
470 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
471 AssertIDA(status);
472
473 status = IDAGetConsistentIC(ida_mem, yy, yp);
474 AssertIDA(status);
475 }
476 else if (type == AdditionalData::use_y_diff)
477 {
478 status =
479 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
480 AssertIDA(status);
481
482 status = IDAGetConsistentIC(ida_mem, yy, yp);
483 AssertIDA(status);
484 }
485 }
486
487 template <typename VectorType>
488 void
490 {
491 reinit_vector = [](VectorType &) {
492 AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
493 };
494
495 residual = [](const double,
496 const VectorType &,
497 const VectorType &,
498 VectorType &) -> int {
499 int ret = 0;
500 AssertThrow(false, ExcFunctionNotProvided("residual"));
501 return ret;
502 };
503
504
505 output_step = [](const double,
506 const VectorType &,
507 const VectorType &,
508 const unsigned int) { return; };
509
510 solver_should_restart =
511 [](const double, VectorType &, VectorType &) -> bool { return false; };
512
513 differential_components = [&]() -> IndexSet {
515 typename VectorMemory<VectorType>::Pointer v(mem);
516 reinit_vector(*v);
517 return v->locally_owned_elements();
518 };
519 }
520
521 template class IDA<Vector<double>>;
522 template class IDA<BlockVector<double>>;
523
524# ifdef DEAL_II_WITH_MPI
525
526# ifdef DEAL_II_WITH_TRILINOS
529# endif // DEAL_II_WITH_TRILINOS
530
531# ifdef DEAL_II_WITH_PETSC
532# ifndef PETSC_USE_COMPLEX
533 template class IDA<PETScWrappers::MPI::Vector>;
535# endif // PETSC_USE_COMPLEX
536# endif // DEAL_II_WITH_PETSC
537
538# endif // DEAL_II_WITH_MPI
539
540} // namespace SUNDIALS
541
543
544#endif // DEAL_II_WITH_SUNDIALS
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
Definition ida.h:781
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition ida.cc:100
MPI_Comm mpi_communicator
Definition ida.h:902
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
Definition ida.h:701
void set_functions_to_trigger_an_assert()
Definition ida.cc:489
void reset(const double t, const double h, VectorType &y, VectorType &yp)
Definition ida.cc:181
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
Definition ida.h:662
IDA(const AdditionalData &data=AdditionalData())
Definition ida.cc:49
std::function< void(const VectorType &rhs, VectorType &dst)> solve_jacobian_system
Definition ida.h:736
SUNContext ida_ctx
Definition ida.h:895
std::exception_ptr pending_exception
Definition ida.h:914
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition config.h:359
#define DEAL_II_SUNDIALS_VERSION_LT(major, minor, patch)
Definition config.h:366
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & RecoverableUserCallbackError()
#define AssertThrow(cond, exc)
#define AssertIDA(code)
Definition ida.h:57
int call_and_possibly_capture_exception(const F &f, std::exception_ptr &eptr, Args &&...args)
Definition utilities.h:46