Curvilinear FDSBP

examples/unstructured_2d_fdsbp/elixir_euler_convergence.jl

@@ -0,0 +1,64 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the FDSBP approximation operator
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order=1, accuracy_order=4,
+                            xmin=-1.0, xmax=1.0, N=12)
+solver = FDSBP(D_SBP,
+               surface_integral=SurfaceIntegralStrongForm(flux_lax_friedrichs),
+               volume_integral=VolumeIntegralStrongForm())
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh",
+                                       default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity=true)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+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)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol=1.0e-9, reltol=1.0e-9,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary

examples/unstructured_2d_fdsbp/elixir_euler_convergence_upwind.jl

@@ -0,0 +1,111 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Modify the manufactured solution test to use `L = sqrt(2)` in the initial and source terms
+# such that testing works on the flipped mesh
+function initial_condition_convergence_upwind(x, t, equations::CompressibleEulerEquations2D)
+  c = 2
+  A = 0.1
+  L = sqrt(2)
+  f = 1/L
+  ω = 2 * pi * f
+  ini = c + A * sin(ω * (x[1] + x[2] - t))
+
+  rho = ini
+  rho_v1 = ini
+  rho_v2 = ini
+  rho_e = ini^2
+
+  return SVector(rho, rho_v1, rho_v2, rho_e)
+end
+
+@inline function source_terms_convergence_upwind(u, x, t, equations::CompressibleEulerEquations2D)
+  # Same settings as in `initial_condition`
+  c = 2
+  A = 0.1
+  L = sqrt(2)
+  f = 1/L
+  ω = 2 * pi * f
+  γ = equations.gamma
+
+  x1, x2 = x
+  si, co = sincos(ω * (x1 + x2 - t))
+  rho = c + A * si
+  rho_x = ω * A * co
+  # Note that d/dt rho = -d/dx rho = -d/dy rho.
+
+  tmp = (2 * rho - 1) * (γ - 1)
+
+  du1 = rho_x
+  du2 = rho_x * (1 + tmp)
+  du3 = du2
+  du4 = 2 * rho_x * (rho + tmp)
+
+  return SVector(du1, du2, du3, du4)
+end
+
+initial_condition = initial_condition_convergence_upwind
+
+###############################################################################
+# Get the FDSBP approximation operator
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1,
+                         accuracy_order=4,
+                         xmin=-1.0, xmax=1.0,
+                         N=10)
+
+flux_splitting = splitting_lax_friedrichs
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+default_mesh_file = joinpath(@__DIR__, "mesh_multiple_flips.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/9a967f066bc5bf081e77ef2519b3918e40a964ed/mesh_multiple_flips.mesh",
+                                       default_mesh_file)
+
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity=true)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_upwind)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+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)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(), abstol=1.0e-9, reltol=1.0e-9,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary

examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl

