From 2f61c49b9b45761724b7abcc78e32561eef96ed4 Mon Sep 17 00:00:00 2001 From: Spencer Lee Date: Sat, 30 Mar 2024 20:37:24 -0400 Subject: [PATCH] Made new convergence gatherer --- src/Tests/test_convergence.jl | 141 ++++++++++++---------------------- 1 file changed, 50 insertions(+), 91 deletions(-) diff --git a/src/Tests/test_convergence.jl b/src/Tests/test_convergence.jl index 4136404..c8edff2 100644 --- a/src/Tests/test_convergence.jl +++ b/src/Tests/test_convergence.jl @@ -1,5 +1,6 @@ # I should include a manufactured solution in this + """ Double the step size to get more precise solutions, and compare the change in error with the "true" solution (solution using the most steps / smallest step @@ -16,131 +17,89 @@ way more computation than necessary to get the 'true' solution. Richardson extrapolation seems like a good idea for that. """ -function get_history_convergence(prob, control, pcof, N_iterations; - orders=(2, 4, 6, 8), true_history=missing, error_limit=-Inf, n_runs=1, - kwargs... +function get_histories(prob, controls, pcof, N_iterations; orders=(2,4,6,8), + min_error_limit=-Inf, max_error_limit=Inf, base_nsteps=missing, + nsteps_change_factor=2, evolution_kwargs... ) - base_nsteps = prob.nsteps - nsteps_change_factor = 2 - - # Copy problem so we can mutate nsteps without altering input - prob_copy = copy(prob) - - histories = [] - errors_all = Matrix{Float64}(undef, N_iterations, length(orders)) - timing_all = Matrix{Float64}(undef, N_iterations, length(orders)) - timing_stddev_all = Matrix{Float64}(undef, N_iterations, length(orders)) - # Fill with NaN, so that we can break the loop early and still plot, ignoring unfilled values - errors_all .= NaN - timing_all .= NaN - timing_stddev_all .= NaN + prob_original_nsteps = prob.nsteps - step_sizes = (prob.tf/base_nsteps) ./ [nsteps_change_factor^k for k in 0:N_iterations-1] + prob_copy = copy(prob) - for (j, order) in enumerate(orders) - errors = Vector{Float64}(undef, N_iterations) - timing = Vector{Float64}(undef, N_iterations) - timing_stddev = Vector{Float64}(undef, N_iterations) - errors .= NaN - timing .= NaN - timing_stddev .= NaN + # Default to whatever nsteps the problem original had + if ismissing(base_nsteps) + base_nsteps = prob_original_nsteps + end + ret_vec = [] - N_derivatives = div(order, 2) + println("Beginning at ", Dates.now()) - println("========================================") - println("Running Order ", order) - println("========================================") - - histories_this_order = [] - - if ismissing(true_history) - println("----------------------------------------") - println("True history not given, using Richardson extrapolation\n") - println("Calculating solution with base_nsteps=", base_nsteps) - println("----------------------------------------") - prob_copy.nsteps = base_nsteps - base_history = eval_forward( - prob_copy, control, pcof, order=order; - saveEveryNsteps=div(prob_copy.nsteps, base_nsteps), - kwargs... - ) - # Only include state, not derivatives - base_history = base_history[:,1,:,:] - push!(histories_this_order, base_history) - end + for order in orders + println("#"^40, "\n") + println("Order ", order, "\n") + println("#"^40) + nsteps_vec = Int64[] + step_sizes = Float64[] + elapsed_times = Float64[] + histories = Array{Float64, 3}[] + richardson_errors = Float64[] for k in 1:N_iterations - nsteps_multiplier = nsteps_change_factor^k + println("Starting iteration ", k, " at ", Dates.now()) + nsteps_multiplier = nsteps_change_factor^(k-1) prob_copy.nsteps = base_nsteps*nsteps_multiplier - elapsed_times = zeros(n_runs) - - elapsed_times[1] = @elapsed history = eval_forward( - prob_copy, control, pcof, order=order; - saveEveryNsteps=div(prob_copy.nsteps, base_nsteps), - kwargs... - + # Run forward evolution + elapsed_time = @elapsed history = eval_forward( + prob_copy, controls, pcof, order=order; + saveEveryNsteps=nsteps_multiplier, + evolution_kwargs... ) - for rerun_i in 2:n_runs - elapsed_times[rerun_i] = @elapsed history = eval_forward( - prob_copy, control, pcof, order=order; - saveEveryNsteps=div(prob_copy.nsteps, base_nsteps), - kwargs... - ) - end - mean_elapsed_time = sum(elapsed_times)/length(elapsed_times) - stddev_elapsed_time = sum((elapsed_times .- mean_elapsed_time) .^ 2) - stddev_elapsed_time /= length(elapsed_times)-1 - stddev_elapsed_time = sqrt(stddev_elapsed_time) - # Only compare state, not derivatives history = history[:,1,:,:] - push!(histories_this_order, history) - if ismissing(true_history) - history_prev = histories_this_order[k] - error = richardson_extrap_rel_err(history, history_prev, order) - else - error = norm(history - true_history)/norm(true_history) + # Compute Richardson Error + richardson_err = NaN + if (k > 1) + history_prev = histories[k-1] + richardson_err = richardson_extrap_rel_err(history, history_prev, order) end - errors[k] = error - timing[k] = mean_elapsed_time - timing_stddev[k] = stddev_elapsed_time + push!(nsteps_vec, prob_copy.nsteps) + push!(step_sizes, prob_copy.tf / prob_copy.nsteps) + push!(elapsed_times, elapsed_time) + push!(histories, history) + push!(richardson_errors, richardson_err) - println("Nsteps = ", prob_copy.nsteps) - println("Error = ", error) - println("Mean Elapsed Time = ", mean_elapsed_time) - println("StdDev Elapsed Time = ", stddev_elapsed_time) + println("Finished iteration ", k, " at ", Dates.now()) + println("Nsteps = \t", prob_copy.nsteps) + println("Richardson Error = \t", richardson_err) + println("Elapsed Time = \t", elapsed_time) println("----------------------------------------") # Break once we reach high enough precision - if error < error_limit + if richardson_err < min_error_limit break end + # If we are reasonably precise, break if the error increases twice # (numerical saturation) if k > 2 - if (error < 1e-4) && (error > errors[k-1]) && (errors[k-1] > errors[k-2]) + if (error < max_error_limit) && (error > errors[k-1]) && (errors[k-1] > errors[k-2]) break end end end - - push!(histories, histories_this_order) - - errors_all[:,j] .= errors - timing_all[:,j] .= timing - timing_stddev_all[:,j] .= timing_stddev - #errors_all[:,length(orders)+j] .= order_line + push!(ret_vec, (order, nsteps_vec, step_sizes, elapsed_times, histories, richardson_errors)) end + + println("Finished at ", Dates.now()) + println("Returning (order, nsteps_vec, step_sizes, elapsed_times, histories, richardson_errors) for each order.") - #return step_sizes, errors_all, timing_all, timing_stddev_all - return step_sizes, errors_all, timing_all, timing_stddev_all, histories + return ret_vec end