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