dirty implementation of Jin Xin Relaxation for Compressible Euler 2D

examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl

@@ -48,14 +48,14 @@ solver = DGSEM(basis, surface_flux, volume_integral)
 coordinates_min = (-1.0, -1.0)
 coordinates_max = ( 1.0,  1.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level=5,
+                initial_refinement_level=4,
                 n_cells_max=100_000)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 3.0)
+tspan = (0.0, 3.7)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
@@ -70,7 +70,7 @@ save_solution = SaveSolutionCallback(interval=20,
                                      save_final_solution=true,
                                      solution_variables=cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl=1.3)
+stepsize_callback = StepsizeCallback(cfl=0.9)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

examples/tree_2d_dgsem/elixir_jin_xin_euler.jl

@@ -0,0 +1,93 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+epsilon_relaxation = 1.0e-5
+a1 = a2 = a3 = a4 = 15.0
+b1 = b2 = b3 = b4 = 15.0
+
+equations_relaxation = CompressibleEulerEquations2D(1.4)
+equations = JinXinCompressibleEulerEquations2D(epsilon_relaxation, a1, a2, a3, a4, b1, b2, b3, b4,equations_relaxation)
+
+function initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
+  # change discontinuity to tanh
+  # typical resolution 128^2, 256^2
+  # domain size is [-1,+1]^2
+  slope = 15
+  amplitude = 0.02
+  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(2 * pi * x[1])
+  p = 1.0
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end 
+
+#initial_condition = initial_condition_constant
+initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability)
+solver = DGSEM(polydeg=3, surface_flux=Trixi.flux_upwind)
+
+#surface_flux = Trixi.flux_upwind
+#volume_flux  = flux_central
+#basis = LobattoLegendreBasis(7)
+#limiter_idp = SubcellLimiterIDP(equations, basis;
+#                                positivity_variables_cons=[1],
+#                                positivity_correction_factor=0.5)
+#volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+#                                                volume_flux_dg=volume_flux,
+#                                                volume_flux_fv=surface_flux)
+#solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = ( 1.0,  1.0)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=1_000_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)#,source_terms=source_terms_JinXin_Relaxation)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.7)
+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,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=0.1)
+
+collision_callback = LBMCollisionCallback()  
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,stepsize_callback)#,collision_callback)
+
+###############################################################################
+# run the simulation
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
+                                                     variables=(Trixi.density, pressure))
+
+#sol = solve(ode, CarpenterKennedy2N54(stage_limiter!,williamson_condition=false),
+sol = solve(ode, SSPRK43(stage_limiter!),
+#sol = solve(ode, SSPRK33(stage_limiter!),
+#sol = solve(ode, RDPK3SpFSAL49(),
+#sol = solve(ode, AutoTsit5(Rosenbrock23()),
+            dt=1.0e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks,maxiters=1e7);
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -150,7 +150,8 @@ export AcousticPerturbationEquations2D,
        ShallowWaterEquations1D, ShallowWaterEquations2D,
        ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
        ShallowWaterEquationsQuasi1D,
-       LinearizedEulerEquations2D
+       LinearizedEulerEquations2D,
+       JinXinCompressibleEulerEquations2D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
@@ -199,7 +200,7 @@ export boundary_condition_do_nothing,
        BoundaryConditionCoupled
 
 export initial_condition_convergence_test, source_terms_convergence_test
-export source_terms_harmonic
+export source_terms_harmonic, source_terms_JinXin_Relaxation
 export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic,
        boundary_condition_poisson_nonperiodic
 export initial_condition_eoc_test_coupled_euler_gravity,

src/callbacks_stage/positivity_zhang_shu.jl

@@ -33,8 +33,9 @@ function (limiter!::PositivityPreservingLimiterZhangShu)(u_ode, integrator,
                                                          t)
     u = wrap_array(u_ode, semi)
     @trixi_timeit timer() "positivity-preserving limiter" begin
-        limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables,
-                           mesh_equations_solver_cache(semi)...)
+        #limiter_zhang_shu!(u, limiter!.thresholds, limiter!.variables,
+        #                   mesh_equations_solver_cache(semi)...)
+        perform_idp_correction!(u, integrator.dt, mesh_equations_solver_cache(semi)...)  
     end
 end
 

src/callbacks_stage/positivity_zhang_shu_dg2d.jl

@@ -42,4 +42,40 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
     return nothing
 end
