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