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.
mabarnes · Nov 28, 2024 · 869b08c · 869b08c
1 parent 1ff204e
commit 869b08c
Show file tree

Hide file tree

Showing 5 changed files with 102 additions and 47 deletions.
diff --git a/moment_kinetics/src/input_structs.jl b/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

diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/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,

diff --git a/moment_kinetics/src/runge_kutta.jl b/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

diff --git a/moment_kinetics/src/time_advance.jl b/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

diff --git a/moment_kinetics/test/recycling_fraction_tests.jl b/moment_kinetics/test/recycling_fraction_tests.jl
@@ -169,6 +169,9 @@ 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 split2 case
+test_input_adaptive_split2["timestepping"]["exact_output_times"] = true
+
 """
 Run a test for a single set of parameters
 """
@@ -296,17 +299,17 @@ function runtests()
                       -0.03231930055546638])
         end
 
-        fullf_expected_output = [-0.04372543535228032, -0.02233515082616229,
-                                 -0.012793688037658377, -0.010786492944264052,
-                                 -0.007051439902278702, -0.0001605908774545327,
-                                 0.005982619745890949, 0.0118094191749825,
-                                 0.01954207152061524, 0.02978202423468538,
-                                 0.039384279904624404, 0.042446003403153604,
-                                 0.03181914367119813, 0.021111423438351817,
-                                 0.015103049638495273, 0.009135485828230407,
-                                 0.0010369322036392606, -0.005949066066045502,
-                                 -0.00942148866222427, -0.011607485576226423,
-                                 -0.020871221194795328, -0.03762871759968933]
+        fullf_expected_output = [-0.04374608549468029, -0.02233533054697219,
+                                 -0.012765638772788954, -0.010755864687704396,
+                                 -0.0070362685354843696, -0.00016505331520153655,
+                                 0.0059678980942991025, 0.011786668520847633,
+                                 0.019512193511411577, 0.029755929374471413,
+                                 0.039368696435893843, 0.04243414264798086,
+                                 0.03179501689779317, 0.02108115144969106,
+                                 0.015076445791375681, 0.00911630168608071,
+                                 0.0010301551995991337, -0.0059379253621724346,
+                                 -0.009396663595722667, -0.011574763104943042,
+                                 -0.02087042396890708, -0.03764968016116834]
         @testset "Adaptive timestep - full-f" begin
             test_input_adaptive["output"]["base_directory"] = test_output_directory
             run_test(test_input_adaptive,
@@ -315,14 +318,14 @@ function runtests()
         @testset "Adaptive timestep - split 1" begin
             test_input_adaptive_split1["output"]["base_directory"] = test_output_directory
             run_test(test_input_adaptive_split1,
-                     [-0.04375862714017892, -0.022363510973059945, -0.012739964397542611,
-                      -0.010806509398868007, -0.007052551067569563,
-                      -0.0001618866835357178, 0.005980921838191561, 0.011808361372364367,
-                      0.019540868336503224, 0.02978014755372564, 0.03938085813395519,
-                      0.042446888380863836, 0.031821059258512106, 0.021109010112552534,
-                      0.015101702015235266, 0.009134407186439548, 0.0010347434646523774,
-                      -0.005951302261109976, -0.009412276056941643, -0.011636393512121094,
-                      -0.020739923046188418, -0.03769486232955374], rtol=6.0e-4,
+                     [-0.043762149596224965, -0.022363280144398198, -0.01273491805603744,
+                      -0.010801103203158597, -0.00704987643976939,
+                      -0.00016267906218628162, 0.005978320579621679, 0.011804343061782241,
+                      0.019535591473459873, 0.02977553787462779, 0.03937810441243841,
+                      0.042444795477791036, 0.03181679728649564, 0.021103662728488865,
+                      0.015097003117537804, 0.009131018917786756, 0.0010335446855131265,
+                      -0.005949340355649067, -0.009407891192708668, -0.011630821571091444,
+                      -0.020739593081699037, -0.03769843453544615], rtol=6.0e-4,
                      atol=2.0e-12)
         end
         @testset "Adaptive timestep - split 2" begin
@@ -341,14 +344,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