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