deal.II version GIT relicensing-1936-ge55124254b 2024-10-04 21: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\}}\)
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 long int n_steps;
181 status = IDAGetNumSteps(ida_mem, &n_steps);
182 AssertIDA(status);
183 return n_steps;
184 }
185
186
187
188 template <typename VectorType>
189 void
190 IDA<VectorType>::reset(const double current_time,
191 const double current_time_step,
192 VectorType &solution,
193 VectorType &solution_dot)
194 {
195 bool first_step = (current_time == data.initial_time);
196
197 int status;
198 (void)status;
199
200# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
201 status = SUNContext_Free(&ida_ctx);
202 AssertIDA(status);
203
204 // Same comment applies as in class constructor:
205 status =
206 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? SUN_COMM_NULL :
207 mpi_communicator,
208 &ida_ctx);
209 AssertIDA(status);
210# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
211 status = SUNContext_Free(&ida_ctx);
212 AssertIDA(status);
213
214 // Same comment applies as in class constructor:
215 status =
216 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
217 &mpi_communicator,
218 &ida_ctx);
219 AssertIDA(status);
220# endif
221
222 if (ida_mem)
223 {
224 IDAFree(&ida_mem);
225 // Initialization is version-dependent: do that in a moment
226 }
227
228# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
229 ida_mem = IDACreate();
230# else
231 ida_mem = IDACreate(ida_ctx);
232# endif
233
234 auto yy = internal::make_nvector_view(solution
235# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
236 ,
237 ida_ctx
238# endif
239 );
240 auto yp = internal::make_nvector_view(solution_dot
241# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
242 ,
243 ida_ctx
244# endif
245 );
246
247 status = IDAInit(
248 ida_mem,
249 [](SUNDIALS::realtype tt,
250 N_Vector yy,
251 N_Vector yp,
252 N_Vector rr,
253 void *user_data) -> int {
254 IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
255
256 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
257 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
258 auto *residual = internal::unwrap_nvector<VectorType>(rr);
259
261 solver.residual,
262 solver.pending_exception,
263 tt,
264 *src_yy,
265 *src_yp,
266 *residual);
267 },
268 current_time,
269 yy,
270 yp);
271 AssertIDA(status);
272 if (get_local_tolerances)
273 {
274 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
275# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
276 ,
277 ida_ctx
278# endif
279 );
280 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tols);
281 AssertIDA(status);
282 }
283 else
284 {
285 status = IDASStolerances(ida_mem,
286 data.relative_tolerance,
287 data.absolute_tolerance);
288 AssertIDA(status);
289 }
290
291 status = IDASetInitStep(ida_mem, current_time_step);
292 AssertIDA(status);
293
294 status = IDASetUserData(ida_mem, this);
295 AssertIDA(status);
296
297 if (data.ic_type == AdditionalData::use_y_diff ||
298 data.reset_type == AdditionalData::use_y_diff ||
299 data.ignore_algebraic_terms_for_errors)
300 {
301 VectorType diff_comp_vector(solution);
302 diff_comp_vector = 0.0;
303 for (const auto &component : differential_components())
304 diff_comp_vector[component] = 1.0;
305 diff_comp_vector.compress(VectorOperation::insert);
306
307 const auto diff_id = internal::make_nvector_view(diff_comp_vector
308# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
309 ,
310 ida_ctx
311# endif
312 );
313 status = IDASetId(ida_mem, diff_id);
314 AssertIDA(status);
315 }
316
317 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
318 AssertIDA(status);
319
320 // status = IDASetMaxNumSteps(ida_mem, max_steps);
321 status = IDASetStopTime(ida_mem, data.final_time);
322 AssertIDA(status);
323
324 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
325 AssertIDA(status);
326
327 // Initialize solver
328 SUNMatrix J = nullptr;
329 SUNLinearSolver LS = nullptr;
330
331 // and attach it to the SUNLinSol object. The functions that will get
332 // called do not actually receive the IDAMEM object, just the LS
333 // object, so we have to store a pointer to the current
334 // object in the LS object
335# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
336 LS = SUNLinSolNewEmpty();
337# else
338 LS = SUNLinSolNewEmpty(ida_ctx);
339# endif
340
341 LS->content = this;
342
343 LS->ops->gettype = [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
344 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
345 };
346
347 LS->ops->free = [](SUNLinearSolver LS) -> int {
348 if (LS->content)
349 {
350 LS->content = nullptr;
351 }
352 if (LS->ops)
353 {
354 free(LS->ops);
355 LS->ops = nullptr;
356 }
357 free(LS);
358 LS = nullptr;
359 return 0;
360 };
361
362 AssertThrow(solve_with_jacobian,
363 ExcFunctionNotProvided("solve_with_jacobian"));
364 LS->ops->solve = [](SUNLinearSolver LS,
365 SUNMatrix /*ignored*/,
366 N_Vector x,
367 N_Vector b,
368 SUNDIALS::realtype tol) -> int {
369 IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(LS->content);
370
371 auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
372 auto *dst_x = internal::unwrap_nvector<VectorType>(x);
374 solver.solve_with_jacobian,
375 solver.pending_exception,
376 *src_b,
377 *dst_x,
378 tol);
379 };
380
381 // When we set an iterative solver IDA requires that resid is provided. From
382 // SUNDIALS docs If an iterative method computes the preconditioned initial
383 // residual and returns with a successful solve without performing any
384 // iterations (i.e., either the initial guess or the preconditioner is
385 // sufficiently accurate), then this optional routine may be called by the
386 // SUNDIALS package. This routine should return the N_Vector containing the
387 // preconditioned initial residual vector.
388 LS->ops->resid = [](SUNLinearSolver /*ignored*/) -> N_Vector {
389 return nullptr;
390 };
391 // When we set an iterative solver IDA requires that last number of
392 // iteration is provided. Since we can't know what kind of solver the user
393 // has provided we set 1. This is clearly suboptimal.
394 LS->ops->numiters = [](SUNLinearSolver /*ignored*/) -> int { return 1; };
395 // Even though we don't use it, IDA still wants us to set some
396 // kind of matrix object for the nonlinear solver. This is because
397 // if we don't set it, it won't call the functions that set up
398 // the matrix object (i.e., the argument to the 'IDASetJacFn'
399 // function below).
400# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
401 J = SUNMatNewEmpty();
402# else
403 J = SUNMatNewEmpty(ida_ctx);
404# endif
405 J->content = this;
406
407 J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
408 return SUNMATRIX_CUSTOM;
409 };
410
411 J->ops->destroy = [](SUNMatrix A) {
412 if (A->content)
413 {
414 A->content = nullptr;
415 }
416 if (A->ops)
417 {
418 free(A->ops);
419 A->ops = nullptr;
420 }
421 free(A);
422 A = nullptr;
423 };
424
425 // Now set the linear system and Jacobian objects in the solver:
426 status = IDASetLinearSolver(ida_mem, LS, J);
427 AssertIDA(status);
428
429 status = IDASetLSNormFactor(ida_mem, data.ls_norm_factor);
430 AssertIDA(status);
431 // Finally tell IDA about
432 // it as well. The manual says that this must happen *after*
433 // calling IDASetLinearSolver
434 status = IDASetJacFn(
435 ida_mem,
436 [](SUNDIALS::realtype tt,
438 N_Vector yy,
439 N_Vector yp,
440 N_Vector /* residual */,
441 SUNMatrix /* ignored */,
442 void *user_data,
443 N_Vector /* tmp1 */,
444 N_Vector /* tmp2 */,
445 N_Vector /* tmp3 */) -> int {
446 Assert(user_data != nullptr, ExcInternalError());
447 IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
448
449 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
450 auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
451
453 solver.setup_jacobian,
454 solver.pending_exception,
455 tt,
456 *src_yy,
457 *src_yp,
458 cj);
459 });
460 AssertIDA(status);
461 status = IDASetMaxOrd(ida_mem, data.maximum_order);
462 AssertIDA(status);
463
465 if (first_step)
466 type = data.ic_type;
467 else
468 type = data.reset_type;
469
470 status =
471 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
472 AssertIDA(status);
473
474 if (type == AdditionalData::use_y_dot)
475 {
476 // (re)initialization of the vectors
477 status =
478 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
479 AssertIDA(status);
480
481 status = IDAGetConsistentIC(ida_mem, yy, yp);
482 AssertIDA(status);
483 }
484 else if (type == AdditionalData::use_y_diff)
485 {
486 status =
487 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
488 AssertIDA(status);
489
490 status = IDAGetConsistentIC(ida_mem, yy, yp);
491 AssertIDA(status);
492 }
493 }
494
495 template <typename VectorType>
496 void
498 {
499 reinit_vector = [](VectorType &) {
500 AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
501 };
502
503 residual = [](const double,
504 const VectorType &,
505 const VectorType &,
506 VectorType &) -> int {
507 int ret = 0;
508 AssertThrow(false, ExcFunctionNotProvided("residual"));
509 return ret;
510 };
511
512
513 output_step = [](const double,
514 const VectorType &,
515 const VectorType &,
516 const unsigned int) { return; };
517
518 solver_should_restart =
519 [](const double, VectorType &, VectorType &) -> bool { return false; };
520
521 differential_components = [&]() -> IndexSet {
523 typename VectorMemory<VectorType>::Pointer v(mem);
524 reinit_vector(*v);
525 return v->locally_owned_elements();
526 };
527 }
528
529 template class IDA<Vector<double>>;
530 template class IDA<BlockVector<double>>;
531
532# ifdef DEAL_II_WITH_MPI
533
534# ifdef DEAL_II_WITH_TRILINOS
537# endif // DEAL_II_WITH_TRILINOS
538
539# ifdef DEAL_II_WITH_PETSC
540# ifndef PETSC_USE_COMPLEX
541 template class IDA<PETScWrappers::MPI::Vector>;
543# endif // PETSC_USE_COMPLEX
544# endif // DEAL_II_WITH_PETSC
545
546# endif // DEAL_II_WITH_MPI
547
548} // namespace SUNDIALS
549
551
552#endif // DEAL_II_WITH_SUNDIALS
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
Definition ida.h:995
unsigned int solve_dae(VectorType &solution, VectorType &solution_dot)
Definition ida.cc:107
MPI_Comm mpi_communicator
Definition ida.h:1116
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, const double alpha)> setup_jacobian
Definition ida.h:950
void set_functions_to_trigger_an_assert()
Definition ida.cc:497
void reset(const double t, const double h, VectorType &y, VectorType &yp)
Definition ida.cc:190
std::function< void(const double t, const VectorType &y, const VectorType &y_dot, VectorType &res)> residual
Definition ida.h:914
IDA(const AdditionalData &data=AdditionalData())
Definition ida.cc:48
SUNContext ida_ctx
Definition ida.h:1109
std::exception_ptr pending_exception
Definition ida.h:1128
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition config.h:402
#define DEAL_II_SUNDIALS_VERSION_LT(major, minor, patch)
Definition config.h:409
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:55
int call_and_possibly_capture_exception(const F &f, std::exception_ptr &eptr, Args &&...args)
Definition utilities.h:45
::sunrealtype realtype