WIP: ES volume integral

Project.toml

@@ -14,12 +14,16 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 

examples/2d/elixir_euler_kelvin_helmholtz_instability.jl

@@ -16,13 +16,13 @@ volume_flux  = flux_chandrashekar
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 indicator_sc = IndicatorHennemannGassner(equations, basis,
-                                         alpha_max=0.002,
-                                         alpha_min=0.0001,
-                                         alpha_smooth=true,
-                                         variable=density_pressure)
+                                        alpha_max=0.002,
+                                        alpha_min=0.0001,
+                                        alpha_smooth=true,
+                                        variable=density_pressure)
 volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
-                                                 volume_flux_dg=volume_flux,
-                                                 volume_flux_fv=surface_flux)
+                                                volume_flux_dg=volume_flux,
+                                                volume_flux_fv=surface_flux)
 solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinates_min = (-0.5, -0.5)
@@ -53,7 +53,7 @@ save_solution = SaveSolutionCallback(interval=20,
 stepsize_callback = StepsizeCallback(cfl=1.3)
 
 callbacks = CallbackSet(summary_callback,
-                        analysis_callback, alive_callback, 
+                        analysis_callback, alive_callback,
                         save_solution,
                         stepsize_callback)
 

examples/2d/elixir_euler_kelvin_helmholtz_instability2.jl

@@ -0,0 +1,137 @@
+
+using Random: seed!
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+seed!(0)
+function initial_condition_khi2(x, t, equations::CompressibleEulerEquations2D)
+  # change discontinuity to tanh
+  # typical resolution 128^2, 256^2
+  # domain size is [-1,+1]^2
+  slope = 15
+  B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
+  rho = 0.5 + 0.75 * B
+  v1 = 0.5 * (B - 1)
+  v2 = 0.1 * sin(1 * pi * x[1]) #+ 0.05 * sin(4 * pi * x[1])
+  p = 1.0
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+
+
+initial_condition = initial_condition_khi2
+
+#surface_flux = flux_lax_friedrichs
+surface_flux = flux_hllc
+#surface_flux = flux_hll
+#surface_flux = flux_shima_etal
+volume_flux  = flux_ranocha
+#volume_flux  = flux_shima_etal
+#volume_flux  = flux_central
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+#indicator_sc = IndicatorHennemannGassner(equations, basis,
+#                                        alpha_max=0.002,
+#                                        alpha_min=0.0001,
+#                                        alpha_smooth=true,
+#                                        variable=density_pressure)
+#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+#                                                volume_flux_dg=volume_flux,
+#                                                volume_flux_fv=surface_flux)
+
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+#volume_integral = VolumeIntegralWeakForm()
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+#surface_flux  = flux_hllc
+#volume_flux  = flux_ranocha
+#solver = DGSEM(7, surface_flux, VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux)))
+
+#surface_flux  = flux_hllc
+#volume_flux = flux_ranocha
+#volume_integral = VolumeIntegralFluxComparison(flux_central, volume_flux)
+#solver = DGSEM(15, surface_flux, volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = ( 1.0,  1.0)
+
+#restart_filename = joinpath("out", "solution_010700.h5")
+#mesh = load_mesh(restart_filename, n_cells_max=800_000)
+
+println("start mesh")
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=800_000)
+println("end mesh")
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+println("end semi")
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.7)
+#tspan = (load_time(restart_filename), 3.7)
+ode = semidiscretize(semi, tspan)
+#ode = semidiscretize(semi, tspan, restart_filename)
+println("end ode")
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
+                                     extra_analysis_errors=(:conservation_error,),
+                                     extra_analysis_integrals=(entropy, energy_total,
+                                                               energy_kinetic, energy_internal))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=1000,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+save_restart = SaveRestartCallback(interval=1000)
+
+stepsize_callback = StepsizeCallback(cfl=0.2)
+
+#amr_indicator = IndicatorLöhner(semi,
+#                                variable=Trixi.density)
+#amr_controller = ControllerThreeLevel(semi, amr_indicator,
+#                                      base_level=1,
+#                                      med_level =0, med_threshold=0.1, # med_level = current level
+#                                      max_level =6, max_threshold=0.5)
+#amr_callback = AMRCallback(semi, amr_controller,
+#                           interval=3,
+#                           adapt_initial_condition=false,
+#                           adapt_initial_condition_only_refine=true)
+
+limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-5, 5.0e-5),
+                                               variables=(Trixi.density, pressure))
+stage_limiter! = limiter!
+step_limiter!  = limiter!
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        save_restart,
+#                        amr_callback, 
+#			stepsize_callback
+			)
+
+println("end callbacks")
+
+###############################################################################
+# run the simulation
+
+#sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+#sol = solve(ode, SSPRK33(stage_limiter!, step_limiter!),
+sol = solve(ode, SSPRK43(stage_limiter!, step_limiter!),
+#sol = solve(ode, SSPRK33(),
+            #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/2d/elixir_euler_kelvin_helmholtz_instability3.jl

@@ -0,0 +1,160 @@
+
+using Random: seed!
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 5/3
+equations = CompressibleEulerEquations2D(gamma)
+
+seed!(0)
+function initial_condition_khi3(x, t, equations::CompressibleEulerEquations2D)
+  # from https://arxiv.org/pdf/2001.01927.pdf
+  # typical resolution 128^2, 256^2
+  # domain size is [0,1]^2
+  L = 0.025
+  U1 = 0.5
+  U2 =-0.5
+  Um = 0.5 * (U1-U2)
+  rho1 = 1.0
+  rho2 = 2.0
+  rhom = 0.5 * (rho1 - rho2)
+
+  v2 = 0.01 * sin(4 * pi * x[1])
+  p = 2.5
+
+  if (0 <= x[2] <= 0.25)
+    rho = rho1 - rhom * exp((x[2] - 0.25) / L)
+    v1  = U1 - Um * exp((x[2] - 0.25) / L)
+  elseif (0.25 < x[2] <= 0.5)
+    rho = rho2 + rhom * exp((-x[2] + 0.25) / L)
+    v1  = U2 + Um * exp((-x[2] + 0.25) / L)
+  elseif (0.5 < x[2] <= 0.75)
+    rho = rho2 + rhom * exp((x[2] - 0.75) / L)
+    v1  = U2 + Um * exp((x[2] - 0.75) / L)
+  elseif (0.75 < x[2] <= 1.0)
+    rho = rho1 - rhom * exp((-x[2] + 0.75) / L)
+    v1  = U1 - Um * exp((-x[2] + 0.75) / L)
+  else
+    println("shyte", x[2])
+  end
+  if (rho < 0.0)
+   println("density", x[2], rho)
+  end
+  if (p < 0.0)
+   println("pressure")
+  end
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+
+
+initial_condition = initial_condition_khi3
+
+#surface_flux = flux_lax_friedrichs
+surface_flux = flux_hllc
+#surface_flux = flux_hll
+#surface_flux = flux_shima_etal
+volume_flux  = flux_ranocha
+#volume_flux  = flux_shima_etal
+#volume_flux  = flux_central
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+#indicator_sc = IndicatorHennemannGassner(equations, basis,
+#                                        alpha_max=0.002,
+#                                        alpha_min=0.0001,
+#                                        alpha_smooth=true,
+#                                        variable=density_pressure)
+#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+#                                                volume_flux_dg=volume_flux,
+#                                                volume_flux_fv=surface_flux)
+
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+#volume_integral = VolumeIntegralWeakForm()
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+#surface_flux  = flux_hllc
+#volume_flux  = flux_ranocha
+#solver = DGSEM(7, surface_flux, VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux)))
+
+#surface_flux  = flux_hllc
+#volume_flux = flux_ranocha
+#volume_integral = VolumeIntegralFluxComparison(flux_central, volume_flux)
+#solver = DGSEM(15, surface_flux, volume_integral)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = ( 1.0,  1.0)
+
+#restart_filename = joinpath("out", "solution_010700.h5")
+#mesh = load_mesh(restart_filename, n_cells_max=800_000)
+
+println("start mesh")
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=7,
+                n_cells_max=4_000_000)
+println("end mesh")
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.7)
+#tspan = (load_time(restart_filename), 3.7)
+ode = semidiscretize(semi, tspan)
+#ode = semidiscretize(semi, tspan, restart_filename)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
+                                     extra_analysis_errors=(:conservation_error,),
+                                     extra_analysis_integrals=(entropy, energy_total,
+                                                               energy_kinetic, energy_internal))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=1000,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+save_restart = SaveRestartCallback(interval=1000)
+
+stepsize_callback = StepsizeCallback(cfl=0.8)
+
+#amr_indicator = IndicatorLöhner(semi,
+#                                variable=Trixi.density)
+#amr_controller = ControllerThreeLevel(semi, amr_indicator,
+#                                      base_level=1,
+#                                      med_level =0, med_threshold=0.1, # med_level = current level
+#                                      max_level =6, max_threshold=0.5)
+#amr_callback = AMRCallback(semi, amr_controller,
+#                           interval=3,
+#                           adapt_initial_condition=false,
+#                           adapt_initial_condition_only_refine=true)
+
+limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-5, 5.0e-5),
+                                               variables=(Trixi.density, pressure))
+stage_limiter! = limiter!
+step_limiter!  = limiter!
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        save_restart,
+#                        amr_callback, 
+#			stepsize_callback
+			)
+
+
+###############################################################################
+# run the simulation
+
+#sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+#sol = solve(ode, SSPRK33(stage_limiter!, step_limiter!),
+sol = solve(ode, SSPRK43(stage_limiter!, step_limiter!),
+#sol = solve(ode, SSPRK33(),
+            #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