diff --git a/Project.toml b/Project.toml index bcc1224a96d..3b0c65a22d1 100644 --- a/Project.toml +++ b/Project.toml @@ -14,12 +14,16 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" diff --git a/examples/2d/elixir_euler_kelvin_helmholtz_instability.jl b/examples/2d/elixir_euler_kelvin_helmholtz_instability.jl index 74194f2f302..9db32c27b60 100644 --- a/examples/2d/elixir_euler_kelvin_helmholtz_instability.jl +++ b/examples/2d/elixir_euler_kelvin_helmholtz_instability.jl @@ -16,13 +16,13 @@ volume_flux = flux_chandrashekar polydeg = 3 basis = LobattoLegendreBasis(polydeg) indicator_sc = IndicatorHennemannGassner(equations, basis, - alpha_max=0.002, - alpha_min=0.0001, - alpha_smooth=true, - variable=density_pressure) + alpha_max=0.002, + alpha_min=0.0001, + alpha_smooth=true, + variable=density_pressure) volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg=volume_flux, - volume_flux_fv=surface_flux) + volume_flux_dg=volume_flux, + volume_flux_fv=surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-0.5, -0.5) @@ -53,7 +53,7 @@ save_solution = SaveSolutionCallback(interval=20, stepsize_callback = StepsizeCallback(cfl=1.3) callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, + analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/2d/elixir_euler_kelvin_helmholtz_instability2.jl b/examples/2d/elixir_euler_kelvin_helmholtz_instability2.jl new file mode 100644 index 00000000000..70dffa185c2 --- /dev/null +++ b/examples/2d/elixir_euler_kelvin_helmholtz_instability2.jl @@ -0,0 +1,137 @@ + +using Random: seed! +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +seed!(0) +function initial_condition_khi2(x, t, equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # typical resolution 128^2, 256^2 + # domain size is [-1,+1]^2 + slope = 15 + 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(1 * pi * x[1]) #+ 0.05 * sin(4 * pi * x[1]) + p = 1.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end + + + +initial_condition = initial_condition_khi2 + +#surface_flux = flux_lax_friedrichs +surface_flux = flux_hllc +#surface_flux = flux_hll +#surface_flux = flux_shima_etal +volume_flux = flux_ranocha +#volume_flux = flux_shima_etal +#volume_flux = flux_central +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +#indicator_sc = IndicatorHennemannGassner(equations, basis, +# alpha_max=0.002, +# alpha_min=0.0001, +# alpha_smooth=true, +# variable=density_pressure) +#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; +# volume_flux_dg=volume_flux, +# volume_flux_fv=surface_flux) + +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +#volume_integral = VolumeIntegralWeakForm() +solver = DGSEM(basis, surface_flux, volume_integral) + +#surface_flux = flux_hllc +#volume_flux = flux_ranocha +#solver = DGSEM(7, surface_flux, VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux))) + +#surface_flux = flux_hllc +#volume_flux = flux_ranocha +#volume_integral = VolumeIntegralFluxComparison(flux_central, volume_flux) +#solver = DGSEM(15, surface_flux, volume_integral) + +coordinates_min = (-1.0, -1.0) +coordinates_max = ( 1.0, 1.0) + +#restart_filename = joinpath("out", "solution_010700.h5") +#mesh = load_mesh(restart_filename, n_cells_max=800_000) + +println("start mesh") +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=6, + n_cells_max=800_000) +println("end mesh") +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +println("end semi") + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.7) +#tspan = (load_time(restart_filename), 3.7) +ode = semidiscretize(semi, tspan) +#ode = semidiscretize(semi, tspan, restart_filename) +println("end ode") + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true, + extra_analysis_errors=(:conservation_error,), + extra_analysis_integrals=(entropy, energy_total, + energy_kinetic, energy_internal)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=1000, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +save_restart = SaveRestartCallback(interval=1000) + +stepsize_callback = StepsizeCallback(cfl=0.2) + +#amr_indicator = IndicatorLöhner(semi, +# variable=Trixi.density) +#amr_controller = ControllerThreeLevel(semi, amr_indicator, +# base_level=1, +# med_level =0, med_threshold=0.1, # med_level = current level +# max_level =6, max_threshold=0.5) +#amr_callback = AMRCallback(semi, amr_controller, +# interval=3, +# adapt_initial_condition=false, +# adapt_initial_condition_only_refine=true) + +limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-5, 5.0e-5), + variables=(Trixi.density, pressure)) +stage_limiter! = limiter! +step_limiter! = limiter! + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + save_restart, +# amr_callback, +# stepsize_callback + ) + +println("end callbacks") + +############################################################################### +# run the simulation + +#sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), +#sol = solve(ode, SSPRK33(stage_limiter!, step_limiter!), +sol = solve(ode, SSPRK43(stage_limiter!, step_limiter!), +#sol = solve(ode, SSPRK33(), + #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/examples/2d/elixir_euler_kelvin_helmholtz_instability3.jl b/examples/2d/elixir_euler_kelvin_helmholtz_instability3.jl new file mode 100644 index 00000000000..ad727ef2d3e --- /dev/null +++ b/examples/2d/elixir_euler_kelvin_helmholtz_instability3.jl @@ -0,0 +1,160 @@ + +using Random: seed! +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 5/3 +equations = CompressibleEulerEquations2D(gamma) + +seed!(0) +function initial_condition_khi3(x, t, equations::CompressibleEulerEquations2D) + # from https://arxiv.org/pdf/2001.01927.pdf + # typical resolution 128^2, 256^2 + # domain size is [0,1]^2 + L = 0.025 + U1 = 0.5 + U2 =-0.5 + Um = 0.5 * (U1-U2) + rho1 = 1.0 + rho2 = 2.0 + rhom = 0.5 * (rho1 - rho2) + + v2 = 0.01 * sin(4 * pi * x[1]) + p = 2.5 + + if (0 <= x[2] <= 0.25) + rho = rho1 - rhom * exp((x[2] - 0.25) / L) + v1 = U1 - Um * exp((x[2] - 0.25) / L) + elseif (0.25 < x[2] <= 0.5) + rho = rho2 + rhom * exp((-x[2] + 0.25) / L) + v1 = U2 + Um * exp((-x[2] + 0.25) / L) + elseif (0.5 < x[2] <= 0.75) + rho = rho2 + rhom * exp((x[2] - 0.75) / L) + v1 = U2 + Um * exp((x[2] - 0.75) / L) + elseif (0.75 < x[2] <= 1.0) + rho = rho1 - rhom * exp((-x[2] + 0.75) / L) + v1 = U1 - Um * exp((-x[2] + 0.75) / L) + else + println("shyte", x[2]) + end + if (rho < 0.0) + println("density", x[2], rho) + end + if (p < 0.0) + println("pressure") + end + return prim2cons(SVector(rho, v1, v2, p), equations) +end + + + +initial_condition = initial_condition_khi3 + +#surface_flux = flux_lax_friedrichs +surface_flux = flux_hllc +#surface_flux = flux_hll +#surface_flux = flux_shima_etal +volume_flux = flux_ranocha +#volume_flux = flux_shima_etal +#volume_flux = flux_central +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +#indicator_sc = IndicatorHennemannGassner(equations, basis, +# alpha_max=0.002, +# alpha_min=0.0001, +# alpha_smooth=true, +# variable=density_pressure) +#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; +# volume_flux_dg=volume_flux, +# volume_flux_fv=surface_flux) + +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +#volume_integral = VolumeIntegralWeakForm() +solver = DGSEM(basis, surface_flux, volume_integral) + +#surface_flux = flux_hllc +#volume_flux = flux_ranocha +#solver = DGSEM(7, surface_flux, VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux))) + +#surface_flux = flux_hllc +#volume_flux = flux_ranocha +#volume_integral = VolumeIntegralFluxComparison(flux_central, volume_flux) +#solver = DGSEM(15, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = ( 1.0, 1.0) + +#restart_filename = joinpath("out", "solution_010700.h5") +#mesh = load_mesh(restart_filename, n_cells_max=800_000) + +println("start mesh") +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=7, + n_cells_max=4_000_000) +println("end mesh") +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.7) +#tspan = (load_time(restart_filename), 3.7) +ode = semidiscretize(semi, tspan) +#ode = semidiscretize(semi, tspan, restart_filename) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true, + extra_analysis_errors=(:conservation_error,), + extra_analysis_integrals=(entropy, energy_total, + energy_kinetic, energy_internal)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=1000, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +save_restart = SaveRestartCallback(interval=1000) + +stepsize_callback = StepsizeCallback(cfl=0.8) + +#amr_indicator = IndicatorLöhner(semi, +# variable=Trixi.density) +#amr_controller = ControllerThreeLevel(semi, amr_indicator, +# base_level=1, +# med_level =0, med_threshold=0.1, # med_level = current level +# max_level =6, max_threshold=0.5) +#amr_callback = AMRCallback(semi, amr_controller, +# interval=3, +# adapt_initial_condition=false, +# adapt_initial_condition_only_refine=true) + +limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-5, 5.0e-5), + variables=(Trixi.density, pressure)) +stage_limiter! = limiter! +step_limiter! = limiter! + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + save_restart, +# amr_callback, +# stepsize_callback + ) + + +############################################################################### +# run the simulation + +#sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), +#sol = solve(ode, SSPRK33(stage_limiter!, step_limiter!), +sol = solve(ode, SSPRK43(stage_limiter!, step_limiter!), +#sol = solve(ode, SSPRK33(), + #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/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison.jl b/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison.jl new file mode 100644 index 00000000000..da06658dc85 --- /dev/null +++ b/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison.jl @@ -0,0 +1,60 @@ + +# TODO: Needs tests +using Random: seed! +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +seed!(0) +initial_condition = initial_condition_khi + +surface_flux = flux_lax_friedrichs +volume_flux = flux_chandrashekar +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +solver = DGSEM(3, surface_flux, VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux))) + +coordinates_min = (-0.5, -0.5) +coordinates_max = ( 0.5, 0.5) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=5, + n_cells_max=100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=20, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl=1.3) + +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/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison_pure_ec.jl b/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison_pure_ec.jl new file mode 100644 index 00000000000..71b0cdfc366 --- /dev/null +++ b/examples/2d/elixir_euler_kelvin_helmholtz_instability_comparison_pure_ec.jl @@ -0,0 +1,62 @@ + +# TODO: Needs tests +using Random: seed! +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +seed!(0) +initial_condition = initial_condition_khi + +surface_flux = flux_lax_friedrichs +volume_flux = flux_chandrashekar +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +volume_integral = VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux), + (;alpha_ec=1, alpha_cen=0), variant=:parametric) +solver = DGSEM(3, surface_flux, volume_integral) + +coordinates_min = (-0.5, -0.5) +coordinates_max = ( 0.5, 0.5) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=5, + n_cells_max=100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=20, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl=1.3) + +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/examples/2d/elixir_euler_source_terms.jl b/examples/2d/elixir_euler_source_terms.jl index 6b36bba9322..b146ca4b8a1 100644 --- a/examples/2d/elixir_euler_source_terms.jl +++ b/examples/2d/elixir_euler_source_terms.jl @@ -8,7 +8,12 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) initial_condition = initial_condition_convergence_test -solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +# TODO: Change back +# solver = DGSEM(polydeg=3, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) +solver = DGSEM(polydeg=3, surface_flux = surface_flux, volume_integral = VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux))) coordinates_min = (0, 0) coordinates_max = (2, 2) diff --git a/examples/3d/elixir_euler_taylor_green_vortex.jl b/examples/3d/elixir_euler_taylor_green_vortex.jl index 33684c46afc..a4612f3956c 100644 --- a/examples/3d/elixir_euler_taylor_green_vortex.jl +++ b/examples/3d/elixir_euler_taylor_green_vortex.jl @@ -9,14 +9,44 @@ equations = CompressibleEulerEquations3D(1.4) initial_condition = initial_condition_taylor_green_vortex +# TODO: Change back + +# TGV run "VolumeIntegralLocalFluxComparison+secret" +# surface_flux = FluxComparedToCentral(flux_ranocha) +# volume_flux = flux_ranocha +# volume_integral = VolumeIntegralLocalFluxComparison(flux_central, volume_flux) + +# TGV run "VolumeIntegralLocalFluxComparison" +surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha -solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, - volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) +volume_integral = VolumeIntegralLocalFluxComparison(flux_central, volume_flux) + +# TGV run "VolumeIntegralFluxComparison" +# surface_flux = flux_lax_friedrichs +# volume_flux = flux_ranocha +# volume_integral = VolumeIntegralFluxComparison(flux_central, volume_flux) + +# TGV run "secret+LLF" +# surface_flux = flux_lax_friedrichs +# volume_flux = FluxComparedToCentral(flux_ranocha) +# volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +# TGV run "secret+secret" +# surface_flux = FluxComparedToCentral(flux_ranocha) +# volume_flux = FluxComparedToCentral(flux_ranocha) +# volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +# TGV run "mod1/VolumeIntegralLocalComparison" +# surface_flux = flux_lax_friedrichs +# volume_flux = flux_ranocha +# volume_integral = VolumeIntegralLocalComparison(VolumeIntegralFluxDifferencing(volume_flux)) + +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, volume_integral = volume_integral) coordinates_min = (-pi, -pi, -pi) coordinates_max = ( pi, pi, pi) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=3, + initial_refinement_level=4, n_cells_max=100_000) @@ -31,8 +61,10 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +analysis_interval = 20 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + save_analysis=true, + extra_analysis_integrals=(Trixi.energy_kinetic,)) alive_callback = AliveCallback(analysis_interval=analysis_interval) diff --git a/src/Trixi.jl b/src/Trixi.jl index 70366329562..3b801b71725 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -88,9 +88,9 @@ export AcousticPerturbationEquations2D, InviscidBurgersEquation1D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D -export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_godunov, +export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_godunov, flux_secret, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_kennedy_gruber, flux_shima_etal, - flux_ec, + flux_ec, FluxComparedToCentral, flux_left, flux_right, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, @@ -145,6 +145,8 @@ export TreeMesh, CurvedMesh export DG, DGSEM, LobattoLegendreBasis, VolumeIntegralWeakForm, VolumeIntegralFluxDifferencing, + VolumeIntegralLocalComparison, VolumeIntegralFluxComparison, + VolumeIntegralLocalFluxComparison, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, MortarL2 diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 4cb31860ddf..c3bd70e22dc 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -765,6 +765,20 @@ Entropy conserving two-point flux by return SVector(f1, f2, f3, f4) end +# TODO: Superseded by FluxComparedToCentral +@inline function flux_secret(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) + f_ranocha = flux_ranocha(u_ll, u_rr, orientation, equations) + f_central = flux_central(u_ll, u_rr, orientation, equations) + w_ll = cons2entropy(u_ll,equations) + w_rr = cons2entropy(u_rr,equations) + delta_entropy = dot((w_rr-w_ll),f_central-f_ranocha) + if (delta_entropy < 0.0) + f = f_central + else + f = 2*f_ranocha - f_central + end + return f +end """ flux_ranocha(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl index e69695e9abb..f13f002558c 100644 --- a/src/equations/inviscid_burgers_1d.jl +++ b/src/equations/inviscid_burgers_1d.jl @@ -72,6 +72,20 @@ end end +@inline function flux_godunov(u_ll, u_rr, orientation, equation::InviscidBurgersEquation1D) + u_L = u_ll[1] + u_R = u_rr[1] + + if u_L < 0 && u_R < 0 + return SVector(0.5 * u_R^2) + elseif u_L > 0 && u_R > 0 + return SVector(0.5 * u_L^2) + else + return zero(u_ll) + end +end + + # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation @inline function max_abs_speed_naive(u_ll, u_rr, orientation, equations::InviscidBurgersEquation1D) u_L = u_ll[1] @@ -86,11 +100,12 @@ end end -function flux_ec(u_ll, u_rr, orientation, equation::InviscidBurgersEquation1D) +@inline function flux_ec(u_ll, u_rr, orientation, equation::InviscidBurgersEquation1D) u_L = u_ll[1] u_R = u_rr[1] + one_sixth = one(typeof(u_L)) / 6 - return SVector((u_L^2 + u_L * u_R + u_R^2) / 6) + return SVector(one_sixth * (u_L^2 + u_L * u_R + u_R^2)) end diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index d9e17e51592..9fdc279dffa 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -18,6 +18,63 @@ DG method (except floating point errors). end +# TODO: Must be documented and might need a better name +struct FluxComparedToCentral{Numflux} + numflux::Numflux +end + +@inline function (f::FluxComparedToCentral{Numflux})(u_ll, u_rr, orientation, equations) where {Numflux} + + f_baseline = f.numflux(u_ll, u_rr, orientation, equations) + f_central = flux_central(u_ll, u_rr, orientation, equations) + w_ll = cons2entropy(u_ll, equations) + w_rr = cons2entropy(u_rr, equations) + # The local entropy production of a numerical flux at an interface is + # dot(w_rr - w_ll, numerical_flux) - (psi_rr - psi_ll), + # see Tadmor (1987). Since the flux potential is the same for both, we can + # omit that part. + delta_entropy = dot(w_rr - w_ll, f_central - f_baseline) + if delta_entropy < 0 + return f_central + else + return 2 * f_baseline - f_central + end +end + + +""" + flux_left(u_ll, u_rr, orientation, equations) = flux(u_ll, orientation, equations) + +This flux can be useful to construct upwind summation by parts (SBP) operators, cf. +- Mattsson (2017), + Diagonal norm upwind SBP operators. + [DOI: 10.1016/j.jcp.2017.01.042](https://doi.org/10.1016/j.jcp.2017.01.042). +- Ranocha, Mitsotakis, Ketcheson (2021), + A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations. + [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119). + [arXiv: 2006.14802 [math.NA]](https://arxiv.org/abs/2006.14802). +""" +@inline function flux_left(u_ll, u_rr, orientation, equations) + flux(u_ll, orientation, equations) +end + +""" + flux_right(u_ll, u_rr, orientation, equations) = flux(u_rr, orientation, equations) + +This flux can be useful to construct upwind summation by parts (SBP) operators, cf. +- Mattsson (2017), + Diagonal norm upwind SBP operators. + [DOI: 10.1016/j.jcp.2017.01.042](https://doi.org/10.1016/j.jcp.2017.01.042). +- Ranocha, Mitsotakis, Ketcheson (2021), + A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations. + [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119). + [arXiv: 2006.14802 [math.NA]](https://arxiv.org/abs/2006.14802). +""" +@inline function flux_right(u_ll, u_rr, orientation, equations) + flux(u_rr, orientation, equations) +end + + """ FluxPlusDissipation(numerical_flux, dissipation) diff --git a/src/solvers/dg/dg.jl b/src/solvers/dg/dg.jl index aa68b42dce0..5c1d22e1587 100644 --- a/src/solvers/dg/dg.jl +++ b/src/solvers/dg/dg.jl @@ -57,6 +57,86 @@ function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxDiffe end +# TODO: It would be really nice if we could convert volume integrals to a purely +# element-based approach. Then, we could have a new volume integral that +# contains two different methods and performs the blending of them locally. +# The choice coded below would then be the special case of blending the +# weak form with flux differencing, but we could also use other choices. +# TODO: This needs a better name... +struct VolumeIntegralLocalComparison{Variant, VolumeFlux, Parameters} <: AbstractVolumeIntegral + volume_integral_flux_differencing::VolumeIntegralFluxDifferencing{VolumeFlux} + parameters::Parameters +end + +function VolumeIntegralLocalComparison(volume_integral_flux_differencing, parameters=NamedTuple(); + variant=:default) + VolumeIntegralLocalComparison{Val{variant}, + typeof(volume_integral_flux_differencing.volume_flux), + typeof(parameters)}(volume_integral_flux_differencing, parameters) +end + +function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalComparison{Variant}) where Variant + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + setup = [ + "variant" => Variant, + "volume flux" => integral.volume_integral_flux_differencing.volume_flux, + "parameters" => integral.parameters + ] + summary_box(io, "VolumeIntegralLocalComparison", setup) + end +end + + +# TODO: It would be really nice if we could convert other volume integrals to +# the flux form required here easily... +# TODO: This needs a better name... +struct VolumeIntegralFluxComparison{VolumeFlux_a, VolumeFlux_b} <: AbstractVolumeIntegral + volume_flux_a::VolumeFlux_a + volume_flux_b::VolumeFlux_b +end + +function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxComparison) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + setup = [ + "volume flux a" => integral.volume_flux_a, + "volume flux b" => integral.volume_flux_b + ] + summary_box(io, "VolumeIntegralFluxComparison", setup) + end +end + + +# TODO: It would be really nice if we could convert other volume integrals to +# the flux form required here easily... +# TODO: This needs a better name... +struct VolumeIntegralLocalFluxComparison{VolumeFluxA, VolumeFluxB} <: AbstractVolumeIntegral + volume_flux_a::VolumeFluxA + volume_flux_b::VolumeFluxB +end + +function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralLocalFluxComparison) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + setup = [ + "volume flux a" => integral.volume_flux_a, + "volume flux b" => integral.volume_flux_b + ] + summary_box(io, "VolumeIntegralLocalFluxComparison", setup) + end +end + + """ VolumeIntegralShockCapturingHG diff --git a/src/solvers/dg/dg_1d.jl b/src/solvers/dg/dg_1d.jl index dd7410f778f..c5a0d48196a 100644 --- a/src/solvers/dg/dg_1d.jl +++ b/src/solvers/dg/dg_1d.jl @@ -22,6 +22,12 @@ function create_cache(mesh::TreeMesh{1}, equations::AbstractEquations{1}, cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...) + # TODO: This should be moved somewhere else but that would require more involved + # changes since we wouldn't know `nelements(dg, cache)`... + du_ec = Array{uEltype, 3}(undef, nvariables(equations), nnodes(dg), nelements(dg, cache)) + du_cen = similar(du_ec) + cache = (;cache..., du_ec, du_cen) + return cache end @@ -44,6 +50,40 @@ end # end +function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralLocalComparison, dg::DG, uEltype) + @unpack volume_integral_flux_differencing = volume_integral + + # TODO: We should allocate the additional temporary storage `du_ec, du_cen` here + cache = create_cache(mesh, equations, volume_integral_flux_differencing, dg, uEltype) + + return cache +end + + +function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 1}(undef, nnodes(dg)) for _ in 1:Threads.nthreads()] + fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + + return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded) +end + + +function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 1}(undef, nnodes(dg)) for _ in 1:Threads.nthreads()] + + return (; w_threaded) +end + + function create_cache(mesh::TreeMesh{1}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -119,24 +159,28 @@ end function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, - nonconservative_terms::Val{false}, equations, + nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) + @threaded for element in eachelement(dg, cache) + weak_form_kernel!(du, u, nonconservative_terms, equations, dg, cache, element) + end +end + +@inline function weak_form_kernel!(du::AbstractArray{<:Any,3}, u, + nonconservative_terms::Val{false}, equations, + dg::DGSEM, cache, element, alpha=true) @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) - for i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, element) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) - flux1 = flux(u_node, 1, equations) - for ii in eachnode(dg) - integral_contribution = derivative_dhat[ii, i] * flux1 - add_to_node_vars!(du, integral_contribution, equations, dg, ii, element) - end + flux1 = flux(u_node, 1, equations) + for ii in eachnode(dg) + integral_contribution = alpha * derivative_dhat[ii, i] * flux1 + add_to_node_vars!(du, integral_contribution, equations, dg, ii, element) end end - - return nothing end @@ -179,6 +223,156 @@ end end +function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, + nonconservative_terms, equations, + volume_integral::VolumeIntegralLocalComparison, + dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack du_ec, du_cen = cache + @unpack volume_flux = volume_integral.volume_integral_flux_differencing + + du_ec .= zero(eltype(du_ec)) + du_cen .= zero(eltype(du_cen)) + + @threaded for element in eachelement(dg, cache) + # compute volume integral with flux, and for comparison with central flux + split_form_kernel!(du_ec, u, nonconservative_terms, equations, volume_flux, dg, cache, element) + weak_form_kernel!(du_cen, u, nonconservative_terms, equations, dg, cache, element) + + # compute entropy production of both volume integrals + delta_entropy = zero(eltype(du)) + for i in eachnode(dg) + du_ec_node = get_node_vars(du_ec, equations, dg, i, element) + du_cen_node = get_node_vars(du_cen, equations, dg, i, element) + w_node = cons2entropy(get_node_vars(u, equations, dg, i, element), equations) + delta_entropy += weights[i] * dot(w_node, du_ec_node - du_cen_node) + end + if delta_entropy < 0 + for i in eachnode(dg) + du_cen_node = get_node_vars(du_cen, equations, dg, i, element) + add_to_node_vars!(du, du_cen_node, equations, dg, i, element) + end + else + for i in eachnode(dg) + du_ec_node = get_node_vars(du_ec, equations, dg, i, element) + du_cen_node = get_node_vars(du_cen, equations, dg, i, element) + add_to_node_vars!(du, 2*du_ec_node-du_cen_node, equations, dg, i, element) + end + end + end +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack inverse_weights = dg.basis + @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + fluxes_a = fluxes_a_threaded[Threads.threadid()] + fluxes_b = fluxes_b_threaded[Threads.threadid()] + + # compute entropy variables + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + w[i] = cons2entropy(u_node, equations) + end + + # compute local high-order fluxes + local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, element) + local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for i in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[i] + flux_b = fluxes_b[i] + # if dot(w[i+1] - w[i], flux_a - flux_b) <= 0 + # fluxes_a[i] = flux_a # flux_a is more dissipative + # # fluxes_a[i] = 2 * flux_a - flux_b # flux_a is more dissipative + # else + # fluxes_a[i] = flux_b # flux_b is more dissipative + # # fluxes_a[i] = 2 * flux_b - flux_a # flux_b is more dissipative + # end + b = dot(w[i+1] - w[i], flux_b - flux_a) + c = 1.0e-7 + # hyp = sqrt(b^2 + c^2) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[i] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, element) + for i in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), element) + end +end + +@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,3}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, element) + @unpack weights, derivative_split = dg.basis + + for i in eachindex(fluxes) + flux1 = zero(eltype(fluxes)) + for iip in i+1:nnodes(dg) + uiip = get_node_vars(u, equations, dg, iip, element) + for iim in 1:i + uiim = get_node_vars(u, equations, dg, iim, element) + flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations) + end + end + fluxes[i] = flux1 + end +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralLocalFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack derivative_split = dg.basis + @unpack w_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + + # compute entropy variables + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + w[i] = cons2entropy(u_node, equations) + end + + # perform flux-differencing in symmetric form + for i in eachnode(dg) + u_i = get_node_vars(u, equations, dg, i, element) + du_i = zero(u_i) + for ii in eachnode(dg) + u_ii = get_node_vars(u, equations, dg, ii, element) + flux_a = volume_flux_a(u_i, u_ii, 1, equations) + flux_b = volume_flux_b(u_i, u_ii, 1, equations) + d_i_ii = derivative_split[i, ii] + b = -sign(d_i_ii) * dot(w[i] - w[ii], flux_b - flux_a) + c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠 + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a)) + end + add_to_node_vars!(du, du_i, equations, dg, i, element) + end + end +end + + # TODO: Taal dimension agnostic function calc_volume_integral!(du::AbstractArray{<:Any,3}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, diff --git a/src/solvers/dg/dg_2d.jl b/src/solvers/dg/dg_2d.jl index 92a1b9a1407..9781bc60ef5 100644 --- a/src/solvers/dg/dg_2d.jl +++ b/src/solvers/dg/dg_2d.jl @@ -24,6 +24,13 @@ function create_cache(mesh::TreeMesh{2}, equations::AbstractEquations{2}, cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...) + # TODO: This should be moved somewhere else but that would require more involved + # changes since we wouldn't know `nelements(dg, cache)`... + # TODO: This does not work with AMR + du_ec = Array{uEltype, 4}(undef, nvariables(equations), nnodes(dg), nnodes(dg), nelements(dg, cache)) + du_cen = similar(du_ec) + cache = (;cache..., du_ec, du_cen) + return cache end @@ -58,6 +65,40 @@ function create_cache(mesh::TreeMesh{2}, nonconservative_terms::Val{true}, equat end +function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralLocalComparison, dg::DG, uEltype) + @unpack volume_integral_flux_differencing = volume_integral + + # TODO: We should allocate the additional temporary storage `du_ec, du_cen` here + cache = create_cache(mesh, equations, volume_integral_flux_differencing, dg, uEltype) + + return cache +end + + +function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 2}(undef, nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + + return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded) +end + + +function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 2}(undef, nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + + return (; w_threaded) +end + + function create_cache(mesh::TreeMesh{2}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -159,30 +200,34 @@ end function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, - nonconservative_terms::Val{false}, equations, + nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) + @threaded for element in eachelement(dg, cache) + weak_form_kernel!(du, u, nonconservative_terms, equations, dg, cache, element) + end +end + +@inline function weak_form_kernel!(du::AbstractArray{<:Any,4}, u, + nonconservative_terms::Val{false}, equations, + dg::DGSEM, cache, element, alpha=true) @unpack derivative_dhat = dg.basis - @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) + 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) - integral_contribution = derivative_dhat[ii, i] * flux1 - add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element) - end + flux1 = flux(u_node, 1, equations) + for ii in eachnode(dg) + integral_contribution = alpha * derivative_dhat[ii, i] * flux1 + add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, element) + end - flux2 = flux(u_node, 2, equations) - for jj in eachnode(dg) - integral_contribution = derivative_dhat[jj, j] * flux2 - add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element) - end + flux2 = flux(u_node, 2, equations) + for jj in eachnode(dg) + integral_contribution = alpha * derivative_dhat[jj, j] * flux2 + add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, element) end end - - return nothing end @@ -241,6 +286,7 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) + @threaded for element in eachelement(dg, cache) split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element) end @@ -321,6 +367,267 @@ end end +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, + nonconservative_terms, equations, + volume_integral::VolumeIntegralLocalComparison{Val{:default}}, + dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack du_ec, du_cen = cache + @unpack volume_flux = volume_integral.volume_integral_flux_differencing + + du_ec .= zero(eltype(du_ec)) + du_cen .= zero(eltype(du_cen)) + + c = 1.0e-12 + + @threaded for element in eachelement(dg, cache) + # compute volume integral with flux, and for comparison with central flux + split_form_kernel!(du_ec, u, nonconservative_terms, equations, volume_flux, dg, cache, element) + weak_form_kernel!(du_cen, u, nonconservative_terms, equations, dg, cache, element) + + # compute entropy production of both volume integrals + delta_entropy = zero(eltype(du)) + for j in eachnode(dg), i in eachnode(dg) + du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element) + w_node = cons2entropy(get_node_vars(u, equations, dg, i, j, element), equations) + delta_entropy += weights[i] * weights[j] * dot(w_node, du_ec_node - du_cen_node) + end + b = - delta_entropy + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + δ = (hyp - b) / hyp # add anti-dissipation as dissipation + #δ = (hyp - b) / 2hyp # just use the more dissipative flux + for j in eachnode(dg), i in eachnode(dg) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element) + du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element) + add_to_node_vars!(du, du_cen_node + δ * (du_ec_node - du_cen_node) , equations, dg, i, j, element) + end + #if delta_entropy < 0 + # for j in eachnode(dg), i in eachnode(dg) + # du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element) + # add_to_node_vars!(du, du_cen_node, equations, dg, i, j, element) + # end + #else + # for j in eachnode(dg), i in eachnode(dg) + # du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element) + # du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element) + # add_to_node_vars!(du, 2*du_ec_node-du_cen_node, equations, dg, i, j, element) + # end + #end + end +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, + nonconservative_terms, equations, + volume_integral::VolumeIntegralLocalComparison{Val{:parametric}}, + dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack du_ec, du_cen = cache + @unpack volume_flux = volume_integral.volume_integral_flux_differencing + @unpack parameters = volume_integral + @unpack alpha_ec = parameters + @unpack alpha_cen = parameters + + du_ec .= zero(eltype(du_ec)) + du_cen .= zero(eltype(du_cen)) + + @threaded for element in eachelement(dg, cache) + # compute volume integral with flux, and for comparison with central flux + split_form_kernel!(du_ec, u, nonconservative_terms, equations, volume_flux, dg, cache, element) + weak_form_kernel!(du_cen, u, nonconservative_terms, equations, dg, cache, element) + + # compute entropy production of both volume integrals + delta_entropy = zero(eltype(du)) + for j in eachnode(dg), i in eachnode(dg) + du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element) + w_node = cons2entropy(get_node_vars(u, equations, dg, i, j, element), equations) + delta_entropy += weights[i] * weights[j] * dot(w_node, du_ec_node - du_cen_node) + end + if delta_entropy < 0 + for j in eachnode(dg), i in eachnode(dg) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element) + add_to_node_vars!(du, du_cen_node, equations, dg, i, j, element) + end + else + for j in eachnode(dg), i in eachnode(dg) + du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element) + add_to_node_vars!(du, alpha_ec*du_ec_node + alpha_cen * du_cen_node, + equations, dg, i, j, element) + end + end + end +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack inverse_weights = dg.basis + @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + fluxes_a = fluxes_a_threaded[Threads.threadid()] + fluxes_b = fluxes_b_threaded[Threads.threadid()] + c = 1.0e-12 + + # compute entropy variables + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + w[i,j] = cons2entropy(u_node, equations) + end + + # x + for j in eachnode(dg) + # compute local high-order fluxes + local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, j, element) + local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, j, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for i in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[i] + flux_b = fluxes_b[i] + b = dot(w[i+1,j] - w[i,j], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[i] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, j, element) + for i in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, j, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), j, element) + end + + # y + for i in eachnode(dg) + # compute local high-order fluxes + local_fluxes_2!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, element) + local_fluxes_2!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for j in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[j] + flux_b = fluxes_b[j] + b = dot(w[i,j+1] - w[i,j], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[j] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, 1, element) + for j in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[j] * (fluxes_a[j] - fluxes_a[j-1]), equations, dg, i, j, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, nnodes(dg), element) + end + + end +end + +@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,4}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, j, element) + @unpack weights, derivative_split = dg.basis + + for i in eachindex(fluxes) + flux1 = zero(eltype(fluxes)) + for iip in i+1:nnodes(dg) + uiip = get_node_vars(u, equations, dg, iip, j, element) + for iim in 1:i + uiim = get_node_vars(u, equations, dg, iim, j, element) + flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations) + end + end + fluxes[i] = flux1 + end +end + +@inline function local_fluxes_2!(fluxes, u::AbstractArray{<:Any,4}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, i, element) + @unpack weights, derivative_split = dg.basis + + for j in eachindex(fluxes) + flux2 = zero(eltype(fluxes)) + for jjp in j+1:nnodes(dg) + ujjp = get_node_vars(u, equations, dg, i, jjp, element) + for jjm in 1:j + ujjm = get_node_vars(u, equations, dg, i, jjm, element) + flux2 += weights[jjm] * derivative_split[jjm,jjp] * volume_flux(ujjm, ujjp, 2, equations) + end + end + fluxes[j] = flux2 + end +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralLocalFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack derivative_split = dg.basis + @unpack w_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + + # compute entropy variables + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + w[i,j] = cons2entropy(u_node, equations) + end + + # perform flux-differencing in symmetric form + for j in eachnode(dg), i in eachnode(dg) + u_i = get_node_vars(u, equations, dg, i, j, element) + du_i = zero(u_i) + c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠 + + # x + for ii in eachnode(dg) + u_ii = get_node_vars(u, equations, dg, ii, j, element) + flux_a = volume_flux_a(u_i, u_ii, 1, equations) + flux_b = volume_flux_b(u_i, u_ii, 1, equations) + d_i_ii = derivative_split[i, ii] + b = -sign(d_i_ii) * dot(w[i,j] - w[ii,j], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a)) + end + + # y + for jj in eachnode(dg) + u_jj = get_node_vars(u, equations, dg, i, jj, element) + flux_a = volume_flux_a(u_i, u_jj, 2, equations) + flux_b = volume_flux_b(u_i, u_jj, 2, equations) + d_j_jj = derivative_split[j, jj] + b = -sign(d_j_jj) * dot(w[i,j] - w[i,jj], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_j_jj * (flux_a + δ * (flux_b - flux_a)) + end + + add_to_node_vars!(du, du_i, equations, dg, i, j, element) + end + + end +end + + # TODO: Taal dimension agnostic function calc_volume_integral!(du::AbstractArray{<:Any,4}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, diff --git a/src/solvers/dg/dg_3d.jl b/src/solvers/dg/dg_3d.jl index ac5d7cbd698..3303c5f34a9 100644 --- a/src/solvers/dg/dg_3d.jl +++ b/src/solvers/dg/dg_3d.jl @@ -24,6 +24,13 @@ function create_cache(mesh::TreeMesh{3}, equations::AbstractEquations{3}, cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...) + # TODO: This should be moved somewhere else but that would require more involved + # changes since we wouldn't know `nelements(dg, cache)`... + # TODO: This does not work with AMR + du_ec = Array{uEltype, 5}(undef, nvariables(equations), nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache)) + du_cen = similar(du_ec) + cache = (;cache..., du_ec, du_cen) + return cache end @@ -77,6 +84,40 @@ function create_cache(mesh::TreeMesh{3}, nonconservative_terms::Val{true}, equat end +function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralLocalComparison, dg::DG, uEltype) + @unpack volume_integral_flux_differencing = volume_integral + + # TODO: We should allocate the additional temporary storage `du_ec, du_cen` here + cache = create_cache(mesh, equations, volume_integral_flux_differencing, dg, uEltype) + + return cache +end + + +function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 3}(undef, nnodes(dg), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + fluxes_a_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + fluxes_b_threaded = [Vector{FluxType}(undef, nnodes(dg) - 1) for _ in 1:Threads.nthreads()] + + return (; w_threaded, fluxes_a_threaded, fluxes_b_threaded) +end + + +function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralLocalFluxComparison, dg::DG, uEltype) + + have_nonconservative_terms(equations) !== Val(false) && error("Are you kidding me?") + + FluxType = SVector{nvariables(equations), uEltype} + w_threaded = [Array{FluxType, 3}(undef, nnodes(dg), nnodes(dg), nnodes(dg)) for _ in 1:Threads.nthreads()] + + return (; w_threaded) +end + + function create_cache(mesh::TreeMesh{3}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -189,36 +230,40 @@ end function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, - nonconservative_terms::Val{false}, equations, + nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) + @threaded for element in eachelement(dg, cache) + weak_form_kernel!(du, u, nonconservative_terms, equations, dg, cache, element) + end +end + +@inline function weak_form_kernel!(du::AbstractArray{<:Any,5}, u, + nonconservative_terms::Val{false}, equations, + dg::DGSEM, cache, element, alpha=true) @unpack derivative_dhat = dg.basis - @threaded for element in eachelement(dg, cache) - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, j, k, element) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) - flux1 = flux(u_node, 1, equations) - for ii in eachnode(dg) - integral_contribution = derivative_dhat[ii, i] * flux1 - add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, k, element) - end + flux1 = flux(u_node, 1, equations) + for ii in eachnode(dg) + integral_contribution = alpha * derivative_dhat[ii, i] * flux1 + add_to_node_vars!(du, integral_contribution, equations, dg, ii, j, k, element) + end - flux2 = flux(u_node, 2, equations) - for jj in eachnode(dg) - integral_contribution = derivative_dhat[jj, j] * flux2 - add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, k, element) - end + flux2 = flux(u_node, 2, equations) + for jj in eachnode(dg) + integral_contribution = alpha * derivative_dhat[jj, j] * flux2 + add_to_node_vars!(du, integral_contribution, equations, dg, i, jj, k, element) + end - flux3 = flux(u_node, 3, equations) - for kk in eachnode(dg) - integral_contribution = derivative_dhat[kk, k] * flux3 - add_to_node_vars!(du, integral_contribution, equations, dg, i, j, kk, element) - end + flux3 = flux(u_node, 3, equations) + for kk in eachnode(dg) + integral_contribution = alpha * derivative_dhat[kk, k] * flux3 + add_to_node_vars!(du, integral_contribution, equations, dg, i, j, kk, element) end end - - return nothing end @@ -288,6 +333,7 @@ function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) + @threaded for element in eachelement(dg, cache) split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element) end @@ -385,6 +431,268 @@ end end +function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, + nonconservative_terms, equations, + volume_integral::VolumeIntegralLocalComparison, + dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack du_ec, du_cen = cache + @unpack volume_flux = volume_integral.volume_integral_flux_differencing + + du_ec .= zero(eltype(du_ec)) + du_cen .= zero(eltype(du_cen)) + + @threaded for element in eachelement(dg, cache) + # compute volume integral with flux, and for comparison with central flux + split_form_kernel!(du_ec, u, nonconservative_terms, equations, volume_flux, dg, cache, element) + weak_form_kernel!(du_cen, u, nonconservative_terms, equations, dg, cache, element) + + # compute entropy production of both volume integrals + delta_entropy = zero(eltype(du)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + du_ec_node = get_node_vars(du_ec, equations, dg, i, j, k, element) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, k, element) + w_node = cons2entropy(get_node_vars(u, equations, dg, i, j, k, element), equations) + delta_entropy += weights[i] * weights[j] * weights[k] * dot(w_node, du_ec_node - du_cen_node) + end + if delta_entropy < 0 + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, k, element) + add_to_node_vars!(du, du_cen_node, equations, dg, i, j, k, element) + end + else + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + du_ec_node = get_node_vars(du_ec, equations, dg, i, j, k, element) + du_cen_node = get_node_vars(du_cen, equations, dg, i, j, k, element) + add_to_node_vars!(du, 2*du_ec_node-du_cen_node, equations, dg, i, j, k, element) + end + end + end +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack inverse_weights = dg.basis + @unpack w_threaded, fluxes_a_threaded, fluxes_b_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + fluxes_a = fluxes_a_threaded[Threads.threadid()] + fluxes_b = fluxes_b_threaded[Threads.threadid()] + c = 1.0e-7 + + # compute entropy variables + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + w[i,j,k] = cons2entropy(u_node, equations) + end + + # x + for k in eachnode(dg), j in eachnode(dg) + # compute local high-order fluxes + local_fluxes_1!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, j, k, element) + local_fluxes_1!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, j, k, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for i in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[i] + flux_b = fluxes_b[i] + b = dot(w[i+1,j,k] - w[i,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[i] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, 1, j, k, element) + for i in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[i] * (fluxes_a[i] - fluxes_a[i-1]), equations, dg, i, j, k, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, nnodes(dg), j, k, element) + end + + # y + for k in eachnode(dg), i in eachnode(dg) + # compute local high-order fluxes + local_fluxes_2!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, k, element) + local_fluxes_2!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, k, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for j in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[j] + flux_b = fluxes_b[j] + b = dot(w[i,j+1,k] - w[i,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[j] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, 1, k, element) + for j in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[j] * (fluxes_a[j] - fluxes_a[j-1]), equations, dg, i, j, k, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, nnodes(dg), k, element) + end + + # z + for j in eachnode(dg), i in eachnode(dg) + # compute local high-order fluxes + local_fluxes_3!(fluxes_a, u, nonconservative_terms, equations, volume_flux_a, dg, i, j, element) + local_fluxes_3!(fluxes_b, u, nonconservative_terms, equations, volume_flux_b, dg, i, j, element) + + # compare entropy production of both fluxes and choose the more dissipative one + for k in eachindex(fluxes_a, fluxes_b) + flux_a = fluxes_a[k] + flux_b = fluxes_b[k] + b = dot(w[i,j,k+1] - w[i,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + fluxes_a[k] = flux_a + δ * (flux_b - flux_a) + end + + # update volume contribution in locally conservative form + add_to_node_vars!(du, inverse_weights[1] * fluxes_a[1], equations, dg, i, j, 1, element) + for k in 2:nnodes(dg)-1 + add_to_node_vars!(du, inverse_weights[k] * (fluxes_a[k] - fluxes_a[k-1]), equations, dg, i, j, k, element) + end + add_to_node_vars!(du, -inverse_weights[end] * fluxes_a[end], equations, dg, i, j, nnodes(dg), element) + end + + end +end + +@inline function local_fluxes_1!(fluxes, u::AbstractArray{<:Any,5}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, j, k, element) + @unpack weights, derivative_split = dg.basis + + for i in eachindex(fluxes) + flux1 = zero(eltype(fluxes)) + for iip in i+1:nnodes(dg) + uiip = get_node_vars(u, equations, dg, iip, j, k, element) + for iim in 1:i + uiim = get_node_vars(u, equations, dg, iim, j, k, element) + flux1 += weights[iim] * derivative_split[iim,iip] * volume_flux(uiim, uiip, 1, equations) + end + end + fluxes[i] = flux1 + end +end + +@inline function local_fluxes_2!(fluxes, u::AbstractArray{<:Any,5}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, i, k, element) + @unpack weights, derivative_split = dg.basis + + for j in eachindex(fluxes) + flux2 = zero(eltype(fluxes)) + for jjp in j+1:nnodes(dg) + ujjp = get_node_vars(u, equations, dg, i, jjp, k, element) + for jjm in 1:j + ujjm = get_node_vars(u, equations, dg, i, jjm, k, element) + flux2 += weights[jjm] * derivative_split[jjm,jjp] * volume_flux(ujjm, ujjp, 2, equations) + end + end + fluxes[j] = flux2 + end +end + +@inline function local_fluxes_3!(fluxes, u::AbstractArray{<:Any,5}, + nonconservative_terms::Val{false}, equations, + volume_flux, dg::DGSEM, i, j, element) + @unpack weights, derivative_split = dg.basis + + for k in eachindex(fluxes) + flux3 = zero(eltype(fluxes)) + for kkp in k+1:nnodes(dg) + ukkp = get_node_vars(u, equations, dg, i, j, kkp, element) + for kkm in 1:k + ukkm = get_node_vars(u, equations, dg, i, j, kkm, element) + flux3 += weights[kkm] * derivative_split[kkm,kkp] * volume_flux(ukkm, ukkp, 3, equations) + end + end + fluxes[k] = flux3 + end +end + + +function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, + nonconservative_terms::Val{false}, equations, + volume_integral::VolumeIntegralLocalFluxComparison, + dg::DGSEM, cache) + @unpack volume_flux_a, volume_flux_b = volume_integral + @unpack derivative_split = dg.basis + @unpack w_threaded = cache + + @threaded for element in eachelement(dg, cache) + w = w_threaded[Threads.threadid()] + + # compute entropy variables + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + w[i,j,k] = cons2entropy(u_node, equations) + end + + # perform flux-differencing in symmetric form + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_i = get_node_vars(u, equations, dg, i, j, k, element) + du_i = zero(u_i) + c = 1.0e-4 # TODO: magic constant determining linear and nonlinear stability 😠 + + # x + for ii in eachnode(dg) + u_ii = get_node_vars(u, equations, dg, ii, j, k, element) + flux_a = volume_flux_a(u_i, u_ii, 1, equations) + flux_b = volume_flux_b(u_i, u_ii, 1, equations) + d_i_ii = derivative_split[i, ii] + b = -sign(d_i_ii) * dot(w[i,j,k] - w[ii,j,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_i_ii * (flux_a + δ * (flux_b - flux_a)) + end + + # y + for jj in eachnode(dg) + u_jj = get_node_vars(u, equations, dg, i, jj, k, element) + flux_a = volume_flux_a(u_i, u_jj, 2, equations) + flux_b = volume_flux_b(u_i, u_jj, 2, equations) + d_j_jj = derivative_split[j, jj] + b = -sign(d_j_jj) * dot(w[i,j,k] - w[i,jj,k], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_j_jj * (flux_a + δ * (flux_b - flux_a)) + end + + # z + for kk in eachnode(dg) + u_kk = get_node_vars(u, equations, dg, i, j, kk, element) + flux_a = volume_flux_a(u_i, u_kk, 3, equations) + flux_b = volume_flux_b(u_i, u_kk, 3, equations) + d_k_kk = derivative_split[k, kk] + b = -sign(d_k_kk) * dot(w[i,j,k] - w[i,j,kk], flux_b - flux_a) + hyp = hypot(b, c) # sqrt(b^2 + c^2) computed in a numerically stable way + # δ = (hyp - b) / hyp # add anti-dissipation as dissipation + δ = (hyp - b) / 2hyp # just use the more dissipative flux + du_i += d_k_kk * (flux_a + δ * (flux_b - flux_a)) + end + + add_to_node_vars!(du, du_i, equations, dg, i, j, k, element) + end + + end +end + + # TODO: Taal dimension agnostic function calc_volume_integral!(du::AbstractArray{<:Any,5}, u, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, diff --git a/test/test_examples_1d_burgers.jl b/test/test_examples_1d_burgers.jl index 2c75d58ce92..9c24bf79b0e 100644 --- a/test/test_examples_1d_burgers.jl +++ b/test/test_examples_1d_burgers.jl @@ -15,6 +15,13 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "1d") linf = [0.00016152468882624227]) end + @testset "elixir_burgers_basic.jl with flux_godunov" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_basic.jl"), + surface_flux = flux_godunov, + l2 = [2.967494855204592e-5], + linf = [0.0001615256989919711]) + end + @testset "elixir_burgers_linear_stability.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_linear_stability.jl"), l2 = [0.5660569881106876],