From c6027b6f3a277d9ccb06f43c9bd205e748abe3e1 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Nov 2024 22:53:38 +0000 Subject: [PATCH] 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 | 1 + moment_kinetics/src/moment_kinetics_input.jl | 1 + moment_kinetics/src/runge_kutta.jl | 65 +++++++++++++++---- moment_kinetics/src/time_advance.jl | 25 +++++-- .../test/braginskii_electrons_imex_tests.jl | 1 + .../test/recycling_fraction_tests.jl | 21 +++--- 6 files changed, 86 insertions(+), 28 deletions(-) diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 682830365..dfc80b3f4 100644 --- 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 index e1008bf0e..c6ff98c85 100644 --- 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 index 0c87369b7..772cda8ba 100644 --- 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 index 8f2d7095b..d2c0e295a 100644 --- 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/braginskii_electrons_imex_tests.jl b/moment_kinetics/test/braginskii_electrons_imex_tests.jl index 1104271f3..fb91dfc3c 100644 --- a/moment_kinetics/test/braginskii_electrons_imex_tests.jl +++ b/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, diff --git a/moment_kinetics/test/recycling_fraction_tests.jl b/moment_kinetics/test/recycling_fraction_tests.jl index 04d289637..77b14fc2b 100644 --- a/moment_kinetics/test/recycling_fraction_tests.jl +++ b/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