diff --git a/examples/tree_3d_dgsem/elixir_euler_amr_visualization.jl b/examples/tree_3d_dgsem/elixir_euler_amr_visualization.jl new file mode 100644 index 0000000000..643ea465db --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_euler_amr_visualization.jl @@ -0,0 +1,87 @@ + +using OrdinaryDiffEq +using Trixi +using Plots +using GLMakie + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +""" + initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D) + +A Gaussian pulse in the density with constant velocity and pressure; reduces the +compressible Euler equations to the linear advection equations. +""" +function initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D) + rho = 1 + exp(-(x[1]^2 + x[2]^2 + x[3]^2)) / 2 + v1 = 1 + v2 = 1 + v3 = 1 + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_v3 = rho * v3 + p = 1 + rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2 + v3^2) + return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e) +end +initial_condition = initial_condition_density_pulse +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 = 4, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.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_restart = SaveRestartCallback(interval = 100, + save_final_restart = true) + +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) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 1.05, + max_level = 3, max_threshold = 1.3) +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_restart, 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 diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index 8ebe69e0ca..e03128ef72 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -14,6 +14,7 @@ mutable struct VisualizationCallback{SolutionVariables, VariableNames, PlotDataC plot_data_creator::PlotDataCreator plot_creator::PlotCreator figure + axis plot_arguments::Dict{Symbol, Any} end @@ -99,7 +100,7 @@ function VisualizationCallback(; interval = 0, visualization_callback = VisualizationCallback(interval, solution_variables, variable_names, show_mesh, - plot_data_creator, plot_creator, nothing, + plot_data_creator, plot_creator, nothing, [], Dict{Symbol, Any}(plot_arguments)) DiscreteCallback(visualization_callback, visualization_callback, # the first one is the condition, the second the affect! @@ -137,7 +138,7 @@ end function (visualization_callback::VisualizationCallback)(integrator) u_ode = integrator.u semi = integrator.p - @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure = visualization_callback + @unpack plot_arguments, solution_variables, variable_names, show_mesh, plot_data_creator, plot_creator, figure, axis = visualization_callback mesh, equations, solver, cache = mesh_equations_solver_cache(integrator.p) n = Trixi.ndims(mesh) @@ -152,7 +153,7 @@ function (visualization_callback::VisualizationCallback)(integrator) # Create plot plot_creator(n, visualization_callback, plot_data, variable_names; show_mesh = show_mesh, plot_arguments = plot_arguments, - time = integrator.t, timestep = integrator.stats.naccept, figure = figure) + time = integrator.t, timestep = integrator.stats.naccept, figure = figure, axis = axis) # avoid re-evaluating possible FSAL stages u_modified!(integrator, false) @@ -178,7 +179,7 @@ See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) """ function show_plot(ndims, visualization_callback, plot_data, variable_names; show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) + time = nothing, timestep = nothing, figure = nothing, axis = nothing) # Gather subplots plots = [] for v in variable_names @@ -289,30 +290,25 @@ See also: [`VisualizationCallback`](@ref), [`save_plot`](@ref) """ function show_plot_makie(ndims, visualization_callback, plot_data, variable_names; show_mesh = true, plot_arguments = Dict{Symbol, Any}(), - time = nothing, timestep = nothing, figure = nothing) - + time = nothing, timestep = nothing, figure = nothing, axis = []) if figure === nothing @warn "Creating new figure" visualization_callback.figure = GLMakie.Figure() figure = visualization_callback.figure - if ndims == 2 - for v in 1:size(variable_names)[1] - GLMakie.heatmap(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v], axis=(; title = variable_names[v])) - end - else - for v in 1:size(variable_names)[1] - GLMakie.volume(figure[makieLayoutHelper(v)...], plot_data.data[v], axis=(; title = variable_names[v])) - end + for v in 1:size(variable_names)[1] + push!(axis, (ndims == 2) ? GLMakie.Axis(figure[makieLayoutHelper(v)...], title = variable_names[v]) : + GLMakie.Axis3(figure[makieLayoutHelper(v)...], aspect=:equal, title = variable_names[v])) end + visualization_callback.axis = axis GLMakie.display(figure) else if ndims == 2 for v in 1:size(variable_names)[1] - GLMakie.heatmap!(figure[makieLayoutHelper(v)...], plot_data.x, plot_data.y, plot_data.data[v]) + GLMakie.heatmap!(axis[v], plot_data.x, plot_data.y, plot_data.data[v]) end else for v in 1:size(variable_names)[1] - GLMakie.volume!(figure[makieLayoutHelper(v)...], plot_data.data[v]) + GLMakie.volume!(axis[v], plot_data.data[v]) end end end