From bff694dd718e594784510b7394c99d67f5e23cc8 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 19 Dec 2023 10:29:40 +0100 Subject: [PATCH] Update elixir_navierstokes_3d_blast_wave_amr.jl --- .../elixir_navierstokes_3d_blast_wave_amr.jl | 176 +++++++----------- 1 file changed, 72 insertions(+), 104 deletions(-) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl index 935ff8d59c5..cb28987143c 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl @@ -1,3 +1,4 @@ + using OrdinaryDiffEq using Trixi @@ -9,36 +10,36 @@ prandtl_number() = 0.72 mu() = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = - CompressibleNavierStokesDiffusion3D(equations, mu = mu(), Prandtl = prandtl_number()) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), + Prandtl=prandtl_number()) function initial_condition_3d_blast_wave(x, t, equations::CompressibleEulerEquations3D) - - rho_c = 1.0 - p_c = 1.0 - u_c = 0.0 - - rho_o = 0.125 - p_o = 0.1 - u_o = 0.0 - - rc = 0.5 - r = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - if r < rc - rho = rho_c - v1 = u_c - v2 = u_c - v3 = u_c - p = p_c - else - rho = rho_o - v1 = u_o - v2 = u_o - v3 = u_o - p = p_o - end - - return prim2cons(SVector(rho, v1, v2, v3, p), equations) + + rho_c = 1.0 + p_c = 1.0 + u_c = 0.0 + + rho_o = 0.125 + p_o = 0.1 + u_o = 0.0 + + rc = 0.5 + r = sqrt(x[1]^2+x[2]^2+x[3]^2) + if r < rc + rho = rho_c + v1 = u_c + v2 = u_c + v3 = u_c + p = p_c + else + rho = rho_o + v1 = u_o + v2 = u_o + v3 = u_o + p = p_o + end + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) end initial_condition = initial_condition_3d_blast_wave @@ -46,44 +47,29 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha polydeg = 3 basis = LobattoLegendreBasis(polydeg) -indicator_sc = IndicatorHennemannGassner( - equations, - basis, - alpha_max = 1.0, - 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(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max=1.0, + 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(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral) coordinates_min = (-1.0, -1.0, -1.0) .* pi -coordinates_max = (1.0, 1.0, 1.0) .* pi +coordinates_max = ( 1.0, 1.0, 1.0) .* pi trees_per_dimension = (4, 4, 4) -mesh = P4estMesh( - trees_per_dimension, - polydeg = 3, - coordinates_min = coordinates_min, - coordinates_max = coordinates_max, - periodicity = (true, true, true), - initial_refinement_level = 1, -) +mesh = P4estMesh(trees_per_dimension, polydeg=3, + coordinates_min=coordinates_min, coordinates_max=coordinates_max, + periodicity=(true, true, true), initial_refinement_level=1) -semi = SemidiscretizationHyperbolicParabolic( - mesh, - (equations, equations_parabolic), - initial_condition, - solver, -) +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) ############################################################################### # ODE solvers, callbacks etc. @@ -94,54 +80,36 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() analysis_interval = 50 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -save_solution = SaveSolutionCallback( - interval = 20, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim, -) - -amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) - -amr_controller = ControllerThreeLevel( - semi, - amr_indicator, - base_level = 0, - med_level = 1, - med_threshold = 0.05, - max_level = 4, - max_threshold = 0.1, -) -amr_callback = AMRCallback( - semi, - amr_controller, - interval = 5, - adapt_initial_condition = true, - adapt_initial_condition_only_refine = true, -) - - -alive_callback = AliveCallback(analysis_interval = analysis_interval) - -callbacks = CallbackSet( - summary_callback, - analysis_callback, - alive_callback, - amr_callback, - save_solution, -) +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +save_solution = SaveSolutionCallback(interval=20, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +amr_indicator = IndicatorLöhner(semi, variable=Trixi.density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level=0, + med_level=1, med_threshold=0.05, + max_level=4, max_threshold=0.1) +amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + + +alive_callback = AliveCallback(analysis_interval=analysis_interval,) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + save_solution) ############################################################################### # run the simulation time_int_tol = 1e-8 -sol = solve( - ode, - RDPK3SpFSAL49(); - abstol = time_int_tol, - reltol = time_int_tol, - ode_default_options()..., - callback = callbacks, -) +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) summary_callback() # print the timer summary