17 #include <deal.II/base/config.h> 19 #include <deal.II/sundials/arkode.h> 21 #ifdef DEAL_II_WITH_SUNDIALS 23 # include <deal.II/base/utilities.h> 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> 30 # ifdef DEAL_II_WITH_PETSC 31 # include <deal.II/lac/petsc_block_vector.h> 32 # include <deal.II/lac/petsc_vector.h> 34 # include <deal.II/base/utilities.h> 36 # include <deal.II/sundials/copy.h> 38 # include <arkode/arkode_impl.h> 39 # include <sundials/sundials_config.h> 45 const auto &SundialsARKode = ARKode;
47 DEAL_II_NAMESPACE_OPEN
55 template <
typename VectorType>
57 t_arkode_explicit_function(realtype tt,
62 ARKode<VectorType> &solver =
63 *
static_cast<ARKode<VectorType> *
>(user_data);
67 solver.reinit_vector(*src_yy);
70 solver.reinit_vector(*dst_yp);
74 int err = solver.explicit_function(tt, *src_yy, *dst_yp);
83 template <
typename VectorType>
85 t_arkode_implicit_function(realtype tt,
90 ARKode<VectorType> &solver =
91 *
static_cast<ARKode<VectorType> *
>(user_data);
95 solver.reinit_vector(*src_yy);
98 solver.reinit_vector(*dst_yp);
102 int err = solver.implicit_function(tt, *src_yy, *dst_yp);
111 template <
typename VectorType>
113 t_arkode_setup_jacobian(ARKodeMem arkode_mem,
117 booleantype *jcurPtr,
122 ARKode<VectorType> &solver =
123 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
127 solver.reinit_vector(*src_ypred);
130 solver.reinit_vector(*src_fpred);
132 copy(*src_ypred, ypred);
133 copy(*src_fpred, fpred);
135 int err = solver.setup_jacobian(convfail,
137 arkode_mem->ark_gamma,
147 template <
typename VectorType>
149 t_arkode_solve_jacobian(ARKodeMem arkode_mem,
151 #
if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
157 ARKode<VectorType> &solver =
158 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
162 solver.reinit_vector(*src);
165 solver.reinit_vector(*src_ycur);
168 solver.reinit_vector(*src_fcur);
171 solver.reinit_vector(*dst);
174 copy(*src_ycur, ycur);
175 copy(*src_fcur, fcur);
177 int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
178 arkode_mem->ark_gamma,
190 template <
typename VectorType>
192 t_arkode_setup_mass(ARKodeMem arkode_mem, N_Vector, N_Vector, N_Vector)
194 ARKode<VectorType> &solver =
195 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
196 int err = solver.setup_mass(arkode_mem->ark_tn);
202 template <
typename VectorType>
204 t_arkode_solve_mass(ARKodeMem arkode_mem,
205 #
if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
213 ARKode<VectorType> &solver =
214 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
218 solver.reinit_vector(*src);
221 solver.reinit_vector(*dst);
225 int err = solver.solve_mass_system(*src, *dst);
232 template <
typename VectorType>
234 const MPI_Comm mpi_comm)
236 , arkode_mem(nullptr)
241 Utilities::MPI::duplicate_communicator(mpi_comm))
246 template <
typename VectorType>
250 ARKodeFree(&arkode_mem);
251 # ifdef DEAL_II_WITH_MPI 263 template <
typename VectorType>
267 unsigned int system_size = solution.size();
271 unsigned int step_number = 0;
279 # ifdef DEAL_II_WITH_MPI 282 const IndexSet is = solution.locally_owned_elements();
283 const size_t local_system_size = is.
n_elements();
285 yy = N_VNew_Parallel(
communicator, local_system_size, system_size);
288 N_VNew_Parallel(
communicator, local_system_size, system_size);
295 "Trying to use a serial code with a parallel vector."));
296 yy = N_VNew_Serial(system_size);
310 status = SundialsARKode(arkode_mem, next_time, yy, &t, ARK_NORMAL);
312 AssertARKode(status);
314 status = ARKodeGetLastStep(arkode_mem, &h);
315 AssertARKode(status);
320 reset(t, h, solution);
329 # ifdef DEAL_II_WITH_MPI 332 N_VDestroy_Parallel(yy);
338 N_VDestroy_Serial(yy);
345 template <
typename VectorType>
348 const double & current_time_step,
349 const VectorType &solution)
351 unsigned int system_size;
354 ARKodeFree(&arkode_mem);
356 arkode_mem = ARKodeCreate();
361 # ifdef DEAL_II_WITH_MPI 364 N_VDestroy_Parallel(yy);
370 N_VDestroy_Serial(yy);
377 system_size = solution.size();
378 # ifdef DEAL_II_WITH_MPI 381 const IndexSet is = solution.locally_owned_elements();
382 const size_t local_system_size = is.
n_elements();
384 yy = N_VNew_Parallel(
communicator, local_system_size, system_size);
387 N_VNew_Parallel(
communicator, local_system_size, system_size);
392 yy = N_VNew_Serial(system_size);
399 ExcFunctionNotProvided(
"explicit_function || implicit_function"));
407 AssertARKode(status);
414 AssertARKode(status);
418 status = ARKodeSStolerances(arkode_mem,
421 AssertARKode(status);
424 status = ARKodeSetInitStep(arkode_mem, current_time_step);
425 AssertARKode(status);
427 status = ARKodeSetUserData(arkode_mem, (
void *)
this);
428 AssertARKode(status);
431 AssertARKode(status);
435 AssertARKode(status);
438 ARKodeMem ARKode_mem = (ARKodeMem)arkode_mem;
442 status = ARKodeSetNewton(arkode_mem);
443 AssertARKode(status);
446 status = ARKodeSetLinear(
448 AssertARKode(status);
452 ARKode_mem->ark_lsolve = t_arkode_solve_jacobian<VectorType>;
455 ARKode_mem->ark_lsetup = t_arkode_setup_jacobian<VectorType>;
456 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) 457 ARKode_mem->ark_setupNonNull =
true;
465 AssertARKode(status);
471 ARKode_mem->ark_msolve = t_arkode_solve_mass<VectorType>;
475 ARKode_mem->ark_msetup = t_arkode_setup_mass<VectorType>;
476 # if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) 477 ARKode_mem->ark_MassSetupNonNull =
true;
483 AssertARKode(status);
486 template <
typename VectorType>
491 AssertThrow(
false, ExcFunctionNotProvided(
"reinit_vector"));
502 # ifdef DEAL_II_WITH_MPI 504 # ifdef DEAL_II_WITH_TRILINOS 507 # endif // DEAL_II_WITH_TRILINOS 509 # ifdef DEAL_II_WITH_PETSC 510 # ifndef PETSC_USE_COMPLEX 513 # endif // PETSC_USE_COMPLEX 514 # endif // DEAL_II_WITH_PETSC 516 # endif // DEAL_II_WITH_MPI 520 DEAL_II_NAMESPACE_CLOSE
bool implicit_function_is_time_independent
#define AssertNothrow(cond, exc)
size_type n_elements() const
std::function< int(const VectorType &rhs, VectorType &dst)> solve_mass_system
double relative_tolerance
#define AssertThrow(cond, exc)
std::function< void(const double t, const VectorType &sol, const unsigned int step_number)> output_step
void reset(const double &t, const double &h, const VectorType &y)
std::function< VectorType &()> get_local_tolerances
unsigned int maximum_non_linear_iterations
std::function< int(const double t)> setup_mass
std::function< bool(const double t, VectorType &sol)> solver_should_restart
std::function< int(const double t, const double gamma, const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType &dst)> solve_jacobian_system
#define Assert(cond, exc)
std::function< int(const double t, const VectorType &y, VectorType &res)> implicit_function
bool implicit_function_is_linear
unsigned int maximum_order
std::function< int(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function
std::function< void(VectorType &)> reinit_vector
unsigned int solve_ode(VectorType &solution)
std::function< int(const int convfail, const double t, const double gamma, const VectorType &ypred, const VectorType &fpred, bool &j_is_current)> setup_jacobian
void set_functions_to_trigger_an_assert()
ARKode(const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD)
static::ExceptionBase & ExcInternalError()
double absolute_tolerance