adjust the test and add a seed to the randomized initial guess for re…

…producibility
DanielDoehring · Jun 4, 2024 · 7bf717a · 7bf717a
1 parent 1332c28
commit 7bf717a
Show file tree

Hide file tree

Showing 5 changed files with 55 additions and 52 deletions.
diff --git 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"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
diff --git 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
@@ -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
@@ -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
@@ -1679,50 +1679,43 @@ 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
     dt_opt = 0.004485771991312504
     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