WIP: In-situ visualization #1: Make Makie available in VisualizationCallback

Project.toml

@@ -83,7 +83,7 @@ LinearAlgebra = "1"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.151"
 MPI = "0.20.6"
-Makie = "0.19, 0.20"
+Makie = "0.22"
 MuladdMacro = "0.2.2"
 NLsolve = "4.5.1"
 Octavian = "0.3.21"
@@ -121,4 +121,4 @@ julia = "1.8"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
-NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"

examples/structured_2d_dgsem/elixir_advection_basic_visualization.jl

@@ -0,0 +1,67 @@
+# The same setup as tree_2d_dgsem/elixir_advection_basic.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEq
+using Trixi
+using Plots
+using GLMakie
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+cells_per_dimension = (16, 16)
+
+# Create curved mesh with 16 x 16 elements
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 1.0))
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
+# Enable in-situ visualization with a new plot generated every 20 time steps
+# and additional plotting options passed as keyword arguments
+visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_creator = Trixi.show_plot_makie, show_mesh = true)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        visualization, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

examples/tree_2d_dgsem/elixir_advection_amr_visualization.jl

@@ -2,6 +2,7 @@
 using OrdinaryDiffEq
 using Trixi
 using Plots
+using GLMakie
 
 ###############################################################################
 # semidiscretization of the linear advection equation
@@ -24,15 +25,15 @@ solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
 coordinates_min = (-5.0, -5.0)
 coordinates_max = (5.0, 5.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 3,
+                initial_refinement_level = 4,
                 n_cells_max = 30_000)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 20.0)
+tspan = (0.0, 100.0)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
@@ -50,12 +51,12 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 # Enable in-situ visualization with a new plot generated every 20 time steps
 # and additional plotting options passed as keyword arguments
-visualization = VisualizationCallback(interval = 20, clims = (0, 1))
+visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_creator = Trixi.show_plot_makie, show_mesh = true)
 
 amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
-                                      base_level = 3,
-                                      med_level = 4, med_threshold = 0.1,
-                                      max_level = 5, max_threshold = 0.6)
+                                      base_level = 4,
+                                      med_level = 5, med_threshold = 0.1,
+                                      max_level = 6, max_threshold = 0.6)
 amr_callback = AMRCallback(semi, amr_controller,
                            interval = 5,
                            adapt_initial_condition = true,

examples/tree_2d_dgsem/elixir_advection_amr_visualization_makie.jl

@@ -0,0 +1,77 @@
+using OrdinaryDiffEq
+using Trixi
+using GLMakie
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+function initial_condition_gauss_largedomain(x, t,
+                                             equation::LinearScalarAdvectionEquation2D)
+    # Store translated coordinate for easy use of exact solution
+    domain_length = SVector(10, 10)
+    x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t, domain_length)
+
+    return SVector(exp(-(x_trans[1]^2 + x_trans[2]^2)))
+end
+initial_condition = initial_condition_gauss_largedomain
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-5.0, -5.0)
+coordinates_max = (5.0, 5.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 20.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+# Enable in-situ visualization with a new plot generated every 20 time steps
+# and additional plotting options passed as keyword arguments
+visualization = VisualizationCallback(interval = 100, clims = (0, 1),
+                                      plot_creator = Trixi.show_plot_makie)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 3,
+                                      med_level = 4, med_threshold = 0.1,
+                                      max_level = 5, max_threshold = 0.6)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution, visualization,
+                        amr_callback, stepsize_callback);
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_2d_dgsem/elixir_euler_blast_wave_amr_visualization.jl

@@ -0,0 +1,109 @@
+
+using OrdinaryDiffEq
+using Trixi
+using GLMakie
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+"""
+    initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+
+A medium blast wave taken from
+- Sebastian Hennemann, Gregor J. Gassner (2020)
+  A provably entropy stable subcell shock capturing approach for high order split form DG
+  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
+"""
+function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+    phi = atan(y_norm, x_norm)
+    sin_phi, cos_phi = sincos(phi)
+
+    # Calculate primitive variables
+    rho = r > 0.5 ? 1.0 : 1.1691
+    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
+    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
+    p = r > 0.5 ? 1.0E-3 : 1.245
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-2.0, -2.0)
+coordinates_max = (2.0, 2.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 12.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_indicator = IndicatorHennemannGassner(semi,
+                                          alpha_max = 0.5,
+                                          alpha_min = 0.001,
+                                          alpha_smooth = true,
+                                          variable = density_pressure)
+
+
+visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_creator = Trixi.show_plot_makie)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 4,
+                                      max_level = 6, max_threshold = 0.01)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        visualization,
+                        amr_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_3d_dgsem/elixir_advection_amr_visualization.jl

@@ -0,0 +1,72 @@
+
+using OrdinaryDiffEq
+using Trixi
+using Plots
+using GLMakie
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (20, -7, 5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_gauss
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-5.0, -5.0, -5.0)
+coordinates_max = (5.0, 5.0, 5.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 2,
+                n_cells_max = 30_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 300)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+
+# visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_data_creator = Trixi.PlotData3D, plot_creator = Trixi.show_plot_makie)
+visualization = VisualizationCallback(interval = 20, clims = (0, 1), plot_data_creator = Trixi.PlotData3D, plot_creator = Trixi.show_plot_makie_slicing)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 2,
+                                      med_level = 3, med_threshold = 0.1,
+                                      max_level = 4, max_threshold = 0.6)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        visualization,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary