From: Reza Rastak Date: Mon, 27 Apr 2020 22:28:22 +0000 (-0700) Subject: arkode ode solver internally uses DiscreteTime X-Git-Tag: v9.3.0-rc1~1224^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9973%2Fhead;p=dealii.git arkode ode solver internally uses DiscreteTime --- diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h index 648411dd41..c327ff5249 100644 --- a/include/deal.II/sundials/arkode.h +++ b/include/deal.II/sundials/arkode.h @@ -325,7 +325,7 @@ namespace SUNDIALS * @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 + * @param output_period Desired time interval between each output * * Running parameters: * @@ -475,7 +475,8 @@ namespace SUNDIALS unsigned int maximum_order; /** - * Time period between each output. + * Desired time period between each output. The actual output time period + * may be adjusted by Arkode. */ double output_period; diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc index cec5de071f..64be7b34b2 100644 --- a/source/sundials/arkode.cc +++ b/source/sundials/arkode.cc @@ -20,6 +20,7 @@ #ifdef DEAL_II_WITH_SUNDIALS +# include # include # include @@ -270,18 +271,9 @@ namespace SUNDIALS unsigned int ARKode::solve_ode(VectorType &solution) { - unsigned int system_size = solution.size(); + const unsigned int system_size = solution.size(); - double t = data.initial_time; - double h = data.initial_step_size; - unsigned int step_number = 0; - - int status; - (void)status; - - // The solution is stored in - // solution. Here we take only a - // view of it. + // The solution is stored in solution. Here we take only a view of it. # ifdef DEAL_II_WITH_MPI if (is_serial_vector::value == false) { @@ -302,33 +294,38 @@ namespace SUNDIALS yy = N_VNew_Serial(system_size); abs_tolls = N_VNew_Serial(system_size); } - reset(data.initial_time, data.initial_step_size, solution); - double next_time = data.initial_time; + 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(0, solution, 0); + output_step(time.get_current_time(), solution, time.get_step_number()); - while (t < data.final_time) + while (time.is_at_end() == false) { - next_time += data.output_period; - - status = SundialsARKode(arkode_mem, next_time, yy, &t, ARK_NORMAL); - - AssertARKode(status); - - status = ARKodeGetLastStep(arkode_mem, &h); + time.set_desired_next_step_size(data.output_period); + double actual_next_time; + const auto status = SundialsARKode( + arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL); + (void)status; AssertARKode(status); copy(solution, yy); - while (solver_should_restart(t, solution)) - reset(t, h, solution); + time.set_next_step_size(actual_next_time - time.get_current_time()); + time.advance_time(); - step_number++; + while (solver_should_restart(time.get_current_time(), solution)) + reset(time.get_current_time(), + time.get_previous_step_size(), + solution); if (output_step) - output_step(t, solution, step_number); + output_step(time.get_current_time(), + solution, + time.get_step_number()); } // Free the vectors which are no longer used. @@ -345,7 +342,7 @@ namespace SUNDIALS N_VDestroy_Serial(abs_tolls); } - return step_number; + return time.get_step_number(); } template diff --git a/tests/sundials/harmonic_oscillator_05.output b/tests/sundials/harmonic_oscillator_05.output index 507044c3b2..04f8401db1 100644 --- a/tests/sundials/harmonic_oscillator_05.output +++ b/tests/sundials/harmonic_oscillator_05.output @@ -99,4 +99,3 @@ 9.8 0.991567 2.63469 2.49996 9.9 1.02585 2.61958 2.50009 10 1.06497 2.59596 2.49997 -10 1.06497 2.59596 2.49997 diff --git a/tests/sundials/harmonic_oscillator_06.output b/tests/sundials/harmonic_oscillator_06.output index 507044c3b2..04f8401db1 100644 --- a/tests/sundials/harmonic_oscillator_06.output +++ b/tests/sundials/harmonic_oscillator_06.output @@ -99,4 +99,3 @@ 9.8 0.991567 2.63469 2.49996 9.9 1.02585 2.61958 2.50009 10 1.06497 2.59596 2.49997 -10 1.06497 2.59596 2.49997