From 7ec342318c15fb89048bb54b01671c5400f421ca Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Tue, 27 Oct 2020 17:09:17 +0100 Subject: [PATCH] SUNDIALS: update ARKode to version 5.4.0 Provide new std::functions to interface with SUNDIALS arkode module IDA and KINSOL not updated --- cmake/configure/configure_sundials.cmake | 20 - doc/news/changes/major/20201211Proell | 3 + include/deal.II/sundials/arkode.h | 731 ++++++++++++++-- include/deal.II/sundials/ida.h | 74 +- include/deal.II/sundials/kinsol.h | 35 +- include/deal.II/sundials/sunlinsol_newempty.h | 104 +++ source/sundials/arkode.cc | 806 +++++++++++++++++- source/sundials/ida.cc | 80 +- source/sundials/kinsol.cc | 90 +- ...harmonic_oscillator_02.cc => arkode_01.cc} | 8 +- ...rmonic_oscillator_02.prm => arkode_01.prm} | 0 .../arkode_01.with_sundials.geq.5.4.0.output | 128 +++ .../arkode_01.with_sundials.le.5.4.0.output | 128 +++ ...harmonic_oscillator_03.cc => arkode_02.cc} | 10 +- .../arkode_02.with_sundials.geq.5.4.0.output | 14 + .../arkode_02.with_sundials.le.4.0.0.output | 14 + .../arkode_02.with_sundials=4.0.0.output | 14 + ...harmonic_oscillator_04.cc => arkode_03.cc} | 10 +- ...rmonic_oscillator_04.prm => arkode_03.prm} | 0 .../arkode_03.with_sundials.geq.5.4.0.output | 12 + .../arkode_03.with_sundials.le.5.4.0.output | 12 + ...harmonic_oscillator_05.cc => arkode_04.cc} | 9 +- ...rmonic_oscillator_05.prm => arkode_04.prm} | 0 .../arkode_04.with_sundials.le.4.0.0.output | 102 +++ ...harmonic_oscillator_06.cc => arkode_05.cc} | 9 +- ...rmonic_oscillator_06.prm => arkode_05.prm} | 0 .../arkode_05.with_sundials.le.4.0.0.output | 102 +++ tests/sundials/arkode_06.cc | 159 ++++ tests/sundials/arkode_06.prm | 15 + .../arkode_06.with_sundials.geq.5.4.0.output | 102 +++ .../arkode_06.with_sundials=4.0.0.output | 102 +++ tests/sundials/arkode_07.cc | 183 ++++ tests/sundials/arkode_07.prm | 15 + .../arkode_07.with_sundials.geq.5.4.0.output | 102 +++ .../arkode_07.with_sundials=4.0.0.output | 102 +++ tests/sundials/arkode_08.cc | 162 ++++ tests/sundials/arkode_08.prm | 15 + .../arkode_08.with_sundials.geq.5.4.0.output | 279 ++++++ .../arkode_08.with_sundials=4.0.0.output | 89 ++ tests/sundials/harmonic_oscillator_02.output | 127 --- tests/sundials/harmonic_oscillator_03.output | 13 - tests/sundials/harmonic_oscillator_04.output | 11 - tests/sundials/harmonic_oscillator_05.output | 101 --- tests/sundials/harmonic_oscillator_06.output | 101 --- .../{harmonic_oscillator_01.cc => ida_01.cc} | 4 +- ...{harmonic_oscillator_01.prm => ida_01.prm} | 0 ...t => ida_01.with_sundials.le.4.0.0.output} | 0 ...> kinsol_01.with_sundials.le.4.0.0.output} | 0 tests/sundials/serial_sundials.cc | 3 +- ...al_sundials.with_sundials.geq.4.0.0.output | 2 + ...al_sundials.with_sundials.le.4.0.0.output} | 0 51 files changed, 3563 insertions(+), 639 deletions(-) create mode 100644 doc/news/changes/major/20201211Proell create mode 100644 include/deal.II/sundials/sunlinsol_newempty.h rename tests/sundials/{harmonic_oscillator_02.cc => arkode_01.cc} (90%) rename tests/sundials/{harmonic_oscillator_02.prm => arkode_01.prm} (100%) create mode 100644 tests/sundials/arkode_01.with_sundials.geq.5.4.0.output create mode 100644 tests/sundials/arkode_01.with_sundials.le.5.4.0.output rename tests/sundials/{harmonic_oscillator_03.cc => arkode_02.cc} (90%) create mode 100644 tests/sundials/arkode_02.with_sundials.geq.5.4.0.output create mode 100644 tests/sundials/arkode_02.with_sundials.le.4.0.0.output create mode 100644 tests/sundials/arkode_02.with_sundials=4.0.0.output rename tests/sundials/{harmonic_oscillator_04.cc => arkode_03.cc} (91%) rename tests/sundials/{harmonic_oscillator_04.prm => arkode_03.prm} (100%) create mode 100644 tests/sundials/arkode_03.with_sundials.geq.5.4.0.output create mode 100644 tests/sundials/arkode_03.with_sundials.le.5.4.0.output rename tests/sundials/{harmonic_oscillator_05.cc => arkode_04.cc} (93%) rename tests/sundials/{harmonic_oscillator_05.prm => arkode_04.prm} (100%) create mode 100644 tests/sundials/arkode_04.with_sundials.le.4.0.0.output rename tests/sundials/{harmonic_oscillator_06.cc => arkode_05.cc} (93%) rename tests/sundials/{harmonic_oscillator_06.prm => arkode_05.prm} (100%) create mode 100644 tests/sundials/arkode_05.with_sundials.le.4.0.0.output create mode 100644 tests/sundials/arkode_06.cc create mode 100644 tests/sundials/arkode_06.prm create mode 100644 tests/sundials/arkode_06.with_sundials.geq.5.4.0.output create mode 100644 tests/sundials/arkode_06.with_sundials=4.0.0.output create mode 100644 tests/sundials/arkode_07.cc create mode 100644 tests/sundials/arkode_07.prm create mode 100644 tests/sundials/arkode_07.with_sundials.geq.5.4.0.output create mode 100644 tests/sundials/arkode_07.with_sundials=4.0.0.output create mode 100644 tests/sundials/arkode_08.cc create mode 100644 tests/sundials/arkode_08.prm create mode 100644 tests/sundials/arkode_08.with_sundials.geq.5.4.0.output create mode 100644 tests/sundials/arkode_08.with_sundials=4.0.0.output delete mode 100644 tests/sundials/harmonic_oscillator_02.output delete mode 100644 tests/sundials/harmonic_oscillator_03.output delete mode 100644 tests/sundials/harmonic_oscillator_04.output delete mode 100644 tests/sundials/harmonic_oscillator_05.output delete mode 100644 tests/sundials/harmonic_oscillator_06.output rename tests/sundials/{harmonic_oscillator_01.cc => ida_01.cc} (96%) rename tests/sundials/{harmonic_oscillator_01.prm => ida_01.prm} (100%) rename tests/sundials/{harmonic_oscillator_01.output => ida_01.with_sundials.le.4.0.0.output} (100%) rename tests/sundials/{kinsol_01.output => kinsol_01.with_sundials.le.4.0.0.output} (100%) create mode 100644 tests/sundials/serial_sundials.with_sundials.geq.4.0.0.output rename tests/sundials/{serial_sundials.output => serial_sundials.with_sundials.le.4.0.0.output} (100%) diff --git a/cmake/configure/configure_sundials.cmake b/cmake/configure/configure_sundials.cmake index 87893e3806..62522094ce 100644 --- a/cmake/configure/configure_sundials.cmake +++ b/cmake/configure/configure_sundials.cmake @@ -22,26 +22,6 @@ MACRO(FEATURE_SUNDIALS_FIND_EXTERNAL var) IF(SUNDIALS_FOUND) SET(${var} TRUE) - - # - # We don't support version 4.0.0 or later yet. - # - SET(_first_unsupported_sundials_version 4.0.0) - IF(NOT SUNDIALS_VERSION VERSION_LESS ${_first_unsupported_sundials_version}) - MESSAGE(STATUS - "Insufficient SUNDIALS installation found: " - "version ${_first_unsupported_sundials_version} " - "or later is not yet supported, " - "but version ${SUNDIALS_VERSION} was found." - ) - SET(SUNDIALS_ADDITIONAL_ERROR_STRING - "Insufficient SUNDIALS installation found!\n" - "Version ${_first_unsupported_sundials_version} " - "or later is not yet supported, " - "but version ${SUNDIALS_VERSION} was found.\n" - ) - SET(${var} FALSE) - ENDIF() ENDIF() ENDMACRO() diff --git a/doc/news/changes/major/20201211Proell b/doc/news/changes/major/20201211Proell new file mode 100644 index 0000000000..aac3d03065 --- /dev/null +++ b/doc/news/changes/major/20201211Proell @@ -0,0 +1,3 @@ +Improved: Update SUNDIALS ARKODE interface to support versions > 4.0.0 and add preconditioner support +
+(Sebastian Proell, 2020/12/11) diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index 86435ad991..5f565e3640 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -34,13 +34,16 @@ # include # include -# include +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) +# include +# endif # include # ifdef DEAL_II_WITH_MPI # include # endif # include +# include # include # include @@ -49,6 +52,24 @@ DEAL_II_NAMESPACE_OPEN +// Forward declarations +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) +# ifndef DOXYGEN +namespace SUNDIALS +{ + // Forward declaration + template + struct SundialsOperator; + + template + struct SundialsPreconditioner; + + template + class SundialsLinearSolverWrapper; +} // namespace SUNDIALS +# endif +# endif + // Shorthand notation for ARKODE error codes. # define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code)) @@ -57,6 +78,41 @@ DEAL_II_NAMESPACE_OPEN */ namespace SUNDIALS { +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) + /** + * Type of function objects to interface with SUNDIALS linear solvers + * + * This function type encapsulates the action of solving $P^{-1}Ax=P^{-1}b$. + * The LinearOperator @p op encapsulates the matrix vector product $Ax$ and + * the LinearOperator @p prec encapsulates the application of the + * preconditioner $P^{-1}z$. + * The user can specify function objects of this type to attach custom linear + * solver routines to SUNDIALS. The two LinearOperators @p op and @p prec are + * built internally by SUNDIALS based on user settings. The parameters are + * interpreted as follows: + * + * @param[in] op A LinearOperator that applies the matrix vector product + * @param[in] prec A LinearOperator that applies the preconditioner + * @param[out] x The output solution vector + * @param[in] b The right-hand side + * @param[in] tol Tolerance for the iterative solver + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + template + using LinearSolveFunction = + std::function & op, + SundialsPreconditioner &prec, + VectorType & x, + const VectorType & b, + double tol)>; +# endif + /** * Interface to SUNDIALS additive Runge-Kutta methods (ARKode). * @@ -192,9 +248,10 @@ namespace SUNDIALS * $f_I(t, y)$ depends linearly on $y$, and if the Newton-based nonlinear * solver is chosen, then the system will be solved using only a single * Newton iteration. Notice that in order for the Newton solver to be used, - * at least the solve_jacobian_system() function should be supplied. If this - * function is not supplied, then only the fixed-point iteration will be - * supported, and the `implicit_function_is_linear` setting is ignored. + * at least the jacobian_times_vector() function (or solve_jacobian_system() + * for SUNDIALS version > 4.0.0) should be supplied. If this function is not + * supplied, then only the fixed-point iteration will be supported, and the + *`implicit_function_is_linear` setting is ignored. * * The optimal solver (Newton vs fixed-point) is highly problem-dependent. * Since fixed-point solvers do not require the solution of any linear @@ -221,30 +278,49 @@ namespace SUNDIALS * $z_j$ where $j < i$). Additional information on the specific predictor * algorithms implemented in ARKode is provided in ARKode documentation. * - * The user has to provide the implementation of the following std::function: - * - reinit_vector; + * The user has to provide the implementation of the following + *`std::function`s: + * - reinit_vector() * and either one or both of - * - implicit_function; - * - explicit_function; + * - implicit_function() + * - explicit_function() * * If the mass matrix is different from the identity, the user should supply - * - solve_mass_system; - * and, optionally, - * - setup_mass; + * - mass_times_vector() (or solve_mass_system() for SUNDIALS version + * < 4.0.0) and, optionally, + * - mass_times_setup() (or setup_mass() for SUNDIALS version < 4.0.0) * * If the use of a Newton method is desired, then the user should also supply - * - solve_jacobian_system; - * and optionally - * - setup_jacobian; + * - jacobian_times_vector (or solve_jacobian_system() for SUNDIALS version + * < 4.0.0) + * - optional: jacobian_times_setup() (or setup_jacobian() for SUNDIALS + * version < 4.0.0) + * + * @note Although SUNDIALS can provide a difference quotient approximation + * of the Jacobian, this is currently not supported through this wrapper. + * + * Only for SUNDIALS version > 4.0.0: A SUNDIALS default solver (SPGMR) is + * used to solve the linear systems. To use a custom linear solver for the + * mass matrix and/or Jacobian, set: + * - solve_mass() and/or + * - solve_jacobian() + * + * Only for SUNDIALS version > 4.0.0: To use a custom preconditioner with + * either a default or custom linear solver, set: + * - jacobian_preconditioner_solve() and/or mass_preconditioner_solve() + * and, optionally, + * - jacobian_preconditioner_setup() and/or mass_preconditioner_setup() * * Also the following functions could be rewritten. By default * they do nothing, or are not required. - * - solver_should_restart; - * - get_local_tolerances; + * - solver_should_restart() + * - get_local_tolerances() * - * To produce output at fixed steps, overload the function - * - output_step; + * To produce output at fixed steps, set the function + * - output_step() * + * Any other custom settings of the ARKODE object can be specified in + * - custom_setup() * * To provide a simple example, consider the harmonic oscillator problem: * \f[ @@ -332,11 +408,15 @@ namespace SUNDIALS * @param minimum_step_size Minimum step size * @param maximum_order Maximum ARK order * @param maximum_non_linear_iterations Maximum number of nonlinear - * iterations + * iterations * @param implicit_function_is_linear Specifies that the implicit portion - * of the problem is linear + * of the problem is linear * @param implicit_function_is_time_independent Specifies that the - * implicit portion of the problem is linear and time independent + * implicit portion of the problem is linear and time independent + * @param mass_is_time_independent Specifies that the mass pre-factor is + * independent of time + * @param anderson_acceleration_subspace The number of vectors to use for + * Anderson acceleration within the packaged SUNDIALS solver. * * Error parameters: * @@ -355,6 +435,8 @@ namespace SUNDIALS const unsigned int maximum_non_linear_iterations = 10, const bool implicit_function_is_linear = false, const bool implicit_function_is_time_independent = false, + const bool mass_is_time_independent = false, + const int anderson_acceleration_subspace = 3, // Error parameters const double absolute_tolerance = 1e-6, const double relative_tolerance = 1e-5) @@ -370,6 +452,8 @@ namespace SUNDIALS , implicit_function_is_linear(implicit_function_is_linear) , implicit_function_is_time_independent( implicit_function_is_time_independent) + , mass_is_time_independent(mass_is_time_independent) + , anderson_acceleration_subspace(anderson_acceleration_subspace) {} /** @@ -377,34 +461,6 @@ namespace SUNDIALS * object. When the parameters are parsed from a file, the internal * parameters are automatically updated. * - * The following parameters are declared: - * - * @code - * set Final time = 1.000000 - * set Initial time = 0.000000 - * set Time interval between each output = 0.2 - * subsection Error control - * set Absolute error tolerance = 0.000001 - * set Ignore algebraic terms for error computations = true - * set Relative error tolerance = 0.00001 - * set Use local tolerances = false - * end - * subsection Initial condition correction parameters - * set Correction type at initial time = none - * set Correction type after restart = none - * set Maximum number of nonlinear iterations = 5 - * end - * subsection Running parameters - * set Initial step size = 0.1 - * set Maximum number of nonlinear iterations = 10 - * set Maximum order of ARK = 5 - * set Minimum step size = 0.000001 - * end - * @endcode - * - * These are one-to-one with the options you can pass at construction - * time. - * * The options you pass at construction time are set as default values in * the ParameterHandler object `prm`. You can later modify them by parsing * a parameter file using `prm`. The values of the parameter will be @@ -431,6 +487,9 @@ namespace SUNDIALS implicit_function_is_linear); prm.add_parameter("Implicit function is time independent", implicit_function_is_time_independent); + prm.add_parameter("Mass is time independent", mass_is_time_independent); + prm.add_parameter("Anderson-acceleration subspace", + anderson_acceleration_subspace); prm.leave_subsection(); prm.enter_subsection("Error control"); @@ -496,6 +555,18 @@ namespace SUNDIALS * independent. */ bool implicit_function_is_time_independent; + + /** + * Specify whether the mass pre-factor is time independent. Has no effect + * if no mass is specified. + */ + bool mass_is_time_independent; + + /** + * Number of subspace vectors to use for Anderson acceleration. Only + * meaningful if the packaged SUNDIALS fixed-point solver is used. + */ + int anderson_acceleration_subspace; }; /** @@ -541,6 +612,44 @@ namespace SUNDIALS void reset(const double t, const double h, const VectorType &y); + /*! + * Create a new SUNDIALS vector from a given template. + * + * @note Vectors created this way should be freed with free_vector(). + * + * @param template_vector The vector to use as a template for the layout of + * a new vector. + */ + N_Vector + create_vector(const VectorType &template_vector) const; + + /** + * Free a SUNDIALS vector created with create_vector(). + * + * @param vector The vector to free + */ + void + free_vector(N_Vector vector) const; + + /** + * Provides user access to the internally used ARKODE memory. + * + * This functionality is intended for users who wish to query additional + * information directly from the ARKODE integrator, refer to the ARKODE + * manual for the various `ARKStepGet...` functions. The `ARKStepSet...` + * functions should not be called since this might lead to conflicts with + * various settings that are performed by this ARKode object. + * + * @note If custom settings of ARKODE functionality (that are not achievable + * via the interface of this class) are required, the function + * custom_setup() should be used. + * + * @return pointer to the ARKODE memory block that can be passed to SUNDIALS + * functions + */ + void * + get_arkode_memory() const; + /** * A function object that users need to supply and that is intended to * reinit the given vector. @@ -558,10 +667,10 @@ namespace SUNDIALS * * 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. + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. */ std::function< int(const double t, const VectorType &y, VectorType &explicit_f)> @@ -578,14 +687,15 @@ namespace SUNDIALS * * 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. + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. */ 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 @@ -668,7 +778,7 @@ namespace SUNDIALS * - >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. + * assertion will be thrown. */ std::function0: 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. + * assertion will be thrown. */ std::function0: 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. + * assertion will be thrown. */ std::function setup_mass; @@ -782,10 +892,325 @@ namespace SUNDIALS * - >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. + * 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 + * will be called by ARKode (possibly several times) after + * mass_times_setup() has been called at least once. ARKode tries to do its + * best to call mass_times_setup() the minimum amount of times. + * + * A call to this function should store in `Mv` the result of $M$ + * applied to `v`. + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function + mass_times_vector; + + /** + * 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 guaranteed to be called by ARKode at least once, before + * any call to mass_times_vector(). + * + * 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 mass_times_vector() function, then the + * identity is used. If the mass_times_setup() function is not provided, + * then mass_times_vector() 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 should be called to + * assemble the matrix. Subsequent calls (possibly more than one) to + * mass_times_vector() can assume that this function has been called at + * least once. + * + * @note No assumption is made by this interface on what the user + * should do in this function. ARKode only assumes that after a call to + * mass_times_setup() it is possible to call mass_times_vector(). + * + * @param t The current evaluation time + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function mass_times_setup; + + /** + * A function object that users may supply and that is intended to compute + * the product of the Jacobian matrix with a given vector `v`. The Jacobian + * here refers to $J=\frac{\partial f_I}{\partial y}$, i.e., the Jacobian of + * the user-specified implicit_function. + * + * A call to this function should store in `Jv` the result of $J$ + * applied to `v`. + * + * Arguments to the function are + * + * @param[in] v The vector to be mulitplied by the Jacobian + * @param[out] Jv The vector to be filled with the product J*v + * @param[in] t The current time + * @param[in] y The current $y$ vector for the current ARKode internal + * step + * @param[in] fy The current value of the implicit right-hand side at y, + * $f_I (t_n, y)$. + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function + jacobian_times_vector; + + /** + * A function object that users may supply and that is intended to set up + * all data necessary for the application of jacobian_times_setup(). This + * function is called by ARKode any time a Jacobian update is required. + * The user should compute the Jacobian (or update all the variables that + * allow the application of Jacobian). This function is guaranteed to + * be called by ARKode at least once, before any call to + * jacobian_times_vector(). + * + * If the jacobian_times_setup() function is not provided, then + * jacobian_times_vector() should do all the work by itself. + * + * If the user uses a matrix based computation of the Jacobian, then this is + * the right place where an assembly routine should be called to assemble + * the matrix. Subsequent calls (possibly more than one) to + * jacobian_times_vector() can assume that this function has been called at + * least once. + * + * @note No assumption is made by this interface on what the user + * should do in this function. ARKode only assumes that after a call to + * jacobian_times_setup() it is possible to call jacobian_times_vector(). + * + * @param t The current time + * @param y The current ARKode internal solution vector $y$ + * @param fy The implicit right-hand side function evaluated at the + * current time $t$ and state $y$, i.e., $f_I(y,t)$ + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function + jacobian_times_setup; + + /** + * A LinearSolveFunction object that users may supply and that is intended + * to solve the linearized system $Ax=b$, where $A = M-\gamma J$ is the + * Jacobian of the nonlinear residual. The application fo the mass matrix + * $M$ and Jacobian $J$ are known through the functions mass_times_vector() + * and jacobian_times_vector() and $\gamma$ is a factor provided by + * SUNDIALS. The matrix-vector product $Ax$ is encoded in the supplied + * SundialsOperator. If a preconditioner was set through + * jacobian_preconditioner_solve(), it is encoded in the + * SundialsPreconditioner. If no preconditioner was supplied this way, the + * preconditioner is the identity matrix, i.e., no preconditioner. The user + * is free to use a custom preconditioner in this function object that is + * not supplied through SUNDIALS. + * + * If you do not specify a solve_linearized_system() function, then a + * SUNDIALS packaged SPGMR solver with default settings is used. + * + * For more details on the function type refer to LinearSolveFunction. + */ + LinearSolveFunction solve_linearized_system; + + /** + * A LinearSolveFunction object that users may supply and that is intended + * to solve the mass system $Mx=b$. The matrix-vector product $Mx$ is + * encoded in the supplied SundialsOperator. If a preconditioner was set + * through mass_preconditioner_solve(), it is encoded in the + * SundialsPreconditioner. If no preconditioner was supplied this way, the + * preconditioner is the identity matrix, i.e., no preconditioner. The user + * is free to use a custom preconditioner in this function object that is + * not supplied through SUNDIALS. + * + * The user must specify this function if a non-identity mass matrix is used + * and applied in mass_times_vector(). + * + * For more details on the function type refer to LinearSolveFunction. + */ + LinearSolveFunction solve_mass; + + + /** + * A function object that users may supply to either pass a preconditioner + * to a SUNDIALS built-in solver or to apply a custom preconditioner within + * the user's own linear solve specified in solve_linearized_system(). + * + * This function should compute the solution to the preconditioner equation + * $Pz=r$ and store it in @p z. In this equation $P$ should approximate the + * Jacobian $M-\gamma J$ of the nonlinear system. + * + * @param[in] t The current time + * @param[in] y The current $y$ vector for the current ARKode internal + * step + * @param[in] fy The current value of the implicit right-hand side at y, + * $f_I (t_n, y)$. + * @param[in] r The right-hand side of the preconditioner equation + * @param[out] z The solution of applying the preconditioner, i.e., solving + * $Pz=r$ + * @param[in] gamma The value $\gamma$ in the preconditioner equation + * @param[in] tol The tolerance up to which the system should be solved + * @param[in] lr An input flag indicating whether the preconditioner solve + * is to use the left preconditioner (lr = 1) or the right preconditioner + * (lr = 2). Only relevant if used with a SUNDIALS packaged solver. If + * used with a custom solve_mass() function this will be set to zero. + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function + jacobian_preconditioner_solve; + + /** + * A function object that users may supply to set up a preconditioner + * specified in jacobian_preconditioner_solve(). + * + * This function should prepare the solution of the preconditioner equation + * $Pz=r$. In this equation $P$ should approximate the Jacobian $M-\gamma J$ + * of the nonlinear system. + * + * If the jacobian_preconditioner_setup() function is not provided, then + * jacobian_preconditioner_solve() should do all the work by itself. + * + * @note No assumption is made by this interface on what the user + * should do in this function. ARKode only assumes that after a call to + * jacobian_preconditioner_setup() it is possible to call + * jacobian_preconditioner_solve(). + * + * @param[in] t The current time + * @param[in] y The current $y$ vector for the current ARKode internal + * step + * @param[in] fy The current value of the implicit right-hand side at y, + * $f_I (t_n, y)$. + * @param[in] jok An input flag indicating whether the Jacobian-related + * data needs to be updated. The jok argument provides for the reuse of + * Jacobian data in the preconditioner solve function. When jok = + * SUNFALSE, the Jacobian-related data should be recomputed from scratch. + * When jok = SUNTRUE the Jacobian data, if saved from the previous call + * to this function, can be reused (with the current value of gamma). A + * call with jok = SUNTRUE can only occur after a call with jok = + * SUNFALSE. + * @param[out] jcur On output this should be set to SUNTRUE if Jacobian data + * was recomputed, or set to SUNFALSE if Jacobian data was not recomputed, + * but saved data was still reused. + * @param[in] gamma The value $\gamma$ in $M-\gamma J$. The preconditioner + * should approximate the inverse of this matrix. + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function + jacobian_preconditioner_setup; + + /** + * A function object that users may supply to either pass a preconditioner + * to a SUNDIALS built-in solver or to apply a custom preconditioner within + * the user's own linear solve specified in solve_mass(). + * + * This function should compute the solution to the preconditioner equation + * $Pz=r$ and store it in @p z. In this equation $P$ should approximate the + * mass matrix $M$. + * + * @param[in] t The current time + * @param[in] r The right-hand side of the preconditioner equation + * @param[out] z The solution of applying the preconditioner, i.e., solving + * $Pz=r$ + * @param[in] gamma The value $\gamma$ in the preconditioner equation + * @param[in] tol The tolerance up to which the system should be solved + * @param[in] lr An input flag indicating whether the preconditioner solve + * is to use the left preconditioner (lr = 1) or the right preconditioner + * (lr = 2). Only relevant if used with a SUNDIALS packaged solver. If + * used with a custom solve_mass() function this will be set to zero. + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function< + int(double t, const VectorType &r, VectorType &z, double tol, int lr)> + mass_preconditioner_solve; + + /** + * A function object that users may supply to set up a preconditioner + * specified in mass_preconditioner_solve(). + * + * This function should prepare the solution of the preconditioner equation + * $Pz=r$. In this equation $P$ should approximate the mass matrix $M$. + * + * If the mass_preconditioner_setup() function is not provided, then + * mass_preconditioner_solve() should do all the work by itself. + * + * @note No assumption is made by this interface on what the user + * should do in this function. ARKode only assumes that after a call to + * mass_preconditioner_setup() it is possible to call + * mass_preconditioner_solve(). + * + * @param[in] t The current time + * + * This function should return: + * - 0: Success + * - >0: Recoverable error, ARKode will reattempt the solution and call this + * function again. + * - <0: Unrecoverable error, the computation will be aborted and an + * assertion will be thrown. + */ + std::function mass_preconditioner_setup; +# endif /** * A function object that users may supply and that is intended to @@ -794,12 +1219,12 @@ namespace SUNDIALS * polynomial interpolation of the solution, computed using the current ARK * order and the (internally stored) previously computed solution steps. * - * Notice that it is well possible that internally ARKode computes a time - * step which is much larger than the `output_period` step, and therefore - * calls this function consecutively several times by simply performing all - * intermediate interpolations. There is no relationship between how many - * times this function is called and how many time steps have actually been - * computed. + * @note It is well possible that internally ARKode computes a time + * step which is much larger than the `output_period` step, and therefore + * calls this function consecutively several times by simply performing + * all intermediate interpolations. There is no relationship between how + * many times this function is called and how many time steps have + * actually been computed. */ std::function get_local_tolerances; /** - * Handle ARKode exceptions. + * A function object that users may supply and which is intended to perform + * custom settings on the supplied @p arkode_mem object. Refer to the + * SUNDIALS documentation for valid options. + * + * For instance, the following code attaches two files for diagnostic and + * error output of the internal ARKODE implementation: + * + * @code + * ode.custom_setup = [&](void *arkode_mem) { + * ARKStepSetErrFile(arkode_mem, errfile); + * ARKStepSetDiagnostics(arkode_mem, diagnostics_file); + * }; + * @endcode + * + * @note This function will be called at the end of all other set up right + * before the actual time evloution is started or continued with + * solve_ode(). This function is also called when the solver is restarted, + * see solver_should_restart(). Consult the SUNDIALS manual to see which + * options are still available at this point. + * + * @param arkode_mem pointer to the ARKODE memory block which can be used + * for custom calls to `ARKStepSet...` methods. */ - DeclException1(ExcARKodeError, - int, - << "One of the SUNDIALS ARKode internal functions " - << " returned a negative error code: " << arg1 - << ". Please consult SUNDIALS manual."); - + std::function custom_setup; private: /** @@ -853,6 +1294,26 @@ namespace SUNDIALS << "Please provide an implementation for the function \"" << arg1 << "\""); +# 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. + * @param solution The solution vector whihc is used as a template to create + * new vectors. + */ + void + setup_system_solver(const VectorType &solution); + + /** + * Set up the solver and preconditioner for a non-identity mass matrix in + * the ARKODE memory object based on the user-specified functions. + */ + void + setup_mass_solver(); + +# endif + /** * This function is executed at construction time to set the * std::function above to trigger an assert if they are not @@ -893,6 +1354,11 @@ namespace SUNDIALS */ GrowingVectorMemory mem; +# 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 static_assert(!std::is_same::value, @@ -907,6 +1373,117 @@ namespace SUNDIALS # endif // DEAL_II_WITH_PETSC }; +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) + + /** + * A linear operator that wraps SUNDIALS functionality. + */ + template + struct SundialsOperator + { + /** + * Apply this LinearOperator to @p src and store the result in @p dst. + */ + void + vmult(VectorType &dst, const VectorType &src) const; + + /** + * Constructor. + * + * @param solver The ARKode solver that uses this operator + * @param A_data Data required by @p a_times_fn + * @param a_times_fn A function pointer to the function that computes A*v + */ + SundialsOperator(ARKode &solver, + void * A_data, + ATimesFn a_times_fn); + + private: + /** + * Reference to the ARKode object that uses this SundialsOperator. + */ + ARKode &solver; + /** + * Data necessary to evaluate a_times_fn. + */ + void *A_data; + + /** + * Function pointer declared by SUNDIALS to evaluate the matrix vector + * product. + */ + ATimesFn a_times_fn; + }; + + + + /** + * A linear operator that wraps preconditioner functionality as specified by + * SUNDIALS. The vmult() function solves the preconditioner equation $Px=b$, + * i.e., it computes $x=P^{-1}b$. + */ + template + struct SundialsPreconditioner + { + /** + * Apply the wrapped preconditioner, i.e., solve $Px=b$ where $x$ is the + * @p dst vector and $b$ the @p src vector. + * + * @param dst Result vector of the preconditioner application + * @param src Target vector of the preconditioner application + */ + void + vmult(VectorType &dst, const VectorType &src) const; + + /** + * Constructor. + * + * @param solver The ARKode solver that uses this operator + * @param P_data Data required by @p p_solve_fn + * @param p_solve_fn A function pointer to the function that computes A*v + * @param tol Tolerance, that an iterative solver should use to judge + * convergence + */ + SundialsPreconditioner(ARKode &solver, + void * P_data, + PSolveFn p_solve_fn, + double tol); + + private: + /** + * Reference to the ARKode object that uses this SundialsPreconditioner. + */ + ARKode &solver; + + /** + * Data necessary to calls p_solve_fn + */ + void *P_data; + + /** + * Function pointer to a function that computes the preconditioner + * application. + */ + PSolveFn p_solve_fn; + + /** + * Potential tolerance to use in the internal solve of the preconditioner + * equation. + */ + double tol; + }; + +# endif + + /** + * Handle ARKode exceptions. + */ + DeclException1(ExcARKodeError, + int, + << "One of the SUNDIALS ARKode internal functions " + << " returned a negative error code: " << arg1 + << ". Please consult SUNDIALS manual."); + } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index d6cfca4df9..268bb55ca7 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -20,43 +20,44 @@ #include #ifdef DEAL_II_WITH_SUNDIALS - -# include -# include -# include -# include -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif -# include -# include - -# ifdef DEAL_II_SUNDIALS_WITH_IDAS -# include -# else -# include -# endif - -# include -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) -# include -# include -# include -# endif -# include - -# include -# include -# include - -# include +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) + +# include +# include +# include +# include +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif +# include +# include + +# ifdef DEAL_II_SUNDIALS_WITH_IDAS +# include +# else +# include +# endif + +# include +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# include +# include +# include +# endif +# include + +# include +# include +# include + +# include DEAL_II_NAMESPACE_OPEN // Shorthand notation for IDA error codes. -# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code)) +# define AssertIDA(code) Assert(code >= 0, ExcIDAError(code)) namespace SUNDIALS { @@ -833,8 +834,8 @@ namespace SUNDIALS */ GrowingVectorMemory mem; -# ifdef DEAL_II_WITH_PETSC -# ifdef PETSC_USE_COMPLEX +# ifdef DEAL_II_WITH_PETSC +# ifdef PETSC_USE_COMPLEX static_assert(!std::is_same::value, "Sundials does not support complex scalar types, " "but PETSc is configured to use a complex scalar type!"); @@ -843,14 +844,15 @@ namespace SUNDIALS !std::is_same::value, "Sundials does not support complex scalar types, " "but PETSc is configured to use a complex scalar type!"); -# endif // PETSC_USE_COMPLEX -# endif // DEAL_II_WITH_PETSC +# endif // PETSC_USE_COMPLEX +# endif // DEAL_II_WITH_PETSC }; } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE +# endif #endif // DEAL_II_WITH_SUNDIALS #endif diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h index 88793da27a..fbf3c0a584 100644 --- a/include/deal.II/sundials/kinsol.h +++ b/include/deal.II/sundials/kinsol.h @@ -17,32 +17,35 @@ #define dealii_sundials_kinsol_h #include + #ifdef DEAL_II_WITH_SUNDIALS +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) + -# include -# include -# include -# include -# include +# include +# include +# include +# include +# include -# include -# include +# include +# include -# include +# include -# include -# include -# include -# include -# include +# include +# include +# include +# include +# include -# include +# include DEAL_II_NAMESPACE_OPEN // Shorthand notation for KINSOL error codes. -# define AssertKINSOL(code) Assert(code >= 0, ExcKINSOLError(code)) +# define AssertKINSOL(code) Assert(code >= 0, ExcKINSOLError(code)) namespace SUNDIALS { @@ -677,7 +680,7 @@ namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE +# endif #endif - #endif diff --git a/include/deal.II/sundials/sunlinsol_newempty.h b/include/deal.II/sundials/sunlinsol_newempty.h new file mode 100644 index 0000000000..b992d66c23 --- /dev/null +++ b/include/deal.II/sundials/sunlinsol_newempty.h @@ -0,0 +1,104 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2020 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`. + */ + 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 + */ + 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/source/sundials/arkode.cc b/source/sundials/arkode.cc index b0dd8c672d..0a52dfeced 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -35,14 +35,24 @@ # include -# 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 +# 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 @@ -59,6 +69,7 @@ namespace SUNDIALS N_Vector yp, void * user_data) { + Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); GrowingVectorMemory mem; @@ -87,6 +98,7 @@ namespace SUNDIALS N_Vector yp, void * user_data) { + Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); GrowingVectorMemory mem; @@ -108,6 +120,7 @@ namespace SUNDIALS +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) template int t_arkode_setup_jacobian(ARKodeMem arkode_mem, @@ -119,6 +132,7 @@ namespace SUNDIALS N_Vector, N_Vector) { + Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(arkode_mem->ark_user_data); GrowingVectorMemory mem; @@ -140,11 +154,11 @@ namespace SUNDIALS *src_ypred, *src_fpred, jcurPtr_tmp); -# if DEAL_II_SUNDIALS_VERSION_GTE(2, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(2, 0, 0) *jcurPtr = jcurPtr_tmp ? SUNTRUE : SUNFALSE; -# else +# else *jcurPtr = jcurPtr_tmp ? TRUE : FALSE; -# endif +# endif return err; } @@ -155,12 +169,13 @@ namespace SUNDIALS int t_arkode_solve_jacobian(ARKodeMem arkode_mem, N_Vector b, -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) N_Vector, -# endif +# endif N_Vector ycur, N_Vector fcur) { + Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(arkode_mem->ark_user_data); GrowingVectorMemory mem; @@ -198,6 +213,7 @@ namespace SUNDIALS int t_arkode_setup_mass(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); @@ -209,14 +225,15 @@ namespace SUNDIALS template int t_arkode_solve_mass(ARKodeMem arkode_mem, -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) N_Vector b, N_Vector -# else +# else N_Vector b -# endif +# endif ) { + Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(arkode_mem->ark_user_data); GrowingVectorMemory mem; @@ -234,8 +251,380 @@ namespace SUNDIALS return err; } +# else + + template + int + t_arkode_jac_times_vec_function(N_Vector v, + N_Vector Jv, + realtype t, + N_Vector y, + N_Vector fy, + void * user_data, + N_Vector) + { + Assert(user_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_v(mem); + solver.reinit_vector(*src_v); + typename VectorMemory::Pointer dst_Jv(mem); + solver.reinit_vector(*dst_Jv); + typename VectorMemory::Pointer src_y(mem); + solver.reinit_vector(*src_y); + typename VectorMemory::Pointer src_fy(mem); + solver.reinit_vector(*src_fy); + copy(*src_v, v); + copy(*src_y, y); + copy(*src_fy, fy); + + int err = + solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy); + copy(Jv, *dst_Jv); + return err; + } + + + + template + int + t_arkode_jac_times_setup_function(realtype t, + N_Vector y, + N_Vector fy, + void * user_data) + { + Assert(user_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(user_data); + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_y(mem); + solver.reinit_vector(*src_y); + typename VectorMemory::Pointer src_fy(mem); + solver.reinit_vector(*src_fy); + copy(*src_y, y); + copy(*src_fy, fy); + int err = solver.jacobian_times_setup(t, *src_y, *src_fy); + return err; + } + + + + template + int + t_arkode_prec_solve_function(realtype t, + N_Vector y, + N_Vector fy, + N_Vector r, + N_Vector z, + realtype gamma, + realtype delta, + int lr, + void * user_data) + { + Assert(user_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(user_data); + + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_y(mem); + solver.reinit_vector(*src_y); + typename VectorMemory::Pointer src_fy(mem); + solver.reinit_vector(*src_fy); + typename VectorMemory::Pointer src_r(mem); + solver.reinit_vector(*src_r); + typename VectorMemory::Pointer dst_z(mem); + solver.reinit_vector(*dst_z); + copy(*src_y, y); + copy(*src_fy, fy); + copy(*src_r, r); + + int status = solver.jacobian_preconditioner_solve( + t, *src_y, *src_fy, *src_r, *dst_z, gamma, delta, lr); + copy(z, *dst_z); + return status; + } + + + + template + int + t_arkode_prec_setup_function(realtype t, + N_Vector y, + N_Vector fy, + booleantype jok, + booleantype *jcurPtr, + realtype gamma, + void * user_data) + { + Assert(user_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(user_data); + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_y(mem); + solver.reinit_vector(*src_y); + typename VectorMemory::Pointer src_fy(mem); + solver.reinit_vector(*src_fy); + copy(*src_y, y); + copy(*src_fy, fy); + return solver.jacobian_preconditioner_setup( + t, *src_y, *src_fy, jok, *jcurPtr, gamma); + } + + + + template + int + t_arkode_mass_times_vec_function(N_Vector v, + N_Vector Mv, + realtype t, + void * mtimes_data) + { + Assert(mtimes_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(mtimes_data); + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_v(mem); + solver.reinit_vector(*src_v); + typename VectorMemory::Pointer dst_Mv(mem); + solver.reinit_vector(*dst_Mv); + + copy(*src_v, v); + int err = solver.mass_times_vector(t, *src_v, *dst_Mv); + copy(Mv, *dst_Mv); + + return err; + } + + + + template + int + t_arkode_mass_times_setup_function(realtype t, void *mtimes_data) + { + Assert(mtimes_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(mtimes_data); + + return solver.mass_times_setup(t); + } + + + + template + int + t_arkode_mass_prec_solve_function(realtype t, + N_Vector r, + N_Vector z, + realtype delta, + int lr, + void * user_data) + { + Assert(user_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(user_data); + + GrowingVectorMemory mem; + typename VectorMemory::Pointer src_r(mem); + solver.reinit_vector(*src_r); + typename VectorMemory::Pointer dst_z(mem); + solver.reinit_vector(*dst_z); + copy(*src_r, r); + + int status = + solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr); + copy(z, *dst_z); + return status; + } + + + + template + int + t_arkode_mass_prec_setup_function(realtype t, void *user_data) + { + Assert(user_data != nullptr, ExcInternalError()); + ARKode &solver = + *static_cast *>(user_data); + + return solver.mass_preconditioner_setup(t); + } + + + + /** + * storage for internal content of the linear solver wrapper + */ + template + struct LinearSolverContent + { + ATimesFn a_times_fn; + PSetupFn preconditioner_setup; + PSolveFn preconditioner_solve; + + LinearSolveFunction lsolve; + + //! solver reference to forward all solver related calls to user-specified + //! functions + ARKode *solver; + void * P_data; + void * A_data; + }; + + + + /** + * Access our LinearSolverContent from the generic content of the + * SUNLinearSolver @p ls. + */ + template + LinearSolverContent * + access_content(SUNLinearSolver ls) + { + Assert(ls->content != nullptr, ExcInternalError()); + return static_cast *>(ls->content); + } + + + + SUNLinearSolver_Type arkode_linsol_get_type(SUNLinearSolver) + { + return SUNLINEARSOLVER_ITERATIVE; + } + + + + template + int + arkode_linsol_solve(SUNLinearSolver LS, + SUNMatrix, + N_Vector x, + N_Vector b, + realtype tol) + { + auto content = access_content(LS); + ARKode &solver = *(content->solver); + + GrowingVectorMemory mem; + typename VectorMemory::Pointer rhs(mem); + solver.reinit_vector(*rhs); + typename VectorMemory::Pointer dst(mem); + solver.reinit_vector(*dst); + + copy(*rhs, b); + + SundialsOperator op(solver, + content->A_data, + content->a_times_fn); + SundialsPreconditioner preconditioner( + solver, content->P_data, content->preconditioner_solve, tol); + + int err = content->lsolve(op, preconditioner, *dst, *rhs, tol); + copy(x, *dst); + return err; + } + + + + template + int + arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix) + { + auto content = access_content(LS); + if (content->preconditioner_setup) + return content->preconditioner_setup(content->P_data); + return 0; + } + + + + template + int arkode_linsol_initialize(SUNLinearSolver) + { + // this method is currently only provided because SUNDIALS 4.0.0 requires + // it - no user-set action is implemented so far + return 0; + } + + + + template + int + arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes) + { + auto content = access_content(LS); + content->A_data = A_data; + content->a_times_fn = ATimes; + return 0; + } + + + + template + int + arkode_linsol_set_preconditioner(SUNLinearSolver LS, + void * P_data, + PSetupFn p_setup, + PSolveFn p_solve) + { + auto content = access_content(LS); + content->P_data = P_data; + content->preconditioner_setup = p_setup; + content->preconditioner_solve = p_solve; + return 0; + } +# endif + } // namespace +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) + + /*! + * Attach wrapper functions to SUNDIALS' linear solver interface. We pretend + * that the user-supplied linear solver is matrix-free, even though it can + * be matrix-based. This way SUNDIALS does not need to understand our matrix + * types. + */ + template + class SundialsLinearSolverWrapper + { + public: + SundialsLinearSolverWrapper(ARKode & solver, + LinearSolveFunction lsolve) + { + sun_linear_solver = SUNLinSolNewEmpty(); + sun_linear_solver->ops->gettype = arkode_linsol_get_type; + sun_linear_solver->ops->solve = arkode_linsol_solve; + sun_linear_solver->ops->setup = arkode_linsol_setup; + sun_linear_solver->ops->initialize = arkode_linsol_initialize; + sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times; + sun_linear_solver->ops->setpreconditioner = + arkode_linsol_set_preconditioner; + + content.solver = &solver; + content.lsolve = lsolve; + sun_linear_solver->content = &content; + } + + ~SundialsLinearSolverWrapper() + { + SUNLinSolFreeEmpty(sun_linear_solver); + } + + SUNLinearSolver + get_wrapped_solver() + { + return sun_linear_solver; + } + + private: + SUNLinearSolver sun_linear_solver; + LinearSolverContent content; + }; + +# endif + template ARKode::ARKode(const AdditionalData &data, const MPI_Comm & mpi_comm) @@ -254,7 +643,12 @@ namespace SUNDIALS ARKode::~ARKode() { if (arkode_mem) +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) ARKodeFree(&arkode_mem); +# else + ARKStepFree(&arkode_mem); +# endif + # ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { @@ -306,9 +700,14 @@ namespace SUNDIALS while (time.is_at_end() == false) { time.set_desired_next_step_size(data.output_period); - double actual_next_time; + double actual_next_time; +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) const auto status = SundialsARKode( arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL); +# else + const auto status = ARKStepEvolve( + arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL); +# endif (void)status; AssertARKode(status); @@ -345,6 +744,7 @@ namespace SUNDIALS return time.get_step_number(); } +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) template void ARKode::reset(const double current_time, @@ -361,24 +761,14 @@ namespace SUNDIALS // Free the vectors which are no longer used. if (yy) { -# ifdef DEAL_II_WITH_MPI - if (is_serial_vector::value == false) - { - N_VDestroy_Parallel(yy); - N_VDestroy_Parallel(abs_tolls); - } - else -# endif - { - N_VDestroy_Serial(yy); - N_VDestroy_Serial(abs_tolls); - } + free_vector(yy); + free_vector(abs_tolls); } int status; (void)status; system_size = solution.size(); -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const IndexSet is = solution.locally_owned_elements(); @@ -390,7 +780,7 @@ namespace SUNDIALS N_VNew_Parallel(communicator, local_system_size, system_size); } else -# endif +# endif { yy = N_VNew_Serial(system_size); abs_tolls = N_VNew_Serial(system_size); @@ -456,9 +846,9 @@ namespace SUNDIALS if (setup_jacobian) { ARKode_mem->ark_lsetup = t_arkode_setup_jacobian; -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) ARKode_mem->ark_setupNonNull = true; -# endif +# endif } } else @@ -476,16 +866,223 @@ namespace SUNDIALS if (setup_mass) { ARKode_mem->ark_msetup = t_arkode_setup_mass; -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) ARKode_mem->ark_MassSetupNonNull = true; -# endif +# endif } } 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) + { + if (arkode_mem) + ARKStepFree(&arkode_mem); + + // Free the vectors which are no longer used. + if (yy) + { + free_vector(yy); + free_vector(abs_tolls); + } + + int status; + (void)status; + yy = create_vector(solution); + abs_tolls = create_vector(solution); + copy(yy, solution); + + Assert(explicit_function || implicit_function, + ExcFunctionNotProvided("explicit_function || implicit_function")); + + arkode_mem = ARKStepCreate( + explicit_function ? &t_arkode_explicit_function : nullptr, + implicit_function ? &t_arkode_implicit_function : nullptr, + current_time, + yy); + + if (get_local_tolerances) + { + copy(abs_tolls, get_local_tolerances()); + status = + ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls); + AssertARKode(status); + } + else + { + status = ARKStepSStolerances(arkode_mem, + data.relative_tolerance, + data.absolute_tolerance); + AssertARKode(status); + } + + // for version 4.0.0 this call must be made before solver settings + status = ARKStepSetUserData(arkode_mem, this); + AssertARKode(status); + + setup_system_solver(solution); + + setup_mass_solver(); + + status = ARKStepSetInitStep(arkode_mem, current_time_step); + AssertARKode(status); + + status = ARKStepSetStopTime(arkode_mem, data.final_time); + AssertARKode(status); + + status = ARKStepSetOrder(arkode_mem, data.maximum_order); + AssertARKode(status); + + if (custom_setup) + custom_setup(arkode_mem); } + + + template + void + ARKode::setup_system_solver(const VectorType &solution) + { + // no point in setting up a solver when there is no linear system + if (!implicit_function) + return; + + int status; + (void)status; + // Initialize solver + // currently only iterative linear solvers are supported + if (jacobian_times_vector) + { + SUNLinearSolver sun_linear_solver; + if (solve_linearized_system) + { + linear_solver = + std::make_unique>( + *this, solve_linearized_system); + sun_linear_solver = linear_solver->get_wrapped_solver(); + } + else + { + // use default solver from SUNDIALS + // TODO give user options + sun_linear_solver = + SUNLinSol_SPGMR(yy, + PREC_NONE, + 0 /*krylov subvectors, 0 uses default*/); + } + status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver, nullptr); + AssertARKode(status); + status = + ARKStepSetJacTimes(arkode_mem, + jacobian_times_setup ? + t_arkode_jac_times_setup_function : + nullptr, + t_arkode_jac_times_vec_function); + AssertARKode(status); + if (jacobian_preconditioner_solve) + { + status = ARKStepSetPreconditioner( + arkode_mem, + jacobian_preconditioner_setup ? + t_arkode_prec_setup_function : + nullptr, + t_arkode_prec_solve_function); + AssertARKode(status); + } + if (data.implicit_function_is_linear) + { + status = ARKStepSetLinear( + arkode_mem, data.implicit_function_is_time_independent ? 0 : 1); + AssertARKode(status); + } + } + else + { + N_Vector y_template = create_vector(solution); + + SUNNonlinearSolver fixed_point_solver = + SUNNonlinSol_FixedPoint(y_template, + data.anderson_acceleration_subspace); + + free_vector(y_template); + status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver); + AssertARKode(status); + } + + status = + ARKStepSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations); + AssertARKode(status); + } + + + + template + void + ARKode::setup_mass_solver() + { + int status; + (void)status; + + if (mass_times_vector) + { + SUNLinearSolver sun_mass_linear_solver; + if (solve_mass) + { + mass_solver = + std::make_unique>( + *this, solve_mass); + sun_mass_linear_solver = mass_solver->get_wrapped_solver(); + } + else + { + sun_mass_linear_solver = + SUNLinSol_SPGMR(yy, + PREC_NONE, + 0 /*krylov subvectors, 0 uses default*/); + } + booleantype mass_time_dependent = + data.mass_is_time_independent ? SUNFALSE : SUNTRUE; + status = ARKStepSetMassLinearSolver(arkode_mem, + sun_mass_linear_solver, + nullptr, + mass_time_dependent); + AssertARKode(status); + status = + ARKStepSetMassTimes(arkode_mem, + mass_times_setup ? + t_arkode_mass_times_setup_function : + nullptr, + t_arkode_mass_times_vec_function, + this); + AssertARKode(status); + + if (mass_preconditioner_solve) + { + status = ARKStepSetMassPreconditioner( + arkode_mem, + mass_preconditioner_setup ? + t_arkode_mass_prec_setup_function : + nullptr, + t_arkode_mass_prec_solve_function); + AssertARKode(status); + } + } + } +# endif + + + template void ARKode::set_functions_to_trigger_an_assert() @@ -499,20 +1096,169 @@ namespace SUNDIALS }; } + + + template + N_Vector + ARKode::create_vector(const VectorType &template_vector) const + { + N_Vector new_vector; + auto system_size = template_vector.size(); + +# ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + { + new_vector = + N_VNew_Parallel(communicator, + template_vector.locally_owned_elements().n_elements(), + system_size); + } + else +# endif + { + new_vector = N_VNew_Serial(system_size); + } + return new_vector; + } + + template + void + ARKode::free_vector(N_Vector vector) const + { +# ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + N_VDestroy_Parallel(vector); + else +# endif + N_VDestroy_Serial(vector); + } + + + + template + void * + ARKode::get_arkode_memory() const + { + return arkode_mem; + } + +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) + + template + SundialsOperator::SundialsOperator(ARKode &solver, + void * A_data, + ATimesFn a_times_fn) + : solver(solver) + , A_data(A_data) + , a_times_fn(a_times_fn) + + { + Assert(a_times_fn != nullptr, ExcInternalError()); + } + + + + template + void + SundialsOperator::vmult(VectorType & dst, + const VectorType &src) const + { + N_Vector sun_dst = solver.create_vector(dst); + N_Vector sun_src = solver.create_vector(src); + copy(sun_src, src); + int status = a_times_fn(A_data, sun_src, sun_dst); + (void)status; + AssertARKode(status); + + copy(dst, sun_dst); + + solver.free_vector(sun_dst); + solver.free_vector(sun_src); + } + + + + template + SundialsPreconditioner::SundialsPreconditioner( + ARKode &solver, + void * P_data, + PSolveFn p_solve_fn, + double tol) + : solver(solver) + , P_data(P_data) + , p_solve_fn(p_solve_fn) + , tol(tol) + {} + + + + template + void + SundialsPreconditioner::vmult(VectorType & dst, + const VectorType &src) const + { + // apply identity preconditioner if nothing else specified + if (!p_solve_fn) + { + dst = src; + return; + } + + N_Vector sun_dst = solver.create_vector(dst); + N_Vector sun_src = solver.create_vector(src); + copy(sun_src, src); + // for custom preconditioners no distinction between left and right + // preconditioning is made + int status = + p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/); + (void)status; + AssertARKode(status); + + copy(dst, sun_dst); + + solver.free_vector(sun_dst); + solver.free_vector(sun_src); + } +# endif + template class ARKode>; template class ARKode>; +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) + template struct SundialsOperator>; + template struct SundialsOperator>; + + template struct SundialsPreconditioner>; + template struct SundialsPreconditioner>; +# endif + # ifdef DEAL_II_WITH_MPI # ifdef DEAL_II_WITH_TRILINOS template class ARKode; template class ARKode; + +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) + template struct SundialsOperator; + template struct SundialsOperator; + + template struct SundialsPreconditioner; + template struct SundialsPreconditioner; +# endif # endif // DEAL_II_WITH_TRILINOS # ifdef DEAL_II_WITH_PETSC # ifndef PETSC_USE_COMPLEX template class ARKode; template class ARKode; + +# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) + template class SundialsOperator; + template class SundialsOperator; + + template class SundialsPreconditioner; + template class SundialsPreconditioner; +# endif # endif // PETSC_USE_COMPLEX # endif // DEAL_II_WITH_PETSC diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index 2e8454aebd..1fe5fc9061 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -20,28 +20,30 @@ #ifdef DEAL_II_WITH_SUNDIALS -# include +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -# include -# ifdef DEAL_II_WITH_TRILINOS -# include -# include -# endif -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif +# include + +# include +# ifdef DEAL_II_WITH_TRILINOS +# include +# include +# endif +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif -# include +# include -# ifdef DEAL_II_SUNDIALS_WITH_IDAS -# include -# else -# include -# endif +# ifdef DEAL_II_SUNDIALS_WITH_IDAS +# include +# else +# include +# endif -# include -# include +# include +# include DEAL_II_NAMESPACE_OPEN @@ -172,14 +174,14 @@ namespace SUNDIALS { if (ida_mem) IDAFree(&ida_mem); -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const int ierr = MPI_Comm_free(&communicator); (void)ierr; AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr)); } -# endif +# endif } @@ -200,7 +202,7 @@ namespace SUNDIALS // The solution is stored in // solution. Here we take only a // view of it. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const IndexSet is = solution.locally_owned_elements(); @@ -216,7 +218,7 @@ namespace SUNDIALS N_VNew_Parallel(communicator, local_system_size, system_size); } else -# endif +# endif { Assert(is_serial_vector::value, ExcInternalError( @@ -254,7 +256,7 @@ namespace SUNDIALS } // Free the vectors which are no longer used. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { N_VDestroy_Parallel(yy); @@ -263,7 +265,7 @@ namespace SUNDIALS N_VDestroy_Parallel(diff_id); } else -# endif +# endif { N_VDestroy_Serial(yy); N_VDestroy_Serial(yp); @@ -293,7 +295,7 @@ namespace SUNDIALS // Free the vectors which are no longer used. if (yy) { -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { N_VDestroy_Parallel(yy); @@ -302,7 +304,7 @@ namespace SUNDIALS N_VDestroy_Parallel(diff_id); } else -# endif +# endif { N_VDestroy_Serial(yy); N_VDestroy_Serial(yp); @@ -314,7 +316,7 @@ namespace SUNDIALS int status; (void)status; system_size = solution.size(); -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const IndexSet is = solution.locally_owned_elements(); @@ -330,7 +332,7 @@ namespace SUNDIALS N_VNew_Parallel(communicator, local_system_size, system_size); } else -# endif +# endif { yy = N_VNew_Serial(system_size); yp = N_VNew_Serial(system_size); @@ -394,9 +396,9 @@ namespace SUNDIALS IDA_mem->ida_lsetup = t_dae_lsetup; IDA_mem->ida_lsolve = t_dae_solve; -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) IDA_mem->ida_setupNonNull = true; -# endif +# endif status = IDASetMaxOrd(ida_mem, data.maximum_order); AssertIDA(status); @@ -490,24 +492,26 @@ namespace SUNDIALS template class IDA>; template class IDA>; -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI -# ifdef DEAL_II_WITH_TRILINOS +# ifdef DEAL_II_WITH_TRILINOS template class IDA; template class IDA; -# 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 class IDA; template class IDA; -# 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 // DEAL_II_WITH_SUNDIALS diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc index ed2f26844b..0bccafdbbc 100644 --- a/source/sundials/kinsol.cc +++ b/source/sundials/kinsol.cc @@ -20,31 +20,33 @@ #ifdef DEAL_II_WITH_SUNDIALS -# include +# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) -# include -# ifdef DEAL_II_WITH_TRILINOS -# include -# include -# endif -# ifdef DEAL_II_WITH_PETSC -# include -# include -# endif +# include + +# include +# ifdef DEAL_II_WITH_TRILINOS +# include +# include +# endif +# ifdef DEAL_II_WITH_PETSC +# include +# include +# endif -# include +# include -# include -# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) -# include -# include -# include -# else -# include -# endif +# include +# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) +# include +# include +# include +# else +# include +# endif -# include -# include +# include +# include DEAL_II_NAMESPACE_OPEN @@ -171,14 +173,14 @@ namespace SUNDIALS { if (kinsol_mem) KINFree(&kinsol_mem); -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const int ierr = MPI_Comm_free(&communicator); (void)ierr; AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr)); } -# endif +# endif } @@ -192,7 +194,7 @@ namespace SUNDIALS // The solution is stored in // solution. Here we take only a // view of it. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { const IndexSet is = initial_guess_and_solution.locally_owned_elements(); @@ -208,7 +210,7 @@ namespace SUNDIALS N_VConst_Parallel(1.e0, f_scale); } else -# endif +# endif { Assert(is_serial_vector::value, ExcInternalError( @@ -267,10 +269,10 @@ namespace SUNDIALS status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error); AssertKINSOL(status); -# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) SUNMatrix J = nullptr; SUNLinearSolver LS = nullptr; -# endif +# endif if (solve_jacobian_system) { @@ -279,20 +281,20 @@ namespace SUNDIALS if (setup_jacobian) { KIN_mem->kin_lsetup = t_kinsol_setup_jacobian; -# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0) KIN_mem->kin_setupNonNull = true; -# endif +# endif } } else { -# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) J = SUNDenseMatrix(system_size, system_size); LS = SUNDenseLinearSolver(u_scale, J); status = KINDlsSetLinearSolver(kinsol_mem, LS, J); -# else +# else status = KINDense(kinsol_mem, system_size); -# endif +# endif AssertKINSOL(status); } @@ -311,7 +313,7 @@ namespace SUNDIALS copy(initial_guess_and_solution, solution); // Free the vectors which are no longer used. -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { N_VDestroy_Parallel(solution); @@ -319,7 +321,7 @@ namespace SUNDIALS N_VDestroy_Parallel(f_scale); } else -# endif +# endif { N_VDestroy_Serial(solution); N_VDestroy_Serial(u_scale); @@ -330,10 +332,10 @@ namespace SUNDIALS status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters); AssertKINSOL(status); -# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) +# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0) SUNMatDestroy(J); SUNLinSolFree(LS); -# endif +# endif KINFree(&kinsol_mem); return static_cast(nniters); @@ -351,24 +353,24 @@ namespace SUNDIALS template class KINSOL>; template class KINSOL>; -# ifdef DEAL_II_WITH_MPI +# ifdef DEAL_II_WITH_MPI -# ifdef DEAL_II_WITH_TRILINOS +# ifdef DEAL_II_WITH_TRILINOS template class KINSOL; template class KINSOL; -# endif +# endif -# ifdef DEAL_II_WITH_PETSC -# ifndef PETSC_USE_COMPLEX +# ifdef DEAL_II_WITH_PETSC +# ifndef PETSC_USE_COMPLEX template class KINSOL; template class KINSOL; +# endif # endif -# endif -# endif +# endif } // namespace SUNDIALS DEAL_II_NAMESPACE_CLOSE - +# endif #endif diff --git a/tests/sundials/harmonic_oscillator_02.cc b/tests/sundials/arkode_01.cc similarity index 90% rename from tests/sundials/harmonic_oscillator_02.cc rename to tests/sundials/arkode_01.cc index 349131e2c8..28c49b0378 100644 --- a/tests/sundials/harmonic_oscillator_02.cc +++ b/tests/sundials/arkode_01.cc @@ -52,7 +52,7 @@ int main(int argc, char **argv) { - std::ofstream out("output"); + initlog(); Utilities::MPI::MPI_InitFinalize mpi_initialization( argc, argv, numbers::invalid_unsigned_int); @@ -66,12 +66,12 @@ main(int argc, char **argv) // Set to true to reset input file. if (false) { - std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_02.prm"); + std::ofstream ofile(SOURCE_DIR "/arkode_01.prm"); prm.print_parameters(ofile, ParameterHandler::ShortText); ofile.close(); } - std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_02.prm"); + std::ifstream ifile(SOURCE_DIR "/arkode_01.prm"); prm.parse_input(ifile); SUNDIALS::ARKode ode(data); @@ -90,7 +90,7 @@ main(int argc, char **argv) ode.output_step = [&](const double t, const VectorType & sol, const unsigned int step_number) -> int { - out << t << " " << sol[0] << " " << sol[1] << std::endl; + deallog << t << " " << sol[0] << " " << sol[1] << std::endl; return 0; }; diff --git a/tests/sundials/harmonic_oscillator_02.prm b/tests/sundials/arkode_01.prm similarity index 100% rename from tests/sundials/harmonic_oscillator_02.prm rename to tests/sundials/arkode_01.prm diff --git a/tests/sundials/arkode_01.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_01.with_sundials.geq.5.4.0.output new file mode 100644 index 0000000000..40a95d1ef5 --- /dev/null +++ b/tests/sundials/arkode_01.with_sundials.geq.5.4.0.output @@ -0,0 +1,128 @@ + +DEAL::0 0 1 +DEAL::0.05 0.0499792 0.99875 +DEAL::0.1 0.0998334 0.995004 +DEAL::0.15 0.149438 0.988771 +DEAL::0.2 0.198669 0.980067 +DEAL::0.25 0.247404 0.968912 +DEAL::0.3 0.295521 0.955336 +DEAL::0.35 0.3429 0.939372 +DEAL::0.4 0.38942 0.92106 +DEAL::0.45 0.434967 0.900447 +DEAL::0.5 0.479426 0.877583 +DEAL::0.55 0.522687 0.852525 +DEAL::0.6 0.564643 0.825335 +DEAL::0.65 0.605188 0.796083 +DEAL::0.7 0.644219 0.764841 +DEAL::0.75 0.68164 0.731688 +DEAL::0.8 0.717356 0.696707 +DEAL::0.85 0.751281 0.659983 +DEAL::0.9 0.78333 0.621606 +DEAL::0.95 0.813421 0.581675 +DEAL::1 0.841478 0.540292 +DEAL::1.05 0.867431 0.49756 +DEAL::1.1 0.891214 0.453587 +DEAL::1.15 0.912769 0.408482 +DEAL::1.2 0.932042 0.362356 +DEAL::1.25 0.948986 0.315322 +DEAL::1.3 0.96356 0.267498 +DEAL::1.35 0.975725 0.219002 +DEAL::1.4 0.985452 0.169958 +DEAL::1.45 0.992716 0.120492 +DEAL::1.5 0.997498 0.0707265 +DEAL::1.55 0.999787 0.0207865 +DEAL::1.6 0.999576 -0.0292043 +DEAL::1.65 0.996867 -0.0791226 +DEAL::1.7 0.991667 -0.128845 +DEAL::1.75 0.983988 -0.178249 +DEAL::1.8 0.973849 -0.227207 +DEAL::1.85 0.961276 -0.275598 +DEAL::1.9 0.946301 -0.323297 +DEAL::1.95 0.928961 -0.370187 +DEAL::2 0.909299 -0.416151 +DEAL::2.05 0.887365 -0.461075 +DEAL::2.1 0.863211 -0.504848 +DEAL::2.15 0.836899 -0.547361 +DEAL::2.2 0.808496 -0.588507 +DEAL::2.25 0.778072 -0.62818 +DEAL::2.3 0.745704 -0.666282 +DEAL::2.35 0.711473 -0.702718 +DEAL::2.4 0.675464 -0.737397 +DEAL::2.45 0.637766 -0.770234 +DEAL::2.5 0.598473 -0.801147 +DEAL::2.55 0.557682 -0.830058 +DEAL::2.6 0.515498 -0.856895 +DEAL::2.65 0.472027 -0.881588 +DEAL::2.7 0.427377 -0.904078 +DEAL::2.75 0.38166 -0.924307 +DEAL::2.8 0.334989 -0.942226 +DEAL::2.85 0.287478 -0.957791 +DEAL::2.9 0.239247 -0.970962 +DEAL::2.95 0.190418 -0.981707 +DEAL::3 0.141114 -0.989997 +DEAL::3.05 0.0914593 -0.995813 +DEAL::3.1 0.041577 -0.99914 +DEAL::3.15 -0.00840899 -0.999969 +DEAL::3.2 -0.0583749 -0.998299 +DEAL::3.25 -0.108197 -0.994133 +DEAL::3.3 -0.157751 -0.987483 +DEAL::3.35 -0.206909 -0.978364 +DEAL::3.4 -0.255548 -0.966801 +DEAL::3.45 -0.303548 -0.952821 +DEAL::3.5 -0.350788 -0.93646 +DEAL::3.55 -0.397151 -0.917758 +DEAL::3.6 -0.442523 -0.896762 +DEAL::3.65 -0.486791 -0.873523 +DEAL::3.7 -0.529842 -0.848102 +DEAL::3.75 -0.571568 -0.82056 +DEAL::3.8 -0.611864 -0.790969 +DEAL::3.85 -0.65063 -0.759401 +DEAL::3.9 -0.68777 -0.725935 +DEAL::3.95 -0.723192 -0.690654 +DEAL::4 -0.756808 -0.653645 +DEAL::4.05 -0.788532 -0.615002 +DEAL::4.1 -0.818284 -0.574823 +DEAL::4.15 -0.84599 -0.533208 +DEAL::4.2 -0.871582 -0.490261 +DEAL::4.25 -0.894994 -0.446089 +DEAL::4.3 -0.916171 -0.4008 +DEAL::4.35 -0.935058 -0.354508 +DEAL::4.4 -0.951608 -0.30733 +DEAL::4.45 -0.965779 -0.259385 +DEAL::4.5 -0.977536 -0.210793 +DEAL::4.55 -0.98685 -0.161675 +DEAL::4.6 -0.993697 -0.112152 +DEAL::4.65 -0.99806 -0.062348 +DEAL::4.7 -0.999929 -0.0123862 +DEAL::4.75 -0.999298 0.0376067 +DEAL::4.8 -0.99617 0.0875044 +DEAL::4.85 -0.990552 0.137182 +DEAL::4.9 -0.982458 0.186516 +DEAL::4.95 -0.971908 0.235384 +DEAL::5 -0.958929 0.283665 +DEAL::5.05 -0.943553 0.331238 +DEAL::5.1 -0.925818 0.377984 +DEAL::5.15 -0.90577 0.423783 +DEAL::5.2 -0.883458 0.468523 +DEAL::5.25 -0.858939 0.512091 +DEAL::5.3 -0.832272 0.554378 +DEAL::5.35 -0.803524 0.595282 +DEAL::5.4 -0.772767 0.634699 +DEAL::5.45 -0.740079 0.672529 +DEAL::5.5 -0.705541 0.708677 +DEAL::5.55 -0.669241 0.743054 +DEAL::5.6 -0.631269 0.775572 +DEAL::5.65 -0.591718 0.806152 +DEAL::5.7 -0.550688 0.834719 +DEAL::5.75 -0.50828 0.861199 +DEAL::5.8 -0.464601 0.885527 +DEAL::5.85 -0.419763 0.907641 +DEAL::5.9 -0.373876 0.927486 +DEAL::5.95 -0.327055 0.945012 +DEAL::6 -0.279416 0.960177 +DEAL::6.05 -0.231078 0.972942 +DEAL::6.1 -0.182162 0.983275 +DEAL::6.15 -0.132791 0.991151 +DEAL::6.2 -0.0830885 0.996549 +DEAL::6.25 -0.0331785 0.999456 +DEAL::6.28 -0.0031845 1 \ No newline at end of file diff --git a/tests/sundials/arkode_01.with_sundials.le.5.4.0.output b/tests/sundials/arkode_01.with_sundials.le.5.4.0.output new file mode 100644 index 0000000000..1f47ae5c6d --- /dev/null +++ b/tests/sundials/arkode_01.with_sundials.le.5.4.0.output @@ -0,0 +1,128 @@ + +DEAL::0 0 1 +DEAL::0.05 0.0499791 0.99875 +DEAL::0.1 0.0998334 0.995004 +DEAL::0.15 0.149438 0.988771 +DEAL::0.2 0.198669 0.980066 +DEAL::0.25 0.247404 0.968911 +DEAL::0.3 0.295516 0.955326 +DEAL::0.35 0.34289 0.939353 +DEAL::0.4 0.38941 0.92104 +DEAL::0.45 0.43496 0.900433 +DEAL::0.5 0.479424 0.877579 +DEAL::0.55 0.522687 0.852524 +DEAL::0.6 0.564637 0.825329 +DEAL::0.65 0.605176 0.796071 +DEAL::0.7 0.644206 0.764828 +DEAL::0.75 0.681631 0.73168 +DEAL::0.8 0.717355 0.696705 +DEAL::0.85 0.751278 0.659982 +DEAL::0.9 0.783303 0.621595 +DEAL::0.95 0.813363 0.581652 +DEAL::1 0.841396 0.540258 +DEAL::1.05 0.867339 0.497522 +DEAL::1.1 0.891131 0.453553 +DEAL::1.15 0.912711 0.408458 +DEAL::1.2 0.932015 0.362345 +DEAL::1.25 0.948983 0.315321 +DEAL::1.3 0.963554 0.267498 +DEAL::1.35 0.975694 0.219003 +DEAL::1.4 0.985391 0.169961 +DEAL::1.45 0.992634 0.120495 +DEAL::1.5 0.997413 0.07073 +DEAL::1.55 0.999716 0.0207894 +DEAL::1.6 0.999533 -0.0292026 +DEAL::1.65 0.996854 -0.0791221 +DEAL::1.7 0.991667 -0.128845 +DEAL::1.75 0.983976 -0.178243 +DEAL::1.8 0.973816 -0.227193 +DEAL::1.85 0.961228 -0.275575 +DEAL::1.9 0.946248 -0.323273 +DEAL::1.95 0.928918 -0.370167 +DEAL::2 0.909275 -0.41614 +DEAL::2.05 0.88736 -0.461072 +DEAL::2.1 0.86321 -0.504847 +DEAL::2.15 0.836886 -0.547347 +DEAL::2.2 0.808467 -0.588477 +DEAL::2.25 0.778033 -0.62814 +DEAL::2.3 0.745665 -0.666242 +DEAL::2.35 0.711444 -0.702689 +DEAL::2.4 0.675451 -0.737384 +DEAL::2.45 0.637765 -0.770233 +DEAL::2.5 0.598471 -0.801142 +DEAL::2.55 0.557673 -0.830034 +DEAL::2.6 0.515481 -0.856851 +DEAL::2.65 0.472005 -0.881535 +DEAL::2.7 0.427357 -0.904029 +DEAL::2.75 0.381646 -0.924275 +DEAL::2.8 0.334983 -0.942215 +DEAL::2.85 0.287479 -0.957791 +DEAL::2.9 0.239248 -0.970952 +DEAL::2.95 0.190419 -0.981674 +DEAL::3 0.141115 -0.989946 +DEAL::3.05 0.091459 -0.995756 +DEAL::3.1 0.0415761 -0.999091 +DEAL::3.15 -0.00841 -0.999941 +DEAL::3.2 -0.0583754 -0.998292 +DEAL::3.25 -0.108196 -0.994133 +DEAL::3.3 -0.157743 -0.987468 +DEAL::3.35 -0.206893 -0.978328 +DEAL::3.4 -0.255528 -0.96675 +DEAL::3.45 -0.303527 -0.952769 +DEAL::3.5 -0.350773 -0.93642 +DEAL::3.55 -0.397144 -0.917738 +DEAL::3.6 -0.442522 -0.89676 +DEAL::3.65 -0.486788 -0.873522 +DEAL::3.7 -0.529827 -0.848088 +DEAL::3.75 -0.571542 -0.820533 +DEAL::3.8 -0.611833 -0.790935 +DEAL::3.85 -0.650604 -0.759371 +DEAL::3.9 -0.687755 -0.725916 +DEAL::3.95 -0.723188 -0.690649 +DEAL::4 -0.756806 -0.653646 +DEAL::4.05 -0.788517 -0.614997 +DEAL::4.1 -0.818254 -0.574809 +DEAL::4.15 -0.84595 -0.533188 +DEAL::4.2 -0.871543 -0.490241 +DEAL::4.25 -0.894968 -0.446074 +DEAL::4.3 -0.916161 -0.400795 +DEAL::4.35 -0.935058 -0.35451 +DEAL::4.4 -0.951598 -0.307332 +DEAL::4.45 -0.965751 -0.259384 +DEAL::4.5 -0.977494 -0.210788 +DEAL::4.55 -0.986805 -0.161669 +DEAL::4.6 -0.993662 -0.112147 +DEAL::4.65 -0.998044 -0.0623455 +DEAL::4.7 -0.999928 -0.0123877 +DEAL::4.75 -0.999294 0.0376024 +DEAL::4.8 -0.996148 0.0874965 +DEAL::4.85 -0.990514 0.137172 +DEAL::4.9 -0.982413 0.186506 +DEAL::4.95 -0.97187 0.235376 +DEAL::5 -0.958908 0.28366 +DEAL::5.05 -0.94355 0.331236 +DEAL::5.1 -0.925819 0.37798 +DEAL::5.15 -0.905761 0.423774 +DEAL::5.2 -0.883439 0.468508 +DEAL::5.25 -0.858915 0.512075 +DEAL::5.3 -0.832253 0.554367 +DEAL::5.35 -0.803517 0.595277 +DEAL::5.4 -0.772769 0.634698 +DEAL::5.45 -0.740077 0.672522 +DEAL::5.5 -0.705531 0.708661 +DEAL::5.55 -0.669225 0.743031 +DEAL::5.6 -0.631253 0.775551 +DEAL::5.65 -0.591709 0.80614 +DEAL::5.7 -0.550686 0.834716 +DEAL::5.75 -0.508281 0.861197 +DEAL::5.8 -0.464599 0.885514 +DEAL::5.85 -0.419756 0.907617 +DEAL::5.9 -0.373868 0.927459 +DEAL::5.95 -0.327048 0.944992 +DEAL::6 -0.279413 0.960169 +DEAL::6.05 -0.231078 0.972942 +DEAL::6.1 -0.182163 0.983273 +DEAL::6.15 -0.132791 0.991144 +DEAL::6.2 -0.0830885 0.996543 +DEAL::6.25 -0.0331784 0.999455 +DEAL::6.28 -0.00318445 1 \ No newline at end of file diff --git a/tests/sundials/harmonic_oscillator_03.cc b/tests/sundials/arkode_02.cc similarity index 90% rename from tests/sundials/harmonic_oscillator_03.cc rename to tests/sundials/arkode_02.cc index 80961bca48..afae0d877b 100644 --- a/tests/sundials/harmonic_oscillator_03.cc +++ b/tests/sundials/arkode_02.cc @@ -52,7 +52,7 @@ int main(int argc, char **argv) { - std::ofstream out("output"); + initlog(); Utilities::MPI::MPI_InitFinalize mpi_initialization( argc, argv, numbers::invalid_unsigned_int); @@ -63,8 +63,8 @@ main(int argc, char **argv) SUNDIALS::ARKode::AdditionalData data; data.add_parameters(prm); - // Use the same parameters of test 2. - std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_02.prm"); + // Use the same parameters of test 1. + std::ifstream ifile(SOURCE_DIR "/arkode_01.prm"); prm.parse_input(ifile); SUNDIALS::ARKode ode(data); @@ -86,8 +86,8 @@ main(int argc, char **argv) // limit the output to every 10th step and increase the precision to make // the test more robust if (step_number % 10 == 0) - out << t << " " << std::setprecision(7) << sol[0] << " " << sol[1] - << std::endl; + deallog << t << " " << std::setprecision(7) << sol[0] << " " << sol[1] + << std::endl; return 0; }; diff --git a/tests/sundials/arkode_02.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_02.with_sundials.geq.5.4.0.output new file mode 100644 index 0000000000..8100d2ca42 --- /dev/null +++ b/tests/sundials/arkode_02.with_sundials.geq.5.4.0.output @@ -0,0 +1,14 @@ + +DEAL::0 0 1 +DEAL::0.5 0.4794368 0.877573 +DEAL::1 0.8414828 0.5402389 +DEAL::1.5 0.9974885 0.07066302 +DEAL::2 0.9091538 -0.416307 +DEAL::2.5 0.5984422 -0.8010988 +DEAL::3 0.1410421 -0.9899475 +DEAL::3.5 -0.3508891 -0.9363109 +DEAL::4 -0.756723 -0.6536002 +DEAL::4.5 -0.9774401 -0.2107849 +DEAL::5 -0.9588134 0.2836821 +DEAL::5.5 -0.7054196 0.7086219 +DEAL::6 -0.2794164 0.9600536 \ No newline at end of file diff --git a/tests/sundials/arkode_02.with_sundials.le.4.0.0.output b/tests/sundials/arkode_02.with_sundials.le.4.0.0.output new file mode 100644 index 0000000000..67768818ec --- /dev/null +++ b/tests/sundials/arkode_02.with_sundials.le.4.0.0.output @@ -0,0 +1,14 @@ + +DEAL::0 0 1 +DEAL::0.5 0.4793775 0.8775139 +DEAL::1 0.8412554 0.540225 +DEAL::1.5 0.9970436 0.07063589 +DEAL::2 0.9084636 -0.4157026 +DEAL::2.5 0.5984372 -0.8010479 +DEAL::3 0.1409945 -0.9894479 +DEAL::3.5 -0.350428 -0.9357101 +DEAL::4 -0.7566852 -0.6535909 +DEAL::4.5 -0.9772954 -0.21075 +DEAL::5 -0.9583971 0.2834953 +DEAL::5.5 -0.7052272 0.7082618 +DEAL::6 -0.279415 0.9600362 \ No newline at end of file diff --git a/tests/sundials/arkode_02.with_sundials=4.0.0.output b/tests/sundials/arkode_02.with_sundials=4.0.0.output new file mode 100644 index 0000000000..643ce6ea39 --- /dev/null +++ b/tests/sundials/arkode_02.with_sundials=4.0.0.output @@ -0,0 +1,14 @@ + +DEAL::0.00000 0.000000 1.000000 +DEAL::0.5000000 0.4793775 0.8775139 +DEAL::1.000000 0.8412555 0.5402250 +DEAL::1.500000 0.9970430 0.07063579 +DEAL::2.000000 0.9084635 -0.4157024 +DEAL::2.500000 0.5984374 -0.8010482 +DEAL::3.000000 0.1409943 -0.9894467 +DEAL::3.500000 -0.3504279 -0.9357104 +DEAL::4.000000 -0.7566856 -0.6535912 +DEAL::4.500000 -0.9772941 -0.2107495 +DEAL::5.000000 -0.9583959 0.2834948 +DEAL::5.500000 -0.7052279 0.7082622 +DEAL::6.000000 -0.2794151 0.9600365 diff --git a/tests/sundials/harmonic_oscillator_04.cc b/tests/sundials/arkode_03.cc similarity index 91% rename from tests/sundials/harmonic_oscillator_04.cc rename to tests/sundials/arkode_03.cc index 882a95bff6..db0a27ce60 100644 --- a/tests/sundials/harmonic_oscillator_04.cc +++ b/tests/sundials/arkode_03.cc @@ -49,7 +49,7 @@ int main(int argc, char **argv) { - std::ofstream out("output"); + initlog(); Utilities::MPI::MPI_InitFinalize mpi_initialization( argc, argv, numbers::invalid_unsigned_int); @@ -62,12 +62,12 @@ main(int argc, char **argv) if (false) { - std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_04.prm"); + std::ofstream ofile(SOURCE_DIR "/arkode_03.prm"); prm.print_parameters(ofile, ParameterHandler::ShortText); ofile.close(); } - std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_04.prm"); + std::ifstream ifile(SOURCE_DIR "/arkode_03.prm"); prm.parse_input(ifile); SUNDIALS::ARKode ode(data); @@ -103,8 +103,8 @@ main(int argc, char **argv) // limit the output to every 10th step and increase the precision to make // the test more robust if (step_number % 10 == 0) - out << t << " " << std::setprecision(10) << sol[0] << " " << sol[1] << " " - << sol[2] << std::endl; + deallog << t << " " << std::setprecision(10) << sol[0] << " " << sol[1] + << " " << sol[2] << std::endl; return 0; }; diff --git a/tests/sundials/harmonic_oscillator_04.prm b/tests/sundials/arkode_03.prm similarity index 100% rename from tests/sundials/harmonic_oscillator_04.prm rename to tests/sundials/arkode_03.prm diff --git a/tests/sundials/arkode_03.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_03.with_sundials.geq.5.4.0.output new file mode 100644 index 0000000000..9e2ed6e8a6 --- /dev/null +++ b/tests/sundials/arkode_03.with_sundials.geq.5.4.0.output @@ -0,0 +1,12 @@ + +DEAL::0 3.9 1.1 2.8 +DEAL::1 2.075274578 1.059601971 2.499947283 +DEAL::2 1.100460018 1.723610111 2.499972388 +DEAL::3 0.7896202372 2.326668043 2.499979776 +DEAL::4 0.8117726433 2.725583965 2.499987084 +DEAL::5 1.187002916 2.595541878 2.499973736 +DEAL::6 1.887311324 1.494098057 2.499952698 +DEAL::7 1.323271341 1.619016372 2.499965777 +DEAL::8 0.9217874641 2.130219386 2.499970637 +DEAL::9 0.8495425771 2.539475069 2.499991366 +DEAL::10 1.064969831 2.595955794 2.499973383 \ No newline at end of file diff --git a/tests/sundials/arkode_03.with_sundials.le.5.4.0.output b/tests/sundials/arkode_03.with_sundials.le.5.4.0.output new file mode 100644 index 0000000000..16be7ac6e6 --- /dev/null +++ b/tests/sundials/arkode_03.with_sundials.le.5.4.0.output @@ -0,0 +1,12 @@ + +DEAL::0 3.9 1.1 2.8 +DEAL::1 2.075274578 1.059601971 2.499948149 +DEAL::2 1.100460018 1.723610111 2.499972517 +DEAL::3 0.7896202372 2.326668043 2.499980207 +DEAL::4 0.8117726432 2.725583965 2.499979833 +DEAL::5 1.187002916 2.595541878 2.499970301 +DEAL::6 1.887311324 1.494098057 2.499952687 +DEAL::7 1.323271341 1.619016372 2.499966912 +DEAL::8 0.9217874641 2.130219386 2.499976935 +DEAL::9 0.849542577 2.539475069 2.499978894 +DEAL::10 1.064969831 2.595955794 2.499973382 \ No newline at end of file diff --git a/tests/sundials/harmonic_oscillator_05.cc b/tests/sundials/arkode_04.cc similarity index 93% rename from tests/sundials/harmonic_oscillator_05.cc rename to tests/sundials/arkode_04.cc index bb32408f62..3198c8a57e 100644 --- a/tests/sundials/harmonic_oscillator_05.cc +++ b/tests/sundials/arkode_04.cc @@ -49,7 +49,7 @@ int main(int argc, char **argv) { - std::ofstream out("output"); + initlog(); Utilities::MPI::MPI_InitFinalize mpi_initialization( argc, argv, numbers::invalid_unsigned_int); @@ -62,12 +62,12 @@ main(int argc, char **argv) if (false) { - std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_05.prm"); + std::ofstream ofile(SOURCE_DIR "/arkode_04.prm"); prm.print_parameters(ofile, ParameterHandler::ShortText); ofile.close(); } - std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_05.prm"); + std::ifstream ifile(SOURCE_DIR "/arkode_04.prm"); prm.parse_input(ifile); SUNDIALS::ARKode ode(data); @@ -117,7 +117,8 @@ main(int argc, char **argv) ode.output_step = [&](const double t, const VectorType & sol, const unsigned int step_number) -> int { - out << t << " " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl; + deallog << t << " " << sol[0] << " " << sol[1] << " " << sol[2] + << std::endl; return 0; }; diff --git a/tests/sundials/harmonic_oscillator_05.prm b/tests/sundials/arkode_04.prm similarity index 100% rename from tests/sundials/harmonic_oscillator_05.prm rename to tests/sundials/arkode_04.prm diff --git a/tests/sundials/arkode_04.with_sundials.le.4.0.0.output b/tests/sundials/arkode_04.with_sundials.le.4.0.0.output new file mode 100644 index 0000000000..890ce32ab9 --- /dev/null +++ b/tests/sundials/arkode_04.with_sundials.le.4.0.0.output @@ -0,0 +1,102 @@ + +DEAL::0 3.9 1.1 2.8 +DEAL::0.1 3.99995 0.720431 2.49989 +DEAL::0.2 3.78894 0.660967 2.49996 +DEAL::0.3 3.52275 0.681501 2.49984 +DEAL::0.4 3.26217 0.72296 2.49995 +DEAL::0.5 3.01935 0.771864 2.49995 +DEAL::0.6 2.79594 0.82467 2.49988 +DEAL::0.7 2.59116 0.880241 2.49998 +DEAL::0.8 2.4037 0.938097 2.49993 +DEAL::0.9 2.23217 0.997959 2.49992 +DEAL::1 2.07527 1.0596 2.49999 +DEAL::1.1 1.93182 1.12281 2.49992 +DEAL::1.2 1.80073 1.18737 2.49995 +DEAL::1.3 1.68106 1.25304 2.49999 +DEAL::1.4 1.57194 1.31959 2.49993 +DEAL::1.5 1.47261 1.38677 2.49997 +DEAL::1.6 1.38237 1.45434 2.49999 +DEAL::1.7 1.3006 1.52203 2.49994 +DEAL::1.8 1.2267 1.58963 2.49998 +DEAL::1.9 1.16015 1.65689 2.49998 +DEAL::2 1.10046 1.72361 2.49996 +DEAL::2.1 1.04715 1.78959 2.49998 +DEAL::2.2 0.999783 1.85466 2.49997 +DEAL::2.3 0.957942 1.91866 2.49997 +DEAL::2.4 0.921233 1.98145 2.49999 +DEAL::2.5 0.889283 2.04291 2.49997 +DEAL::2.6 0.861743 2.10294 2.49998 +DEAL::2.7 0.838291 2.16142 2.49999 +DEAL::2.8 0.818628 2.21827 2.49997 +DEAL::2.9 0.802485 2.27338 2.49998 +DEAL::3 0.78962 2.32667 2.49998 +DEAL::3.1 0.779824 2.37802 2.49998 +DEAL::3.2 0.772919 2.42731 2.49998 +DEAL::3.3 0.768758 2.47441 2.49999 +DEAL::3.4 0.767229 2.51916 2.49998 +DEAL::3.5 0.768257 2.56138 2.5 +DEAL::3.6 0.7718 2.60085 2.49995 +DEAL::3.7 0.777859 2.63733 2.50005 +DEAL::3.8 0.786474 2.67052 2.49999 +DEAL::3.9 0.797733 2.70008 2.49988 +DEAL::4 0.811773 2.72558 2.5001 +DEAL::4.1 0.828789 2.74657 2.5001 +DEAL::4.2 0.84904 2.76245 2.50009 +DEAL::4.3 0.872862 2.77257 2.50012 +DEAL::4.4 0.900673 2.77611 2.49986 +DEAL::4.5 0.932987 2.77216 2.50009 +DEAL::4.6 0.970423 2.75959 2.49999 +DEAL::4.7 1.01371 2.73716 2.50002 +DEAL::4.8 1.06366 2.70339 2.50006 +DEAL::4.9 1.12115 2.65673 2.49999 +DEAL::5 1.187 2.59554 2.49993 +DEAL::5.1 1.26179 2.51839 2.50005 +DEAL::5.2 1.34554 2.42435 2.49986 +DEAL::5.3 1.43719 2.31362 2.50002 +DEAL::5.4 1.53404 2.18823 2.50001 +DEAL::5.5 1.6313 2.05269 2.4999 +DEAL::5.6 1.72208 1.91414 2.49987 +DEAL::5.7 1.79854 1.7815 2.49991 +DEAL::5.8 1.85371 1.66352 2.49995 +DEAL::5.9 1.88345 1.56671 2.49996 +DEAL::6 1.88731 1.4941 2.49995 +DEAL::6.1 1.86803 1.44543 2.49992 +DEAL::6.2 1.83015 1.41827 2.49989 +DEAL::6.3 1.77867 1.40922 2.50001 +DEAL::6.4 1.71817 1.41482 2.50005 +DEAL::6.5 1.65242 1.43201 2.50005 +DEAL::6.6 1.58434 1.45824 2.50005 +DEAL::6.7 1.51607 1.49149 2.49991 +DEAL::6.8 1.44916 1.53017 2.4999 +DEAL::6.9 1.38465 1.57301 2.50006 +DEAL::7 1.32327 1.61902 2.49995 +DEAL::7.1 1.26548 1.6674 2.49989 +DEAL::7.2 1.21156 1.7175 2.50005 +DEAL::7.3 1.16164 1.76879 2.49998 +DEAL::7.4 1.11578 1.82081 2.49991 +DEAL::7.5 1.07396 1.87319 2.50002 +DEAL::7.6 1.03609 1.92558 2.5 +DEAL::7.7 1.00208 1.97772 2.49993 +DEAL::7.8 0.971784 2.02935 2.49999 +DEAL::7.9 0.945072 2.08025 2.5 +DEAL::8 0.921787 2.13022 2.49995 +DEAL::8.1 0.901781 2.17907 2.49997 +DEAL::8.2 0.884906 2.22664 2.5 +DEAL::8.3 0.871025 2.27275 2.49997 +DEAL::8.4 0.860013 2.31723 2.49997 +DEAL::8.5 0.851761 2.35992 2.49999 +DEAL::8.6 0.84618 2.40062 2.49998 +DEAL::8.7 0.843199 2.43916 2.49998 +DEAL::8.8 0.842774 2.47531 2.49999 +DEAL::8.9 0.844885 2.50883 2.49997 +DEAL::9 0.849543 2.53948 2.50004 +DEAL::9.1 0.856788 2.56694 2.49994 +DEAL::9.2 0.866696 2.59088 2.49999 +DEAL::9.3 0.879382 2.61091 2.50009 +DEAL::9.4 0.895002 2.6266 2.49992 +DEAL::9.5 0.913755 2.63743 2.49985 +DEAL::9.6 0.935893 2.64284 2.49987 +DEAL::9.7 0.961713 2.64217 2.50002 +DEAL::9.8 0.991567 2.63469 2.49996 +DEAL::9.9 1.02585 2.61958 2.50009 +DEAL::10 1.06497 2.59596 2.49997 \ No newline at end of file diff --git a/tests/sundials/harmonic_oscillator_06.cc b/tests/sundials/arkode_05.cc similarity index 93% rename from tests/sundials/harmonic_oscillator_06.cc rename to tests/sundials/arkode_05.cc index 5fb4acbda3..40339066bd 100644 --- a/tests/sundials/harmonic_oscillator_06.cc +++ b/tests/sundials/arkode_05.cc @@ -49,7 +49,7 @@ int main(int argc, char **argv) { - std::ofstream out("output"); + initlog(); Utilities::MPI::MPI_InitFinalize mpi_initialization( argc, argv, numbers::invalid_unsigned_int); @@ -62,12 +62,12 @@ main(int argc, char **argv) if (false) { - std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_06.prm"); + std::ofstream ofile(SOURCE_DIR "/arkode_05.prm"); prm.print_parameters(ofile, ParameterHandler::ShortText); ofile.close(); } - std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_06.prm"); + std::ifstream ifile(SOURCE_DIR "/arkode_05.prm"); prm.parse_input(ifile); SUNDIALS::ARKode ode(data); @@ -128,7 +128,8 @@ main(int argc, char **argv) ode.output_step = [&](const double t, const VectorType & sol, const unsigned int step_number) -> int { - out << t << " " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl; + deallog << t << " " << sol[0] << " " << sol[1] << " " << sol[2] + << std::endl; return 0; }; diff --git a/tests/sundials/harmonic_oscillator_06.prm b/tests/sundials/arkode_05.prm similarity index 100% rename from tests/sundials/harmonic_oscillator_06.prm rename to tests/sundials/arkode_05.prm diff --git a/tests/sundials/arkode_05.with_sundials.le.4.0.0.output b/tests/sundials/arkode_05.with_sundials.le.4.0.0.output new file mode 100644 index 0000000000..890ce32ab9 --- /dev/null +++ b/tests/sundials/arkode_05.with_sundials.le.4.0.0.output @@ -0,0 +1,102 @@ + +DEAL::0 3.9 1.1 2.8 +DEAL::0.1 3.99995 0.720431 2.49989 +DEAL::0.2 3.78894 0.660967 2.49996 +DEAL::0.3 3.52275 0.681501 2.49984 +DEAL::0.4 3.26217 0.72296 2.49995 +DEAL::0.5 3.01935 0.771864 2.49995 +DEAL::0.6 2.79594 0.82467 2.49988 +DEAL::0.7 2.59116 0.880241 2.49998 +DEAL::0.8 2.4037 0.938097 2.49993 +DEAL::0.9 2.23217 0.997959 2.49992 +DEAL::1 2.07527 1.0596 2.49999 +DEAL::1.1 1.93182 1.12281 2.49992 +DEAL::1.2 1.80073 1.18737 2.49995 +DEAL::1.3 1.68106 1.25304 2.49999 +DEAL::1.4 1.57194 1.31959 2.49993 +DEAL::1.5 1.47261 1.38677 2.49997 +DEAL::1.6 1.38237 1.45434 2.49999 +DEAL::1.7 1.3006 1.52203 2.49994 +DEAL::1.8 1.2267 1.58963 2.49998 +DEAL::1.9 1.16015 1.65689 2.49998 +DEAL::2 1.10046 1.72361 2.49996 +DEAL::2.1 1.04715 1.78959 2.49998 +DEAL::2.2 0.999783 1.85466 2.49997 +DEAL::2.3 0.957942 1.91866 2.49997 +DEAL::2.4 0.921233 1.98145 2.49999 +DEAL::2.5 0.889283 2.04291 2.49997 +DEAL::2.6 0.861743 2.10294 2.49998 +DEAL::2.7 0.838291 2.16142 2.49999 +DEAL::2.8 0.818628 2.21827 2.49997 +DEAL::2.9 0.802485 2.27338 2.49998 +DEAL::3 0.78962 2.32667 2.49998 +DEAL::3.1 0.779824 2.37802 2.49998 +DEAL::3.2 0.772919 2.42731 2.49998 +DEAL::3.3 0.768758 2.47441 2.49999 +DEAL::3.4 0.767229 2.51916 2.49998 +DEAL::3.5 0.768257 2.56138 2.5 +DEAL::3.6 0.7718 2.60085 2.49995 +DEAL::3.7 0.777859 2.63733 2.50005 +DEAL::3.8 0.786474 2.67052 2.49999 +DEAL::3.9 0.797733 2.70008 2.49988 +DEAL::4 0.811773 2.72558 2.5001 +DEAL::4.1 0.828789 2.74657 2.5001 +DEAL::4.2 0.84904 2.76245 2.50009 +DEAL::4.3 0.872862 2.77257 2.50012 +DEAL::4.4 0.900673 2.77611 2.49986 +DEAL::4.5 0.932987 2.77216 2.50009 +DEAL::4.6 0.970423 2.75959 2.49999 +DEAL::4.7 1.01371 2.73716 2.50002 +DEAL::4.8 1.06366 2.70339 2.50006 +DEAL::4.9 1.12115 2.65673 2.49999 +DEAL::5 1.187 2.59554 2.49993 +DEAL::5.1 1.26179 2.51839 2.50005 +DEAL::5.2 1.34554 2.42435 2.49986 +DEAL::5.3 1.43719 2.31362 2.50002 +DEAL::5.4 1.53404 2.18823 2.50001 +DEAL::5.5 1.6313 2.05269 2.4999 +DEAL::5.6 1.72208 1.91414 2.49987 +DEAL::5.7 1.79854 1.7815 2.49991 +DEAL::5.8 1.85371 1.66352 2.49995 +DEAL::5.9 1.88345 1.56671 2.49996 +DEAL::6 1.88731 1.4941 2.49995 +DEAL::6.1 1.86803 1.44543 2.49992 +DEAL::6.2 1.83015 1.41827 2.49989 +DEAL::6.3 1.77867 1.40922 2.50001 +DEAL::6.4 1.71817 1.41482 2.50005 +DEAL::6.5 1.65242 1.43201 2.50005 +DEAL::6.6 1.58434 1.45824 2.50005 +DEAL::6.7 1.51607 1.49149 2.49991 +DEAL::6.8 1.44916 1.53017 2.4999 +DEAL::6.9 1.38465 1.57301 2.50006 +DEAL::7 1.32327 1.61902 2.49995 +DEAL::7.1 1.26548 1.6674 2.49989 +DEAL::7.2 1.21156 1.7175 2.50005 +DEAL::7.3 1.16164 1.76879 2.49998 +DEAL::7.4 1.11578 1.82081 2.49991 +DEAL::7.5 1.07396 1.87319 2.50002 +DEAL::7.6 1.03609 1.92558 2.5 +DEAL::7.7 1.00208 1.97772 2.49993 +DEAL::7.8 0.971784 2.02935 2.49999 +DEAL::7.9 0.945072 2.08025 2.5 +DEAL::8 0.921787 2.13022 2.49995 +DEAL::8.1 0.901781 2.17907 2.49997 +DEAL::8.2 0.884906 2.22664 2.5 +DEAL::8.3 0.871025 2.27275 2.49997 +DEAL::8.4 0.860013 2.31723 2.49997 +DEAL::8.5 0.851761 2.35992 2.49999 +DEAL::8.6 0.84618 2.40062 2.49998 +DEAL::8.7 0.843199 2.43916 2.49998 +DEAL::8.8 0.842774 2.47531 2.49999 +DEAL::8.9 0.844885 2.50883 2.49997 +DEAL::9 0.849543 2.53948 2.50004 +DEAL::9.1 0.856788 2.56694 2.49994 +DEAL::9.2 0.866696 2.59088 2.49999 +DEAL::9.3 0.879382 2.61091 2.50009 +DEAL::9.4 0.895002 2.6266 2.49992 +DEAL::9.5 0.913755 2.63743 2.49985 +DEAL::9.6 0.935893 2.64284 2.49987 +DEAL::9.7 0.961713 2.64217 2.50002 +DEAL::9.8 0.991567 2.63469 2.49996 +DEAL::9.9 1.02585 2.61958 2.50009 +DEAL::10 1.06497 2.59596 2.49997 \ No newline at end of file diff --git a/tests/sundials/arkode_06.cc b/tests/sundials/arkode_06.cc new file mode 100644 index 0000000000..8bd82b01ba --- /dev/null +++ b/tests/sundials/arkode_06.cc @@ -0,0 +1,159 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +//----------------------------------------------------------- + +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + + +// Test implicit-explicit time stepper. jac_times_vector/setup + +// custom linear solver + +/** + * This test problem is called "brusselator", and is a typical benchmark for + * ODE solvers. This problem has 3 dependent variables u, v and w, that depend + * on the independent variable t via the IVP system + * + * du/dt = a − (w + 1)u + v u^2 + * dv/dt = w u − v u^2 + * dw/dt = (b − w)/eps -w u + * + * We integrate over the interval 0 ≤ t ≤ 10, with the initial conditions + * + * u(0) = 3.9, v(0) = 1.1, w(0) = 2.8, + * + * and parameters + * + * a = 1.2, b = 2.5, and eps = 10−5 + * + * The implicit part only contains the stiff part of the problem (the part with + * eps in right hand side of the third equation). + */ +int +main(int argc, char **argv) +{ + initlog(); + // restrict output to highest level + deallog.depth_file(1); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + typedef Vector VectorType; + + ParameterHandler prm; + SUNDIALS::ARKode::AdditionalData data; + data.add_parameters(prm); + + if (false) + { + std::ofstream ofile(SOURCE_DIR "/arkode_06.prm"); + prm.print_parameters(ofile, ParameterHandler::ShortText); + ofile.close(); + } + + std::ifstream ifile(SOURCE_DIR "/arkode_06.prm"); + prm.parse_input(ifile); + + SUNDIALS::ARKode ode(data); + + ode.reinit_vector = [&](VectorType &v) { + // Three independent variables + v.reinit(3); + }; + + // Parameters + double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b = 2.5, eps = 1e-5; + // Explicit jacobian. + FullMatrix J(3, 3); + + ode.implicit_function = + [&](double, const VectorType &y, VectorType &ydot) -> int { + ydot[0] = 0; + ydot[1] = 0; + ydot[2] = (b - y[2]) / eps; + return 0; + }; + + + ode.explicit_function = + [&](double, const VectorType &y, VectorType &ydot) -> int { + ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; + ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; + ydot[2] = -y[2] * y[0]; + return 0; + }; + + + ode.jacobian_times_setup = + [&](realtype t, const VectorType &y, const VectorType &fy) -> int { + J = 0; + J(2, 2) = -1.0 / eps; + return 0; + }; + + ode.jacobian_times_vector = [&](const VectorType &v, + VectorType & Jv, + double t, + const VectorType &y, + const VectorType &fy) -> int { + J.vmult(Jv, v); + return 0; + }; + + ode.solve_linearized_system = + [&](SUNDIALS::SundialsOperator &op, + SUNDIALS::SundialsPreconditioner &, + VectorType & x, + const VectorType &b, + double tol) -> int { + ReductionControl control; + SolverCG solver_cg(control); + solver_cg.solve(op, x, b, PreconditionIdentity()); + return 0; + }; + + ode.output_step = [&](const double t, + const VectorType & sol, + const unsigned int step_number) -> int { + deallog << t << " " << sol[0] << " " << sol[1] << " " << sol[2] + << std::endl; + return 0; + }; + + // after 5.2.0 a special interpolation mode should be used for stiff problems +#if DEAL_II_SUNDIALS_VERSION_GTE(5, 2, 0) + ode.custom_setup = [&](void *arkode_mem) { + ARKStepSetInterpolantType(arkode_mem, ARK_INTERP_LAGRANGE); + }; +#endif + + Vector y(3); + y[0] = u0; + y[1] = v0; + y[2] = w0; + ode.solve_ode(y); + return 0; +} diff --git a/tests/sundials/arkode_06.prm b/tests/sundials/arkode_06.prm new file mode 100644 index 0000000000..55a91a0b7c --- /dev/null +++ b/tests/sundials/arkode_06.prm @@ -0,0 +1,15 @@ +set Final time = 10. +set Initial time = 0. +set Time interval between each output = 0.1 +subsection Error control + set Absolute error tolerance = 0.000001 + set Relative error tolerance = 0.000010 +end +subsection Running parameters + set Implicit function is linear = true + set Implicit function is time independent = true + set Initial step size = 0.010000 + set Maximum number of nonlinear iterations = 10 + set Maximum order of ARK = 5 + set Minimum step size = 0.000001 +end diff --git a/tests/sundials/arkode_06.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_06.with_sundials.geq.5.4.0.output new file mode 100644 index 0000000000..7c605fb0b8 --- /dev/null +++ b/tests/sundials/arkode_06.with_sundials.geq.5.4.0.output @@ -0,0 +1,102 @@ + +DEAL::0.00000 3.90000 1.10000 2.80000 +DEAL::0.100000 3.99995 0.720432 2.49990 +DEAL::0.200000 3.78894 0.660967 2.49990 +DEAL::0.300000 3.52275 0.681501 2.49991 +DEAL::0.400000 3.26217 0.722960 2.49992 +DEAL::0.500000 3.01935 0.771864 2.49992 +DEAL::0.600000 2.79594 0.824670 2.49993 +DEAL::0.700000 2.59116 0.880241 2.49993 +DEAL::0.800000 2.40370 0.938097 2.49994 +DEAL::0.900000 2.23217 0.997959 2.49994 +DEAL::1.00000 2.07527 1.05960 2.49995 +DEAL::1.10000 1.93182 1.12281 2.49995 +DEAL::1.20000 1.80073 1.18737 2.49995 +DEAL::1.30000 1.68106 1.25304 2.49996 +DEAL::1.40000 1.57194 1.31959 2.49996 +DEAL::1.50000 1.47261 1.38677 2.49996 +DEAL::1.60000 1.38237 1.45434 2.49997 +DEAL::1.70000 1.30060 1.52203 2.49997 +DEAL::1.80000 1.22670 1.58963 2.49997 +DEAL::1.90000 1.16015 1.65689 2.49997 +DEAL::2.00000 1.10046 1.72361 2.49997 +DEAL::2.10000 1.04715 1.78959 2.49997 +DEAL::2.20000 0.999783 1.85466 2.49997 +DEAL::2.30000 0.957942 1.91866 2.49998 +DEAL::2.40000 0.921233 1.98145 2.49998 +DEAL::2.50000 0.889283 2.04291 2.49998 +DEAL::2.60000 0.861743 2.10294 2.49998 +DEAL::2.70000 0.838291 2.16142 2.49998 +DEAL::2.80000 0.818628 2.21827 2.49998 +DEAL::2.90000 0.802485 2.27338 2.49998 +DEAL::3.00000 0.789620 2.32667 2.49998 +DEAL::3.10000 0.779824 2.37802 2.49998 +DEAL::3.20000 0.772919 2.42731 2.49998 +DEAL::3.30000 0.768758 2.47441 2.49999 +DEAL::3.40000 0.767229 2.51916 2.49998 +DEAL::3.50000 0.768257 2.56138 2.49998 +DEAL::3.60000 0.771800 2.60085 2.49998 +DEAL::3.70000 0.777859 2.63733 2.49998 +DEAL::3.80000 0.786474 2.67052 2.49998 +DEAL::3.90000 0.797733 2.70008 2.49998 +DEAL::4.00000 0.811773 2.72558 2.49998 +DEAL::4.10000 0.828789 2.74657 2.49998 +DEAL::4.20000 0.849040 2.76245 2.49998 +DEAL::4.30000 0.872862 2.77257 2.49998 +DEAL::4.40000 0.900673 2.77611 2.49998 +DEAL::4.50000 0.932987 2.77216 2.49998 +DEAL::4.60000 0.970423 2.75959 2.49998 +DEAL::4.70000 1.01371 2.73716 2.49998 +DEAL::4.80000 1.06366 2.70339 2.49997 +DEAL::4.90000 1.12115 2.65673 2.49997 +DEAL::5.00000 1.18700 2.59554 2.49997 +DEAL::5.10000 1.26179 2.51839 2.49997 +DEAL::5.20000 1.34554 2.42435 2.49997 +DEAL::5.30000 1.43719 2.31362 2.49996 +DEAL::5.40000 1.53404 2.18823 2.49996 +DEAL::5.50000 1.63130 2.05269 2.49996 +DEAL::5.60000 1.72208 1.91414 2.49996 +DEAL::5.70000 1.79854 1.78150 2.49996 +DEAL::5.80000 1.85371 1.66352 2.49995 +DEAL::5.90000 1.88345 1.56671 2.49995 +DEAL::6.00000 1.88731 1.49410 2.49995 +DEAL::6.10000 1.86803 1.44543 2.49995 +DEAL::6.20000 1.83015 1.41827 2.49995 +DEAL::6.30000 1.77867 1.40922 2.49995 +DEAL::6.40000 1.71817 1.41482 2.49996 +DEAL::6.50000 1.65242 1.43201 2.49996 +DEAL::6.60000 1.58434 1.45824 2.49996 +DEAL::6.70000 1.51607 1.49149 2.49996 +DEAL::6.80000 1.44916 1.53017 2.49996 +DEAL::6.90000 1.38465 1.57301 2.49997 +DEAL::7.00000 1.32327 1.61902 2.49997 +DEAL::7.10000 1.26548 1.66740 2.49997 +DEAL::7.20000 1.21156 1.71750 2.49997 +DEAL::7.30000 1.16164 1.76879 2.49997 +DEAL::7.40000 1.11578 1.82081 2.49997 +DEAL::7.50000 1.07396 1.87319 2.49997 +DEAL::7.60000 1.03609 1.92558 2.49997 +DEAL::7.70000 1.00208 1.97772 2.49997 +DEAL::7.80000 0.971784 2.02935 2.49998 +DEAL::7.90000 0.945072 2.08025 2.49998 +DEAL::8.00000 0.921787 2.13022 2.49998 +DEAL::8.10000 0.901781 2.17907 2.49998 +DEAL::8.20000 0.884906 2.22664 2.49998 +DEAL::8.30000 0.871025 2.27275 2.49998 +DEAL::8.40000 0.860013 2.31723 2.49998 +DEAL::8.50000 0.851761 2.35992 2.49998 +DEAL::8.60000 0.846180 2.40062 2.49998 +DEAL::8.70000 0.843199 2.43916 2.49998 +DEAL::8.80000 0.842774 2.47531 2.49998 +DEAL::8.90000 0.844885 2.50883 2.49998 +DEAL::9.00000 0.849543 2.53948 2.49998 +DEAL::9.10000 0.856788 2.56694 2.49998 +DEAL::9.20000 0.866696 2.59088 2.49998 +DEAL::9.30000 0.879382 2.61091 2.49998 +DEAL::9.40000 0.895002 2.62660 2.49998 +DEAL::9.50000 0.913755 2.63743 2.49998 +DEAL::9.60000 0.935893 2.64284 2.49998 +DEAL::9.70000 0.961713 2.64217 2.49998 +DEAL::9.80000 0.991567 2.63469 2.49998 +DEAL::9.90000 1.02585 2.61958 2.49998 +DEAL::10.0000 1.06497 2.59596 2.49997 diff --git a/tests/sundials/arkode_06.with_sundials=4.0.0.output b/tests/sundials/arkode_06.with_sundials=4.0.0.output new file mode 100644 index 0000000000..32dccf75f3 --- /dev/null +++ b/tests/sundials/arkode_06.with_sundials=4.0.0.output @@ -0,0 +1,102 @@ + +DEAL::0.00000 3.90000 1.10000 2.80000 +DEAL::0.100000 3.99995 0.720432 2.49990 +DEAL::0.200000 3.78894 0.660967 2.49991 +DEAL::0.300000 3.52275 0.681501 2.49992 +DEAL::0.400000 3.26217 0.722960 2.49993 +DEAL::0.500000 3.01935 0.771864 2.49994 +DEAL::0.600000 2.79594 0.824670 2.49995 +DEAL::0.700000 2.59116 0.880241 2.49996 +DEAL::0.800000 2.40370 0.938097 2.49996 +DEAL::0.900000 2.23217 0.997959 2.49997 +DEAL::1.00000 2.07527 1.05960 2.49997 +DEAL::1.10000 1.93182 1.12281 2.49997 +DEAL::1.20000 1.80073 1.18737 2.49998 +DEAL::1.30000 1.68106 1.25304 2.49998 +DEAL::1.40000 1.57194 1.31959 2.49998 +DEAL::1.50000 1.47261 1.38677 2.49998 +DEAL::1.60000 1.38237 1.45434 2.49998 +DEAL::1.70000 1.30060 1.52203 2.49999 +DEAL::1.80000 1.22670 1.58963 2.49999 +DEAL::1.90000 1.16015 1.65689 2.49999 +DEAL::2.00000 1.10046 1.72361 2.49999 +DEAL::2.10000 1.04715 1.78959 2.49999 +DEAL::2.20000 0.999783 1.85466 2.49999 +DEAL::2.30000 0.957942 1.91866 2.49999 +DEAL::2.40000 0.921233 1.98145 2.49999 +DEAL::2.50000 0.889283 2.04291 2.49999 +DEAL::2.60000 0.861743 2.10294 2.49999 +DEAL::2.70000 0.838291 2.16142 2.49999 +DEAL::2.80000 0.818628 2.21827 2.49999 +DEAL::2.90000 0.802485 2.27338 2.49998 +DEAL::3.00000 0.789620 2.32667 2.49998 +DEAL::3.10000 0.779824 2.37802 2.49998 +DEAL::3.20000 0.772919 2.42731 2.49998 +DEAL::3.30000 0.768758 2.47441 2.49998 +DEAL::3.40000 0.767229 2.51916 2.49998 +DEAL::3.50000 0.768257 2.56138 2.49997 +DEAL::3.60000 0.771800 2.60085 2.50001 +DEAL::3.70000 0.777859 2.63733 2.50006 +DEAL::3.80000 0.786474 2.67052 2.49990 +DEAL::3.90000 0.797733 2.70008 2.50014 +DEAL::4.00000 0.811773 2.72558 2.49986 +DEAL::4.10000 0.828789 2.74657 2.50000 +DEAL::4.20000 0.849040 2.76245 2.50012 +DEAL::4.30000 0.872862 2.77257 2.50006 +DEAL::4.40000 0.900673 2.77611 2.49984 +DEAL::4.50000 0.932987 2.77216 2.50011 +DEAL::4.60000 0.970423 2.75959 2.50002 +DEAL::4.70000 1.01371 2.73716 2.50005 +DEAL::4.80000 1.06366 2.70339 2.50003 +DEAL::4.90000 1.12115 2.65673 2.50001 +DEAL::5.00000 1.18700 2.59554 2.49996 +DEAL::5.10000 1.26179 2.51839 2.50006 +DEAL::5.20000 1.34554 2.42435 2.49986 +DEAL::5.30000 1.43719 2.31362 2.49999 +DEAL::5.40000 1.53404 2.18823 2.49998 +DEAL::5.50000 1.63130 2.05269 2.49988 +DEAL::5.60000 1.72208 1.91414 2.49987 +DEAL::5.70000 1.79854 1.78150 2.49992 +DEAL::5.80000 1.85371 1.66352 2.49996 +DEAL::5.90000 1.88345 1.56671 2.49997 +DEAL::6.00000 1.88731 1.49410 2.49995 +DEAL::6.10000 1.86803 1.44543 2.49990 +DEAL::6.20000 1.83015 1.41827 2.49998 +DEAL::6.30000 1.77867 1.40922 2.49995 +DEAL::6.40000 1.71817 1.41482 2.50007 +DEAL::6.50000 1.65242 1.43201 2.50007 +DEAL::6.60000 1.58434 1.45824 2.49999 +DEAL::6.70000 1.51607 1.49149 2.49987 +DEAL::6.80000 1.44916 1.53017 2.50000 +DEAL::6.90000 1.38465 1.57301 2.50003 +DEAL::7.00000 1.32327 1.61902 2.49990 +DEAL::7.10000 1.26548 1.66740 2.49995 +DEAL::7.20000 1.21156 1.71750 2.50004 +DEAL::7.30000 1.16164 1.76879 2.49994 +DEAL::7.40000 1.11578 1.82081 2.49993 +DEAL::7.50000 1.07396 1.87319 2.50003 +DEAL::7.60000 1.03609 1.92558 2.49996 +DEAL::7.70000 1.00208 1.97772 2.49993 +DEAL::7.80000 0.971784 2.02935 2.50002 +DEAL::7.90000 0.945072 2.08025 2.49998 +DEAL::8.00000 0.921787 2.13022 2.49994 +DEAL::8.10000 0.901781 2.17907 2.50000 +DEAL::8.20000 0.884906 2.22664 2.49999 +DEAL::8.30000 0.871025 2.27275 2.49996 +DEAL::8.40000 0.860013 2.31723 2.49998 +DEAL::8.50000 0.851761 2.35992 2.49989 +DEAL::8.60000 0.846180 2.40062 2.49997 +DEAL::8.70000 0.843199 2.43916 2.49998 +DEAL::8.80000 0.842774 2.47531 2.49999 +DEAL::8.90000 0.844885 2.50883 2.49996 +DEAL::9.00000 0.849543 2.53948 2.50003 +DEAL::9.10000 0.856788 2.56694 2.50000 +DEAL::9.20000 0.866696 2.59088 2.49989 +DEAL::9.30000 0.879382 2.61091 2.50011 +DEAL::9.40000 0.895002 2.62660 2.50004 +DEAL::9.50000 0.913755 2.63743 2.49995 +DEAL::9.60000 0.935893 2.64284 2.49998 +DEAL::9.70000 0.961713 2.64217 2.50011 +DEAL::9.80000 0.991567 2.63469 2.49985 +DEAL::9.90000 1.02585 2.61958 2.50009 +DEAL::10.0000 1.06497 2.59596 2.49997 diff --git a/tests/sundials/arkode_07.cc b/tests/sundials/arkode_07.cc new file mode 100644 index 0000000000..cbbb08378a --- /dev/null +++ b/tests/sundials/arkode_07.cc @@ -0,0 +1,183 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +//----------------------------------------------------------- + +#include + +#include +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + + +// Test implicit-explicit time stepper. jac_times_vector/setup + +// custom linear solver + custom preconditioner supplied through SUNDIALS + +/** + * This test problem is called "brusselator", and is a typical benchmark for + * ODE solvers. This problem has 3 dependent variables u, v and w, that depend + * on the independent variable t via the IVP system + * + * du/dt = a − (w + 1)u + v u^2 + * dv/dt = w u − v u^2 + * dw/dt = (b − w)/eps -w u + * + * We integrate over the interval 0 ≤ t ≤ 10, with the initial conditions + * + * u(0) = 3.9, v(0) = 1.1, w(0) = 2.8, + * + * and parameters + * + * a = 1.2, b = 2.5, and eps = 10−5 + * + * The implicit part only contains the stiff part of the problem (the part with + * eps in right hand side of the third equation). + */ +int +main(int argc, char **argv) +{ + initlog(); + // restrict output to highest level + deallog.depth_file(1); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + typedef Vector VectorType; + + ParameterHandler prm; + SUNDIALS::ARKode::AdditionalData data; + data.add_parameters(prm); + + if (false) + { + std::ofstream ofile(SOURCE_DIR "/arkode_07.prm"); + prm.print_parameters(ofile, ParameterHandler::ShortText); + ofile.close(); + } + + std::ifstream ifile(SOURCE_DIR "/arkode_07.prm"); + prm.parse_input(ifile); + + SUNDIALS::ARKode ode(data); + + ode.reinit_vector = [&](VectorType &v) { + // Three independent variables + v.reinit(3); + }; + + // Parameters + double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b = 2.5, eps = 1e-5; + // Explicit jacobian. + FullMatrix J(3, 3); + + ode.implicit_function = + [&](double, const VectorType &y, VectorType &ydot) -> int { + ydot[0] = 0; + ydot[1] = 0; + ydot[2] = (b - y[2]) / eps; + return 0; + }; + + + ode.explicit_function = + [&](double, const VectorType &y, VectorType &ydot) -> int { + ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0]; + ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0]; + ydot[2] = -y[2] * y[0]; + return 0; + }; + + + ode.jacobian_times_setup = + [&](realtype t, const VectorType &y, const VectorType &fy) -> int { + J = 0; + J(2, 2) = -1.0 / eps; + return 0; + }; + + ode.jacobian_times_vector = [&](const VectorType &v, + VectorType & Jv, + double t, + const VectorType &y, + const VectorType &fy) -> int { + J.vmult(Jv, v); + return 0; + }; + + ode.solve_linearized_system = + [&](SUNDIALS::SundialsOperator & op, + SUNDIALS::SundialsPreconditioner &prec, + VectorType & x, + const VectorType & b, + double tol) -> int { + ReductionControl control; + SolverCG solver_cg(control); + solver_cg.solve(op, x, b, prec); + return 0; + }; + + ode.jacobian_preconditioner_setup = [&](double t, + const VectorType &y, + const VectorType &fy, + int jok, + int & jcur, + double gamma) -> int { + deallog << "jacobian_preconditioner_setup called\n"; + return 0; + }; + + ode.jacobian_preconditioner_solve = [&](double t, + const VectorType &y, + const VectorType &fy, + const VectorType &r, + VectorType & z, + double gamma, + double delta, + int lr) -> int { + deallog << "jacobian_preconditioner_solve called\n"; + z = r; + return 0; + }; + + + ode.output_step = [&](const double t, + const VectorType & sol, + const unsigned int step_number) -> int { + deallog << t << " " << sol[0] << " " << sol[1] << " " << sol[2] + << std::endl; + return 0; + }; + + // after 5.2.0 a special interpolation mode should be used for stiff problems +#if DEAL_II_SUNDIALS_VERSION_GTE(5, 2, 0) + ode.custom_setup = [&](void *arkode_mem) { + ARKStepSetInterpolantType(arkode_mem, ARK_INTERP_LAGRANGE); + }; +#endif + + Vector y(3); + y[0] = u0; + y[1] = v0; + y[2] = w0; + ode.solve_ode(y); + return 0; +} diff --git a/tests/sundials/arkode_07.prm b/tests/sundials/arkode_07.prm new file mode 100644 index 0000000000..55a91a0b7c --- /dev/null +++ b/tests/sundials/arkode_07.prm @@ -0,0 +1,15 @@ +set Final time = 10. +set Initial time = 0. +set Time interval between each output = 0.1 +subsection Error control + set Absolute error tolerance = 0.000001 + set Relative error tolerance = 0.000010 +end +subsection Running parameters + set Implicit function is linear = true + set Implicit function is time independent = true + set Initial step size = 0.010000 + set Maximum number of nonlinear iterations = 10 + set Maximum order of ARK = 5 + set Minimum step size = 0.000001 +end diff --git a/tests/sundials/arkode_07.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_07.with_sundials.geq.5.4.0.output new file mode 100644 index 0000000000..7c605fb0b8 --- /dev/null +++ b/tests/sundials/arkode_07.with_sundials.geq.5.4.0.output @@ -0,0 +1,102 @@ + +DEAL::0.00000 3.90000 1.10000 2.80000 +DEAL::0.100000 3.99995 0.720432 2.49990 +DEAL::0.200000 3.78894 0.660967 2.49990 +DEAL::0.300000 3.52275 0.681501 2.49991 +DEAL::0.400000 3.26217 0.722960 2.49992 +DEAL::0.500000 3.01935 0.771864 2.49992 +DEAL::0.600000 2.79594 0.824670 2.49993 +DEAL::0.700000 2.59116 0.880241 2.49993 +DEAL::0.800000 2.40370 0.938097 2.49994 +DEAL::0.900000 2.23217 0.997959 2.49994 +DEAL::1.00000 2.07527 1.05960 2.49995 +DEAL::1.10000 1.93182 1.12281 2.49995 +DEAL::1.20000 1.80073 1.18737 2.49995 +DEAL::1.30000 1.68106 1.25304 2.49996 +DEAL::1.40000 1.57194 1.31959 2.49996 +DEAL::1.50000 1.47261 1.38677 2.49996 +DEAL::1.60000 1.38237 1.45434 2.49997 +DEAL::1.70000 1.30060 1.52203 2.49997 +DEAL::1.80000 1.22670 1.58963 2.49997 +DEAL::1.90000 1.16015 1.65689 2.49997 +DEAL::2.00000 1.10046 1.72361 2.49997 +DEAL::2.10000 1.04715 1.78959 2.49997 +DEAL::2.20000 0.999783 1.85466 2.49997 +DEAL::2.30000 0.957942 1.91866 2.49998 +DEAL::2.40000 0.921233 1.98145 2.49998 +DEAL::2.50000 0.889283 2.04291 2.49998 +DEAL::2.60000 0.861743 2.10294 2.49998 +DEAL::2.70000 0.838291 2.16142 2.49998 +DEAL::2.80000 0.818628 2.21827 2.49998 +DEAL::2.90000 0.802485 2.27338 2.49998 +DEAL::3.00000 0.789620 2.32667 2.49998 +DEAL::3.10000 0.779824 2.37802 2.49998 +DEAL::3.20000 0.772919 2.42731 2.49998 +DEAL::3.30000 0.768758 2.47441 2.49999 +DEAL::3.40000 0.767229 2.51916 2.49998 +DEAL::3.50000 0.768257 2.56138 2.49998 +DEAL::3.60000 0.771800 2.60085 2.49998 +DEAL::3.70000 0.777859 2.63733 2.49998 +DEAL::3.80000 0.786474 2.67052 2.49998 +DEAL::3.90000 0.797733 2.70008 2.49998 +DEAL::4.00000 0.811773 2.72558 2.49998 +DEAL::4.10000 0.828789 2.74657 2.49998 +DEAL::4.20000 0.849040 2.76245 2.49998 +DEAL::4.30000 0.872862 2.77257 2.49998 +DEAL::4.40000 0.900673 2.77611 2.49998 +DEAL::4.50000 0.932987 2.77216 2.49998 +DEAL::4.60000 0.970423 2.75959 2.49998 +DEAL::4.70000 1.01371 2.73716 2.49998 +DEAL::4.80000 1.06366 2.70339 2.49997 +DEAL::4.90000 1.12115 2.65673 2.49997 +DEAL::5.00000 1.18700 2.59554 2.49997 +DEAL::5.10000 1.26179 2.51839 2.49997 +DEAL::5.20000 1.34554 2.42435 2.49997 +DEAL::5.30000 1.43719 2.31362 2.49996 +DEAL::5.40000 1.53404 2.18823 2.49996 +DEAL::5.50000 1.63130 2.05269 2.49996 +DEAL::5.60000 1.72208 1.91414 2.49996 +DEAL::5.70000 1.79854 1.78150 2.49996 +DEAL::5.80000 1.85371 1.66352 2.49995 +DEAL::5.90000 1.88345 1.56671 2.49995 +DEAL::6.00000 1.88731 1.49410 2.49995 +DEAL::6.10000 1.86803 1.44543 2.49995 +DEAL::6.20000 1.83015 1.41827 2.49995 +DEAL::6.30000 1.77867 1.40922 2.49995 +DEAL::6.40000 1.71817 1.41482 2.49996 +DEAL::6.50000 1.65242 1.43201 2.49996 +DEAL::6.60000 1.58434 1.45824 2.49996 +DEAL::6.70000 1.51607 1.49149 2.49996 +DEAL::6.80000 1.44916 1.53017 2.49996 +DEAL::6.90000 1.38465 1.57301 2.49997 +DEAL::7.00000 1.32327 1.61902 2.49997 +DEAL::7.10000 1.26548 1.66740 2.49997 +DEAL::7.20000 1.21156 1.71750 2.49997 +DEAL::7.30000 1.16164 1.76879 2.49997 +DEAL::7.40000 1.11578 1.82081 2.49997 +DEAL::7.50000 1.07396 1.87319 2.49997 +DEAL::7.60000 1.03609 1.92558 2.49997 +DEAL::7.70000 1.00208 1.97772 2.49997 +DEAL::7.80000 0.971784 2.02935 2.49998 +DEAL::7.90000 0.945072 2.08025 2.49998 +DEAL::8.00000 0.921787 2.13022 2.49998 +DEAL::8.10000 0.901781 2.17907 2.49998 +DEAL::8.20000 0.884906 2.22664 2.49998 +DEAL::8.30000 0.871025 2.27275 2.49998 +DEAL::8.40000 0.860013 2.31723 2.49998 +DEAL::8.50000 0.851761 2.35992 2.49998 +DEAL::8.60000 0.846180 2.40062 2.49998 +DEAL::8.70000 0.843199 2.43916 2.49998 +DEAL::8.80000 0.842774 2.47531 2.49998 +DEAL::8.90000 0.844885 2.50883 2.49998 +DEAL::9.00000 0.849543 2.53948 2.49998 +DEAL::9.10000 0.856788 2.56694 2.49998 +DEAL::9.20000 0.866696 2.59088 2.49998 +DEAL::9.30000 0.879382 2.61091 2.49998 +DEAL::9.40000 0.895002 2.62660 2.49998 +DEAL::9.50000 0.913755 2.63743 2.49998 +DEAL::9.60000 0.935893 2.64284 2.49998 +DEAL::9.70000 0.961713 2.64217 2.49998 +DEAL::9.80000 0.991567 2.63469 2.49998 +DEAL::9.90000 1.02585 2.61958 2.49998 +DEAL::10.0000 1.06497 2.59596 2.49997 diff --git a/tests/sundials/arkode_07.with_sundials=4.0.0.output b/tests/sundials/arkode_07.with_sundials=4.0.0.output new file mode 100644 index 0000000000..f9684c5bee --- /dev/null +++ b/tests/sundials/arkode_07.with_sundials=4.0.0.output @@ -0,0 +1,102 @@ + +DEAL::0.00000 3.90000 1.10000 2.80000 +DEAL::0.100000 3.99995 0.720432 2.49990 +DEAL::0.200000 3.78894 0.660967 2.49991 +DEAL::0.300000 3.52275 0.681501 2.49992 +DEAL::0.400000 3.26217 0.722960 2.49993 +DEAL::0.500000 3.01935 0.771864 2.49994 +DEAL::0.600000 2.79594 0.824670 2.49995 +DEAL::0.700000 2.59116 0.880241 2.49996 +DEAL::0.800000 2.40370 0.938097 2.49996 +DEAL::0.900000 2.23217 0.997959 2.49997 +DEAL::1.00000 2.07527 1.05960 2.49997 +DEAL::1.10000 1.93182 1.12281 2.49997 +DEAL::1.20000 1.80073 1.18737 2.49998 +DEAL::1.30000 1.68106 1.25304 2.49998 +DEAL::1.40000 1.57194 1.31959 2.49998 +DEAL::1.50000 1.47261 1.38677 2.49998 +DEAL::1.60000 1.38237 1.45434 2.49998 +DEAL::1.70000 1.30060 1.52203 2.49999 +DEAL::1.80000 1.22670 1.58963 2.49999 +DEAL::1.90000 1.16015 1.65689 2.49999 +DEAL::2.00000 1.10046 1.72361 2.49999 +DEAL::2.10000 1.04715 1.78959 2.49999 +DEAL::2.20000 0.999783 1.85466 2.49999 +DEAL::2.30000 0.957942 1.91866 2.49999 +DEAL::2.40000 0.921233 1.98145 2.49999 +DEAL::2.50000 0.889283 2.04291 2.49999 +DEAL::2.60000 0.861743 2.10294 2.49999 +DEAL::2.70000 0.838291 2.16142 2.49999 +DEAL::2.80000 0.818628 2.21827 2.49999 +DEAL::2.90000 0.802485 2.27338 2.49998 +DEAL::3.00000 0.789620 2.32667 2.49998 +DEAL::3.10000 0.779824 2.37802 2.49998 +DEAL::3.20000 0.772919 2.42731 2.49998 +DEAL::3.30000 0.768758 2.47441 2.49998 +DEAL::3.40000 0.767229 2.51916 2.49998 +DEAL::3.50000 0.768257 2.56138 2.49997 +DEAL::3.60000 0.771800 2.60085 2.50001 +DEAL::3.70000 0.777859 2.63733 2.50006 +DEAL::3.80000 0.786474 2.67052 2.49990 +DEAL::3.90000 0.797733 2.70008 2.50014 +DEAL::4.00000 0.811773 2.72558 2.49986 +DEAL::4.10000 0.828789 2.74657 2.50000 +DEAL::4.20000 0.849040 2.76245 2.50012 +DEAL::4.30000 0.872862 2.77257 2.50006 +DEAL::4.40000 0.900673 2.77611 2.49984 +DEAL::4.50000 0.932987 2.77216 2.50011 +DEAL::4.60000 0.970423 2.75959 2.50002 +DEAL::4.70000 1.01371 2.73716 2.50005 +DEAL::4.80000 1.06366 2.70339 2.50003 +DEAL::4.90000 1.12115 2.65673 2.50001 +DEAL::5.00000 1.18700 2.59554 2.49996 +DEAL::5.10000 1.26179 2.51839 2.50006 +DEAL::5.20000 1.34554 2.42435 2.49986 +DEAL::5.30000 1.43719 2.31362 2.49999 +DEAL::5.40000 1.53404 2.18823 2.49998 +DEAL::5.50000 1.63130 2.05269 2.49988 +DEAL::5.60000 1.72208 1.91414 2.49987 +DEAL::5.70000 1.79854 1.78150 2.49992 +DEAL::5.80000 1.85371 1.66352 2.49996 +DEAL::5.90000 1.88345 1.56671 2.49997 +DEAL::6.00000 1.88731 1.49410 2.49995 +DEAL::6.10000 1.86803 1.44543 2.49990 +DEAL::6.20000 1.83015 1.41827 2.49998 +DEAL::6.30000 1.77867 1.40922 2.49995 +DEAL::6.40000 1.71817 1.41482 2.50007 +DEAL::6.50000 1.65242 1.43201 2.50007 +DEAL::6.60000 1.58434 1.45824 2.49999 +DEAL::6.70000 1.51607 1.49149 2.49987 +DEAL::6.80000 1.44916 1.53017 2.50000 +DEAL::6.90000 1.38465 1.57301 2.50003 +DEAL::7.00000 1.32327 1.61902 2.49990 +DEAL::7.10000 1.26548 1.66740 2.49995 +DEAL::7.20000 1.21156 1.71750 2.50004 +DEAL::7.30000 1.16164 1.76879 2.49994 +DEAL::7.40000 1.11578 1.82081 2.49993 +DEAL::7.50000 1.07396 1.87319 2.50003 +DEAL::7.60000 1.03609 1.92558 2.49996 +DEAL::7.70000 1.00208 1.97772 2.49993 +DEAL::7.80000 0.971784 2.02935 2.50002 +DEAL::7.90000 0.945072 2.08025 2.49998 +DEAL::8.00000 0.921787 2.13022 2.49994 +DEAL::8.10000 0.901781 2.17907 2.50000 +DEAL::8.20000 0.884906 2.22664 2.49999 +DEAL::8.30000 0.871025 2.27275 2.49996 +DEAL::8.40000 0.860013 2.31723 2.49998 +DEAL::8.50000 0.851761 2.35992 2.49989 +DEAL::8.60000 0.846180 2.40062 2.49997 +DEAL::8.70000 0.843199 2.43916 2.49998 +DEAL::8.80000 0.842774 2.47531 2.49999 +DEAL::8.90000 0.844885 2.50883 2.49996 +DEAL::9.00000 0.849543 2.53948 2.50003 +DEAL::9.10000 0.856788 2.56694 2.50000 +DEAL::9.20000 0.866696 2.59088 2.49989 +DEAL::9.30000 0.879382 2.61091 2.50011 +DEAL::9.40000 0.895002 2.62660 2.50003 +DEAL::9.50000 0.913755 2.63743 2.49994 +DEAL::9.60000 0.935893 2.64284 2.49997 +DEAL::9.70000 0.961713 2.64217 2.50011 +DEAL::9.80000 0.991567 2.63469 2.49985 +DEAL::9.90000 1.02585 2.61958 2.50009 +DEAL::10.0000 1.06497 2.59596 2.49997 diff --git a/tests/sundials/arkode_08.cc b/tests/sundials/arkode_08.cc new file mode 100644 index 0000000000..97f8a863b2 --- /dev/null +++ b/tests/sundials/arkode_08.cc @@ -0,0 +1,162 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +//----------------------------------------------------------- + +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + + +// Test implicit-explicit time stepper. solve_jacobian + custom linear solver + +// custom mass solver + +/** + * Test problem inspired by linear 1D FE problem with two unknowns u = [u1 u2]: + * + * M*u' = K*u + f + * + * where the mass matrix M = [2/3 1/3; 1/3 2/3] + * and the stiffness matrix K = [1/2 -1/2; 1/2 -1/2] + * and the force vector f = [1,2]. + * The initial values is chosen as u0 = [1,0]. + * + * The implicit part is taken as K*u while f is explicit. + * + * Note that this is not the classical second-order mechanical system and thus + * the term "mass matrix" may be misleading. It is still used here to be in line + * with SUNDIALS nomenclature. + */ +int +main(int argc, char **argv) +{ + initlog(); + + typedef Vector VectorType; + + ParameterHandler prm; + SUNDIALS::ARKode::AdditionalData data; + data.add_parameters(prm); + + if (false) + { + std::ofstream ofile(SOURCE_DIR "/arkode_08.prm"); + prm.print_parameters(ofile, ParameterHandler::ShortText); + ofile.close(); + } + + std::ifstream ifile(SOURCE_DIR "/arkode_08.prm"); + prm.parse_input(ifile); + + SUNDIALS::ARKode ode(data); + + ode.reinit_vector = [&](VectorType &v) { v.reinit(2); }; + + // Explicit jacobian = stiffness matrix + FullMatrix K(2, 2); + K(0, 0) = K(1, 1) = 0.5; + K(1, 0) = K(0, 1) = -0.5; + + // mass matrix + FullMatrix M(2, 2); + M(0, 0) = M(1, 1) = 2.0 / 3; + M(1, 0) = M(0, 1) = 1.0 / 3; + + ode.implicit_function = + [&](double, const VectorType &y, VectorType &ydot) -> int { + K.vmult(ydot, y); + return 0; + }; + + + ode.explicit_function = + [&](double, const VectorType &y, VectorType &ydot) -> int { + ydot[0] = 1; + ydot[1] = 2; + return 0; + }; + + ode.jacobian_times_vector = [&](const VectorType &v, + VectorType & Jv, + double t, + const VectorType &y, + const VectorType &fy) -> int { + K.vmult(Jv, v); + return 0; + }; + + const auto solve_function = + [&](SUNDIALS::SundialsOperator & op, + SUNDIALS::SundialsPreconditioner &prec, + VectorType & x, + const VectorType & b, + double tol) -> int { + ReductionControl control; + SolverCG solver_cg(control); + solver_cg.solve(op, x, b, prec); + return 0; + }; + + ode.solve_linearized_system = solve_function; + + ode.solve_mass = solve_function; + + FullMatrix M_inv(2, 2); + + ode.mass_preconditioner_solve = [&](double t, + const VectorType &r, + VectorType & z, + double gamma, + int lr) -> int { + LogStream::Prefix prefix("mass_preconditioner_solve"); + deallog << "applied" << std::endl; + M_inv.vmult(z, r); + return 0; + }; + + ode.mass_preconditioner_setup = [&](double t) -> int { + LogStream::Prefix prefix("mass_preconditioner_setup"); + deallog << "applied" << std::endl; + M_inv.invert(M); + return 0; + }; + + ode.mass_times_vector = + [&](const double t, const VectorType &v, VectorType &Mv) -> int { + M.vmult(Mv, v); + return 0; + }; + + + ode.output_step = [&](const double t, + const VectorType & sol, + const unsigned int step_number) -> int { + deallog << t << " " << sol[0] << " " << sol[1] << std::endl; + return 0; + }; + + Vector y(2); + y[0] = 1; + y[1] = 0; + ode.solve_ode(y); + return 0; +} diff --git a/tests/sundials/arkode_08.prm b/tests/sundials/arkode_08.prm new file mode 100644 index 0000000000..f7d62b16cd --- /dev/null +++ b/tests/sundials/arkode_08.prm @@ -0,0 +1,15 @@ +set Final time = 1. +set Initial time = 0. +set Time interval between each output = 0.1 +subsection Error control + set Absolute error tolerance = 0.000001 + set Relative error tolerance = 0.000010 +end +subsection Running parameters + set Implicit function is linear = true + set Implicit function is time independent = true + set Initial step size = 0.010000 + set Maximum number of nonlinear iterations = 10 + set Maximum order of ARK = 5 + set Minimum step size = 0.000001 +end diff --git a/tests/sundials/arkode_08.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_08.with_sundials.geq.5.4.0.output new file mode 100644 index 0000000000..1783c03d3c --- /dev/null +++ b/tests/sundials/arkode_08.with_sundials.geq.5.4.0.output @@ -0,0 +1,279 @@ + +DEAL::0.00000 1.00000 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.00869741 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.00551394 +DEAL:cg::Convergence step 1 value 6.13317e-19 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 7.85046e-17 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.00420341 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.0195161 +DEAL:cg::Convergence step 1 value 1.35486e-17 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.00509117 +DEAL:cg::Convergence step 1 value 1.30104e-18 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.0127279 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 7.85046e-17 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.0212132 +DEAL:cg::Convergence step 1 value 5.48568e-18 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.235321 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.149188 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.113729 +DEAL:cg::Convergence step 1 value 1.38778e-17 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.528038 +DEAL:cg::Convergence step 1 value 5.55112e-17 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.137749 +DEAL:cg::Convergence step 1 value 3.10317e-17 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.344373 +DEAL:cg::Convergence step 1 value 2.16778e-16 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.573955 +DEAL:cg::Convergence step 1 value 7.85046e-17 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.100000 1.15000 0.150000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.200000 1.30000 0.300000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.625723 +DEAL:cg::Convergence step 1 value 1.11022e-16 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.396692 +DEAL:cg::Convergence step 1 value 5.55112e-17 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.302408 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 1.40406 +DEAL:cg::Convergence step 1 value 1.11022e-16 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.366277 +DEAL:cg::Convergence step 1 value 7.85046e-17 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.915692 +DEAL:cg::Convergence step 1 value 3.14018e-16 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 1.52615 +DEAL:cg::Convergence step 1 value 4.96507e-16 +DEAL:cg::Starting value 0.707107 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 2.23607 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.300000 1.45000 0.450000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.400000 1.60000 0.600000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.500000 1.75000 0.750000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.600000 1.90000 0.900000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.700000 2.05000 1.05000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.800000 2.20000 1.20000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.900000 2.35000 1.35000 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::1.00000 2.50000 1.50000 diff --git a/tests/sundials/arkode_08.with_sundials=4.0.0.output b/tests/sundials/arkode_08.with_sundials=4.0.0.output new file mode 100644 index 0000000000..1f5cf89bca --- /dev/null +++ b/tests/sundials/arkode_08.with_sundials=4.0.0.output @@ -0,0 +1,89 @@ + +DEAL::0.00000 1.00000 0.00000 +DEAL:mass_preconditioner_setup::applied +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Convergence step 0 value 0.00000 +DEAL:cg::Starting value 0.00869741 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.00551394 +DEAL:cg::Convergence step 1 value 4.33681e-19 +DEAL:cg::Starting value 0.00420341 +DEAL:cg::Convergence step 1 value 1.37142e-18 +DEAL:cg::Starting value 0.0195161 +DEAL:cg::Convergence step 1 value 9.81308e-18 +DEAL:cg::Starting value 0.00509117 +DEAL:cg::Convergence step 1 value 8.67362e-19 +DEAL:cg::Starting value 0.0127279 +DEAL:cg::Convergence step 1 value 3.46945e-18 +DEAL:cg::Starting value 0.0212132 +DEAL:cg::Convergence step 1 value 2.45327e-18 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.0212132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 2.45327e-18 +DEAL:cg::Starting value 2.74317e-18 +DEAL:cg::Convergence step 0 value 2.74317e-18 +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Convergence step 0 value 0.00000 +DEAL:cg::Starting value 0.235321 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.149188 +DEAL:cg::Convergence step 1 value 1.96262e-17 +DEAL:cg::Starting value 0.113729 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 0.528038 +DEAL:cg::Convergence step 1 value 5.55112e-17 +DEAL:cg::Starting value 0.137749 +DEAL:cg::Convergence step 1 value 2.77556e-17 +DEAL:cg::Starting value 0.344373 +DEAL:cg::Convergence step 1 value 8.77708e-17 +DEAL:cg::Starting value 0.573955 +DEAL:cg::Convergence step 1 value 1.75542e-16 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 0.573955 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 1.24127e-16 +DEAL:cg::Starting value 7.36006e-17 +DEAL:cg::Convergence step 0 value 7.36006e-17 +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.100000 1.15000 0.150000 +DEAL::0.200000 1.30000 0.300000 +DEAL:cg::Convergence step 0 value 0.00000 +DEAL:cg::Starting value 0.625723 +DEAL:cg::Convergence step 1 value 1.75542e-16 +DEAL:cg::Starting value 0.396692 +DEAL:cg::Convergence step 1 value 1.24127e-16 +DEAL:cg::Starting value 0.302408 +DEAL:cg::Convergence step 1 value 6.20634e-17 +DEAL:cg::Starting value 1.40406 +DEAL:cg::Convergence step 1 value 3.14018e-16 +DEAL:cg::Starting value 0.366277 +DEAL:cg::Convergence step 1 value 5.55112e-17 +DEAL:cg::Starting value 0.915692 +DEAL:cg::Convergence step 1 value 1.11022e-16 +DEAL:cg::Starting value 1.52615 +DEAL:cg::Convergence step 1 value 2.22045e-16 +DEAL:mass_preconditioner_setup::applied +DEAL:cg::Starting value 1.52615 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:cg::Starting value 1.98900e-16 +DEAL:cg::Convergence step 0 value 1.98900e-16 +DEAL:cg::Starting value 2.12132 +DEAL:cg:mass_preconditioner_solve::applied +DEAL:cg::Convergence step 1 value 0.00000 +DEAL::0.300000 1.45000 0.450000 +DEAL::0.400000 1.60000 0.600000 +DEAL::0.500000 1.75000 0.750000 +DEAL::0.600000 1.90000 0.900000 +DEAL::0.700000 2.05000 1.05000 +DEAL::0.800000 2.20000 1.20000 +DEAL::0.900000 2.35000 1.35000 +DEAL::1.00000 2.50000 1.50000 diff --git a/tests/sundials/harmonic_oscillator_02.output b/tests/sundials/harmonic_oscillator_02.output deleted file mode 100644 index 3571d90f5a..0000000000 --- a/tests/sundials/harmonic_oscillator_02.output +++ /dev/null @@ -1,127 +0,0 @@ -0 0 1 -0.05 0.0499791 0.99875 -0.1 0.0998334 0.995004 -0.15 0.149438 0.988771 -0.2 0.198669 0.980066 -0.25 0.247404 0.968911 -0.3 0.295516 0.955326 -0.35 0.34289 0.939353 -0.4 0.38941 0.92104 -0.45 0.43496 0.900433 -0.5 0.479424 0.877579 -0.55 0.522687 0.852524 -0.6 0.564637 0.825329 -0.65 0.605176 0.796071 -0.7 0.644206 0.764828 -0.75 0.681631 0.73168 -0.8 0.717355 0.696705 -0.85 0.751278 0.659982 -0.9 0.783303 0.621595 -0.95 0.813363 0.581652 -1 0.841396 0.540258 -1.05 0.867339 0.497522 -1.1 0.891131 0.453553 -1.15 0.912711 0.408458 -1.2 0.932015 0.362345 -1.25 0.948983 0.315321 -1.3 0.963554 0.267498 -1.35 0.975694 0.219003 -1.4 0.985391 0.169961 -1.45 0.992634 0.120495 -1.5 0.997413 0.07073 -1.55 0.999716 0.0207894 -1.6 0.999533 -0.0292026 -1.65 0.996854 -0.0791221 -1.7 0.991667 -0.128845 -1.75 0.983976 -0.178243 -1.8 0.973816 -0.227193 -1.85 0.961228 -0.275575 -1.9 0.946248 -0.323273 -1.95 0.928918 -0.370167 -2 0.909275 -0.41614 -2.05 0.88736 -0.461072 -2.1 0.86321 -0.504847 -2.15 0.836886 -0.547347 -2.2 0.808467 -0.588477 -2.25 0.778033 -0.62814 -2.3 0.745665 -0.666242 -2.35 0.711444 -0.702689 -2.4 0.675451 -0.737384 -2.45 0.637765 -0.770233 -2.5 0.598471 -0.801142 -2.55 0.557673 -0.830034 -2.6 0.515481 -0.856851 -2.65 0.472005 -0.881535 -2.7 0.427357 -0.904029 -2.75 0.381646 -0.924275 -2.8 0.334983 -0.942215 -2.85 0.287479 -0.957791 -2.9 0.239248 -0.970952 -2.95 0.190419 -0.981674 -3 0.141115 -0.989946 -3.05 0.091459 -0.995756 -3.1 0.0415761 -0.999091 -3.15 -0.00841 -0.999941 -3.2 -0.0583754 -0.998292 -3.25 -0.108196 -0.994133 -3.3 -0.157743 -0.987468 -3.35 -0.206893 -0.978328 -3.4 -0.255528 -0.96675 -3.45 -0.303527 -0.952769 -3.5 -0.350773 -0.93642 -3.55 -0.397144 -0.917738 -3.6 -0.442522 -0.89676 -3.65 -0.486788 -0.873522 -3.7 -0.529827 -0.848088 -3.75 -0.571542 -0.820533 -3.8 -0.611833 -0.790935 -3.85 -0.650604 -0.759371 -3.9 -0.687755 -0.725916 -3.95 -0.723188 -0.690649 -4 -0.756806 -0.653646 -4.05 -0.788517 -0.614997 -4.1 -0.818254 -0.574809 -4.15 -0.84595 -0.533188 -4.2 -0.871543 -0.490241 -4.25 -0.894968 -0.446074 -4.3 -0.916161 -0.400795 -4.35 -0.935058 -0.35451 -4.4 -0.951598 -0.307332 -4.45 -0.965751 -0.259384 -4.5 -0.977494 -0.210788 -4.55 -0.986805 -0.161669 -4.6 -0.993662 -0.112147 -4.65 -0.998044 -0.0623455 -4.7 -0.999928 -0.0123877 -4.75 -0.999294 0.0376024 -4.8 -0.996148 0.0874965 -4.85 -0.990514 0.137172 -4.9 -0.982413 0.186506 -4.95 -0.97187 0.235376 -5 -0.958908 0.28366 -5.05 -0.94355 0.331236 -5.1 -0.925819 0.37798 -5.15 -0.905761 0.423774 -5.2 -0.883439 0.468508 -5.25 -0.858915 0.512075 -5.3 -0.832253 0.554367 -5.35 -0.803517 0.595277 -5.4 -0.772769 0.634698 -5.45 -0.740077 0.672522 -5.5 -0.705531 0.708661 -5.55 -0.669225 0.743031 -5.6 -0.631253 0.775551 -5.65 -0.591709 0.80614 -5.7 -0.550686 0.834716 -5.75 -0.508281 0.861197 -5.8 -0.464599 0.885514 -5.85 -0.419756 0.907617 -5.9 -0.373868 0.927459 -5.95 -0.327048 0.944992 -6 -0.279413 0.960169 -6.05 -0.231078 0.972942 -6.1 -0.182163 0.983273 -6.15 -0.132791 0.991144 -6.2 -0.0830885 0.996543 -6.25 -0.0331784 0.999455 -6.28 -0.00318445 1 diff --git a/tests/sundials/harmonic_oscillator_03.output b/tests/sundials/harmonic_oscillator_03.output deleted file mode 100644 index 2b90296f5a..0000000000 --- a/tests/sundials/harmonic_oscillator_03.output +++ /dev/null @@ -1,13 +0,0 @@ -0 0 1 -0.5 0.4793775 0.8775139 -1 0.8412554 0.540225 -1.5 0.9970436 0.07063589 -2 0.9084636 -0.4157026 -2.5 0.5984372 -0.8010479 -3 0.1409945 -0.9894479 -3.5 -0.350428 -0.9357101 -4 -0.7566852 -0.6535909 -4.5 -0.9772954 -0.21075 -5 -0.9583971 0.2834953 -5.5 -0.7052272 0.7082618 -6 -0.279415 0.9600362 diff --git a/tests/sundials/harmonic_oscillator_04.output b/tests/sundials/harmonic_oscillator_04.output deleted file mode 100644 index ae732b9f8c..0000000000 --- a/tests/sundials/harmonic_oscillator_04.output +++ /dev/null @@ -1,11 +0,0 @@ -0 3.9 1.1 2.8 -1 2.075274578 1.059601971 2.499948149 -2 1.100460018 1.723610111 2.499972517 -3 0.7896202372 2.326668043 2.499980207 -4 0.8117726432 2.725583965 2.499979833 -5 1.187002916 2.595541878 2.499970301 -6 1.887311324 1.494098057 2.499952687 -7 1.323271341 1.619016372 2.499966912 -8 0.9217874641 2.130219386 2.499976935 -9 0.849542577 2.539475069 2.499978894 -10 1.064969831 2.595955794 2.499973382 diff --git a/tests/sundials/harmonic_oscillator_05.output b/tests/sundials/harmonic_oscillator_05.output deleted file mode 100644 index 04f8401db1..0000000000 --- a/tests/sundials/harmonic_oscillator_05.output +++ /dev/null @@ -1,101 +0,0 @@ -0 3.9 1.1 2.8 -0.1 3.99995 0.720431 2.49989 -0.2 3.78894 0.660967 2.49996 -0.3 3.52275 0.681501 2.49984 -0.4 3.26217 0.72296 2.49995 -0.5 3.01935 0.771864 2.49995 -0.6 2.79594 0.82467 2.49988 -0.7 2.59116 0.880241 2.49998 -0.8 2.4037 0.938097 2.49993 -0.9 2.23217 0.997959 2.49992 -1 2.07527 1.0596 2.49999 -1.1 1.93182 1.12281 2.49992 -1.2 1.80073 1.18737 2.49995 -1.3 1.68106 1.25304 2.49999 -1.4 1.57194 1.31959 2.49993 -1.5 1.47261 1.38677 2.49997 -1.6 1.38237 1.45434 2.49999 -1.7 1.3006 1.52203 2.49994 -1.8 1.2267 1.58963 2.49998 -1.9 1.16015 1.65689 2.49998 -2 1.10046 1.72361 2.49996 -2.1 1.04715 1.78959 2.49998 -2.2 0.999783 1.85466 2.49997 -2.3 0.957942 1.91866 2.49997 -2.4 0.921233 1.98145 2.49999 -2.5 0.889283 2.04291 2.49997 -2.6 0.861743 2.10294 2.49998 -2.7 0.838291 2.16142 2.49999 -2.8 0.818628 2.21827 2.49997 -2.9 0.802485 2.27338 2.49998 -3 0.78962 2.32667 2.49998 -3.1 0.779824 2.37802 2.49998 -3.2 0.772919 2.42731 2.49998 -3.3 0.768758 2.47441 2.49999 -3.4 0.767229 2.51916 2.49998 -3.5 0.768257 2.56138 2.5 -3.6 0.7718 2.60085 2.49995 -3.7 0.777859 2.63733 2.50005 -3.8 0.786474 2.67052 2.49999 -3.9 0.797733 2.70008 2.49988 -4 0.811773 2.72558 2.5001 -4.1 0.828789 2.74657 2.5001 -4.2 0.84904 2.76245 2.50009 -4.3 0.872862 2.77257 2.50012 -4.4 0.900673 2.77611 2.49986 -4.5 0.932987 2.77216 2.50009 -4.6 0.970423 2.75959 2.49999 -4.7 1.01371 2.73716 2.50002 -4.8 1.06366 2.70339 2.50006 -4.9 1.12115 2.65673 2.49999 -5 1.187 2.59554 2.49993 -5.1 1.26179 2.51839 2.50005 -5.2 1.34554 2.42435 2.49986 -5.3 1.43719 2.31362 2.50002 -5.4 1.53404 2.18823 2.50001 -5.5 1.6313 2.05269 2.4999 -5.6 1.72208 1.91414 2.49987 -5.7 1.79854 1.7815 2.49991 -5.8 1.85371 1.66352 2.49995 -5.9 1.88345 1.56671 2.49996 -6 1.88731 1.4941 2.49995 -6.1 1.86803 1.44543 2.49992 -6.2 1.83015 1.41827 2.49989 -6.3 1.77867 1.40922 2.50001 -6.4 1.71817 1.41482 2.50005 -6.5 1.65242 1.43201 2.50005 -6.6 1.58434 1.45824 2.50005 -6.7 1.51607 1.49149 2.49991 -6.8 1.44916 1.53017 2.4999 -6.9 1.38465 1.57301 2.50006 -7 1.32327 1.61902 2.49995 -7.1 1.26548 1.6674 2.49989 -7.2 1.21156 1.7175 2.50005 -7.3 1.16164 1.76879 2.49998 -7.4 1.11578 1.82081 2.49991 -7.5 1.07396 1.87319 2.50002 -7.6 1.03609 1.92558 2.5 -7.7 1.00208 1.97772 2.49993 -7.8 0.971784 2.02935 2.49999 -7.9 0.945072 2.08025 2.5 -8 0.921787 2.13022 2.49995 -8.1 0.901781 2.17907 2.49997 -8.2 0.884906 2.22664 2.5 -8.3 0.871025 2.27275 2.49997 -8.4 0.860013 2.31723 2.49997 -8.5 0.851761 2.35992 2.49999 -8.6 0.84618 2.40062 2.49998 -8.7 0.843199 2.43916 2.49998 -8.8 0.842774 2.47531 2.49999 -8.9 0.844885 2.50883 2.49997 -9 0.849543 2.53948 2.50004 -9.1 0.856788 2.56694 2.49994 -9.2 0.866696 2.59088 2.49999 -9.3 0.879382 2.61091 2.50009 -9.4 0.895002 2.6266 2.49992 -9.5 0.913755 2.63743 2.49985 -9.6 0.935893 2.64284 2.49987 -9.7 0.961713 2.64217 2.50002 -9.8 0.991567 2.63469 2.49996 -9.9 1.02585 2.61958 2.50009 -10 1.06497 2.59596 2.49997 diff --git a/tests/sundials/harmonic_oscillator_06.output b/tests/sundials/harmonic_oscillator_06.output deleted file mode 100644 index 04f8401db1..0000000000 --- a/tests/sundials/harmonic_oscillator_06.output +++ /dev/null @@ -1,101 +0,0 @@ -0 3.9 1.1 2.8 -0.1 3.99995 0.720431 2.49989 -0.2 3.78894 0.660967 2.49996 -0.3 3.52275 0.681501 2.49984 -0.4 3.26217 0.72296 2.49995 -0.5 3.01935 0.771864 2.49995 -0.6 2.79594 0.82467 2.49988 -0.7 2.59116 0.880241 2.49998 -0.8 2.4037 0.938097 2.49993 -0.9 2.23217 0.997959 2.49992 -1 2.07527 1.0596 2.49999 -1.1 1.93182 1.12281 2.49992 -1.2 1.80073 1.18737 2.49995 -1.3 1.68106 1.25304 2.49999 -1.4 1.57194 1.31959 2.49993 -1.5 1.47261 1.38677 2.49997 -1.6 1.38237 1.45434 2.49999 -1.7 1.3006 1.52203 2.49994 -1.8 1.2267 1.58963 2.49998 -1.9 1.16015 1.65689 2.49998 -2 1.10046 1.72361 2.49996 -2.1 1.04715 1.78959 2.49998 -2.2 0.999783 1.85466 2.49997 -2.3 0.957942 1.91866 2.49997 -2.4 0.921233 1.98145 2.49999 -2.5 0.889283 2.04291 2.49997 -2.6 0.861743 2.10294 2.49998 -2.7 0.838291 2.16142 2.49999 -2.8 0.818628 2.21827 2.49997 -2.9 0.802485 2.27338 2.49998 -3 0.78962 2.32667 2.49998 -3.1 0.779824 2.37802 2.49998 -3.2 0.772919 2.42731 2.49998 -3.3 0.768758 2.47441 2.49999 -3.4 0.767229 2.51916 2.49998 -3.5 0.768257 2.56138 2.5 -3.6 0.7718 2.60085 2.49995 -3.7 0.777859 2.63733 2.50005 -3.8 0.786474 2.67052 2.49999 -3.9 0.797733 2.70008 2.49988 -4 0.811773 2.72558 2.5001 -4.1 0.828789 2.74657 2.5001 -4.2 0.84904 2.76245 2.50009 -4.3 0.872862 2.77257 2.50012 -4.4 0.900673 2.77611 2.49986 -4.5 0.932987 2.77216 2.50009 -4.6 0.970423 2.75959 2.49999 -4.7 1.01371 2.73716 2.50002 -4.8 1.06366 2.70339 2.50006 -4.9 1.12115 2.65673 2.49999 -5 1.187 2.59554 2.49993 -5.1 1.26179 2.51839 2.50005 -5.2 1.34554 2.42435 2.49986 -5.3 1.43719 2.31362 2.50002 -5.4 1.53404 2.18823 2.50001 -5.5 1.6313 2.05269 2.4999 -5.6 1.72208 1.91414 2.49987 -5.7 1.79854 1.7815 2.49991 -5.8 1.85371 1.66352 2.49995 -5.9 1.88345 1.56671 2.49996 -6 1.88731 1.4941 2.49995 -6.1 1.86803 1.44543 2.49992 -6.2 1.83015 1.41827 2.49989 -6.3 1.77867 1.40922 2.50001 -6.4 1.71817 1.41482 2.50005 -6.5 1.65242 1.43201 2.50005 -6.6 1.58434 1.45824 2.50005 -6.7 1.51607 1.49149 2.49991 -6.8 1.44916 1.53017 2.4999 -6.9 1.38465 1.57301 2.50006 -7 1.32327 1.61902 2.49995 -7.1 1.26548 1.6674 2.49989 -7.2 1.21156 1.7175 2.50005 -7.3 1.16164 1.76879 2.49998 -7.4 1.11578 1.82081 2.49991 -7.5 1.07396 1.87319 2.50002 -7.6 1.03609 1.92558 2.5 -7.7 1.00208 1.97772 2.49993 -7.8 0.971784 2.02935 2.49999 -7.9 0.945072 2.08025 2.5 -8 0.921787 2.13022 2.49995 -8.1 0.901781 2.17907 2.49997 -8.2 0.884906 2.22664 2.5 -8.3 0.871025 2.27275 2.49997 -8.4 0.860013 2.31723 2.49997 -8.5 0.851761 2.35992 2.49999 -8.6 0.84618 2.40062 2.49998 -8.7 0.843199 2.43916 2.49998 -8.8 0.842774 2.47531 2.49999 -8.9 0.844885 2.50883 2.49997 -9 0.849543 2.53948 2.50004 -9.1 0.856788 2.56694 2.49994 -9.2 0.866696 2.59088 2.49999 -9.3 0.879382 2.61091 2.50009 -9.4 0.895002 2.6266 2.49992 -9.5 0.913755 2.63743 2.49985 -9.6 0.935893 2.64284 2.49987 -9.7 0.961713 2.64217 2.50002 -9.8 0.991567 2.63469 2.49996 -9.9 1.02585 2.61958 2.50009 -10 1.06497 2.59596 2.49997 diff --git a/tests/sundials/harmonic_oscillator_01.cc b/tests/sundials/ida_01.cc similarity index 96% rename from tests/sundials/harmonic_oscillator_01.cc rename to tests/sundials/ida_01.cc index 172d117530..6eaacbd6a7 100644 --- a/tests/sundials/harmonic_oscillator_01.cc +++ b/tests/sundials/ida_01.cc @@ -145,11 +145,11 @@ main(int argc, char **argv) ParameterHandler prm; data.add_parameters(prm); - // std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_01.prm"); + // std::ofstream ofile(SOURCE_DIR "/ida_01.prm"); // prm.print_parameters(ofile, ParameterHandler::ShortText); // ofile.close(); - std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_01.prm"); + std::ifstream ifile(SOURCE_DIR "/ida_01.prm"); prm.parse_input(ifile); diff --git a/tests/sundials/harmonic_oscillator_01.prm b/tests/sundials/ida_01.prm similarity index 100% rename from tests/sundials/harmonic_oscillator_01.prm rename to tests/sundials/ida_01.prm diff --git a/tests/sundials/harmonic_oscillator_01.output b/tests/sundials/ida_01.with_sundials.le.4.0.0.output similarity index 100% rename from tests/sundials/harmonic_oscillator_01.output rename to tests/sundials/ida_01.with_sundials.le.4.0.0.output diff --git a/tests/sundials/kinsol_01.output b/tests/sundials/kinsol_01.with_sundials.le.4.0.0.output similarity index 100% rename from tests/sundials/kinsol_01.output rename to tests/sundials/kinsol_01.with_sundials.le.4.0.0.output diff --git a/tests/sundials/serial_sundials.cc b/tests/sundials/serial_sundials.cc index 180548d36f..8a421ba1c5 100644 --- a/tests/sundials/serial_sundials.cc +++ b/tests/sundials/serial_sundials.cc @@ -29,7 +29,7 @@ int main() { initlog(); - +#if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) SUNDIALS::IDA> ida_solver; (void)ida_solver; deallog << "IDA OK" << std::endl; @@ -37,6 +37,7 @@ main() SUNDIALS::KINSOL> kinsol_solver; (void)kinsol_solver; deallog << "KINSOL OK" << std::endl; +#endif SUNDIALS::ARKode> arkode_solver; (void)arkode_solver; diff --git a/tests/sundials/serial_sundials.with_sundials.geq.4.0.0.output b/tests/sundials/serial_sundials.with_sundials.geq.4.0.0.output new file mode 100644 index 0000000000..048c4e5268 --- /dev/null +++ b/tests/sundials/serial_sundials.with_sundials.geq.4.0.0.output @@ -0,0 +1,2 @@ + +DEAL::ARKODE OK diff --git a/tests/sundials/serial_sundials.output b/tests/sundials/serial_sundials.with_sundials.le.4.0.0.output similarity index 100% rename from tests/sundials/serial_sundials.output rename to tests/sundials/serial_sundials.with_sundials.le.4.0.0.output -- 2.39.5