From 7bf717ad12157a0ea853696dba714c1c838e28ed Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 4 Jun 2024 20:44:23 +0200 Subject: [PATCH] adjust the test and add a seed to the randomized initial guess for reproducibility --- Project.toml | 9 ++-- ext/TrixiNLsolveExt.jl | 11 ++++- .../methods_PERK3.jl | 30 ++++++------ test/test_structured_1d.jl | 8 +-- test/test_unit.jl | 49 ++++++++----------- 5 files changed, 55 insertions(+), 52 deletions(-) diff --git a/Project.toml b/Project.toml index 6ae18fb29a4..9abf4d36b5d 100644 --- a/Project.toml +++ b/Project.toml @@ -28,6 +28,7 @@ Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" @@ -50,14 +51,14 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [weakdeps] -Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Convex = "f65535da-76fb-5f13-bab9-19810c17039a" ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" [extensions] -TrixiMakieExt = "Makie" TrixiConvexECOSExt = ["Convex", "ECOS"] +TrixiMakieExt = "Makie" TrixiNLsolveExt = "NLsolve" [compat] @@ -112,7 +113,7 @@ UUIDs = "1.6" julia = "1.8" [extras] -Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Convex = "f65535da-76fb-5f13-bab9-19810c17039a" ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" -NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" \ No newline at end of file +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" diff --git a/ext/TrixiNLsolveExt.jl b/ext/TrixiNLsolveExt.jl index 4d2c1108ea7..cbe022b1f3a 100644 --- a/ext/TrixiNLsolveExt.jl +++ b/ext/TrixiNLsolveExt.jl @@ -9,6 +9,9 @@ else using ..NLsolve: nlsolve end +# Use other necessary libraries +using Random: seed! + # Use functions and additional symbols that are not exported using Trixi: Trixi, PairedExplicitRK3_butcher_tableau_objective_function, @muladd @@ -36,8 +39,13 @@ function Trixi.solve_a_unknown!(a_unknown, num_stages, monomial_coeffs, c_s2, c; c_s2) end + # Set the seed for reproducibility of the initial guess of a_unknown + seed!(5555) + while !is_sol_valid + # Initialize initial guess + # The nonlinear system may have multiple valid solutions, so a reproducible initial guess is important x0 = 0.1 .* rand(num_stages) x0[1] = 0.0 x0[2] = c[2] @@ -47,7 +55,8 @@ function Trixi.solve_a_unknown!(a_unknown, num_stages, monomial_coeffs, c_s2, c; a_unknown = sol.zero - # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) and subsequently c[i] - a[i, i-1] >= 0.0 + # Check if the values a[i, i-1] >= 0.0 (which stem from the nonlinear solver) + # and subsequently c[i] - a[i, i-1] >= 0.0 is_sol_valid = all(x -> !isnan(x) && x >= 0, a_unknown[3:end]) && all(x -> !isnan(x) && x >= 0, c[3:end] .- a_unknown[3:end]) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 6c6057f4bc3..ae2042872af 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -94,7 +94,8 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, num_stages) # Solve the nonlinear system of equations from monomial coefficient and - # Butcher array abscissae c to find Butcher matrix A + # Butcher array abscissae c to find Butcher matrix A + # This function is extended in TrixiNLsolveExt.jl a_unknown = solve_a_unknown!(a_unknown, num_stages, monomial_coeffs, cS2, c; verbose) end @@ -107,9 +108,8 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, return a_matrix, c, dt_opt end -#TODO: Correct this and make it actually import A matrix # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3 -# using provided monomial coefficients file +# using provided values of coefficients a in a matrix of Butcher tableau function compute_PairedExplicitRK3_butcher_tableau(num_stages, base_path_a_coeffs::AbstractString; cS2) @@ -118,26 +118,26 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, c = compute_c_coeffs_SSP33(num_stages, cS2) # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency) - coeffs_max = num_stages - 2 + a_coeffs_max = num_stages - 2 - a_matrix = zeros(coeffs_max, 2) + a_matrix = zeros(a_coeffs_max, 2) a_matrix[:, 1] = c[3:end] path_a_coeffs = joinpath(base_path_a_coeffs, - "a_" * string(num_stages) * "_" * string(num_stages) * ".txt") + "a_" * string(num_stages) * "_" * string(num_stages) * + ".txt") @assert isfile(path_a_coeffs) "Couldn't find file" - A = readdlm(path_a_coeffs, Float64) - num_a_coeffs = size(A, 1) + a_coeffs = readdlm(path_a_coeffs, Float64) + num_a_coeffs = size(a_coeffs, 1) - @assert num_a_coeffs == coeffs_max - a_matrix[:, 1] -= A - a_matrix[:, 2] = A + @assert num_a_coeffs == a_coeffs_max + a_matrix[:, 1] -= a_coeffs + a_matrix[:, 2] = a_coeffs return a_matrix, c end -#TODO: explain dt_opt and also explain base_path_a_coeffs in the first constructor @doc raw""" PairedExplicitRK3(num_stages, base_path_a_coeffs::AbstractString, dt_opt; @@ -149,9 +149,9 @@ end Parameters: - `num_stages` (`Int`): Number of stages in the PERK method. - - `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing - monomial coefficients of the stability polynomial of PERK method. - The coefficients should be stored in a text file at `joinpath(base_path_monomial_coeffs, "gamma_$(num_stages).txt")` and separated by line breaks. + - `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the matrix A in + the Butcher tableau of the Runge Kutta method. + The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages)_.$(num_stages)txt")` and separated by line breaks. - `tspan`: Time span of the simulation. - `semi` (`AbstractSemidiscretization`): Semidiscretization setup. - `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index 3401b314ceb..7b1c28aa1d6 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -47,7 +47,7 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_shockcapturing.jl"), l2=[0.08015029105233593], linf=[0.610709468736576], - atol=1.0e-5) + atol=1.0e-12) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -58,11 +58,11 @@ end end end - @trixi_testset "elixir_burgers_perk3.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_perk3.jl"), - l2=[8.120426329330083e-8], - linf=[4.906376900315479e-7]) + l2=[8.120426320528704e-8], + linf=[4.906376875890572e-7], + atol=1.0e-5) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_unit.jl b/test/test_unit.jl index eecbe315bc8..62a04496c07 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1679,18 +1679,9 @@ end end @testset "PERK Single p3 Constructors" begin - # We are only testing from the second row to the end of `ode_algorithm.a_matrix` due to the nature of the system of equations - # used to find the `a_matrix` in the Butcher tableau. This system can have multiple valid solutions, which can cause slight - # variations in the results of the tests. These variations don't affect the validity of the solutions, but they can cause - # the tests to fail due to small differences in the expected and actual values. - - # The first row of `ode_algorithm.a_matrix` is particularly susceptible to these variations, while the rest of the rows - # remain consistent across different runs. Therefore, to make the tests more robust and less prone to false negatives, - # we are excluding the first row from the test. This approach allows us to verify the correctness of the majority of the - # `a_matrix` while avoiding the issue of test instability caused by the multiple valid solutions of the system of equations. path_coeff_file = mktempdir() Trixi.download("https://gist.githubusercontent.com/warisa-r/0796db36abcd5abe735ac7eebf41b973/raw/32889062fd5dcf7f450748f4f5f0797c8155a18d/a_8_8.txt", - joinpath(path_coeff_file, "a_8_8.txt")) + joinpath(path_coeff_file, "a_8_8.txt")) # Value of dt_opt obtained from running the simulation in elixir_burgers_perk3 # The value plays no role in the result but added so that the constructor can be called @@ -1698,31 +1689,33 @@ end ode_algorithm = Trixi.PairedExplicitRK3(8, path_coeff_file, dt_opt) #TODO: adjust this value according to the result in the test pipeline - @test isapprox(ode_algorithm.a_matrix[2:end, :], - [0.496535 0.103465 - 0.649689 0.150311 - 0.789172 0.210828 - 0.752297 0.247703 - 0.311926 0.188074], atol = 1e-13) - - #TODO: make this a p3 Constructors + println(ode_algorithm.a_matrix) # Value in CI differs slightly from what I get locally + @test isapprox(ode_algorithm.a_matrix, + [0.3355167784195604 0.06448322158043965 + 0.4965349205803965 0.10346507941960345 + 0.6496890792935297 0.15031092070647037 + 0.789172498521197 0.21082750147880308 + 0.7522972036571336 0.2477027963428664 + 0.31192569908571666 0.18807430091428337], atol = 1e-13) + Trixi.download("https://gist.githubusercontent.com/warisa-r/8d93f6a3ae0635e13b9f51ee32ab7fff/raw/54dc5b14be9288e186b745facb5bbcb04d1476f8/EigenvalueList_Refined2.txt", joinpath(path_coeff_file, "spectrum.txt")) eig_vals = readdlm(joinpath(path_coeff_file, "spectrum.txt"), ComplexF64) tspan = (0.0, 1.0) - ode_algorithm = Trixi.PairedExplicitRK3(10, tspan, vec(eig_vals)) + ode_algorithm = Trixi.PairedExplicitRK3(13, tspan, vec(eig_vals)) #TODO: adjust this value according to the result in the test pipeline - display(ode_algorithm.a_matrix) # Value in CI differs slightly from what I get locally - @test isapprox(ode_algorithm.a_matrix[2:end, :], - [0.406023 0.0225489 - 0.534288 0.0371408 - 0.654943 0.0593431 - 0.76216 0.0949827 - 0.844659 0.155341 - 0.771998 0.228002 - 0.307001 0.192999], atol = 1e-13) + println(ode_algorithm.a_matrix) # Value in CI differs slightly from what I get locally + @test isapprox(ode_algorithm.a_matrix, + [0.27321088155198703 0.01250340416229867 + 0.4060225225166573 0.022548906054771216 + 0.534287756076577 0.03714081535199439 + 0.6549425779583463 0.05934313632736803 + 0.7621601562844809 0.09498270085837623 + 0.8446587253087918 0.1553412746912082 + 0.7719976108598626 0.22800238914013735 + 0.30700059728503437 0.1929994027149656], atol = 1e-13) end @testset "Sutherlands Law" begin