]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Require SUNDIALS 5.4 or later.
authorDavid Wells <drwells@email.unc.edu>
Mon, 23 May 2022 13:59:33 +0000 (09:59 -0400)
committerDavid Wells <drwells@email.unc.edu>
Fri, 27 May 2022 12:19:48 +0000 (08:19 -0400)
12 files changed:
cmake/configure/configure_50_sundials.cmake
include/deal.II/sundials/arkode.h
include/deal.II/sundials/ida.h
include/deal.II/sundials/kinsol.h
include/deal.II/sundials/n_vector.templates.h
include/deal.II/sundials/sundials_backport.h [deleted file]
include/deal.II/sundials/sunlinsol_newempty.h [deleted file]
include/deal.II/sundials/sunlinsol_wrapper.h
source/sundials/arkode.cc
source/sundials/ida.cc
source/sundials/kinsol.cc
source/sundials/sunlinsol_wrapper.cc

index b8042db44f7ec5795909b082efdab24a42101322..fa22c76ef9a1b427e4baee48015c452796c8b395 100644 (file)
@@ -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}, "
index 0b2bd9117c19601745bacc6a72c048262d7152a9..1d8fc96e86c40853fee460654f2fcc1971eaebd6 100644 (file)
@@ -34,9 +34,6 @@
 #  include <deal.II/lac/vector_memory.h>
 
 #  include <arkode/arkode.h>
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-#    include <arkode/arkode_impl.h>
-#  endif
 #  include <nvector/nvector_serial.h>
 #  ifdef DEAL_II_WITH_MPI
 #    include <nvector/nvector_parallel.h>
@@ -621,209 +618,6 @@ namespace SUNDIALS
     std::function<int(const double t, const VectorType &y, VectorType &res)>
       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<int(const int         convfail,
-                      const double      t,
-                      const double      gamma,
-                      const VectorType &ypred,
-                      const VectorType &fpred,
-                      bool &            j_is_current)>
-      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<int(const double      t,
-                      const double      gamma,
-                      const VectorType &ycur,
-                      const VectorType &fcur,
-                      const VectorType &rhs,
-                      VectorType &      dst)>
-      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<int(const double t)> 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<int(const VectorType &rhs, VectorType &dst)>
-      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<int(double t)> 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<internal::LinearSolverWrapper<VectorType>> linear_solver;
     std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
-#  endif
 
 #  ifdef DEAL_II_WITH_PETSC
 #    ifdef PETSC_USE_COMPLEX
index 82161698d6b016a5e5c1359e6af32fa0202f42c3..7e374851da6348226fe3d195f4630e53274c2ea2 100644 (file)
 #    include <ida/ida.h>
 #  endif
 
-#  include <sundials/sundials_config.h>
-#  if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
-#    include <ida/ida_spbcgs.h>
-#    include <ida/ida_spgmr.h>
-#    include <ida/ida_sptfqmr.h>
-#  endif
 #  include <deal.II/sundials/sunlinsol_wrapper.h>
 
 #  include <boost/signals2.hpp>
 
 #  include <nvector/nvector_serial.h>
+#  include <sundials/sundials_config.h>
 #  include <sundials/sundials_math.h>
 #  include <sundials/sundials_types.h>
 
index ff22bb541d9c649d2ece48fa469f80cf80ecdcb0..5b2c50cd194c25df9252a90d5d69e284e28acc3b 100644 (file)
@@ -33,9 +33,6 @@
 #  include <boost/signals2.hpp>
 
 #  include <kinsol/kinsol.h>
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
-#    include <kinsol/kinsol_impl.h>
-#  endif
 #  include <nvector/nvector_serial.h>
 #  include <sundials/sundials_math.h>
 #  include <sundials/sundials_types.h>
index a030c20d73ebeb83d7c3b48d76f66516c8595229..436961963f948e5307c7a7d507df1f8943fbcf45 100644 (file)
 #  include <deal.II/lac/trilinos_vector.h>
 #  include <deal.II/lac/vector_memory.h>
 
-#  if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-#    include <deal.II/sundials/sundials_backport.h>
-#  endif
-
 DEAL_II_NAMESPACE_OPEN
 
 namespace SUNDIALS
