From 7548300cb8c629ff8ca44a0ed5c177496781c2bf Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Mon, 9 Oct 2023 09:34:16 +0200 Subject: [PATCH 01/20] dirty implementation of Jin Xin Relaxation for Compressible Euler 2D --- ...ixir_euler_kelvin_helmholtz_instability.jl | 6 +- .../tree_2d_dgsem/elixir_jin_xin_euler.jl | 93 +++++++++ src/Trixi.jl | 5 +- src/callbacks_stage/positivity_zhang_shu.jl | 5 +- .../positivity_zhang_shu_dg2d.jl | 36 ++++ .../subcell_limiter_idp_correction_2d.jl | 36 ++++ src/callbacks_step/lbm_collision_dg2d.jl | 12 ++ src/equations/equations.jl | 5 + .../jin_xin_compressible_euler_2d.jl | 195 ++++++++++++++++++ 9 files changed, 386 insertions(+), 7 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_jin_xin_euler.jl create mode 100644 src/equations/jin_xin_compressible_euler_2d.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl index 4fdedf516ef..cbae69163d3 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl +++ b/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, diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl new file mode 100644 index 00000000000..af5519b9329 --- /dev/null +++ b/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 diff --git a/src/Trixi.jl b/src/Trixi.jl index b65d03e7975..feed878b075 100644 --- a/src/Trixi.jl +++ b/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, diff --git a/src/callbacks_stage/positivity_zhang_shu.jl b/src/callbacks_stage/positivity_zhang_shu.jl index 92141c4b26e..5cd2bf4ea34 100644 --- a/src/callbacks_stage/positivity_zhang_shu.jl +++ b/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 diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index b37ed9c49d5..e48a22a1452 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/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 diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index f6b91444578..c9d636d9a4a 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/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 diff --git a/src/callbacks_step/lbm_collision_dg2d.jl b/src/callbacks_step/lbm_collision_dg2d.jl index 932edfd61f6..887c4ca12f2 100644 --- a/src/callbacks_step/lbm_collision_dg2d.jl +++ b/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) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 9bae563d85f..bdfe1e54661 100644 --- a/src/equations/equations.jl +++ b/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 diff --git a/src/equations/jin_xin_compressible_euler_2d.jl b/src/equations/jin_xin_compressible_euler_2d.jl new file mode 100644 index 00000000000..8472e295005 --- /dev/null +++ b/src/equations/jin_xin_compressible_euler_2d.jl @@ -0,0 +1,195 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + JinXinCompressibleEulerEquations2D{RealT <: Real} <: + AbstractJinXinEquations{2, 12} +""" +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 +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)) + + collision_op = collision_bgk + + JinXinCompressibleEulerEquations2D(eps_relaxation, 1.0/eps_relaxation, a, sa, sai, b, sb, sbi, equations_relaxation, collision_op) +end + +function varnames(::typeof(cons2cons), equations::JinXinCompressibleEulerEquations2D) + ("rho", "rho_v1", "rho_v2", "rho_e", "vu1", "vu2", "vu3", "vu4", "wu1", "wu2", "wu3", "wu4") +end +function varnames(::typeof(cons2prim), equations::JinXinCompressibleEulerEquations2D) + 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]) +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 + + +# Calculate 1D flux in for a single point +@inline function flux(u, orientation::Integer, equations::JinXinCompressibleEulerEquations2D) + 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]) + end + return flux +end + + +@inline function flux_upwind(u_ll, u_rr, orientation::Integer, + equations::JinXinCompressibleEulerEquations2D) + 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])) + end + return 0.5 * (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]) +end + +# not correct yet!! +# Convert conservative variables to primitive +@inline cons2prim(u, equations::JinXinCompressibleEulerEquations2D) = u + +# Convert conservative variables to entropy variables +@inline cons2entropy(u, equations::JinXinCompressibleEulerEquations2D) = u +end # @muladd From 86a063ad87c945bcc08f82cba51bcbdc5647c5cf Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 10 Oct 2023 13:39:12 +0200 Subject: [PATCH 02/20] fixed some double definition and reduced resolution in elixir --- .../tree_2d_dgsem/elixir_jin_xin_euler.jl | 8 ++- .../subcell_limiter_idp_correction_2d.jl | 72 +++++++++---------- 2 files changed, 41 insertions(+), 39 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index af5519b9329..a4234b596d6 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -6,8 +6,8 @@ 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 +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) @@ -28,6 +28,7 @@ end #initial_condition = initial_condition_constant initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability) +#initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave) solver = DGSEM(polydeg=3, surface_flux=Trixi.flux_upwind) #surface_flux = Trixi.flux_upwind @@ -45,7 +46,7 @@ coordinates_min = (-1.0, -1.0) coordinates_max = ( 1.0, 1.0) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=6, + initial_refinement_level=5, n_cells_max=1_000_000) @@ -56,6 +57,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)# # ODE solvers, callbacks etc. tspan = (0.0, 3.7) +#tspan = (0.0, 2.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl index c9d636d9a4a..1c1133b0138 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl @@ -41,40 +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 +#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 From 48391591c06372caa6cca664d3649a9b350d9593 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Wed, 17 Apr 2024 14:10:07 +0200 Subject: [PATCH 03/20] small changes --- .../elixir_euler_kelvin_helmholtz_instability.jl | 12 ++++++------ examples/tree_2d_dgsem/elixir_jin_xin_euler.jl | 10 ++++------ 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl index cbae69163d3..5e6194f5be4 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl @@ -32,12 +32,12 @@ end initial_condition = initial_condition_kelvin_helmholtz_instability surface_flux = flux_lax_friedrichs -volume_flux = flux_ranocha +volume_flux = flux_central polydeg = 3 basis = LobattoLegendreBasis(polydeg) indicator_sc = IndicatorHennemannGassner(equations, basis, - alpha_max=0.002, - alpha_min=0.0001, + alpha_max=0.000, + alpha_min=0.0000, alpha_smooth=true, variable=density_pressure) volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; @@ -48,7 +48,7 @@ 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=4, + initial_refinement_level=6, n_cells_max=100_000) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) @@ -60,12 +60,12 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval) alive_callback = AliveCallback(analysis_interval=analysis_interval) -save_solution = SaveSolutionCallback(interval=20, +save_solution = SaveSolutionCallback(interval=1000, save_initial_solution=true, save_final_solution=true, solution_variables=cons2prim) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index a4234b596d6..65676e25a15 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -5,7 +5,7 @@ using Trixi ############################################################################### # semidiscretization of the compressible Euler equations -epsilon_relaxation = 1.0e-5 +epsilon_relaxation = 5.0e-6 a1 = a2 = a3 = a4 = 30.0 b1 = b2 = b3 = b4 = 30.0 @@ -46,18 +46,16 @@ 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=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) -#tspan = (0.0, 2.0) +#tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -72,7 +70,7 @@ save_solution = SaveSolutionCallback(interval=1000, save_final_solution=true, solution_variables=cons2prim) -stepsize_callback = StepsizeCallback(cfl=0.1) +stepsize_callback = StepsizeCallback(cfl=0.5) collision_callback = LBMCollisionCallback() From c4b4803dea83501ca987e27fe218af957ef91295 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 28 May 2024 10:13:19 +0100 Subject: [PATCH 04/20] working towards projection --- src/auxiliary/precompile.jl | 6 +- .../positivity_zhang_shu_dg2d.jl | 30 +++++++++ src/solvers/dgsem/basis_lobatto_legendre.jl | 61 +++++++++++++++---- 3 files changed, 82 insertions(+), 15 deletions(-) diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 4d5399b5ba3..70ebf923b87 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -186,6 +186,8 @@ function _precompile_manual_() Matrix{RealT}, # DerivativeMatrix #StaticArrays.SArray{Tuple{nnodes_,nnodes_},RealT,2,nnodes_^2}, + Matrix{RealT}, + Matrix{RealT}, Matrix{RealT}} end @@ -302,9 +304,9 @@ function _precompile_manual_() # Constructors of the basis are inherently type-unstable since we pass integers # and use their values as parameters of static arrays. # Nevertheless, we can still precompile methods used to construct the bases. - Base.precompile(Tuple{Type{LobattoLegendreBasis}, Int}) + #Base.precompile(Tuple{Type{LobattoLegendreBasis}, Int}) for RealT in (Float64,) - Base.precompile(Tuple{Type{LobattoLegendreBasis}, RealT, Int}) + # Base.precompile(Tuple{Type{LobattoLegendreBasis}, RealT, Int}) @assert Base.precompile(Tuple{typeof(Trixi.calc_dhat), Vector{RealT}, Vector{RealT}}) @assert Base.precompile(Tuple{typeof(Trixi.calc_dsplit), Vector{RealT}, diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index e5a656f8b5d..ca97bd9cc33 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -55,12 +55,42 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr factor =1.0/ (eps + dt_) eq_relax = equations.equations_relaxation + # prepare local storage for projection + @unpack interpolate_N_to_M, project_M_to_N = dg.basis + nnodes_,nnodes_projection = size(project_M_to_N) + nVars = nvariables(eq_relax) + RealT = real(dg) + u_N = zeros(RealT, nVars, nnodes_, nnodes_) + f_N = zeros(RealT, nVars, nnodes_, nnodes_) + g_N = zeros(RealT, nVars, nnodes_, nnodes_) + u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + @threaded 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) + for v in eachvariable(eq_relax) + u_N[v,i,j] = u_node[v] + end + end + # compute projection of f,g + for j in eachnode(dg), i in eachnode(dg) + u_cons = u_N[:,i,j] + f_N[:,i,j] = flux(u_cons,1,eq_relax) + g_N[:,i,j] = flux(u_cons,2,eq_relax) + end + + 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 = f_N[:,i,j] #flux(u_cons,1,eq_relax) + #wu = g_N[:,i,j] #flux(u_cons,2,eq_relax) vu = flux(u_cons,1,eq_relax) wu = flux(u_cons,2,eq_relax) # compute relaxation terms diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index cac1dba9c74..d1c9b909c6e 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -18,7 +18,9 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, VectorT <: AbstractVector{RealT}, InverseVandermondeLegendre <: AbstractMatrix{RealT}, BoundaryMatrix <: AbstractMatrix{RealT}, - DerivativeMatrix <: AbstractMatrix{RealT}} <: + DerivativeMatrix <: AbstractMatrix{RealT}, + Matrix_MxN <: AbstractMatrix{RealT}, + Matrix_NxM <: AbstractMatrix{RealT}} <: AbstractBasisSBP{RealT} nodes::VectorT weights::VectorT @@ -32,9 +34,12 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, derivative_split_transpose::DerivativeMatrix # transpose of `derivative_split` derivative_dhat::DerivativeMatrix # weak form matrix "dhat", # negative adjoint wrt the SBP dot product + # L2 projection operators + interpolate_N_to_M::Matrix_MxN # interpolates from N nodes to M nodes + project_M_to_N::Matrix_NxM # compute projection via Legendre modes, cut off modes at N, back to N nodes end -function LobattoLegendreBasis(RealT, polydeg::Integer) +function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) nnodes_ = polydeg + 1 # compute everything using `Float64` by default @@ -67,20 +72,33 @@ function LobattoLegendreBasis(RealT, polydeg::Integer) derivative_split_transpose = Matrix{RealT}(derivative_split_transpose_) derivative_dhat = Matrix{RealT}(derivative_dhat_) + # L2 projection operators + nnodes_projection = polydeg_projection + 1 + nodes_projection, weights_projection = gauss_lobatto_nodes_weights(nnodes_projection) + interpolate_N_to_M_ = polynomial_interpolation_matrix(nodes_, nodes_projection) + interpolate_N_to_M = Matrix{RealT}(interpolate_N_to_M_) + + project_M_to_N_ = calc_projection_matrix(nodes_projection, nodes_) + project_M_to_N = Matrix{RealT}(project_M_to_N_) + return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes), typeof(inverse_vandermonde_legendre), typeof(boundary_interpolation), - typeof(derivative_matrix)}(nodes, weights, - inverse_weights, - inverse_vandermonde_legendre, - boundary_interpolation, - derivative_matrix, - derivative_split, - derivative_split_transpose, - derivative_dhat) -end - -LobattoLegendreBasis(polydeg::Integer) = LobattoLegendreBasis(Float64, polydeg) + typeof(derivative_matrix), + typeof(interpolate_N_to_M), + typeof(project_M_to_N)}(nodes, weights, + inverse_weights, + inverse_vandermonde_legendre, + boundary_interpolation, + derivative_matrix, + derivative_split, + derivative_split_transpose, + derivative_dhat, + interpolate_N_to_M, + project_M_to_N) +end + +LobattoLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) = LobattoLegendreBasis(Float64, polydeg; polydeg_projection) function Base.show(io::IO, basis::LobattoLegendreBasis) @nospecialize basis # reduce precompilation time @@ -780,4 +798,21 @@ function vandermonde_legendre(nodes, N) return vandermonde, inverse_vandermonde end vandermonde_legendre(nodes) = vandermonde_legendre(nodes, length(nodes) - 1) + +function calc_projection_matrix(nodes_in,nodes_out) + # nodes_in are size M>N + nnodes_in=length(nodes_in) + polydeg_in = nnodes_in - 1 + # nodes_out are size N + nnodes_out = length(nodes_out) + vandermonde_in,inverse_vandermonde_in = vandermonde_legendre(nodes_in,polydeg_in) + filter_matrix = zeros(nnodes_in,nnodes_in) + for j in 1:nnodes_out + filter_matrix[j,j] = 1 + end + interpolate_M_to_N = polynomial_interpolation_matrix(nodes_in, nodes_out) + projection_matrix = interpolate_M_to_N * vandermonde_in * filter_matrix * inverse_vandermonde_in + return projection_matrix +end + end # @muladd From 8941d537527e2c2e5edffaa8146858b2c6ff9167 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 28 May 2024 12:20:48 +0100 Subject: [PATCH 05/20] standard L2 projection implemented --- .../tree_2d_dgsem/elixir_jin_xin_euler.jl | 7 +++- .../positivity_zhang_shu_dg2d.jl | 41 ++++++++++++------- 2 files changed, 31 insertions(+), 17 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index 65676e25a15..a00a13f79b1 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -5,7 +5,7 @@ using Trixi ############################################################################### # semidiscretization of the compressible Euler equations -epsilon_relaxation = 5.0e-6 +epsilon_relaxation = 5.0e-7 a1 = a2 = a3 = a4 = 30.0 b1 = b2 = b3 = b4 = 30.0 @@ -29,7 +29,10 @@ end #initial_condition = initial_condition_constant initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability) #initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave) -solver = DGSEM(polydeg=3, surface_flux=Trixi.flux_upwind) +polydeg = 3 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) +solver = DGSEM(basis, Trixi.flux_upwind) +#solver = DGSEM(polydeg = 3, surface_flux = Trixi.flux_upwind) #surface_flux = Trixi.flux_upwind #volume_flux = flux_central diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index ca97bd9cc33..59f2f4b0f8f 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -67,32 +67,43 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - @threaded for element in eachelement(dg, cache) - + #@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) for v in eachvariable(eq_relax) - u_N[v,i,j] = u_node[v] + u_N[v,i,j] = u_node[v] end end - - # compute projection of f,g - for j in eachnode(dg), i in eachnode(dg) - u_cons = u_N[:,i,j] - f_N[:,i,j] = flux(u_cons,1,eq_relax) - g_N[:,i,j] = flux(u_cons,2,eq_relax) + # 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 conservative f,g 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) + f_cons = flux(u_cons,1,eq_relax) + set_node_vars!(f_M,f_cons,eq_relax,dg,i,j) + g_cons = flux(u_cons,2,eq_relax) + set_node_vars!(g_M,g_cons,eq_relax,dg,i,j) end - + # 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) + #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) + # set_node_vars!(f_N,f_cons,eq_relax,dg,i,j) + # g_cons = flux(u_cons,2,eq_relax) + # set_node_vars!(g_N,g_cons,eq_relax,dg,i,j) + #end + 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 = f_N[:,i,j] #flux(u_cons,1,eq_relax) - #wu = g_N[:,i,j] #flux(u_cons,2,eq_relax) - vu = flux(u_cons,1,eq_relax) - wu = flux(u_cons,2,eq_relax) + vu = get_node_vars(f_N,eq_relax,dg,i,j) + wu = get_node_vars(g_N,eq_relax,dg,i,j) # compute relaxation terms du1 = u_node[1] du2 = u_node[2] From accbecbc0c70698121592fdbc4b614bf07d6a667 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 28 May 2024 12:36:12 +0100 Subject: [PATCH 06/20] first version of entropy projection --- examples/tree_2d_dgsem/elixir_jin_xin_euler.jl | 2 +- .../positivity_zhang_shu_dg2d.jl | 17 ++++++++++++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index a00a13f79b1..d2bcc88e3f8 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -5,7 +5,7 @@ using Trixi ############################################################################### # semidiscretization of the compressible Euler equations -epsilon_relaxation = 5.0e-7 +epsilon_relaxation = 5.0e-6 a1 = a2 = a3 = a4 = 30.0 b1 = b2 = b3 = b4 = 30.0 diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index 59f2f4b0f8f..e56b900131a 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -61,9 +61,11 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr nVars = nvariables(eq_relax) RealT = real(dg) u_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) f_N = zeros(RealT, nVars, nnodes_, nnodes_) g_N = zeros(RealT, nVars, nnodes_, nnodes_) u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) @@ -79,9 +81,22 @@ 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 conservative f,g on the M grid + + # 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) + set_node_vars!(w_M,w_ij,eq_relax,dg,i,j) + end + # compute projection of w with M values down to N + multiply_dimensionwise!(w_N,project_M_to_N,w_M) + # use w now to compute new u on grid M + multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) + + # compute nodal values of conservative f,g on the M grid + for j in 1:nnodes_projection, i in 1:nnodes_projection + w_ij = get_node_vars(w_M, eq_relax, dg, i, j) + u_cons = entropy2cons(w_ij, eq_relax) f_cons = flux(u_cons,1,eq_relax) set_node_vars!(f_M,f_cons,eq_relax,dg,i,j) g_cons = flux(u_cons,2,eq_relax) From 669dcaffc18c90b56e7be7ae9bce530a109df1cc Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 28 May 2024 15:19:57 +0100 Subject: [PATCH 07/20] Add VolumeIntegralWeakFormProjection --- .../tree_2d_dgsem/elixir_jin_xin_euler.jl | 3 +- src/Trixi.jl | 1 + src/solvers/dg.jl | 10 ++++ src/solvers/dgsem_tree/dg_2d.jl | 57 +++++++++++++++++++ 4 files changed, 70 insertions(+), 1 deletion(-) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index d2bcc88e3f8..6ee8d11df98 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -31,7 +31,8 @@ initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmho #initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave) polydeg = 3 basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) -solver = DGSEM(basis, Trixi.flux_upwind) +# solver = DGSEM(basis, Trixi.flux_upwind) +solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakFormProjection()) #solver = DGSEM(polydeg = 3, surface_flux = Trixi.flux_upwind) #surface_flux = Trixi.flux_upwind diff --git a/src/Trixi.jl b/src/Trixi.jl index fa5693b52d7..d0ee4b6f2a0 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -233,6 +233,7 @@ export DG, DGSEM, LobattoLegendreBasis, FDSBP, VolumeIntegralWeakForm, VolumeIntegralStrongForm, + VolumeIntegralWeakFormProjection, VolumeIntegralFluxDifferencing, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 0ab947e697a..768ae5f6973 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -50,6 +50,15 @@ struct VolumeIntegralWeakForm <: AbstractVolumeIntegral end create_cache(mesh, equations, ::VolumeIntegralWeakForm, dg, uEltype) = NamedTuple() +""" + VolumeIntegralWeakFormProjection() + +A weak form volume integral type for DG methods with projection of entropy variables. +""" +struct VolumeIntegralWeakFormProjection <: AbstractVolumeIntegral end + +create_cache(mesh, equations, ::VolumeIntegralWeakFormProjection, dg, uEltype) = NamedTuple() + """ VolumeIntegralFluxDifferencing(volume_flux) @@ -416,6 +425,7 @@ function Base.show(io::IO, mime::MIME"text/plain", dg::DG) show(increment_indent(io), mime, dg.surface_integral) summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof) if !(dg.volume_integral isa VolumeIntegralWeakForm) && + !(dg.volume_integral isa VolumeIntegralWeakFormProjection) && !(dg.volume_integral isa VolumeIntegralStrongForm) show(increment_indent(io), mime, dg.volume_integral) end diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 6b66c2d4bfa..f3a8b35af70 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -229,6 +229,63 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 return nothing end +function calc_volume_integral!(du, u, + mesh::Union{TreeMesh{2}}, + nonconservative_terms, equations, + volume_integral::VolumeIntegralWeakFormProjection, + dg::DGSEM, cache) + # prepare local storage for projection + @unpack interpolate_N_to_M, project_M_to_N = dg.basis + nnodes_,nnodes_projection = size(project_M_to_N) + nVars = size(u, 1) + RealT = real(dg) + u_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + f_N = zeros(RealT, nVars, nnodes_, nnodes_) + g_N = zeros(RealT, nVars, nnodes_, nnodes_) + u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + cache_projection = (; u_N, w_N, f_N, g_N, u_M, w_M, f_M, g_M) + + @threaded for element in eachelement(dg, cache) + weak_form_kernel_projection!(du, u, element, mesh, + nonconservative_terms, equations, + dg, cache, cache_projection) + end + + return nothing +end + +@inline function weak_form_kernel_projection!(du, u, + element, mesh::TreeMesh{2}, + nonconservative_terms::False, equations, + dg::DGSEM, cache, cache_projection) + # true * [some floating point value] == [exactly the same floating point value] + # This can (hopefully) be optimized away due to constant propagation. + @unpack derivative_dhat = dg.basis + + # Calculate volume terms in one element + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + flux1 = flux(u_node, 1, equations) + for ii in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[ii, i], flux1, + equations, dg, ii, j, element) + end + + flux2 = flux(u_node, 2, equations) + for jj in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[jj, j], flux2, + equations, dg, i, jj, element) + end + end + + return nothing +end + # flux differencing volume integral. For curved meshes averaging of the # mapping terms, stored in `cache.elements.contravariant_vectors`, is peeled apart # from the evaluation of the physical fluxes in each Cartesian direction From 8f3c9550bcc1af259eb8a4c4233f1fbe07f094aa Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 28 May 2024 15:42:45 +0100 Subject: [PATCH 08/20] current status of entropy projection --- src/auxiliary/precompile.jl | 6 +- .../positivity_zhang_shu_dg2d.jl | 23 ++++-- src/solvers/dgsem/basis_lobatto_legendre.jl | 16 ++-- src/solvers/dgsem_tree/dg_2d.jl | 73 +++++++++++++++---- 4 files changed, 86 insertions(+), 32 deletions(-) diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 70ebf923b87..860dc95a9d9 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -340,9 +340,9 @@ function _precompile_manual_() for RealT in (Float64,), polydeg in 1:7 nnodes_ = polydeg + 1 basis_type = basis_type_dgsem(RealT, nnodes_) - @assert Base.precompile(Tuple{typeof(Trixi.MortarL2), basis_type}) - @assert Base.precompile(Tuple{Type{Trixi.SolutionAnalyzer}, basis_type}) - @assert Base.precompile(Tuple{Type{Trixi.AdaptorL2}, basis_type}) +# @assert Base.precompile(Tuple{typeof(Trixi.MortarL2), basis_type}) +# @assert Base.precompile(Tuple{Type{Trixi.SolutionAnalyzer}, basis_type}) +# @assert Base.precompile(Tuple{Type{Trixi.AdaptorL2}, basis_type}) end # Constructors: callbacks diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index e56b900131a..a9aebea8353 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -56,7 +56,7 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr eq_relax = equations.equations_relaxation # prepare local storage for projection - @unpack interpolate_N_to_M, project_M_to_N = dg.basis + @unpack interpolate_N_to_M, project_M_to_N, filter_modal_to_N = dg.basis nnodes_,nnodes_projection = size(project_M_to_N) nVars = nvariables(eq_relax) RealT = real(dg) @@ -65,10 +65,15 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr f_N = zeros(RealT, nVars, nnodes_, nnodes_) g_N = zeros(RealT, nVars, nnodes_, nnodes_) u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M_raw = zeros(RealT, nVars, nnodes_projection, nnodes_projection) w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + tmp_MxM = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) + tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) + #@threaded for element in eachelement(dg, cache) for element in eachelement(dg, cache) @@ -80,18 +85,20 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr end end # bring elemtn u_N to grid (M+1)x(M+1) - multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N) + multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) # 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) - set_node_vars!(w_M,w_ij,eq_relax,dg,i,j) + set_node_vars!(w_M_raw,w_ij,eq_relax,dg,i,j) end # compute projection of w with M values down to N - multiply_dimensionwise!(w_N,project_M_to_N,w_M) - # use w now to compute new u on grid M - multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) + multiply_dimensionwise!(w_M,filter_modal_to_N,w_M_raw,tmp_MxM) + + #multiply_dimensionwise!(w_N,project_M_to_N,w_M) + #multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) + # compute nodal values of conservative f,g on the M grid for j in 1:nnodes_projection, i in 1:nnodes_projection @@ -103,8 +110,8 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr set_node_vars!(g_M,g_cons,eq_relax,dg,i,j) end # 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) + multiply_dimensionwise!(f_N,project_M_to_N,f_M,tmp_NxM) + multiply_dimensionwise!(g_N,project_M_to_N,g_M,tmp_NxM) #@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) diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index d1c9b909c6e..aadc60f18eb 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -20,7 +20,8 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, BoundaryMatrix <: AbstractMatrix{RealT}, DerivativeMatrix <: AbstractMatrix{RealT}, Matrix_MxN <: AbstractMatrix{RealT}, - Matrix_NxM <: AbstractMatrix{RealT}} <: + Matrix_NxM <: AbstractMatrix{RealT}, + Matrix_MxM <: AbstractMatrix{RealT}} <: AbstractBasisSBP{RealT} nodes::VectorT weights::VectorT @@ -37,6 +38,7 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, # L2 projection operators interpolate_N_to_M::Matrix_MxN # interpolates from N nodes to M nodes project_M_to_N::Matrix_NxM # compute projection via Legendre modes, cut off modes at N, back to N nodes + filter_modal_to_N::Matrix_MxM # compute modal filter via Legendre modes, cut off modes at N, leave it at M nodes end function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) @@ -78,15 +80,17 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ interpolate_N_to_M_ = polynomial_interpolation_matrix(nodes_, nodes_projection) interpolate_N_to_M = Matrix{RealT}(interpolate_N_to_M_) - project_M_to_N_ = calc_projection_matrix(nodes_projection, nodes_) + project_M_to_N_,filter_modal_to_N_ = calc_projection_matrix(nodes_projection, nodes_) project_M_to_N = Matrix{RealT}(project_M_to_N_) + filter_modal_to_N = Matrix{RealT}(filter_modal_to_N_) return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes), typeof(inverse_vandermonde_legendre), typeof(boundary_interpolation), typeof(derivative_matrix), typeof(interpolate_N_to_M), - typeof(project_M_to_N)}(nodes, weights, + typeof(project_M_to_N), + typeof(filter_modal_to_N)}(nodes, weights, inverse_weights, inverse_vandermonde_legendre, boundary_interpolation, @@ -95,7 +99,8 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ derivative_split_transpose, derivative_dhat, interpolate_N_to_M, - project_M_to_N) + project_M_to_N, + filter_modal_to_N) end LobattoLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) = LobattoLegendreBasis(Float64, polydeg; polydeg_projection) @@ -811,8 +816,9 @@ function calc_projection_matrix(nodes_in,nodes_out) filter_matrix[j,j] = 1 end interpolate_M_to_N = polynomial_interpolation_matrix(nodes_in, nodes_out) + filter_modal = vandermonde_in * filter_matrix * inverse_vandermonde_in projection_matrix = interpolate_M_to_N * vandermonde_in * filter_matrix * inverse_vandermonde_in - return projection_matrix + return projection_matrix, filter_modal end end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index f3a8b35af70..1ea05fbe97a 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -237,31 +237,72 @@ function calc_volume_integral!(du, u, # prepare local storage for projection @unpack interpolate_N_to_M, project_M_to_N = dg.basis nnodes_,nnodes_projection = size(project_M_to_N) - nVars = size(u, 1) + nVars = nvariables(equations) RealT = real(dg) - u_N = zeros(RealT, nVars, nnodes_, nnodes_) - w_N = zeros(RealT, nVars, nnodes_, nnodes_) - f_N = zeros(RealT, nVars, nnodes_, nnodes_) - g_N = zeros(RealT, nVars, nnodes_, nnodes_) - u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - cache_projection = (; u_N, w_N, f_N, g_N, u_M, w_M, f_M, g_M) + u_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + f_N = zeros(RealT, nVars, nnodes_, nnodes_) + g_N = zeros(RealT, nVars, nnodes_, nnodes_) + u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M_raw = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + + tmp_MxM = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) + tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) @threaded for element in eachelement(dg, cache) - weak_form_kernel_projection!(du, u, element, mesh, + # get element u_N + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + for v in eachvariable(equations) + u_N[v,i,j] = u_node[v] + end + end + # bring elemtn u_N to grid (M+1)x(M+1) + multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) + + # 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, equations, dg, i, j) + w_ij = cons2entropy(u_cons,equations) + set_node_vars!(w_M_raw,w_ij,equations,dg,i,j) + end + # compute projection of w with M values down to N + multiply_dimensionwise!(w_M,filter_modal_to_N,w_M_raw,tmp_MxM) + + #multiply_dimensionwise!(w_N,project_M_to_N,w_M) + #multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) + + + # compute nodal values of conservative f,g on the M grid + for j in 1:nnodes_projection, i in 1:nnodes_projection + w_ij = get_node_vars(w_M, equations, dg, i, j) + u_cons = entropy2cons(w_ij, equations) + f_cons = flux(u_cons,1,equations) + set_node_vars!(f_M,f_cons,equations,dg,i,j) + g_cons = flux(u_cons,2,equations) + set_node_vars!(g_M,g_cons,equations,dg,i,j) + end + # compute projection of f with M values down to N, same for g + multiply_dimensionwise!(f_N,project_M_to_N,f_M,tmp_NxM) + multiply_dimensionwise!(g_N,project_M_to_N,g_M,tmp_NxM) + + + weak_form_kernel_projection!(du, u,f_N, g_N, element, mesh, nonconservative_terms, equations, - dg, cache, cache_projection) + dg, cache) end return nothing end -@inline function weak_form_kernel_projection!(du, u, +@inline function weak_form_kernel_projection!(du, u, f_N, g_N, element, mesh::TreeMesh{2}, nonconservative_terms::False, equations, - dg::DGSEM, cache, cache_projection) + dg::DGSEM, cache) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @unpack derivative_dhat = dg.basis @@ -270,13 +311,13 @@ end for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, element) - flux1 = flux(u_node, 1, equations) + flux1 = get_node_vars(f_N, equations, dg, i,j) for ii in eachnode(dg) multiply_add_to_node_vars!(du, derivative_dhat[ii, i], flux1, equations, dg, ii, j, element) end - flux2 = flux(u_node, 2, equations) + flux2 = get_node_vars(g_N, equations, dg, i,j) for jj in eachnode(dg) multiply_add_to_node_vars!(du, derivative_dhat[jj, j], flux2, equations, dg, i, jj, element) From 811c9637f2471258cdfa4965670ed877ad00ea20 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 28 May 2024 16:10:58 +0100 Subject: [PATCH 09/20] first draft of projection dgsem --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 12 ++- ...kelvin_helmholtz_instability_projection.jl | 88 +++++++++++++++++++ src/solvers/dgsem_tree/dg_2d.jl | 2 +- 3 files changed, 98 insertions(+), 4 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index e634a383cdf..c32aae817a7 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -8,9 +8,15 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition = initial_condition_weak_blast_wave -volume_flux = flux_ranocha -solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) +#volume_flux = flux_ranocha +#solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, +# volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) +# +surface_flux = flux_ranocha +polydeg = 1 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 4 * polydeg) +volume_integral = VolumeIntegralWeakFormProjection() +solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-2.0, -2.0) coordinates_max = (2.0, 2.0) diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl new file mode 100644 index 00000000000..9232c12c6b0 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl @@ -0,0 +1,88 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D) + +A version of the classical Kelvin-Helmholtz instability based on +- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021) + A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations + of the Euler Equations + [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017) +""" +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_kelvin_helmholtz_instability + +surface_flux = flux_lax_friedrichs +polydeg = 3 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) +#indicator_sc = IndicatorHennemannGassner(equations, basis, +# alpha_max = 0.000, +# alpha_min = 0.0000, +# alpha_smooth = true, +# variable = density_pressure) +# +#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; +# volume_flux_dg = volume_flux, +# volume_flux_fv = surface_flux) +volume_integral = VolumeIntegralWeakFormProjection() +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 = 100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# 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.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + #save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + 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 diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 1ea05fbe97a..7b8208039b1 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -235,7 +235,7 @@ function calc_volume_integral!(du, u, volume_integral::VolumeIntegralWeakFormProjection, dg::DGSEM, cache) # prepare local storage for projection - @unpack interpolate_N_to_M, project_M_to_N = dg.basis + @unpack interpolate_N_to_M, project_M_to_N, filter_modal_to_N = dg.basis nnodes_,nnodes_projection = size(project_M_to_N) nVars = nvariables(equations) RealT = real(dg) From fdc414ee410d6a6d92685c592a8cc662f9a35c31 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 28 May 2024 18:01:07 +0100 Subject: [PATCH 10/20] current state --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 4 +- src/callbacks_step/analysis_dg2d.jl | 9 ++- src/solvers/dgsem_tree/dg_2d.jl | 73 +++++++++++++++++------ 3 files changed, 66 insertions(+), 20 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index c32aae817a7..2b273091b6a 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -46,11 +46,11 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 1.0) +stepsize_callback = StepsizeCallback(cfl = 0.1) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_solution, + #save_solution, stepsize_callback) ############################################################################### diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index de6b9a2a4a6..5249ba74ea4 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -239,13 +239,20 @@ function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) + u_original = similar(u) + u_original .= u + calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ - integrate_via_indices(u, mesh, equations, dg, cache, + result = integrate_via_indices(u, mesh, equations, dg, cache, du) do u, i, j, element, equations, dg, du u_node = get_node_vars(u, equations, dg, i, j, element) du_node = get_node_vars(du, equations, dg, i, j, element) dot(cons2entropy(u_node, equations), du_node) end + + u .= u_original + return result end function analyze(::Val{:l2_divb}, du, u, t, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 7b8208039b1..5776a2f4d5f 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -113,6 +113,12 @@ function rhs!(du, u, t, mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} + # Copy u for safekeeping + u_original = similar(u) + u_original .= u + + calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # Reset du @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) @@ -175,9 +181,55 @@ function rhs!(du, u, t, calc_sources!(du, u, t, source_terms, equations, dg, cache) end + u .= u_original + return nothing end +function calc_entropy_projection!(u_projected, u, mesh, equations, dg, cache) + # prepare local storage for projection + @unpack interpolate_N_to_M, project_M_to_N = dg.basis + nnodes_,nnodes_projection = size(project_M_to_N) + nVars = nvariables(equations) + RealT = real(dg) + u_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + + tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) + tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) + + @threaded 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) + for v in eachvariable(equations) + u_N[v,i,j] = u_node[v] + end + end + # bring elemtn u_N to grid (M+1)x(M+1) + multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) + + # 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, equations, dg, i, j) + w_ij = cons2entropy(u_cons,equations) + set_node_vars!(w_M,w_ij,equations,dg,i,j) + end + + # compute projection of w with M values down to N + multiply_dimensionwise!(w_N, project_M_to_N ,w_M, tmp_NxM) + + # compute nodal values of conservative variables from the projected entropy variables + for j in eachnode(dg), i in eachnode(dg) + w_ij = get_node_vars(w_N, equations, dg, i, j) + u_cons = entropy2cons(w_ij, equations) + set_node_vars!(u_projected, u_cons, equations, dg, i, j, element) + end + end +end + function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, @@ -257,26 +309,14 @@ function calc_volume_integral!(du, u, # get element u_N for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, element) + w_ij = cons2entropy(u_node, equations) for v in eachvariable(equations) - u_N[v,i,j] = u_node[v] + w_N[v,i,j] = w_ij[v] end end # bring elemtn u_N to grid (M+1)x(M+1) - multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) - - # 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, equations, dg, i, j) - w_ij = cons2entropy(u_cons,equations) - set_node_vars!(w_M_raw,w_ij,equations,dg,i,j) - end - # compute projection of w with M values down to N - multiply_dimensionwise!(w_M,filter_modal_to_N,w_M_raw,tmp_MxM) - - #multiply_dimensionwise!(w_N,project_M_to_N,w_M) - #multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) - - + multiply_dimensionwise!(w_M, interpolate_N_to_M, w_N, tmp_MxN) + # compute nodal values of conservative f,g on the M grid for j in 1:nnodes_projection, i in 1:nnodes_projection w_ij = get_node_vars(w_M, equations, dg, i, j) @@ -290,7 +330,6 @@ function calc_volume_integral!(du, u, multiply_dimensionwise!(f_N,project_M_to_N,f_M,tmp_NxM) multiply_dimensionwise!(g_N,project_M_to_N,g_M,tmp_NxM) - weak_form_kernel_projection!(du, u,f_N, g_N, element, mesh, nonconservative_terms, equations, dg, cache) From 559dfa296a930f1719f0c7a3b8b7d17398213b0b Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 28 May 2024 19:30:34 +0100 Subject: [PATCH 11/20] Add modal cutoff filter --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 7 +-- ...kelvin_helmholtz_instability_projection.jl | 11 +++-- src/callbacks_step/analysis_dg2d.jl | 8 +-- src/solvers/dgsem/basis_lobatto_legendre.jl | 13 +++-- src/solvers/dgsem_tree/dg_2d.jl | 49 ++++++++++++++++++- 5 files changed, 72 insertions(+), 16 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index 2b273091b6a..5e98f1ce785 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -13,9 +13,10 @@ initial_condition = initial_condition_weak_blast_wave # volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # surface_flux = flux_ranocha -polydeg = 1 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 4 * polydeg) -volume_integral = VolumeIntegralWeakFormProjection() +polydeg = 6 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 1) +# volume_integral = VolumeIntegralWeakFormProjection() +volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-2.0, -2.0) diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl index 9232c12c6b0..77c9572f717 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl @@ -33,8 +33,9 @@ end initial_condition = initial_condition_kelvin_helmholtz_instability surface_flux = flux_lax_friedrichs -polydeg = 3 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) +# polydeg = 3 +polydeg = 10 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) #indicator_sc = IndicatorHennemannGassner(equations, basis, # alpha_max = 0.000, # alpha_min = 0.0000, @@ -44,7 +45,8 @@ basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) #volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; # volume_flux_dg = volume_flux, # volume_flux_fv = surface_flux) -volume_integral = VolumeIntegralWeakFormProjection() +# volume_integral = VolumeIntegralWeakFormProjection() +volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-1.0, -1.0) @@ -72,7 +74,8 @@ save_solution = SaveSolutionCallback(interval=1000, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.5) +# stepsize_callback = StepsizeCallback(cfl = 0.5) +stepsize_callback = StepsizeCallback(cfl = 1.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 5249ba74ea4..53f586938b1 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -239,9 +239,9 @@ function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) - u_original = similar(u) - u_original .= u - calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # u_original = similar(u) + # u_original .= u + # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ result = integrate_via_indices(u, mesh, equations, dg, cache, @@ -251,7 +251,7 @@ function analyze(::typeof(entropy_timederivative), du, u, t, dot(cons2entropy(u_node, equations), du_node) end - u .= u_original + # u .= u_original return result end diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index aadc60f18eb..1bd3d9efa18 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -39,9 +39,10 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, interpolate_N_to_M::Matrix_MxN # interpolates from N nodes to M nodes project_M_to_N::Matrix_NxM # compute projection via Legendre modes, cut off modes at N, back to N nodes filter_modal_to_N::Matrix_MxM # compute modal filter via Legendre modes, cut off modes at N, leave it at M nodes + filter_modal_to_cutoff::DerivativeMatrix # compute modal cutoff filter via Legendre modes, cut off modes at polydeg_cutoff end -function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) +function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg, polydeg_cutoff::Integer = div(polydeg + 1, 2) - 1) nnodes_ = polydeg + 1 # compute everything using `Float64` by default @@ -84,6 +85,11 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ project_M_to_N = Matrix{RealT}(project_M_to_N_) filter_modal_to_N = Matrix{RealT}(filter_modal_to_N_) + nnodes_cutoff = polydeg_cutoff + 1 + nodes_cutoff, weights_cutoff = gauss_lobatto_nodes_weights(nnodes_cutoff) + _, filter_modal_to_cutoff_ = calc_projection_matrix(nodes_, nodes_cutoff) + filter_modal_to_cutoff = Matrix{RealT}(filter_modal_to_cutoff_) + return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes), typeof(inverse_vandermonde_legendre), typeof(boundary_interpolation), @@ -100,10 +106,11 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ derivative_dhat, interpolate_N_to_M, project_M_to_N, - filter_modal_to_N) + filter_modal_to_N, + filter_modal_to_cutoff) end -LobattoLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) = LobattoLegendreBasis(Float64, polydeg; polydeg_projection) +LobattoLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg, polydeg_cutoff::Integer = div(polydeg + 1, 2) - 1) = LobattoLegendreBasis(Float64, polydeg; polydeg_projection, polydeg_cutoff) function Base.show(io::IO, basis::LobattoLegendreBasis) @nospecialize basis # reduce precompilation time diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 5776a2f4d5f..899fda82747 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -116,8 +116,20 @@ function rhs!(du, u, t, # Copy u for safekeeping u_original = similar(u) u_original .= u + u_filter_cons = similar(u) + u_filter_prim = similar(u) - calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # calc_filter!(u, u_original, cons2entropy, entropy2cons, mesh, equations, dg, cache) + # calc_filter!(u, u_original, cons2cons, cons2cons, mesh, equations, dg, cache) + # calc_filter!(u, u_original, cons2prim, prim2cons, mesh, equations, dg, cache) + calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) + calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) + + u .= u_filter_cons + @. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] + @. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) + # @autoinfiltrate # Reset du @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) @@ -181,7 +193,7 @@ function rhs!(du, u, t, calc_sources!(du, u, t, source_terms, equations, dg, cache) end - u .= u_original + # u .= u_original return nothing end @@ -230,6 +242,39 @@ function calc_entropy_projection!(u_projected, u, mesh, equations, dg, cache) end end +function calc_filter!(u_filtered, u, cons2filter, filter2cons, mesh, equations, dg, cache) + # prepare local storage for projection + @unpack filter_modal_to_cutoff = dg.basis + nnodes_ = nnodes(dg) + nVars = nvariables(equations) + RealT = real(dg) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N_filtered = zeros(RealT, nVars, nnodes_, nnodes_) + + tmp_NxN = zeros(RealT, nVars, nnodes_, nnodes_) + + @threaded for element in eachelement(dg, cache) + # convert u to entropy variables + for j in eachnode(dg), i in eachnode(dg) + u_cons = get_node_vars(u, equations, dg, i, j, element) + w_ij = cons2filter(u_cons, equations) + for v in eachvariable(equations) + w_N[v,i,j] = w_ij[v] + end + end + + # filter entropy variables + multiply_dimensionwise!(w_N_filtered, filter_modal_to_cutoff, w_N, tmp_NxN) + + # compute nodal values of conservative variables from the projected entropy variables + for j in eachnode(dg), i in eachnode(dg) + w_ij = get_node_vars(w_N_filtered, equations, dg, i, j) + u_cons = filter2cons(w_ij, equations) + set_node_vars!(u_filtered, u_cons, equations, dg, i, j, element) + end + end +end + function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, From 7df151648d1819f9511cf34b44ff94c428258f46 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 28 May 2024 23:48:41 +0100 Subject: [PATCH 12/20] Beginning of using JinXin + filtering --- examples/tree_2d_dgsem/elixir_jin_xin_euler.jl | 8 +++++--- src/solvers/dgsem_tree/dg_2d.jl | 8 ++++---- 2 files changed, 9 insertions(+), 7 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..e06318a500e 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -29,10 +29,12 @@ end #initial_condition = initial_condition_constant initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability) #initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave) -polydeg = 3 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) +polydeg = 5 +polydeg_cutoff = 3 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 1 * polydeg, polydeg_cutoff = polydeg_cutoff) # solver = DGSEM(basis, Trixi.flux_upwind) -solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakFormProjection()) +#solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakFormProjection()) +solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakForm()) #solver = DGSEM(polydeg = 3, surface_flux = Trixi.flux_upwind) #surface_flux = Trixi.flux_upwind diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 899fda82747..1bb0ce91df0 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -117,18 +117,18 @@ function rhs!(du, u, t, u_original = similar(u) u_original .= u u_filter_cons = similar(u) - u_filter_prim = similar(u) + # u_filter_prim = similar(u) # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) # calc_filter!(u, u_original, cons2entropy, entropy2cons, mesh, equations, dg, cache) # calc_filter!(u, u_original, cons2cons, cons2cons, mesh, equations, dg, cache) # calc_filter!(u, u_original, cons2prim, prim2cons, mesh, equations, dg, cache) calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) - calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) + #calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) u .= u_filter_cons - @. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] - @. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) + #@. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] + #@. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) # @autoinfiltrate # Reset du From aabec9c14fac346165be71311d9b7a12808e94f3 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Wed, 29 May 2024 12:37:19 +0100 Subject: [PATCH 13/20] filter testing stuff --- .../elixir_euler_density_wave.jl | 39 ++++++++++++++++--- examples/tree_2d_dgsem/elixir_euler_ec.jl | 18 +++++++-- src/solvers/dgsem_tree/dg_2d.jl | 18 ++++----- 3 files changed, 57 insertions(+), 18 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_density_wave.jl b/examples/tree_2d_dgsem/elixir_euler_density_wave.jl index 0f9e232fa14..e7736252dcf 100644 --- a/examples/tree_2d_dgsem/elixir_euler_density_wave.jl +++ b/examples/tree_2d_dgsem/elixir_euler_density_wave.jl @@ -1,15 +1,32 @@ using OrdinaryDiffEq using Trixi +using LinearAlgebra ############################################################################### # semidiscretization of the compressible Euler equations gamma = 1.4 equations = CompressibleEulerEquations2D(gamma) -initial_condition = initial_condition_density_wave - -solver = DGSEM(polydeg = 5, surface_flux = flux_central) +function initial_condition_density_wave2(x, t, equations::CompressibleEulerEquations2D) + v1 = 0.1 + v2 = 0.2 + rho = 1 + 0.98 * sinpi(2 * (x[1] + x[2] - t * (v1 + v2))) + rho_v1 = rho * v1 + rho_v2 = rho * v2 + p = 20 + rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2) + return SVector(rho, rho_v1, rho_v2, rho_e) +end +initial_condition = initial_condition_density_wave2 + +surface_flux = flux_central +polydeg = 17 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 5) +volume_integral = VolumeIntegralWeakForm() +solver = DGSEM(basis, surface_flux, volume_integral) + +# solver = DGSEM(polydeg = 5, surface_flux = flux_central) coordinates_min = (-1.0, -1.0) coordinates_max = (1.0, 1.0) @@ -17,8 +34,20 @@ mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 2, n_cells_max = 30_000) +@info "Create semi..." semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +@info "Compute Jacobian..." +J = jacobian_ad_forward(semi) + +@info "Compute eigenvalues..." +λ = eigvals(J) + +@info "max real part" maximum(real.(λ)) +# @info "Plot spectrum..." +# scatter(real.(λ), imag.(λ), label="central flux") +wololo + ############################################################################### # ODE solvers, callbacks etc. @@ -27,7 +56,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 100 * 20 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) @@ -41,7 +70,7 @@ stepsize_callback = StepsizeCallback(cfl = 1.6) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_solution, + # save_solution, stepsize_callback) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index 5e98f1ce785..ced8a5221a2 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -1,20 +1,30 @@ using OrdinaryDiffEq using Trixi - +using Random: seed! ############################################################################### # semidiscretization of the compressible Euler equations equations = CompressibleEulerEquations2D(1.4) +seed!(1) +function initial_condition_random_field(x, t, equations::CompressibleEulerEquations2D) +amplitude = 1.5 +rho = 2 + amplitude * rand() +v1 = -3.1 + amplitude * rand() +v2 = 1.3 + amplitude * rand() +p = 7.54 + amplitude * rand() +return prim2cons(SVector(rho, v1, v2, p), equations) +end initial_condition = initial_condition_weak_blast_wave +# initial_condition = initial_condition_random_field #volume_flux = flux_ranocha #solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, # volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # surface_flux = flux_ranocha -polydeg = 6 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 1) +polydeg = 11 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) # volume_integral = VolumeIntegralWeakFormProjection() volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) @@ -47,7 +57,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.1) +stepsize_callback = StepsizeCallback(cfl = 0.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 1bb0ce91df0..50b8e7eac35 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -114,19 +114,19 @@ function rhs!(du, u, t, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Copy u for safekeeping - u_original = similar(u) - u_original .= u + #u_original = similar(u) + #u_original .= u u_filter_cons = similar(u) - # u_filter_prim = similar(u) + u_filter_prim = similar(u) # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) - # calc_filter!(u, u_original, cons2entropy, entropy2cons, mesh, equations, dg, cache) - # calc_filter!(u, u_original, cons2cons, cons2cons, mesh, equations, dg, cache) - # calc_filter!(u, u_original, cons2prim, prim2cons, mesh, equations, dg, cache) - calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) + #calc_filter!(u, u, cons2entropy, entropy2cons, mesh, equations, dg, cache) + # calc_filter!(u, u, cons2cons, cons2cons, mesh, equations, dg, cache) + calc_filter!(u, u, cons2prim, prim2cons, mesh, equations, dg, cache) + #calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) #calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) - u .= u_filter_cons + #u .= u_filter_cons #@. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] #@. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) # @autoinfiltrate @@ -247,7 +247,7 @@ function calc_filter!(u_filtered, u, cons2filter, filter2cons, mesh, equations, @unpack filter_modal_to_cutoff = dg.basis nnodes_ = nnodes(dg) nVars = nvariables(equations) - RealT = real(dg) + RealT = eltype(u) w_N = zeros(RealT, nVars, nnodes_, nnodes_) w_N_filtered = zeros(RealT, nVars, nnodes_, nnodes_) From 306d47eb4f07f064c554620ac5773afcbf3f9079 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 29 May 2024 15:15:28 +0100 Subject: [PATCH 14/20] Implement the Gauss-Legendre basis for DGSEM --- .../elixir_euler_convergence_pure_fv.jl | 14 +- src/Trixi.jl | 2 +- src/solvers/dgsem/basis_gauss_legendre.jl | 240 ++++++++++++++++++ src/solvers/dgsem/dgsem.jl | 11 +- src/solvers/dgsem_tree/dg_2d.jl | 86 ++++++- 5 files changed, 344 insertions(+), 9 deletions(-) create mode 100644 src/solvers/dgsem/basis_gauss_legendre.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl b/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl index 96184b5ba47..ac4b3709465 100644 --- a/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl +++ b/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl @@ -11,14 +11,18 @@ initial_condition = initial_condition_convergence_test surface_flux = flux_hllc -basis = LobattoLegendreBasis(3) -volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc) +# basis = LobattoLegendreBasis(3) +# volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc) +# solver = DGSEM(basis, surface_flux, volume_integral) +polydeg = 4 +basis = GaussLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 4) +volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (2.0, 2.0) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 4, + initial_refinement_level = 2, n_cells_max = 10_000) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, @@ -32,7 +36,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 100 * 10 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) @@ -46,7 +50,7 @@ stepsize_callback = StepsizeCallback(cfl = 0.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_solution, + # save_solution, stepsize_callback) ############################################################################### diff --git a/src/Trixi.jl b/src/Trixi.jl index d0ee4b6f2a0..392138b0f8d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -230,7 +230,7 @@ export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMe T8codeMesh export DG, - DGSEM, LobattoLegendreBasis, + DGSEM, LobattoLegendreBasis, GaussLegendreBasis, FDSBP, VolumeIntegralWeakForm, VolumeIntegralStrongForm, VolumeIntegralWeakFormProjection, diff --git a/src/solvers/dgsem/basis_gauss_legendre.jl b/src/solvers/dgsem/basis_gauss_legendre.jl new file mode 100644 index 00000000000..b9aafeb9da5 --- /dev/null +++ b/src/solvers/dgsem/basis_gauss_legendre.jl @@ -0,0 +1,240 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +""" + GaussLegendreBasis([RealT=Float64,] polydeg::Integer) + +Create a nodal Gauss-Legendre basis for polynomials of degree `polydeg`. +""" +struct GaussLegendreBasis{RealT <: Real, NNODES, + VectorT <: AbstractVector{RealT}, + InverseVandermondeLegendre <: AbstractMatrix{RealT}, + BoundaryMatrix <: AbstractMatrix{RealT}, + DerivativeMatrix <: AbstractMatrix{RealT}, + Matrix_MxN <: AbstractMatrix{RealT}, + Matrix_NxM <: AbstractMatrix{RealT}, + Matrix_MxM <: AbstractMatrix{RealT}} <: + AbstractBasisSBP{RealT} + nodes::VectorT + weights::VectorT + inverse_weights::VectorT + + inverse_vandermonde_legendre::InverseVandermondeLegendre + boundary_interpolation::BoundaryMatrix # lhat + + derivative_matrix::DerivativeMatrix # strong form derivative matrix + derivative_split::DerivativeMatrix # strong form derivative matrix minus boundary terms + derivative_split_transpose::DerivativeMatrix # transpose of `derivative_split` + derivative_dhat::DerivativeMatrix # weak form matrix "dhat", + # negative adjoint wrt the SBP dot product + # L2 projection operators + interpolate_N_to_M::Matrix_MxN # interpolates from N nodes to M nodes + project_M_to_N::Matrix_NxM # compute projection via Legendre modes, cut off modes at N, back to N nodes + filter_modal_to_N::Matrix_MxM # compute modal filter via Legendre modes, cut off modes at N, leave it at M nodes + filter_modal_to_cutoff::DerivativeMatrix # compute modal cutoff filter via Legendre modes, cut off modes at polydeg_cutoff +end + +function GaussLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg, polydeg_cutoff::Integer = div(polydeg + 1, 2) - 1) + nnodes_ = polydeg + 1 + + # compute everything using `Float64` by default + nodes_, weights_ = gauss_nodes_weights(nnodes_) + inverse_weights_ = inv.(weights_) + + _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_) + + boundary_interpolation_ = zeros(nnodes_, 2) + boundary_interpolation_[:, 1] = calc_lhat(-1.0, nodes_, weights_) + boundary_interpolation_[:, 2] = calc_lhat(1.0, nodes_, weights_) + + derivative_matrix_ = polynomial_derivative_matrix(nodes_) + derivative_split_ = calc_dsplit(nodes_, weights_) + derivative_split_transpose_ = Matrix(derivative_split_') + derivative_dhat_ = calc_dhat(nodes_, weights_) + + # type conversions to get the requested real type and enable possible + # optimizations of runtime performance and latency + nodes = SVector{nnodes_, RealT}(nodes_) + weights = SVector{nnodes_, RealT}(weights_) + inverse_weights = SVector{nnodes_, RealT}(inverse_weights_) + + inverse_vandermonde_legendre = convert.(RealT, inverse_vandermonde_legendre_) + boundary_interpolation = convert.(RealT, boundary_interpolation_) + + # Usually as fast as `SMatrix` (when using `let` in the volume integral/`@threaded`) + derivative_matrix = Matrix{RealT}(derivative_matrix_) + derivative_split = Matrix{RealT}(derivative_split_) + derivative_split_transpose = Matrix{RealT}(derivative_split_transpose_) + derivative_dhat = Matrix{RealT}(derivative_dhat_) + + # L2 projection operators + nnodes_projection = polydeg_projection + 1 + nodes_projection, weights_projection = gauss_nodes_weights(nnodes_projection) + interpolate_N_to_M_ = polynomial_interpolation_matrix(nodes_, nodes_projection) + interpolate_N_to_M = Matrix{RealT}(interpolate_N_to_M_) + + project_M_to_N_,filter_modal_to_N_ = calc_projection_matrix(nodes_projection, nodes_) + project_M_to_N = Matrix{RealT}(project_M_to_N_) + filter_modal_to_N = Matrix{RealT}(filter_modal_to_N_) + + nnodes_cutoff = polydeg_cutoff + 1 + nodes_cutoff, weights_cutoff = gauss_nodes_weights(nnodes_cutoff) + _, filter_modal_to_cutoff_ = calc_projection_matrix(nodes_, nodes_cutoff) + filter_modal_to_cutoff = Matrix{RealT}(filter_modal_to_cutoff_) + + return GaussLegendreBasis{RealT, nnodes_, typeof(nodes), + typeof(inverse_vandermonde_legendre), + typeof(boundary_interpolation), + typeof(derivative_matrix), + typeof(interpolate_N_to_M), + typeof(project_M_to_N), + typeof(filter_modal_to_N)}(nodes, weights, + inverse_weights, + inverse_vandermonde_legendre, + boundary_interpolation, + derivative_matrix, + derivative_split, + derivative_split_transpose, + derivative_dhat, + interpolate_N_to_M, + project_M_to_N, + filter_modal_to_N, + filter_modal_to_cutoff) +end + +GaussLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg, polydeg_cutoff::Integer = div(polydeg + 1, 2) - 1) = GaussLegendreBasis(Float64, polydeg; polydeg_projection, polydeg_cutoff) + +function Base.show(io::IO, basis::GaussLegendreBasis) + @nospecialize basis # reduce precompilation time + + print(io, "GaussLegendreBasis{", real(basis), "}(polydeg=", polydeg(basis), ")") +end +function Base.show(io::IO, ::MIME"text/plain", basis::GaussLegendreBasis) + @nospecialize basis # reduce precompilation time + + print(io, "GaussLegendreBasis{", real(basis), "} with polynomials of degree ", + polydeg(basis)) +end + +function Base.:(==)(b1::GaussLegendreBasis, b2::GaussLegendreBasis) + if typeof(b1) != typeof(b2) + return false + end + + for field in fieldnames(typeof(b1)) + if getfield(b1, field) != getfield(b2, field) + return false + end + end + + return true +end + +@inline Base.real(basis::GaussLegendreBasis{RealT}) where {RealT} = RealT + +@inline function nnodes(basis::GaussLegendreBasis{RealT, NNODES}) where {RealT, NNODES + } + NNODES +end + +""" + eachnode(basis::GaussLegendreBasis) + +Return an iterator over the indices that specify the location in relevant data structures +for the nodes in `basis`. +In particular, not the nodes themselves are returned. +""" +@inline eachnode(basis::GaussLegendreBasis) = Base.OneTo(nnodes(basis)) + +@inline polydeg(basis::GaussLegendreBasis) = nnodes(basis) - 1 + +@inline get_nodes(basis::GaussLegendreBasis) = basis.nodes + +""" + integrate(f, u, basis::GaussLegendreBasis) + +Map the function `f` to the coefficients `u` and integrate with respect to the +quadrature rule given by `basis`. +""" +function integrate(f, u, basis::GaussLegendreBasis) + @unpack weights = basis + + res = zero(f(first(u))) + for i in eachindex(u, weights) + res += f(u[i]) * weights[i] + end + return res +end + +# Return the first/last weight of the quadrature associated with `basis`. +# Since the mass matrix for nodal Gauss-Legendre bases is diagonal, +# these weights are the only coefficients necessary for the scaling of +# surface terms/integrals in DGSEM. +left_boundary_weight(basis::GaussLegendreBasis) = first(basis.weights) +right_boundary_weight(basis::GaussLegendreBasis) = last(basis.weights) + +struct GaussLegendreAnalyzer{RealT <: Real, NNODES, + VectorT <: AbstractVector{RealT}, + Vandermonde <: AbstractMatrix{RealT}} <: + SolutionAnalyzer{RealT} + nodes::VectorT + weights::VectorT + vandermonde::Vandermonde +end + +function SolutionAnalyzer(basis::GaussLegendreBasis; + analysis_polydeg = 2 * polydeg(basis)) + RealT = real(basis) + nnodes_ = analysis_polydeg + 1 + + # compute everything using `Float64` by default + nodes_, weights_ = gauss_nodes_weights(nnodes_) + vandermonde_ = polynomial_interpolation_matrix(get_nodes(basis), nodes_) + + # type conversions to get the requested real type and enable possible + # optimizations of runtime performance and latency + nodes = SVector{nnodes_, RealT}(nodes_) + weights = SVector{nnodes_, RealT}(weights_) + + vandermonde = Matrix{RealT}(vandermonde_) + + return GaussLegendreAnalyzer{RealT, nnodes_, typeof(nodes), typeof(vandermonde)}(nodes, + weights, + vandermonde) +end + +function Base.show(io::IO, analyzer::GaussLegendreAnalyzer) + @nospecialize analyzer # reduce precompilation time + + print(io, "GaussLegendreAnalyzer{", real(analyzer), "}(polydeg=", + polydeg(analyzer), ")") +end +function Base.show(io::IO, ::MIME"text/plain", analyzer::GaussLegendreAnalyzer) + @nospecialize analyzer # reduce precompilation time + + print(io, "GaussLegendreAnalyzer{", real(analyzer), + "} with polynomials of degree ", polydeg(analyzer)) +end + +@inline Base.real(analyzer::GaussLegendreAnalyzer{RealT}) where {RealT} = RealT + +@inline function nnodes(analyzer::GaussLegendreAnalyzer{RealT, NNODES}) where {RealT, + NNODES} + NNODES +end +""" + eachnode(analyzer::GaussLegendreAnalyzer) + +Return an iterator over the indices that specify the location in relevant data structures +for the nodes in `analyzer`. +In particular, not the nodes themselves are returned. +""" +@inline eachnode(analyzer::GaussLegendreAnalyzer) = Base.OneTo(nnodes(analyzer)) + +@inline polydeg(analyzer::GaussLegendreAnalyzer) = nnodes(analyzer) - 1 + +end # @muladd diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl index 27caad4d2dc..d3d7b9f9310 100644 --- a/src/solvers/dgsem/dgsem.jl +++ b/src/solvers/dgsem/dgsem.jl @@ -9,6 +9,7 @@ include("interpolation.jl") include("l2projection.jl") include("basis_lobatto_legendre.jl") +include("basis_gauss_legendre.jl") """ DGSEM(; RealT=Float64, polydeg::Integer, @@ -20,7 +21,7 @@ include("basis_lobatto_legendre.jl") Create a discontinuous Galerkin spectral element method (DGSEM) using a [`LobattoLegendreBasis`](@ref) with polynomials of degree `polydeg`. """ -const DGSEM = DG{Basis} where {Basis <: LobattoLegendreBasis} +const DGSEM = DG{Basis} where {Basis <: AbstractBasisSBP} # TODO: Deprecated in v0.3 (no longer documented) function DGSEM(basis::LobattoLegendreBasis, @@ -31,6 +32,14 @@ function DGSEM(basis::LobattoLegendreBasis, return DG{typeof(basis), typeof(mortar), typeof(surface_integral), typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral) end +function DGSEM(basis::GaussLegendreBasis, + surface_flux = flux_central, + volume_integral = VolumeIntegralWeakForm(), + mortar = nothing) + surface_integral = SurfaceIntegralWeakForm(surface_flux) + return DG{typeof(basis), typeof(mortar), typeof(surface_integral), + typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral) +end # TODO: Deprecated in v0.3 (no longer documented) function DGSEM(basis::LobattoLegendreBasis, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 50b8e7eac35..d087e4dcd9d 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -720,7 +720,7 @@ end end function prolong2interfaces!(cache, u, - mesh::TreeMesh{2}, equations, surface_integral, dg::DG) + mesh::TreeMesh{2}, equations, surface_integral, dg::DG{<:LobattoLegendreBasis}) @unpack interfaces = cache @unpack orientations, neighbor_ids = interfaces interfaces_u = interfaces.u @@ -747,6 +747,43 @@ function prolong2interfaces!(cache, u, return nothing end +function prolong2interfaces!(cache, u, + mesh::TreeMesh{2}, equations, surface_integral, dg::DG{<:GaussLegendreBasis}) + @unpack interfaces = cache + @unpack orientations, neighbor_ids = interfaces + @unpack boundary_interpolation, weights = dg.basis + interfaces_u = interfaces.u + + @threaded for interface in eachinterface(dg, cache) + left_element = neighbor_ids[1, interface] + right_element = neighbor_ids[2, interface] + + if orientations[interface] == 1 + # interface in x-direction + for j in eachnode(dg), v in eachvariable(equations) + interfaces_u[1, v, j, interface] = 0 + interfaces_u[2, v, j, interface] = 0 + for ii in eachnode(dg) + interfaces_u[1, v, j, interface] += u[v, ii, j, left_element] * weights[ii] * boundary_interpolation[ii, 2] + interfaces_u[2, v, j, interface] += u[v, ii, j, right_element] * weights[ii] * boundary_interpolation[ii, 1] + end + end + else # if orientations[interface] == 2 + # interface in y-direction + for i in eachnode(dg), v in eachvariable(equations) + interfaces_u[1, v, i, interface] = 0 + interfaces_u[2, v, i, interface] = 0 + for jj in eachnode(dg) + interfaces_u[1, v, i, interface] += u[v, i, jj, left_element] * weights[jj] * boundary_interpolation[jj, 2] + interfaces_u[2, v, i, interface] += u[v, i, jj, right_element] * weights[jj] * boundary_interpolation[jj, 1] + end + end + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, nonconservative_terms::False, equations, @@ -1271,7 +1308,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DG, cache) + dg::DG{<:LobattoLegendreBasis}, cache) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements @@ -1310,6 +1347,51 @@ function calc_surface_integral!(du, u, return nothing end +function calc_surface_integral!(du, u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DG{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation = dg.basis + @unpack surface_flux_values = cache.elements + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # Access the factors only once before beginning the loop to increase performance. + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + for v in eachvariable(equations) + for ii in eachnode(dg) + # surface at -x + du[v, ii, l, element] = (du[v, ii, l, element] - + surface_flux_values[v, l, 1, element] * + boundary_interpolation[ii, 1]) + + # surface at +x + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_values[v, l, 2, element] * + boundary_interpolation[ii, 2]) + end + + for jj in eachnode(dg) + # surface at -y + du[v, l, jj, element] = (du[v, l, jj, element] - + surface_flux_values[v, l, 3, element] * + boundary_interpolation[jj, 1]) + + # surface at +y + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_values[v, l, 4, element] * + boundary_interpolation[jj, 2]) + end + end + end + end + + return nothing +end + function apply_jacobian!(du, mesh::TreeMesh{2}, equations, dg::DG, cache) @unpack inverse_jacobian = cache.elements From fad39b128e65818cf35502cf9b56ce4761603e74 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 29 May 2024 16:05:37 +0100 Subject: [PATCH 15/20] Fix whitespace --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 6 +++--- src/solvers/dgsem_tree/dg_2d.jl | 26 +++++++++++------------ 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index ced8a5221a2..94629fcc596 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -15,15 +15,15 @@ v2 = 1.3 + amplitude * rand() p = 7.54 + amplitude * rand() return prim2cons(SVector(rho, v1, v2, p), equations) end -initial_condition = initial_condition_weak_blast_wave -# initial_condition = initial_condition_random_field +# initial_condition = initial_condition_weak_blast_wave +initial_condition = initial_condition_random_field #volume_flux = flux_ranocha #solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, # volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # surface_flux = flux_ranocha -polydeg = 11 +polydeg = 31 basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) # volume_integral = VolumeIntegralWeakFormProjection() volume_integral = VolumeIntegralWeakForm() diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index d087e4dcd9d..c16434e890c 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -336,19 +336,19 @@ function calc_volume_integral!(du, u, nnodes_,nnodes_projection = size(project_M_to_N) nVars = nvariables(equations) RealT = real(dg) - u_N = zeros(RealT, nVars, nnodes_, nnodes_) - w_N = zeros(RealT, nVars, nnodes_, nnodes_) - f_N = zeros(RealT, nVars, nnodes_, nnodes_) - g_N = zeros(RealT, nVars, nnodes_, nnodes_) - u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - w_M_raw = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - - tmp_MxM = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) - tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) + u_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + f_N = zeros(RealT, nVars, nnodes_, nnodes_) + g_N = zeros(RealT, nVars, nnodes_, nnodes_) + u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M_raw = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + + tmp_MxM = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) + tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) @threaded for element in eachelement(dg, cache) # get element u_N From 1e2aaa3857ac6828ae4495574d7d31589692b3a3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 29 May 2024 18:27:57 +0200 Subject: [PATCH 16/20] draft a generalized version of Jin Xin (#1961) Co-authored-by: Michael Schlottke-Lakemper --- .../tree_2d_dgsem/elixir_jin_xin_euler.jl | 11 +- src/Trixi.jl | 4 +- .../positivity_zhang_shu_dg2d.jl | 6 +- .../subcell_limiter_idp_correction_2d.jl | 24 +- .../jin_xin_compressible_euler_2d.jl | 300 +++++++++--------- 5 files changed, 173 insertions(+), 172 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index e06318a500e..2f993f6a90f 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) @@ -78,11 +79,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 392138b0f8d..acabd30867d 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 a9aebea8353..567217292a3 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, filter_modal_to_N = dg.basis @@ -76,7 +76,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) 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 From 513557406907a9f29af67f297fbd47295d57cb7e Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 29 May 2024 17:34:52 +0100 Subject: [PATCH 17/20] Revert JinXin system behavior to before the generalization --- src/equations/jin_xin_compressible_euler_2d.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/equations/jin_xin_compressible_euler_2d.jl b/src/equations/jin_xin_compressible_euler_2d.jl index edbc202291e..7636ccc7314 100644 --- a/src/equations/jin_xin_compressible_euler_2d.jl +++ b/src/equations/jin_xin_compressible_euler_2d.jl @@ -166,7 +166,8 @@ end 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))..., + # (sqrt_velocities .* (f_rr_base + f_ll_base))..., @ranocha: is this correct? + (sqrt_velocities .* (f_rr_base - f_ll_base))..., zero(u_ll_base)...) else # orientation == 2 sqrt_velocities = equations.sqrt_velocities[2] @@ -174,7 +175,8 @@ end 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))...) + # (sqrt_velocities .* (f_rr_base + f_ll_base))..., # @ranocha: is this correct? + (sqrt_velocities .* (f_rr_base - f_ll_base))...) end return 0.5f0 * (flux(u_ll, orientation, equations) + @@ -184,7 +186,8 @@ end @inline function max_abs_speeds(u, equations::JinXinEquations{2}) return ntuple(Val(ndims(equations))) do n - maximum(equations.sqrt_velocities_inv[n]) + # maximum(equations.sqrt_velocities_inv[n]) @ranocha: is this correct? + maximum(equations.sqrt_velocities[n]) end end From 508e8d685aa40413737cce90b992612ff59e176f Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Thu, 30 May 2024 11:51:12 +0100 Subject: [PATCH 18/20] Add modal filter stage callback --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 13 +++-- src/Trixi.jl | 2 +- src/callbacks_stage/callbacks_stage.jl | 1 + src/callbacks_stage/modal_filter.jl | 62 +++++++++++++++++++++++ src/callbacks_stage/modal_filter_dg2d.jl | 40 +++++++++++++++ src/solvers/dgsem_tree/dg_2d.jl | 3 +- 6 files changed, 116 insertions(+), 5 deletions(-) create mode 100644 src/callbacks_stage/modal_filter.jl create mode 100644 src/callbacks_stage/modal_filter_dg2d.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index 94629fcc596..ef55781b0f7 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -23,8 +23,9 @@ initial_condition = initial_condition_random_field # volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # surface_flux = flux_ranocha -polydeg = 31 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) +polydeg = 19 +# basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) +basis = GaussLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) # volume_integral = VolumeIntegralWeakFormProjection() volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) @@ -67,7 +68,13 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +# Create modal filter and filter initial condition +modal_filter = ModalFilter(solver; polydeg_cutoff = 3, + cons2filter = cons2prim, filter2cons = prim2cons) +modal_filter(ode.u0, semi) + +# sol = solve(ode, CarpenterKennedy2N54(; williamson_condition = false), +sol = solve(ode, CarpenterKennedy2N54(; stage_limiter! = modal_filter, williamson_condition = false), 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 diff --git a/src/Trixi.jl b/src/Trixi.jl index acabd30867d..e154d34c5fb 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -274,7 +274,7 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt, export ControllerThreeLevel, ControllerThreeLevelCombined, IndicatorLöhner, IndicatorLoehner, IndicatorMax -export PositivityPreservingLimiterZhangShu +export PositivityPreservingLimiterZhangShu, ModalFilter export trixi_include, examples_dir, get_examples, default_example, default_example_unstructured, ode_default_options diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index d5abc1d227d..52ac2aa10dc 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -5,6 +5,7 @@ @muladd begin #! format: noindent +include("modal_filter.jl") include("positivity_zhang_shu.jl") include("subcell_limiter_idp_correction.jl") include("subcell_bounds_check.jl") diff --git a/src/callbacks_stage/modal_filter.jl b/src/callbacks_stage/modal_filter.jl new file mode 100644 index 00000000000..86094e8cdb2 --- /dev/null +++ b/src/callbacks_stage/modal_filter.jl @@ -0,0 +1,62 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +""" + ModalFilter(; polydeg_cutoff::Integer, + cons2filter = cons2cons, + filter2cons = cons2cons) + +A modal filter that will cut-off all modes of the solution `u` higher than `polydeg_cutoff`. +The filtering will be done in the filter variables, for which a forward (`cons2filter`) and +reverse (`filter2cons`) transformation is required. By default, the identity transformation +(`cons2cons`) is used. +""" +struct ModalFilter{RealT, Cons2Filter, Filter2Cons} + cons2filter::Cons2Filter + filter2cons::Filter2Cons + filter_matrix::Matrix{RealT} +end + +function ModalFilter(dg; polydeg_cutoff::Integer, cons2filter = cons2cons, + filter2cons = cons2cons) + RealT = real(dg) + if dg.basis isa LobattoLegendreBasis + nodes_cutoff, _ = gauss_lobatto_nodes_weights(polydeg_cutoff) + elseif dg.basis isa GaussLegendreBasis + nodes_cutoff, _ = gauss_nodes_weights(polydeg_cutoff) + else + throw(ArgumentError("Only LobattoLegendreBasis and GaussLegendreBasis are supported.")) + end + _, filter_matrix_ = calc_projection_matrix(dg.basis.nodes, nodes_cutoff) + filter_matrix = Matrix{RealT}(filter_matrix_) + + modal_filter = ModalFilter{RealT, + typeof(cons2filter), + typeof(filter2cons)}(cons2filter, filter2cons, filter_matrix) +end + +# Main function that applies the actual, mesh- and solver-specific filter +function (modal_filter::ModalFilter)(u_ode, semi::AbstractSemidiscretization) + u = wrap_array(u_ode, semi) + @unpack cons2filter, filter2cons, filter_matrix = modal_filter + + @trixi_timeit timer() "modal filter" begin + # println("Applying modal filter") + apply_modal_filter!(u, u, cons2filter, filter2cons, filter_matrix, + mesh_equations_solver_cache(semi)...) + end + + return nothing +end + +# This version is called as the stage limiter version of the filter +function (modal_filter::ModalFilter)(u_ode, integrator, semi::AbstractSemidiscretization, t) + modal_filter(u_ode, semi) +end + +include("modal_filter_dg2d.jl") +end # @muladd diff --git a/src/callbacks_stage/modal_filter_dg2d.jl b/src/callbacks_stage/modal_filter_dg2d.jl new file mode 100644 index 00000000000..1eb05bfec92 --- /dev/null +++ b/src/callbacks_stage/modal_filter_dg2d.jl @@ -0,0 +1,40 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function apply_modal_filter!(u_filtered, u, cons2filter, filter2cons, filter_matrix, + mesh::TreeMesh{2}, equations, dg, cache) + nnodes_ = nnodes(dg) + nvars = nvariables(equations) + RealT = eltype(u) + + u_element = zeros(RealT, nvars, nnodes_, nnodes_) + u_element_filtered = zeros(RealT, nvars, nnodes_, nnodes_) + tmp_NxN = zeros(RealT, nvars, nnodes_, nnodes_) + + @threaded for element in eachelement(dg, cache) + # convert u to filter variables + for j in eachnode(dg), i in eachnode(dg) + u_node_cons = get_node_vars(u, equations, dg, i, j, element) + u_node_filter = cons2filter(u_node_cons, equations) + for v in eachvariable(equations) + u_element[v,i,j] = u_node_filter[v] + end + end + + # Apply modal filter + multiply_dimensionwise!(u_element_filtered, filter_matrix, u_element, tmp_NxN) + + # compute nodal values of conservative variables from the projected entropy variables + for j in eachnode(dg), i in eachnode(dg) + u_node_filter = get_node_vars(u_element_filtered, equations, dg, i, j) + u_node_cons = filter2cons(u_node_filter, equations) + set_node_vars!(u_filtered, u_node_cons, equations, dg, i, j, element) + end + end +end + +end # @muladd \ No newline at end of file diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index c16434e890c..9d2917b9667 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -118,11 +118,12 @@ function rhs!(du, u, t, #u_original .= u u_filter_cons = similar(u) u_filter_prim = similar(u) + # println("rhs! called") # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) #calc_filter!(u, u, cons2entropy, entropy2cons, mesh, equations, dg, cache) # calc_filter!(u, u, cons2cons, cons2cons, mesh, equations, dg, cache) - calc_filter!(u, u, cons2prim, prim2cons, mesh, equations, dg, cache) + # calc_filter!(u, u, cons2prim, prim2cons, mesh, equations, dg, cache) #calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) #calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) From 43605a8e8699491217f8ed8ef80e520cabf7f260 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Thu, 30 May 2024 15:59:44 +0100 Subject: [PATCH 19/20] current status --- ...ixir_euler_kelvin_helmholtz_instability.jl | 11 +++++----- .../elixir_euler_density_wave.jl | 20 +++++++++---------- examples/tree_2d_dgsem/elixir_euler_ec.jl | 8 ++++---- ...ixir_euler_kelvin_helmholtz_instability.jl | 10 +++++----- ...kelvin_helmholtz_instability_projection.jl | 19 +++++++++++------- ...kelvin_helmholtz_instability_sc_subcell.jl | 4 ++-- .../tree_2d_dgsem/elixir_jin_xin_euler.jl | 8 ++++---- src/visualization/utilities.jl | 1 + 8 files changed, 44 insertions(+), 37 deletions(-) diff --git a/examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability.jl b/examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability.jl index 14de0bf0e8b..6d1af59ac08 100644 --- a/examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability.jl +++ b/examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability.jl @@ -1,6 +1,6 @@ using Trixi, OrdinaryDiffEq -dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(), +dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(), surface_integral = SurfaceIntegralWeakForm(FluxLaxFriedrichs()), volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha)) @@ -31,7 +31,7 @@ function initial_condition_kelvin_helmholtz_instability(x, t, end initial_condition = initial_condition_kelvin_helmholtz_instability -cells_per_dimension = (32, 32) +cells_per_dimension = (64, 64) mesh = DGMultiMesh(dg, cells_per_dimension; periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg) @@ -40,8 +40,8 @@ tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -alive_callback = AliveCallback(alive_interval = 10) -analysis_interval = 100 +alive_callback = AliveCallback(alive_interval = 100) +analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg)) callbacks = CallbackSet(summary_callback, analysis_callback, @@ -50,7 +50,8 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +#sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +sol = solve(ode, SSPRK43(), dt = estimate_dt(mesh, dg), save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_euler_density_wave.jl b/examples/tree_2d_dgsem/elixir_euler_density_wave.jl index e7736252dcf..29e37f4a043 100644 --- a/examples/tree_2d_dgsem/elixir_euler_density_wave.jl +++ b/examples/tree_2d_dgsem/elixir_euler_density_wave.jl @@ -34,19 +34,19 @@ mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 2, n_cells_max = 30_000) -@info "Create semi..." -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +# @info "Create semi..." +# semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) -@info "Compute Jacobian..." -J = jacobian_ad_forward(semi) +# @info "Compute Jacobian..." +# J = jacobian_ad_forward(semi) -@info "Compute eigenvalues..." -λ = eigvals(J) +# @info "Compute eigenvalues..." +# λ = eigvals(J) -@info "max real part" maximum(real.(λ)) -# @info "Plot spectrum..." -# scatter(real.(λ), imag.(λ), label="central flux") -wololo +# @info "max real part" maximum(real.(λ)) +# # @info "Plot spectrum..." +# # scatter(real.(λ), imag.(λ), label="central flux") +# wololo ############################################################################### # ODE solvers, callbacks etc. diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index ced8a5221a2..db1b65e207f 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -23,10 +23,10 @@ initial_condition = initial_condition_weak_blast_wave # volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # surface_flux = flux_ranocha -polydeg = 11 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) -# volume_integral = VolumeIntegralWeakFormProjection() -volume_integral = VolumeIntegralWeakForm() +polydeg = 3 +basis = GaussLegendreBasis(polydeg; polydeg_projection = 1 * polydeg, polydeg_cutoff = 3) +volume_integral = VolumeIntegralWeakFormProjection() +#volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-2.0, -2.0) diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl index 689763fa706..8835f5cccf3 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability.jl @@ -32,15 +32,15 @@ function initial_condition_kelvin_helmholtz_instability(x, t, end initial_condition = initial_condition_kelvin_helmholtz_instability -surface_flux = flux_lax_friedrichs -volume_flux = flux_central +surface_flux = flux_hllc +volume_flux = flux_ranocha polydeg = 3 basis = LobattoLegendreBasis(polydeg) indicator_sc = IndicatorHennemannGassner(equations, basis, - alpha_max = 0.000, - alpha_min = 0.0000, + alpha_max = 0.2, + alpha_min = 0.0001, alpha_smooth = true, - variable = density_pressure) + variable = Trixi.density) volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg = volume_flux, diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl index 77c9572f717..5906f094e9b 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl @@ -32,10 +32,10 @@ function initial_condition_kelvin_helmholtz_instability(x, t, end initial_condition = initial_condition_kelvin_helmholtz_instability -surface_flux = flux_lax_friedrichs +surface_flux = FluxLMARS(1.7) # polydeg = 3 -polydeg = 10 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) +polydeg = 3 +basis = GaussLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) #indicator_sc = IndicatorHennemannGassner(equations, basis, # alpha_max = 0.000, # alpha_min = 0.0000, @@ -45,7 +45,7 @@ basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_ #volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; # volume_flux_dg = volume_flux, # volume_flux_fv = surface_flux) -# volume_integral = VolumeIntegralWeakFormProjection() +#volume_integral = VolumeIntegralWeakFormProjection() volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) @@ -75,7 +75,10 @@ save_solution = SaveSolutionCallback(interval=1000, solution_variables = cons2prim) # stepsize_callback = StepsizeCallback(cfl = 0.5) -stepsize_callback = StepsizeCallback(cfl = 1.5) +stepsize_callback = StepsizeCallback(cfl = 0.5) + +#stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-4, 5.0e-4), +# variables = (Trixi.density, pressure)) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, @@ -85,7 +88,9 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), - dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback +#sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +#sol = solve(ode, SSPRK43(stage_limiter!), +sol = solve(ode, SSPRK43(), + dt = 1.0e-3, # 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 diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl index 9e9fb45e7d1..81e55f09848 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl @@ -32,7 +32,7 @@ function initial_condition_kelvin_helmholtz_instability(x, t, end initial_condition = initial_condition_kelvin_helmholtz_instability -surface_flux = flux_lax_friedrichs +surface_flux = flux_hllc volume_flux = flux_ranocha polydeg = 3 basis = LobattoLegendreBasis(polydeg) @@ -48,7 +48,7 @@ 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 = 6, n_cells_max = 100_000) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index e06318a500e..05527f4a69a 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -5,7 +5,7 @@ using Trixi ############################################################################### # semidiscretization of the compressible Euler equations -epsilon_relaxation = 5.0e-6 +epsilon_relaxation = 1.0e-6 a1 = a2 = a3 = a4 = 30.0 b1 = b2 = b3 = b4 = 30.0 @@ -29,9 +29,9 @@ end #initial_condition = initial_condition_constant initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability) #initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave) -polydeg = 5 -polydeg_cutoff = 3 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 1 * polydeg, polydeg_cutoff = polydeg_cutoff) +polydeg = 3 +polydeg_cutoff = 3 +basis = GaussLegendreBasis(polydeg; polydeg_projection = polydeg, polydeg_cutoff = polydeg_cutoff) # solver = DGSEM(basis, Trixi.flux_upwind) #solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakFormProjection()) solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakForm()) diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl index c1108128c92..56bac978d0a 100644 --- a/src/visualization/utilities.jl +++ b/src/visualization/utilities.jl @@ -1179,6 +1179,7 @@ function unstructured2structured(unstructured_data, normalized_coordinates, # Get node coordinates for DG locations on reference element nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in) + #nodes_in, _ = gauss_nodes_weights(n_nodes_in) # Calculate interpolation vandermonde matrices for each level max_level = length(nvisnodes_per_level) - 1 From 20c60d9095146e128c781500a55a0953b4d34b8a Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 7 Jun 2024 20:08:38 +0200 Subject: [PATCH 20/20] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- src/equations/jin_xin_compressible_euler_2d.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/equations/jin_xin_compressible_euler_2d.jl b/src/equations/jin_xin_compressible_euler_2d.jl index 7636ccc7314..da43b3f6b0b 100644 --- a/src/equations/jin_xin_compressible_euler_2d.jl +++ b/src/equations/jin_xin_compressible_euler_2d.jl @@ -166,7 +166,6 @@ end 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))..., @ranocha: is this correct? (sqrt_velocities .* (f_rr_base - f_ll_base))..., zero(u_ll_base)...) else # orientation == 2 @@ -175,7 +174,6 @@ end 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))..., # @ranocha: is this correct? (sqrt_velocities .* (f_rr_base - f_ll_base))...) end @@ -186,7 +184,6 @@ end @inline function max_abs_speeds(u, equations::JinXinEquations{2}) return ntuple(Val(ndims(equations))) do n - # maximum(equations.sqrt_velocities_inv[n]) @ranocha: is this correct? maximum(equations.sqrt_velocities[n]) end end