From 17bd91ee62e859d5b3c60cdf69655eadddaa1295 Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Thu, 25 Feb 2021 09:28:58 +0100 Subject: [PATCH] ARKode: single step solve option --- include/deal.II/sundials/arkode.h | 70 +++++++--- source/sundials/arkode.cc | 102 ++++++++++---- tests/sundials/arkode_repeated_solve.cc | 129 ++++++++++++++++++ tests/sundials/arkode_repeated_solve.prm | 13 ++ ...eated_solve.with_sundials.geq.5.4.0.output | 121 ++++++++++++++++ 5 files changed, 392 insertions(+), 43 deletions(-) create mode 100644 tests/sundials/arkode_repeated_solve.cc create mode 100644 tests/sundials/arkode_repeated_solve.prm create mode 100644 tests/sundials/arkode_repeated_solve.with_sundials.geq.5.4.0.output diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index a6b5918312..27a91311ca 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -41,6 +41,8 @@ # ifdef DEAL_II_WITH_MPI # include # endif +# include + # include # include @@ -535,10 +537,44 @@ namespace SUNDIALS /** * Integrate the initial value problem. This function returns the final * number of computed steps. + * + * @param solution On input, this vector contains the initial condition. On + * output, it contains the solution at the final time. */ unsigned int solve_ode(VectorType &solution); + /** + * Integrate the initial value problem. Compared to the function above, this + * function allows to specify an @p intermediate_time for the next solution. + * Repeated calls of this function must use monotonously increasing values + * for @p intermediate_time. The last solution state is saved internally + * along with the @p intermediate_time and will be reused as initial + * condition for the next call. + * + * Users may find this function useful when integrating ARKode into an outer + * time loop of their own, especially when output_step() is too restrictive. + * + * @note @p intermediate_time may be larger than AdditionalData::final_time, + * which is ignored by this function. + * + * @param solution The final solution. If the solver restarts, either + * because it is the first ever solve or the flag @p reset_solver is + * set, the vector is also used as initial condition. + * @param intermediate_time The time for the incremental solution step. Must + * be greater than the last time that was used in a previous call to this + * function. + * @param reset_solver Optional flag to recreate all internal objects which + * may be desirable for spatial adaptivity methods. If set to `true`, + * reset() is called before solving the ODE, which sets @p solution as + * initial condition. This will *not* reset the stored time from previous + * calls to this function. + */ + unsigned int + solve_ode_incrementally(VectorType & solution, + const double intermediate_time, + const bool reset_solver = false); + /** * Clear internal memory and start with clean objects. This function is * called when the simulation starts and when the user returns true to a @@ -549,12 +585,12 @@ namespace SUNDIALS * a different function to solver_should_restart() that performs all mesh * changes, transfers the solution to the new mesh, and returns true. * - * @param[in] t The new starting time - * @param[in] h The new starting time step - * @param[in,out] y The new initial solution + * @param t The new starting time + * @param h The new starting time step + * @param y The new initial solution */ void - reset(const double t, const double h, VectorType &y); + reset(const double t, const double h, const VectorType &y); /** * Provides user access to the internally used ARKODE memory. @@ -1224,6 +1260,14 @@ namespace SUNDIALS << "Please provide an implementation for the function \"" << arg1 << "\""); + /** + * Internal routine to call ARKode repeatedly. + */ + int + do_evolve_time(VectorType & solution, + dealii::DiscreteTime &time, + const bool do_reset); + # if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) /** @@ -1238,9 +1282,11 @@ namespace SUNDIALS /** * Set up the solver and preconditioner for a non-identity mass matrix in * the ARKODE memory object based on the user-specified functions. + * @param solution The solution vector which is used as a template to create + * new vectors. */ void - setup_mass_solver(); + setup_mass_solver(const VectorType &solution); # endif @@ -1262,16 +1308,6 @@ namespace SUNDIALS */ void *arkode_mem; - /** - * ARKode solution vector. - */ - internal::NVectorView yy; - - /** - * ARKode absolute tolerances vector. - */ - internal::NVectorView abs_tolls; - /** * MPI communicator. SUNDIALS solver runs happily in * parallel. Note that if the library is compiled without MPI @@ -1280,9 +1316,9 @@ namespace SUNDIALS MPI_Comm communicator; /** - * Memory pool of vectors. + * The final time in the last call to solve_ode(). */ - GrowingVectorMemory mem; + double last_end_time; # if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0) std::unique_ptr> linear_solver; diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index 31d970927c..f143a2d43b 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -226,7 +226,6 @@ namespace SUNDIALS Assert(user_data != nullptr, ExcInternalError()); ARKode &solver = *static_cast *>(user_data); - GrowingVectorMemory mem; auto *src_v = internal::unwrap_nvector_const(v); auto *src_y = internal::unwrap_nvector_const(y); @@ -385,6 +384,7 @@ namespace SUNDIALS , communicator(is_serial_vector::value ? MPI_COMM_SELF : Utilities::MPI::duplicate_communicator(mpi_comm)) + , last_end_time(data.initial_time) { set_functions_to_trigger_an_assert(); } @@ -418,21 +418,64 @@ namespace SUNDIALS DiscreteTime time(data.initial_time, data.final_time, data.initial_step_size); - reset(time.get_current_time(), time.get_next_step_size(), solution); - if (output_step) - output_step(time.get_current_time(), solution, time.get_step_number()); + return do_evolve_time(solution, time, /* force_solver_restart = */ true); + } + - while (time.is_at_end() == false) + + template + unsigned int + ARKode::solve_ode_incrementally(VectorType & solution, + const double intermediate_time, + const bool reset_solver) + { + AssertThrow( + intermediate_time > last_end_time, + dealii::ExcMessage( + "The requested intermediate time is smaller than the last requested " + "intermediate time.")); + + const bool do_reset = reset_solver || arkode_mem == nullptr; + DiscreteTime time(last_end_time, intermediate_time, data.initial_step_size); + return do_evolve_time(solution, time, do_reset); + } + + + + template + int + ARKode::do_evolve_time(VectorType & solution, + DiscreteTime &time, + const bool do_reset) + { + if (do_reset) + { + reset(time.get_current_time(), time.get_next_step_size(), solution); + if (output_step) + output_step(time.get_current_time(), + solution, + time.get_step_number()); + } + + auto solution_nvector = internal::make_nvector_view(solution); + + while (!time.is_at_end()) { time.set_desired_next_step_size(data.output_period); 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); + const auto status = SundialsARKode(arkode_mem, + time.get_next_time(), + solution_nvector, + &actual_next_time, + ARK_NORMAL); # else - const auto status = ARKStepEvolve( - arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL); + const auto status = ARKStepEvolve(arkode_mem, + time.get_next_time(), + solution_nvector, + &actual_next_time, + ARK_NORMAL); # endif (void)status; AssertARKode(status); @@ -450,17 +493,20 @@ namespace SUNDIALS solution, time.get_step_number()); } - + last_end_time = time.get_current_time(); return time.get_step_number(); } + + # if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0) template void - ARKode::reset(const double current_time, - const double current_time_step, - VectorType & solution) + ARKode::reset(const double current_time, + const double current_time_step, + const VectorType &solution) { + last_end_time = current_time; if (arkode_mem) ARKodeFree(&arkode_mem); @@ -474,21 +520,21 @@ namespace SUNDIALS // just a view on the memory in solution, all write operations on yy by // ARKODE will automatically be mirrored to solution - yy = internal::make_nvector_view(solution); + auto initial_condition_nvector = internal::make_nvector_view(solution); status = ARKodeInit( arkode_mem, explicit_function ? &t_arkode_explicit_function : nullptr, implicit_function ? &t_arkode_implicit_function : nullptr, current_time, - yy); + initial_condition_nvector); AssertARKode(status); if (get_local_tolerances) { - abs_tolls = make_nvector_view(get_local_tolerances()); + const auto abs_tols = make_nvector_view(get_local_tolerances()); status = - ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls); + ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols); AssertARKode(status); } else @@ -570,8 +616,9 @@ namespace SUNDIALS void ARKode::reset(const double current_time, const double current_time_step, - VectorType &solution) + const VectorType &solution) { + last_end_time = current_time; if (arkode_mem) ARKStepFree(&arkode_mem); @@ -580,7 +627,7 @@ namespace SUNDIALS // just a view on the memory in solution, all write operations on yy by // ARKODE will automatically be mirrored to solution - yy = internal::make_nvector_view(solution); + auto initial_condition_nvector = internal::make_nvector_view(solution); Assert(explicit_function || implicit_function, ExcFunctionNotProvided("explicit_function || implicit_function")); @@ -589,15 +636,16 @@ namespace SUNDIALS explicit_function ? &t_arkode_explicit_function : nullptr, implicit_function ? &t_arkode_implicit_function : nullptr, current_time, - yy); + initial_condition_nvector); Assert(arkode_mem != nullptr, ExcInternalError()); if (get_local_tolerances) { - abs_tolls = internal::make_nvector_view(get_local_tolerances()); + const auto abs_tols = + internal::make_nvector_view(get_local_tolerances()); status = - ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls); + ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols); AssertARKode(status); } else @@ -614,7 +662,7 @@ namespace SUNDIALS setup_system_solver(solution); - setup_mass_solver(); + setup_mass_solver(solution); status = ARKStepSetInitStep(arkode_mem, current_time_step); AssertARKode(status); @@ -657,8 +705,9 @@ namespace SUNDIALS { // use default solver from SUNDIALS // TODO give user options + auto y_template = internal::make_nvector_view(solution); sun_linear_solver = - SUNLinSol_SPGMR(yy, + SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/); } @@ -709,7 +758,7 @@ namespace SUNDIALS template void - ARKode::setup_mass_solver() + ARKode::setup_mass_solver(const VectorType &solution) { int status; (void)status; @@ -726,8 +775,9 @@ namespace SUNDIALS } else { + auto y_template = internal::make_nvector_view(solution); sun_mass_linear_solver = - SUNLinSol_SPGMR(yy, + SUNLinSol_SPGMR(y_template, PREC_NONE, 0 /*krylov subvectors, 0 uses default*/); } diff --git a/tests/sundials/arkode_repeated_solve.cc b/tests/sundials/arkode_repeated_solve.cc new file mode 100644 index 0000000000..ad4a62939f --- /dev/null +++ b/tests/sundials/arkode_repeated_solve.cc @@ -0,0 +1,129 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +//----------------------------------------------------------- + +#include + +#include +#include + +#include + +#include "../tests.h" + +// Test repeated calls to solve_ode with harmonic oscillator problem + +using VectorType = Vector; + +double kappa = 1.0; + +std::unique_ptr> +create_solver() +{ + ParameterHandler prm; + SUNDIALS::ARKode::AdditionalData data; + data.add_parameters(prm); + + std::ifstream ifile(SOURCE_DIR "/arkode_repeated_solve.prm"); + prm.parse_input(ifile); + + auto ode = std::make_unique>(data); + + // will yield analytic solution y[0] = sin(kappa*t); y[1] = kappa*cos(kappa*t) + ode->explicit_function = + [&](double, const VectorType &y, VectorType &ydot) -> int { + ydot[0] = y[1]; + ydot[1] = -kappa * kappa * y[0]; + 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; + }; + return ode; +} + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, numbers::invalid_unsigned_int); + + auto ode = create_solver(); + { + Vector y(2); + y[0] = 0; + y[1] = kappa; + + // solve the ode in one go + ode->solve_ode(y); + + deallog << "First solve done" << std::endl; + + // solve again -> will call reset automatically + y[0] = 0; + y[1] = kappa; + ode->solve_ode(y); + deallog << "Second solve done" << std::endl; + + // Solve for a time greater than the originally requested final time + ode->solve_ode_incrementally(y, 3.0); + deallog << "Solve further than final time done" << std::endl; + } + + // solve in two steps with new object + ode = create_solver(); + + { + Vector y(2); + y[0] = 0; + y[1] = kappa; + + ode->solve_ode_incrementally(y, 1.0); + deallog << "Split solve - part 1 done" << std::endl; + } + + { + // N.B. initial conditions are not necessary as the internal state is reused + Vector y(2); + + ode->solve_ode_incrementally(y, 2.0); + deallog << "Split solve - part 2 done" << std::endl; + } + + { + Vector y(2); + y[0] = 0; + y[1] = kappa; + // solving with single-argument call reruns from start + ode->solve_ode(y); + deallog << "Rerun with single argument call done" << std::endl; + } + + { + Vector y0(2), y(2); + y0[0] = 0; + y0[1] = kappa; + // manually reset solver + ode->reset(0.0, 0.01, y0); + ode->solve_ode_incrementally(y, 2.0); + } + + return 0; +} diff --git a/tests/sundials/arkode_repeated_solve.prm b/tests/sundials/arkode_repeated_solve.prm new file mode 100644 index 0000000000..2bdb712fba --- /dev/null +++ b/tests/sundials/arkode_repeated_solve.prm @@ -0,0 +1,13 @@ +set Final time = 2.0 +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 Initial step size = 0.010000 + set Maximum number of nonlinear iterations = 10 + set Maximum order of ARK = 5 + set Minimum step size = 0.000001 +end diff --git a/tests/sundials/arkode_repeated_solve.with_sundials.geq.5.4.0.output b/tests/sundials/arkode_repeated_solve.with_sundials.geq.5.4.0.output new file mode 100644 index 0000000000..3f6116f868 --- /dev/null +++ b/tests/sundials/arkode_repeated_solve.with_sundials.geq.5.4.0.output @@ -0,0 +1,121 @@ + +DEAL::0.00000 0.00000 1.00000 +DEAL::0.100000 0.0998334 0.995004 +DEAL::0.200000 0.198669 0.980067 +DEAL::0.300000 0.295521 0.955336 +DEAL::0.400000 0.389420 0.921060 +DEAL::0.500000 0.479426 0.877583 +DEAL::0.600000 0.564643 0.825335 +DEAL::0.700000 0.644219 0.764841 +DEAL::0.800000 0.717356 0.696707 +DEAL::0.900000 0.783330 0.621606 +DEAL::1.00000 0.841478 0.540292 +DEAL::1.10000 0.891214 0.453587 +DEAL::1.20000 0.932042 0.362356 +DEAL::1.30000 0.963560 0.267498 +DEAL::1.40000 0.985452 0.169958 +DEAL::1.50000 0.997498 0.0707265 +DEAL::1.60000 0.999576 -0.0292043 +DEAL::1.70000 0.991667 -0.128845 +DEAL::1.80000 0.973849 -0.227205 +DEAL::1.90000 0.946302 -0.323292 +DEAL::2.00000 0.909300 -0.416148 +DEAL::First solve done +DEAL::0.00000 0.00000 1.00000 +DEAL::0.100000 0.0998334 0.995004 +DEAL::0.200000 0.198669 0.980067 +DEAL::0.300000 0.295521 0.955336 +DEAL::0.400000 0.389420 0.921060 +DEAL::0.500000 0.479426 0.877583 +DEAL::0.600000 0.564643 0.825335 +DEAL::0.700000 0.644219 0.764841 +DEAL::0.800000 0.717356 0.696707 +DEAL::0.900000 0.783330 0.621606 +DEAL::1.00000 0.841478 0.540292 +DEAL::1.10000 0.891214 0.453587 +DEAL::1.20000 0.932042 0.362356 +DEAL::1.30000 0.963560 0.267498 +DEAL::1.40000 0.985452 0.169958 +DEAL::1.50000 0.997498 0.0707265 +DEAL::1.60000 0.999576 -0.0292043 +DEAL::1.70000 0.991667 -0.128845 +DEAL::1.80000 0.973849 -0.227205 +DEAL::1.90000 0.946302 -0.323292 +DEAL::2.00000 0.909300 -0.416148 +DEAL::Second solve done +DEAL::2.10000 0.863210 -0.504849 +DEAL::2.20000 0.808497 -0.588504 +DEAL::2.30000 0.745707 -0.666278 +DEAL::2.40000 0.675463 -0.737397 +DEAL::2.50000 0.598472 -0.801147 +DEAL::2.60000 0.515502 -0.856891 +DEAL::2.70000 0.427379 -0.904075 +DEAL::2.80000 0.334987 -0.942226 +DEAL::2.90000 0.239250 -0.970961 +DEAL::3.00000 0.141118 -0.989995 +DEAL::Solve further than final time done +DEAL::0.00000 0.00000 1.00000 +DEAL::0.100000 0.0998334 0.995004 +DEAL::0.200000 0.198669 0.980067 +DEAL::0.300000 0.295521 0.955336 +DEAL::0.400000 0.389420 0.921060 +DEAL::0.500000 0.479426 0.877583 +DEAL::0.600000 0.564643 0.825335 +DEAL::0.700000 0.644219 0.764841 +DEAL::0.800000 0.717356 0.696707 +DEAL::0.900000 0.783330 0.621606 +DEAL::1.00000 0.841478 0.540292 +DEAL::Split solve - part 1 done +DEAL::1.10000 0.891214 0.453587 +DEAL::1.20000 0.932042 0.362356 +DEAL::1.30000 0.963560 0.267498 +DEAL::1.40000 0.985452 0.169958 +DEAL::1.50000 0.997498 0.0707265 +DEAL::1.60000 0.999576 -0.0292043 +DEAL::1.70000 0.991667 -0.128845 +DEAL::1.80000 0.973849 -0.227205 +DEAL::1.90000 0.946302 -0.323292 +DEAL::2.00000 0.909300 -0.416148 +DEAL::Split solve - part 2 done +DEAL::0.00000 0.00000 1.00000 +DEAL::0.100000 0.0998334 0.995004 +DEAL::0.200000 0.198669 0.980067 +DEAL::0.300000 0.295521 0.955336 +DEAL::0.400000 0.389420 0.921060 +DEAL::0.500000 0.479426 0.877583 +DEAL::0.600000 0.564643 0.825335 +DEAL::0.700000 0.644219 0.764841 +DEAL::0.800000 0.717356 0.696707 +DEAL::0.900000 0.783330 0.621606 +DEAL::1.00000 0.841478 0.540292 +DEAL::1.10000 0.891214 0.453587 +DEAL::1.20000 0.932042 0.362356 +DEAL::1.30000 0.963560 0.267498 +DEAL::1.40000 0.985452 0.169958 +DEAL::1.50000 0.997498 0.0707265 +DEAL::1.60000 0.999576 -0.0292043 +DEAL::1.70000 0.991667 -0.128845 +DEAL::1.80000 0.973849 -0.227205 +DEAL::1.90000 0.946302 -0.323292 +DEAL::2.00000 0.909300 -0.416148 +DEAL::Rerun with single argument call done +DEAL::0.100000 0.0998334 0.995004 +DEAL::0.200000 0.198669 0.980067 +DEAL::0.300000 0.295521 0.955336 +DEAL::0.400000 0.389420 0.921060 +DEAL::0.500000 0.479426 0.877583 +DEAL::0.600000 0.564643 0.825335 +DEAL::0.700000 0.644219 0.764841 +DEAL::0.800000 0.717356 0.696707 +DEAL::0.900000 0.783330 0.621606 +DEAL::1.00000 0.841478 0.540292 +DEAL::1.10000 0.891214 0.453587 +DEAL::1.20000 0.932042 0.362356 +DEAL::1.30000 0.963560 0.267498 +DEAL::1.40000 0.985452 0.169958 +DEAL::1.50000 0.997498 0.0707265 +DEAL::1.60000 0.999576 -0.0292043 +DEAL::1.70000 0.991667 -0.128845 +DEAL::1.80000 0.973849 -0.227205 +DEAL::1.90000 0.946302 -0.323292 +DEAL::2.00000 0.909300 -0.416148 -- 2.39.5