From f2211b37bd545c94606928da3e70fd3b3f6ce2c1 Mon Sep 17 00:00:00 2001 From: David Wells Date: Mon, 23 May 2022 09:59:33 -0400 Subject: [PATCH] Require SUNDIALS 5.4 or later. --- cmake/configure/configure_50_sundials.cmake | 4 +- include/deal.II/sundials/arkode.h | 213 ------------ include/deal.II/sundials/ida.h | 7 +- include/deal.II/sundials/kinsol.h | 3 - include/deal.II/sundials/n_vector.templates.h | 6 - include/deal.II/sundials/sundials_backport.h | 241 -------------- include/deal.II/sundials/sunlinsol_newempty.h | 104 ------ include/deal.II/sundials/sunlinsol_wrapper.h | 16 +- source/sundials/arkode.cc | 306 ++---------------- source/sundials/ida.cc | 105 +----- source/sundials/kinsol.cc | 89 +---- source/sundials/sunlinsol_wrapper.cc | 76 ++--- 12 files changed, 94 insertions(+), 1076 deletions(-) delete mode 100644 include/deal.II/sundials/sundials_backport.h delete mode 100644 include/deal.II/sundials/sunlinsol_newempty.h diff --git a/cmake/configure/configure_50_sundials.cmake b/cmake/configure/configure_50_sundials.cmake index b8042db44f..fa22c76ef9 100644 --- a/cmake/configure/configure_50_sundials.cmake +++ b/cmake/configure/configure_50_sundials.cmake @@ -24,9 +24,9 @@ MACRO(FEATURE_SUNDIALS_FIND_EXTERNAL var) SET(${var} TRUE) # - # We require at least sundials 3.0.0 + # We require at least sundials 5.4.0 # - SET(_version_required 3.0.0) + SET(_version_required 5.4.0) IF(SUNDIALS_VERSION VERSION_LESS ${_version_required}) MESSAGE(STATUS "Could not find a sufficient Sundials installation: " "deal.II requires at least version ${_version_required}, " diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index 0b2bd9117c..1d8fc96e86 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -34,9 +34,6 @@ # include # include -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -# include -# endif # include # ifdef DEAL_II_WITH_MPI # include @@ -621,209 +618,6 @@ namespace SUNDIALS std::function implicit_function; -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - /** - * A function object that users may supply and that is intended to - * prepare the linear solver for subsequent calls to - * solve_jacobian_system(). - * - * Make sure that after a call to this function, we know how to compute - * solutions of systems $A x = b$, where $A$ is some approximation to the - * Newton matrix, $M - \gamma \partial f_I/\partial y$. This function is - * optional. If the user does not provide it, then solve_jacobian_system() - * is assumed to also perform the setup internally. - * - * The setup_jacobian() function may call a user-supplied function to - * compute needed data related to the Jacobian matrix. Alternatively, it may - * choose to retrieve and use stored values of this data. In either case, - * setup_jacobian() may also preprocess that data as needed for - * solve_jacobian_system(), which may involve calling a generic function - * (such as for LU factorization). - * - * This data may be intended either for direct use (in a direct linear - * solver) or for use in a preconditioner (in a preconditioned iterative - * linear solver). The setup_jacobian() function is not called at every - * stage solve (or even every time step), but only as frequently as the - * solver determines that it is appropriate to perform the setup task. In - * this way, Jacobian-related data generated by setup_jacobian() is - * expected to be used over a number of time steps. - * - * If the user uses a matrix based computation of the Jacobian, then this - * is the right place where an assembly routine shoulde be called to - * assemble both a matrix and a preconditioner for the Jacobian system. - * Subsequent calls (possibly more than one) to solve_jacobian_system() can - * assume that this function has been called at least once. - * - * Notice that no assumption is made by this interface on what the user - * should do in this function. ARKode only assumes that after a call to - * setup_jacobian() it is possible to call solve_jacobian_system(), to - * obtain a solution $x$ to the system $J x = b$. If this function is not - * provided, then it is never called. - * - * Arguments to the function are - * - * @param[in] t the current time - * @param[in] gamma the current factor to use in the jacobian computation - * @param[in] ypred is the predicted $y$ vector for the current ARKode - * internal step - * @param[in] fpred is the value of the implicit right-hand side at ypred, - * $f_I (t_n, ypred)$. - * - * @param[in] convfail Input flag used to indicate any problem that - * occurred during the solution of the nonlinear equation on the current - * time step for which the linear solver is being used. This flag can be - * used to help decide whether the Jacobian data kept by a linear solver - * needs to be updated or not. Its possible values are: - * - * - ARK_NO_FAILURES: this value is passed if either this is the first - * call for this step, or the local error test failed on the previous - * attempt at this step (but the Newton iteration converged). - * - * - ARK_FAIL_BAD_J: this value is passed if (a) the previous Newton - * corrector iteration did not converge and the linear solver's setup - * function indicated that its Jacobian-related data is not current, or - * (b) during the previous Newton corrector iteration, the linear solver's - * solve function failed in a recoverable manner and the linear solver's - * setup function indicated that its Jacobian-related data is not - * current. - * - * - ARK_FAIL_OTHER: this value is passed if during the current internal - * step try, the previous Newton iteration failed to converge even - * though the linear solver was using current Jacobian-related data. - * - * @param[out] j_is_current: a boolean to be filled in by setup_jacobian(). - * The value should be set to `true` if the Jacobian data is current after - * the call, and should be set to `false` if its Jacobian data is not - * current. If setup_jacobian() calls for re-evaluation of Jacobian data - * (based on convfail and ARKode state data), then it should set - * `j_is_current` to `true` unconditionally, otherwise an infinite loop can - * result. - * - * This function should return: - * - 0: Success - * - >0: Recoverable error (ARKodeReinit will be called if this happens, and - * then last function will be attempted again - * - <0: Unrecoverable error the computation will be aborted and an - * assertion will be thrown. - */ - std::function - setup_jacobian; - - /** - * A function object that users may supply and that is intended to solve - * the Jacobian linear system. This function will be called by ARKode - * (possibly several times) after setup_jacobian() has been called at least - * once. ARKode tries to do its best to call setup_jacobian() the minimum - * amount of times. If convergence can be achieved without updating the - * Jacobian, then ARKode does not call setup_jacobian() again. If, on the - * contrary, internal ARKode convergence tests fail, then ARKode calls - * again setup_jacobian() with updated vectors and coefficients so that - * successive calls to solve_jacobian_systems() lead to better convergence - * in the Newton process. - * - * If you do not specify a solve_jacobian_system() function, then a fixed - * point iteration is used instead of a Newton method. Notice that this may - * not converge, or may converge very slowly. - * - * The jacobian $J$ should be (an approximation of) the system Jacobian - * \f[ - * J = M - \gamma \frac{\partial f_I}{\partial y} - * \f] - * evaluated at `t`, `ycur`. `fcur` is $f_I(t,ycur)$. - * - * A call to this function should store in `dst` the result of $J^{-1}$ - * applied to `src`, i.e., `J*dst = src`. It is the users responsibility to - * set up proper solvers and preconditioners inside this function. - * - * - * Arguments to the function are - * - * @param[in] t the current time - * @param[in] gamma the current factor to use in the jacobian computation - * @param[in] ycur is the current $y$ vector for the current ARKode - * internal step - * @param[in] fcur is the current value of the implicit right-hand side at - * ycur, $f_I (t_n, ypred)$. - * - * - * This function should return: - * - 0: Success - * - >0: Recoverable error (ARKodeReinit will be called if this happens, and - * then last function will be attempted again - * - <0: Unrecoverable error the computation will be aborted and an - * assertion will be thrown. - */ - std::function - solve_jacobian_system; - - - /** - * A function object that users may supply and that is intended to set up - * the mass matrix. This function is called by ARKode any time a mass - * matrix update is required. The user should compute the mass matrix (or - * update all the variables that allow the application of the mass matrix). - * This function is called by ARKode once, before any call to - * solve_mass_system(). - * - * ARKode supports the case where the mass matrix may depend on time, but - * not the case where the mass matrix depends on the solution itself. - * - * If the user does not provide a solve_mass_matrix() function, then the - * identity is used. If the setup_mass() function is not provided, then - * solve_mass_system() should do all the work by itself. - * - * If the user uses a matrix based computation of the mass matrix, then - * this is the right place where an assembly routine shoulde be called to - * assemble both a matrix and a preconditioner. Subsequent calls (possibly - * more than one) to solve_mass_system() can assume that this function - * has been called at least once. - * - * Notice that no assumption is made by this interface on what the user - * should do in this function. ARKode only assumes that after a call to - * setup_mass() it is possible to call solve_mass_system(), to - * obtain a solution $x$ to the system $M x = b$. - * - * This function should return: - * - 0: Success - * - >0: Recoverable error (ARKodeReinit will be called if this happens, and - * then last function will be attempted again - * - <0: Unrecoverable error the computation will be aborted and an - * assertion will be thrown. - */ - std::function setup_mass; - - /** - * A function object that users may supply and that is intended to solve - * the mass matrix linear system. This function will be called by ARKode - * (possibly several times) after setup_mass() has been called at least - * once. ARKode tries to do its best to call setup_mass() the minimum - * amount of times. - * - * A call to this function should store in `dst` the result of $M^{-1}$ - * applied to `src`, i.e., `M*dst = src`. It is the users responsibility to - * set up proper solvers and preconditioners inside this function. - * - * This function should return: - * - 0: Success - * - >0: Recoverable error (ARKodeReinit will be called if this happens, and - * then last function will be attempted again - * - <0: Unrecoverable error the computation will be aborted and an - * assertion will be thrown. - */ - std::function - solve_mass_system; -# else - /** * A function object that users may supply and that is intended to compute * the product of the mass matrix with a given vector `v`. This function @@ -1135,7 +929,6 @@ namespace SUNDIALS * assertion will be thrown. */ std::function mass_preconditioner_setup; -# endif /** * A function object that users may supply and that is intended to @@ -1227,8 +1020,6 @@ namespace SUNDIALS dealii::DiscreteTime &time, const bool do_reset); -# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) - /** * Set up the (non)linear solver and preconditioners in the ARKODE memory * object based on the user-specified functions. @@ -1247,8 +1038,6 @@ namespace SUNDIALS void setup_mass_solver(const VectorType &solution); -# endif - /** * This function is executed at construction time to set the * std::function above to trigger an assert if they are not @@ -1286,10 +1075,8 @@ namespace SUNDIALS */ double last_end_time; -# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) std::unique_ptr> linear_solver; std::unique_ptr> mass_solver; -# endif # ifdef DEAL_II_WITH_PETSC # ifdef PETSC_USE_COMPLEX diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index 82161698d6..7e374851da 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -37,17 +37,12 @@ # include # endif -# include -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) -# include -# include -# include -# endif # include # include # include +# include # include # include diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index ff22bb541d..5b2c50cd19 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -33,9 +33,6 @@ # include # include -# if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0) -# include -# endif # include # include # include diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index a030c20d73..436961963f 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -34,10 +34,6 @@ # include # include -# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) -# include -# endif - DEAL_II_NAMESPACE_OPEN namespace SUNDIALS @@ -1000,11 +996,9 @@ namespace SUNDIALS v->ops->nvcloneempty = &NVectorOperations::clone_empty; v->ops->nvdestroy = &NVectorOperations::destroy; // v->ops->nvspace = undef; -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0) v->ops->nvgetcommunicator = &NVectorOperations::get_communicator_as_void_ptr; v->ops->nvgetlength = &NVectorOperations::get_global_length; -# endif /* standard vector operations */ v->ops->nvlinearsum = &NVectorOperations::linear_sum; diff --git a/include/deal.II/sundials/sundials_backport.h b/include/deal.II/sundials/sundials_backport.h deleted file mode 100644 index 39d528df6c..0000000000 --- a/include/deal.II/sundials/sundials_backport.h +++ /dev/null @@ -1,241 +0,0 @@ -//----------------------------------------------------------- -// -// Copyright (C) 2020 - 2021 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. -// -//----------------------------------------------------------- - -/* - * The functions in this file are based on an implementation distributed within - * the SUNDIALS package, see the license here: - * https://computing.llnl.gov/projects/sundials/license. - * ----------------------------------------------------------------- - * Programmer(s): Daniel Reynolds @ SMU - * David J. Gardner, Carol S. Woodward, and - * Slaven Peles @ LLNL - * -----------------------------------------------------------------*/ - -#ifndef dealii_sundials_backport_h -#define dealii_sundials_backport_h - -#include - -#ifdef DEAL_II_WITH_SUNDIALS - -# include - -DEAL_II_NAMESPACE_OPEN -namespace SUNDIALS -{ - namespace internal - { - N_Vector - N_VNewEmpty() - { - N_Vector v = new _generic_N_Vector; - N_Vector_Ops ops = new _generic_N_Vector_Ops; - - /* initialize operations to nullptr */ - - /* constructors, destructors, and utility operations */ - ops->nvgetvectorid = nullptr; - ops->nvclone = nullptr; - ops->nvcloneempty = nullptr; - ops->nvdestroy = nullptr; - ops->nvspace = nullptr; - ops->nvgetarraypointer = nullptr; - ops->nvsetarraypointer = nullptr; -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0) - ops->nvgetcommunicator = nullptr; - ops->nvgetlength = nullptr; -# endif - - /* standard vector operations */ - ops->nvlinearsum = nullptr; - ops->nvconst = nullptr; - ops->nvprod = nullptr; - ops->nvdiv = nullptr; - ops->nvscale = nullptr; - ops->nvabs = nullptr; - ops->nvinv = nullptr; - ops->nvaddconst = nullptr; - ops->nvdotprod = nullptr; - ops->nvmaxnorm = nullptr; - ops->nvwrmsnormmask = nullptr; - ops->nvwrmsnorm = nullptr; - ops->nvmin = nullptr; - ops->nvwl2norm = nullptr; - ops->nvl1norm = nullptr; - ops->nvcompare = nullptr; - ops->nvinvtest = nullptr; - ops->nvconstrmask = nullptr; - ops->nvminquotient = nullptr; - -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0) - /* fused vector operations (optional) */ - ops->nvlinearcombination = nullptr; - ops->nvscaleaddmulti = nullptr; - ops->nvdotprodmulti = nullptr; - - /* vector array operations (optional) */ - ops->nvlinearsumvectorarray = nullptr; - ops->nvscalevectorarray = nullptr; - ops->nvconstvectorarray = nullptr; - ops->nvwrmsnormvectorarray = nullptr; - ops->nvwrmsnormmaskvectorarray = nullptr; - ops->nvscaleaddmultivectorarray = nullptr; - ops->nvlinearcombinationvectorarray = nullptr; - - /* local reduction operations (optional) */ - ops->nvdotprodlocal = nullptr; - ops->nvmaxnormlocal = nullptr; - ops->nvminlocal = nullptr; - ops->nvl1normlocal = nullptr; - ops->nvinvtestlocal = nullptr; - ops->nvconstrmasklocal = nullptr; - ops->nvminquotientlocal = nullptr; - ops->nvwsqrsumlocal = nullptr; - ops->nvwsqrsummasklocal = nullptr; - -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 4, 0) - /* XBraid interface operations */ - ops->nvbufsize = nullptr; - ops->nvbufpack = nullptr; - ops->nvbufunpack = nullptr; -# endif - -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 3, 0) - /* debugging functions (called when SUNDIALS_DEBUG_PRINTVEC is defined) */ - ops->nvprint = nullptr; - ops->nvprintfile = nullptr; -# endif -# endif - - /* attach ops and initialize content to nullptr */ - v->ops = ops; - v->content = nullptr; - - return v; - } - - - - void - N_VFreeEmpty(N_Vector v) - { - if (v == nullptr) - return; - - /* free non-nullptr ops structure */ - if (v->ops) - delete v->ops; - v->ops = nullptr; - - /* free overall N_Vector object and return */ - delete v; - } - - - - int - N_VCopyOps(N_Vector w, N_Vector v) - { - /* Check that ops structures exist */ - if (w == nullptr || v == nullptr) - return (-1); - if (w->ops == nullptr || v->ops == nullptr) - return (-1); - - /* Copy ops from w to v */ - - /* constructors, destructors, and utility operations */ - v->ops->nvgetvectorid = w->ops->nvgetvectorid; - v->ops->nvclone = w->ops->nvclone; - v->ops->nvcloneempty = w->ops->nvcloneempty; - v->ops->nvdestroy = w->ops->nvdestroy; - v->ops->nvspace = w->ops->nvspace; - v->ops->nvgetarraypointer = w->ops->nvgetarraypointer; - v->ops->nvsetarraypointer = w->ops->nvsetarraypointer; -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0) - v->ops->nvgetcommunicator = w->ops->nvgetcommunicator; - v->ops->nvgetlength = w->ops->nvgetlength; -# endif - - /* standard vector operations */ - v->ops->nvlinearsum = w->ops->nvlinearsum; - v->ops->nvconst = w->ops->nvconst; - v->ops->nvprod = w->ops->nvprod; - v->ops->nvdiv = w->ops->nvdiv; - v->ops->nvscale = w->ops->nvscale; - v->ops->nvabs = w->ops->nvabs; - v->ops->nvinv = w->ops->nvinv; - v->ops->nvaddconst = w->ops->nvaddconst; - v->ops->nvdotprod = w->ops->nvdotprod; - v->ops->nvmaxnorm = w->ops->nvmaxnorm; - v->ops->nvwrmsnormmask = w->ops->nvwrmsnormmask; - v->ops->nvwrmsnorm = w->ops->nvwrmsnorm; - v->ops->nvmin = w->ops->nvmin; - v->ops->nvwl2norm = w->ops->nvwl2norm; - v->ops->nvl1norm = w->ops->nvl1norm; - v->ops->nvcompare = w->ops->nvcompare; - v->ops->nvinvtest = w->ops->nvinvtest; - v->ops->nvconstrmask = w->ops->nvconstrmask; - v->ops->nvminquotient = w->ops->nvminquotient; - -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0) - /* fused vector operations */ - v->ops->nvlinearcombination = w->ops->nvlinearcombination; - v->ops->nvscaleaddmulti = w->ops->nvscaleaddmulti; - v->ops->nvdotprodmulti = w->ops->nvdotprodmulti; - - /* vector array operations */ - v->ops->nvlinearsumvectorarray = w->ops->nvlinearsumvectorarray; - v->ops->nvscalevectorarray = w->ops->nvscalevectorarray; - v->ops->nvconstvectorarray = w->ops->nvconstvectorarray; - v->ops->nvwrmsnormvectorarray = w->ops->nvwrmsnormvectorarray; - v->ops->nvwrmsnormmaskvectorarray = w->ops->nvwrmsnormmaskvectorarray; - v->ops->nvscaleaddmultivectorarray = w->ops->nvscaleaddmultivectorarray; - v->ops->nvlinearcombinationvectorarray = - w->ops->nvlinearcombinationvectorarray; - - /* local reduction operations */ - v->ops->nvdotprodlocal = w->ops->nvdotprodlocal; - v->ops->nvmaxnormlocal = w->ops->nvmaxnormlocal; - v->ops->nvminlocal = w->ops->nvminlocal; - v->ops->nvl1normlocal = w->ops->nvl1normlocal; - v->ops->nvinvtestlocal = w->ops->nvinvtestlocal; - v->ops->nvconstrmasklocal = w->ops->nvconstrmasklocal; - v->ops->nvminquotientlocal = w->ops->nvminquotientlocal; - v->ops->nvwsqrsumlocal = w->ops->nvwsqrsumlocal; - v->ops->nvwsqrsummasklocal = w->ops->nvwsqrsummasklocal; - -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 4, 0) - /* XBraid interface operations */ - v->ops->nvbufsize = w->ops->nvbufsize; - v->ops->nvbufpack = w->ops->nvbufpack; - v->ops->nvbufunpack = w->ops->nvbufunpack; -# endif - -# if DEAL_II_SUNDIALS_VERSION_GTE(5, 3, 0) - /* debugging functions (called when SUNDIALS_DEBUG_PRINTVEC is defined) */ - v->ops->nvprint = w->ops->nvprint; - v->ops->nvprintfile = w->ops->nvprintfile; -# endif -# endif - - return (0); - } - } // namespace internal -} // namespace SUNDIALS -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_WITH_SUNDIALS -#endif // dealii_sundials_sunlinsol_newempty_h diff --git a/include/deal.II/sundials/sunlinsol_newempty.h b/include/deal.II/sundials/sunlinsol_newempty.h deleted file mode 100644 index 66b8bc73e0..0000000000 --- a/include/deal.II/sundials/sunlinsol_newempty.h +++ /dev/null @@ -1,104 +0,0 @@ -//----------------------------------------------------------- -// -// Copyright (C) 2020 - 2021 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. -// -//----------------------------------------------------------- - -/* - * The functions in this file are based on an implementation distributed within - * the SUNDIALS package, see the license here: - * https://computing.llnl.gov/projects/sundials/license. - * ----------------------------------------------------------------- - * Programmer(s): Daniel Reynolds @ SMU - * David J. Gardner, Carol S. Woodward, and - * Slaven Peles @ LLNL - * -----------------------------------------------------------------*/ - -#ifndef dealii_sundials_sunlinsol_newempty_h -#define dealii_sundials_sunlinsol_newempty_h - -#include - -#ifdef DEAL_II_WITH_SUNDIALS - -# include - -DEAL_II_NAMESPACE_OPEN -namespace SUNDIALS -{ - namespace internal - { - /** - * Create a new SUNLinearSolver structure without any content and - * operations set to `nullptr`. - */ - inline SUNLinearSolver - SUNLinSolNewEmpty() - { - /* create linear solver object */ - SUNLinearSolver LS = new _generic_SUNLinearSolver; - - /* create linear solver ops structure */ - SUNLinearSolver_Ops ops = new _generic_SUNLinearSolver_Ops; - - /* initialize operations to nullptr */ - ops->gettype = nullptr; - ops->setatimes = nullptr; - ops->setpreconditioner = nullptr; - ops->setscalingvectors = nullptr; - ops->initialize = nullptr; - ops->setup = nullptr; - ops->solve = nullptr; - ops->numiters = nullptr; - ops->resnorm = nullptr; - ops->resid = nullptr; - ops->lastflag = nullptr; - ops->space = nullptr; - ops->free = nullptr; - - /* attach ops and initialize content to nullptr */ - LS->ops = ops; - LS->content = nullptr; - - return (LS); - } - - /** - * Free the memory associated with @p solver which was previously allocated - * with a call to SUNLinSolNewEmpty(). - * - * @note A call to this function does not deallocate the `content` field. - * - * @param solver The solver memory to free - */ - inline void - SUNLinSolFreeEmpty(SUNLinearSolver solver) - { - if (solver == nullptr) - return; - - /* free non-nullptr ops structure */ - if (solver->ops) - delete solver->ops; - solver->ops = nullptr; - - /* free overall linear solver object */ - delete solver; - return; - } - - } // namespace internal -} // namespace SUNDIALS -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_WITH_SUNDIALS -#endif // dealii_sundials_sunlinsol_newempty_h diff --git a/include/deal.II/sundials/sunlinsol_wrapper.h b/include/deal.II/sundials/sunlinsol_wrapper.h index d5f90e3dba..df9f8aa7d6 100644 --- a/include/deal.II/sundials/sunlinsol_wrapper.h +++ b/include/deal.II/sundials/sunlinsol_wrapper.h @@ -19,25 +19,24 @@ #include #ifdef DEAL_II_WITH_SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) -# include +# include -# include -# include +# include +# include DEAL_II_NAMESPACE_OPEN namespace SUNDIALS { -# ifndef DOXYGEN +# ifndef DOXYGEN // forward declarations namespace internal { template struct LinearSolverContent; } -# endif +# endif /** * A linear operator that wraps SUNDIALS functionality. @@ -167,10 +166,10 @@ namespace SUNDIALS { public: explicit LinearSolverWrapper(LinearSolveFunction lsolve -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , SUNContext &linsol_ctx -# endif +# endif ); ~LinearSolverWrapper(); @@ -189,6 +188,5 @@ namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE -# endif #endif #endif diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index fda97a4a05..6c48c6adec 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -38,25 +38,12 @@ # include # include -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -# include -# include -# else -# include -# include -# include -# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) -# include -# endif -# endif +# include +# include +# include # include -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -// Make sure we know how to call sundials own ARKode() function -const auto &SundialsARKode = ARKode; -# endif - DEAL_II_NAMESPACE_OPEN namespace SUNDIALS @@ -101,105 +88,6 @@ namespace SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - template - int - setup_jacobian_callback(ARKodeMem arkode_mem, - int convfail, - N_Vector ypred, - N_Vector fpred, - booleantype *jcurPtr, - N_Vector, - N_Vector, - N_Vector) - { - Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); - ARKode &solver = - *static_cast *>(arkode_mem->ark_user_data); - - auto *src_ypred = internal::unwrap_nvector_const(ypred); - auto *src_fpred = internal::unwrap_nvector_const(fpred); - // avoid reinterpret_cast - bool jcurPtr_tmp = false; - int err = solver.setup_jacobian(convfail, - arkode_mem->ark_tn, - arkode_mem->ark_gamma, - *src_ypred, - *src_fpred, - jcurPtr_tmp); - *jcurPtr = jcurPtr_tmp ? SUNTRUE : SUNFALSE; - - return err; - } - - - - template - int - solve_with_jacobian_callback(ARKodeMem arkode_mem, - N_Vector b, - N_Vector ycur, - N_Vector fcur) - { - Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); - ARKode &solver = - *static_cast *>(arkode_mem->ark_user_data); - - auto *dst = internal::unwrap_nvector(b); - auto *src_ycur = internal::unwrap_nvector_const(ycur); - auto *src_fcur = internal::unwrap_nvector_const(fcur); - - // make a temporary copy to work on in the user call - VectorType src = *dst; - - int err = solver.solve_jacobian_system(arkode_mem->ark_tn, - arkode_mem->ark_gamma, - *src_ycur, - *src_fcur, - src, - *dst); - - - return err; - } - - - - template - int - setup_mass_matrix_callback(ARKodeMem arkode_mem, - N_Vector, - N_Vector, - N_Vector) - { - Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); - ARKode &solver = - *static_cast *>(arkode_mem->ark_user_data); - int err = solver.setup_mass(arkode_mem->ark_tn); - return err; - } - - - - template - int - solve_with_mass_matrix_callback(ARKodeMem arkode_mem, N_Vector b) - { - Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); - ARKode &solver = - *static_cast *>(arkode_mem->ark_user_data); - - auto *dst = internal::unwrap_nvector(b); - - // make a temporary copy to work on in the user call - VectorType src = *dst; - - int err = solver.solve_mass_system(src, *dst); - - return err; - } -# else - template int jacobian_times_vector_callback(N_Vector v, @@ -358,8 +246,6 @@ namespace SUNDIALS return solver.mass_preconditioner_setup(t); } -# endif - } // namespace @@ -383,11 +269,7 @@ namespace SUNDIALS { if (arkode_mem) { -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - ARKodeFree(&arkode_mem); -# else ARKStepFree(&arkode_mem); -# endif # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) const int status = SUNContext_Free(&arkode_ctx); @@ -463,20 +345,12 @@ namespace SUNDIALS time.set_desired_next_step_size(data.output_period); // Let ARKode advance time by one period: - double actual_next_time; -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - const auto status = SundialsARKode(arkode_mem, - time.get_next_time(), - solution_nvector, - &actual_next_time, - ARK_NORMAL); -# else + double actual_next_time; const auto status = ARKStepEvolve(arkode_mem, time.get_next_time(), solution_nvector, &actual_next_time, ARK_NORMAL); -# endif (void)status; AssertARKode(status); @@ -502,7 +376,6 @@ namespace SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) template void ARKode::reset(const double current_time, @@ -510,126 +383,10 @@ namespace SUNDIALS const VectorType &solution) { last_end_time = current_time; - if (arkode_mem) - ARKodeFree(&arkode_mem); - - arkode_mem = ARKodeCreate(); - int status; (void)status; - Assert(explicit_function || implicit_function, - ExcFunctionNotProvided("explicit_function || implicit_function")); - - // just a view on the memory in solution, all write operations on yy by - // ARKODE will automatically be mirrored to solution - auto initial_condition_nvector = internal::make_nvector_view(solution -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - arkode_ctx -# endif - ); - - status = ARKodeInit( - arkode_mem, - explicit_function ? &explicit_function_callback : nullptr, - implicit_function ? &implicit_function_callback : nullptr, - current_time, - initial_condition_nvector); - AssertARKode(status); - - if (get_local_tolerances) - { - const auto abs_tols = internal::make_nvector_view(get_local_tolerances() -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) - , - arkode_ctx -# endif - ); - status = - ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols); - AssertARKode(status); - } - else - { - status = ARKodeSStolerances(arkode_mem, - data.relative_tolerance, - data.absolute_tolerance); - AssertARKode(status); - } - - status = ARKodeSetInitStep(arkode_mem, current_time_step); - AssertARKode(status); - - status = ARKodeSetUserData(arkode_mem, this); - AssertARKode(status); - - status = ARKodeSetStopTime(arkode_mem, data.final_time); - AssertARKode(status); - - status = - ARKodeSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations); - AssertARKode(status); - - // Initialize solver - auto ARKode_mem = static_cast(arkode_mem); - - if (solve_jacobian_system) - { - status = ARKodeSetNewton(arkode_mem); - AssertARKode(status); - if (data.implicit_function_is_linear) - { - status = ARKodeSetLinear( - arkode_mem, data.implicit_function_is_time_independent ? 0 : 1); - AssertARKode(status); - } - - - ARKode_mem->ark_lsolve = solve_with_jacobian_callback; - if (setup_jacobian) - { - ARKode_mem->ark_lsetup = setup_jacobian_callback; - } - } - else - { - status = - ARKodeSetFixedPoint(arkode_mem, data.maximum_non_linear_iterations); - AssertARKode(status); - } - - - if (solve_mass_system) - { - ARKode_mem->ark_msolve = solve_with_mass_matrix_callback; - - if (setup_mass) - { - ARKode_mem->ark_msetup = setup_mass_matrix_callback; - } - } - - status = ARKodeSetOrder(arkode_mem, data.maximum_order); - AssertARKode(status); - - if (custom_setup) - custom_setup(arkode_mem); - } - -# else - - template - void - ARKode::reset(const double current_time, - const double current_time_step, - const VectorType &solution) - { - last_end_time = current_time; - int status; - (void)status; - -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) if (arkode_ctx) { status = SUNContext_Free(&arkode_ctx); @@ -637,7 +394,7 @@ namespace SUNDIALS } status = SUNContext_Create(&mpi_communicator, &arkode_ctx); AssertARKode(status); -# endif +# endif if (arkode_mem) { @@ -648,10 +405,10 @@ namespace SUNDIALS // just a view on the memory in solution, all write operations on yy by // ARKODE will automatically be mirrored to solution auto initial_condition_nvector = internal::make_nvector_view(solution -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); Assert(explicit_function || implicit_function, @@ -662,20 +419,20 @@ namespace SUNDIALS implicit_function ? &implicit_function_callback : nullptr, current_time, initial_condition_nvector -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); Assert(arkode_mem != nullptr, ExcInternalError()); if (get_local_tolerances) { const auto abs_tols = internal::make_nvector_view(get_local_tolerances() -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); status = ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols); @@ -732,10 +489,10 @@ namespace SUNDIALS linear_solver = std::make_unique>( solve_linearized_system -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); sun_linear_solver = *linear_solver; } @@ -744,23 +501,23 @@ namespace SUNDIALS // use default solver from SUNDIALS // TODO give user options auto y_template = internal::make_nvector_view(solution -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) sun_linear_solver = SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/); -# else +# else sun_linear_solver = SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/, arkode_ctx); -# endif +# endif } status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr); AssertARKode(status); @@ -791,22 +548,22 @@ namespace SUNDIALS else { auto y_template = internal::make_nvector_view(solution -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) SUNNonlinearSolver fixed_point_solver = SUNNonlinSol_FixedPoint(y_template, data.anderson_acceleration_subspace); -# else +# else SUNNonlinearSolver fixed_point_solver = SUNNonlinSol_FixedPoint(y_template, data.anderson_acceleration_subspace, arkode_ctx); -# endif +# endif status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver); AssertARKode(status); @@ -834,33 +591,33 @@ namespace SUNDIALS mass_solver = std::make_unique>( solve_mass -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); sun_mass_linear_solver = *mass_solver; } else { auto y_template = internal::make_nvector_view(solution -# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) , arkode_ctx -# endif +# endif ); -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) sun_mass_linear_solver = SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/); -# else +# else sun_mass_linear_solver = SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/, arkode_ctx); -# endif +# endif } booleantype mass_time_dependent = data.mass_is_time_independent ? SUNFALSE : SUNTRUE; @@ -890,7 +647,6 @@ namespace SUNDIALS } } } -# endif diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index 801e1c0380..c07aa31977 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -36,17 +36,6 @@ # include -# include -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -# ifdef DEAL_II_SUNDIALS_WITH_IDAS -# include -# else -# include -# endif -# endif -# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) -# include -# endif # include # include @@ -77,68 +66,6 @@ namespace SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - template - int - t_dae_lsetup(IDAMem IDA_mem, - N_Vector yy, - N_Vector yp, - N_Vector resp, - N_Vector tmp1, - N_Vector tmp2, - N_Vector tmp3) - { - (void)tmp1; - (void)tmp2; - (void)tmp3; - (void)resp; - IDA &solver = - *static_cast *>(IDA_mem->ida_user_data); - - auto *src_yy = internal::unwrap_nvector_const(yy); - auto *src_yp = internal::unwrap_nvector_const(yp); - - int err = solver.setup_jacobian(IDA_mem->ida_tn, - *src_yy, - *src_yp, - IDA_mem->ida_cj); - - return err; - } - - - - template - int - solve_with_jacobian_callback(IDAMem IDA_mem, - N_Vector b, - N_Vector weight, - N_Vector yy, - N_Vector yp, - N_Vector resp) - { - (void)weight; - (void)yy; - (void)yp; - (void)resp; - IDA &solver = - *static_cast *>(IDA_mem->ida_user_data); - GrowingVectorMemory mem; - - typename VectorMemory::Pointer dst(mem); - solver.reinit_vector(*dst); - - auto *src = internal::unwrap_nvector(b); - - int err = solver.solve_jacobian_system(*src, *dst); - *src = *dst; - - return err; - } - - - -# else template int setup_jacobian_callback(realtype tt, @@ -188,7 +115,6 @@ namespace SUNDIALS return err; } -# endif } // namespace @@ -315,7 +241,7 @@ namespace SUNDIALS status = SUNContext_Create(&mpi_communicator, &ida_ctx); AssertIDA(status); - ida_mem = IDACreate(ida_ctx); + ida_mem = IDACreate(ida_ctx); # endif auto yy = internal::make_nvector_view(solution @@ -391,16 +317,6 @@ namespace SUNDIALS AssertIDA(status); // Initialize solver -# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) - auto IDA_mem = static_cast(ida_mem); - - IDA_mem->ida_lsetup = t_dae_lsetup; - - if (solve_jacobian_system) - IDA_mem->ida_lsolve = solve_with_jacobian_callback; - else - AssertThrow(false, ExcFunctionNotProvided("solve_jacobian_system")); -# else SUNMatrix J = nullptr; SUNLinearSolver LS = nullptr; @@ -408,11 +324,11 @@ namespace SUNDIALS // called do not actually receive the IDAMEM object, just the LS // object, so we have to store a pointer to the current // object in the LS object -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) LS = SUNLinSolNewEmpty(); -# else - LS = SUNLinSolNewEmpty(ida_ctx); -# endif +# else + LS = SUNLinSolNewEmpty(ida_ctx); +# endif LS->content = this; @@ -459,11 +375,11 @@ namespace SUNDIALS // if we don't set it, it won't call the functions that set up // the matrix object (i.e., the argument to the 'IDASetJacFn' // function below). -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) - J = SUNMatNewEmpty(); -# else - J = SUNMatNewEmpty(ida_ctx); -# endif +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + J = SUNMatNewEmpty(); +# else + J = SUNMatNewEmpty(ida_ctx); +# endif J->content = this; J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID { @@ -495,7 +411,6 @@ namespace SUNDIALS // calling IDASetLinearSolver status = IDASetJacFn(ida_mem, &setup_jacobian_callback); AssertIDA(status); -# endif status = IDASetMaxOrd(ida_mem, data.maximum_order); AssertIDA(status); diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index 97869c426b..66b40dfc84 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -170,55 +170,6 @@ namespace SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0) - template - int - setup_jacobian_callback(KINMem kinsol_mem) - { - KINSOL &solver = - *static_cast *>(kinsol_mem->kin_user_data); - - auto *src_ycur = - internal::unwrap_nvector_const(kinsol_mem->kin_uu); - auto *src_fcur = - internal::unwrap_nvector_const(kinsol_mem->kin_fval); - - int err = solver.setup_jacobian(*src_ycur, *src_fcur); - return err; - } - - - - template - int - solve_with_jacobian_callback(KINMem kinsol_mem, - N_Vector x, - N_Vector b, - realtype *sJpnorm, - realtype *sFdotJp) - { - KINSOL &solver = - *static_cast *>(kinsol_mem->kin_user_data); - - auto *src_ycur = - internal::unwrap_nvector_const(kinsol_mem->kin_uu); - auto *src_fcur = - internal::unwrap_nvector_const(kinsol_mem->kin_fval); - auto *src = internal::unwrap_nvector_const(b); - auto *dst = internal::unwrap_nvector(x); - - int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src, *dst); - - *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale); - N_VProd(b, kinsol_mem->kin_fscale, b); - N_VProd(b, kinsol_mem->kin_fscale, b); - *sFdotJp = N_VDotProd(kinsol_mem->kin_fval, b); - - return err; - } - -# else // SUNDIALS 5.0 or later - template int setup_jacobian_callback(N_Vector u, @@ -297,7 +248,6 @@ namespace SUNDIALS return err; } } -# endif } // namespace @@ -479,38 +429,16 @@ namespace SUNDIALS solve_with_jacobian) // user assigned a function object to the solver // slot { -/* interface up to and including 4.0 */ -# if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0) - auto KIN_mem = static_cast(kinsol_mem); - // Old version only works with solve_jacobian_system - Assert(solve_jacobian_system, - ExcFunctionNotProvided("solve_jacobian_system")) - KIN_mem->kin_lsolve = solve_with_jacobian_callback; - if (setup_jacobian) // user assigned a function object to the Jacobian - // set-up slot - KIN_mem->kin_lsetup = setup_jacobian_callback; - -/* interface up to and including 4.1 */ -# elif DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) - - // deal.II does not currently have support for KINSOL in - // SUNDIALS 4.1. One could write this and update this section, - // but it does not seem worthwhile spending the time to - // interface with an old version of SUNDIAL given that the - // code below supports modern SUNDIAL versions just fine. - Assert(false, ExcNotImplemented()); - -# else /* interface starting with SUNDIALS 5.0 */ // Set the operations we care for in the sun_linear_solver object // and attach it to the KINSOL object. The functions that will get // called do not actually receive the KINSOL object, just the LS // object, so we have to store a pointer to the current // object in the LS object -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) - LS = SUNLinSolNewEmpty(); -# else +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + LS = SUNLinSolNewEmpty(); +# else LS = SUNLinSolNewEmpty(kinsol_ctx); -# endif +# endif LS->content = this; LS->ops->gettype = @@ -540,11 +468,11 @@ namespace SUNDIALS // if we don't set it, it won't call the functions that set up // the matrix object (i.e., the argument to the 'KINSetJacFn' // function below). -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) - J = SUNMatNewEmpty(); -# else +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) + J = SUNMatNewEmpty(); +# else J = SUNMatNewEmpty(kinsol_ctx); -# endif +# endif J->content = this; J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID { @@ -578,7 +506,6 @@ namespace SUNDIALS }; status = KINSetJacFn(kinsol_mem, &setup_jacobian_callback); AssertKINSOL(status); -# endif } // call to KINSol diff --git a/source/sundials/sunlinsol_wrapper.cc b/source/sundials/sunlinsol_wrapper.cc index 9062d2f828..49f70c810e 100644 --- a/source/sundials/sunlinsol_wrapper.cc +++ b/source/sundials/sunlinsol_wrapper.cc @@ -18,28 +18,23 @@ #include #ifdef DEAL_II_WITH_SUNDIALS -# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) -# include +# include -# include -# include -# include -# include -# ifdef DEAL_II_WITH_TRILINOS -# include -# include -# endif -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif - -# include +# include +# include +# include +# include +# ifdef DEAL_II_WITH_TRILINOS +# include +# include +# endif +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif -# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0) -# include -# endif +# include DEAL_II_NAMESPACE_OPEN @@ -53,8 +48,8 @@ namespace SUNDIALS << " functions returned a negative error code: " << arg1 << ". Please consult SUNDIALS manual."); -# define AssertSundialsSolver(code) \ - Assert(code >= 0, ExcSundialsSolverError(code)) +# define AssertSundialsSolver(code) \ + Assert(code >= 0, ExcSundialsSolverError(code)) namespace internal { @@ -177,18 +172,18 @@ namespace SUNDIALS template internal::LinearSolverWrapper::LinearSolverWrapper( LinearSolveFunction lsolve -# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if !DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) , SUNContext &linsol_ctx -# endif +# endif ) : content(std::make_unique>()) { -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) sun_linear_solver = SUNLinSolNewEmpty(); -# else +# else sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx); -# endif +# endif sun_linear_solver->ops->gettype = arkode_linsol_get_type; sun_linear_solver->ops->solve = arkode_linsol_solve; @@ -240,13 +235,13 @@ namespace SUNDIALS SundialsOperator::vmult(VectorType & dst, const VectorType &src) const { -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) auto sun_dst = internal::make_nvector_view(dst); auto sun_src = internal::make_nvector_view(src); int status = a_times_fn(A_data, sun_src, sun_dst); (void)status; AssertSundialsSolver(status); -# else +# else // We don't currently know how to implement this: the // internal::make_nvector_view() function called above requires a // SUNContext object, but we don't actually have access to such an @@ -254,7 +249,7 @@ namespace SUNDIALS (void)dst; (void)src; Assert(false, ExcNotImplemented()); -# endif +# endif } @@ -283,7 +278,7 @@ namespace SUNDIALS return; } -# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0) auto sun_dst = internal::make_nvector_view(dst); auto sun_src = internal::make_nvector_view(src); // for custom preconditioners no distinction between left and right @@ -292,7 +287,7 @@ namespace SUNDIALS p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/); (void)status; AssertSundialsSolver(status); -# else +# else // We don't currently know how to implement this: the // internal::make_nvector_view() function called above requires a // SUNContext object, but we don't actually have access to such an @@ -300,7 +295,7 @@ namespace SUNDIALS (void)dst; (void)src; Assert(false, ExcNotImplemented()); -# endif +# endif } template struct SundialsOperator>; @@ -323,8 +318,8 @@ namespace SUNDIALS template class internal::LinearSolverWrapper< LinearAlgebra::distributed::BlockVector>; -# ifdef DEAL_II_WITH_MPI -# ifdef DEAL_II_WITH_TRILINOS +# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_TRILINOS template struct SundialsOperator; template struct SundialsOperator; @@ -334,10 +329,10 @@ namespace SUNDIALS template class internal::LinearSolverWrapper; template class internal::LinearSolverWrapper< TrilinosWrappers::MPI::BlockVector>; -# endif // DEAL_II_WITH_TRILINOS +# endif // DEAL_II_WITH_TRILINOS -# ifdef DEAL_II_WITH_PETSC -# ifndef PETSC_USE_COMPLEX +# ifdef DEAL_II_WITH_PETSC +# ifndef PETSC_USE_COMPLEX template struct SundialsOperator; template struct SundialsOperator; @@ -347,12 +342,11 @@ namespace SUNDIALS template class internal::LinearSolverWrapper; template class internal::LinearSolverWrapper; -# endif // PETSC_USE_COMPLEX -# endif // DEAL_II_WITH_PETSC +# endif // PETSC_USE_COMPLEX +# endif // DEAL_II_WITH_PETSC -# endif // DEAL_II_WITH_MPI +# endif // DEAL_II_WITH_MPI } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE -# endif #endif -- 2.39.5