Add tests for alpha analysis

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)
 
 ###############################################################################

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)
 
 ###############################################################################

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)
 
 ###############################################################################

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)
 
 ###############################################################################

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,15 +156,17 @@ 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,
     alpha_mean, alpha_mean_pressure, alpha_mean_entropy = limiter.cache.subcell_limiter_coefficients
 
     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)

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

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

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

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