draft a generalized version of Jin Xin

trixi-framework · May 28, 2024 · 5f4fdc7 · 5f4fdc7
1 parent 669dcaf
commit 5f4fdc7
Show file tree

Hide file tree

Showing 5 changed files with 176 additions and 175 deletions.
diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl
@@ -9,8 +9,9 @@ epsilon_relaxation = 5.0e-6
 a1 = a2 = a3 = a4 = 30.0
 b1 = b2 = b3 = b4 = 30.0
 
-equations_relaxation = CompressibleEulerEquations2D(1.4)
-equations = JinXinCompressibleEulerEquations2D(epsilon_relaxation, a1, a2, a3, a4, b1, b2, b3, b4,equations_relaxation)
+equations_base = CompressibleEulerEquations2D(1.4)
+velocities = (SVector(a1, a2, a3, a4), SVector(b1, b2, b3, b4))
+equations = JinXinEquations(equations_base, epsilon_relaxation, velocities)
 
 function initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
   # change discontinuity to tanh
@@ -24,7 +25,7 @@ function initial_condition_kelvin_helmholtz_instability(x, t, equations::Compres
   v2 = 0.1 * sin(2 * pi * x[1])
   p = 1.0
   return prim2cons(SVector(rho, v1, v2, p), equations)
-end 
+end
 
 #initial_condition = initial_condition_constant
 initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability)
@@ -76,11 +77,9 @@ save_solution = SaveSolutionCallback(interval=1000,
 
 stepsize_callback = StepsizeCallback(cfl=0.5)
 
-collision_callback = LBMCollisionCallback()  
-
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,stepsize_callback)#,collision_callback)
+                        save_solution,stepsize_callback)
 
 ###############################################################################
 # run the simulation

diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -163,7 +163,7 @@ export AcousticPerturbationEquations2D,
        LinearizedEulerEquations1D, LinearizedEulerEquations2D, LinearizedEulerEquations3D,
        PolytropicEulerEquations2D,
        TrafficFlowLWREquations1D,
-       JinXinCompressibleEulerEquations2D
+       JinXinEquations
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
@@ -210,7 +210,7 @@ export boundary_condition_do_nothing,
        BoundaryConditionCoupled
 
 export initial_condition_convergence_test, source_terms_convergence_test
-export source_terms_harmonic, source_terms_JinXin_Relaxation
+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,

diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl
@@ -47,13 +47,13 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
 
     return nothing
 end
-function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompressibleEulerEquations2D, dg, cache)
+function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinEquations, dg, cache)
 
             # relaxation parameter
             eps = equations.eps_relaxation
             dt_ = dt
             factor =1.0/ (eps + dt_)
-            eq_relax = equations.equations_relaxation
+            eq_relax = equations.equations_base
 
             # prepare local storage for projection
             @unpack interpolate_N_to_M, project_M_to_N = dg.basis
@@ -71,7 +71,7 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr
 
  #@threaded for element in eachelement(dg, cache)
   for element in eachelement(dg, cache)
-       
+
         # get element u_N
         for j in eachnode(dg), i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, j, element)
@@ -81,11 +81,11 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr
         end
         # bring elemtn u_N to grid (M+1)x(M+1)
         multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N)
-        
+
         # compute nodal values of entropy variables w on the M grid
         for j in 1:nnodes_projection, i in 1:nnodes_projection
             u_cons = get_node_vars(u_M, eq_relax, dg, i, j)
-            w_ij   = cons2entropy(u_cons,eq_relax) 
+            w_ij   = cons2entropy(u_cons,eq_relax)
             set_node_vars!(w_M,w_ij,eq_relax,dg,i,j)
         end
         # compute projection of w with M values down to N
@@ -105,7 +105,7 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr
         # compute projection of f with M values down to N, same for g
         multiply_dimensionwise!(f_N,project_M_to_N,f_M)
         multiply_dimensionwise!(g_N,project_M_to_N,g_M)
-        #@assert nnodes_projection == nnodes(dg) 
+        #@assert nnodes_projection == nnodes(dg)
         #for j in 1:nnodes_projection, i in 1:nnodes_projection
         #    u_cons = get_node_vars(u_N, eq_relax, dg, i, j)
         #    f_cons = flux(u_cons,1,eq_relax)

diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl
@@ -44,20 +44,20 @@ function perform_idp_correction!(u, dt,
 
     return nothing
 end
-#function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompressibleEulerEquations2D, dg, cache)
+#function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinEquations, 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)                                 
-#            
+# @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
+#            eq_relax = equations.equations_base
 #            vu = flux(u_cons,1,eq_relax)
 #            wu = flux(u_cons,2,eq_relax)
 #            # compute relaxation terms
@@ -74,10 +74,10 @@ end
 #            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       
+#            set_node_vars!(u, new_u, equations, dg, i, j, element)
+#        end
+#    end
+#
+#    return nothing
+#end
 end # @muladd