From 5f4fdc7e24640d73d10ccbb763d120b6a038ae3a Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 28 May 2024 17:17:52 +0200 Subject: [PATCH] draft a generalized version of Jin Xin --- .../tree_2d_dgsem/elixir_jin_xin_euler.jl | 11 +- src/Trixi.jl | 4 +- .../positivity_zhang_shu_dg2d.jl | 12 +- .../subcell_limiter_idp_correction_2d.jl | 24 +- .../jin_xin_compressible_euler_2d.jl | 300 +++++++++--------- 5 files changed, 176 insertions(+), 175 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index 6ee8d11df98..2272713a7c8 100644 --- 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 index d0ee4b6f2a0..ce53173fc9e 100644 --- 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 index e56b900131a..fc17bea620f 100644 --- 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 index a99fd0c4670..e3ddf36ae74 100644 --- 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 diff --git a/src/equations/jin_xin_compressible_euler_2d.jl b/src/equations/jin_xin_compressible_euler_2d.jl index 8472e295005..edbc202291e 100644 --- a/src/equations/jin_xin_compressible_euler_2d.jl +++ b/src/equations/jin_xin_compressible_euler_2d.jl @@ -6,190 +6,192 @@ #! format: noindent @doc raw""" - JinXinCompressibleEulerEquations2D{RealT <: Real} <: - AbstractJinXinEquations{2, 12} + JinXinEquations + +TODO: Write a proper docstring """ -struct JinXinCompressibleEulerEquations2D{RealT <: Real,CollisionOp} <: - AbstractJinXinEquations{2, 12} - eps_relaxation ::RealT # relaxation parameter epsilon with the tendency epsilon -> 0 - eps_relaxation_inv::RealT # relaxation parameter epsilon with the tendency epsilon -> 0 - - a::SVector{4, RealT} # discrete molecular velocity components in x-direction - sqrt_a::SVector{4, RealT} # discrete molecular velocity components in x-direction - sqrt_a_inv::SVector{4, RealT} # discrete molecular velocity components in x-direction - b::SVector{4, RealT} # discrete molecular velocity components in y-direction - sqrt_b::SVector{4, RealT} # discrete molecular velocity components in x-direction - sqrt_b_inv::SVector{4, RealT} # discrete molecular velocity components in x-direction - equations_relaxation::CompressibleEulerEquations2D{RealT} - collision_op::CollisionOp # collision operator for the collision kernel +struct JinXinEquations{NDIMS, NVARS, NVARS_BASE, EquationsBase <: AbstractEquations{NDIMS, NVARS_BASE}, RealT <: Real} <: + AbstractJinXinEquations{NDIMS, NVARS} + equations_base::EquationsBase + + # relaxation parameter of the Jin Xin relaxation model + # The relaxed equations should converge to the original equations + # as the relaxation parameter epsilon -> 0 + eps_relaxation ::RealT + eps_relaxation_inv::RealT + + # velocity parameters of the Jin Xin relaxation model + # They need to satisfy the subcharacteristic condition, i.e., they + # must be upper bounds of the characteristic speeds of the original + # equations + velocities::NTuple{NDIMS, SVector{NVARS_BASE, RealT}} + sqrt_velocities::NTuple{NDIMS, SVector{NVARS_BASE, RealT}} + sqrt_velocities_inv::NTuple{NDIMS, SVector{NVARS_BASE, RealT}} end -function JinXinCompressibleEulerEquations2D(eps_relaxation,a1,a2,a3,a4,b1,b2,b3,b4,equations_relaxation) - - a = SVector(a1,a2,a3,a4) - sa = SVector(sqrt(a1),sqrt(a2),sqrt(a3),sqrt(a4)) - sai = SVector(1.0/sqrt(a1),1.0/sqrt(a2),1.0/sqrt(a3),1.0/sqrt(a4)) - - b = SVector(b1,b2,b3,b4) - sb = SVector(sqrt(b1),sqrt(b2),sqrt(b3),sqrt(b4)) - sbi = SVector(1.0/sqrt(b1),1.0/sqrt(b2),1.0/sqrt(b3),1.0/sqrt(b4)) +function JinXinEquations(equations_base, eps_relaxation, velocities) - collision_op = collision_bgk + sqrt_velocities = map(velocities) do v + sqrt_v = sqrt.(v) + return sqrt_v + end + sqrt_velocities_inv = map(sqrt_velocities) do sv + sqrt_vi = inv.(sv) + return sqrt_vi + end - JinXinCompressibleEulerEquations2D(eps_relaxation, 1.0/eps_relaxation, a, sa, sai, b, sb, sbi, equations_relaxation, collision_op) + NDIMS = ndims(equations_base) + NVARS_BASE = nvariables(equations_base) + RealT = promote_type(typeof(eps_relaxation), eltype(eltype(velocities))) + JinXinEquations{NDIMS, (NDIMS + 1) * NVARS_BASE, NVARS_BASE, typeof(equations_base), RealT}(equations_base, eps_relaxation, inv(eps_relaxation), + velocities, sqrt_velocities, sqrt_velocities_inv) end -function varnames(::typeof(cons2cons), equations::JinXinCompressibleEulerEquations2D) - ("rho", "rho_v1", "rho_v2", "rho_e", "vu1", "vu2", "vu3", "vu4", "wu1", "wu2", "wu3", "wu4") +function varnames(func::typeof(cons2cons), equations::JinXinEquations) + basenames = varnames(func, equations.equations_base) + + if ndims(equations) == 1 + velocities_1 = ntuple(n -> "v1_" * string(n), + Val(nvariables(equations.equations_base))) + return (basenames..., velocities_1...) + elseif ndims(equations) == 2 + velocities_1 = ntuple(n -> "v1_" * string(n), + Val(nvariables(equations.equations_base))) + velocities_2 = ntuple(n -> "v2_" * string(n), + Val(nvariables(equations.equations_base))) + return (basenames..., velocities_1..., velocities_2...) + elseif ndims(equations) == 3 + velocities_1 = ntuple(n -> "v1_" * string(n), + Val(nvariables(equations.equations_base))) + velocities_2 = ntuple(n -> "v2_" * string(n), + Val(nvariables(equations.equations_base))) + velocities_3 = ntuple(n -> "v3_" * string(n), + Val(nvariables(equations.equations_base))) + return (basenames..., velocities_1..., velocities_2..., velocities_3...) + else + throw(ArgumentError("Number of dimensions $(ndims(equations)) not supported")) + end end -function varnames(::typeof(cons2prim), equations::JinXinCompressibleEulerEquations2D) + +function varnames(::typeof(cons2prim), equations::JinXinEquations) + # TODO: Jin Xin varnames(cons2cons, equations) end # Set initial conditions at physical location `x` for time `t` -""" - initial_condition_constant(x, t, equations::JinXinCompressibleEulerEquations2D) - -A constant initial condition to test free-stream preservation. -""" -function initial_condition_constant(x, t, equations::JinXinCompressibleEulerEquations2D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 - rho_e = 10.0 - - u_cons = (rho, rho_v1, rho_v2, rho_e) - eq_relax = equations.equations_relaxation - vu = flux(u_cons,1,eq_relax) - wu = flux(u_cons,2,eq_relax) - - return SVector(rho, rho_v1, rho_v2, rho_e, vu[1], vu[2], vu[3], vu[4], wu[1], wu[2], wu[3], wu[4]) -end - struct InitialConditionJinXin{IC} initial_condition::IC end -@inline function (ic::InitialConditionJinXin)(x,t,equations) - - eq_relax = equations.equations_relaxation - u = ic.initial_condition(x,t,eq_relax) - vu = flux(u,1,eq_relax) - wu = flux(u,2,eq_relax) - - return SVector(u[1], u[2], u[3], u[4], vu[1], vu[2], vu[3], vu[4], wu[1], wu[2], wu[3], wu[4]) +@inline function (ic::InitialConditionJinXin)(x, t, equations::JinXinEquations) + eq_base = equations.equations_base + u = ic.initial_condition(x, t, eq_base) + + if ndims(equations) == 1 + v1 = flux(u, 1, eq_base) + return SVector(u..., v1...) + elseif ndims(equations) == 2 + v1 = flux(u, 1, eq_base) + v2 = flux(u, 2, eq_base) + return SVector(u..., v1..., v2...) + else + v1 = flux(u, 1, eq_base) + v2 = flux(u, 2, eq_base) + v3 = flux(u, 3, eq_base) + return SVector(u..., v1..., v2..., v3...) + end end # Pre-defined source terms should be implemented as -# function source_terms_WHATEVER(u, x, t, equations::JinXinCompressibleEulerEquations2D) -@inline function source_terms_JinXin_Relaxation(u, x, t, - equations::JinXinCompressibleEulerEquations2D) - - # relaxation parameter - eps_inv= equations.eps_relaxation_inv - - # compute compressible Euler fluxes - u_cons = SVector(u[1], u[2], u[3], u[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 = -eps_inv * (u[5] - vu[1]) - du6 = -eps_inv * (u[6] - vu[2]) - du7 = -eps_inv * (u[7] - vu[3]) - du8 = -eps_inv * (u[8] - vu[4]) - du9 = -eps_inv * (u[9] - wu[1]) - du10= -eps_inv * (u[10]- wu[2]) - du11= -eps_inv * (u[11]- wu[3]) - du12= -eps_inv * (u[12]- wu[4]) - - return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, du10, du11, du12) -end +# @inline function source_terms_JinXin_Relaxation(u, x, t, +# equations::JinXinEquations) + +# # relaxation parameter +# eps_inv = equations.eps_relaxation_inv + +# # compute compressible Euler fluxes +# u_cons = SVector(u[1], u[2], u[3], u[4]) +# eq_relax = equations.equations_base +# 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 = -eps_inv * (u[5] - vu[1]) +# du6 = -eps_inv * (u[6] - vu[2]) +# du7 = -eps_inv * (u[7] - vu[3]) +# du8 = -eps_inv * (u[8] - vu[4]) +# du9 = -eps_inv * (u[9] - wu[1]) +# du10= -eps_inv * (u[10]- wu[2]) +# du11= -eps_inv * (u[11]- wu[3]) +# du12= -eps_inv * (u[12]- wu[4]) + +# return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, du10, du11, du12) +# end + + +function get_block_components(u, n, equations::JinXinEquations) + nvars_base = nvariables(equations.equations_base) + return SVector(ntuple(i -> u[i + (n - 1) * nvars_base], Val(nvars_base))) +end # Calculate 1D flux in for a single point -@inline function flux(u, orientation::Integer, equations::JinXinCompressibleEulerEquations2D) +# TODO: Implement 1D and 3D +@inline function flux(u, orientation::Integer, equations::JinXinEquations{2}) if orientation == 1 - a = equations.a - flux = SVector(u[5],u[6],u[7],u[8],a[1]*u[1],a[2]*u[2],a[3]*u[3],a[4]*u[4],0,0,0,0) - else - b = equations.b - flux = SVector(u[9],u[10],u[11],u[12],0,0,0,0,b[1]*u[1],b[2]*u[2],b[3]*u[3],b[4]*u[4]) + u_base = get_block_components(u, 1, equations) + fluxes = get_block_components(u, 2, equations) + velocities = equations.velocities[1] + return SVector(fluxes..., (velocities .* u_base)..., zero(u_base)...) + else # orientation == 2 + u_base = get_block_components(u, 1, equations) + fluxes = get_block_components(u, 3, equations) + velocities = equations.velocities[2] + return SVector(fluxes..., zero(u_base)..., (velocities .* u_base)...) end - return flux end - +# TODO: Implement 1D and 3D @inline function flux_upwind(u_ll, u_rr, orientation::Integer, - equations::JinXinCompressibleEulerEquations2D) + equations::JinXinEquations{2}) + u_ll_base = get_block_components(u_ll, 1, equations) + u_rr_base = get_block_components(u_rr, 1, equations) + if orientation == 1 - sai = equations.sqrt_a_inv - sa = equations.sqrt_a - #dissipation = SVector(sai[1]*(u_rr[5] - u_ll[5]), sai[2]*(u_rr[6] - u_ll[6]), sai[3]*(u_rr[7] - u_ll[7]), sai[4]*(u_rr[8] - u_ll[8]), sa[1]*(u_rr[1] - u_ll[1]), sa[2]*(u_rr[2] - u_ll[2]), sa[3]*(u_rr[3] - u_ll[3]), sa[4]*(u_rr[4] - u_ll[4]),0.0,0.0,0.0,0.0) - dissipation = SVector(sa[1]*(u_rr[1] - u_ll[1]), sa[2]*(u_rr[2] - u_ll[2]), sa[3]*(u_rr[3] - u_ll[3]), sa[4]*(u_rr[4] - u_ll[4]), sa[1]*(u_rr[5] - u_ll[5]), sa[2]*(u_rr[6] - u_ll[6]), sa[3]*(u_rr[7] - u_ll[7]), sa[4]*(u_rr[8] - u_ll[8]),0.0,0.0,0.0,0.0) - else - sbi = equations.sqrt_b_inv - sb = equations.sqrt_b - #dissipation = SVector(sbi[1]*(u_rr[9] - u_ll[9]), sbi[2]*(u_rr[10] - u_ll[10]), sbi[3]*(u_rr[11] - u_ll[11]), sbi[4]*(u_rr[12] - u_ll[12]),0.0,0.0,0.0,0.0, sb[1]*(u_rr[1] - u_ll[1]), sb[2]*(u_rr[2] - u_ll[2]), sb[3]*(u_rr[3] - u_ll[3]), sb[4]*(u_rr[4] - u_ll[4])) - dissipation = SVector(sb[1]*(u_rr[1] - u_ll[1]), sb[2]*(u_rr[2] - u_ll[2]), sb[3]*(u_rr[3] - u_ll[3]), sb[4]*(u_rr[4] - u_ll[4]),0.0,0.0,0.0,0.0,sb[1]*(u_rr[9] - u_ll[9]), sb[2]*(u_rr[10] - u_ll[10]), sb[3]*(u_rr[11] - u_ll[11]), sb[4]*(u_rr[12] - u_ll[12])) + sqrt_velocities = equations.sqrt_velocities[1] + f_ll_base = get_block_components(u_ll, 2, equations) + f_rr_base = get_block_components(u_rr, 2, equations) + dissipation = SVector((sqrt_velocities .* (u_rr_base - u_ll_base))..., + (sqrt_velocities .* (f_rr_base + f_ll_base))..., + zero(u_ll_base)...) + else # orientation == 2 + sqrt_velocities = equations.sqrt_velocities[2] + f_ll_base = get_block_components(u_ll, 3, equations) + f_rr_base = get_block_components(u_rr, 3, equations) + dissipation = SVector((sqrt_velocities .* (u_rr_base - u_ll_base))..., + zero(u_ll_base)..., + (sqrt_velocities .* (f_rr_base + f_ll_base))...) end - return 0.5 * (flux(u_ll,orientation,equations) + flux(u_rr,orientation,equations) - dissipation) + + return 0.5f0 * (flux(u_ll, orientation, equations) + + flux(u_rr, orientation, equations) - dissipation) end -""" - collision_bgk(u, dt, equations::LatticeBoltzmannEquations2D) - -Collision operator for the Bhatnagar, Gross, and Krook (BGK) model. -""" -@inline function collision_bgk(u, dt, equations::JinXinCompressibleEulerEquations2D) - # relaxation parameter - eps = equations.eps_relaxation - - # compute compressible Euler fluxes - u_cons = SVector(u[1], u[2], u[3], u[4]) - eq_relax = equations.equations_relaxation - vu = flux(u_cons,1,eq_relax) - wu = flux(u_cons,2,eq_relax) - - dt_ = dt * 1.0 - factor =1.0/ (eps + dt_) - - # compute relaxation terms - du1 = 0.0 - du2 = 0.0 - du3 = 0.0 - du4 = 0.0 - du5 = -u[5] + factor * (eps * u[5] + dt_ * vu[1]) - du6 = -u[6] + factor * (eps * u[6] + dt_ * vu[2]) - du7 = -u[7] + factor * (eps * u[7] + dt_ * vu[3]) - du8 = -u[8] + factor * (eps * u[8] + dt_ * vu[4]) - du9 = -u[9] + factor * (eps * u[9] + dt_ * wu[1]) - du10= -u[10]+ factor * (eps * u[10]+ dt_ * wu[2]) - du11= -u[11]+ factor * (eps * u[11]+ dt_ * wu[3]) - du12= -u[12]+ factor * (eps * u[12]+ dt_ * wu[4]) - - return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, du10, du11, du12) -end - - - -@inline function max_abs_speeds(u,equations::JinXinCompressibleEulerEquations2D) - sa = equations.sqrt_a - sb = equations.sqrt_b - - return max(sa[1],sa[2],sa[3],sa[4]), max(sb[1],sb[2],sb[3],sb[4]) +@inline function max_abs_speeds(u, equations::JinXinEquations{2}) + return ntuple(Val(ndims(equations))) do n + maximum(equations.sqrt_velocities_inv[n]) + end end -# not correct yet!! +# TODO: not correct yet!! # Convert conservative variables to primitive -@inline cons2prim(u, equations::JinXinCompressibleEulerEquations2D) = u +@inline cons2prim(u, equations::JinXinEquations) = u # Convert conservative variables to entropy variables -@inline cons2entropy(u, equations::JinXinCompressibleEulerEquations2D) = u +@inline cons2entropy(u, equations::JinXinEquations) = u end # @muladd