@@ -1000,11 +996,9 @@ namespace SUNDIALS
       v->ops->nvcloneempty  = &NVectorOperations::clone_empty;
       v->ops->nvdestroy     = &NVectorOperations::destroy<VectorType>;
       //  v->ops->nvspace           = undef;
-#  if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
       v->ops->nvgetcommunicator =
         &NVectorOperations::get_communicator_as_void_ptr<VectorType>;
       v->ops->nvgetlength = &NVectorOperations::get_global_length<VectorType>;
-#  endif
 
       /* standard vector operations */
       v->ops->nvlinearsum = &NVectorOperations::linear_sum<VectorType>;
diff --git a/include/deal.II/sundials/sundials_backport.h b/include/deal.II/sundials/sundials_backport.h
deleted file mode 100644 (file)
index 39d528d..0000000
+++ /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 <deal.II/base/config.h>
-
-#ifdef DEAL_II_WITH_SUNDIALS
-
-#  include <sundials/sundials_nvector.h>
-
-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 (file)
index 66b8bc7..0000000
+++ /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 <deal.II/base/config.h>
-
-#ifdef DEAL_II_WITH_SUNDIALS
-
-#  include <sundials/sundials_linearsolver.h>
-
-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
index d5f90e3dba8d3bf8d7ea955187770a3996b06312..df9f8aa7d68b508869eb9861583a00f882f31699 100644 (file)
 #include <deal.II/base/config.h>
 
 #ifdef DEAL_II_WITH_SUNDIALS
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
 
-#    include <sundials/sundials_linearsolver.h>
+#  include <sundials/sundials_linearsolver.h>
 
