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