@@ -0,0 +1,76 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Free-stream initial condition
+initial_condition = initial_condition_constant
+
+# Boundary conditions for free-stream testing
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict( :Body    => boundary_condition_free_stream,
+                            :Button1 => boundary_condition_free_stream,
+                            :Button2 => boundary_condition_free_stream,
+                            :Eye1    => boundary_condition_free_stream,
+                            :Eye2    => boundary_condition_free_stream,
+                            :Smile   => boundary_condition_free_stream,
+                            :Bowtie  => boundary_condition_free_stream )
+
+###############################################################################
+# Get the FDSBP approximation space
+
+D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
+                            derivative_order=1, accuracy_order=4,
+                            xmin=-1.0, xmax=1.0, N=9)
+solver = FDSBP(D_SBP,
+               surface_integral=SurfaceIntegralStrongForm(flux_hll),
+               volume_integral=VolumeIntegralStrongForm())
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+default_mesh_file = joinpath(@__DIR__, "mesh_gingerbread_man.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh",
+                                       default_mesh_file)
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+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)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# set small tolerances for the free-stream preservation test
+sol = solve(ode, SSPRK43(), abstol=1.0e-12, reltol=1.0e-12,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary

examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl

@@ -0,0 +1,78 @@
+
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Free-stream initial condition
+initial_condition = initial_condition_constant
+
+# Boundary conditions for free-stream testing
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+
+boundary_conditions = Dict( :Top    => boundary_condition_free_stream,
+                            :Bottom => boundary_condition_free_stream,
+                            :Right  => boundary_condition_free_stream,
+                            :Left   => boundary_condition_free_stream )
+
+###############################################################################
+# Get the Upwind FDSBP approximation space
+
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order=1,
+                         accuracy_order=4,
+                         xmin=-1.0, xmax=1.0,
+                         N=9)
+
+flux_splitting = splitting_lax_friedrichs
+solver = FDSBP(D_upw,
+               surface_integral=SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
+               volume_integral=VolumeIntegralUpwind(flux_splitting))
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+default_mesh_file = joinpath(@__DIR__, "mesh_multiple_flips.mesh")
+isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/9a967f066bc5bf081e77ef2519b3918e40a964ed/mesh_multiple_flips.mesh",
+                                       default_mesh_file)
+
+mesh_file = default_mesh_file
+
+mesh = UnstructuredMesh2D(mesh_file)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=1000,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback,
+                        alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# set small tolerances for the free-stream preservation test
+sol = solve(ode, SSPRK43(), abstol=1.0e-12, reltol=1.0e-12,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary

src/equations/compressible_euler_2d.jl

@@ -780,7 +780,7 @@ end
     rho_2gamma = 0.5 * rho / equations.gamma
     f1m = rho_2gamma * alpha_m
     f2m = rho_2gamma * alpha_m * v1
-    f3m = rho_2gamma * (alpha_m * v2 + a * (lambda2_m-lambda3_m))
+    f3m = rho_2gamma * (alpha_m * v2 + a * (lambda2_m - lambda3_m))
     f4m = rho_2gamma * (alpha_m * 0.5 * (v1^2 + v2^2) + a * v2 * (lambda2_m - lambda3_m)
                         + a^2 * (lambda2_m + lambda3_m) * equations.inv_gamma_minus_one)
   end
@@ -970,6 +970,46 @@ end
   return SVector(f1m, f2m, f3m, f4m)
 end
 
+@inline function splitting_lax_friedrichs(u, ::Val{:plus}, normal_direction::AbstractVector,
+                                          equations::CompressibleEulerEquations2D)
+  rho, rho_v1, rho_v2, rho_e = u
+  v1 = rho_v1 / rho
+  v2 = rho_v2 / rho
+  v_dot_n = normal_direction[1] * v1 + normal_direction[2] * v2
+  p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+
+  a = sqrt(equations.gamma * p / rho)
+  H = (rho_e + p) / rho
+  lambda = 0.5 * (v_dot_n + a * norm(normal_direction))
+
+  f1p = 0.5 * rho * v_dot_n + lambda * u[1]
+  f2p = 0.5 * rho * v_dot_n * v1 + 0.5 * p * normal_direction[1] + lambda * u[2]
+  f3p = 0.5 * rho * v_dot_n * v2 + 0.5 * p * normal_direction[2] + lambda * u[3]
+  f4p = 0.5 * rho * v_dot_n * H + lambda * u[4]
+
+  return SVector(f1p, f2p, f3p, f4p)
+end
+
+@inline function splitting_lax_friedrichs(u, ::Val{:minus}, normal_direction::AbstractVector,
+                                          equations::CompressibleEulerEquations2D)
+  rho, rho_v1, rho_v2, rho_e = u
+  v1 = rho_v1 / rho
+  v2 = rho_v2 / rho
+  v_dot_n = normal_direction[1] * v1 + normal_direction[2] * v2
+  p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+
+  a = sqrt(equations.gamma * p / rho)
+  H = (rho_e + p) / rho
+  lambda = 0.5 * (v_dot_n + a * norm(normal_direction))
+
+  f1m = 0.5 * rho * v_dot_n - lambda * u[1]
+  f2m = 0.5 * rho * v_dot_n * v1 + 0.5 * p * normal_direction[1] - lambda * u[2]
+  f3m = 0.5 * rho * v_dot_n * v2 + 0.5 * p * normal_direction[2] - lambda * u[3]
+  f4m = 0.5 * rho * v_dot_n * H - lambda * u[4]
+
+  return SVector(f1m, f2m, f3m, f4m)
+end
+
 
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
 # maximum velocity magnitude plus the maximum speed of sound