From f71b09376bf689e1095983ab882e8aa01fbc8951 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Fri, 8 Sep 2017 16:45:02 +0200 Subject: [PATCH] Use GrowingVectorMemory in IDA + example in doc. --- include/deal.II/sundials/ida.h | 130 +++++++++++++++++++++-- source/sundials/ida.cc | 66 ++++++++---- tests/quick_tests/sundials-ida.cc | 4 +- tests/sundials/harmonic_oscillator_01.cc | 15 +-- 4 files changed, 174 insertions(+), 41 deletions(-) diff --git a/include/deal.II/sundials/ida.h b/include/deal.II/sundials/ida.h index a5e1130d30..ed00b4b198 100644 --- a/include/deal.II/sundials/ida.h +++ b/include/deal.II/sundials/ida.h @@ -26,7 +26,7 @@ #include #include #include - +#include #include #include @@ -56,14 +56,14 @@ namespace SUNDIALS * systems of Differential-Algebraic Equations (DAEs). * * The user has to provide the implementation of the following std::functions: - * - create_new_vector; + * - reinit_vector; * - residual; * - setup_jacobian; * - solve_jacobian_system; * - * Optionally, also the following functions could be implemented. By default + * Optionally, also the following functions could be rewritten. By default * they do nothing, or are not required. If you call the constructor in a way - * that requires any of the not-implemented functions, an Assertion will be + * that requires a not-implemented function, an Assertion will be * thrown. * - solver_should_restart; * - differential_components; @@ -125,6 +125,99 @@ namespace SUNDIALS * scalar $\alpha$ changes whenever the step size or method order * changes. * + * To provide a simple example, consider the following harmonic oscillator problem: + * \f[ + * \begin{split} + * u'' & = -k^2 u \\ + * u (0) & = 0 \\ + * u'(0) & = k + * \end{split} + * \f] + * + * We write it in terms of a first order ode: + *\f[ + * \begin{matrix} + * y_0' & -y_1 & = 0 \\ + * y_1' & + k^2 y_0 & = 0 + * \end{matrix} + * \f] + * + * That is $F(y', y, t) = y' + A y = 0 $ + * + * where + * \f[ + * \begin{matrix} + * 0 & -1 \\ + * k^2 &0 + * \end{matrix} + * \f] + * + * 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)$ + * + * The Jacobian to assemble is the following: $J = \alpha I + A$ + * + * This is achieved by the following snippet of code: + * @code + * typedef Vector VectorType; + * + * VectorType y(2); + * VectorType y_dot(2); + * + * FullMatrix A(2,2); + * A(0,1) = -1.0; + * A(1,0) = kappa*kappa; + * + * FullMatrix J(2,2); + * FullMatrix Jinv(2,2); + * + * double kappa = 1.0; + * + * IDA time_stepper; + * + * time_stepper.reinit_vector = [&] (VectorType&v) + * { + * v.reinit(2); + * }; + * + * time_stepper.residual = [&](const double t, + * const VectorType &y, + * const VectorType &y_dot, + * VectorType &res) ->int + * { + * res = y_dot; + * A.vmult_add(res, y); + * return 0; + * }; + * + * time_stepper.setup_jacobian = [&](const double , + * const VectorType &, + * const VectorType &, + * const double alpha) ->int + * { + * J = A; + * + * J(0,0) = alpha; + * J(1,1) = alpha; + * + * Jinv.invert(J); + * return 0; + * }; + * + * time_stepper.solve_jacobian_system = [&](const VectorType &src, + * VectorType &dst) ->int + * { + * Jinv.vmult(dst,src); + * return 0; + * }; + * + * y[1] = kappa; + * y_dot[0] = kappa; + * time_stepper.solve_dae(y,y_dot); + * @endcode + * * @author Luca Heltai, Alberto Sartori, 2017. */ template > @@ -149,6 +242,11 @@ namespace SUNDIALS * algebraic components in the function get_differential_components. * - use_y_dot: compute all components of y, given y_dot. * + * By default, this class assumes that all components are differential, and that you want + * to solve a standard ode. In this case, the initial component type is set to `use_y_diff`, + * so that the `y_dot` at time t=`initial_time` is computed by solving the nonlinear problem + * $F(y_dot, y(t0), t0) = 0$ in the variable `y_dot`. + * * Notice that a Newton solver is used for this computation. The Newton solver parameters * can be tweaked by acting on `ic_alpha` and `ic_max_iter`. * @@ -265,9 +363,9 @@ namespace SUNDIALS VectorType &yp); /** - * Return a newly allocated unique_ptr. + * Reinit vector to have the right size, MPI communicator, etc. */ - std::function()> create_new_vector; + std::function reinit_vector; /** * Compute residual. @@ -331,6 +429,9 @@ namespace SUNDIALS * t. This function should then return true, and change the dimension of * both sol and sol_dot to reflect the new dimension. Since IDA does not * know about the new dimension, an internal reset is necessary. + * + * The default implementation simply returns `false`, i.e., no restart is + * performed during the evolution. */ std::function differential_components; @@ -350,6 +452,9 @@ namespace SUNDIALS * the vector norm. The implementation of this function is optional, and it * is used only when `use_local_tolerances` is set to true at construction * time, or through the parameter file. + * + * If you do not overwrite this function and set `use_local_tolerances` to + * true, an exception will be thrown when trying to start the time stepper. */ std::function get_local_tolerances; @@ -531,6 +636,11 @@ namespace SUNDIALS * output verbose information about time steps. */ ConditionalOStream pcout; + + /** + * Memory pool of vectors. + */ + GrowingVectorMemory mem; }; } diff --git a/source/sundials/ida.cc b/source/sundials/ida.cc index e5b5c312db..aaa4c6afe3 100644 --- a/source/sundials/ida.cc +++ b/source/sundials/ida.cc @@ -49,10 +49,16 @@ namespace SUNDIALS N_Vector rr, void *user_data) { IDA &solver = *static_cast *>(user_data); + GrowingVectorMemory mem; - std::shared_ptr src_yy = solver.create_new_vector(); - std::shared_ptr src_yp = solver.create_new_vector(); - std::shared_ptr residual = solver.create_new_vector(); + VectorType *src_yy = mem.alloc(); + solver.reinit_vector(*src_yy); + + VectorType *src_yp = mem.alloc(); + solver.reinit_vector(*src_yp); + + VectorType *residual = mem.alloc(); + solver.reinit_vector(*residual); copy(*src_yy, yy); copy(*src_yp, yp); @@ -61,6 +67,10 @@ namespace SUNDIALS copy(rr, *residual); + mem.free(src_yy); + mem.free(src_yp); + mem.free(residual); + return err; } @@ -80,9 +90,13 @@ namespace SUNDIALS (void) tmp3; (void) resp; IDA &solver = *static_cast *>(IDA_mem->ida_user_data); + GrowingVectorMemory mem; + + VectorType *src_yy = mem.alloc(); + solver.reinit_vector(*src_yy); - std::shared_ptr src_yy = solver.create_new_vector(); - std::shared_ptr src_yp = solver.create_new_vector(); + VectorType *src_yp = mem.alloc(); + solver.reinit_vector(*src_yp); copy(*src_yy, yy); copy(*src_yp, yp); @@ -91,6 +105,10 @@ namespace SUNDIALS *src_yy, *src_yp, IDA_mem->ida_cj); + + mem.free(src_yy); + mem.free(src_yp); + return err; } @@ -108,14 +126,22 @@ namespace SUNDIALS (void) yp; (void) resp; IDA &solver = *static_cast *>(IDA_mem->ida_user_data); + GrowingVectorMemory mem; + + VectorType *src = mem.alloc(); + solver.reinit_vector(*src); - std::shared_ptr dst = solver.create_new_vector(); - std::shared_ptr src = solver.create_new_vector(); + VectorType *dst = mem.alloc(); + solver.reinit_vector(*dst); copy(*src, b); int err = solver.solve_jacobian_system(*src,*dst); copy(b, *dst); + + mem.free(src); + mem.free(dst); + return err; } @@ -297,9 +323,10 @@ namespace SUNDIALS << std::endl; } status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL); + AssertIDA(status); status = IDAGetLastStep(ida_mem, &h); - AssertThrow(status == 0, ExcMessage("Error in IDA Solver")); + AssertIDA(status); copy(solution, yy); copy(solution_dot, yp); @@ -415,11 +442,13 @@ namespace SUNDIALS if (use_local_tolerances) { copy(abs_tolls, get_local_tolerances()); - status += IDASVtolerances(ida_mem, rel_tol, abs_tolls); + status = IDASVtolerances(ida_mem, rel_tol, abs_tolls); + AssertIDA(status); } else { - status += IDASStolerances(ida_mem, rel_tol, abs_tol); + status = IDASStolerances(ida_mem, rel_tol, abs_tol); + AssertIDA(status); } status = IDASetInitStep(ida_mem, current_time_step); @@ -517,11 +546,9 @@ namespace SUNDIALS void IDA::set_functions_to_trigger_an_assert() { - create_new_vector = []() ->std::unique_ptr + reinit_vector = [](VectorType &) { - std::unique_ptr p; - AssertThrow(false, ExcFunctionNotProvided("create_new_vector")); - return p; + AssertThrow(false, ExcFunctionNotProvided("reinit_vector")); }; residual = [](const double, @@ -567,11 +594,14 @@ namespace SUNDIALS return false; }; - differential_components = []() ->IndexSet + differential_components = [&]() ->IndexSet { - IndexSet i; - AssertThrow(false, ExcFunctionNotProvided("differential_components")); - return i; + GrowingVectorMemory mem; + auto v = mem.alloc(); + reinit_vector(*v); + unsigned int size = v->size(); + mem.free(v); + return complete_index_set(size); }; get_local_tolerances = []() ->VectorType & diff --git a/tests/quick_tests/sundials-ida.cc b/tests/quick_tests/sundials-ida.cc index de6ad8b982..75a15f51cc 100644 --- a/tests/quick_tests/sundials-ida.cc +++ b/tests/quick_tests/sundials-ida.cc @@ -68,9 +68,9 @@ public: diff[0] = 1.0; diff[1] = 1.0; - time_stepper.create_new_vector = [&] () -> std::unique_ptr > + time_stepper.reinit_vector = [&] (Vector &v) { - return std::unique_ptr>(new Vector(2)); + v.reinit(2); }; diff --git a/tests/sundials/harmonic_oscillator_01.cc b/tests/sundials/harmonic_oscillator_01.cc index f048672014..e8eda44b5f 100644 --- a/tests/sundials/harmonic_oscillator_01.cc +++ b/tests/sundials/harmonic_oscillator_01.cc @@ -64,14 +64,14 @@ public: kappa(_kappa), out("output") { - time_stepper.create_new_vector = [&] () -> std::unique_ptr > + typedef Vector VectorType; + + time_stepper.reinit_vector = [&] (VectorType&v) { - return std::unique_ptr>(new Vector(2)); + v.reinit(2); }; - typedef Vector VectorType; - time_stepper.residual = [&](const double t, const VectorType &y, const VectorType &y_dot, @@ -115,13 +115,6 @@ public: << sol[0] << " " << sol[1] << " " << sol_dot[0] << " " << sol_dot[1] << std::endl; return 0; }; - - time_stepper.solver_should_restart = [](const double , - VectorType &, - VectorType &) ->bool - { - return false; - }; } void run() -- 2.39.5