Reference documentation for deal.II version Git 191d06ed00 2021-05-11 21:15:49 -0400
\(\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\}}\)
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 
19 #include <deal.II/sundials/ida.h>
20 
21 #ifdef DEAL_II_WITH_SUNDIALS
22 
23 # if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
24 
25 # include <deal.II/base/utilities.h>
26 
28 # ifdef DEAL_II_WITH_TRILINOS
31 # endif
32 # ifdef DEAL_II_WITH_PETSC
35 # endif
36 
37 # include <deal.II/sundials/copy.h>
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 
50 namespace 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
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1528
void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src)
Definition: copy.cc:136
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1013
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:160
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
size_type n_elements() const
Definition: index_set.h:1832
static ::ExceptionBase & ExcInternalError()