Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
arkode.cc
Go to the documentation of this file.
1//-----------------------------------------------------------
2//
3// Copyright (C) 2017 - 2021 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
20
21#ifdef DEAL_II_WITH_SUNDIALS
22
25
29# ifdef DEAL_II_WITH_TRILINOS
32# endif
33# ifdef DEAL_II_WITH_PETSC
36# endif
37
40
41# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
42# include <arkode/arkode_impl.h>
43# include <sundials/sundials_config.h>
44# else
45# include <arkode/arkode_arkstep.h>
46# include <sunlinsol/sunlinsol_spgmr.h>
47# include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
48# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
50# endif
51# endif
52
53# include <iostream>
54
55# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
56// Make sure we know how to call sundials own ARKode() function
57const auto &SundialsARKode = ARKode;
58# endif
59
61
62namespace SUNDIALS
63{
64 using namespace internal;
65
66 namespace
67 {
68 template <typename VectorType>
69 int
70 t_arkode_explicit_function(realtype tt,
71 N_Vector yy,
72 N_Vector yp,
73 void * user_data)
74 {
75 Assert(user_data != nullptr, ExcInternalError());
76 ARKode<VectorType> &solver =
77 *static_cast<ARKode<VectorType> *>(user_data);
78
79 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
80 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
81
82 return solver.explicit_function(tt, *src_yy, *dst_yp);
83 }
84
85
86
87 template <typename VectorType>
88 int
89 t_arkode_implicit_function(realtype tt,
90 N_Vector yy,
91 N_Vector yp,
92 void * user_data)
93 {
94 Assert(user_data != nullptr, ExcInternalError());
95 ARKode<VectorType> &solver =
96 *static_cast<ARKode<VectorType> *>(user_data);
97
98 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
99 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
100
101 return solver.implicit_function(tt, *src_yy, *dst_yp);
102 }
103
104
105
106# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
107 template <typename VectorType>
108 int
109 t_arkode_setup_jacobian(ARKodeMem arkode_mem,
110 int convfail,
111 N_Vector ypred,
112 N_Vector fpred,
113 booleantype *jcurPtr,
114 N_Vector,
115 N_Vector,
116 N_Vector)
117 {
118 Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
119 ARKode<VectorType> &solver =
120 *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
121
122 auto *src_ypred = internal::unwrap_nvector_const<VectorType>(ypred);
123 auto *src_fpred = internal::unwrap_nvector_const<VectorType>(fpred);
124 // avoid reinterpret_cast
125 bool jcurPtr_tmp = false;
126 int err = solver.setup_jacobian(convfail,
127 arkode_mem->ark_tn,
128 arkode_mem->ark_gamma,
129 *src_ypred,
130 *src_fpred,
131 jcurPtr_tmp);
132# if DEAL_II_SUNDIALS_VERSION_GTE(2, 0, 0)
133 *jcurPtr = jcurPtr_tmp ? SUNTRUE : SUNFALSE;
134# else
135 *jcurPtr = jcurPtr_tmp ? TRUE : FALSE;
136# endif
137
138 return err;
139 }
140
141
142
143 template <typename VectorType>
144 int
145 t_arkode_solve_jacobian(ARKodeMem arkode_mem,
146 N_Vector b,
147# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
148 N_Vector,
149# endif
150 N_Vector ycur,
151 N_Vector fcur)
152 {
153 Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
154 ARKode<VectorType> &solver =
155 *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
156
157 auto *dst = internal::unwrap_nvector<VectorType>(b);
158 auto *src_ycur = internal::unwrap_nvector_const<VectorType>(ycur);
159 auto *src_fcur = internal::unwrap_nvector_const<VectorType>(fcur);
160
161 // make a temporary copy to work on in the user call
162 VectorType src = *dst;
163
164 int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
165 arkode_mem->ark_gamma,
166 *src_ycur,
167 *src_fcur,
168 src,
169 *dst);
170
171
172 return err;
173 }
174
175
176
177 template <typename VectorType>
178 int
179 t_arkode_setup_mass(ARKodeMem arkode_mem, N_Vector, N_Vector, N_Vector)
180 {
181 Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
182 ARKode<VectorType> &solver =
183 *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
184 int err = solver.setup_mass(arkode_mem->ark_tn);
185 return err;
186 }
187
188
189
190 template <typename VectorType>
191 int
192 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 Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
202 ARKode<VectorType> &solver =
203 *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
204
205 auto *dst = internal::unwrap_nvector<VectorType>(b);
206
207 // make a temporary copy to work on in the user call
208 VectorType src = *dst;
209
210 int err = solver.solve_mass_system(src, *dst);
211
212 return err;
213 }
214# else
215
216 template <typename VectorType>
217 int
218 t_arkode_jac_times_vec_function(N_Vector v,
219 N_Vector Jv,
220 realtype t,
221 N_Vector y,
222 N_Vector fy,
223 void * user_data,
224 N_Vector)
225 {
226 Assert(user_data != nullptr, ExcInternalError());
227 ARKode<VectorType> &solver =
228 *static_cast<ARKode<VectorType> *>(user_data);
229
230 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
231 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
232 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
233
234 auto *dst_Jv = internal::unwrap_nvector<VectorType>(Jv);
235
236 return solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy);
237 }
238
239
240
241 template <typename VectorType>
242 int
243 t_arkode_jac_times_setup_function(realtype t,
244 N_Vector y,
245 N_Vector fy,
246 void * user_data)
247 {
248 Assert(user_data != nullptr, ExcInternalError());
249 ARKode<VectorType> &solver =
250 *static_cast<ARKode<VectorType> *>(user_data);
251
252 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
253 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
254
255 return solver.jacobian_times_setup(t, *src_y, *src_fy);
256 }
257
258
259
260 template <typename VectorType>
261 int
262 t_arkode_prec_solve_function(realtype t,
263 N_Vector y,
264 N_Vector fy,
265 N_Vector r,
266 N_Vector z,
267 realtype gamma,
268 realtype delta,
269 int lr,
270 void * user_data)
271 {
272 Assert(user_data != nullptr, ExcInternalError());
273 ARKode<VectorType> &solver =
274 *static_cast<ARKode<VectorType> *>(user_data);
275
276 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
277 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
278 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
279
280 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
281
282 return solver.jacobian_preconditioner_solve(
283 t, *src_y, *src_fy, *src_r, *dst_z, gamma, delta, lr);
284 }
285
286
287
288 template <typename VectorType>
289 int
290 t_arkode_prec_setup_function(realtype t,
291 N_Vector y,
292 N_Vector fy,
293 booleantype jok,
294 booleantype *jcurPtr,
295 realtype gamma,
296 void * user_data)
297 {
298 Assert(user_data != nullptr, ExcInternalError());
299 ARKode<VectorType> &solver =
300 *static_cast<ARKode<VectorType> *>(user_data);
301
302 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
303 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
304
305 return solver.jacobian_preconditioner_setup(
306 t, *src_y, *src_fy, jok, *jcurPtr, gamma);
307 }
308
309
310
311 template <typename VectorType>
312 int
313 t_arkode_mass_times_vec_function(N_Vector v,
314 N_Vector Mv,
315 realtype t,
316 void * mtimes_data)
317 {
318 Assert(mtimes_data != nullptr, ExcInternalError());
319 ARKode<VectorType> &solver =
320 *static_cast<ARKode<VectorType> *>(mtimes_data);
321
322 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
323 auto *dst_Mv = internal::unwrap_nvector<VectorType>(Mv);
324
325 return solver.mass_times_vector(t, *src_v, *dst_Mv);
326 }
327
328
329
330 template <typename VectorType>
331 int
332 t_arkode_mass_times_setup_function(realtype t, void *mtimes_data)
333 {
334 Assert(mtimes_data != nullptr, ExcInternalError());
335 ARKode<VectorType> &solver =
336 *static_cast<ARKode<VectorType> *>(mtimes_data);
337
338 return solver.mass_times_setup(t);
339 }
340
341
342
343 template <typename VectorType>
344 int
345 t_arkode_mass_prec_solve_function(realtype t,
346 N_Vector r,
347 N_Vector z,
348 realtype delta,
349 int lr,
350 void * user_data)
351 {
352 Assert(user_data != nullptr, ExcInternalError());
353 ARKode<VectorType> &solver =
354 *static_cast<ARKode<VectorType> *>(user_data);
355
356 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
357 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
358
359 return solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr);
360 }
361
362
363
364 template <typename VectorType>
365 int
366 t_arkode_mass_prec_setup_function(realtype t, void *user_data)
367 {
368 Assert(user_data != nullptr, ExcInternalError());
369 ARKode<VectorType> &solver =
370 *static_cast<ARKode<VectorType> *>(user_data);
371
372 return solver.mass_preconditioner_setup(t);
373 }
374# endif
375
376 } // namespace
377
378
379 template <typename VectorType>
381 const MPI_Comm & mpi_comm)
382 : data(data)
383 , arkode_mem(nullptr)
384 , communicator(is_serial_vector<VectorType>::value ?
385 MPI_COMM_SELF :
386 Utilities::MPI::duplicate_communicator(mpi_comm))
387 , last_end_time(data.initial_time)
388 {
390 }
391
392 template <typename VectorType>
394 {
395 if (arkode_mem)
396# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
397 ARKodeFree(&arkode_mem);
398# else
399 ARKStepFree(&arkode_mem);
400# endif
401
402# ifdef DEAL_II_WITH_MPI
404 {
405 const int ierr = MPI_Comm_free(&communicator);
406 (void)ierr;
407 AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
408 }
409# endif
410 }
411
412
413
414 template <typename VectorType>
415 unsigned int
416 ARKode<VectorType>::solve_ode(VectorType &solution)
417 {
418 DiscreteTime time(data.initial_time,
419 data.final_time,
420 data.initial_step_size);
421
422 return do_evolve_time(solution, time, /* force_solver_restart = */ true);
423 }
424
425
426
427 template <typename VectorType>
428 unsigned int
430 const double intermediate_time,
431 const bool reset_solver)
432 {
434 intermediate_time > last_end_time,
436 "The requested intermediate time is smaller than the last requested "
437 "intermediate time."));
438
439 const bool do_reset = reset_solver || arkode_mem == nullptr;
440 DiscreteTime time(last_end_time, intermediate_time, data.initial_step_size);
441 return do_evolve_time(solution, time, do_reset);
442 }
443
444
445
446 template <typename VectorType>
447 int
449 DiscreteTime &time,
450 const bool do_reset)
451 {
452 if (do_reset)
453 {
454 reset(time.get_current_time(), time.get_next_step_size(), solution);
455 if (output_step)
456 output_step(time.get_current_time(),
457 solution,
458 time.get_step_number());
459 }
460
461 auto solution_nvector = internal::make_nvector_view(solution);
462
463 while (!time.is_at_end())
464 {
465 time.set_desired_next_step_size(data.output_period);
466 double actual_next_time;
467# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
468 const auto status = SundialsARKode(arkode_mem,
469 time.get_next_time(),
470 solution_nvector,
471 &actual_next_time,
472 ARK_NORMAL);
473# else
474 const auto status = ARKStepEvolve(arkode_mem,
475 time.get_next_time(),
476 solution_nvector,
477 &actual_next_time,
478 ARK_NORMAL);
479# endif
480 (void)status;
481 AssertARKode(status);
482
483 time.set_next_step_size(actual_next_time - time.get_current_time());
484 time.advance_time();
485
486 while (solver_should_restart(time.get_current_time(), solution))
487 reset(time.get_current_time(),
489 solution);
490
491 if (output_step)
492 output_step(time.get_current_time(),
493 solution,
494 time.get_step_number());
495 }
496 last_end_time = time.get_current_time();
497 return time.get_step_number();
498 }
499
500
501
502# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
503 template <typename VectorType>
504 void
505 ARKode<VectorType>::reset(const double current_time,
506 const double current_time_step,
507 const VectorType &solution)
508 {
509 last_end_time = current_time;
510 if (arkode_mem)
511 ARKodeFree(&arkode_mem);
512
513 arkode_mem = ARKodeCreate();
514
515 int status;
516 (void)status;
517
518 Assert(explicit_function || implicit_function,
519 ExcFunctionNotProvided("explicit_function || implicit_function"));
520
521 // just a view on the memory in solution, all write operations on yy by
522 // ARKODE will automatically be mirrored to solution
523 auto initial_condition_nvector = internal::make_nvector_view(solution);
524
525 status = ARKodeInit(
526 arkode_mem,
527 explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
528 implicit_function ? &t_arkode_implicit_function<VectorType> : nullptr,
529 current_time,
530 initial_condition_nvector);
531 AssertARKode(status);
532
533 if (get_local_tolerances)
534 {
535 const auto abs_tols = make_nvector_view(get_local_tolerances());
536 status =
537 ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
538 AssertARKode(status);
539 }
540 else
541 {
542 status = ARKodeSStolerances(arkode_mem,
543 data.relative_tolerance,
544 data.absolute_tolerance);
545 AssertARKode(status);
546 }
547
548 status = ARKodeSetInitStep(arkode_mem, current_time_step);
549 AssertARKode(status);
550
551 status = ARKodeSetUserData(arkode_mem, this);
552 AssertARKode(status);
553
554 status = ARKodeSetStopTime(arkode_mem, data.final_time);
555 AssertARKode(status);
556
557 status =
558 ARKodeSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
559 AssertARKode(status);
560
561 // Initialize solver
562 auto ARKode_mem = static_cast<ARKodeMem>(arkode_mem);
563
564 if (solve_jacobian_system)
565 {
566 status = ARKodeSetNewton(arkode_mem);
567 AssertARKode(status);
568 if (data.implicit_function_is_linear)
569 {
570 status = ARKodeSetLinear(
571 arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
572 AssertARKode(status);
573 }
574
575
576 ARKode_mem->ark_lsolve = t_arkode_solve_jacobian<VectorType>;
577 if (setup_jacobian)
578 {
579 ARKode_mem->ark_lsetup = t_arkode_setup_jacobian<VectorType>;
580# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
581 ARKode_mem->ark_setupNonNull = true;
582# endif
583 }
584 }
585 else
586 {
587 status =
588 ARKodeSetFixedPoint(arkode_mem, data.maximum_non_linear_iterations);
589 AssertARKode(status);
590 }
591
592
593 if (solve_mass_system)
594 {
595 ARKode_mem->ark_msolve = t_arkode_solve_mass<VectorType>;
596
597 if (setup_mass)
598 {
599 ARKode_mem->ark_msetup = t_arkode_setup_mass<VectorType>;
600# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
601 ARKode_mem->ark_MassSetupNonNull = true;
602# endif
603 }
604 }
605
606 status = ARKodeSetOrder(arkode_mem, data.maximum_order);
607 AssertARKode(status);
608
609 if (custom_setup)
610 custom_setup(arkode_mem);
611 }
612
613# else
614
615 template <typename VectorType>
616 void
617 ARKode<VectorType>::reset(const double current_time,
618 const double current_time_step,
619 const VectorType &solution)
620 {
621 last_end_time = current_time;
622 if (arkode_mem)
623 ARKStepFree(&arkode_mem);
624
625 int status;
626 (void)status;
627
628 // just a view on the memory in solution, all write operations on yy by
629 // ARKODE will automatically be mirrored to solution
630 auto initial_condition_nvector = internal::make_nvector_view(solution);
631
632 Assert(explicit_function || implicit_function,
633 ExcFunctionNotProvided("explicit_function || implicit_function"));
634
635 arkode_mem = ARKStepCreate(
636 explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
637 implicit_function ? &t_arkode_implicit_function<VectorType> : nullptr,
638 current_time,
639 initial_condition_nvector);
640
641 Assert(arkode_mem != nullptr, ExcInternalError());
642
643 if (get_local_tolerances)
644 {
645 const auto abs_tols =
646 internal::make_nvector_view(get_local_tolerances());
647 status =
648 ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
649 AssertARKode(status);
650 }
651 else
652 {
653 status = ARKStepSStolerances(arkode_mem,
654 data.relative_tolerance,
655 data.absolute_tolerance);
656 AssertARKode(status);
657 }
658
659 // for version 4.0.0 this call must be made before solver settings
660 status = ARKStepSetUserData(arkode_mem, this);
661 AssertARKode(status);
662
663 setup_system_solver(solution);
664
665 setup_mass_solver(solution);
666
667 status = ARKStepSetInitStep(arkode_mem, current_time_step);
668 AssertARKode(status);
669
670 status = ARKStepSetStopTime(arkode_mem, data.final_time);
671 AssertARKode(status);
672
673 status = ARKStepSetOrder(arkode_mem, data.maximum_order);
674 AssertARKode(status);
675
676 if (custom_setup)
677 custom_setup(arkode_mem);
678 }
679
680
681
682 template <typename VectorType>
683 void
684 ARKode<VectorType>::setup_system_solver(const VectorType &solution)
685 {
686 // no point in setting up a solver when there is no linear system
687 if (!implicit_function)
688 return;
689
690 int status;
691 (void)status;
692 // Initialize solver
693 // currently only iterative linear solvers are supported
694 if (jacobian_times_vector)
695 {
696 SUNLinearSolver sun_linear_solver;
697 if (solve_linearized_system)
698 {
699 linear_solver =
700 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
701 solve_linearized_system);
702 sun_linear_solver = *linear_solver;
703 }
704 else
705 {
706 // use default solver from SUNDIALS
707 // TODO give user options
708 auto y_template = internal::make_nvector_view(solution);
709 sun_linear_solver =
710 SUNLinSol_SPGMR(y_template,
711 PREC_NONE,
712 0 /*krylov subvectors, 0 uses default*/);
713 }
714 status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr);
715 AssertARKode(status);
716 status =
717 ARKStepSetJacTimes(arkode_mem,
718 jacobian_times_setup ?
719 t_arkode_jac_times_setup_function<VectorType> :
720 nullptr,
721 t_arkode_jac_times_vec_function<VectorType>);
722 AssertARKode(status);
723 if (jacobian_preconditioner_solve)
724 {
725 status = ARKStepSetPreconditioner(
726 arkode_mem,
727 jacobian_preconditioner_setup ?
728 t_arkode_prec_setup_function<VectorType> :
729 nullptr,
730 t_arkode_prec_solve_function<VectorType>);
731 AssertARKode(status);
732 }
733 if (data.implicit_function_is_linear)
734 {
735 status = ARKStepSetLinear(
736 arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
737 AssertARKode(status);
738 }
739 }
740 else
741 {
742 auto y_template = internal::make_nvector_view(solution);
743
744 SUNNonlinearSolver fixed_point_solver =
745 SUNNonlinSol_FixedPoint(y_template,
746 data.anderson_acceleration_subspace);
747
748 status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
749 AssertARKode(status);
750 }
751
752 status =
753 ARKStepSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
754 AssertARKode(status);
755 }
756
757
758
759 template <typename VectorType>
760 void
761 ARKode<VectorType>::setup_mass_solver(const VectorType &solution)
762 {
763 int status;
764 (void)status;
765
766 if (mass_times_vector)
767 {
768 SUNLinearSolver sun_mass_linear_solver;
769 if (solve_mass)
770 {
771 mass_solver =
772 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
773 solve_mass);
774 sun_mass_linear_solver = *mass_solver;
775 }
776 else
777 {
778 auto y_template = internal::make_nvector_view(solution);
779 sun_mass_linear_solver =
780 SUNLinSol_SPGMR(y_template,
781 PREC_NONE,
782 0 /*krylov subvectors, 0 uses default*/);
783 }
784 booleantype mass_time_dependent =
785 data.mass_is_time_independent ? SUNFALSE : SUNTRUE;
786 status = ARKStepSetMassLinearSolver(arkode_mem,
787 sun_mass_linear_solver,
788 nullptr,
789 mass_time_dependent);
790 AssertARKode(status);
791 status =
792 ARKStepSetMassTimes(arkode_mem,
793 mass_times_setup ?
794 t_arkode_mass_times_setup_function<VectorType> :
795 nullptr,
796 t_arkode_mass_times_vec_function<VectorType>,
797 this);
798 AssertARKode(status);
799
800 if (mass_preconditioner_solve)
801 {
802 status = ARKStepSetMassPreconditioner(
803 arkode_mem,
804 mass_preconditioner_setup ?
805 t_arkode_mass_prec_setup_function<VectorType> :
806 nullptr,
807 t_arkode_mass_prec_solve_function<VectorType>);
808 AssertARKode(status);
809 }
810 }
811 }
812# endif
813
814
815
816 template <typename VectorType>
817 void
819 {
820 reinit_vector = [](VectorType &) {
821 AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
822 };
823
824 solver_should_restart = [](const double, VectorType &) -> bool {
825 return false;
826 };
827 }
828
829
830
831 template <typename VectorType>
832 void *
834 {
835 return arkode_mem;
836 }
837
838 template class ARKode<Vector<double>>;
839 template class ARKode<BlockVector<double>>;
840
843
844# ifdef DEAL_II_WITH_MPI
845
846# ifdef DEAL_II_WITH_TRILINOS
849
850# endif // DEAL_II_WITH_TRILINOS
851
852# ifdef DEAL_II_WITH_PETSC
853# ifndef PETSC_USE_COMPLEX
856# endif // PETSC_USE_COMPLEX
857# endif // DEAL_II_WITH_PETSC
858
859# endif // DEAL_II_WITH_MPI
860
861} // namespace SUNDIALS
862
864
865#endif
const auto & SundialsARKode
Definition: arkode.cc:57
#define AssertARKode(code)
Definition: arkode.h:62
double get_next_step_size() const
void set_desired_next_step_size(const double time_step_size)
unsigned int get_step_number() const
bool is_at_end() const
double get_current_time() const
void set_next_step_size(const double time_step_size)
double get_next_time() const
double get_previous_step_size() const
void advance_time()
void setup_mass_solver(const VectorType &solution)
Definition: arkode.cc:761
void * get_arkode_memory() const
Definition: arkode.cc:833
ARKode(const AdditionalData &data=AdditionalData(), const MPI_Comm &mpi_comm=MPI_COMM_WORLD)
Definition: arkode.cc:380
int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
Definition: arkode.cc:448
void reset(const double t, const double h, const VectorType &y)
Definition: arkode.cc:505
void set_functions_to_trigger_an_assert()
Definition: arkode.cc:818
unsigned int solve_ode(VectorType &solution)
Definition: arkode.cc:416
void setup_system_solver(const VectorType &solution)
Definition: arkode.cc:684
unsigned int solve_ode_incrementally(VectorType &solution, const double intermediate_time, const bool reset_solver=false)
Definition: arkode.cc:429
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_SUNDIALS_VERSION_LT(major, minor, patch)
Definition: config.h:314
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1528
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:160
long double gamma(const unsigned int n)