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