From: Daniel Arndt Date: Sun, 19 Nov 2023 15:34:38 +0000 (-0500) Subject: Introduce sundials_types.h X-Git-Tag: relicensing~289^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16271%2Fhead;p=dealii.git Introduce sundials_types.h --- diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index 7e254952cb..f969eb80e7 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -43,13 +43,13 @@ # include # include +# include # include # include # include # include -# include # include # include diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 45dcc0d141..9e6acf22b2 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -39,6 +39,7 @@ # include # endif +# include # include # include @@ -46,7 +47,6 @@ # include # include # include -# include # include diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 1c980951da..3f098c7d77 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -31,12 +31,13 @@ # include # include +# include + # include # include # include # include -# include # include # include diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index 1fe1c0cf9a..12455a9c22 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -20,6 +20,7 @@ #include #include +#include #ifdef DEAL_II_WITH_SUNDIALS @@ -176,12 +177,6 @@ namespace SUNDIALS */ namespace NVectorOperations { -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - using realtype = ::sunrealtype; -# else - using realtype = ::realtype; -# endif - N_Vector_ID get_vector_id(N_Vector v); @@ -202,18 +197,22 @@ namespace SUNDIALS template void - linear_sum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z); + linear_sum(SUNDIALS::realtype a, + N_Vector x, + SUNDIALS::realtype b, + N_Vector y, + N_Vector z); template - realtype + SUNDIALS::realtype dot_product(N_Vector x, N_Vector y); template - realtype + SUNDIALS::realtype weighted_l2_norm(N_Vector x, N_Vector y); template - realtype + SUNDIALS::realtype l1_norm(N_Vector x); template @@ -251,47 +250,47 @@ namespace SUNDIALS elementwise_abs(N_Vector x, N_Vector z); template - realtype + SUNDIALS::realtype weighted_rms_norm(N_Vector x, N_Vector w); template - realtype + SUNDIALS::realtype weighted_rms_norm_mask(N_Vector x, N_Vector w, N_Vector mask); template - realtype + SUNDIALS::realtype max_norm(N_Vector x); template ::value, int> = 0> - realtype + SUNDIALS::realtype min_element(N_Vector x); template ::value && !IsBlockVector::value, int> = 0> - realtype + SUNDIALS::realtype min_element(N_Vector x); template ::value && IsBlockVector::value, int> = 0> - realtype + SUNDIALS::realtype min_element(N_Vector x); template void - scale(realtype c, N_Vector x, N_Vector z); + scale(SUNDIALS::realtype c, N_Vector x, N_Vector z); template void - set_constant(realtype c, N_Vector v); + set_constant(SUNDIALS::realtype c, N_Vector v); template void - add_constant(N_Vector x, realtype b, N_Vector z); + add_constant(N_Vector x, SUNDIALS::realtype b, N_Vector z); template const MPI_Comm & @@ -652,7 +651,11 @@ namespace SUNDIALS template void - linear_sum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z) + linear_sum(SUNDIALS::realtype a, + N_Vector x, + SUNDIALS::realtype b, + N_Vector y, + N_Vector z) { auto *x_dealii = unwrap_nvector_const(x); auto *y_dealii = unwrap_nvector_const(y); @@ -672,7 +675,7 @@ namespace SUNDIALS template - realtype + SUNDIALS::realtype dot_product(N_Vector x, N_Vector y) { return *unwrap_nvector_const(x) * @@ -682,7 +685,7 @@ namespace SUNDIALS template - realtype + SUNDIALS::realtype weighted_l2_norm(N_Vector x, N_Vector w) { // TODO copy can be avoided by a custom kernel @@ -695,7 +698,7 @@ namespace SUNDIALS template - realtype + SUNDIALS::realtype l1_norm(N_Vector x) { return unwrap_nvector_const(x)->l1_norm(); @@ -705,7 +708,7 @@ namespace SUNDIALS template void - set_constant(realtype c, N_Vector v) + set_constant(SUNDIALS::realtype c, N_Vector v) { auto *v_dealii = unwrap_nvector(v); *v_dealii = c; @@ -715,7 +718,7 @@ namespace SUNDIALS template void - add_constant(N_Vector x, realtype b, N_Vector z) + add_constant(N_Vector x, SUNDIALS::realtype b, N_Vector z) { auto *x_dealii = unwrap_nvector_const(x); auto *z_dealii = unwrap_nvector(z); @@ -752,7 +755,7 @@ namespace SUNDIALS template - realtype + SUNDIALS::realtype weighted_rms_norm(N_Vector x, N_Vector w) { // TODO copy can be avoided by a custom kernel @@ -766,7 +769,7 @@ namespace SUNDIALS template - realtype + SUNDIALS::realtype weighted_rms_norm_mask(N_Vector x, N_Vector w, N_Vector mask) { // TODO copy can be avoided by a custom kernel @@ -782,7 +785,7 @@ namespace SUNDIALS template - realtype + SUNDIALS::realtype max_norm(N_Vector x) { return unwrap_nvector_const(x)->linfty_norm(); @@ -792,7 +795,7 @@ namespace SUNDIALS template ::value, int>> - realtype + SUNDIALS::realtype min_element(N_Vector x) { auto *vector = unwrap_nvector_const(x); @@ -805,7 +808,7 @@ namespace SUNDIALS std::enable_if_t::value && !IsBlockVector::value, int>> - realtype + SUNDIALS::realtype min_element(N_Vector x) { auto *vector = unwrap_nvector_const(x); @@ -831,7 +834,7 @@ namespace SUNDIALS std::enable_if_t::value && IsBlockVector::value, int>> - realtype + SUNDIALS::realtype min_element(N_Vector x) { auto *vector = unwrap_nvector_const(x); @@ -869,7 +872,7 @@ namespace SUNDIALS template void - scale(realtype c, N_Vector x, N_Vector z) + scale(SUNDIALS::realtype c, N_Vector x, N_Vector z) { auto *x_dealii = unwrap_nvector_const(x); auto *z_dealii = unwrap_nvector(z); diff --git a/include/deal.II/sundials/sundials_types.h b/include/deal.II/sundials/sundials_types.h new file mode 100644 index 0000000000..66db338bb6 --- /dev/null +++ b/include/deal.II/sundials/sundials_types.h @@ -0,0 +1,41 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_sundials_types_h +#define dealii_sundials_types_h + +#include + +#ifdef DEAL_II_WITH_SUNDIALS +# include + +DEAL_II_NAMESPACE_OPEN + +namespace SUNDIALS +{ +/** + * Alias for the real type used by SUNDIALS. + */ +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) + using realtype = ::sunrealtype; +# else + using realtype = ::realtype; +# endif +} // namespace SUNDIALS + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_WITH_SUNDIALS +#endif // dealii_sundials_types_h diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index a27a52a06e..abb320af76 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -269,15 +269,10 @@ namespace SUNDIALS Assert(explicit_function || implicit_function, ExcFunctionNotProvided("explicit_function || implicit_function")); - auto explicit_function_callback = []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype tt, -# else - realtype tt, -# endif - N_Vector yy, - N_Vector yp, - void *user_data) -> int { + auto explicit_function_callback = [](SUNDIALS::realtype tt, + N_Vector yy, + N_Vector yp, + void *user_data) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); @@ -294,15 +289,10 @@ namespace SUNDIALS }; - auto implicit_function_callback = []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype tt, -# else - realtype tt, -# endif - N_Vector yy, - N_Vector yp, - void *user_data) -> int { + auto implicit_function_callback = [](SUNDIALS::realtype tt, + N_Vector yy, + N_Vector yp, + void *user_data) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); @@ -428,16 +418,12 @@ namespace SUNDIALS status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr); AssertARKode(status); - auto jacobian_times_vector_callback = [](N_Vector v, - N_Vector Jv, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif - N_Vector y, - N_Vector fy, - void *user_data, + auto jacobian_times_vector_callback = [](N_Vector v, + N_Vector Jv, + SUNDIALS::realtype t, + N_Vector y, + N_Vector fy, + void *user_data, N_Vector) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = @@ -459,15 +445,10 @@ namespace SUNDIALS *src_fy); }; - auto jacobian_times_vector_setup_callback = []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif - N_Vector y, - N_Vector fy, - void *user_data) -> int { + auto jacobian_times_vector_setup_callback = [](SUNDIALS::realtype t, + N_Vector y, + N_Vector fy, + void *user_data) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); @@ -490,25 +471,15 @@ namespace SUNDIALS AssertARKode(status); if (jacobian_preconditioner_solve) { - auto solve_with_jacobian_callback = []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif - N_Vector y, - N_Vector fy, - N_Vector r, - N_Vector z, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype gamma, - sunrealtype delta, -# else - realtype gamma, - realtype delta, -# endif - int lr, - void *user_data) -> int { + auto solve_with_jacobian_callback = [](SUNDIALS::realtype t, + N_Vector y, + N_Vector fy, + N_Vector r, + N_Vector z, + SUNDIALS::realtype gamma, + SUNDIALS::realtype delta, + int lr, + void *user_data) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); @@ -532,22 +503,13 @@ namespace SUNDIALS lr); }; - auto jacobian_solver_setup_callback = []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif - N_Vector y, - N_Vector fy, - booleantype jok, - booleantype *jcurPtr, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype gamma, -# else - realtype gamma, -# endif - void *user_data) -> int { + auto jacobian_solver_setup_callback = [](SUNDIALS::realtype t, + N_Vector y, + N_Vector fy, + booleantype jok, + booleantype *jcurPtr, + SUNDIALS::realtype gamma, + void *user_data) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); @@ -665,13 +627,7 @@ namespace SUNDIALS AssertARKode(status); auto mass_matrix_times_vector_setup_callback = - []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif - void *mtimes_data) -> int { + [](SUNDIALS::realtype t, void *mtimes_data) -> int { Assert(mtimes_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(mtimes_data); @@ -680,13 +636,9 @@ namespace SUNDIALS solver.mass_times_setup, solver.pending_exception, t); }; - auto mass_matrix_times_vector_callback = [](N_Vector v, - N_Vector Mv, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif + auto mass_matrix_times_vector_callback = [](N_Vector v, + N_Vector Mv, + SUNDIALS::realtype t, void *mtimes_data) -> int { Assert(mtimes_data != nullptr, ExcInternalError()); ARKode &solver = @@ -713,13 +665,8 @@ namespace SUNDIALS if (mass_preconditioner_solve) { - auto mass_matrix_solver_setup_callback = []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif - void *user_data) -> int { + auto mass_matrix_solver_setup_callback = + [](SUNDIALS::realtype t, void *user_data) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); @@ -728,21 +675,12 @@ namespace SUNDIALS solver.mass_preconditioner_setup, solver.pending_exception, t); }; - auto solve_with_mass_matrix_callback = []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -# else - realtype t, -# endif - N_Vector r, - N_Vector z, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype delta, -# else - realtype delta, -# endif - int lr, - void *user_data) -> int { + auto solve_with_mass_matrix_callback = [](SUNDIALS::realtype t, + N_Vector r, + N_Vector z, + SUNDIALS::realtype delta, + int lr, + void *user_data) -> int { Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index ee26b25fa6..6d57319cb7 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -227,16 +227,11 @@ namespace SUNDIALS status = IDAInit( ida_mem, - []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype tt, -# else - realtype tt, -# endif - N_Vector yy, - N_Vector yp, - N_Vector rr, - void *user_data) -> int { + [](SUNDIALS::realtype tt, + N_Vector yy, + N_Vector yp, + N_Vector rr, + void *user_data) -> int { IDA &solver = *static_cast *>(user_data); auto *src_yy = internal::unwrap_nvector_const(yy); @@ -321,7 +316,7 @@ namespace SUNDIALS # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) LS = SUNLinSolNewEmpty(); # else - LS = SUNLinSolNewEmpty(ida_ctx); + LS = SUNLinSolNewEmpty(ida_ctx); # endif LS->content = this; @@ -349,14 +344,9 @@ namespace SUNDIALS ExcFunctionNotProvided("solve_with_jacobian")); LS->ops->solve = [](SUNLinearSolver LS, SUNMatrix /*ignored*/, - N_Vector x, - N_Vector b, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype tol -# else - realtype tol -# endif - ) -> int { + N_Vector x, + N_Vector b, + SUNDIALS::realtype tol) -> int { IDA &solver = *static_cast *>(LS->content); auto *src_b = internal::unwrap_nvector_const(b); @@ -391,7 +381,7 @@ namespace SUNDIALS # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) J = SUNMatNewEmpty(); # else - J = SUNMatNewEmpty(ida_ctx); + J = SUNMatNewEmpty(ida_ctx); # endif J->content = this; @@ -424,22 +414,16 @@ namespace SUNDIALS // calling IDASetLinearSolver status = IDASetJacFn( ida_mem, - []( -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype tt, - sunrealtype cj, -# else - realtype tt, - realtype cj, -# endif - N_Vector yy, - N_Vector yp, - N_Vector /* residual */, - SUNMatrix /* ignored */, - void *user_data, - N_Vector /* tmp1 */, - N_Vector /* tmp2 */, - N_Vector /* tmp3 */) -> int { + [](SUNDIALS::realtype tt, + SUNDIALS::realtype cj, + N_Vector yy, + N_Vector yp, + N_Vector /* residual */, + SUNMatrix /* ignored */, + void *user_data, + N_Vector /* tmp1 */, + N_Vector /* tmp2 */, + N_Vector /* tmp3 */) -> int { Assert(user_data != nullptr, ExcInternalError()); IDA &solver = *static_cast *>(user_data); diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index f80fe1d73c..483844536c 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -362,14 +362,9 @@ namespace SUNDIALS LS->ops->solve = [](SUNLinearSolver LS, SUNMatrix /*ignored*/, - N_Vector x, - N_Vector b, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype tol -# else - realtype tol -# endif - ) -> int { + N_Vector x, + N_Vector b, + SUNDIALS::realtype tol) -> int { // Receive the object that describes the linear solver and // unpack the pointer to the KINSOL object from which we can then // get the 'reinit' and 'solve' functions. @@ -399,7 +394,7 @@ namespace SUNDIALS # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) J = SUNMatNewEmpty(); # else - J = SUNMatNewEmpty(kinsol_ctx); + J = SUNMatNewEmpty(kinsol_ctx); # endif J->content = this; diff --git a/source/sundials/sunlinsol_wrapper.cc b/source/sundials/sunlinsol_wrapper.cc index 45101dcf5b..d3877f20ca 100644 --- a/source/sundials/sunlinsol_wrapper.cc +++ b/source/sundials/sunlinsol_wrapper.cc @@ -36,6 +36,7 @@ # endif # include +# include DEAL_II_NAMESPACE_OPEN @@ -195,14 +196,9 @@ namespace SUNDIALS int arkode_linsol_solve(SUNLinearSolver LS, SUNMatrix /*ignored*/, - N_Vector x, - N_Vector b, -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype tol -# else - realtype tol -# endif - ) + N_Vector x, + N_Vector b, + SUNDIALS::realtype tol) { LinearSolverContent *content = access_content(LS); diff --git a/tests/sundials/arkode_06.cc b/tests/sundials/arkode_06.cc index 701cb3829b..ab54be76c3 100644 --- a/tests/sundials/arkode_06.cc +++ b/tests/sundials/arkode_06.cc @@ -95,17 +95,11 @@ main() }; - ode.jacobian_times_setup = [&]( -#if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -#else - realtype t, -#endif - const VectorType &y, - const VectorType &fy) { - J = 0; - J(2, 2) = -1.0 / eps; - }; + ode.jacobian_times_setup = + [&](SUNDIALS::realtype t, const VectorType &y, const VectorType &fy) { + J = 0; + J(2, 2) = -1.0 / eps; + }; ode.jacobian_times_vector = [&](const VectorType &v, VectorType &Jv, diff --git a/tests/sundials/arkode_07.cc b/tests/sundials/arkode_07.cc index 1e715c40f8..8d9745bd41 100644 --- a/tests/sundials/arkode_07.cc +++ b/tests/sundials/arkode_07.cc @@ -95,17 +95,11 @@ main() }; - ode.jacobian_times_setup = [&]( -#if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - sunrealtype t, -#else - realtype t, -#endif - const VectorType &y, - const VectorType &fy) { - J = 0; - J(2, 2) = -1.0 / eps; - }; + ode.jacobian_times_setup = + [&](SUNDIALS::realtype t, const VectorType &y, const VectorType &fy) { + J = 0; + J(2, 2) = -1.0 / eps; + }; ode.jacobian_times_vector = [&](const VectorType &v, VectorType &Jv,