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