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