-#    include <functional>
-#    include <memory>
+#  include <functional>
+#  include <memory>
 
 DEAL_II_NAMESPACE_OPEN
 
 namespace SUNDIALS
 {
-#    ifndef DOXYGEN
+#  ifndef DOXYGEN
   // forward declarations
   namespace internal
   {
     template <typename VectorType>
     struct LinearSolverContent;
   }
-#    endif
+#  endif
 
   /**
    * A linear operator that wraps SUNDIALS functionality.
@@ -167,10 +166,10 @@ namespace SUNDIALS
     {
     public:
       explicit LinearSolverWrapper(LinearSolveFunction<VectorType> 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
index fda97a4a05698e4ee4f84a994bfd00153f59ed76..6c48c6adec08a0d148ddc7883e52b8738f7cdb21 100644 (file)
 #  include <deal.II/sundials/n_vector.h>
 #  include <deal.II/sundials/sunlinsol_wrapper.h>
 
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-#    include <arkode/arkode_impl.h>
-#    include <sundials/sundials_config.h>
-#  else
-#    include <arkode/arkode_arkstep.h>
-#    include <sunlinsol/sunlinsol_spgmr.h>
-#    include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
-#    if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-#      include <deal.II/sundials/sunlinsol_newempty.h>
-#    endif
-#  endif
+#  include <arkode/arkode_arkstep.h>
+#  include <sunlinsol/sunlinsol_spgmr.h>
+#  include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
 
 #  include <iostream>
 
-#  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 <typename VectorType>
-    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<VectorType> &solver =
-        *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-
-      auto *src_ypred = internal::unwrap_nvector_const<VectorType>(ypred);
-      auto *src_fpred = internal::unwrap_nvector_const<VectorType>(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 <typename VectorType>
-    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<VectorType> &solver =
-        *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-
-      auto *dst      = internal::unwrap_nvector<VectorType>(b);
-      auto *src_ycur = internal::unwrap_nvector_const<VectorType>(ycur);
-      auto *src_fcur = internal::unwrap_nvector_const<VectorType>(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 <typename VectorType>
-    int
-    setup_mass_matrix_callback(ARKodeMem arkode_mem,
-                               N_Vector,
-                               N_Vector,
-                               N_Vector)
-    {
-      Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
-      ARKode<VectorType> &solver =
-        *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-      int err = solver.setup_mass(arkode_mem->ark_tn);
-      return err;
-    }
-
-
-
-    template <typename VectorType>
-    int
-    solve_with_mass_matrix_callback(ARKodeMem arkode_mem, N_Vector b)
-    {
-      Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
-      ARKode<VectorType> &solver =
-        *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
-
-      auto *dst = internal::unwrap_nvector<VectorType>(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 <typename VectorType>
     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 <typename VectorType>
   void
   ARKode<VectorType>::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<VectorType> : nullptr,
-      implicit_function ? &implicit_function_callback<VectorType> : 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<ARKodeMem>(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<VectorType>;
-        if (setup_jacobian)
-          {
-            ARKode_mem->ark_lsetup = setup_jacobian_callback<VectorType>;
-          }
-      }
-    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<VectorType>;
-
-        if (setup_mass)
-          {
-            ARKode_mem->ark_msetup = setup_mass_matrix_callback<VectorType>;
-          }
-      }
-
-    status = ARKodeSetOrder(arkode_mem, data.maximum_order);
-    AssertARKode(status);
-
-    if (custom_setup)
-      custom_setup(arkode_mem);
-  }
-
-#  else
-
-  template <typename VectorType>
-  void
-  ARKode<VectorType>::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<VectorType> : 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<internal::LinearSolverWrapper<VectorType>>(
                 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<internal::LinearSolverWrapper<VectorType>>(
                 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
 
 
 
index 801e1c0380aca90d411f0ca4e149fc20034ace4b..c07aa31977d36889f1816c04adc1bf6bcdd06ea2 100644 (file)
 
 #  include <deal.II/sundials/n_vector.h>
 
-#  include <sundials/sundials_config.h>
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-#    ifdef DEAL_II_SUNDIALS_WITH_IDAS
-#      include <idas/idas_impl.h>
-#    else
-#      include <ida/ida_impl.h>
-#    endif
-#  endif
-#  if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-#    include <deal.II/sundials/sunlinsol_newempty.h>
-#  endif
 #  include <iomanip>
 #  include <iostream>
 
@@ -77,68 +66,6 @@ namespace SUNDIALS
 
 
 
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-    template <typename VectorType>
-    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<VectorType> &solver =
-        *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
-
-      auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
-      auto *src_yp = internal::unwrap_nvector_const<VectorType>(yp);
-
-      int err = solver.setup_jacobian(IDA_mem->ida_tn,
-                                      *src_yy,
-                                      *src_yp,
-                                      IDA_mem->ida_cj);
-
-      return err;
-    }
-
-
-
-    template <typename VectorType>
-    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<VectorType> &solver =
-        *static_cast<IDA<VectorType> *>(IDA_mem->ida_user_data);
-      GrowingVectorMemory<VectorType> mem;
-
-      typename VectorMemory<VectorType>::Pointer dst(mem);
-      solver.reinit_vector(*dst);
-
-      auto *src = internal::unwrap_nvector<VectorType>(b);
-
-      int err = solver.solve_jacobian_system(*src, *dst);
-      *src    = *dst;
-
-      return err;
-    }
-
-
-
-#  else
     template <typename VectorType>
     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<IDAMem>(ida_mem);
-
-    IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
-
-    if (solve_jacobian_system)
-      IDA_mem->ida_lsolve = solve_with_jacobian_callback<VectorType>;
-    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<VectorType>);
     AssertIDA(status);
-#  endif
     status = IDASetMaxOrd(ida_mem, data.maximum_order);
     AssertIDA(status);
 
index 97869c426b68ce2b985a4ee929566ac4f88b87f4..66b40dfc846dda4c1f111a0e11224a4b9397dac4 100644 (file)
@@ -170,55 +170,6 @@ namespace SUNDIALS
 
 
 
-#  if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
-    template <typename VectorType>
-    int
-    setup_jacobian_callback(KINMem kinsol_mem)
-    {
-      KINSOL<VectorType> &solver =
-        *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
-
-      auto *src_ycur =
-        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_uu);
-      auto *src_fcur =
-        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_fval);
-
-      int err = solver.setup_jacobian(*src_ycur, *src_fcur);
-      return err;
-    }
-
-
-
-    template <typename VectorType>
-    int
-    solve_with_jacobian_callback(KINMem    kinsol_mem,
-                                 N_Vector  x,
-                                 N_Vector  b,
-                                 realtype *sJpnorm,
-                                 realtype *sFdotJp)
-    {
-      KINSOL<VectorType> &solver =
-        *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
-
-      auto *src_ycur =
-        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_uu);
-      auto *src_fcur =
-        internal::unwrap_nvector_const<VectorType>(kinsol_mem->kin_fval);
-      auto *src = internal::unwrap_nvector_const<VectorType>(b);
-      auto *dst = internal::unwrap_nvector<VectorType>(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 <typename VectorType>
     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<KINMem>(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<VectorType>;
-        if (setup_jacobian) // user assigned a function object to the Jacobian
-          // set-up slot
-          KIN_mem->kin_lsetup = setup_jacobian_callback<VectorType>;
-
-/* 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<VectorType>);
         AssertKINSOL(status);
-#  endif
       }
 
     // call to KINSol
index 9062d2f828c844f894e17994e60a94116eba2813..49f70c810e9a7b6015137562ae95be17ad942e92 100644 (file)
 #include <deal.II/sundials/sunlinsol_wrapper.h>
 
 #ifdef DEAL_II_WITH_SUNDIALS
-#  if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
 
-#    include <deal.II/base/exceptions.h>
+#  include <deal.II/base/exceptions.h>
 
-#    include <deal.II/lac/block_vector.h>
-#    include <deal.II/lac/la_parallel_block_vector.h>
-#    include <deal.II/lac/la_parallel_vector.h>
-#    include <deal.II/lac/vector.h>
-#    ifdef DEAL_II_WITH_TRILINOS
-#      include <deal.II/lac/trilinos_parallel_block_vector.h>
-#      include <deal.II/lac/trilinos_vector.h>
-#    endif
-#    ifdef DEAL_II_WITH_PETSC
-#      include <deal.II/lac/petsc_block_vector.h>
-#      include <deal.II/lac/petsc_vector.h>
-#    endif
-
-#    include <deal.II/sundials/n_vector.h>
+#  include <deal.II/lac/block_vector.h>
+#  include <deal.II/lac/la_parallel_block_vector.h>
+#  include <deal.II/lac/la_parallel_vector.h>
+#  include <deal.II/lac/vector.h>
+#  ifdef DEAL_II_WITH_TRILINOS
+#    include <deal.II/lac/trilinos_parallel_block_vector.h>
+#    include <deal.II/lac/trilinos_vector.h>
+#  endif
+#  ifdef DEAL_II_WITH_PETSC
+#    include <deal.II/lac/petsc_block_vector.h>
+#    include <deal.II/lac/petsc_vector.h>
+#  endif
 
-#    if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
-#      include <deal.II/sundials/sunlinsol_newempty.h>
-#    endif
+#  include <deal.II/sundials/n_vector.h>
 
 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 <typename VectorType>
   internal::LinearSolverWrapper<VectorType>::LinearSolverWrapper(
     LinearSolveFunction<VectorType> 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<LinearSolverContent<VectorType>>())
   {
-#    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<VectorType>;
@@ -240,13 +235,13 @@ namespace SUNDIALS
   SundialsOperator<VectorType>::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<Vector<double>>;
@@ -323,8 +318,8 @@ namespace SUNDIALS
   template class internal::LinearSolverWrapper<
     LinearAlgebra::distributed::BlockVector<double>>;
 
-#    ifdef DEAL_II_WITH_MPI
-#      ifdef DEAL_II_WITH_TRILINOS
+#  ifdef DEAL_II_WITH_MPI
+#    ifdef DEAL_II_WITH_TRILINOS
   template struct SundialsOperator<TrilinosWrappers::MPI::Vector>;
   template struct SundialsOperator<TrilinosWrappers::MPI::BlockVector>;
 
@@ -334,10 +329,10 @@ namespace SUNDIALS
   template class internal::LinearSolverWrapper<TrilinosWrappers::MPI::Vector>;
   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<PETScWrappers::MPI::Vector>;
   template struct SundialsOperator<PETScWrappers::MPI::BlockVector>;
@@ -347,12 +342,11 @@ namespace SUNDIALS
 
   template class internal::LinearSolverWrapper<PETScWrappers::MPI::Vector>;
   template class internal::LinearSolverWrapper<PETScWrappers::MPI::BlockVector>;
-#        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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.