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()
--- /dev/null
+Improved: Update SUNDIALS ARKODE interface to support versions > 4.0.0 and add preconditioner support
+<br>
+(Sebastian Proell, 2020/12/11)
# include <deal.II/lac/vector_memory.h>
# include <arkode/arkode.h>
-# include <arkode/arkode_impl.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
+# include <arkode/arkode_impl.h>
+# endif
# include <nvector/nvector_serial.h>
# ifdef DEAL_II_WITH_MPI
# include <nvector/nvector_parallel.h>
# endif
# include <boost/signals2.hpp>
+# include <sundials/sundials_linearsolver.h>
# include <sundials/sundials_math.h>
# include <sundials/sundials_types.h>
DEAL_II_NAMESPACE_OPEN
+// Forward declarations
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+# ifndef DOXYGEN
+namespace SUNDIALS
+{
+ // Forward declaration
+ template <typename VectorType>
+ struct SundialsOperator;
+
+ template <typename VectorType>
+ struct SundialsPreconditioner;
+
+ template <typename VectorType>
+ class SundialsLinearSolverWrapper;
+} // namespace SUNDIALS
+# endif
+# endif
+
// Shorthand notation for ARKODE error codes.
# define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code))
*/
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 <typename VectorType>
+ using LinearSolveFunction =
+ std::function<int(SundialsOperator<VectorType> & op,
+ SundialsPreconditioner<VectorType> &prec,
+ VectorType & x,
+ const VectorType & b,
+ double tol)>;
+# endif
+
/**
* Interface to SUNDIALS additive Runge-Kutta methods (ARKode).
*
* $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
* $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[
* @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:
*
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)
, 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)
{}
/**
* 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
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");
* 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;
};
/**
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.
*
* 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)>
*
* 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 &res)>
implicit_function;
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
/**
* A function object that users may supply and that is intended to
* prepare the linear solver for subsequent calls to
* - >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<int(const int convfail,
const double t,
* @param[in] t the current time
* @param[in] gamma the current factor to use in the jacobian computation
* @param[in] ycur is the current $y$ vector for the current ARKode
- * internal step
+ * internal step
* @param[in] fcur is the current value of the implicit right-hand side at
- * ycur, $f_I (t_n, ypred)$.
+ * ycur, $f_I (t_n, ypred)$.
*
*
* This function should return:
* - >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<int(const double t,
const double gamma,
/**
- * A function object that users may supply and that is intended to setup
+ * 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).
* - >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<int(const double t)> setup_mass;
* - >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<int(const VectorType &rhs, VectorType &dst)>
solve_mass_system;
+# else
+
+ /**
+ * A function object that users may supply and that is intended to compute
+ * the product of the mass matrix with a given vector `v`. This function
+ * 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<int(double t, const VectorType &v, VectorType &Mv)>
+ 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<int(const double t)> 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<int(const VectorType &v,
+ VectorType & Jv,
+ double t,
+ const VectorType &y,
+ const VectorType &fy)>
+ 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<int(realtype t, const VectorType &y, const VectorType &fy)>
+ 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<VectorType> 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<VectorType> 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<int(double t,
+ const VectorType &y,
+ const VectorType &fy,
+ const VectorType &r,
+ VectorType & z,
+ double gamma,
+ double tol,
+ int lr)>
+ 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<int(double t,
+ const VectorType &y,
+ const VectorType &fy,
+ int jok,
+ int & jcur,
+ double gamma)>
+ 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<int(double t)> mass_preconditioner_setup;
+# endif
/**
* A function object that users may supply and that is intended to
* 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<void(const double t,
const VectorType & sol,
std::function<VectorType &()> 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<void(void *arkode_mem)> custom_setup;
private:
/**
<< "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
*/
GrowingVectorMemory<VectorType> mem;
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+ std::unique_ptr<SundialsLinearSolverWrapper<VectorType>> linear_solver;
+ std::unique_ptr<SundialsLinearSolverWrapper<VectorType>> mass_solver;
+# endif
+
# ifdef DEAL_II_WITH_PETSC
# ifdef PETSC_USE_COMPLEX
static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
# endif // DEAL_II_WITH_PETSC
};
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+
+ /**
+ * A linear operator that wraps SUNDIALS functionality.
+ */
+ template <typename VectorType>
+ 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<VectorType> &solver,
+ void * A_data,
+ ATimesFn a_times_fn);
+
+ private:
+ /**
+ * Reference to the ARKode object that uses this SundialsOperator.
+ */
+ ARKode<VectorType> &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 <typename VectorType>
+ 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<VectorType> &solver,
+ void * P_data,
+ PSolveFn p_solve_fn,
+ double tol);
+
+ private:
+ /**
+ * Reference to the ARKode object that uses this SundialsPreconditioner.
+ */
+ ARKode<VectorType> &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
#include <deal.II/base/mpi.h>
#ifdef DEAL_II_WITH_SUNDIALS
-
-# include <deal.II/base/conditional_ostream.h>
-# include <deal.II/base/exceptions.h>
-# include <deal.II/base/logstream.h>
-# include <deal.II/base/parameter_handler.h>
-# ifdef DEAL_II_WITH_PETSC
-# include <deal.II/lac/petsc_block_vector.h>
-# include <deal.II/lac/petsc_vector.h>
-# endif
-# include <deal.II/lac/vector.h>
-# include <deal.II/lac/vector_memory.h>
-
-# ifdef DEAL_II_SUNDIALS_WITH_IDAS
-# include <idas/idas.h>
-# else
-# include <ida/ida.h>
-# endif
-
-# include <sundials/sundials_config.h>
-# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
-# include <ida/ida_spbcgs.h>
-# include <ida/ida_spgmr.h>
-# include <ida/ida_sptfqmr.h>
-# endif
-# include <boost/signals2.hpp>
-
-# include <nvector/nvector_serial.h>
-# include <sundials/sundials_math.h>
-# include <sundials/sundials_types.h>
-
-# include <memory>
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
+
+# include <deal.II/base/conditional_ostream.h>
+# include <deal.II/base/exceptions.h>
+# include <deal.II/base/logstream.h>
+# include <deal.II/base/parameter_handler.h>
+# ifdef DEAL_II_WITH_PETSC
+# include <deal.II/lac/petsc_block_vector.h>
+# include <deal.II/lac/petsc_vector.h>
+# endif
+# include <deal.II/lac/vector.h>
+# include <deal.II/lac/vector_memory.h>
+
+# ifdef DEAL_II_SUNDIALS_WITH_IDAS
+# include <idas/idas.h>
+# else
+# include <ida/ida.h>
+# endif
+
+# include <sundials/sundials_config.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
+# include <ida/ida_spbcgs.h>
+# include <ida/ida_spgmr.h>
+# include <ida/ida_sptfqmr.h>
+# endif
+# include <boost/signals2.hpp>
+
+# include <nvector/nvector_serial.h>
+# include <sundials/sundials_math.h>
+# include <sundials/sundials_types.h>
+
+# include <memory>
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
{
*/
GrowingVectorMemory<VectorType> mem;
-# ifdef DEAL_II_WITH_PETSC
-# ifdef PETSC_USE_COMPLEX
+# ifdef DEAL_II_WITH_PETSC
+# ifdef PETSC_USE_COMPLEX
static_assert(!std::is_same<VectorType, PETScWrappers::MPI::Vector>::value,
"Sundials does not support complex scalar types, "
"but PETSc is configured to use a complex scalar type!");
!std::is_same<VectorType, PETScWrappers::MPI::BlockVector>::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
#define dealii_sundials_kinsol_h
#include <deal.II/base/config.h>
+
#ifdef DEAL_II_WITH_SUNDIALS
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
+
-# include <deal.II/base/conditional_ostream.h>
-# include <deal.II/base/exceptions.h>
-# include <deal.II/base/logstream.h>
-# include <deal.II/base/mpi.h>
-# include <deal.II/base/parameter_handler.h>
+# include <deal.II/base/conditional_ostream.h>
+# include <deal.II/base/exceptions.h>
+# include <deal.II/base/logstream.h>
+# include <deal.II/base/mpi.h>
+# include <deal.II/base/parameter_handler.h>
-# include <deal.II/lac/vector.h>
-# include <deal.II/lac/vector_memory.h>
+# include <deal.II/lac/vector.h>
+# include <deal.II/lac/vector_memory.h>
-# include <boost/signals2.hpp>
+# include <boost/signals2.hpp>
-# include <kinsol/kinsol.h>
-# include <kinsol/kinsol_impl.h>
-# include <nvector/nvector_serial.h>
-# include <sundials/sundials_math.h>
-# include <sundials/sundials_types.h>
+# include <kinsol/kinsol.h>
+# include <kinsol/kinsol_impl.h>
+# include <nvector/nvector_serial.h>
+# include <sundials/sundials_math.h>
+# include <sundials/sundials_types.h>
-# include <memory>
+# include <memory>
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
{
DEAL_II_NAMESPACE_CLOSE
+# endif
#endif
-
#endif
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+
+# include <sundials/sundials_linearsolver.h>
+
+DEAL_II_NAMESPACE_OPEN
+namespace SUNDIALS
+{
+ namespace internal
+ {
+ /**
+ * Create a new SUNLinearSolver structure without any content and
+ * operations set to `nullptr`.
+ */
+ 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
# include <deal.II/sundials/copy.h>
-# include <arkode/arkode_impl.h>
-# include <sundials/sundials_config.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
+# include <arkode/arkode_impl.h>
+# include <sundials/sundials_config.h>
+# else
+# include <arkode/arkode_arkstep.h>
+# include <sunlinsol/sunlinsol_spgmr.h>
+# include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
+# include <deal.II/sundials/sunlinsol_newempty.h>
+# endif
+# endif
-# include <iomanip>
# include <iostream>
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
// Make sure we know how to call sundials own ARKode() function
const auto &SundialsARKode = ARKode;
+# endif
DEAL_II_NAMESPACE_OPEN
N_Vector yp,
void * user_data)
{
+ Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
GrowingVectorMemory<VectorType> mem;
N_Vector yp,
void * user_data)
{
+ Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
GrowingVectorMemory<VectorType> mem;
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
template <typename VectorType>
int
t_arkode_setup_jacobian(ARKodeMem arkode_mem,
N_Vector,
N_Vector)
{
+ Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
GrowingVectorMemory<VectorType> mem;
*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;
}
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<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
GrowingVectorMemory<VectorType> mem;
int
t_arkode_setup_mass(ARKodeMem arkode_mem, N_Vector, N_Vector, N_Vector)
{
+ Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
int err = solver.setup_mass(arkode_mem->ark_tn);
template <typename VectorType>
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<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
GrowingVectorMemory<VectorType> mem;
return err;
}
+# else
+
+ template <typename VectorType>
+ 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<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(user_data);
+ GrowingVectorMemory<VectorType> mem;
+
+ typename VectorMemory<VectorType>::Pointer src_v(mem);
+ solver.reinit_vector(*src_v);
+ typename VectorMemory<VectorType>::Pointer dst_Jv(mem);
+ solver.reinit_vector(*dst_Jv);
+ typename VectorMemory<VectorType>::Pointer src_y(mem);
+ solver.reinit_vector(*src_y);
+ typename VectorMemory<VectorType>::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 <typename VectorType>
+ int
+ t_arkode_jac_times_setup_function(realtype t,
+ N_Vector y,
+ N_Vector fy,
+ void * user_data)
+ {
+ Assert(user_data != nullptr, ExcInternalError());
+ ARKode<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(user_data);
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer src_y(mem);
+ solver.reinit_vector(*src_y);
+ typename VectorMemory<VectorType>::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 <typename VectorType>
+ 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<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(user_data);
+
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer src_y(mem);
+ solver.reinit_vector(*src_y);
+ typename VectorMemory<VectorType>::Pointer src_fy(mem);
+ solver.reinit_vector(*src_fy);
+ typename VectorMemory<VectorType>::Pointer src_r(mem);
+ solver.reinit_vector(*src_r);
+ typename VectorMemory<VectorType>::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 <typename VectorType>
+ 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<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(user_data);
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer src_y(mem);
+ solver.reinit_vector(*src_y);
+ typename VectorMemory<VectorType>::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 <typename VectorType>
+ int
+ t_arkode_mass_times_vec_function(N_Vector v,
+ N_Vector Mv,
+ realtype t,
+ void * mtimes_data)
+ {
+ Assert(mtimes_data != nullptr, ExcInternalError());
+ ARKode<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(mtimes_data);
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer src_v(mem);
+ solver.reinit_vector(*src_v);
+ typename VectorMemory<VectorType>::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 <typename VectorType>
+ int
+ t_arkode_mass_times_setup_function(realtype t, void *mtimes_data)
+ {
+ Assert(mtimes_data != nullptr, ExcInternalError());
+ ARKode<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(mtimes_data);
+
+ return solver.mass_times_setup(t);
+ }
+
+
+
+ template <typename VectorType>
+ 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<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(user_data);
+
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer src_r(mem);
+ solver.reinit_vector(*src_r);
+ typename VectorMemory<VectorType>::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 <typename VectorType>
+ int
+ t_arkode_mass_prec_setup_function(realtype t, void *user_data)
+ {
+ Assert(user_data != nullptr, ExcInternalError());
+ ARKode<VectorType> &solver =
+ *static_cast<ARKode<VectorType> *>(user_data);
+
+ return solver.mass_preconditioner_setup(t);
+ }
+
+
+
+ /**
+ * storage for internal content of the linear solver wrapper
+ */
+ template <typename VectorType>
+ struct LinearSolverContent
+ {
+ ATimesFn a_times_fn;
+ PSetupFn preconditioner_setup;
+ PSolveFn preconditioner_solve;
+
+ LinearSolveFunction<VectorType> lsolve;
+
+ //! solver reference to forward all solver related calls to user-specified
+ //! functions
+ ARKode<VectorType> *solver;
+ void * P_data;
+ void * A_data;
+ };
+
+
+
+ /**
+ * Access our LinearSolverContent from the generic content of the
+ * SUNLinearSolver @p ls.
+ */
+ template <typename VectorType>
+ LinearSolverContent<VectorType> *
+ access_content(SUNLinearSolver ls)
+ {
+ Assert(ls->content != nullptr, ExcInternalError());
+ return static_cast<LinearSolverContent<VectorType> *>(ls->content);
+ }
+
+
+
+ SUNLinearSolver_Type arkode_linsol_get_type(SUNLinearSolver)
+ {
+ return SUNLINEARSOLVER_ITERATIVE;
+ }
+
+
+
+ template <typename VectorType>
+ int
+ arkode_linsol_solve(SUNLinearSolver LS,
+ SUNMatrix,
+ N_Vector x,
+ N_Vector b,
+ realtype tol)
+ {
+ auto content = access_content<VectorType>(LS);
+ ARKode<VectorType> &solver = *(content->solver);
+
+ GrowingVectorMemory<VectorType> mem;
+ typename VectorMemory<VectorType>::Pointer rhs(mem);
+ solver.reinit_vector(*rhs);
+ typename VectorMemory<VectorType>::Pointer dst(mem);
+ solver.reinit_vector(*dst);
+
+ copy(*rhs, b);
+
+ SundialsOperator<VectorType> op(solver,
+ content->A_data,
+ content->a_times_fn);
+ SundialsPreconditioner<VectorType> preconditioner(
+ solver, content->P_data, content->preconditioner_solve, tol);
+
+ int err = content->lsolve(op, preconditioner, *dst, *rhs, tol);
+ copy(x, *dst);
+ return err;
+ }
+
+
+
+ template <typename VectorType>
+ int
+ arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix)
+ {
+ auto content = access_content<VectorType>(LS);
+ if (content->preconditioner_setup)
+ return content->preconditioner_setup(content->P_data);
+ return 0;
+ }
+
+
+
+ template <typename VectorType>
+ 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 <typename VectorType>
+ int
+ arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes)
+ {
+ auto content = access_content<VectorType>(LS);
+ content->A_data = A_data;
+ content->a_times_fn = ATimes;
+ return 0;
+ }
+
+
+
+ template <typename VectorType>
+ int
+ arkode_linsol_set_preconditioner(SUNLinearSolver LS,
+ void * P_data,
+ PSetupFn p_setup,
+ PSolveFn p_solve)
+ {
+ auto content = access_content<VectorType>(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 <typename VectorType>
+ class SundialsLinearSolverWrapper
+ {
+ public:
+ SundialsLinearSolverWrapper(ARKode<VectorType> & solver,
+ LinearSolveFunction<VectorType> lsolve)
+ {
+ sun_linear_solver = SUNLinSolNewEmpty();
+ sun_linear_solver->ops->gettype = arkode_linsol_get_type;
+ sun_linear_solver->ops->solve = arkode_linsol_solve<VectorType>;
+ sun_linear_solver->ops->setup = arkode_linsol_setup<VectorType>;
+ sun_linear_solver->ops->initialize = arkode_linsol_initialize<VectorType>;
+ sun_linear_solver->ops->setatimes = arkode_linsol_set_a_times<VectorType>;
+ sun_linear_solver->ops->setpreconditioner =
+ arkode_linsol_set_preconditioner<VectorType>;
+
+ 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<VectorType> content;
+ };
+
+# endif
+
template <typename VectorType>
ARKode<VectorType>::ARKode(const AdditionalData &data,
const MPI_Comm & mpi_comm)
ARKode<VectorType>::~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<VectorType>::value == false)
{
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);
return time.get_step_number();
}
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
template <typename VectorType>
void
ARKode<VectorType>::reset(const double current_time,
// Free the vectors which are no longer used.
if (yy)
{
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::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<VectorType>::value == false)
{
const IndexSet is = solution.locally_owned_elements();
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);
if (setup_jacobian)
{
ARKode_mem->ark_lsetup = t_arkode_setup_jacobian<VectorType>;
-# 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
if (setup_mass)
{
ARKode_mem->ark_msetup = t_arkode_setup_mass<VectorType>;
-# 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 <typename VectorType>
+ void
+ ARKode<VectorType>::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<VectorType> : nullptr,
+ implicit_function ? &t_arkode_implicit_function<VectorType> : 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 <typename VectorType>
+ void
+ ARKode<VectorType>::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<SundialsLinearSolverWrapper<VectorType>>(
+ *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<VectorType> :
+ nullptr,
+ t_arkode_jac_times_vec_function<VectorType>);
+ AssertARKode(status);
+ if (jacobian_preconditioner_solve)
+ {
+ status = ARKStepSetPreconditioner(
+ arkode_mem,
+ jacobian_preconditioner_setup ?
+ t_arkode_prec_setup_function<VectorType> :
+ nullptr,
+ t_arkode_prec_solve_function<VectorType>);
+ 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 <typename VectorType>
+ void
+ ARKode<VectorType>::setup_mass_solver()
+ {
+ int status;
+ (void)status;
+
+ if (mass_times_vector)
+ {
+ SUNLinearSolver sun_mass_linear_solver;
+ if (solve_mass)
+ {
+ mass_solver =
+ std::make_unique<SundialsLinearSolverWrapper<VectorType>>(
+ *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<VectorType> :
+ nullptr,
+ t_arkode_mass_times_vec_function<VectorType>,
+ this);
+ AssertARKode(status);
+
+ if (mass_preconditioner_solve)
+ {
+ status = ARKStepSetMassPreconditioner(
+ arkode_mem,
+ mass_preconditioner_setup ?
+ t_arkode_mass_prec_setup_function<VectorType> :
+ nullptr,
+ t_arkode_mass_prec_solve_function<VectorType>);
+ AssertARKode(status);
+ }
+ }
+ }
+# endif
+
+
+
template <typename VectorType>
void
ARKode<VectorType>::set_functions_to_trigger_an_assert()
};
}
+
+
+ template <typename VectorType>
+ N_Vector
+ ARKode<VectorType>::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<VectorType>::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 <typename VectorType>
+ void
+ ARKode<VectorType>::free_vector(N_Vector vector) const
+ {
+# ifdef DEAL_II_WITH_MPI
+ if (is_serial_vector<VectorType>::value == false)
+ N_VDestroy_Parallel(vector);
+ else
+# endif
+ N_VDestroy_Serial(vector);
+ }
+
+
+
+ template <typename VectorType>
+ void *
+ ARKode<VectorType>::get_arkode_memory() const
+ {
+ return arkode_mem;
+ }
+
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+
+ template <typename VectorType>
+ SundialsOperator<VectorType>::SundialsOperator(ARKode<VectorType> &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 <typename VectorType>
+ void
+ SundialsOperator<VectorType>::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 <typename VectorType>
+ SundialsPreconditioner<VectorType>::SundialsPreconditioner(
+ ARKode<VectorType> &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 <typename VectorType>
+ void
+ SundialsPreconditioner<VectorType>::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<Vector<double>>;
template class ARKode<BlockVector<double>>;
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+ template struct SundialsOperator<Vector<double>>;
+ template struct SundialsOperator<BlockVector<double>>;
+
+ template struct SundialsPreconditioner<Vector<double>>;
+ template struct SundialsPreconditioner<BlockVector<double>>;
+# endif
+
# ifdef DEAL_II_WITH_MPI
# ifdef DEAL_II_WITH_TRILINOS
template class ARKode<TrilinosWrappers::MPI::Vector>;
template class ARKode<TrilinosWrappers::MPI::BlockVector>;
+
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+ template struct SundialsOperator<TrilinosWrappers::MPI::Vector>;
+ template struct SundialsOperator<TrilinosWrappers::MPI::BlockVector>;
+
+ template struct SundialsPreconditioner<TrilinosWrappers::MPI::Vector>;
+ template struct SundialsPreconditioner<TrilinosWrappers::MPI::BlockVector>;
+# endif
# endif // DEAL_II_WITH_TRILINOS
# ifdef DEAL_II_WITH_PETSC
# ifndef PETSC_USE_COMPLEX
template class ARKode<PETScWrappers::MPI::Vector>;
template class ARKode<PETScWrappers::MPI::BlockVector>;
+
+# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
+ template class SundialsOperator<PETScWrappers::MPI::Vector>;
+ template class SundialsOperator<PETScWrappers::MPI::BlockVector>;
+
+ template class SundialsPreconditioner<PETScWrappers::MPI::Vector>;
+ template class SundialsPreconditioner<PETScWrappers::MPI::BlockVector>;
+# endif
# endif // PETSC_USE_COMPLEX
# endif // DEAL_II_WITH_PETSC
#ifdef DEAL_II_WITH_SUNDIALS
-# include <deal.II/base/utilities.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-# include <deal.II/lac/block_vector.h>
-# ifdef DEAL_II_WITH_TRILINOS
-# include <deal.II/lac/trilinos_parallel_block_vector.h>
-# include <deal.II/lac/trilinos_vector.h>
-# endif
-# ifdef DEAL_II_WITH_PETSC
-# include <deal.II/lac/petsc_block_vector.h>
-# include <deal.II/lac/petsc_vector.h>
-# endif
+# include <deal.II/base/utilities.h>
+
+# include <deal.II/lac/block_vector.h>
+# ifdef DEAL_II_WITH_TRILINOS
+# include <deal.II/lac/trilinos_parallel_block_vector.h>
+# include <deal.II/lac/trilinos_vector.h>
+# endif
+# ifdef DEAL_II_WITH_PETSC
+# include <deal.II/lac/petsc_block_vector.h>
+# include <deal.II/lac/petsc_vector.h>
+# endif
-# include <deal.II/sundials/copy.h>
+# include <deal.II/sundials/copy.h>
-# ifdef DEAL_II_SUNDIALS_WITH_IDAS
-# include <idas/idas_impl.h>
-# else
-# include <ida/ida_impl.h>
-# endif
+# ifdef DEAL_II_SUNDIALS_WITH_IDAS
+# include <idas/idas_impl.h>
+# else
+# include <ida/ida_impl.h>
+# endif
-# include <iomanip>
-# include <iostream>
+# include <iomanip>
+# include <iostream>
DEAL_II_NAMESPACE_OPEN
{
if (ida_mem)
IDAFree(&ida_mem);
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
const int ierr = MPI_Comm_free(&communicator);
(void)ierr;
AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
}
-# endif
+# endif
}
// 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<VectorType>::value == false)
{
const IndexSet is = solution.locally_owned_elements();
N_VNew_Parallel(communicator, local_system_size, system_size);
}
else
-# endif
+# endif
{
Assert(is_serial_vector<VectorType>::value,
ExcInternalError(
}
// Free the vectors which are no longer used.
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
N_VDestroy_Parallel(yy);
N_VDestroy_Parallel(diff_id);
}
else
-# endif
+# endif
{
N_VDestroy_Serial(yy);
N_VDestroy_Serial(yp);
// Free the vectors which are no longer used.
if (yy)
{
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
N_VDestroy_Parallel(yy);
N_VDestroy_Parallel(diff_id);
}
else
-# endif
+# endif
{
N_VDestroy_Serial(yy);
N_VDestroy_Serial(yp);
int status;
(void)status;
system_size = solution.size();
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
const IndexSet is = solution.locally_owned_elements();
N_VNew_Parallel(communicator, local_system_size, system_size);
}
else
-# endif
+# endif
{
yy = N_VNew_Serial(system_size);
yp = N_VNew_Serial(system_size);
IDA_mem->ida_lsetup = t_dae_lsetup<VectorType>;
IDA_mem->ida_lsolve = t_dae_solve<VectorType>;
-# 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);
template class IDA<Vector<double>>;
template class IDA<BlockVector<double>>;
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
-# ifdef DEAL_II_WITH_TRILINOS
+# ifdef DEAL_II_WITH_TRILINOS
template class IDA<TrilinosWrappers::MPI::Vector>;
template class IDA<TrilinosWrappers::MPI::BlockVector>;
-# endif // DEAL_II_WITH_TRILINOS
+# endif // DEAL_II_WITH_TRILINOS
-# ifdef DEAL_II_WITH_PETSC
-# ifndef PETSC_USE_COMPLEX
+# ifdef DEAL_II_WITH_PETSC
+# ifndef PETSC_USE_COMPLEX
template class IDA<PETScWrappers::MPI::Vector>;
template class IDA<PETScWrappers::MPI::BlockVector>;
-# endif // PETSC_USE_COMPLEX
-# endif // DEAL_II_WITH_PETSC
+# endif // PETSC_USE_COMPLEX
+# endif // DEAL_II_WITH_PETSC
-# endif // DEAL_II_WITH_MPI
+# endif // DEAL_II_WITH_MPI
} // namespace SUNDIALS
DEAL_II_NAMESPACE_CLOSE
+# endif
+
#endif // DEAL_II_WITH_SUNDIALS
#ifdef DEAL_II_WITH_SUNDIALS
-# include <deal.II/base/utilities.h>
+# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
-# include <deal.II/lac/block_vector.h>
-# ifdef DEAL_II_WITH_TRILINOS
-# include <deal.II/lac/trilinos_parallel_block_vector.h>
-# include <deal.II/lac/trilinos_vector.h>
-# endif
-# ifdef DEAL_II_WITH_PETSC
-# include <deal.II/lac/petsc_block_vector.h>
-# include <deal.II/lac/petsc_vector.h>
-# endif
+# include <deal.II/base/utilities.h>
+
+# include <deal.II/lac/block_vector.h>
+# ifdef DEAL_II_WITH_TRILINOS
+# include <deal.II/lac/trilinos_parallel_block_vector.h>
+# include <deal.II/lac/trilinos_vector.h>
+# endif
+# ifdef DEAL_II_WITH_PETSC
+# include <deal.II/lac/petsc_block_vector.h>
+# include <deal.II/lac/petsc_vector.h>
+# endif
-# include <deal.II/sundials/copy.h>
+# include <deal.II/sundials/copy.h>
-# include <sundials/sundials_config.h>
-# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0)
-# include <kinsol/kinsol_direct.h>
-# include <sunlinsol/sunlinsol_dense.h>
-# include <sunmatrix/sunmatrix_dense.h>
-# else
-# include <kinsol/kinsol_dense.h>
-# endif
+# include <sundials/sundials_config.h>
+# if DEAL_II_SUNDIALS_VERSION_GTE(3, 0, 0)
+# include <kinsol/kinsol_direct.h>
+# include <sunlinsol/sunlinsol_dense.h>
+# include <sunmatrix/sunmatrix_dense.h>
+# else
+# include <kinsol/kinsol_dense.h>
+# endif
-# include <iomanip>
-# include <iostream>
+# include <iomanip>
+# include <iostream>
DEAL_II_NAMESPACE_OPEN
{
if (kinsol_mem)
KINFree(&kinsol_mem);
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
if (is_serial_vector<VectorType>::value == false)
{
const int ierr = MPI_Comm_free(&communicator);
(void)ierr;
AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
}
-# endif
+# endif
}
// 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<VectorType>::value == false)
{
const IndexSet is = initial_guess_and_solution.locally_owned_elements();
N_VConst_Parallel(1.e0, f_scale);
}
else
-# endif
+# endif
{
Assert(is_serial_vector<VectorType>::value,
ExcInternalError(
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)
{
if (setup_jacobian)
{
KIN_mem->kin_lsetup = t_kinsol_setup_jacobian<VectorType>;
-# 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);
}
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<VectorType>::value == false)
{
N_VDestroy_Parallel(solution);
N_VDestroy_Parallel(f_scale);
}
else
-# endif
+# endif
{
N_VDestroy_Serial(solution);
N_VDestroy_Serial(u_scale);
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<unsigned int>(nniters);
template class KINSOL<Vector<double>>;
template class KINSOL<BlockVector<double>>;
-# ifdef DEAL_II_WITH_MPI
+# ifdef DEAL_II_WITH_MPI
-# ifdef DEAL_II_WITH_TRILINOS
+# ifdef DEAL_II_WITH_TRILINOS
template class KINSOL<TrilinosWrappers::MPI::Vector>;
template class KINSOL<TrilinosWrappers::MPI::BlockVector>;
-# endif
+# endif
-# ifdef DEAL_II_WITH_PETSC
-# ifndef PETSC_USE_COMPLEX
+# ifdef DEAL_II_WITH_PETSC
+# ifndef PETSC_USE_COMPLEX
template class KINSOL<PETScWrappers::MPI::Vector>;
template class KINSOL<PETScWrappers::MPI::BlockVector>;
+# endif
# endif
-# endif
-# endif
+# endif
} // namespace SUNDIALS
DEAL_II_NAMESPACE_CLOSE
-
+# endif
#endif
int
main(int argc, char **argv)
{
- std::ofstream out("output");
+ initlog();
Utilities::MPI::MPI_InitFinalize mpi_initialization(
argc, argv, numbers::invalid_unsigned_int);
// 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<VectorType> ode(data);
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;
};
--- /dev/null
+
+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
--- /dev/null
+
+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
int
main(int argc, char **argv)
{
- std::ofstream out("output");
+ initlog();
Utilities::MPI::MPI_InitFinalize mpi_initialization(
argc, argv, numbers::invalid_unsigned_int);
SUNDIALS::ARKode<VectorType>::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<VectorType> ode(data);
// 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;
};
--- /dev/null
+
+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
--- /dev/null
+
+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
--- /dev/null
+
+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
int
main(int argc, char **argv)
{
- std::ofstream out("output");
+ initlog();
Utilities::MPI::MPI_InitFinalize mpi_initialization(
argc, argv, numbers::invalid_unsigned_int);
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<VectorType> ode(data);
// 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;
};
--- /dev/null
+
+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
--- /dev/null
+
+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
int
main(int argc, char **argv)
{
- std::ofstream out("output");
+ initlog();
Utilities::MPI::MPI_InitFinalize mpi_initialization(
argc, argv, numbers::invalid_unsigned_int);
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<VectorType> ode(data);
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;
};
--- /dev/null
+
+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
int
main(int argc, char **argv)
{
- std::ofstream out("output");
+ initlog();
Utilities::MPI::MPI_InitFinalize mpi_initialization(
argc, argv, numbers::invalid_unsigned_int);
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<VectorType> ode(data);
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;
};
--- /dev/null
+
+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
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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 <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/arkode.h>
+
+#include <arkode/arkode_arkstep.h>
+
+#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<double> VectorType;
+
+ ParameterHandler prm;
+ SUNDIALS::ARKode<VectorType>::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<VectorType> 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<double> 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<VectorType> &op,
+ SUNDIALS::SundialsPreconditioner<VectorType> &,
+ VectorType & x,
+ const VectorType &b,
+ double tol) -> int {
+ ReductionControl control;
+ SolverCG<VectorType> 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<double> y(3);
+ y[0] = u0;
+ y[1] = v0;
+ y[2] = w0;
+ ode.solve_ode(y);
+ return 0;
+}
--- /dev/null
+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
--- /dev/null
+
+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
--- /dev/null
+
+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
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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 <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/arkode.h>
+
+#include <arkode/arkode_arkstep.h>
+
+#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<double> VectorType;
+
+ ParameterHandler prm;
+ SUNDIALS::ARKode<VectorType>::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<VectorType> 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<double> 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<VectorType> & op,
+ SUNDIALS::SundialsPreconditioner<VectorType> &prec,
+ VectorType & x,
+ const VectorType & b,
+ double tol) -> int {
+ ReductionControl control;
+ SolverCG<VectorType> 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<double> y(3);
+ y[0] = u0;
+ y[1] = v0;
+ y[2] = w0;
+ ode.solve_ode(y);
+ return 0;
+}
--- /dev/null
+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
--- /dev/null
+
+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
--- /dev/null
+
+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
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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 <deal.II/base/logstream.h>
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/arkode.h>
+
+#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<double> VectorType;
+
+ ParameterHandler prm;
+ SUNDIALS::ARKode<VectorType>::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<VectorType> ode(data);
+
+ ode.reinit_vector = [&](VectorType &v) { v.reinit(2); };
+
+ // Explicit jacobian = stiffness matrix
+ FullMatrix<double> K(2, 2);
+ K(0, 0) = K(1, 1) = 0.5;
+ K(1, 0) = K(0, 1) = -0.5;
+
+ // mass matrix
+ FullMatrix<double> 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<VectorType> & op,
+ SUNDIALS::SundialsPreconditioner<VectorType> &prec,
+ VectorType & x,
+ const VectorType & b,
+ double tol) -> int {
+ ReductionControl control;
+ SolverCG<VectorType> solver_cg(control);
+ solver_cg.solve(op, x, b, prec);
+ return 0;
+ };
+
+ ode.solve_linearized_system = solve_function;
+
+ ode.solve_mass = solve_function;
+
+ FullMatrix<double> 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<double> y(2);
+ y[0] = 1;
+ y[1] = 0;
+ ode.solve_ode(y);
+ return 0;
+}
--- /dev/null
+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
--- /dev/null
+
+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
--- /dev/null
+
+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
+++ /dev/null
-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
+++ /dev/null
-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
+++ /dev/null
-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
+++ /dev/null
-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
+++ /dev/null
-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
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);
main()
{
initlog();
-
+#if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
SUNDIALS::IDA<Vector<double>> ida_solver;
(void)ida_solver;
deallog << "IDA OK" << std::endl;
SUNDIALS::KINSOL<Vector<double>> kinsol_solver;
(void)kinsol_solver;
deallog << "KINSOL OK" << std::endl;
+#endif
SUNDIALS::ARKode<Vector<double>> arkode_solver;
(void)arkode_solver;
--- /dev/null
+
+DEAL::ARKODE OK
projects / dealii.git / commitdiff