When using adaptive timestep, don't force fixed output times by default

When running kinetic electron simulations, it can cause problems to take
a very short ion timestep. When writing outputs at exactly fixed output
times, this can happen if the previous timestep happened to end just
before the output time. To avoid the very short step default to writing
output at whatever time the end of the timestep is that exceeds the set
output time. There is an option to force the previous behaviour of a
decreased timestep so that output is written exactly at the nominal
output time.

moment_kinetics/src/input_structs.jl

@@ -57,6 +57,7 @@ struct time_info{Terrorsum <: Real, T_debug_output, T_electron, Trkimp, Timpzero
     limit_caused_by::Vector{mk_int}
     nwrite_moments::mk_int
     nwrite_dfns::mk_int
+    exact_output_times::Bool
     moments_output_times::Vector{mk_float}
     dfns_output_times::Vector{mk_float}
     type::String

moment_kinetics/src/moment_kinetics_input.jl

@@ -146,6 +146,7 @@ function mk_input(input_dict=OptionsDict(); save_inputs_to_txt=false, ignore_MPI
         CFL_prefactor=-1.0,
         nwrite=1,
         nwrite_dfns=-1,
+        exact_output_times=false,
         type="SSPRK4",
         split_operators=false,
         steady_state_residual=false,

moment_kinetics/src/runge_kutta.jl

@@ -1166,10 +1166,25 @@ function adaptive_timestep_update_t_params!(t_params, CFL_limits, error_norms,
             t_params.failure_caused_by[end] += 1
         end
 
-        # If we were trying to take a step to the output timestep, dt will be smaller on
-        # the re-try, so will not reach the output time.
-        t_params.step_to_moments_output[] = false
-        t_params.step_to_dfns_output[] = false
+        if t_params.exact_output_times
+            # If we were trying to take a step to the output timestep, dt will be smaller on
+            # the re-try, so will not reach the output time.
+            t_params.step_to_moments_output[] = false
+            t_params.step_to_dfns_output[] = false
+        else
+            # If with the reduced dt the step will not pass the next output time,
+            # deactivate step_to_*_output[].
+            if (t_params.step_to_moments_output[]
+                && t_params.t[] + t_params.previous_dt[] + t_params.dt[] <
+                   t_params.moments_output_times[t_params.moments_output_counter[]])
+                t_params.step_to_moments_output[] = false
+            end
+            if (t_params.step_to_dfns_output[]
+                && t_params.t[] + t_params.previous_dt[] + t_params.dt[] <
+                   t_params.dfns_output_times[t_params.dfns_output_counter[]])
+                t_params.step_to_dfns_output[] = false
+            end
+        end
     elseif (error_norm[] > 1.0 || isnan(error_norm[])) && t_params.dt[] > t_params.minimum_dt * (1.0 + 1.0e-13)
         # (1.0 + 1.0e-13) fudge factor accounts for possible rounding errors when
         # t+dt=next_output_time.
@@ -1199,10 +1214,25 @@ function adaptive_timestep_update_t_params!(t_params, CFL_limits, error_norms,
             t_params.failure_caused_by[max_error_variable_index] += 1
         end
 
-        # If we were trying to take a step to the output timestep, dt will be smaller on
-        # the re-try, so will not reach the output time.
-        t_params.step_to_moments_output[] = false
-        t_params.step_to_dfns_output[] = false
+        if t_params.exact_output_times
+            # If we were trying to take a step to the output timestep, dt will be smaller on
+            # the re-try, so will not reach the output time.
+            t_params.step_to_moments_output[] = false
+            t_params.step_to_dfns_output[] = false
+        else
+            # If with the reduced dt the step will not pass the next output time,
+            # deactivate step_to_*_output[].
+            if (t_params.step_to_moments_output[]
+                && t_params.t[] + t_params.previous_dt[] + t_params.dt[] <
+                   t_params.moments_output_times[t_params.moments_output_counter[]])
+                t_params.step_to_moments_output[] = false
+            end
+            if (t_params.step_to_dfns_output[]
+                && t_params.t[] + t_params.previous_dt[] + t_params.dt[] <
+                   t_params.dfns_output_times[t_params.dfns_output_counter[]])
+                t_params.step_to_dfns_output[] = false
+            end
+        end
 
         #println("t=$t, timestep failed, error_norm=$(error_norm[]), error_norms=$error_norms, decreasing timestep to ", t_params.dt[])
     else
@@ -1211,9 +1241,11 @@ function adaptive_timestep_update_t_params!(t_params, CFL_limits, error_norms,
         t_params.previous_dt[] = t_params.dt[]
 
         if t_params.step_to_moments_output[] || t_params.step_to_dfns_output[]
-            # Completed an output step, reset dt to what it was before it was reduced to reach
-            # the output time
-            t_params.dt[] = t_params.dt_before_output[]
+            if !t_params.exact_output_times
+                # Completed an output step, reset dt to what it was before it was reduced to reach
+                # the output time
+                t_params.dt[] = t_params.dt_before_output[]
+            end
 
             if t_params.step_to_moments_output[]
                 t_params.step_to_moments_output[] = false
@@ -1227,7 +1259,8 @@ function adaptive_timestep_update_t_params!(t_params, CFL_limits, error_norms,
             if t_params.dt[] > CFL_limit[]
                 t_params.dt[] = CFL_limit[]
             end
-        else
+        end
+        if !t_params.exact_output_times || !(t_params.write_moments_output[] || t_params.write_dfns_output[])
             # Adjust timestep according to Fehlberg's suggestion
             # (https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method).
             # `step_update_prefactor` is a constant numerical factor to make the estimate
@@ -1337,7 +1370,9 @@ function adaptive_timestep_update_t_params!(t_params, CFL_limits, error_norms,
         && (current_time + t_params.dt[] >= t_params.moments_output_times[t_params.moments_output_counter[]]))
 
         t_params.dt_before_output[] = current_dt
-        t_params.dt[] = t_params.moments_output_times[t_params.moments_output_counter[]] - current_time
+        if t_params.exact_output_times
+            t_params.dt[] = t_params.moments_output_times[t_params.moments_output_counter[]] - current_time
+        end
         t_params.step_to_moments_output[] = true
 
         if t_params.dt[] < 0.0
@@ -1352,7 +1387,9 @@ function adaptive_timestep_update_t_params!(t_params, CFL_limits, error_norms,
         && (current_time + t_params.dt[] >= t_params.dfns_output_times[t_params.dfns_output_counter[]]))
 
         t_params.dt_before_output[] = current_dt
-        t_params.dt[] = t_params.dfns_output_times[t_params.dfns_output_counter[]] - current_time
+        if t_params.exact_output_times
+            t_params.dt[] = t_params.dfns_output_times[t_params.dfns_output_counter[]] - current_time
+        end
         t_params.step_to_dfns_output[] = true
 
         if t_params.dt[] < 0.0

moment_kinetics/src/time_advance.jl

@@ -355,7 +355,7 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload,
 
     end_time = mk_float(code_time + t_input["dt"] * t_input["nstep"])
     epsilon = 1.e-11
-    if adaptive || t_input["write_after_fixed_step_count"]
+    if adaptive && !t_input["write_after_fixed_step_count"]
         if t_input["nwrite"] == 0
             moments_output_times = [end_time]
         else
@@ -483,11 +483,12 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload,
                      step_to_moments_output, step_to_dfns_output, write_moments_output,
                      write_dfns_output, Ref(0), Ref(0), Ref{mk_float}(0.0), Ref(0),
                      Ref(0), Ref(0), mk_int[], mk_int[], t_input["nwrite"],
-                     t_input["nwrite_dfns"], moments_output_times, dfns_output_times,
-                     t_input["type"], rk_coefs, rk_coefs_implicit,
-                     implicit_coefficient_is_zero, n_rk_stages, rk_order, adaptive,
-                     low_storage, mk_float(t_input["rtol"]), mk_float(t_input["atol"]),
-                     mk_float(t_input["atol_upar"]),
+                     t_input["nwrite_dfns"],
+                     electron !== nothing && t_input["exact_output_times"],
+                     moments_output_times, dfns_output_times, t_input["type"], rk_coefs,
+                     rk_coefs_implicit, implicit_coefficient_is_zero, n_rk_stages,
+                     rk_order, adaptive, low_storage, mk_float(t_input["rtol"]),
+                     mk_float(t_input["atol"]), mk_float(t_input["atol_upar"]),
                      mk_float(t_input["step_update_prefactor"]),
                      mk_float(t_input["max_increase_factor"]),
                      mk_float(t_input["max_increase_factor_near_last_fail"]),
@@ -1878,9 +1879,21 @@ function  time_advance!(pdf, scratch, scratch_implicit, scratch_electron, t_para
             end
             if write_moments
                 t_params.moments_output_counter[] += 1
+                if !t_params.exact_output_times
+                    while (t_params.moments_output_counter[] ≤ length(t_params.moments_output_times)
+                           && t_params.moments_output_times[t_params.moments_output_counter[]] ≤ t_params.t[])
+                        t_params.moments_output_counter[] += 1
+                    end
+                end
             end
             if write_dfns
                 t_params.dfns_output_counter[] += 1
+                if !t_params.exact_output_times
+                    while (t_params.dfns_output_counter[] ≤ length(t_params.dfns_output_times)
+                           && t_params.dfns_output_times[t_params.dfns_output_counter[]] ≤ t_params.t[])
+                        t_params.dfns_output_counter[] += 1
+                    end
+                end
             end
 
             if write_moments || write_dfns || finish_now

moment_kinetics/test/braginskii_electrons_imex_tests.jl

@@ -68,6 +68,7 @@ test_input = OptionsDict( "composition" => OptionsDict("n_ion_species" => 1,
                                                      "minimum_dt" => 1.e-7,
                                                      "rtol" => 1.0e-7,
                                                      "nwrite" => 10000,
+                                                     "exact_output_times" => true,
                                                      "high_precision_error_sum" => true),
                   "nonlinear_solver" => OptionsDict("nonlinear_max_iterations" => 100),
                   "r" => OptionsDict("ngrid" => 1,

moment_kinetics/test/recycling_fraction_tests.jl

@@ -169,6 +169,11 @@ test_input_adaptive_split3["timestepping"] = recursive_merge(test_input_adaptive
                                                                          "minimum_dt" => 1.0e-7,
                                                                          "step_update_prefactor" => 0.064))
 
+# Test exact_output_times option in full-f/split1/split2 cases
+test_input_adaptive["timestepping"]["exact_output_times"] = true
+test_input_adaptive_split1["timestepping"]["exact_output_times"] = true
+test_input_adaptive_split2["timestepping"]["exact_output_times"] = true
+
 """
 Run a test for a single set of parameters
 """
@@ -341,14 +346,14 @@ function runtests()
         @testset "Adaptive timestep - split 3" begin
             test_input_adaptive_split3["output"]["base_directory"] = test_output_directory
             run_test(test_input_adaptive_split3,
-                     [-0.034623352735472034, -0.03200541773193755, -0.02714032291656093,
-                      -0.020924986472905527, -0.01015057042512689, 0.0027893133203071574,
-                      0.012837899470698978, 0.022096372980618853, 0.0330348469665054,
-                      0.041531828755231016, 0.045382106043818246, 0.046246244563868354,
-                      0.042551970615727366, 0.034815169767529956, 0.027080688565416164,
-                      0.017886490800418996, 0.004784403555306537, -0.007762152788142663,
-                      -0.01629330539573498, -0.02413421820486561, -0.0315621379076817,
-                      -0.03416924694766477], rtol=6.0e-4, atol=2.0e-12)
+                     [-0.0346196925024167, -0.03200201693849987, -0.02713764319615098,
+                      -0.02092311349672712, -0.010150026206894121, 0.0027883420935253572,
+                      0.012835791449600767, 0.02209326318113659, 0.03303078703903627,
+                      0.04152829640863164, 0.04538051487359227, 0.04624543438581702,
+                      0.04254876799453081, 0.03481104153755928, 0.027077084096581314,
+                      0.01788382934269672, 0.00478320487966262, -0.0077618876322877485,
+                      -0.016292009420807548, -0.024131976958124225, -0.031559093785483404,
+                      -0.0341657304695615], rtol=6.0e-4, atol=2.0e-12)
         end
 
         @long @testset "Check other timestep - $type" for