From 5895118c51f730fbd7880e359facfd519a3d1bc3 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Tue, 22 Oct 2024 14:51:40 +0200 Subject: [PATCH] Add tests for alpha analysis (#132) * Add tests for alpha analysis * Add comments; Fix two tests * Fix test * Update alpha tests * Rm not needed code * Simplify test --- ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 4 ++ .../elixir_euler_double_mach.jl | 4 ++ .../elixir_euler_double_mach_MCL.jl | 4 ++ .../elixir_euler_blast_wave_MCL.jl | 6 ++- src/callbacks_step/limiting_analysis.jl | 8 +-- src/callbacks_step/limiting_analysis_2d.jl | 15 ++++-- test/test_p4est_2d.jl | 17 +++++++ test/test_structured_2d.jl | 49 +++++++++++++++++++ test/test_tree_2d_euler.jl | 46 +++++++++++++++++ 9 files changed, 144 insertions(+), 9 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index 6fcb3b0fa92..f346f6ad59f 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -86,10 +86,14 @@ save_solution = SaveSolutionCallback(interval = 300, stepsize_callback = StepsizeCallback(cfl = 0.5) +limiting_analysis_callback = LimitingAnalysisCallback(output_directory = "out", + interval = 1) + callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + limiting_analysis_callback, stepsize_callback) ############################################################################### diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl index 33be05c615c..12bc936b1cc 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach.jl @@ -153,9 +153,13 @@ save_solution = SaveSolutionCallback(interval = 1000, stepsize_callback = StepsizeCallback(cfl = 0.9) +limiting_analysis_callback = LimitingAnalysisCallback(output_directory = "out", + interval = 1) + callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback, + limiting_analysis_callback, save_solution) ############################################################################### diff --git a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl index 68290ac31ff..9c2a03a26f0 100644 --- a/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl +++ b/examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl @@ -155,9 +155,13 @@ save_solution = SaveSolutionCallback(interval = 1000, stepsize_callback = StepsizeCallback(cfl = 0.9) +limiting_analysis_callback = LimitingAnalysisCallback(output_directory = "out", + interval = 1) + callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback, + limiting_analysis_callback, save_solution) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_MCL.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_MCL.jl index ba5f27d9f9f..78a4a72ed6e 100644 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_MCL.jl +++ b/examples/tree_2d_dgsem/elixir_euler_blast_wave_MCL.jl @@ -48,7 +48,7 @@ limiter_mcl = SubcellLimiterMCL(equations, basis; positivity_limiter_pressure_exact = false, entropy_limiter_semidiscrete = true, smoothness_indicator = true, - Plotting = false) + Plotting = true) volume_integral = VolumeIntegralSubcellLimiting(limiter_mcl; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) @@ -82,9 +82,13 @@ save_solution = SaveSolutionCallback(interval = 500, stepsize_callback = StepsizeCallback(cfl = 0.9) +limiting_analysis_callback = LimitingAnalysisCallback(output_directory = "out", + interval = 1) + callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + limiting_analysis_callback, stepsize_callback) ############################################################################### diff --git a/src/callbacks_step/limiting_analysis.jl b/src/callbacks_step/limiting_analysis.jl index 5bbf8c8b801..5dce041b549 100644 --- a/src/callbacks_step/limiting_analysis.jl +++ b/src/callbacks_step/limiting_analysis.jl @@ -145,7 +145,7 @@ end @unpack output_directory = limiting_analysis_callback @unpack alpha = limiter.cache.subcell_limiter_coefficients - alpha_avg = analyze_coefficient_IDP(mesh, equations, dg, cache, limiter) + alpha_avg = analyze_coefficient(mesh, equations, dg, cache, limiter) open("$output_directory/alphas.txt", "a") do f println(f, iter, ", ", time, ", ", maximum(alpha), ", ", alpha_avg) @@ -156,6 +156,8 @@ end dg, cache, limiter::SubcellLimiterMCL, time, iter) + @assert limiter.Plotting "Parameter `Plotting` needs to be activated for analysis of limiting factor with `LimitingAnalysisCallback`" + @unpack output_directory = limiting_analysis_callback @unpack weights = dg.basis @unpack alpha, alpha_pressure, alpha_entropy, @@ -163,8 +165,8 @@ end n_vars = nvariables(equations) - alpha_min_avg, alpha_mean_avg = analyze_coefficient_MCL(mesh, equations, dg, cache, - limiter) + alpha_min_avg, alpha_mean_avg = analyze_coefficient(mesh, equations, dg, cache, + limiter) open("$output_directory/alphas_min.txt", "a") do f print(f, iter, ", ", time) diff --git a/src/callbacks_step/limiting_analysis_2d.jl b/src/callbacks_step/limiting_analysis_2d.jl index ef547c01967..2b77257228d 100644 --- a/src/callbacks_step/limiting_analysis_2d.jl +++ b/src/callbacks_step/limiting_analysis_2d.jl @@ -5,7 +5,8 @@ @muladd begin #! format: noindent -function analyze_coefficient_IDP(mesh::TreeMesh2D, equations, dg, cache, limiter) +function analyze_coefficient(mesh::TreeMesh2D, equations, dg, cache, + limiter::SubcellLimiterIDP) @unpack weights = dg.basis @unpack alpha = limiter.cache.subcell_limiter_coefficients @@ -22,7 +23,9 @@ function analyze_coefficient_IDP(mesh::TreeMesh2D, equations, dg, cache, limiter return alpha_avg / total_volume end -function analyze_coefficient_IDP(mesh::StructuredMesh{2}, equations, dg, cache, limiter) +function analyze_coefficient(mesh::Union{StructuredMesh{2}, P4estMesh{2}}, + equations, dg, cache, + limiter::SubcellLimiterIDP) @unpack weights = dg.basis @unpack alpha = limiter.cache.subcell_limiter_coefficients @@ -39,7 +42,8 @@ function analyze_coefficient_IDP(mesh::StructuredMesh{2}, equations, dg, cache, return alpha_avg / total_volume end -function analyze_coefficient_MCL(mesh::TreeMesh2D, equations, dg, cache, limiter) +function analyze_coefficient(mesh::TreeMesh2D, equations, dg, cache, + limiter::SubcellLimiterMCL) @unpack weights = dg.basis @unpack alpha, alpha_mean, alpha_pressure, alpha_mean_pressure, alpha_entropy, alpha_mean_entropy = limiter.cache.subcell_limiter_coefficients @@ -83,8 +87,9 @@ function analyze_coefficient_MCL(mesh::TreeMesh2D, equations, dg, cache, limiter return alpha_avg ./ total_volume, alpha_mean_avg ./ total_volume end -function analyze_coefficient_MCL(mesh::StructuredMesh{2}, equations, dg, cache, - limiter) +function analyze_coefficient(mesh::Union{StructuredMesh{2}, P4estMesh{2}}, + equations, dg, cache, + limiter::SubcellLimiterMCL) @unpack weights = dg.basis @unpack alpha, alpha_mean, alpha_pressure, alpha_mean_pressure, alpha_entropy, alpha_mean_entropy = limiter.cache.subcell_limiter_coefficients diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index a8c20ee3bd1..c3c844b9ee7 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -335,6 +335,7 @@ end end @trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl" begin + rm(joinpath("out", "alphas.txt"), force = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ @@ -350,6 +351,22 @@ end 6.268843623142863 ], tspan=(0.0, 0.3)) + # Test alphas.txt + lines = readlines(joinpath("out", "alphas.txt")) + @test lines[1] == + "# iter, simu_time, alpha_max, alpha_avg" + cmd = string(Base.julia_cmd()) + coverage = occursin("--code-coverage", cmd) && + !occursin("--code-coverage=none", cmd) + if coverage + # Run with coverage takes 1 time step. + @test occursin(r"1, 0.00404[0-9]*, 1.0, 0.96795", lines[end]) + else + # Run without coverage takes 85 time steps. + @test startswith(lines[end], "85, 0.3, 1.0, 0.57771") + end + @test count(",", lines[end]) == 3 + @test !any(occursin.(r"NaN", lines)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 8254712cf3f..1e797e2015c 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -719,6 +719,7 @@ end end @trixi_testset "elixir_euler_double_mach.jl" begin + rm(joinpath("out", "alphas.txt"), force = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_double_mach.jl"), l2=[ 0.87417841433288, @@ -734,6 +735,22 @@ end ], initial_refinement_level=2, tspan=(0.0, 0.05)) + # Test alphas.txt + lines = readlines(joinpath("out", "alphas.txt")) + @test lines[1] == + "# iter, simu_time, alpha_max, alpha_avg" + cmd = string(Base.julia_cmd()) + coverage = occursin("--code-coverage", cmd) && + !occursin("--code-coverage=none", cmd) + if coverage + # Run with coverage takes 1 time steps. + @test occursin(r"1, 0.0002[0-9]*, 1.0, 0.9809", lines[end]) + else + # Run without coverage takes 193 time steps. + @test startswith(lines[end], "193, 0.05, 1.0, 0.3160") + end + @test count(",", lines[end]) == 3 + @test !any(occursin.(r"NaN", lines)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -749,6 +766,8 @@ end end @trixi_testset "elixir_euler_double_mach_MCL.jl" begin + rm(joinpath("out", "alphas_mean.txt"), force = true) + rm(joinpath("out", "alphas_min.txt"), force = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_double_mach_MCL.jl"), l2=[ 0.8887316108902574, @@ -764,6 +783,36 @@ end ], initial_refinement_level=2, tspan=(0.0, 0.05)) + # Test alphas_mean.txt + lines = readlines(joinpath("out", "alphas_mean.txt")) + @test lines[1] == + "# iter, simu_time, alpha_min_rho, alpha_avg_rho, alpha_min_rho_v1, alpha_avg_rho_v1, alpha_min_rho_v2, alpha_avg_rho_v2, alpha_min_rho_e, alpha_avg_rho_e" + cmd = string(Base.julia_cmd()) + coverage = occursin("--code-coverage", cmd) && + !occursin("--code-coverage=none", cmd) + if coverage + # Run with coverage takes 1 time steps. + @test startswith(lines[end], "1, 0.0002") + else + # Run without coverage takes 191 time steps. + @test startswith(lines[end], "191, 0.05, 3.70") + end + @test count(",", lines[end]) == 9 + @test !any(occursin.(r"NaN", lines)) + + # Test alphas_min.txt + lines = readlines(joinpath("out", "alphas_min.txt")) + @test lines[1] == + "# iter, simu_time, alpha_min_rho, alpha_avg_rho, alpha_min_rho_v1, alpha_avg_rho_v1, alpha_min_rho_v2, alpha_avg_rho_v2, alpha_min_rho_e, alpha_avg_rho_e" + if coverage + # Run with coverage takes 1 time steps. + @test startswith(lines[end], "1, 0.0002") # TODO + else + # Run without coverage takes 191 time steps. + @test startswith(lines[end], "191, 0.05, -0.0, 0.7216") + end + @test count(",", lines[end]) == 9 + @test !any(occursin.(r"NaN", lines)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 65448da0e3e..2874439b47b 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -532,6 +532,7 @@ end @trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl" begin rm(joinpath("out", "deviations.txt"), force = true) + rm(joinpath("out", "alphas.txt"), force = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_sc_subcell.jl"), l2=[ @@ -563,6 +564,20 @@ end # Run without coverage takes 381 time steps. @test startswith(lines[end], "381") end + + # Test alphas.txt + lines = readlines(joinpath("out", "alphas.txt")) + @test lines[1] == + "# iter, simu_time, alpha_max, alpha_avg" + if coverage + # Run with coverage takes 6 time steps. + @test occursin(r"6, 0.014[0-9]*, 1.0, 0.953", lines[end]) + else + # Run without coverage takes 381 time steps. + @test startswith(lines[end], "381, 1.0, 1.0, 0.544") + end + @test count(",", lines[end]) == 3 + @test !any(occursin.(r"NaN", lines)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -579,6 +594,8 @@ end @trixi_testset "elixir_euler_sedov_blast_wave_MCL.jl" begin rm(joinpath("out", "deviations.txt"), force = true) + rm(joinpath("out", "alphas_mean.txt"), force = true) + rm(joinpath("out", "alphas_min.txt"), force = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_MCL.jl"), l2=[ 0.4740321851943766, @@ -596,6 +613,7 @@ end initial_refinement_level=4, coverage_override=(maxiters = 6,), save_errors=true) + # Test deviations.txt lines = readlines(joinpath("out", "deviations.txt")) @test lines[1] == "# iter, simu_time, rho_min, rho_max, rho_v1_min, rho_v1_max, rho_v2_min, rho_v2_max, rho_e_min, rho_e_max, pressure_min" @@ -609,6 +627,34 @@ end # Run without coverage takes 349 time steps. @test startswith(lines[end], "349") end + + # Test alphas_mean.txt + lines = readlines(joinpath("out", "alphas_mean.txt")) + @test lines[1] == + "# iter, simu_time, alpha_min_rho, alpha_avg_rho, alpha_min_rho_v1, alpha_avg_rho_v1, alpha_min_rho_v2, alpha_avg_rho_v2, alpha_min_rho_e, alpha_avg_rho_e, alpha_min_pressure, alpha_avg_pressure, alpha_min_entropy, alpha_avg_entropy" + if coverage + # Run with coverage takes 6 time steps. + @test startswith(lines[end], "6, 0.014") + else + # Run without coverage takes 349 time steps. + @test startswith(lines[end], "349, 1.0, 0.0002") + end + @test count(",", lines[end]) == 13 + @test !any(occursin.(r"NaN", lines)) + + # Test alphas_min.txt + lines = readlines(joinpath("out", "alphas_min.txt")) + @test lines[1] == + "# iter, simu_time, alpha_min_rho, alpha_avg_rho, alpha_min_rho_v1, alpha_avg_rho_v1, alpha_min_rho_v2, alpha_avg_rho_v2, alpha_min_rho_e, alpha_avg_rho_e, alpha_min_pressure, alpha_avg_pressure, alpha_min_entropy, alpha_avg_entropy" + if coverage + # Run with coverage takes 6 time steps. + @test startswith(lines[end], "6, 0.014") + else + # Run without coverage takes 349 time steps. + @test startswith(lines[end], "349, 1.0, -0.0, 0.773") + end + @test count(",", lines[end]) == 13 + @test !any(occursin.(r"NaN", lines)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let