+function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompressibleEulerEquations2D, dg, cache)
+
+            # relaxation parameter
+            eps = equations.eps_relaxation
+            dt_ = dt
+            factor =1.0/ (eps + dt_)
+            eq_relax = equations.equations_relaxation
+
+ @threaded for element in eachelement(dg, cache)
+        for j in eachnode(dg), i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, j, element)
+
+            # compute compressible Euler fluxes
+            u_cons = SVector(u_node[1], u_node[2], u_node[3], u_node[4])
+            vu = flux(u_cons,1,eq_relax)
+            wu = flux(u_cons,2,eq_relax)
+            # compute relaxation terms
+            du1 = u_node[1]
+            du2 = u_node[2]
+            du3 = u_node[3]
+            du4 = u_node[4]
+            du5 = factor * (eps * u_node[5] + dt_ * vu[1])
+            du6 = factor * (eps * u_node[6] + dt_ * vu[2])
+            du7 = factor * (eps * u_node[7] + dt_ * vu[3])
+            du8 = factor * (eps * u_node[8] + dt_ * vu[4])
+            du9 = factor * (eps * u_node[9] + dt_ * wu[1])
+            du10= factor * (eps * u_node[10]+ dt_ * wu[2])
+            du11= factor * (eps * u_node[11]+ dt_ * wu[3])
+            du12= factor * (eps * u_node[12]+ dt_ * wu[4])
+            new_u = SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, du10, du11, du12)
+            set_node_vars!(u, new_u, equations, dg, i, j, element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -41,4 +41,40 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache)
 
     return nothing
 end
+function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompressibleEulerEquations2D, dg, cache)
+
+            # relaxation parameter
+            eps = equations.eps_relaxation
+            dt_ = dt * 1.0
+            factor =1.0/ (eps + dt_)
+
+ @threaded for element in eachelement(dg, cache)                                                 
+        for j in eachnode(dg), i in eachnode(dg)                                                    
+            u_node = get_node_vars(u, equations, dg, i, j, element)                                 
+
+            # compute compressible Euler fluxes
+            u_cons = SVector(u_node[1], u_node[2], u_node[3], u_node[4])
+            eq_relax = equations.equations_relaxation
+            vu = flux(u_cons,1,eq_relax)
+            wu = flux(u_cons,2,eq_relax)
+            # compute relaxation terms
+            du1 = 0.0
+            du2 = 0.0
+            du3 = 0.0
+            du4 = 0.0
+            du5 = factor * (eps * u[5] + dt_ * vu[1])
+            du6 = factor * (eps * u[6] + dt_ * vu[2])
+            du7 = factor * (eps * u[7] + dt_ * vu[3])
+            du8 = factor * (eps * u[8] + dt_ * vu[4])
+            du9 = factor * (eps * u[9] + dt_ * wu[1])
+            du10= factor * (eps * u[10]+ dt_ * wu[2])
+            du11= factor * (eps * u[11]+ dt_ * wu[3])
+            du12= factor * (eps * u[12]+ dt_ * wu[4])
+            new_u = SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, du10, du11, du12)
+            set_node_vars!(u, new_u, equations, dg, i, j, element)                                            
+        end                                                                                         
+    end                                                                                             
+
+    return nothing                                                                                  
+end       
 end # @muladd

src/callbacks_step/lbm_collision_dg2d.jl

@@ -5,6 +5,18 @@
 @muladd begin
 #! format: noindent
 
+#function apply_collision!(u, dt, collision_op,
+#                          mesh::AbstractMesh{2}, equations, dg::DG, cache)
+#    @threaded for element in eachelement(dg, cache)
+#        for j in eachnode(dg), i in eachnode(dg)
+#            u_node = get_node_vars(u, equations, dg, i, j, element)
+#            update = collision_op(u_node, dt, equations)
+#            add_to_node_vars!(u, update, equations, dg, i, j, element)
+#        end
+#    end
+#
+#    return nothing
+#end
 function apply_collision!(u, dt, collision_op,
                           mesh::AbstractMesh{2}, equations, dg::DG, cache)
     @threaded for element in eachelement(dg, cache)

src/equations/equations.jl

@@ -450,4 +450,9 @@ include("linearized_euler_2d.jl")
 
 abstract type AbstractEquationsParabolic{NDIMS, NVARS} <:
               AbstractEquations{NDIMS, NVARS} end
+
+# Jin Xin Relaxation type equation (Compressible Euler variant)
+abstract type AbstractJinXinEquations{NDIMS, NVARS} <:
+              AbstractEquations{NDIMS, NVARS} end
+include("jin_xin_compressible_euler_2d.jl")
 end # @muladd