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