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\}}\)
arkode.cc
Go to the documentation of this file.
1 //-----------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 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/arkode.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 # include <arkode/arkode_impl.h>
38 # include <sundials/sundials_config.h>
39 
40 # include <iomanip>
41 # include <iostream>
42 
43 // Make sure we know how to call sundials own ARKode() function
44 const auto &SundialsARKode = ARKode;
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_arkode_explicit_function(realtype tt,
57  N_Vector yy,
58  N_Vector yp,
59  void * user_data)
60  {
61  ARKode<VectorType> &solver =
62  *static_cast<ARKode<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 dst_yp(mem);
69  solver.reinit_vector(*dst_yp);
70 
71  copy(*src_yy, yy);
72 
73  int err = solver.explicit_function(tt, *src_yy, *dst_yp);
74 
75  copy(yp, *dst_yp);
76 
77  return err;
78  }
79 
80 
81 
82  template <typename VectorType>
83  int
84  t_arkode_implicit_function(realtype tt,
85  N_Vector yy,
86  N_Vector yp,
87  void * user_data)
88  {
89  ARKode<VectorType> &solver =
90  *static_cast<ARKode<VectorType> *>(user_data);
91  GrowingVectorMemory<VectorType> mem;
92 
93  typename VectorMemory<VectorType>::Pointer src_yy(mem);
94  solver.reinit_vector(*src_yy);
95 
96  typename VectorMemory<VectorType>::Pointer dst_yp(mem);
97  solver.reinit_vector(*dst_yp);
98 
99  copy(*src_yy, yy);
100 
101  int err = solver.implicit_function(tt, *src_yy, *dst_yp);
102 
103  copy(yp, *dst_yp);
104 
105  return err;
106  }
107 
108 
109 
110  template <typename VectorType>
111  int
112  t_arkode_setup_jacobian(ARKodeMem arkode_mem,
113  int convfail,
114  N_Vector ypred,
115  N_Vector fpred,
116  booleantype *jcurPtr,
117  N_Vector,
118  N_Vector,
119  N_Vector)
120  {
121  ARKode<VectorType> &solver =
122  *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
123  GrowingVectorMemory<VectorType> mem;
124 
125  typename VectorMemory<VectorType>::Pointer src_ypred(mem);
126  solver.reinit_vector(*src_ypred);
127 
128  typename VectorMemory<VectorType>::Pointer src_fpred(mem);
129  solver.reinit_vector(*src_fpred);
130 
131  copy(*src_ypred, ypred);
132  copy(*src_fpred, fpred);
133 
134  // avoid reinterpret_cast
135  bool jcurPtr_tmp = false;
136  int err = solver.setup_jacobian(convfail,
137  arkode_mem->ark_tn,
138  arkode_mem->ark_gamma,
139  *src_ypred,
140  *src_fpred,
141  jcurPtr_tmp);
142 # if DEAL_II_SUNDIALS_VERSION_GTE(2, 0, 0)
143  *jcurPtr = jcurPtr_tmp ? SUNTRUE : SUNFALSE;
144 # else
145  *jcurPtr = jcurPtr_tmp ? TRUE : FALSE;
146 # endif
147 
148  return err;
149  }
150 
151 
152 
153  template <typename VectorType>
154  int
155  t_arkode_solve_jacobian(ARKodeMem arkode_mem,
156  N_Vector b,
157 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
158  N_Vector,
159 # endif
160  N_Vector ycur,
161  N_Vector fcur)
162  {
163  ARKode<VectorType> &solver =
164  *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
165  GrowingVectorMemory<VectorType> mem;
166 
167  typename VectorMemory<VectorType>::Pointer src(mem);
168  solver.reinit_vector(*src);
169 
170  typename VectorMemory<VectorType>::Pointer src_ycur(mem);
171  solver.reinit_vector(*src_ycur);
172 
173  typename VectorMemory<VectorType>::Pointer src_fcur(mem);
174  solver.reinit_vector(*src_fcur);
175 
176  typename VectorMemory<VectorType>::Pointer dst(mem);
177  solver.reinit_vector(*dst);
178 
179  copy(*src, b);
180  copy(*src_ycur, ycur);
181  copy(*src_fcur, fcur);
182 
183  int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
184  arkode_mem->ark_gamma,
185  *src_ycur,
186  *src_fcur,
187  *src,
188  *dst);
189  copy(b, *dst);
190 
191  return err;
192  }
193 
194 
195 
196  template <typename VectorType>
197  int
198  t_arkode_setup_mass(ARKodeMem arkode_mem, N_Vector, N_Vector, N_Vector)
199  {
200  ARKode<VectorType> &solver =
201  *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
202  int err = solver.setup_mass(arkode_mem->ark_tn);
203  return err;
204  }
205 
206 
207 
208  template <typename VectorType>
209  int
210  t_arkode_solve_mass(ARKodeMem arkode_mem,
211 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
212  N_Vector b,
213  N_Vector
214 # else
215  N_Vector b
216 # endif
217  )
218  {
219  ARKode<VectorType> &solver =
220  *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
221  GrowingVectorMemory<VectorType> mem;
222 
223  typename VectorMemory<VectorType>::Pointer src(mem);
224  solver.reinit_vector(*src);
225 
226  typename VectorMemory<VectorType>::Pointer dst(mem);
227  solver.reinit_vector(*dst);
228 
229  copy(*src, b);
230 
231  int err = solver.solve_mass_system(*src, *dst);
232  copy(b, *dst);
233 
234  return err;
235  }
236  } // namespace
237 
238  template <typename VectorType>
239  ARKode<VectorType>::ARKode(const AdditionalData &data,
240  const MPI_Comm mpi_comm)
241  : data(data)
242  , arkode_mem(nullptr)
243  , yy(nullptr)
244  , abs_tolls(nullptr)
245  , communicator(is_serial_vector<VectorType>::value ?
246  MPI_COMM_SELF :
247  Utilities::MPI::duplicate_communicator(mpi_comm))
248  {
249  set_functions_to_trigger_an_assert();
250  }
251 
252  template <typename VectorType>
253  ARKode<VectorType>::~ARKode()
254  {
255  if (arkode_mem)
256  ARKodeFree(&arkode_mem);
257 # ifdef DEAL_II_WITH_MPI
258  if (is_serial_vector<VectorType>::value == false)
259  {
260  const int ierr = MPI_Comm_free(&communicator);
261  (void)ierr;
262  AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
263  }
264 # endif
265  }
266 
267 
268 
269  template <typename VectorType>
270  unsigned int
271  ARKode<VectorType>::solve_ode(VectorType &solution)
272  {
273  unsigned int system_size = solution.size();
274 
275  double t = data.initial_time;
276  double h = data.initial_step_size;
277  unsigned int step_number = 0;
278 
279  int status;
280  (void)status;
281 
282  // The solution is stored in
283  // solution. Here we take only a
284  // view of it.
285 # ifdef DEAL_II_WITH_MPI
286  if (is_serial_vector<VectorType>::value == false)
287  {
288  const IndexSet is = solution.locally_owned_elements();
289  const std::size_t local_system_size = is.n_elements();
290 
291  yy = N_VNew_Parallel(communicator, local_system_size, system_size);
292 
293  abs_tolls =
294  N_VNew_Parallel(communicator, local_system_size, system_size);
295  }
296  else
297 # endif
298  {
299  Assert(is_serial_vector<VectorType>::value,
300  ExcInternalError(
301  "Trying to use a serial code with a parallel vector."));
302  yy = N_VNew_Serial(system_size);
303  abs_tolls = N_VNew_Serial(system_size);
304  }
305  reset(data.initial_time, data.initial_step_size, solution);
306 
307  double next_time = data.initial_time;
308 
309  if (output_step)
310  output_step(0, solution, 0);
311 
312  while (t < data.final_time)
313  {
314  next_time += data.output_period;
315 
316  status = SundialsARKode(arkode_mem, next_time, yy, &t, ARK_NORMAL);
317 
318  AssertARKode(status);
319 
320  status = ARKodeGetLastStep(arkode_mem, &h);
321  AssertARKode(status);
322 
323  copy(solution, yy);
324 
325  while (solver_should_restart(t, solution))
326  reset(t, h, solution);
327 
328  step_number++;
329 
330  if (output_step)
331  output_step(t, solution, step_number);
332  }
333 
334  // Free the vectors which are no longer used.
335 # ifdef DEAL_II_WITH_MPI
336  if (is_serial_vector<VectorType>::value == false)
337  {
338  N_VDestroy_Parallel(yy);
339  N_VDestroy_Parallel(abs_tolls);
340  }
341  else
342 # endif
343  {
344  N_VDestroy_Serial(yy);
345  N_VDestroy_Serial(abs_tolls);
346  }
347 
348  return step_number;
349  }
350 
351  template <typename VectorType>
352  void
353  ARKode<VectorType>::reset(const double current_time,
354  const double current_time_step,
355  const VectorType &solution)
356  {
357  unsigned int system_size;
358 
359  if (arkode_mem)
360  ARKodeFree(&arkode_mem);
361 
362  arkode_mem = ARKodeCreate();
363 
364  // Free the vectors which are no longer used.
365  if (yy)
366  {
367 # ifdef DEAL_II_WITH_MPI
368  if (is_serial_vector<VectorType>::value == false)
369  {
370  N_VDestroy_Parallel(yy);
371  N_VDestroy_Parallel(abs_tolls);
372  }
373  else
374 # endif
375  {
376  N_VDestroy_Serial(yy);
377  N_VDestroy_Serial(abs_tolls);
378  }
379  }
380 
381  int status;
382  (void)status;
383  system_size = solution.size();
384 # ifdef DEAL_II_WITH_MPI
385  if (is_serial_vector<VectorType>::value == false)
386  {
387  const IndexSet is = solution.locally_owned_elements();
388  const std::size_t local_system_size = is.n_elements();
389 
390  yy = N_VNew_Parallel(communicator, local_system_size, system_size);
391 
392  abs_tolls =
393  N_VNew_Parallel(communicator, local_system_size, system_size);
394  }
395  else
396 # endif
397  {
398  yy = N_VNew_Serial(system_size);
399  abs_tolls = N_VNew_Serial(system_size);
400  }
401 
402  copy(yy, solution);
403 
404  Assert(explicit_function || implicit_function,
405  ExcFunctionNotProvided("explicit_function || implicit_function"));
406 
407  status = ARKodeInit(
408  arkode_mem,
409  explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
410  implicit_function ? &t_arkode_implicit_function<VectorType> : nullptr,
411  current_time,
412  yy);
413  AssertARKode(status);
414 
415  if (get_local_tolerances)
416  {
417  copy(abs_tolls, get_local_tolerances());
418  status =
419  ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
420  AssertARKode(status);
421  }
422  else
423  {
424  status = ARKodeSStolerances(arkode_mem,
425  data.relative_tolerance,
426  data.absolute_tolerance);
427  AssertARKode(status);
428  }
429 
430  status = ARKodeSetInitStep(arkode_mem, current_time_step);
431  AssertARKode(status);
432 
433  status = ARKodeSetUserData(arkode_mem, this);
434  AssertARKode(status);
435 
436  status = ARKodeSetStopTime(arkode_mem, data.final_time);
437  AssertARKode(status);
438 
439  status =
440  ARKodeSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
441  AssertARKode(status);
442 
443  // Initialize solver
444  auto ARKode_mem = static_cast<ARKodeMem>(arkode_mem);
445 
446  if (solve_jacobian_system)
447  {
448  status = ARKodeSetNewton(arkode_mem);
449  AssertARKode(status);
450  if (data.implicit_function_is_linear)
451  {
452  status = ARKodeSetLinear(
453  arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
454  AssertARKode(status);
455  }
456 
457 
458  ARKode_mem->ark_lsolve = t_arkode_solve_jacobian<VectorType>;
459  if (setup_jacobian)
460  {
461  ARKode_mem->ark_lsetup = t_arkode_setup_jacobian<VectorType>;
462 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
463  ARKode_mem->ark_setupNonNull = true;
464 # endif
465  }
466  }
467  else
468  {
469  status =
470  ARKodeSetFixedPoint(arkode_mem, data.maximum_non_linear_iterations);
471  AssertARKode(status);
472  }
473 
474 
475  if (solve_mass_system)
476  {
477  ARKode_mem->ark_msolve = t_arkode_solve_mass<VectorType>;
478 
479  if (setup_mass)
480  {
481  ARKode_mem->ark_msetup = t_arkode_setup_mass<VectorType>;
482 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
483  ARKode_mem->ark_MassSetupNonNull = true;
484 # endif
485  }
486  }
487 
488  status = ARKodeSetOrder(arkode_mem, data.maximum_order);
489  AssertARKode(status);
490  }
491 
492  template <typename VectorType>
493  void
494  ARKode<VectorType>::set_functions_to_trigger_an_assert()
495  {
496  reinit_vector = [](VectorType &) {
497  AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
498  };
499 
500  solver_should_restart = [](const double, VectorType &) -> bool {
501  return false;
502  };
503  }
504 
505  template class ARKode<Vector<double>>;
506  template class ARKode<BlockVector<double>>;
507 
508 # ifdef DEAL_II_WITH_MPI
509 
510 # ifdef DEAL_II_WITH_TRILINOS
511  template class ARKode<TrilinosWrappers::MPI::Vector>;
512  template class ARKode<TrilinosWrappers::MPI::BlockVector>;
513 # endif // DEAL_II_WITH_TRILINOS
514 
515 # ifdef DEAL_II_WITH_PETSC
516 # ifndef PETSC_USE_COMPLEX
517  template class ARKode<PETScWrappers::MPI::Vector>;
518  template class ARKode<PETScWrappers::MPI::BlockVector>;
519 # endif // PETSC_USE_COMPLEX
520 # endif // DEAL_II_WITH_PETSC
521 
522 # endif // DEAL_II_WITH_MPI
523 
524 } // namespace SUNDIALS
525 
526 DEAL_II_NAMESPACE_CLOSE
527 
528 #endif
IndexSet
Definition: index_set.h:74
SundialsARKode
const auto & SundialsARKode
Definition: arkode.cc:44
SUNDIALS::ARKode
Definition: arkode.h:313
utilities.h
SUNDIALS::ARKode::reset
void reset(const double t, const double h, const VectorType &y)
Definition: arkode.cc:353
SUNDIALS::internal::copy
void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src)
Definition: copy.cc:65
AssertARKode
#define AssertARKode(code)
Definition: arkode.h:53
VectorType
MPI_Comm
SUNDIALS
Definition: arkode.h:58
copy.h
SUNDIALS::ARKode::~ARKode
~ARKode()
Definition: arkode.cc:253
GrowingVectorMemory
Definition: vector_memory.h:320
IndexSet::n_elements
size_type n_elements() const
Definition: index_set.h:1833
SUNDIALS::ARKode::solve_ode
unsigned int solve_ode(VectorType &solution)
Definition: arkode.cc:271
is_serial_vector
Definition: vector_type_traits.h:49
SUNDIALS::ARKode::AdditionalData
Definition: arkode.h:319
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)
petsc_block_vector.h
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
Utilities::MPI::duplicate_communicator
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:139
SUNDIALS::ARKode::ARKode
ARKode(const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
Definition: arkode.cc:239
SUNDIALS::ARKode::set_functions_to_trigger_an_assert
void set_functions_to_trigger_an_assert()
Definition: arkode.cc:494
DEAL_II_SUNDIALS_VERSION_LT
#define DEAL_II_SUNDIALS_VERSION_LT(major, minor, patch)
Definition: config.h:272
arkode.h
StandardExceptions::ExcMPI
Definition: exceptions.h:1168
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
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
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
Generated by   doxygen 1.8.17