Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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 - 2020 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
17#include <deal.II/base/config.h>
18
20
21#ifdef DEAL_II_WITH_SUNDIALS
22
23# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
24
26
28# ifdef DEAL_II_WITH_TRILINOS
31# endif
32# ifdef DEAL_II_WITH_PETSC
35# endif
36
38
39# ifdef DEAL_II_SUNDIALS_WITH_IDAS
40# include <idas/idas_impl.h>
41# else
42# include <ida/ida_impl.h>
43# endif
44
45# include <iomanip>
46# include <iostream>
47
49
50namespace SUNDIALS
51{
52 using namespace internal;
53
54 namespace
55 {
56 template <typename VectorType>
57 int
58 t_dae_residual(realtype tt,
59 N_Vector yy,
60 N_Vector yp,
61 N_Vector rr,
62 void * user_data)
63 {
64 IDA<VectorType> &solver = *static_cast<IDA<VectorType> *>(user_data);
66
67 typename VectorMemory<VectorType>::Pointer src_yy(mem);
68 solver.reinit_vector(*src_yy);
69
70 typename VectorMemory<VectorType>::Pointer src_yp(mem);
71 solver.reinit_vector(*src_yp);
72
73 typename VectorMemory<VectorType>::Pointer residual(mem);
74 solver.reinit_vector(*residual);
75
76 copy(*src_yy, yy);
77 copy(*src_yp, yp);
78
79 int err = solver.residual(tt, *src_yy, *src_yp, *residual);
80
81 copy(rr, *residual);
82
83 return err;
84 }
85
86
87
88 template <typename VectorType>
89 int
90 t_dae_lsetup(IDAMem IDA_mem,
91 N_Vector yy,
92 N_Vector yp,
93 N_Vector resp,
94 N_Vector tmp1,
95 N_Vector tmp2,
96 N_Vector tmp3)
97 {
98 (void)tmp1;
99 (void)tmp2;
100 (void)tmp3;
101 (void)resp;
102 IDA<VectorType> &solver =
103 *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
105
106 typename VectorMemory<VectorType>::Pointer src_yy(mem);
107 solver.reinit_vector(*src_yy);
108
109 typename VectorMemory<VectorType>::Pointer src_yp(mem);
110 solver.reinit_vector(*src_yp);
111
112 copy(*src_yy, yy);
113 copy(*src_yp, yp);
114
115 int err = solver.setup_jacobian(IDA_mem->ida_tn,
116 *src_yy,
117 *src_yp,
118 IDA_mem->ida_cj);
119
120 return err;
121 }
122
123
124 template <typename VectorType>
125 int
126 t_dae_solve(IDAMem IDA_mem,
127 N_Vector b,
128 N_Vector weight,
129 N_Vector yy,
130 N_Vector yp,
131 N_Vector resp)
132 {
133 (void)weight;
134 (void)yy;
135 (void)yp;
136 (void)resp;
137 IDA<VectorType> &solver =
138 *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
140
141 typename VectorMemory<VectorType>::Pointer src(mem);
142 solver.reinit_vector(*src);
143
144 typename VectorMemory<VectorType>::Pointer dst(mem);
145 solver.reinit_vector(*dst);
146
147 copy(*src, b);
148
149 int err = solver.solve_jacobian_system(*src, *dst);
150 copy(b, *dst);
151
152 return err;
153 }
154
155 } // namespace
156
157 template <typename VectorType>
158 IDA<VectorType>::IDA(const AdditionalData &data, const MPI_Comm &mpi_comm)
159 : data(data)
160 , ida_mem(nullptr)
161 , yy(nullptr)
162 , yp(nullptr)
163 , abs_tolls(nullptr)
164 , diff_id(nullptr)
165 , communicator(is_serial_vector<VectorType>::value ?
166 MPI_COMM_SELF :
167 Utilities::MPI::duplicate_communicator(mpi_comm))
168 {
169 set_functions_to_trigger_an_assert();
170 }
171
172 template <typename VectorType>
173 IDA<VectorType>::~IDA()
174 {
175 if (ida_mem)
176 IDAFree(&ida_mem);
177# ifdef DEAL_II_WITH_MPI
179 {
180 const int ierr = MPI_Comm_free(&communicator);
181 (void)ierr;
182 AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
183 }
184# endif
185 }
186
187
188
189 template <typename VectorType>
190 unsigned int
191 IDA<VectorType>::solve_dae(VectorType &solution, VectorType &solution_dot)
192 {
193 unsigned int system_size = solution.size();
194
195 double t = data.initial_time;
196 double h = data.initial_step_size;
197 unsigned int step_number = 0;
198
199 int status;
200 (void)status;
201
202 // The solution is stored in
203 // solution. Here we take only a
204 // view of it.
205# ifdef DEAL_II_WITH_MPI
207 {
208 const IndexSet is = solution.locally_owned_elements();
209 const std::size_t local_system_size = is.n_elements();
210
211 yy = N_VNew_Parallel(communicator, local_system_size, system_size);
212
213 yp = N_VNew_Parallel(communicator, local_system_size, system_size);
214
215 diff_id = N_VNew_Parallel(communicator, local_system_size, system_size);
216
217 abs_tolls =
218 N_VNew_Parallel(communicator, local_system_size, system_size);
219 }
220 else
221# endif
222 {
225 "Trying to use a serial code with a parallel vector."));
226 yy = N_VNew_Serial(system_size);
227 yp = N_VNew_Serial(system_size);
228 diff_id = N_VNew_Serial(system_size);
229 abs_tolls = N_VNew_Serial(system_size);
230 }
231 reset(data.initial_time, data.initial_step_size, solution, solution_dot);
232
233 double next_time = data.initial_time;
234
235 output_step(0, solution, solution_dot, 0);
236
237 while (t < data.final_time)
238 {
239 next_time += data.output_period;
240
241 status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL);
242 AssertIDA(status);
243
244 status = IDAGetLastStep(ida_mem, &h);
245 AssertIDA(status);
246
247 copy(solution, yy);
248 copy(solution_dot, yp);
249
250 while (solver_should_restart(t, solution, solution_dot))
251 reset(t, h, solution, solution_dot);
252
253 step_number++;
254
255 output_step(t, solution, solution_dot, step_number);
256 }
257
258 // Free the vectors which are no longer used.
259# ifdef DEAL_II_WITH_MPI
261 {
262 N_VDestroy_Parallel(yy);
263 N_VDestroy_Parallel(yp);
264 N_VDestroy_Parallel(abs_tolls);
265 N_VDestroy_Parallel(diff_id);
266 }
267 else
268# endif
269 {
270 N_VDestroy_Serial(yy);
271 N_VDestroy_Serial(yp);
272 N_VDestroy_Serial(abs_tolls);
273 N_VDestroy_Serial(diff_id);
274 }
275
276 return step_number;
277 }
278
279 template <typename VectorType>
280 void
281 IDA<VectorType>::reset(const double current_time,
282 const double current_time_step,
283 VectorType & solution,
284 VectorType & solution_dot)
285 {
286 unsigned int system_size;
287 bool first_step = (current_time == data.initial_time);
288
289 if (ida_mem)
290 IDAFree(&ida_mem);
291
292 ida_mem = IDACreate();
293
294
295 // Free the vectors which are no longer used.
296 if (yy)
297 {
298# ifdef DEAL_II_WITH_MPI
300 {
301 N_VDestroy_Parallel(yy);
302 N_VDestroy_Parallel(yp);
303 N_VDestroy_Parallel(abs_tolls);
304 N_VDestroy_Parallel(diff_id);
305 }
306 else
307# endif
308 {
309 N_VDestroy_Serial(yy);
310 N_VDestroy_Serial(yp);
311 N_VDestroy_Serial(abs_tolls);
312 N_VDestroy_Serial(diff_id);
313 }
314 }
315
316 int status;
317 (void)status;
318 system_size = solution.size();
319# ifdef DEAL_II_WITH_MPI
321 {
322 const IndexSet is = solution.locally_owned_elements();
323 const std::size_t local_system_size = is.n_elements();
324
325 yy = N_VNew_Parallel(communicator, local_system_size, system_size);
326
327 yp = N_VNew_Parallel(communicator, local_system_size, system_size);
328
329 diff_id = N_VNew_Parallel(communicator, local_system_size, system_size);
330
331 abs_tolls =
332 N_VNew_Parallel(communicator, local_system_size, system_size);
333 }
334 else
335# endif
336 {
337 yy = N_VNew_Serial(system_size);
338 yp = N_VNew_Serial(system_size);
339 diff_id = N_VNew_Serial(system_size);
340 abs_tolls = N_VNew_Serial(system_size);
341 }
342
343 copy(yy, solution);
344 copy(yp, solution_dot);
345
346 status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
347 AssertIDA(status);
348
349 if (get_local_tolerances)
350 {
351 copy(abs_tolls, get_local_tolerances());
352 status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tolls);
353 AssertIDA(status);
354 }
355 else
356 {
357 status = IDASStolerances(ida_mem,
358 data.relative_tolerance,
359 data.absolute_tolerance);
360 AssertIDA(status);
361 }
362
363 status = IDASetInitStep(ida_mem, current_time_step);
364 AssertIDA(status);
365
366 status = IDASetUserData(ida_mem, this);
367 AssertIDA(status);
368
369 if (data.ic_type == AdditionalData::use_y_diff ||
370 data.reset_type == AdditionalData::use_y_diff ||
371 data.ignore_algebraic_terms_for_errors)
372 {
373 VectorType diff_comp_vector(solution);
374 diff_comp_vector = 0.0;
375 auto dc = differential_components();
376 for (auto i = dc.begin(); i != dc.end(); ++i)
377 diff_comp_vector[*i] = 1.0;
378
379 copy(diff_id, diff_comp_vector);
380 status = IDASetId(ida_mem, diff_id);
381 AssertIDA(status);
382 }
383
384 status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors);
385 AssertIDA(status);
386
387 // status = IDASetMaxNumSteps(ida_mem, max_steps);
388 status = IDASetStopTime(ida_mem, data.final_time);
389 AssertIDA(status);
390
391 status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations);
392 AssertIDA(status);
393
394 // Initialize solver
395 auto IDA_mem = static_cast<IDAMem>(ida_mem);
396
397 IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
398 IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
399# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
400 IDA_mem->ida_setupNonNull = true;
401# endif
402
403 status = IDASetMaxOrd(ida_mem, data.maximum_order);
404 AssertIDA(status);
405
406 typename AdditionalData::InitialConditionCorrection type;
407 if (first_step)
408 type = data.ic_type;
409 else
410 type = data.reset_type;
411
412 status =
413 IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic);
414 AssertIDA(status);
415
416 if (type == AdditionalData::use_y_dot)
417 {
418 // (re)initialization of the vectors
419 status =
420 IDACalcIC(ida_mem, IDA_Y_INIT, current_time + current_time_step);
421 AssertIDA(status);
422
423 status = IDAGetConsistentIC(ida_mem, yy, yp);
424 AssertIDA(status);
425
426 copy(solution, yy);
427 copy(solution_dot, yp);
428 }
429 else if (type == AdditionalData::use_y_diff)
430 {
431 status =
432 IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time + current_time_step);
433 AssertIDA(status);
434
435 status = IDAGetConsistentIC(ida_mem, yy, yp);
436 AssertIDA(status);
437
438 copy(solution, yy);
439 copy(solution_dot, yp);
440 }
441 }
442
443 template <typename VectorType>
444 void
445 IDA<VectorType>::set_functions_to_trigger_an_assert()
446 {
447 reinit_vector = [](VectorType &) {
448 AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
449 };
450
451 residual = [](const double,
452 const VectorType &,
453 const VectorType &,
454 VectorType &) -> int {
455 int ret = 0;
456 AssertThrow(false, ExcFunctionNotProvided("residual"));
457 return ret;
458 };
459
460 setup_jacobian = [](const double,
461 const VectorType &,
462 const VectorType &,
463 const double) -> int {
464 int ret = 0;
465 AssertThrow(false, ExcFunctionNotProvided("setup_jacobian"));
466 return ret;
467 };
468
469 solve_jacobian_system = [](const VectorType &, VectorType &) -> int {
470 int ret = 0;
471 AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system"));
472 return ret;
473 };
474
475 output_step = [](const double,
476 const VectorType &,
477 const VectorType &,
478 const unsigned int) { return; };
479
480 solver_should_restart =
481 [](const double, VectorType &, VectorType &) -> bool { return false; };
482
483 differential_components = [&]() -> IndexSet {
485 typename VectorMemory<VectorType>::Pointer v(mem);
486 reinit_vector(*v);
487 const unsigned int size = v->size();
488 return complete_index_set(size);
489 };
490 }
491
492 template class IDA<Vector<double>>;
493 template class IDA<BlockVector<double>>;
494
495# ifdef DEAL_II_WITH_MPI
496
497# ifdef DEAL_II_WITH_TRILINOS
498 template class IDA<TrilinosWrappers::MPI::Vector>;
499 template class IDA<TrilinosWrappers::MPI::BlockVector>;
500# endif // DEAL_II_WITH_TRILINOS
501
502# ifdef DEAL_II_WITH_PETSC
503# ifndef PETSC_USE_COMPLEX
504 template class IDA<PETScWrappers::MPI::Vector>;
505 template class IDA<PETScWrappers::MPI::BlockVector>;
506# endif // PETSC_USE_COMPLEX
507# endif // DEAL_II_WITH_PETSC
508
509# endif // DEAL_II_WITH_MPI
510
511} // namespace SUNDIALS
512
514
515# endif
516
517#endif // DEAL_II_WITH_SUNDIALS
size_type n_elements() const
Definition: index_set.h:1832
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1013
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1528
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:160
void copy(const T *begin, const T *end, U *dest)