From 371407d69c18467ca7772105e73c5acdab509638 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 12 Sep 2017 13:43:43 +0200 Subject: [PATCH] Moved options to AdditionalData. Fixed doc. --- include/deal.II/sundials/ida.h | 585 ++++++++++++++-------- source/sundials/ida.cc | 167 ++---- tests/sundials/harmonic_oscillator_01.cc | 10 +- tests/sundials/harmonic_oscillator_01.prm | 36 +- 4 files changed, 433 insertions(+), 365 deletions(-) diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index ed00b4b198..a3a7c3a8b4 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -37,6 +37,7 @@ #include #include +#include #include @@ -49,11 +50,11 @@ namespace SUNDIALS { /** - * Interface to SUNDIALS IDA (Implicit Differential-Algebraic) solver. + * Interface to SUNDIALS Implicit Differential-Algebraic (IDA) solver. * - * The class IDAInterface is a wrapper to the Implicit - * Differential-Algebraic solver which is a general purpose solver for - * systems of Differential-Algebraic Equations (DAEs). + * The class IDA is a wrapper to SUNDIALS Implicit Differential-Algebraic + * solver which is a general purpose solver for systems of + * Differential-Algebraic Equations (DAEs). * * The user has to provide the implementation of the following std::functions: * - reinit_vector; @@ -67,9 +68,11 @@ namespace SUNDIALS * thrown. * - solver_should_restart; * - differential_components; - * - output_step; * - get_local_tolerances; * + * To output steps, connect a function to the signal + * - output_step; + * * Citing from the SUNDIALS documentation: * * Consider a system of Differential-Algebraic Equations written in the @@ -143,7 +146,6 @@ namespace SUNDIALS * \f] * * That is $F(y', y, t) = y' + A y = 0 $ - * * where * \f[ * \begin{matrix} @@ -151,13 +153,12 @@ namespace SUNDIALS * k^2 &0 * \end{matrix} * \f] - * - * and $y(0)=(0, k)$, $y'(0) = (k, 0)$ + * and $y(0)=(0, k)$, $y'(0) = (k, 0)$. * * The exact solution is $y_0(t) = \sin(k t)$, $y_1(t) = y_0'(t) = k \cos(k t)$, - * $y_1'(t) = -k^2 \sin(k t)$ + * $y_1'(t) = -k^2 \sin(k t)$. * - * The Jacobian to assemble is the following: $J = \alpha I + A$ + * The Jacobian to assemble is the following: $J = \alpha I + A$. * * This is achieved by the following snippet of code: * @code @@ -226,8 +227,308 @@ namespace SUNDIALS public: /** - * Constructor. It is possible to fine tune the SUNDIALS IDA solver by tweaking some of - * the input parameters. + * Additional parameters that can be passed to the IDA class. + */ + class AdditionalData + { + public: + /** + * IDA is a Differential Algebraic solver. As such, it requires initial + * conditions also for the first order derivatives. If you do not provide + * consistent initial conditions, (i.e., conditions for which F(y_dot(0), + * y(0), 0) = 0), you can ask SUNDIALS to compute initial conditions for + * you by specifying InitialConditionCorrection for the initial + * conditions both at the `initial_time` (`ic_type`) and after a reset + * has occurred (`reset_type`). + */ + enum InitialConditionCorrection + { + /** + * Do not try to make initial conditions consistent. + */ + none = 0, + + /** + * Compute the algebraic components of y and differential + * components of y_dot, given the differential components of y. + * This option requires that the user specifies differential and + * algebraic components in the function get_differential_components. + */ + use_y_diff = 1, + + /** + * Compute all components of y, given y_dot. + */ + use_y_dot = 2 + }; + + /** + * Initialization parameters for IDA. + * + * Global parameters: + * + * @param initial_time Initial time + * @param final_time Final time + * @param initial_step_size Initial step size + * @param output_period Time interval between each output + * + * Running parameters: + * + * @param minimum_step_size Minimum step size + * @param maximum_order Maximum BDF order + * @param maximum_non_linear_iterations Maximum number of nonlinear iterations + * + * Error parameters: + * + * @param absolute_tolerance Absolute error tolerance + * @param relative_tolerance Relative error tolerance + * @param use_local_tolerances Use local tolerances when computing errors + * @param ignore_algebraic_terms_for_errors Ignore algebraic terms for error computations + * + * Initial condition correction parameters: + * + * @param ic_type Initial condition correction type + * @param reset_type Initial condition correction type after restart + * @param maximum_non_linear_iterations_ic Initial condition Newton max iterations + */ + AdditionalData(// Initial parameters + const double &initial_time = 0.0, + const double &final_time = 1.0, + const double &initial_step_size = 1e-2, + const double &output_period = 1e-1, + // Running parameters + const double &minimum_step_size = 1e-6, + const unsigned int &maximum_order = 5, + const unsigned int &maximum_non_linear_iterations = 10, + // Error parameters + const double &absolute_tolerance = 1e-6, + const double &relative_tolerance = 1e-5, + const bool &use_local_tolerances = false, + const bool &ignore_algebraic_terms_for_errors = true, + // Initial conditions parameters + const InitialConditionCorrection &ic_type = use_y_diff, + const InitialConditionCorrection &reset_type = use_y_diff, + const unsigned int &maximum_non_linear_iterations_ic = 5) : + initial_time(initial_time), + final_time(final_time), + initial_step_size(initial_step_size), + minimum_step_size(minimum_step_size), + absolute_tolerance(absolute_tolerance), + relative_tolerance(relative_tolerance), + maximum_order(maximum_order), + output_period(output_period), + ignore_algebraic_terms_for_errors(ignore_algebraic_terms_for_errors), + ic_type(ic_type), + reset_type(reset_type), + maximum_non_linear_iterations_ic(maximum_non_linear_iterations_ic), + maximum_non_linear_iterations(maximum_non_linear_iterations), + use_local_tolerances(use_local_tolerances) + {}; + + /** + * Add all AdditionalData() parameters to the given ParameterHandler + * 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 BDF = 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 updated + * whenever the content of `prm` is updated. + * + * Make sure that this class lives longer than `prm`. Undefined behaviour + * will occurr if you destroy this class, and then parse a parameter file + * using `prm`. + */ + void add_parameters(ParameterHandler &prm) + { + prm.add_parameter("Initial time", initial_time); + prm.add_parameter("Final time", final_time); + prm.add_parameter("Time interval between each output", output_period); + + prm.enter_subsection("Running parameters"); + prm.add_parameter("Initial step size",initial_step_size); + prm.add_parameter("Minimum step size", minimum_step_size); + prm.add_parameter("Maximum order of BDF", maximum_order); + prm.add_parameter("Maximum number of nonlinear iterations", maximum_non_linear_iterations); + prm.leave_subsection(); + + prm.enter_subsection("Error control"); + prm.add_parameter("Absolute error tolerance", absolute_tolerance); + prm.add_parameter("Relative error tolerance", relative_tolerance); + prm.add_parameter("Ignore algebraic terms for error computations", ignore_algebraic_terms_for_errors, + "Indicate whether or not to suppress algebraic variables " + "in the local error test."); + prm.add_parameter("Use local tolerances", use_local_tolerances); + prm.leave_subsection(); + + prm.enter_subsection("Initial condition correction parameters"); + static std::string ic_type_str="use_y_diff"; + prm.add_parameter("Correction type at initial time", ic_type_str, + "This is one of the following three options for the " + "initial condition calculation. \n" + " none: do not try to make initial conditions consistent. \n" + " use_y_diff: compute the algebraic components of y and differential\n" + " components of y_dot, given the differential components of y. \n" + " This option requires that the user specifies differential and \n" + " algebraic components in the function get_differential_components.\n" + " use_y_dot: compute all components of y, given y_dot.", + Patterns::Selection("none|use_y_diff|use_y_dot")); + prm.add_action("Correction type at initial time",[&](const std::string &value) + { + if (value == "use_y_diff") + ic_type = use_y_diff; + else if (value == "use_y_dot") + ic_type = use_y_dot; + else if (value == "none") + ic_type = none; + else + AssertThrow(false, ExcInternalError()); + }); + + static std::string reset_type_str="use_y_diff"; + prm.add_parameter("Correction type after restart", reset_type_str, + "This is one of the following three options for the " + "initial condition calculation. \n" + " none: do not try to make initial conditions consistent. \n" + " use_y_diff: compute the algebraic components of y and differential\n" + " components of y_dot, given the differential components of y. \n" + " This option requires that the user specifies differential and \n" + " algebraic components in the function get_differential_components.\n" + " use_y_dot: compute all components of y, given y_dot.", + Patterns::Selection("none|use_y_diff|use_y_dot")); + prm.add_action("Correction type after restart",[&](const std::string &value) + { + if (value == "use_y_diff") + reset_type = use_y_diff; + else if (value == "use_y_dot") + reset_type = use_y_dot; + else if (value == "none") + reset_type = none; + else + AssertThrow(false, ExcInternalError()); + }); + prm.add_parameter("Maximum number of nonlinear iterations", maximum_non_linear_iterations_ic); + prm.leave_subsection(); + } + + /** + * Initial time for the DAE. + */ + double initial_time; + + /** + * Final time. + */ + double final_time; + + /** + * Initial step size. + */ + double initial_step_size; + + /** + * Minimum step size. + */ + double minimum_step_size; + + /** + * Absolute error tolerance for adaptive time stepping. + */ + double absolute_tolerance; + + /** + * Relative error tolerance for adaptive time stepping. + */ + double relative_tolerance; + + /** + * Maximum order of BDF. + */ + unsigned int maximum_order; + + /** + * Time period between each output. + */ + double output_period; + + /** + * Ignore algebraic terms for errors. + */ + bool ignore_algebraic_terms_for_errors; + + /** + * Type of correction for initial conditions. + * + * If you do not provide consistent initial conditions, (i.e., conditions + * for which $F(y_dot(0), y(0), 0) = 0$), you can ask SUNDIALS to compute + * initial conditions for you by using the `ic_type` parameter at + * construction time. + * + * Notice that you could in principle use this capabilities to solve for + * stady state problems by setting y_dot to zero, and asking to compute + * $y(0)$ that satisfies $F(0, y(0), 0) = 0$, however the nonlinear solver + * used inside IDA may not be robust enough for complex problems with + * several millions unknowns. + */ + InitialConditionCorrection ic_type; + + /** + * Type of correction for initial conditions to be used after a solver + * restart. + * + * If you do not have consistent initial conditions after a restart, (i.e., + * conditions for which F(y_dot(t_restart), y(t_restart), t_restart) = 0), + * you can ask SUNDIALS to compute the new initial conditions for you by + * using the `reset_type` parameter at construction time. + */ + InitialConditionCorrection reset_type; + + /** + * Maximum number of iterations for Newton method in IC calculation. + */ + unsigned maximum_non_linear_iterations_ic; + + /** + * Maximum number of iterations for Newton method during time advancement. + */ + unsigned int maximum_non_linear_iterations; + + /** + * Use local tolerances when computing absolute tolerance. + */ + bool use_local_tolerances; + }; + + /** + * Constructor. It is possible to fine tune the SUNDIALS IDA solver by passing an + * AdditionalData() object that sets all of the solver parameters. * * IDA is a Differential Algebraic solver. As such, it requires initial conditions also for * the first order derivatives. If you do not provide consistent initial conditions, (i.e., @@ -257,81 +558,19 @@ namespace SUNDIALS * options that you used for the initial conditions in `reset_type`. * * @param mpi_comm MPI communicator - * @param initial_time Initial time - * @param final_time Final time - * @param initial_step_size Initial step size - * @param min_step_size Minimum step size - * @param abs_tol Absolute error tolerance - * @param rel_tol Relative error tolerance - * @param max_order Maximum BDF order - * @param output_period Time interval between each output - * @param ignore_algebraic_terms_for_errors Ignore algebraic terms for error computations - * @param ic_type Initial condition type - * @param reset_type Initial condition type after restart - * @param ic_alpha Initial condition Newton parameter - * @param ic_max_iter Initial condition Newton max iterations - * @param max_non_linear_iterations Maximum number of nonlinear iterations - * @param verbose Show output of time steps - * @param use_local_tolerances Use local tolerances + * @param data IDA configuration data */ IDA(const MPI_Comm mpi_comm = MPI_COMM_WORLD, - const double &initial_time = 0.0, - const double &final_time = 1.0, - const double &initial_step_size = 1e-4, - const double &min_step_size = 1e-6, - const double &abs_tol = 1e-6, - const double &rel_tol = 1e-5, - const unsigned int &max_order = 5, - const double &output_period = .1, - const bool &ignore_algebraic_terms_for_errors = true, - const std::string &ic_type = "none", - const std::string &reset_type = "none", - const double &ic_alpha = .33, - const unsigned int &ic_max_iter = 5, - const unsigned int &max_non_linear_iterations = 10, - const bool &verbose = true, - const bool &use_local_tolerances = false); + const AdditionalData &data=AdditionalData()); /** - * House cleaning. + * Destructor. */ ~IDA(); /** - * Add all internal parameters to the given ParameterHandler object. When - * the paramaters are parsed from a file, the internal parameters are automatically - * updated. - * - * The following parameters are declared: - * - * @code - * set Absolute error tolerance = 0.000001 - * set Final time = 1.000000 - * set Ignore algebraic terms for error computations = true - * set Initial condition Newton max iterations = 5 - * set Initial condition Newton parameter = 0.330000 - * set Initial condition type = use_y_diff - * set Initial condition type after restart = use_y_dot - * set Initial step size = 0.000100 - * set Initial time = 0.000000 - * set Maximum number of nonlinear iterations = 10 - * set Maximum order of BDF = 5 - * set Min step size = 0.000001 - * set Relative error tolerance = 0.000010 - * set Show output of time steps = true - * set Time units between each output = 0.05 - * set Use local tolerances = false - * @endcode - * - * These are one-to-one with the options you can pass at construction time. - * - * Those options are set as default value in the ParameterHandler object - * `prm`. - */ - virtual void add_parameters(ParameterHandler &prm); - - /** - * Evolve. This function returns the final number of steps. + * Integrate differential-algebraic equations. This function returns the + * final number of computed steps. */ unsigned int solve_dae(VectorType &solution, VectorType &solution_dot); @@ -341,11 +580,11 @@ namespace SUNDIALS * called when the simulation start and when the user returns true to a * call to solver_should_restart(). * - * By default solver_should_restart returns false. If the user needs to + * By default solver_should_restart() returns false. If the user needs to * implement, for example, local adaptivity in space, he or she may assign * a different function to solver_should_restart() that performs all mesh - * changes, transfers the solution and the solution dot to the new mesh, and - * returns true. + * changes, transfers the solution and the solution dot to the new mesh, + * and returns true. * * During reset(), both y and yp are checked for consistency, and according to * what was specified as ic_type (if t==initial_time) or reset_type (if @@ -368,7 +607,7 @@ namespace SUNDIALS std::function reinit_vector; /** - * Compute residual. + * Compute residual. Return $F(t, y, \dot y)$. * * This function should return: * - 0: Success @@ -383,7 +622,27 @@ namespace SUNDIALS VectorType &res)> residual; /** - * Compute Jacobian. + * Compute Jacobian. This function is called by IDA any time a Jacobian + * update is required. The user should compute the Jacobian (or update all + * the variables that allow the application of the Jacobian). This function + * is called by IDA once, before any call to solve_jacobian_system(). + * + * The Jacobian $J$ should be a (possibly inexact) computation of + * \f[ + * J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} + + * \alpha \dfrac{\partial F}{\partial \dot y}. + * \f] + * + * If the user uses a matrix based computation of the Jacobian, than this + * is the right place where an assembly routine shoulde be called to + * assemble both a matrix and a preconditioner for the Jacobian system. + * Subsequent calls (possibly more than one) to solve_jacobian_system() can + * assume that this function has been called at least once. + * + * Notice that no assumption is made by this interface on what the user + * should do in this function. IDA only assumes that after a call to + * setup_jacobian() it is possible to call solve_jacobian_system(), to + * obtain a solution $x$ to the system $J x = b$. * * This function should return: * - 0: Success @@ -398,7 +657,25 @@ namespace SUNDIALS const double alpha)> setup_jacobian; /** - * Solve linear system. + * Solve the Jacobian linear system. This function will be called by IDA + * (possibly several times) after setup_jacobian() has been called at least + * once. IDA tries to do its best to call setup_jacobian() the minimum + * amount of times. If convergence can be achieved without updating the + * Jacobian, then IDA does not call setup_jacobian() again. If, on the + * contrary, internal IDA convergence tests fail, then IDA calls again + * setup_jacobian() with updated vectors and coefficents so that successive + * calls to solve_jacobian_systems() lead to better convergence in the + * Newton process. + * + * The jacobian $J$ should be (an approximation of) the system Jacobian + * \f[ + * J=\dfrac{\partial G}{\partial y} = \dfrac{\partial F}{\partial y} + + * \alpha \dfrac{\partial F}{\partial \dot y}. + * \f] + * + * A call to this function should store in `dst` the result of $J^{-1}$ + * applied to `src`, i.e., `J*dst = src`. It is the users responsability to + * set up proper solvers and preconditioners inside this function. * * This function should return: * - 0: Success @@ -410,7 +687,18 @@ namespace SUNDIALS std::function solve_jacobian_system; /** - * Store solutions to file. + * Process solution. This function is called by IDA at fixed time steps, + * every `output_period` seconds, and it is passed a polynomial + * interpolation of the solution and of its time derivative, computed using + * the current BDF order and the (internally stored) previously computed + * solution steps. + * + * Notice that it is well possible that internally IDA computes a time step + * which is much larger than the `output_period` step, and therefore calls + * this function consecutively several times by simply performing all + * intermediate interpolations. There is no relationship between how many + * times this function is called and how many time steps have actually been + * computed. */ std::function get_local_tolerances; - /** - * Set initial time equal to @p t disregarding what is written - * in the parameter file. - */ - void set_initial_time(const double &t); - /** * Handle IDA exceptions. */ @@ -488,116 +770,9 @@ namespace SUNDIALS void set_functions_to_trigger_an_assert(); /** - * Initial time for the DAE. - */ - double initial_time; - - /** - * Final time. - */ - double final_time; - - /** - * Initial step size. - */ - double initial_step_size; - - /** - * Minimum step size. - */ - double min_step_size; - - /** - * Absolute error tolerance for adaptive time stepping. - */ - double abs_tol; - - /** - * Relative error tolerance for adaptive time stepping. - */ - double rel_tol; - - /** - * Maximum order of BDF. - */ - unsigned int max_order; - - /** - * Time period between each output. + * IDA configuration data. */ - double output_period; - - /** - * Ignore algebraic terms for errors. - */ - bool ignore_algebraic_terms_for_errors; - - /** - * Type of initial conditions. - * - * IDA is a Differential Algebraic solver. As such, it requires initial conditions also for - * the first order derivatives. If you do not provide consistent initial conditions, (i.e., - * conditions for which $F(y_dot(0), y(0), 0) = 0$), you can ask SUNDIALS to compute initial - * conditions for you by using the `ic_type` parameter at construction time. - * - * You have three options - * - none: do not try to make initial conditions consistent. - * - use_y_diff: compute the algebraic components of y and differential - * components of y_dot, given the differential components of y. - * This option requires that the user specifies differential and - * algebraic components in the function get_differential_components. - * - use_y_dot: compute all components of y, given y_dot. - * - * Notice that you could in principle use this capabilities to solve for - * stady state problems by setting y_dot to zero, and asking to compute - * $y(0)$ that satisfies $F(0, y(0), 0) = 0$, however the nonlinear solver - * used inside IDA may not be robust enough for complex problems with - * several millions unknowns. - */ - std::string ic_type; - - /** - * Type of conditions to be used after a solver restart. - * - * If you do not have consistent initial conditions after a restart, (i.e., - * conditions for which F(y_dot(t_restart), y(t_restart), t_restart) = 0), - * you can ask SUNDIALS to compute the new initial conditions for you by - * using the `reset_type` parameter at construction time. - * - * You have three options - * - none: do not try to make initial conditions consistent. - * - use_y_diff: compute the algebraic components of y and differential - * components of y_dot, given the differential components of y. - * This option requires that the user specifies differential and - * algebraic components in the function get_differential_components. - * - use_y_dot: compute all components of y, given y_dot. - */ - std::string reset_type; - - /** - * Alpha to use in Newton method for IC calculation. - */ - double ic_alpha; - - /** - * Maximum number of iterations for Newton method in IC calculation. - */ - unsigned ic_max_iter; - - /** - * Maximum number of iterations for Newton method during time advancement. - */ - unsigned int max_non_linear_iterations; - - /** - * Show the progress of the time stepper. - */ - bool verbose; - - /** - * Use local tolerances when computing absolute tolerance. - */ - bool use_local_tolerances; + AdditionalData data; /** * IDA memory object. @@ -631,12 +806,6 @@ namespace SUNDIALS MPI_Comm communicator; #endif - /** - * Output stream. If run in parallel, only processor zero will - * output verbose information about time steps. - */ - ConditionalOStream pcout; - /** * Memory pool of vectors. */ diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index aaa4c6afe3..d4212bff02 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -148,42 +148,10 @@ namespace SUNDIALS } template - IDA::IDA(const MPI_Comm mpi_comm, - const double &initial_time, - const double &final_time, - const double &initial_step_size, - const double &min_step_size, - const double &abs_tol, - const double &rel_tol, - const unsigned int &max_order, - const double &output_period, - const bool &ignore_algebraic_terms_for_errors, - const std::string &ic_type, - const std::string &reset_type, - const double &ic_alpha, - const unsigned int &ic_max_iter, - const unsigned int &max_non_linear_iterations, - const bool &verbose, - const bool &use_local_tolerances) : - initial_time(initial_time), - final_time(final_time), - initial_step_size(initial_step_size), - min_step_size(min_step_size), - abs_tol(abs_tol), - rel_tol(rel_tol), - max_order(max_order), - output_period(output_period), - ignore_algebraic_terms_for_errors(ignore_algebraic_terms_for_errors), - ic_type(ic_type), - reset_type(reset_type), - ic_alpha(ic_alpha), - ic_max_iter(ic_max_iter), - max_non_linear_iterations(max_non_linear_iterations), - verbose(verbose), - use_local_tolerances(use_local_tolerances), + IDA::IDA(const MPI_Comm mpi_comm, const AdditionalData &data) : + data(data), ida_mem(nullptr), - communicator(Utilities::MPI::duplicate_communicator(mpi_comm)), - pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_comm)==0) + communicator(Utilities::MPI::duplicate_communicator(mpi_comm)) { set_functions_to_trigger_an_assert(); } @@ -196,61 +164,6 @@ namespace SUNDIALS MPI_Comm_free(&communicator); } - template - void IDA::add_parameters(ParameterHandler &prm) - { - prm.add_parameter("Initial step size",initial_step_size); - - prm.add_parameter("Minimum step size", min_step_size); - - prm.add_parameter("Absolute error tolerance", abs_tol); - - prm.add_parameter("Relative error tolerance", rel_tol); - - prm.add_parameter("Initial time", initial_time); - - prm.add_parameter("Final time", final_time); - - prm.add_parameter("Time interval between each output", output_period); - - prm.add_parameter("Maximum order of BDF", max_order); - - prm.add_parameter("Maximum number of nonlinear iterations", max_non_linear_iterations); - - prm.add_parameter("Ignore algebraic terms for error computations", ignore_algebraic_terms_for_errors, - "Indicate whether or not to suppress algebraic variables " - "in the local error test."); - - prm.add_parameter("Initial condition type", ic_type, - "This is one of the following three options for the " - "initial condition calculation. \n" - " none: do not try to make initial conditions consistent. \n" - " use_y_diff: compute the algebraic components of y and differential\n" - " components of y_dot, given the differential components of y. \n" - " This option requires that the user specifies differential and \n" - " algebraic components in the function get_differential_components.\n" - " use_y_dot: compute all components of y, given y_dot.", - Patterns::Selection("none|use_y_diff|use_y_dot")); - - prm.add_parameter("Initial condition type after restart", reset_type, - "This is one of the following three options for the " - "initial condition calculation. \n" - " none: do not try to make initial conditions consistent. \n" - " use_y_diff: compute the algebraic components of y and differential\n" - " components of y_dot, given the differential components of y. \n" - " This option requires that the user specifies differential and \n" - " algebraic components in the function get_differential_components.\n" - " use_y_dot: compute all components of y, given y_dot.", - Patterns::Selection("none|use_y_diff|use_y_dot")); - - prm.add_parameter("Initial condition Newton parameter", ic_alpha); - - prm.add_parameter("Initial condition Newton max iterations", ic_max_iter); - - prm.add_parameter("Use local tolerances", use_local_tolerances); - - prm.add_parameter("Show output of time steps", verbose); - } template @@ -261,8 +174,8 @@ namespace SUNDIALS unsigned int system_size = solution.size(); unsigned int local_system_size = system_size; - double t = initial_time; - double h = initial_step_size; + double t = data.initial_time; + double h = data.initial_step_size; unsigned int step_number = 0; int status; @@ -302,26 +215,20 @@ namespace SUNDIALS diff_id = N_VNew_Serial(system_size); abs_tolls = N_VNew_Serial(system_size); } - reset(initial_time, - initial_step_size, + reset(data.initial_time, + data.initial_step_size, solution, solution_dot); - double next_time = initial_time; + double next_time = data.initial_time; output_step( 0, solution, solution_dot, 0); - while (t " - << std::setw(5) << next_time - << std::endl; - } + next_time += data.output_period; + status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL); AssertIDA(status); @@ -339,7 +246,6 @@ namespace SUNDIALS output_step(t, solution, solution_dot, step_number); } - pcout << std::endl; // Free the vectors which are no longer used. #ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) @@ -370,7 +276,7 @@ namespace SUNDIALS unsigned int system_size; unsigned int local_system_size; - bool first_step = (current_time == initial_time); + bool first_step = (current_time == data.initial_time); if (ida_mem) IDAFree(&ida_mem); @@ -439,15 +345,15 @@ namespace SUNDIALS status = IDAInit(ida_mem, t_dae_residual, current_time, yy, yp); AssertIDA(status); - if (use_local_tolerances) + if (data.use_local_tolerances) { copy(abs_tolls, get_local_tolerances()); - status = IDASVtolerances(ida_mem, rel_tol, abs_tolls); + status = IDASVtolerances(ida_mem, data.relative_tolerance, abs_tolls); AssertIDA(status); } else { - status = IDASStolerances(ida_mem, rel_tol, abs_tol); + status = IDASStolerances(ida_mem, data.relative_tolerance, data.absolute_tolerance); AssertIDA(status); } @@ -457,7 +363,9 @@ namespace SUNDIALS status = IDASetUserData(ida_mem, (void *) this); AssertIDA(status); - if (ic_type == "use_y_diff" || reset_type == "use_y_diff" || ignore_algebraic_terms_for_errors) + if (data.ic_type == AdditionalData::use_y_diff || + data.reset_type == AdditionalData::use_y_diff || + data.ignore_algebraic_terms_for_errors) { VectorType diff_comp_vector(solution); diff_comp_vector = 0.0; @@ -470,14 +378,14 @@ namespace SUNDIALS AssertIDA(status); } - status = IDASetSuppressAlg(ida_mem, ignore_algebraic_terms_for_errors); + status = IDASetSuppressAlg(ida_mem, data.ignore_algebraic_terms_for_errors); AssertIDA(status); // status = IDASetMaxNumSteps(ida_mem, max_steps); - status = IDASetStopTime(ida_mem, final_time); + status = IDASetStopTime(ida_mem, data.final_time); AssertIDA(status); - status = IDASetMaxNonlinIters(ida_mem, max_non_linear_iterations); + status = IDASetMaxNonlinIters(ida_mem, data.maximum_non_linear_iterations); AssertIDA(status); // Initialize solver @@ -488,24 +396,19 @@ namespace SUNDIALS IDA_mem->ida_lsolve = t_dae_solve; IDA_mem->ida_setupNonNull = true; - status = IDASetMaxOrd(ida_mem, max_order); + status = IDASetMaxOrd(ida_mem, data.maximum_order); AssertIDA(status); - std::string type; + typename AdditionalData::InitialConditionCorrection type; if (first_step) - type = ic_type; + type = data.ic_type; else - type = reset_type; + type = data.reset_type; - if (verbose) - { - pcout << "computing consistent initial conditions with the option " - << type - << " please be patient." - << std::endl; - } + status = IDASetMaxNumItersIC(ida_mem, data.maximum_non_linear_iterations_ic); + AssertIDA(status); - if (type == "use_y_dot") + if (type == AdditionalData::use_y_dot) { // (re)initialization of the vectors status = IDACalcIC(ida_mem, IDA_Y_INIT, current_time+current_time_step); @@ -517,7 +420,7 @@ namespace SUNDIALS copy(solution, yy); copy(solution_dot, yp); } - else if (type == "use_y_diff") + else if (type == AdditionalData::use_y_diff) { status = IDACalcIC(ida_mem, IDA_YA_YDP_INIT, current_time+current_time_step); AssertIDA(status); @@ -528,18 +431,6 @@ namespace SUNDIALS copy(solution, yy); copy(solution_dot, yp); } - - if (verbose) - { - pcout << "compute initial conditions: done." - << std::endl; - } - } - - template - void IDA::set_initial_time(const double &t) - { - initial_time = t; } template diff --git a/tests/sundials/harmonic_oscillator_01.cc b/tests/sundials/harmonic_oscillator_01.cc index e8eda44b5f..3f0c3be311 100644 --- a/tests/sundials/harmonic_oscillator_01.cc +++ b/tests/sundials/harmonic_oscillator_01.cc @@ -55,7 +55,8 @@ class HarmonicOscillator { public: - HarmonicOscillator(double _kappa=1.0) : + HarmonicOscillator(double _kappa, const typename SUNDIALS::IDA>::AdditionalData &data) : + time_stepper(MPI_COMM_WORLD, data), y(2), y_dot(2), J(2,2), @@ -140,9 +141,9 @@ int main (int argc, char **argv) { Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int); - HarmonicOscillator ode(1.0); + SUNDIALS::IDA>::AdditionalData data; ParameterHandler prm; - ode.time_stepper.add_parameters(prm); + data.add_parameters(prm); // std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_01.prm"); // prm.print_parameters(ofile, ParameterHandler::ShortText); @@ -150,6 +151,9 @@ int main (int argc, char **argv) std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_01.prm"); prm.parse_input(ifile); + + + HarmonicOscillator ode(1.0, data); ode.run(); return 0; } diff --git a/tests/sundials/harmonic_oscillator_01.prm b/tests/sundials/harmonic_oscillator_01.prm index 81776b7c19..36b33edc44 100644 --- a/tests/sundials/harmonic_oscillator_01.prm +++ b/tests/sundials/harmonic_oscillator_01.prm @@ -1,16 +1,20 @@ -set Absolute error tolerance = 0.000001 -set Final time = 6.3 -set Ignore algebraic terms for error computations = false -set Initial condition Newton max iterations = 5 -set Initial condition Newton parameter = 0.330000 -set Initial condition type = none -set Initial condition type after restart = none -set Initial step size = 0.1 -set Initial time = 0 -set Maximum number of nonlinear iterations = 10 -set Maximum order of BDF = 5 -set Minimum step size = 0.000001 -set Relative error tolerance = 0.00001 -set Show output of time steps = true -set Time interval between each output = 0.2 -set Use local tolerances = false +set Final time = 6.3 +set Initial time = 0 +set Time interval between each output = 0.2 +subsection Error control + set Absolute error tolerance = 1e-6 + set Ignore algebraic terms for error computations = true + set Relative error tolerance = 1e-5 + 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 BDF = 5 + set Minimum step size = 0.000001 +end -- 2.39.5