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