From 9acaa9033fc2d062ab0a7e6d219414cb9aa30eeb Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Thu, 9 Nov 2023 08:07:17 +0100 Subject: [PATCH 01/17] Add Trixi2Vtk.jl to downstream tests (#1722) --- .github/workflows/downstream.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml index d9e3a2d8a7e..b40d5d365c0 100644 --- a/.github/workflows/downstream.yml +++ b/.github/workflows/downstream.yml @@ -60,6 +60,7 @@ jobs: arch: - x64 package: + - Trixi2Vtk.jl - TrixiShallowWater.jl steps: - uses: actions/checkout@v4 From 63c625824c3d510611dcebba12a0319d7e9cad9a Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 10 Nov 2023 07:28:45 +0100 Subject: [PATCH 02/17] Fix bug (#1719) --- src/equations/compressible_euler_2d.jl | 2 +- test/test_p4est_2d.jl | 16 ++++++++-------- test/test_unit.jl | 4 ++-- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 6c8f2e1e848..a992f99eaf4 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1351,7 +1351,7 @@ of the numerical flux. v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) * inv_sum_sqrt_rho v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) * inv_sum_sqrt_rho v_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] - v_roe_mag = (v1_roe * normal_direction[1])^2 + (v2_roe * normal_direction[2])^2 + v_roe_mag = v1_roe^2 + v2_roe^2 H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) * inv_sum_sqrt_rho c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * v_roe_mag)) * norm_ diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 546e5bff8a6..07c6d02bbcd 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -203,16 +203,16 @@ end @trixi_testset "elixir_euler_sedov.jl (HLLE)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[ - 0.411541263324004, - 0.2558929632770186, - 0.2558929632770193, - 1.298715766843915, + 0.40853279043747015, + 0.25356771650524296, + 0.2535677165052422, + 1.2984601729572691, ], linf=[ - 1.3457201726152221, - 1.3138961427140758, - 1.313896142714079, - 6.293305112638921, + 1.3840909333784284, + 1.3077772519086124, + 1.3077772519086157, + 6.298798630968632, ], surface_flux=flux_hlle, tspan=(0.0, 0.3)) diff --git a/test/test_unit.jl b/test/test_unit.jl index a73dfab5504..29390161ebe 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1191,7 +1191,7 @@ end u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - flux_hll, FluxHLL(min_max_speed_davis)] + flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] for f_std in fluxes f_rot = FluxRotated(f_std) @@ -1215,7 +1215,7 @@ end SVector(1.5, -0.2, 0.1, 0.2, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, FluxLMARS(340), - flux_hll, FluxHLL(min_max_speed_davis)] + flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] for f_std in fluxes f_rot = FluxRotated(f_std) From 910df76691685040cd46ed661253f7348cb35344 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Fri, 10 Nov 2023 07:44:38 +0100 Subject: [PATCH 03/17] fix IndicatorMax with multi-threading on ARM (#1723) --- src/solvers/dgsem_tree/indicators_1d.jl | 3 ++- src/solvers/dgsem_tree/indicators_2d.jl | 3 ++- src/solvers/dgsem_tree/indicators_3d.jl | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 8b57348861c..40bfd1e98c7 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -289,6 +289,7 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 3}, kwargs...) @unpack alpha, indicator_threaded = indicator_max.cache resize!(alpha, nelements(dg, cache)) + indicator_variable = indicator_max.variable @threaded for element in eachelement(dg, cache) indicator = indicator_threaded[Threads.threadid()] @@ -296,7 +297,7 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 3}, # Calculate indicator variables at Gauss-Lobatto nodes for i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, element) - indicator[i] = indicator_max.variable(u_local, equations) + indicator[i] = indicator_variable(u_local, equations) end alpha[element] = maximum(indicator) diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 2f34e0eb661..da81b2a1d36 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -323,6 +323,7 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 4}, kwargs...) @unpack alpha, indicator_threaded = indicator_max.cache resize!(alpha, nelements(dg, cache)) + indicator_variable = indicator_max.variable @threaded for element in eachelement(dg, cache) indicator = indicator_threaded[Threads.threadid()] @@ -330,7 +331,7 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 4}, # Calculate indicator variables at Gauss-Lobatto nodes for j in eachnode(dg), i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, j, element) - indicator[i, j] = indicator_max.variable(u_local, equations) + indicator[i, j] = indicator_variable(u_local, equations) end alpha[element] = maximum(indicator) diff --git a/src/solvers/dgsem_tree/indicators_3d.jl b/src/solvers/dgsem_tree/indicators_3d.jl index 69041ed1298..40362889397 100644 --- a/src/solvers/dgsem_tree/indicators_3d.jl +++ b/src/solvers/dgsem_tree/indicators_3d.jl @@ -234,6 +234,7 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 5}, kwargs...) @unpack alpha, indicator_threaded = indicator_max.cache resize!(alpha, nelements(dg, cache)) + indicator_variable = indicator_max.variable @threaded for element in eachelement(dg, cache) indicator = indicator_threaded[Threads.threadid()] @@ -241,7 +242,7 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 5}, # Calculate indicator variables at Gauss-Lobatto nodes for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, j, k, element) - indicator[i, j, k] = indicator_max.variable(u_local, equations) + indicator[i, j, k] = indicator_variable(u_local, equations) end alpha[element] = maximum(indicator) From 9a00184c6ca656ce2fab4c8600bc837883a02b06 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 10 Nov 2023 07:46:35 +0100 Subject: [PATCH 04/17] Remove AMR CB from HLLE 2D Tree test (#1720) * remove amr cb * fmt --- test/test_tree_2d_euler.jl | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index db36cb7d79f..93490f8ae09 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -425,19 +425,22 @@ end @trixi_testset "elixir_euler_sedov_blast_wave.jl (HLLE)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"), l2=[ - 0.35267161504176747, - 0.17218309138797958, - 0.17218307467125854, - 0.6236143054619037, + 0.352405949321075, + 0.17207721487429464, + 0.17207721487433883, + 0.6263024434020885, ], linf=[ - 2.77484045816607, - 1.8281111268370718, - 1.8281110470490887, - 6.24263735888126, + 2.760997358628186, + 1.8279186132509326, + 1.8279186132502805, + 6.251573757093399, ], tspan=(0.0, 0.5), - surface_flux=flux_hlle) + callbacks=CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback), + surface_flux=flux_hlle), # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From cb8181f6e5ae2fa74054ab227b2e19a393edcee0 Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Fri, 10 Nov 2023 15:59:07 +0100 Subject: [PATCH 05/17] Fix the wet/dry test for SWE (#1727) * Adjust tolerance to fix wet_dry test * apply formatter --- test/test_tree_1d_shallowwater.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 658f178c941..7ec3089d33a 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -109,7 +109,9 @@ end 2.2447689894899726e-13, 1.9999999999999714, ], - tspan=(0.0, 0.25)) + tspan=(0.0, 0.25), + # Soften the tolerance as test results vary between different CPUs + atol=1000 * eps()) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 6158ba4cc319d7715d1b9419043c7ba55f15c473 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 11 Nov 2023 19:15:16 +0100 Subject: [PATCH 06/17] fix failing `DGMultiMesh` and Compressible Navier-Stokes convergence tests (#1728) (#1732) * fix failing test * more fixes * formatting * fix dropped part of source terms * fix p4est parabolic Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> --- .../p4est_2d_dgsem/elixir_navierstokes_convergence.jl | 1 + .../tree_2d_dgsem/elixir_navierstokes_convergence.jl | 11 ++++++----- test/test_parabolic_2d.jl | 3 ++- test/test_special_elixirs.jl | 6 ++++-- 4 files changed, 13 insertions(+), 8 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl index 7aa14e25750..54ec09d2be8 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl @@ -161,6 +161,7 @@ end v1_yy * v1 * mu_ - v2_xy * v1 * mu_ - v1_y * v1_y * mu_ - + v2_x * v1_y * mu_ - 4.0 / 3.0 * v2_yy * v2 * mu_ + 2.0 / 3.0 * v1_xy * v2 * mu_ - 4.0 / 3.0 * v2_y * v2_y * mu_ + diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl index 57558224854..b0c8678baad 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl @@ -160,15 +160,16 @@ end # stress tensor and temperature gradient terms from y-direction v1_yy * v1 * mu_ - v2_xy * v1 * mu_ - + v1_y * v1_y * mu_ - v2_x * v1_y * mu_ - 4.0 / 3.0 * v2_yy * v2 * mu_ + - 2.0 / 3.0 * v1_xy - + 2.0 / 3.0 * v1_xy * v2 * mu_ - 4.0 / 3.0 * v2_y * v2_y * mu_ + - 2.0 / 3.0 * v1_x * v2_y * mu_ - - - - + 2.0 / 3.0 * v1_x * v2_y * mu_ - + T_const * inv_rho_cubed * (p_yy * rho * rho - - 2.0 * p_y * rho * rho_y + - - + 2.0 * p_y * rho * rho_y + + 2.0 * p * rho_y * rho_y - p * rho * rho_yy) * mu_) return SVector(du1, du2, du3, du4) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 22a5a8b4e31..152ca52a6ca 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -16,7 +16,8 @@ isdir(outdir) && rm(outdir, recursive = true) dg = DGMulti(polydeg = 2, element_type = Quad(), approximation_type = Polynomial(), surface_integral = SurfaceIntegralWeakForm(flux_central), volume_integral = VolumeIntegralWeakForm()) - mesh = DGMultiMesh(dg, cells_per_dimension = (2, 2)) + cells_per_dimension = (2, 2) + mesh = DGMultiMesh(dg, cells_per_dimension) # test with polynomial initial condition x^2 * y # test if we recover the exact second derivative diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl index 4f42414ccbf..85671002ba6 100644 --- a/test/test_special_elixirs.jl +++ b/test/test_special_elixirs.jl @@ -203,7 +203,8 @@ end volume_integral = VolumeIntegralWeakForm()) # DGMultiMesh is on [-1, 1]^ndims by default - mesh = DGMultiMesh(solver, cells_per_dimension = (2, 2), + cells_per_dimension = (2, 2) + mesh = DGMultiMesh(solver, cells_per_dimension, periodicity = (true, true)) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, @@ -225,7 +226,8 @@ end volume_integral = VolumeIntegralFluxDifferencing(flux_central)) # DGMultiMesh is on [-1, 1]^ndims by default - mesh = DGMultiMesh(solver, cells_per_dimension = (2, 2), + cells_per_dimension = (2, 2) + mesh = DGMultiMesh(solver, cells_per_dimension, periodicity = (true, true)) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, From 53456098659a9b85fdd5c0a22184d10b88a3ba83 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Sat, 11 Nov 2023 14:10:55 -0600 Subject: [PATCH 07/17] Fix `flux` for quasi-1D SWE (#1731) * fix flux (remove non-conservative part) * add test * formatting --------- Co-authored-by: Hendrik Ranocha --- ...xir_shallowwater_quasi_1d_discontinuous.jl | 73 +++++++++++++++++++ src/equations/shallow_water_quasi_1d.jl | 9 +-- test/test_tree_1d_shallowwater.jl | 25 +++++++ 3 files changed, 101 insertions(+), 6 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl new file mode 100644 index 00000000000..76c04759389 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl @@ -0,0 +1,73 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the quasi 1d shallow water equations +# See Chan et al. https://doi.org/10.48550/arXiv.2307.12089 for details + +equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81) + +function initial_condition_discontinuity(x, t, equations::ShallowWaterEquationsQuasi1D) + H = 2 + 0.1 * exp(-25 * x[1]^2) + v = 0.0 + + if x[1] > 0 + b = 0.1 + a = 1.0 + else + b = 0.0 + a = 1.1 + end + + return prim2cons(SVector(H, v, b, a), equations) +end + +initial_condition = initial_condition_discontinuity + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = -0.5 +coordinates_max = 0.5 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = true) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.25) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 200, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index 217a764e173..d3935f0e75f 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/src/equations/shallow_water_quasi_1d.jl @@ -137,17 +137,14 @@ as defined in [`initial_condition_convergence_test`](@ref). return SVector(du1, du2, 0.0, 0.0) end -# Calculate 1D flux for a single point +# Calculate 1D conservative flux for a single point # Note, the bottom topography and channel width have no flux @inline function flux(u, orientation::Integer, equations::ShallowWaterEquationsQuasi1D) - a_h, a_h_v, _, a = u - h = waterheight(u, equations) + _, a_h_v, _, _ = u v = velocity(u, equations) - p = 0.5 * a * equations.gravity * h^2 - f1 = a_h_v - f2 = a_h_v * v + p + f2 = a_h_v * v return SVector(f1, f2, zero(eltype(u)), zero(eltype(u))) end diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 7ec3089d33a..fd69a0c1d0e 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -387,6 +387,31 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_shallowwater_quasi_1d_discontinuous.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_quasi_1d_discontinuous.jl"), + l2=[ + 0.02843233740533314, + 0.14083324483705398, + 0.0054554472558998, + 0.005455447255899814, + ], + linf=[ + 0.26095842440037487, + 0.45919004549253795, + 0.09999999999999983, + 0.10000000000000009, + ],) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end end # module From 6c70d18bd3c428662a2f80c7762099414d2ac4e3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 11 Nov 2023 21:11:48 +0100 Subject: [PATCH 08/17] set version to v0.5.48 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 497f4a1add7..311891e3b85 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.48-pre" +version = "0.5.48" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From fbf7d54c1975a031dded027cb92645f44ecbf67d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 11 Nov 2023 21:12:07 +0100 Subject: [PATCH 09/17] set development version to v0.5.49-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 311891e3b85..03973a8c8b4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.48" +version = "0.5.49-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 6d95941e919472e56d2fbace741fb63510f4b2a4 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Sat, 4 Nov 2023 13:05:10 +0100 Subject: [PATCH 10/17] Set version to v0.6.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 03973a8c8b4..e95e5347ee3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.49-pre" +version = "0.6.0" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 726f1a2fa96ef877624887a4bf10afcf75ff08e2 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 7 Nov 2023 07:22:09 +0100 Subject: [PATCH 11/17] Migrate neural network-based indicators to new repository (#1701) * Remove all neural network indicator stuff from `src/` * Migrate neural network tests * Migrate neural network examples * Migrate test dependencies * Update NEWS.md * Fix typo * Remove Requires.jl-based use of Flux.jl * Fix formatting * Add migration of indicators to section with breaking changes --------- Co-authored-by: Hendrik Ranocha --- NEWS.md | 19 + ...blast_wave_neuralnetwork_perssonperaire.jl | 109 ------ ...r_blast_wave_neuralnetwork_rayhesthaven.jl | 109 ------ ...ixir_euler_blast_wave_neuralnetwork_cnn.jl | 107 ------ ...blast_wave_neuralnetwork_perssonperaire.jl | 107 ------ ...r_blast_wave_neuralnetwork_rayhesthaven.jl | 110 ------ ...bility_amr_neuralnetwork_perssonperaire.jl | 134 ------- ...blast_wave_neuralnetwork_perssonperaire.jl | 126 ------- src/Trixi.jl | 8 +- src/solvers/dgsem_tree/indicators.jl | 195 ---------- src/solvers/dgsem_tree/indicators_1d.jl | 218 ----------- src/solvers/dgsem_tree/indicators_2d.jl | 351 ------------------ test/Project.toml | 4 - test/test_tree_1d_euler.jl | 20 - test/test_tree_2d_euler.jl | 114 ------ test/test_unit.jl | 8 - 16 files changed, 20 insertions(+), 1719 deletions(-) delete mode 100644 examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl delete mode 100644 examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl delete mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_cnn.jl delete mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl delete mode 100644 examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl delete mode 100644 examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl delete mode 100644 examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl diff --git a/NEWS.md b/NEWS.md index 6baa7f77d23..fc23e74d1db 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,22 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes when updating to v0.6 from v0.5.x + +#### Added + +#### Changed + +#### Deprecated + +#### Removed + +- The neural network-based shock indicators have been migrated to a new repository + [TrixiSmartShockFinder.jl](https://github.com/trixi-framework/TrixiSmartShockFinder.jl). + To continue using the indicators, you will need to use both Trixi.jl and + TrixiSmartShockFinder.jl, as explained in the latter packages' `README.md`. + + ## Changes in the v0.5 lifecycle #### Added @@ -32,6 +48,9 @@ for human readability. #### Removed +- Migrate neural network-based shock indicators to a new repository + [TrixiSmartShockFinder.jl](https://github.com/trixi-framework/TrixiSmartShockFinder.jl). + ## Changes when updating to v0.5 from v0.4.x diff --git a/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl deleted file mode 100644 index 05f392624fd..00000000000 --- a/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl +++ /dev/null @@ -1,109 +0,0 @@ -using Downloads: download -using Flux -using BSON: load -network = joinpath(@__DIR__, "modelnnpp-0.97-0.0001.bson") -download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.97-0.0001.bson", - network) -model1d = load(network, @__MODULE__)[:model1d] - -using OrdinaryDiffEq -using Trixi - -# This elixir was one of the setups used in the following master thesis: -# - Julia Odenthal (2021) -# Shock capturing with artificial neural networks -# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper -# This motivates the particular choice of fluxes, mesh resolution etc. - -############################################################################### -# semidiscretization of the compressible Euler equations - -equations = CompressibleEulerEquations1D(1.4) - -""" - initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) - -A medium blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) -""" -function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) - # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" - # Set up polar coordinates - inicenter = SVector(0.0) - x_norm = x[1] - inicenter[1] - r = abs(x_norm) - # The following code is equivalent to - # phi = atan(0.0, x_norm) - # cos_phi = cos(phi) - # in 1D but faster - cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) - - # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - p = r > 0.5 ? 1.0E-3 : 1.245 - - return prim2cons(SVector(rho, v1, p), equations) -end -initial_condition = initial_condition_blast_wave - -surface_flux = flux_lax_friedrichs -volume_flux = flux_chandrashekar -basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkPerssonPeraire(), - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = false, - alpha_amr = false, - variable = density_pressure, - network = model1d) -volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral) - -coordinates_min = (-2.0,) -coordinates_max = (2.0,) -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 6, - n_cells_max = 10_000) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 12.5) -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 = 100, - 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/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl deleted file mode 100644 index de2f5134a49..00000000000 --- a/examples/tree_1d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl +++ /dev/null @@ -1,109 +0,0 @@ -using Downloads: download -using Flux -using BSON: load -network = joinpath(@__DIR__, "modelnnrh-0.95-0.009.bson") -download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnrh-0.95-0.009.bson", - network) -model1d = load(network, @__MODULE__)[:model1d] - -using OrdinaryDiffEq -using Trixi - -# This elixir was one of the setups used in the following master thesis: -# - Julia Odenthal (2021) -# Shock capturing with artificial neural networks -# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper -# This motivates the particular choice of fluxes, mesh resolution etc. - -############################################################################### -# semidiscretization of the compressible Euler equations - -equations = CompressibleEulerEquations1D(1.4) - -""" - initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) - -A medium blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) -""" -function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) - # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" - # Set up polar coordinates - inicenter = SVector(0.0) - x_norm = x[1] - inicenter[1] - r = abs(x_norm) - # The following code is equivalent to - # phi = atan(0.0, x_norm) - # cos_phi = cos(phi) - # in 1D but faster - cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) - - # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - p = r > 0.5 ? 1.0E-3 : 1.245 - - return prim2cons(SVector(rho, v1, p), equations) -end -initial_condition = initial_condition_blast_wave - -surface_flux = flux_lax_friedrichs -volume_flux = flux_chandrashekar -basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkRayHesthaven(), - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = false, - alpha_amr = false, - variable = density_pressure, - network = model1d) -volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral) - -coordinates_min = (-2.0,) -coordinates_max = (2.0,) -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 6, - n_cells_max = 10_000) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 12.5) -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 = 100, - 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/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_cnn.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_cnn.jl deleted file mode 100644 index 79d7474dc66..00000000000 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_cnn.jl +++ /dev/null @@ -1,107 +0,0 @@ -using Downloads: download -using Flux -using BSON: load -network = joinpath(@__DIR__, "modelcnn-0.964-0.001.bson") -download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelcnn-0.964-0.001.bson", - network) -model2dcnn = load(network, @__MODULE__)[:model2dcnn] - -using OrdinaryDiffEq -using Trixi - -# This elixir was one of the setups used in the following master thesis: -# - Julia Odenthal (2021) -# Shock capturing with artificial neural networks -# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper -# This motivates the particular choice of fluxes, mesh resolution etc. - -############################################################################### -# semidiscretization of the compressible Euler equations - -equations = CompressibleEulerEquations2D(1.4) - -""" - initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - -A medium blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) -""" -function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" - # Set up polar coordinates - inicenter = SVector(0.0, 0.0) - x_norm = x[1] - inicenter[1] - y_norm = x[2] - inicenter[2] - r = sqrt(x_norm^2 + y_norm^2) - phi = atan(y_norm, x_norm) - sin_phi, cos_phi = sincos(phi) - - # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-3 : 1.245 - - return prim2cons(SVector(rho, v1, v2, p), equations) -end -initial_condition = initial_condition_blast_wave - -surface_flux = flux_lax_friedrichs -volume_flux = flux_chandrashekar -basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkCNN(), - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = false, - variable = density_pressure, - network = model2dcnn) -volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral) - -coordinates_min = (-2.0, -2.0) -coordinates_max = (2.0, 2.0) -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 6, - n_cells_max = 10_000) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 12.5) -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 = 100, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim) - -stepsize_callback = StepsizeCallback(cfl = 0.9) - -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/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl deleted file mode 100644 index 27398593efd..00000000000 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl +++ /dev/null @@ -1,107 +0,0 @@ -using Downloads: download -using Flux -using BSON: load -network = joinpath(@__DIR__, "modelnnpp-0.904-0.0005.bson") -download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.904-0.0005.bson", - network) -model2d = load(network, @__MODULE__)[:model2d] - -using OrdinaryDiffEq -using Trixi - -# This elixir was one of the setups used in the following master thesis: -# - Julia Odenthal (2021) -# Shock capturing with artificial neural networks -# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper -# This motivates the particular choice of fluxes, mesh resolution etc. - -############################################################################### -# semidiscretization of the compressible Euler equations - -equations = CompressibleEulerEquations2D(1.4) - -""" - initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - -A medium blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) -""" -function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" - # Set up polar coordinates - inicenter = SVector(0.0, 0.0) - x_norm = x[1] - inicenter[1] - y_norm = x[2] - inicenter[2] - r = sqrt(x_norm^2 + y_norm^2) - phi = atan(y_norm, x_norm) - sin_phi, cos_phi = sincos(phi) - - # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-3 : 1.245 - - return prim2cons(SVector(rho, v1, v2, p), equations) -end -initial_condition = initial_condition_blast_wave - -surface_flux = flux_lax_friedrichs -volume_flux = flux_chandrashekar -basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkPerssonPeraire(), - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = false, - variable = density_pressure, - network = model2d) -volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral) - -coordinates_min = (-2.0, -2.0) -coordinates_max = (2.0, 2.0) -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 6, - n_cells_max = 10_000) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 12.5) -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 = 100, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim) - -stepsize_callback = StepsizeCallback(cfl = 0.9) - -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/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl deleted file mode 100644 index 6c67f948636..00000000000 --- a/examples/tree_2d_dgsem/elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl +++ /dev/null @@ -1,110 +0,0 @@ -using Downloads: download -using Flux -using Random -using BSON: load -network = joinpath(@__DIR__, "modelnnrhs-0.973-0.001.bson") -download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnrhs-0.973-0.001.bson", - network) -model2d = load(network, @__MODULE__)[:model2d] - -using OrdinaryDiffEq -using Trixi - -# This elixir was one of the setups used in the following master thesis: -# - Julia Odenthal (2021) -# Shock capturing with artificial neural networks -# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper -# This motivates the particular choice of fluxes, mesh resolution etc. - -############################################################################### -# semidiscretization of the compressible Euler equations - -equations = CompressibleEulerEquations2D(1.4) - -""" - initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - -A medium blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) -""" -function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" - # Set up polar coordinates - inicenter = SVector(0.0, 0.0) - x_norm = x[1] - inicenter[1] - y_norm = x[2] - inicenter[2] - r = sqrt(x_norm^2 + y_norm^2) - phi = atan(y_norm, x_norm) - sin_phi, cos_phi = sincos(phi) - - # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi - v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi - p = r > 0.5 ? 1.0E-3 : 1.245 - - return prim2cons(SVector(rho, v1, v2, p), equations) -end -initial_condition = initial_condition_blast_wave - -surface_flux = flux_lax_friedrichs -volume_flux = flux_chandrashekar -basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkRayHesthaven(), - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = false, - variable = density_pressure, - network = model2d) -volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral) - -coordinates_min = (-2.0, -2.0) -coordinates_max = (2.0, 2.0) -refinement_patches = () # To allow for specifying them via `trixi_include` -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 6, - refinement_patches = refinement_patches, - n_cells_max = 10_000) - -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 12.5) -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 = 100, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim) - -stepsize_callback = StepsizeCallback(cfl = 0.9) - -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/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl deleted file mode 100644 index d2cab04223e..00000000000 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl +++ /dev/null @@ -1,134 +0,0 @@ -using Downloads: download -using Flux -using BSON: load -network = joinpath(@__DIR__, "modelnnpp-0.904-0.0005.bson") -download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.904-0.0005.bson", - network) -model2d = load(network, @__MODULE__)[:model2d] - -using Random: seed! -seed!(0) - -using OrdinaryDiffEq -using Trixi - -# This elixir was one of the setups used in the following master thesis: -# - Julia Odenthal (2021) -# Shock capturing with artificial neural networks -# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper -# This motivates the particular choice of fluxes, mesh resolution etc. - -############################################################################### -# 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 -https://rsaa.anu.edu.au/research/established-projects/fyris/2-d-kelvin-helmholtz-test. -""" -function initial_condition_kelvin_helmholtz_instability(x, t, - equations::CompressibleEulerEquations2D) - # change discontinuity to tanh - # typical resolution 128^2, 256^2 - # domain size is [-0.5,0.5]^2 - dens0 = 1.0 # outside density - dens1 = 2.0 # inside density - velx0 = -0.5 # outside velocity - velx1 = 0.5 # inside velocity - slope = 50 # used for tanh instead of discontinuous initial condition - # pressure equilibrium - p = 2.5 - # y velocity v2 is only white noise - v2 = 0.01 * (rand(Float64, 1)[1] - 0.5) - # density - rho = dens0 + - (dens1 - dens0) * 0.5 * - (1 + (tanh(slope * (x[2] + 0.25)) - (tanh(slope * (x[2] - 0.25)) + 1))) - # x velocity is also augmented with noise - v1 = velx0 + - (velx1 - velx0) * 0.5 * - (1 + (tanh(slope * (x[2] + 0.25)) - (tanh(slope * (x[2] - 0.25)) + 1))) + - 0.01 * (rand(Float64, 1)[1] - 0.5) - return prim2cons(SVector(rho, v1, v2, p), equations) -end -initial_condition = initial_condition_kelvin_helmholtz_instability - -surface_flux = flux_lax_friedrichs -volume_flux = flux_chandrashekar -polydeg = 3 -basis = LobattoLegendreBasis(polydeg) -indicator_sc = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkPerssonPeraire(), - alpha_max = 0.002, - alpha_min = 0.0001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = false, - variable = density_pressure, - network = model2d) -volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -solver = DGSEM(basis, 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 = 100, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim) - -amr_indicator = IndicatorNeuralNetwork(semi, - indicator_type = NeuralNetworkPerssonPeraire(), - alpha_max = 1.0, - alpha_min = 0.0001, - alpha_smooth = false, - alpha_continuous = true, - alpha_amr = true, - variable = density_pressure, - network = model2d) -amr_controller = ControllerThreeLevel(semi, amr_indicator, - base_level = 4, - med_level = 6, med_threshold = 0.3, # med_level = current level - max_level = 7, max_threshold = 0.5) -amr_callback = AMRCallback(semi, amr_controller, - interval = 1, - adapt_initial_condition = true, - adapt_initial_condition_only_refine = true) - -stepsize_callback = StepsizeCallback(cfl = 1.3) - -callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, - save_solution, - amr_callback, 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/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl deleted file mode 100644 index 39ea947f872..00000000000 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl +++ /dev/null @@ -1,126 +0,0 @@ -using Downloads: download -using Flux -using BSON: load -network = joinpath(@__DIR__, "modelnnpp-0.904-0.0005.bson") -download("https://github.com/trixi-framework/Trixi_IndicatorNeuralNetwork_networks/raw/main/networks/modelnnpp-0.904-0.0005.bson", - network) -model2d = load(network, @__MODULE__)[:model2d] - -using OrdinaryDiffEq -using Trixi - -# This elixir was one of the setups used in the following master thesis: -# - Julia Odenthal (2021) -# Shock capturing with artificial neural networks -# University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper -# This motivates the particular choice of fluxes, mesh resolution etc. - -############################################################################### -# semidiscretization of the compressible Euler equations -gamma = 1.4 -equations = CompressibleEulerEquations2D(gamma) - -""" - initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) - -The Sedov blast wave setup based on Flash -- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 -""" -function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D) - # Set up polar coordinates - inicenter = SVector(0.0, 0.0) - x_norm = x[1] - inicenter[1] - y_norm = x[2] - inicenter[2] - r = sqrt(x_norm^2 + y_norm^2) - - # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 - r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) - # r0 = 0.5 # = more reasonable setup - E = 1.0 - p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) - p0_outer = 1.0e-5 # = true Sedov setup - # p0_outer = 1.0e-3 # = more reasonable setup - - # Calculate primitive variables - rho = 1.0 - v1 = 0.0 - v2 = 0.0 - p = r > r0 ? p0_outer : p0_inner - - return prim2cons(SVector(rho, v1, v2, p), equations) -end -initial_condition = initial_condition_sedov_blast_wave - -surface_flux = flux_lax_friedrichs -volume_flux = flux_chandrashekar -basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkPerssonPeraire(), - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = false, - variable = density_pressure, - network = model2d) -volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; - volume_flux_dg = volume_flux, - volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral) - -coordinates_min = (-2.0, -2.0) -coordinates_max = (2.0, 2.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, 12.5) -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 = 100, - save_initial_solution = true, - save_final_solution = true, - solution_variables = cons2prim) - -amr_indicator = IndicatorNeuralNetwork(semi, - indicator_type = NeuralNetworkPerssonPeraire(), - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = true, - variable = density_pressure, - network = model2d) -amr_controller = ControllerThreeLevel(semi, amr_indicator, - base_level = 4, - max_level = 6, max_threshold = 0.22) -amr_callback = AMRCallback(semi, amr_controller, - interval = 5, - adapt_initial_condition = true, - adapt_initial_condition_only_refine = true) - -stepsize_callback = StepsizeCallback(cfl = 0.9) - -callbacks = CallbackSet(summary_callback, - analysis_callback, alive_callback, - save_solution, - amr_callback, 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/Trixi.jl b/src/Trixi.jl index 97d518d5b78..79810186d4d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -260,9 +260,7 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt, load_adaptive_time_integrator! export ControllerThreeLevel, ControllerThreeLevelCombined, - IndicatorLöhner, IndicatorLoehner, IndicatorMax, - IndicatorNeuralNetwork, NeuralNetworkPerssonPeraire, NeuralNetworkRayHesthaven, - NeuralNetworkCNN + IndicatorLöhner, IndicatorLoehner, IndicatorMax # TODO: TrixiShallowWater: move new limiter export PositivityPreservingLimiterZhangShu, PositivityPreservingLimiterShallowWater @@ -303,10 +301,6 @@ function __init__() end end - @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin - using .Flux: params - end - # FIXME upstream. This is a hacky workaround for # https://github.com/trixi-framework/Trixi.jl/issues/628 # https://github.com/trixi-framework/Trixi.jl/issues/1185 diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index 4b83e9c1a9e..bb9109f2762 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -332,199 +332,4 @@ function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorMax) summary_box(io, "IndicatorMax", setup) end end - -""" - IndicatorNeuralNetwork - -Artificial neural network based indicator used for shock-capturing or AMR. -Depending on the indicator_type, different input values and corresponding trained networks are used. - -`indicator_type = NeuralNetworkPerssonPeraire()` -- Input: The energies in lower modes as well as nnodes(dg). - -`indicator_type = NeuralNetworkRayHesthaven()` -- 1d Input: Cell average of the cell and its neighboring cells as well as the interface values. -- 2d Input: Linear modal values of the cell and its neighboring cells. - -- Ray, Hesthaven (2018) - "An artificial neural network as a troubled-cell indicator" - [doi:10.1016/j.jcp.2018.04.029](https://doi.org/10.1016/j.jcp.2018.04.029) -- Ray, Hesthaven (2019) - "Detecting troubled-cells on two-dimensional unstructured grids using a neural network" - [doi:10.1016/j.jcp.2019.07.043](https://doi.org/10.1016/j.jcp.2019.07.043) - -`indicator_type = CNN (Only in 2d)` -- Based on convolutional neural network. -- 2d Input: Interpolation of the nodal values of the `indicator.variable` to the 4x4 LGL nodes. - -If `alpha_continuous == true` the continuous network output for troubled cells (`alpha > 0.5`) is considered. -If the cells are good (`alpha < 0.5`), `alpha` is set to `0`. -If `alpha_continuous == false`, the blending factor is set to `alpha = 0` for good cells and -`alpha = 1` for troubled cells. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - -""" -struct IndicatorNeuralNetwork{IndicatorType, RealT <: Real, Variable, Chain, Cache} <: - AbstractIndicator - indicator_type::IndicatorType - alpha_max::RealT - alpha_min::RealT - alpha_smooth::Bool - alpha_continuous::Bool - alpha_amr::Bool - variable::Variable - network::Chain - cache::Cache -end - -# this method is used when the indicator is constructed as for shock-capturing volume integrals -function IndicatorNeuralNetwork(equations::AbstractEquations, basis; - indicator_type, - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = false, - variable, - network) - alpha_max, alpha_min = promote(alpha_max, alpha_min) - IndicatorType = typeof(indicator_type) - cache = create_cache(IndicatorNeuralNetwork{IndicatorType}, equations, basis) - IndicatorNeuralNetwork{IndicatorType, typeof(alpha_max), typeof(variable), - typeof(network), typeof(cache)}(indicator_type, alpha_max, - alpha_min, alpha_smooth, - alpha_continuous, alpha_amr, - variable, - network, cache) -end - -# this method is used when the indicator is constructed as for AMR -function IndicatorNeuralNetwork(semi::AbstractSemidiscretization; - indicator_type, - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - alpha_continuous = true, - alpha_amr = true, - variable, - network) - alpha_max, alpha_min = promote(alpha_max, alpha_min) - IndicatorType = typeof(indicator_type) - cache = create_cache(IndicatorNeuralNetwork{IndicatorType}, semi) - IndicatorNeuralNetwork{IndicatorType, typeof(alpha_max), typeof(variable), - typeof(network), typeof(cache)}(indicator_type, alpha_max, - alpha_min, alpha_smooth, - alpha_continuous, alpha_amr, - variable, - network, cache) -end - -function Base.show(io::IO, indicator::IndicatorNeuralNetwork) - @nospecialize indicator # reduce precompilation time - - print(io, "IndicatorNeuralNetwork(") - print(io, indicator.indicator_type) - print(io, ", alpha_max=", indicator.alpha_max) - print(io, ", alpha_min=", indicator.alpha_min) - print(io, ", alpha_smooth=", indicator.alpha_smooth) - print(io, ", alpha_continuous=", indicator.alpha_continuous) - print(io, indicator.variable) - print(io, ")") -end - -function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorNeuralNetwork) - @nospecialize indicator # reduce precompilation time - - if get(io, :compact, false) - show(io, indicator) - else - setup = [ - "indicator type" => indicator.indicator_type, - "max. α" => indicator.alpha_max, - "min. α" => indicator.alpha_min, - "smooth α" => (indicator.alpha_smooth ? "yes" : "no"), - "continuous α" => (indicator.alpha_continuous ? "yes" : "no"), - "indicator variable" => indicator.variable, - ] - summary_box(io, "IndicatorNeuralNetwork", setup) - end -end - -# Convert probability for troubled cell to indicator value for shockcapturing/AMR -@inline function probability_to_indicator(probability_troubled_cell, alpha_continuous, - alpha_amr, - alpha_min, alpha_max) - # Initialize indicator to zero - alpha_element = zero(probability_troubled_cell) - - if alpha_continuous && !alpha_amr - # Set good cells to 0 and troubled cells to continuous value of the network prediction - if probability_troubled_cell > 0.5 - alpha_element = probability_troubled_cell - else - alpha_element = zero(probability_troubled_cell) - end - - # Take care of the case close to pure FV - if alpha_element > 1 - alpha_min - alpha_element = one(alpha_element) - end - - # Scale the probability for a troubled cell (in [0,1]) to the maximum allowed alpha - alpha_element *= alpha_max - elseif !alpha_continuous && !alpha_amr - # Set good cells to 0 and troubled cells to 1 - if probability_troubled_cell > 0.5 - alpha_element = alpha_max - else - alpha_element = zero(alpha_max) - end - elseif alpha_amr - # The entire continuous output of the neural network is used for AMR - alpha_element = probability_troubled_cell - - # Scale the probability for a troubled cell (in [0,1]) to the maximum allowed alpha - alpha_element *= alpha_max - end - - return alpha_element -end - -""" - NeuralNetworkPerssonPeraire - -Indicator type for creating an `IndicatorNeuralNetwork` indicator. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - -See also: [`IndicatorNeuralNetwork`](@ref) -""" -struct NeuralNetworkPerssonPeraire end - -""" - NeuralNetworkRayHesthaven - -Indicator type for creating an `IndicatorNeuralNetwork` indicator. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - -See also: [`IndicatorNeuralNetwork`](@ref) -""" -struct NeuralNetworkRayHesthaven end - -""" - NeuralNetworkCNN - -Indicator type for creating an `IndicatorNeuralNetwork` indicator. - -!!! warning "Experimental implementation" - This is an experimental feature and may change in future releases. - -See also: [`IndicatorNeuralNetwork`](@ref) -""" -struct NeuralNetworkCNN end end # @muladd diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 40bfd1e98c7..4006932352e 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -305,222 +305,4 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 3}, return alpha end - -# this method is used when the indicator is constructed as for shock-capturing volume integrals -# empty cache is default -function create_cache(::Type{<:IndicatorNeuralNetwork}, - equations::AbstractEquations{1}, basis::LobattoLegendreBasis) - return NamedTuple() -end - -# cache for NeuralNetworkPerssonPeraire-type indicator -function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire}}, - equations::AbstractEquations{1}, basis::LobattoLegendreBasis) - alpha = Vector{real(basis)}() - alpha_tmp = similar(alpha) - A = Array{real(basis), ndims(equations)} - - prototype = A(undef, nnodes(basis)) - indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - modal_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - - return (; alpha, alpha_tmp, indicator_threaded, modal_threaded) -end - -# cache for NeuralNetworkRayHesthaven-type indicator -function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkRayHesthaven}}, - equations::AbstractEquations{1}, basis::LobattoLegendreBasis) - alpha = Vector{real(basis)}() - alpha_tmp = similar(alpha) - A = Array{real(basis), ndims(equations)} - - prototype = A(undef, nnodes(basis)) - indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - neighbor_ids = Vector{Int}(undef, 2) - - return (; alpha, alpha_tmp, indicator_threaded, neighbor_ids) -end - -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{<:IndicatorNeuralNetwork}, - mesh, equations::AbstractEquations{1}, dg::DGSEM, cache) - create_cache(typ, equations, dg.basis) -end - -function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire})(u::AbstractArray{ - <:Any, - 3 - }, - mesh, - equations, - dg::DGSEM, - cache; - kwargs...) - @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann - - @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_ann.cache - # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? - # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` - # or just `resize!` whenever we call the relevant methods as we do now? - resize!(alpha, nelements(dg, cache)) - if alpha_smooth - resize!(alpha_tmp, nelements(dg, cache)) - end - - @threaded for element in eachelement(dg, cache) - indicator = indicator_threaded[Threads.threadid()] - modal = modal_threaded[Threads.threadid()] - - # Calculate indicator variables at Gauss-Lobatto nodes - for i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, element) - indicator[i] = indicator_ann.variable(u_local, equations) - end - - # Convert to modal representation - multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, - indicator) - - # Calculate total energies for all modes, without highest, without two highest - total_energy = zero(eltype(modal)) - for i in 1:nnodes(dg) - total_energy += modal[i]^2 - end - total_energy_clip1 = zero(eltype(modal)) - for i in 1:(nnodes(dg) - 1) - total_energy_clip1 += modal[i]^2 - end - total_energy_clip2 = zero(eltype(modal)) - for i in 1:(nnodes(dg) - 2) - total_energy_clip2 += modal[i]^2 - end - - # Calculate energy in highest modes - X1 = (total_energy - total_energy_clip1) / total_energy - X2 = (total_energy_clip1 - total_energy_clip2) / total_energy_clip1 - - # There are two versions of the network: - # The first one only takes the highest energy modes as input, the second one also the number of - # nodes. Automatically use the correct input by checking the number of inputs of the network. - if size(params(network)[1], 2) == 2 - network_input = SVector(X1, X2) - elseif size(params(network)[1], 2) == 3 - network_input = SVector(X1, X2, nnodes(dg)) - end - - # Scale input data - network_input = network_input / - max(maximum(abs, network_input), one(eltype(network_input))) - probability_troubled_cell = network(network_input)[1] - - # Compute indicator value - alpha[element] = probability_to_indicator(probability_troubled_cell, - alpha_continuous, - alpha_amr, alpha_min, alpha_max) - end - - if alpha_smooth - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end - - return alpha -end - -function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkRayHesthaven})(u::AbstractArray{ - <:Any, - 3 - }, - mesh, - equations, - dg::DGSEM, - cache; - kwargs...) - @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann - - @unpack alpha, alpha_tmp, indicator_threaded, neighbor_ids = indicator_ann.cache - # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? - # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` - # or just `resize!` whenever we call the relevant methods as we do now? - resize!(alpha, nelements(dg, cache)) - if alpha_smooth - resize!(alpha_tmp, nelements(dg, cache)) - end - - c2e = zeros(Int, length(mesh.tree)) - for element in eachelement(dg, cache) - c2e[cache.elements.cell_ids[element]] = element - end - - @threaded for element in eachelement(dg, cache) - indicator = indicator_threaded[Threads.threadid()] - cell_id = cache.elements.cell_ids[element] - - for direction in eachdirection(mesh.tree) - if !has_any_neighbor(mesh.tree, cell_id, direction) - neighbor_ids[direction] = element - continue - end - if has_neighbor(mesh.tree, cell_id, direction) - neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id] - if has_children(mesh.tree, neighbor_cell_id) # Cell has small neighbor - if direction == 1 - neighbor_ids[direction] = c2e[mesh.tree.child_ids[2, - neighbor_cell_id]] - else - neighbor_ids[direction] = c2e[mesh.tree.child_ids[1, - neighbor_cell_id]] - end - else # Cell has same refinement level neighbor - neighbor_ids[direction] = c2e[neighbor_cell_id] - end - else # Cell is small and has large neighbor - parent_id = mesh.tree.parent_ids[cell_id] - neighbor_cell_id = mesh.tree.neighbor_ids[direction, parent_id] - neighbor_ids[direction] = c2e[neighbor_cell_id] - end - end - - # Calculate indicator variables at Gauss-Lobatto nodes - for i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, element) - indicator[i] = indicator_ann.variable(u_local, equations) - end - - # Cell average and interface values of the cell - X2 = sum(indicator) / nnodes(dg) - X4 = indicator[1] - X5 = indicator[end] - - # Calculate indicator variables from left neighboring cell at Gauss-Lobatto nodes - for i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, neighbor_ids[1]) - indicator[i] = indicator_ann.variable(u_local, equations) - end - X1 = sum(indicator) / nnodes(dg) - - # Calculate indicator variables from right neighboring cell at Gauss-Lobatto nodes - for i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, neighbor_ids[2]) - indicator[i] = indicator_ann.variable(u_local, equations) - end - X3 = sum(indicator) / nnodes(dg) - network_input = SVector(X1, X2, X3, X4, X5) - - # Scale input data - network_input = network_input / - max(maximum(abs, network_input), one(eltype(network_input))) - probability_troubled_cell = network(network_input)[1] - - # Compute indicator value - alpha[element] = probability_to_indicator(probability_troubled_cell, - alpha_continuous, - alpha_amr, alpha_min, alpha_max) - end - - if alpha_smooth - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end - - return alpha -end end # @muladd diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index da81b2a1d36..8333bb515d3 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -339,355 +339,4 @@ function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any, 4}, return alpha end - -# this method is used when the indicator is constructed as for shock-capturing volume integrals -# empty cache is default -function create_cache(::Type{IndicatorNeuralNetwork}, - equations::AbstractEquations{2}, basis::LobattoLegendreBasis) - return NamedTuple() -end - -# cache for NeuralNetworkPerssonPeraire-type indicator -function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire}}, - equations::AbstractEquations{2}, basis::LobattoLegendreBasis) - alpha = Vector{real(basis)}() - alpha_tmp = similar(alpha) - A = Array{real(basis), ndims(equations)} - - @assert nnodes(basis)>=4 "Indicator only works for nnodes >= 4 (polydeg > 2)" - - prototype = A(undef, nnodes(basis), nnodes(basis)) - indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - modal_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - modal_tmp1_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - - return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded) -end - -# cache for NeuralNetworkRayHesthaven-type indicator -function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkRayHesthaven}}, - equations::AbstractEquations{2}, basis::LobattoLegendreBasis) - alpha = Vector{real(basis)}() - alpha_tmp = similar(alpha) - A = Array{real(basis), ndims(equations)} - - prototype = A(undef, nnodes(basis), nnodes(basis)) - indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - modal_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - modal_tmp1_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - - network_input = Vector{Float64}(undef, 15) - neighbor_ids = Array{Int64}(undef, 8) - neighbor_mean = Array{Float64}(undef, 4, 3) - - return (; alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded, - network_input, neighbor_ids, neighbor_mean) -end - -# cache for NeuralNetworkCNN-type indicator -function create_cache(::Type{IndicatorNeuralNetwork{NeuralNetworkCNN}}, - equations::AbstractEquations{2}, basis::LobattoLegendreBasis) - alpha = Vector{real(basis)}() - alpha_tmp = similar(alpha) - A = Array{real(basis), ndims(equations)} - - prototype = A(undef, nnodes(basis), nnodes(basis)) - indicator_threaded = [similar(prototype) for _ in 1:Threads.nthreads()] - n_cnn = 4 - nodes, _ = gauss_lobatto_nodes_weights(nnodes(basis)) - cnn_nodes, _ = gauss_lobatto_nodes_weights(n_cnn) - vandermonde = polynomial_interpolation_matrix(nodes, cnn_nodes) - network_input = Array{Float32}(undef, n_cnn, n_cnn, 1, 1) - - return (; alpha, alpha_tmp, indicator_threaded, nodes, cnn_nodes, vandermonde, - network_input) -end - -# this method is used when the indicator is constructed as for AMR -function create_cache(typ::Type{<:IndicatorNeuralNetwork}, - mesh, equations::AbstractEquations{2}, dg::DGSEM, cache) - create_cache(typ, equations, dg.basis) -end - -function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkPerssonPeraire})(u, - mesh::TreeMesh{ - 2 - }, - equations, - dg::DGSEM, - cache; - kwargs...) - @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann - - @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded = indicator_ann.cache - # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? - # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` - # or just `resize!` whenever we call the relevant methods as we do now? - resize!(alpha, nelements(dg, cache)) - if alpha_smooth - resize!(alpha_tmp, nelements(dg, cache)) - end - - @threaded for element in eachelement(dg, cache) - indicator = indicator_threaded[Threads.threadid()] - modal = modal_threaded[Threads.threadid()] - modal_tmp1 = modal_tmp1_threaded[Threads.threadid()] - - # Calculate indicator variables at Gauss-Lobatto nodes - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - indicator[i, j] = indicator_ann.variable(u_local, equations) - end - - # Convert to modal representation - multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, - indicator, modal_tmp1) - - # Calculate total energies for all modes, without highest, without two highest - total_energy = zero(eltype(modal)) - for j in 1:nnodes(dg), i in 1:nnodes(dg) - total_energy += modal[i, j]^2 - end - total_energy_clip1 = zero(eltype(modal)) - for j in 1:(nnodes(dg) - 1), i in 1:(nnodes(dg) - 1) - total_energy_clip1 += modal[i, j]^2 - end - total_energy_clip2 = zero(eltype(modal)) - for j in 1:(nnodes(dg) - 2), i in 1:(nnodes(dg) - 2) - total_energy_clip2 += modal[i, j]^2 - end - total_energy_clip3 = zero(eltype(modal)) - for j in 1:(nnodes(dg) - 3), i in 1:(nnodes(dg) - 3) - total_energy_clip3 += modal[i, j]^2 - end - - # Calculate energy in higher modes and polynomial degree for the network input - X1 = (total_energy - total_energy_clip1) / total_energy - X2 = (total_energy_clip1 - total_energy_clip2) / total_energy_clip1 - X3 = (total_energy_clip2 - total_energy_clip3) / total_energy_clip2 - X4 = nnodes(dg) - network_input = SVector(X1, X2, X3, X4) - - # Scale input data - network_input = network_input / - max(maximum(abs, network_input), one(eltype(network_input))) - probability_troubled_cell = network(network_input)[1] - - # Compute indicator value - alpha[element] = probability_to_indicator(probability_troubled_cell, - alpha_continuous, - alpha_amr, alpha_min, alpha_max) - end - - if alpha_smooth - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end - - return alpha -end - -function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkRayHesthaven})(u, - mesh::TreeMesh{ - 2 - }, - equations, - dg::DGSEM, - cache; - kwargs...) - @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann - - @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded, network_input, neighbor_ids, neighbor_mean = indicator_ann.cache #X, network_input - # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? - # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` - # or just `resize!` whenever we call the relevant methods as we do now? - resize!(alpha, nelements(dg, cache)) - if alpha_smooth - resize!(alpha_tmp, nelements(dg, cache)) - end - - c2e = zeros(Int, length(mesh.tree)) - for element in eachelement(dg, cache) - c2e[cache.elements.cell_ids[element]] = element - end - - X = Array{Float64}(undef, 3, nelements(dg, cache)) - - @threaded for element in eachelement(dg, cache) - indicator = indicator_threaded[Threads.threadid()] - modal = modal_threaded[Threads.threadid()] - modal_tmp1 = modal_tmp1_threaded[Threads.threadid()] - - # Calculate indicator variables at Gauss-Lobatto nodes - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - indicator[i, j] = indicator_ann.variable(u_local, equations) - end - - # Convert to modal representation - multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, - indicator, modal_tmp1) - # Save linear modal coefficients for the network input - X[1, element] = modal[1, 1] - X[2, element] = modal[1, 2] - X[3, element] = modal[2, 1] - end - - @threaded for element in eachelement(dg, cache) - cell_id = cache.elements.cell_ids[element] - - network_input[1] = X[1, element] - network_input[2] = X[2, element] - network_input[3] = X[3, element] - - for direction in eachdirection(mesh.tree) - if direction == 1 # -x - dir = 4 - elseif direction == 2 # +x - dir = 1 - elseif direction == 3 # -y - dir = 3 - elseif direction == 4 # +y - dir = 2 - end - - # Of no neighbor exists and current cell is not small - if !has_any_neighbor(mesh.tree, cell_id, direction) - network_input[3 * dir + 1] = X[1, element] - network_input[3 * dir + 2] = X[2, element] - network_input[3 * dir + 3] = X[3, element] - continue - end - - # Get Input data from neighbors - if has_neighbor(mesh.tree, cell_id, direction) - neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id] - if has_children(mesh.tree, neighbor_cell_id) # Cell has small neighbor - # Mean over 4 neighbor cells - neighbor_ids[1] = mesh.tree.child_ids[1, neighbor_cell_id] - neighbor_ids[2] = mesh.tree.child_ids[2, neighbor_cell_id] - neighbor_ids[3] = mesh.tree.child_ids[3, neighbor_cell_id] - neighbor_ids[4] = mesh.tree.child_ids[4, neighbor_cell_id] - - for i in 1:4 - if has_children(mesh.tree, neighbor_ids[i]) - neighbor_ids5 = c2e[mesh.tree.child_ids[1, neighbor_ids[i]]] - neighbor_ids6 = c2e[mesh.tree.child_ids[2, neighbor_ids[i]]] - neighbor_ids7 = c2e[mesh.tree.child_ids[3, neighbor_ids[i]]] - neighbor_ids8 = c2e[mesh.tree.child_ids[4, neighbor_ids[i]]] - - neighbor_mean[i, 1] = (X[1, neighbor_ids5] + - X[1, neighbor_ids6] + - X[1, neighbor_ids7] + - X[1, neighbor_ids8]) / 4 - neighbor_mean[i, 2] = (X[2, neighbor_ids5] + - X[2, neighbor_ids6] + - X[2, neighbor_ids7] + - X[2, neighbor_ids8]) / 4 - neighbor_mean[i, 3] = (X[3, neighbor_ids5] + - X[3, neighbor_ids6] + - X[3, neighbor_ids7] + - X[3, neighbor_ids8]) / 4 - else - neighbor_id = c2e[neighbor_ids[i]] - neighbor_mean[i, 1] = X[1, neighbor_id] - neighbor_mean[i, 2] = X[2, neighbor_id] - neighbor_mean[i, 3] = X[3, neighbor_id] - end - end - network_input[3 * dir + 1] = (neighbor_mean[1, 1] + - neighbor_mean[2, 1] + - neighbor_mean[3, 1] + - neighbor_mean[4, 1]) / 4 - network_input[3 * dir + 2] = (neighbor_mean[1, 2] + - neighbor_mean[2, 2] + - neighbor_mean[3, 2] + - neighbor_mean[4, 2]) / 4 - network_input[3 * dir + 3] = (neighbor_mean[1, 3] + - neighbor_mean[2, 3] + - neighbor_mean[3, 3] + - neighbor_mean[4, 3]) / 4 - - else # Cell has same refinement level neighbor - neighbor_id = c2e[neighbor_cell_id] - network_input[3 * dir + 1] = X[1, neighbor_id] - network_input[3 * dir + 2] = X[2, neighbor_id] - network_input[3 * dir + 3] = X[3, neighbor_id] - end - else # Cell is small and has large neighbor - parent_id = mesh.tree.parent_ids[cell_id] - neighbor_id = c2e[mesh.tree.neighbor_ids[direction, parent_id]] - - network_input[3 * dir + 1] = X[1, neighbor_id] - network_input[3 * dir + 2] = X[2, neighbor_id] - network_input[3 * dir + 3] = X[3, neighbor_id] - end - end - - # Scale input data - network_input = network_input / - max(maximum(abs, network_input), one(eltype(network_input))) - probability_troubled_cell = network(network_input)[1] - - # Compute indicator value - alpha[element] = probability_to_indicator(probability_troubled_cell, - alpha_continuous, - alpha_amr, alpha_min, alpha_max) - end - - if alpha_smooth - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end - - return alpha -end - -function (indicator_ann::IndicatorNeuralNetwork{NeuralNetworkCNN})(u, mesh::TreeMesh{2}, - equations, dg::DGSEM, - cache; kwargs...) - @unpack indicator_type, alpha_max, alpha_min, alpha_smooth, alpha_continuous, alpha_amr, variable, network = indicator_ann - - @unpack alpha, alpha_tmp, indicator_threaded, nodes, cnn_nodes, vandermonde, network_input = indicator_ann.cache - # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? - # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` - # or just `resize!` whenever we call the relevant methods as we do now? - resize!(alpha, nelements(dg, cache)) - if alpha_smooth - resize!(alpha_tmp, nelements(dg, cache)) - end - - @threaded for element in eachelement(dg, cache) - indicator = indicator_threaded[Threads.threadid()] - - # Calculate indicator variables at Gauss-Lobatto nodes - for j in eachnode(dg), i in eachnode(dg) - u_local = get_node_vars(u, equations, dg, i, j, element) - indicator[i, j] = indicator_ann.variable(u_local, equations) - end - - # Interpolate nodal data to 4x4 LGL nodes - for j in 1:4, i in 1:4 - acc = zero(eltype(indicator)) - for jj in eachnode(dg), ii in eachnode(dg) - acc += vandermonde[i, ii] * indicator[ii, jj] * vandermonde[j, jj] - end - network_input[i, j, 1, 1] = acc - end - - # Scale input data - network_input = network_input / - max(maximum(abs, network_input), one(eltype(network_input))) - probability_troubled_cell = network(network_input)[1] - - # Compute indicator value - alpha[element] = probability_to_indicator(probability_troubled_cell, - alpha_continuous, - alpha_amr, alpha_min, alpha_max) - end - - if alpha_smooth - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end - - return alpha -end end # @muladd diff --git a/test/Project.toml b/test/Project.toml index 83b431e269b..fcb96b9e48f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,9 +1,7 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" -BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" @@ -15,10 +13,8 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.7" -BSON = "0.3.3" CairoMakie = "0.6, 0.7, 0.8, 0.9, 0.10" Downloads = "1" -Flux = "0.13.15, 0.14" ForwardDiff = "0.10" LinearAlgebra = "1" MPI = "0.20" diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 62afc5baee3..6cd7998ab02 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -393,26 +393,6 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end - -@trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"), - l2=[0.21814833203212694, 0.2818328665444332, 0.5528379124720818], - linf=[1.5548653877320868, 1.4474018998129738, 2.071919577393772], - maxiters=30) -end - -@trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), - l2=[0.22054468879127423, 0.2828269190680846, 0.5542369885642424], - linf=[ - 1.5623359741479623, - 1.4290121654488288, - 2.1040405133123072, - ], - maxiters=30) -end end end # module diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 93490f8ae09..4a1dc42772d 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -260,89 +260,6 @@ end end end -@trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"), - l2=[ - 0.4758794741390833, - 0.21045415565179362, - 0.21045325630191866, - 0.7022517958549878, - ], - linf=[ - 1.710832148442441, - 0.9711663578827681, - 0.9703787873632452, - 2.9619758810532653, - ], - initial_refinement_level=4, - maxiters=50) -end - -@trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), - l2=[ - 0.472445774440313, - 0.2090782039442978, - 0.20885558673697927, - 0.700569533591275, - ], - linf=[ - 1.7066492792835155, - 0.9856122336679919, - 0.9784316656930644, - 2.9372978989672873, - ], - initial_refinement_level=4, - maxiters=50) -end - -@trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl with mortars" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"), - l2=[ - 0.016486406327766923, - 0.03097329879894433, - 0.03101012918167401, - 0.15157175775429868, - ], - linf=[ - 0.27688647744873407, - 0.5653724536715139, - 0.565695523611447, - 2.513047611639946, - ], - refinement_patches=((type = "box", - coordinates_min = (-0.25, -0.25), - coordinates_max = (0.25, 0.25)), - (type = "box", - coordinates_min = (-0.125, -0.125), - coordinates_max = (0.125, 0.125))), - initial_refinement_level=4, - maxiters=5) -end - -@trixi_testset "elixir_euler_blast_wave_neuralnetwork_cnn.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_blast_wave_neuralnetwork_cnn.jl"), - l2=[ - 0.4795795496408325, - 0.2125148972465021, - 0.21311260934645868, - 0.7033388737692883, - ], - linf=[ - 1.8295385992182336, - 0.9687795218482794, - 0.9616033072376108, - 2.9513245978047133, - ], - initial_refinement_level=4, - maxiters=50, - rtol=1.0e-7) -end - @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_pure_fv.jl"), l2=[ @@ -451,25 +368,6 @@ end end end -@trixi_testset "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl"), - l2=[ - 0.0845430093623868, - 0.09271459184623232, - 0.09271459184623232, - 0.4377291875101709, - ], - linf=[ - 1.3608553480069898, - 1.6822884847136004, - 1.6822884847135997, - 4.2201475428867035, - ], - maxiters=30, - coverage_override=(maxiters = 6,)) -end - @trixi_testset "elixir_euler_positivity.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_positivity.jl"), l2=[ @@ -627,18 +525,6 @@ end end end -@trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl"), - # This stuff is experimental and annoying to test. In the future, we plan - # to move it to another repository. Thus, we save developer time right now - # and do not run these tests anymore. - # l2 = [0.0009823702998067061, 0.004943231496200673, 0.0048604522073091815, 0.00496983530893294], - # linf = [0.00855717053383187, 0.02087422420794427, 0.017121993783086185, 0.02720703869972585], - maxiters=30, - coverage_override=(maxiters = 2,)) -end - @trixi_testset "elixir_euler_colliding_flow.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"), l2=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index 29390161ebe..609100793ba 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -427,14 +427,6 @@ end indicator_max = IndicatorMax("variable", (; cache = nothing)) @test_nowarn show(stdout, indicator_max) - - equations = CompressibleEulerEquations2D(1.4) - basis = LobattoLegendreBasis(3) - indicator_neuralnetwork = IndicatorNeuralNetwork(equations, basis, - indicator_type = NeuralNetworkPerssonPeraire(), - variable = density, - network = nothing) - @test_nowarn show(stdout, indicator_neuralnetwork) end @timed_testset "LBM 2D constructor" begin From 6eb0b5b4c9e803722abe297cff899595ded72da3 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 8 Nov 2023 12:23:07 +0100 Subject: [PATCH 12/17] Make parabolic terms nonexperimental (#1714) * Make parabolic terms non-experimental * Make NSE a separate item * Add MPI to supported features * Mention that parabolic terms are now officially supported in NEWS.md Co-authored-by: Hendrik Ranocha --- NEWS.md | 3 +++ README.md | 4 +++- docs/src/index.md | 4 +++- docs/src/overview.md | 3 +-- src/equations/compressible_navier_stokes.jl | 6 ------ src/equations/compressible_navier_stokes_1d.jl | 3 --- src/equations/compressible_navier_stokes_2d.jl | 3 --- src/equations/compressible_navier_stokes_3d.jl | 3 --- 8 files changed, 10 insertions(+), 19 deletions(-) diff --git a/NEWS.md b/NEWS.md index fc23e74d1db..c037d47db81 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,9 @@ for human readability. #### Changed +- Parabolic diffusion terms are now officially supported and not marked as experimental + anymore. + #### Deprecated #### Removed diff --git a/README.md b/README.md index 1d52089ae3e..c531ab4d1a4 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@

-**Trixi.jl** is a numerical simulation framework for hyperbolic conservation +**Trixi.jl** is a numerical simulation framework for conservation laws written in [Julia](https://julialang.org). A key objective for the framework is to be useful to both scientists and students. Therefore, next to having an extensible design with a fast implementation, Trixi.jl is @@ -46,6 +46,7 @@ installation and postprocessing procedures. Its features include: * Periodic and weakly-enforced boundary conditions * Multiple governing equations: * Compressible Euler equations + * Compressible Navier-Stokes equations * Magnetohydrodynamics (MHD) equations * Multi-component compressible Euler and MHD equations * Linearized Euler and acoustic perturbation equations @@ -56,6 +57,7 @@ installation and postprocessing procedures. Its features include: * Multi-physics simulations * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics) * Shared-memory parallelization via multithreading +* Multi-node parallelization via MPI * Visualization and postprocessing of the results * In-situ and a posteriori visualization with [Plots.jl](https://github.com/JuliaPlots/Plots.jl) * Interactive visualization with [Makie.jl](https://makie.juliaplots.org/) diff --git a/docs/src/index.md b/docs/src/index.md index fd923348928..fbb4b36b224 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -12,7 +12,7 @@ [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3996439.svg)](https://doi.org/10.5281/zenodo.3996439) [**Trixi.jl**](https://github.com/trixi-framework/Trixi.jl) -is a numerical simulation framework for hyperbolic conservation +is a numerical simulation framework for conservation laws written in [Julia](https://julialang.org). A key objective for the framework is to be useful to both scientists and students. Therefore, next to having an extensible design with a fast implementation, Trixi.jl is @@ -40,6 +40,7 @@ installation and postprocessing procedures. Its features include: * Periodic and weakly-enforced boundary conditions * Multiple governing equations: * Compressible Euler equations + * Compressible Navier-Stokes equations * Magnetohydrodynamics (MHD) equations * Multi-component compressible Euler and MHD equations * Linearized Euler and acoustic perturbation equations @@ -50,6 +51,7 @@ installation and postprocessing procedures. Its features include: * Multi-physics simulations * [Self-gravitating gas dynamics](https://github.com/trixi-framework/paper-self-gravitating-gas-dynamics) * Shared-memory parallelization via multithreading +* Multi-node parallelization via MPI * Visualization and postprocessing of the results * In-situ and a posteriori visualization with [Plots.jl](https://github.com/JuliaPlots/Plots.jl) * Interactive visualization with [Makie.jl](https://makie.juliaplots.org/) diff --git a/docs/src/overview.md b/docs/src/overview.md index 9cd11a5df93..9bc523ca297 100644 --- a/docs/src/overview.md +++ b/docs/src/overview.md @@ -60,10 +60,9 @@ different features on different mesh types. | Flux differencing | ✅ | ✅ | ✅ | ✅ | ✅ | [`VolumeIntegralFluxDifferencing`](@ref) | Shock capturing | ✅ | ✅ | ✅ | ✅ | ❌ | [`VolumeIntegralShockCapturingHG`](@ref) | Nonconservative equations | ✅ | ✅ | ✅ | ✅ | ✅ | e.g., GLM MHD or shallow water equations -| Parabolic termsᵇ | ✅ | ✅ | ❌ | ✅ | ✅ | e.g., [`CompressibleNavierStokesDiffusion2D`](@ref) +| Parabolic terms | ✅ | ✅ | ❌ | ✅ | ✅ | e.g., [`CompressibleNavierStokesDiffusion2D`](@ref) ᵃ: quad = quadrilateral, hex = hexahedron -ᵇ: Parabolic terms do not currently support adaptivity. ## Time integration methods diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl index af7897d4586..3059771197c 100644 --- a/src/equations/compressible_navier_stokes.jl +++ b/src/equations/compressible_navier_stokes.jl @@ -6,9 +6,6 @@ Creates a wall-type boundary conditions for the compressible Navier-Stokes equat The fields `boundary_condition_velocity` and `boundary_condition_heat_flux` are intended to be boundary condition types such as the `NoSlip` velocity boundary condition and the `Adiabatic` or `Isothermal` heat boundary condition. - -!!! warning "Experimental feature" - This is an experimental feature and may change in future releases. """ struct BoundaryConditionNavierStokesWall{V, H} boundary_condition_velocity::V @@ -52,9 +49,6 @@ struct Adiabatic{F} end """ -!!! warning "Experimental code" - This code is experimental and may be changed or removed in any future release. - `GradientVariablesPrimitive` and `GradientVariablesEntropy` are gradient variable type parameters for `CompressibleNavierStokesDiffusion1D`. By default, the gradient variables are set to be `GradientVariablesPrimitive`. Specifying `GradientVariablesEntropy` instead uses the entropy variable diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index 74d672ce7ae..f5ae598f389 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -79,9 +79,6 @@ where ```math w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p} ``` - -!!! warning "Experimental code" - This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{1} diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 80857999017..cc0c580bf10 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -79,9 +79,6 @@ where ```math w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} ``` - -!!! warning "Experimental code" - This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2} diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index de2cad99ea8..e9ad3eb6aad 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -79,9 +79,6 @@ where ```math w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p} ``` - -!!! warning "Experimental code" - This code is experimental and may be changed or removed in any future release. """ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{3} From b860e651895276c261fbe5d19f44531c1e498fc3 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Wed, 8 Nov 2023 09:16:58 -0600 Subject: [PATCH 13/17] Deprecate some `DGMultiMesh` constructors (#1709) * remove previously deprecated functions * fix typo in NEWS.md about deprecation vs removal * fix literate tutorial * removing other deprecation * format * Revert "fix typo in NEWS.md about deprecation vs removal" This reverts commit 6b03020a8c881a86550484891a0f53bca018e447. --- docs/literate/src/files/DGMulti_1.jl | 2 +- src/solvers/dgmulti/types.jl | 12 ------------ 2 files changed, 1 insertion(+), 13 deletions(-) diff --git a/docs/literate/src/files/DGMulti_1.jl b/docs/literate/src/files/DGMulti_1.jl index 5ef577e8eeb..6c9a5aea936 100644 --- a/docs/literate/src/files/DGMulti_1.jl +++ b/docs/literate/src/files/DGMulti_1.jl @@ -168,7 +168,7 @@ meshIO = StartUpDG.triangulate_domain(StartUpDG.RectangularDomainWithHole()); # The pre-defined Triangulate geometry in StartUpDG has integer boundary tags. With [`DGMultiMesh`](@ref) # we assign boundary faces based on these integer boundary tags and create a mesh compatible with Trixi.jl. -mesh = DGMultiMesh(meshIO, dg, Dict(:outer_boundary=>1, :inner_boundary=>2)) +mesh = DGMultiMesh(dg, meshIO, Dict(:outer_boundary=>1, :inner_boundary=>2)) #- boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition) boundary_conditions = (; :outer_boundary => boundary_condition_convergence_test, diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl index fe6510856b0..ae1eed7fd52 100644 --- a/src/solvers/dgmulti/types.jl +++ b/src/solvers/dgmulti/types.jl @@ -433,15 +433,3 @@ function LinearAlgebra.mul!(b_in, A_kronecker::SimpleKronecker{3}, x_in) return nothing end end # @muladd - -# TODO: deprecations introduced in Trixi.jl v0.6 -@deprecate DGMultiMesh(dg::DGMulti{NDIMS}; cells_per_dimension, kwargs...) where {NDIMS} DGMultiMesh(dg, - cells_per_dimension; - kwargs...) - -# TODO: deprecations introduced in Trixi.jl v0.5 -@deprecate DGMultiMesh(vertex_coordinates, EToV, dg::DGMulti{NDIMS}; - kwargs...) where {NDIMS} DGMultiMesh(dg, vertex_coordinates, EToV; - kwargs...) -@deprecate DGMultiMesh(triangulateIO, dg::DGMulti{2, Tri}, boundary_dict::Dict{Symbol, Int}; - kwargs...) DGMultiMesh(dg, triangulateIO, boundary_dict; kwargs...) From 07ef07c4c6a7c930b9e6b29e09fe1d1e79d37ffc Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Wed, 8 Nov 2023 09:17:35 -0600 Subject: [PATCH 14/17] add gradient variable type parameter to `AbstractEquationsParabolic` (#1409) * add gradient variable type parameter * fix parabolic literate test * remove trailing comment * remove unnecessary abstract type * move gradient variable structs * formatting * fix dropped changes * try to fix doc tests * fixing navier stokes 1D * formatting * remove duplicate GradientVariablesPrimitive/Entropy definition * update news --- NEWS.md | 1 + docs/literate/src/files/adding_new_parabolic_terms.jl | 11 +++++++++-- src/Trixi.jl | 2 +- src/equations/compressible_navier_stokes_1d.jl | 2 +- src/equations/compressible_navier_stokes_2d.jl | 2 +- src/equations/compressible_navier_stokes_3d.jl | 2 +- src/equations/equations.jl | 2 +- src/equations/equations_parabolic.jl | 11 ++++++++--- 8 files changed, 23 insertions(+), 10 deletions(-) diff --git a/NEWS.md b/NEWS.md index c037d47db81..5d258fa65bb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -35,6 +35,7 @@ for human readability. - Implementation of the quasi-1D shallow water equations - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh` - AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh` +- Added `GradientVariables` type parameter to `AbstractEquationsParabolic` #### Changed diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl index f5c2b815f33..209ef62c988 100644 --- a/docs/literate/src/files/adding_new_parabolic_terms.jl +++ b/docs/literate/src/files/adding_new_parabolic_terms.jl @@ -18,8 +18,15 @@ equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity); # `ConstantAnisotropicDiffusion2D` has a field for `equations_hyperbolic`. It is useful to have # information about the hyperbolic system available to the parabolic part so that we can reuse # functions defined for hyperbolic equations (such as `varnames`). - -struct ConstantAnisotropicDiffusion2D{E, T} <: Trixi.AbstractEquationsParabolic{2, 1} +# +# The abstract type `Trixi.AbstractEquationsParabolic` has three parameters: `NDIMS` (the spatial dimension, +# e.g., 1D, 2D, or 3D), `NVARS` (the number of variables), and `GradientVariable`, which we set as +# `GradientVariablesConservative`. This indicates that the gradient should be taken with respect to the +# conservative variables (e.g., the same variables used in `equations_hyperbolic`). Users can also take +# the gradient with respect to a different set of variables; see, for example, the implementation of +# [`CompressibleNavierStokesDiffusion2D`](@ref), which can utilize either "primitive" or "entropy" variables. + +struct ConstantAnisotropicDiffusion2D{E, T} <: Trixi.AbstractEquationsParabolic{2, 1, GradientVariablesConservative} diffusivity::T equations_hyperbolic::E end diff --git a/src/Trixi.jl b/src/Trixi.jl index 79810186d4d..2d58a6b72a4 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -157,7 +157,7 @@ export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D, CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D -export GradientVariablesPrimitive, GradientVariablesEntropy +export GradientVariablesConservative, GradientVariablesPrimitive, GradientVariablesEntropy export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index f5ae598f389..73436c99b7c 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -83,7 +83,7 @@ w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p} struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{1} } <: - AbstractCompressibleNavierStokesDiffusion{1, 3} + AbstractCompressibleNavierStokesDiffusion{1, 3, GradientVariables} # TODO: parabolic # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index cc0c580bf10..ad0db001872 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -83,7 +83,7 @@ w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{2} } <: - AbstractCompressibleNavierStokesDiffusion{2, 4} + AbstractCompressibleNavierStokesDiffusion{2, 4, GradientVariables} # TODO: parabolic # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index e9ad3eb6aad..c6a55983b53 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -83,7 +83,7 @@ w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p} struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, E <: AbstractCompressibleEulerEquations{3} } <: - AbstractCompressibleNavierStokesDiffusion{3, 5} + AbstractCompressibleNavierStokesDiffusion{3, 5, GradientVariables} # TODO: parabolic # 1) For now save gamma and inv(gamma-1) again, but could potentially reuse them from the Euler equations # 2) Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 0e77b92e045..78b1b829b06 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -479,6 +479,6 @@ abstract type AbstractLinearizedEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("linearized_euler_2d.jl") -abstract type AbstractEquationsParabolic{NDIMS, NVARS} <: +abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <: AbstractEquations{NDIMS, NVARS} end end # @muladd diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 47a76174cb1..a063e9f2758 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -2,16 +2,21 @@ # specialize this function to compute gradients e.g., of primitive variables instead of conservative gradient_variable_transformation(::AbstractEquationsParabolic) = cons2cons +# By default, the gradients are taken with respect to the conservative variables. +# this is reflected by the type parameter `GradientVariablesConservative` in the abstract +# type `AbstractEquationsParabolic{NDIMS, NVARS, GradientVariablesConservative}`. +struct GradientVariablesConservative end + # Linear scalar diffusion for use in linear scalar advection-diffusion problems abstract type AbstractLaplaceDiffusion{NDIMS, NVARS} <: - AbstractEquationsParabolic{NDIMS, NVARS} end + AbstractEquationsParabolic{NDIMS, NVARS, GradientVariablesConservative} end include("laplace_diffusion_1d.jl") include("laplace_diffusion_2d.jl") include("laplace_diffusion_3d.jl") # Compressible Navier-Stokes equations -abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: - AbstractEquationsParabolic{NDIMS, NVARS} end +abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS, GradientVariables} <: + AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} end include("compressible_navier_stokes.jl") include("compressible_navier_stokes_1d.jl") include("compressible_navier_stokes_2d.jl") From b75cab9859e592f1d1d1514e7e8628e7621b3e3c Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Fri, 10 Nov 2023 07:54:56 +0100 Subject: [PATCH 15/17] Add new non-conservative flux for standard and two-layer shallow water equations (#1508) * Add flux_nonconservative_ersing_etal * Change formulation of flux_nonconservative_ersing * remove old fluxes (fjordholm, wintermeyer) * Update tests for new flux_ersing_etal * Add flux_nonconservative_ersing for standard SWE * Edit comments, change flux in elixir * fix well-balanced test * Update src/Trixi.jl Change order of exporting statements Co-authored-by: Andrew Winters * fix indentation * fix indentation II * Update src/equations/shallow_water_2d.jl Co-authored-by: Andrew Winters * fix energy total function * Apply formatter * add docstring references * remove flux_nonconservative_fjordholm_etal * apply formatter * Format examples --------- Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha --- ...lixir_shallowwater_twolayer_convergence.jl | 4 +- .../elixir_shallowwater_twolayer_dam_break.jl | 4 +- ...xir_shallowwater_twolayer_well_balanced.jl | 4 +- ...lixir_shallowwater_twolayer_convergence.jl | 4 +- ...xir_shallowwater_twolayer_well_balanced.jl | 4 +- ...lixir_shallowwater_twolayer_convergence.jl | 4 +- .../elixir_shallowwater_twolayer_dam_break.jl | 4 +- ...xir_shallowwater_twolayer_well_balanced.jl | 4 +- src/Trixi.jl | 3 +- src/equations/shallow_water_1d.jl | 38 ++ src/equations/shallow_water_2d.jl | 68 +++ src/equations/shallow_water_two_layer_1d.jl | 293 +++++-------- src/equations/shallow_water_two_layer_2d.jl | 391 +++++------------- test/test_tree_1d_shallowwater.jl | 50 +++ test/test_tree_1d_shallowwater_twolayer.jl | 138 +++---- test/test_tree_2d_shallowwater.jl | 58 +++ test/test_tree_2d_shallowwater_twolayer.jl | 162 +++----- test/test_unstructured_2d.jl | 94 ++++- 18 files changed, 630 insertions(+), 697 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl index fc76b4f034b..e6a01849852 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -13,9 +13,9 @@ initial_condition = initial_condition_convergence_test ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 3, - surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl index b2e6a81401b..03b93754d0f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -36,9 +36,9 @@ initial_condition = initial_condition_dam_break ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 3, - surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl index 7236f1697d0..098e3aaf601 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -35,9 +35,9 @@ initial_condition = initial_condition_fjordholm_well_balanced ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 3, - surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal), + surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index f7c8ab3a249..790916e4467 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -13,9 +13,9 @@ initial_condition = initial_condition_convergence_test ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 3, - surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl index 1495e6d8568..264c26390fe 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -31,8 +31,8 @@ initial_condition = initial_condition_well_balanced ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 3, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index 2cab68b1cb5..0b86095663a 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -15,8 +15,8 @@ initial_condition = initial_condition_convergence_test ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 6, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl index 9d70e9287cf..4ad5f7e3201 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -41,8 +41,8 @@ boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_b ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 6, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl index 35b027c3a81..6a727df2502 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -33,8 +33,8 @@ initial_condition = initial_condition_well_balanced ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg = 6, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/src/Trixi.jl b/src/Trixi.jl index 2d58a6b72a4..b8110cf5bdd 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -164,8 +164,9 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, flux_nonconservative_powell_local_symmetric, flux_kennedy_gruber, flux_shima_etal, flux_ec, - flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal, + flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, + flux_es_ersing_etal, flux_nonconservative_ersing_etal, flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, # TODO: TrixiShallowWater: move anything with "chen_noelle" to new file diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 32782d5478c..25ce0fa79fe 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -380,6 +380,44 @@ Further details on the hydrostatic reconstruction and its motivation can be foun z) end +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations1D) + +!!! warning "Experimental code" + This numerical flux is experimental and may change in any future release. + +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref). + +This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy +conservation and well-balancedness in both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations1D) + # Pull the necessary left and right state information + h_ll = waterheight(u_ll, equations) + b_rr = u_rr[3] + b_ll = u_ll[3] + + # Calculate jump + b_jump = b_rr - b_ll + + z = zero(eltype(u_ll)) + + # Bottom gradient nonconservative term: (0, g h b_x, 0) + f = SVector(z, equations.gravity * h_ll * b_jump, z) + + return f +end + """ flux_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterEquations1D) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index a81fddeed49..e75c92a27d0 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -702,6 +702,74 @@ end return SVector(f1, f2, f3, f4) end +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquations2D) + +!!! warning "Experimental code" + This numerical flux is experimental and may change in any future release. + +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref). + +On curvilinear meshes, this nonconservative flux depends on both the +contravariant vector (normal direction) at the current node and the averaged +one. This is different from numerical fluxes used to discretize conservative +terms. + +This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy +conservation and well-balancedness in both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations2D) + # Pull the necessary left and right state information + h_ll = waterheight(u_ll, equations) + b_rr = u_rr[4] + b_ll = u_ll[4] + + # Calculate jump + b_jump = b_rr - b_ll + + z = zero(eltype(u_ll)) + # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) + if orientation == 1 + f = SVector(z, equations.gravity * h_ll * b_jump, z, z) + else # orientation == 2 + f = SVector(z, z, equations.gravity * h_ll * b_jump, z) + end + return f +end + +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquations2D) + # Pull the necessary left and right state information + h_ll = waterheight(u_ll, equations) + b_rr = u_rr[4] + b_ll = u_ll[4] + + # Calculate jump + b_jump = b_rr - b_ll + # Note this routine only uses the `normal_direction_average` and the average of the + # bottom topography to get a quadratic split form DG gradient on curved elements + return SVector(zero(eltype(u_ll)), + normal_direction_average[1] * equations.gravity * h_ll * b_jump, + normal_direction_average[2] * equations.gravity * h_ll * b_jump, + zero(eltype(u_ll))) +end + """ flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterEquations2D) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 4b64481cca3..42ff393593e 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -87,15 +87,15 @@ end have_nonconservative_terms(::ShallowWaterTwoLayerEquations1D) = True() function varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations1D) - ("h_upper", "h_v1_upper", - "h_lower", "h_v1_lower", "b") + ("h_upper", "h_v_upper", + "h_lower", "h_v_lower", "b") end # Note, we use the total water height, H_lower = h_upper + h_lower + b, and first layer total height # H_upper = h_upper + b as the first primitive variable for easier visualization and setting initial # conditions function varnames(::typeof(cons2prim), ::ShallowWaterTwoLayerEquations1D) - ("H_upper", "v1_upper", - "H_lower", "v1_lower", "b") + ("H_upper", "v_upper", + "H_lower", "v_lower", "b") end # Set initial conditions at physical location `x` for time `t` @@ -113,11 +113,11 @@ function initial_condition_convergence_test(x, t, H_lower = 2.0 + 0.1 * sin(ω * x[1] + t) H_upper = 4.0 + 0.1 * cos(ω * x[1] + t) - v1_lower = 1.0 - v1_upper = 0.9 + v_lower = 1.0 + v_upper = 0.9 b = 1.0 + 0.1 * cos(2.0 * ω * x[1]) - return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) + return prim2cons(SVector(H_upper, v_upper, H_lower, v_lower, b), equations) end """ @@ -196,165 +196,77 @@ end # Note, the bottom topography has no flux @inline function flux(u, orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u + h_upper, h_v_upper, h_lower, h_v_lower, _ = u # Calculate velocities - v1_upper, v1_lower = velocity(u, equations) + v_upper, v_lower = velocity(u, equations) # Calculate pressure - p1 = 0.5 * equations.gravity * h_upper^2 - p2 = 0.5 * equations.gravity * h_lower^2 + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 - f1 = h_v1_upper - f2 = h_v1_upper * v1_upper + p1 - f3 = h_v2_lower - f4 = h_v2_lower * v1_lower + p2 + f1 = h_v_upper + f2 = h_v_upper * v_upper + p_upper + f3 = h_v_lower + f4 = h_v_lower * v_lower + p_lower return SVector(f1, f2, f3, f4, zero(eltype(u))) end """ - flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. -Non-symmetric two-point volume flux discretizing the nonconservative (source) term -that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. This is a slightly modified version -to account for the additional source term compared to the standard SWE described in the paper. +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations1D`](@ref) and an +additional term that couples the momentum of both layers. -Further details are available in the paper: -- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) - An entropy stable nodal discontinuous Galerkin method for the two dimensional - shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry - [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy +conservation and well-balancedness in both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ -@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) # Pull the necessary left and right state information h_upper_ll, h_lower_ll = waterheight(u_ll, equations) h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[5] + b_ll = u_ll[5] - z = zero(eltype(u_ll)) - - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, - # 0, g*h_lower*(b+r*h_upper)_x, 0) - f = SVector(z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), - z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), - z) - return f -end - -""" - flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - -!!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. - -Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains -the gradients of the bottom topography and an additional term that couples the momentum of both -layers [`ShallowWaterTwoLayerEquations2D`](@ref). - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Pull the necessary left and right state information - h_upper_ll, _, h_lower_ll, _, b_ll = u_ll - h_upper_rr, _, h_lower_rr, _, b_rr = u_rr - - # Create average and jump values - h_upper_average = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_average = 0.5 * (h_lower_ll + h_lower_rr) - h_upper_jump = h_upper_rr - h_upper_ll - h_lower_jump = h_lower_rr - h_lower_ll - b_jump = b_rr - b_ll - - # Assign variables for constants for better readability - g = equations.gravity + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) z = zero(eltype(u_ll)) # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, # 0, g*h_lower*(b+r*h_upper)_x, 0) f = SVector(z, - g * h_upper_ll * (b_ll + h_lower_ll) + - g * h_upper_average * (b_jump + h_lower_jump), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, - g * h_lower_ll * (b_ll + equations.r * h_upper_ll) + - g * h_lower_average * (b_jump + - equations.r * h_upper_jump), + equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), z) return f end -""" - flux_fjordholm_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations1D) - -Total energy conservative (mathematical entropy for shallow water equations). When the bottom -topography is nonzero this should only be used as a surface flux otherwise the scheme will not be -well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). - -Details are available in Eq. (4.1) in the paper: -- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) - Well-balanced and energy stable schemes for the shallow water equations with discontinuous - topography [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) -and the application to two layers is shown in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Unpack left and right state - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - v1_ll, v2_ll = velocity(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - v1_rr, v2_rr = velocity(u_rr, equations) - - # Average each factor of products in flux - h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p1_avg = 0.25 * equations.gravity * (h_upper_ll^2 + h_upper_rr^2) - p2_avg = 0.25 * equations.gravity * (h_lower_ll^2 + h_lower_rr^2) - - # Calculate fluxes - f1 = h_upper_avg * v1_avg - f2 = f1 * v1_avg + p1_avg - f3 = h_lower_avg * v2_avg - f4 = f3 * v2_avg + p2_avg - - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) -end - """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. -When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. -The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow +When the bottom topography is nonzero this scheme will be well-balanced when used with the +nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the +two-layer shallow water equations the flux that is described in the paper for the normal shallow water equations is used within each layer. Further details are available in Theorem 1 of the paper: @@ -367,51 +279,48 @@ Further details are available in Theorem 1 of the paper: orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) # Unpack left and right state - h_upper_ll, h_v1_upper_ll, h_lower_ll, h_v2_lower_ll, _ = u_ll - h_upper_rr, h_v1_upper_rr, h_lower_rr, h_v2_lower_rr, _ = u_rr + h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll + h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr # Get the velocities on either side - v1_ll, v2_ll = velocity(u_ll, equations) - v1_rr, v2_rr = velocity(u_rr, equations) + v_upper_ll, v_lower_ll = velocity(u_ll, equations) + v_upper_rr, v_lower_rr = velocity(u_rr, equations) # Average each factor of products in flux - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p1_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p2_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + v_upper_avg = 0.5 * (v_upper_ll + v_upper_rr) + v_lower_avg = 0.5 * (v_lower_ll + v_lower_rr) + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr # Calculate fluxes - f1 = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) - f2 = f1 * v1_avg + p1_avg - f3 = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) - f4 = f3 * v2_avg + p2_avg + f1 = 0.5 * (h_v_upper_ll + h_v_upper_rr) + f2 = f1 * v_upper_avg + p_upper_avg + f3 = 0.5 * (h_v_lower_ll + h_v_lower_rr) + f4 = f3 * v_lower_avg + p_lower_avg return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end """ - flux_es_fjordholm_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations1D) - -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy -conservative [`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump -of entropy variables. - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. + flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations1D) +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ -@inline function flux_es_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) +@inline function flux_es_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) # Compute entropy conservative flux but without the bottom topography - f_ec = flux_fjordholm_etal(u_ll, u_rr, - orientation, - equations) + f_ec = flux_wintermeyer_etal(u_ll, u_rr, + orientation, + equations) # Get maximum signal velocity λ = max_abs_speed_naive(u_ll, u_rr, orientation, equations) @@ -474,12 +383,12 @@ end orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) # Unpack left and right state - h_upper_ll, h_v1_upper_ll, h_lower_ll, h_v2_lower_ll, _ = u_ll - h_upper_rr, h_v1_upper_rr, h_lower_rr, h_v2_lower_rr, _ = u_rr + h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll + h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr # Get the averaged velocity - v_m_ll = (h_v1_upper_ll + h_v2_lower_ll) / (h_upper_ll + h_lower_ll) - v_m_rr = (h_v1_upper_rr + h_v2_lower_rr) / (h_upper_rr + h_lower_rr) + v_m_ll = (h_v_upper_ll + h_v_lower_ll) / (h_upper_ll + h_lower_ll) + v_m_rr = (h_v_upper_rr + h_v_lower_rr) / (h_upper_rr + h_lower_rr) # Calculate the wave celerity on the left and right h_upper_ll, h_lower_ll = waterheight(u_ll, equations) @@ -503,10 +412,10 @@ end # Absolute speed of the barotropic mode @inline function max_abs_speeds(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u + h_upper, h_v_upper, h_lower, h_v_lower, _ = u # Calculate averaged velocity of both layers - v_m = (h_v1_upper + h_v2_lower) / (h_upper + h_lower) + v_m = (h_v_upper + h_v_lower) / (h_upper + h_lower) c = sqrt(equations.gravity * (h_upper + h_lower)) return (abs(v_m) + c) @@ -514,11 +423,11 @@ end # Helper function to extract the velocity vector from the conservative variables @inline function velocity(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u + h_upper, h_v_upper, h_lower, h_v_lower, _ = u - v1_upper = h_v1_upper / h_upper - v1_lower = h_v2_lower / h_lower - return SVector(v1_upper, v1_lower) + v_upper = h_v_upper / h_upper + v_lower = h_v_lower / h_lower + return SVector(v_upper, v_lower) end # Convert conservative variables to primitive @@ -527,8 +436,8 @@ end H_lower = h_lower + b H_upper = h_lower + h_upper + b - v1_upper, v1_lower = velocity(u, equations) - return SVector(H_upper, v1_upper, H_lower, v1_lower, b) + v_upper, v_lower = velocity(u, equations) + return SVector(H_upper, v_upper, H_lower, v_lower, b) end # Convert conservative variables to entropy variables @@ -536,26 +445,26 @@ end # bottom topography values for convenience @inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations1D) h_upper, _, h_lower, _, b = u - v1_upper, v1_lower = velocity(u, equations) - - w1 = equations.rho_upper * - (equations.gravity * (h_upper + h_lower + b) - 0.5 * v1_upper^2) - w2 = equations.rho_upper * v1_upper - w3 = equations.rho_lower * - (equations.gravity * (equations.r * h_upper + h_lower + b) - 0.5 * v1_lower^2) - w4 = equations.rho_lower * v1_lower + v_upper, v_lower = velocity(u, equations) + + w1 = (equations.rho_upper * + (equations.gravity * (h_upper + h_lower + b) - 0.5 * v_upper^2)) + w2 = equations.rho_upper * v_upper + w3 = (equations.rho_lower * + (equations.gravity * (equations.r * h_upper + h_lower + b) - 0.5 * v_lower^2)) + w4 = equations.rho_lower * v_lower return SVector(w1, w2, w3, w4, b) end # Convert primitive to conservative variables @inline function prim2cons(prim, equations::ShallowWaterTwoLayerEquations1D) - H_upper, v1_upper, H_lower, v1_lower, b = prim + H_upper, v_upper, H_lower, v_lower, b = prim h_lower = H_lower - b h_upper = H_upper - h_lower - b - h_v1_upper = h_upper * v1_upper - h_v2_lower = h_lower * v1_lower - return SVector(h_upper, h_v1_upper, h_lower, h_v2_lower, b) + h_v_upper = h_upper * v_upper + h_v_lower = h_lower * v_lower + return SVector(h_upper, h_v_upper, h_lower, h_v_lower, b) end @inline function waterheight(u, equations::ShallowWaterTwoLayerEquations1D) @@ -569,23 +478,23 @@ end # Calculate total energy for a conservative state `cons` @inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, b = cons + h_upper, h_v_upper, h_lower, h_v_lower, b = cons # Set new variables for better readability g = equations.gravity rho_upper = equations.rho_upper rho_lower = equations.rho_lower - e = (0.5 * rho_upper * (h_v1_upper^2 / h_upper + g * h_upper^2) + - 0.5 * rho_lower * (h_v2_lower^2 / h_lower + g * h_lower^2) + + e = (0.5 * rho_upper * (h_v_upper^2 / h_upper + g * h_upper^2) + + 0.5 * rho_lower * (h_v_lower^2 / h_lower + g * h_lower^2) + g * rho_lower * h_lower * b + g * rho_upper * h_upper * (h_lower + b)) return e end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u - return 0.5 * equations.rho_upper * h_v1_upper^2 / h_upper + - 0.5 * equations.rho_lower * h_v2_lower^2 / h_lower + h_upper, h_v_upper, h_lower, h_v_lower, _ = u + return (0.5 * equations.rho_upper * h_v_upper^2 / h_upper + + 0.5 * equations.rho_lower * h_v_lower^2 / h_lower) end # Calculate potential energy for a conservative state `cons` diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 87249e91948..a31d881f2ef 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -271,24 +271,24 @@ end v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) # Calculate pressure - p1 = 0.5 * equations.gravity * h_upper^2 - p2 = 0.5 * equations.gravity * h_lower^2 + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 # Calculate fluxes depending on orientation if orientation == 1 f1 = h_v1_upper - f2 = h_v1_upper * v1_upper + p1 + f2 = h_v1_upper * v1_upper + p_upper f3 = h_v1_upper * v2_upper f4 = h_v1_lower - f5 = h_v1_lower * v1_lower + p2 + f5 = h_v1_lower * v1_lower + p_lower f6 = h_v1_lower * v2_lower else f1 = h_v2_upper f2 = h_v2_upper * v1_upper - f3 = h_v2_upper * v2_upper + p1 + f3 = h_v2_upper * v2_upper + p_upper f4 = h_v2_lower f5 = h_v2_lower * v1_lower - f6 = h_v2_lower * v2_lower + p2 + f6 = h_v2_lower * v2_lower + p_lower end return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) end @@ -305,44 +305,57 @@ end h_v_upper_normal = h_upper * v_normal_upper h_v_lower_normal = h_lower * v_normal_lower - p1 = 0.5 * equations.gravity * h_upper^2 - p2 = 0.5 * equations.gravity * h_lower^2 + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 f1 = h_v_upper_normal - f2 = h_v_upper_normal * v1_upper + p1 * normal_direction[1] - f3 = h_v_upper_normal * v2_upper + p1 * normal_direction[2] + f2 = h_v_upper_normal * v1_upper + p_upper * normal_direction[1] + f3 = h_v_upper_normal * v2_upper + p_upper * normal_direction[2] f4 = h_v_lower_normal - f5 = h_v_lower_normal * v1_lower + p2 * normal_direction[1] - f6 = h_v_lower_normal * v2_lower + p2 * normal_direction[2] + f5 = h_v_lower_normal * v1_lower + p_lower * normal_direction[1] + f6 = h_v_lower_normal * v2_lower + p_lower * normal_direction[2] return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) end """ - flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) !!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. + This numerical flux is experimental and may change in any future release. -Non-symmetric two-point volume flux discretizing the nonconservative (source) term +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. This is a slightly modified version -to account for the additional source term compared to the standard SWE described in the paper. +additional term that couples the momentum of both layers. -Further details are available in the paper: -- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) - An entropy stable nodal discontinuous Galerkin method for the two dimensional - shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry - [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy +conservation and well-balancedness in both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ -@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) # Pull the necessary left and right state information h_upper_ll, h_lower_ll = waterheight(u_ll, equations) h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[7] + b_ll = u_ll[7] + + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) z = zero(eltype(u_ll)) @@ -351,268 +364,64 @@ Further details are available in the paper: # g*h_lower*(b + r*h_upper)_y, 0) if orientation == 1 f = SVector(z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), + equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), z, z) else # orientation == 2 f = SVector(z, z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), + equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), z) end return f end -@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) # Pull the necessary left and right state information h_upper_ll, h_lower_ll = waterheight(u_ll, equations) h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[7] + b_ll = u_ll[7] + + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) # Note this routine only uses the `normal_direction_average` and the average of the # bottom topography to get a quadratic split form DG gradient on curved elements return SVector(zero(eltype(u_ll)), normal_direction_average[1] * equations.gravity * h_upper_ll * - (b_rr + h_lower_rr), + (b_jump + h_lower_jump), normal_direction_average[2] * equations.gravity * h_upper_ll * - (b_rr + h_lower_rr), + (b_jump + h_lower_jump), zero(eltype(u_ll)), normal_direction_average[1] * equations.gravity * h_lower_ll * - (b_rr + - equations.r * h_upper_rr), + (b_jump + equations.r * h_upper_jump), normal_direction_average[2] * equations.gravity * h_lower_ll * - (b_rr + - equations.r * h_upper_rr), + (b_jump + equations.r * h_upper_jump), zero(eltype(u_ll))) end -""" - flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - -!!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. - -Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains -the gradients of the bottom topography and an additional term that couples the momentum of both -layers [`ShallowWaterTwoLayerEquations2D`](@ref). - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - # Pull the necessary left and right state information - h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, b_ll = u_ll - h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, b_rr = u_rr - - # Create average and jump values - h_upper_average = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_average = 0.5 * (h_lower_ll + h_lower_rr) - h_upper_jump = h_upper_rr - h_upper_ll - h_lower_jump = h_lower_rr - h_lower_ll - b_jump = b_rr - b_ll - - # Assign variables for constants for better readability - g = equations.gravity - - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, g*h_upper*(b+h_lower)_y, 0, - # g*h_lower*(b+r*h_upper)_x, g*h_lower*(b+r*h_upper)_x, 0) - - # Includes two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid - # cross-averaging across a discontinuous bottom topography - # (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry - z = zero(eltype(u_ll)) - if orientation == 1 - f = SVector(z, - g * h_upper_ll * (b_ll + h_lower_ll) + - g * h_upper_average * (b_jump + h_lower_jump), - z, z, - g * h_lower_ll * (b_ll + equations.r * h_upper_ll) + - g * h_lower_average * (b_jump + - equations.r * h_upper_jump), - z, z) - else # orientation == 2 - f = SVector(z, z, - g * h_upper_ll * (b_ll + h_lower_ll) + - g * h_upper_average * (b_jump + h_lower_jump), - z, z, - g * h_lower_ll * (b_ll + equations.r * h_upper_ll) + - g * h_lower_average * (b_jump + - equations.r * h_upper_jump), - z) - end - - return f -end - -@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - # Pull the necessary left and right state information - h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, b_ll = u_ll - h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, b_rr = u_rr - - # Create average and jump values - h_upper_average = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_average = 0.5 * (h_lower_ll + h_lower_rr) - h_upper_jump = h_upper_rr - h_upper_ll - h_lower_jump = h_lower_rr - h_lower_ll - b_jump = b_rr - b_ll - - # Comes in two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `normal_direction_average` - # but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography - f2 = normal_direction_average[1] * equations.gravity * h_upper_ll * - (b_ll + h_lower_ll) - f3 = normal_direction_average[2] * equations.gravity * h_upper_ll * - (b_ll + h_lower_ll) - f5 = normal_direction_average[1] * equations.gravity * h_lower_ll * - (b_ll + equations.r * h_upper_ll) - f6 = normal_direction_average[2] * equations.gravity * h_lower_ll * - (b_ll + equations.r * h_upper_ll) - # (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump` - # to handle discontinuous bathymetry - f2 += normal_direction_ll[1] * equations.gravity * h_upper_average * - (b_jump + h_lower_jump) - f3 += normal_direction_ll[2] * equations.gravity * h_upper_average * - (b_jump + h_lower_jump) - f5 += normal_direction_ll[1] * equations.gravity * h_lower_average * - (b_jump + - equations.r * h_upper_jump) - f6 += normal_direction_ll[2] * equations.gravity * h_lower_average * - (b_jump + - equations.r * h_upper_jump) - - # Continuity equations do not have a nonconservative flux - f1 = f4 = zero(eltype(u_ll)) - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - -""" - flux_fjordholm_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations2D) - -Total energy conservative (mathematical entropy for two-layer shallow water equations). When the -bottom topography is nonzero this should only be used as a surface flux otherwise the scheme will -not be well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). - -Details are available in Eq. (4.1) in the paper: -- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011) - Well-balanced and energy stable schemes for the shallow water equations with discontinuous - topography [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) -and the application to two layers is shown in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) - - # Average each factor of products in flux - h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) - v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) - v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) - v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) - v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.25 * equations.gravity * (h_upper_ll^2 + h_upper_rr^2) - p2_avg = 0.25 * equations.gravity * (h_lower_ll^2 + h_lower_rr^2) - - # Calculate fluxes depending on orientation - if orientation == 1 - f1 = h_upper_avg * v1_upper_avg - f2 = f1 * v1_upper_avg + p1_avg - f3 = f1 * v2_upper_avg - f4 = h_lower_avg * v1_lower_avg - f5 = f4 * v1_lower_avg + p2_avg - f6 = f4 * v2_lower_avg - else - f1 = h_upper_avg * v2_upper_avg - f2 = f1 * v1_upper_avg - f3 = f1 * v2_upper_avg + p1_avg - f4 = h_lower_avg * v2_lower_avg - f5 = f4 * v1_lower_avg - f6 = f4 * v2_lower_avg + p2_avg - end - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - -@inline function flux_fjordholm_etal(u_ll, u_rr, - normal_direction::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) - - # Compute velocity in normal direction - v_upper_dot_n_ll = v1_upper_ll * normal_direction[1] + - v2_upper_ll * normal_direction[2] - v_upper_dot_n_rr = v1_upper_rr * normal_direction[1] + - v2_upper_rr * normal_direction[2] - v_lower_dot_n_ll = v1_lower_ll * normal_direction[1] + - v2_lower_ll * normal_direction[2] - v_lower_dot_n_rr = v1_lower_rr * normal_direction[1] + - v2_lower_rr * normal_direction[2] - - # Average each factor of products in flux - h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) - v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) - v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) - v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) - v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.25 * equations.gravity * (h_upper_ll^2 + h_upper_rr^2) - p2_avg = 0.25 * equations.gravity * (h_lower_ll^2 + h_lower_rr^2) - v_upper_dot_n_avg = 0.5 * (v_upper_dot_n_ll + v_upper_dot_n_rr) - v_lower_dot_n_avg = 0.5 * (v_lower_dot_n_ll + v_lower_dot_n_rr) - - # Calculate fluxes depending on normal_direction - f1 = h_upper_avg * v_upper_dot_n_avg - f2 = f1 * v1_upper_avg + p1_avg * normal_direction[1] - f3 = f1 * v2_upper_avg + p1_avg * normal_direction[2] - f4 = h_lower_avg * v_lower_dot_n_avg - f5 = f4 * v1_lower_avg + p2_avg * normal_direction[1] - f6 = f4 * v2_lower_avg + p2_avg * normal_direction[2] - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations2D) + flux_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. -When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. -The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow +When the bottom topography is nonzero this scheme will be well-balanced when used with the +nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the +two-layer shallow water equations the flux that is described in the paper for the normal shallow water equations is used within each layer. Further details are available in Theorem 1 of the paper: @@ -637,24 +446,24 @@ Further details are available in Theorem 1 of the paper: v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p2_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr # Calculate fluxes depending on orientation if orientation == 1 f1 = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) - f2 = f1 * v1_upper_avg + p1_avg + f2 = f1 * v1_upper_avg + p_upper_avg f3 = f1 * v2_upper_avg f4 = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) - f5 = f4 * v1_lower_avg + p2_avg + f5 = f4 * v1_lower_avg + p_lower_avg f6 = f4 * v2_lower_avg else f1 = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) f2 = f1 * v1_upper_avg - f3 = f1 * v2_upper_avg + p1_avg + f3 = f1 * v2_upper_avg + p_upper_avg f4 = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) f5 = f4 * v1_lower_avg - f6 = f4 * v2_lower_avg + p2_avg + f6 = f4 * v2_lower_avg + p_lower_avg end return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) @@ -676,8 +485,8 @@ end v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p2_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr h_v1_upper_avg = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) h_v2_upper_avg = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) h_v1_lower_avg = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) @@ -685,38 +494,36 @@ end # Calculate fluxes depending on normal_direction f1 = h_v1_upper_avg * normal_direction[1] + h_v2_upper_avg * normal_direction[2] - f2 = f1 * v1_upper_avg + p1_avg * normal_direction[1] - f3 = f1 * v2_upper_avg + p1_avg * normal_direction[2] + f2 = f1 * v1_upper_avg + p_upper_avg * normal_direction[1] + f3 = f1 * v2_upper_avg + p_upper_avg * normal_direction[2] f4 = h_v1_lower_avg * normal_direction[1] + h_v2_lower_avg * normal_direction[2] - f5 = f4 * v1_lower_avg + p2_avg * normal_direction[1] - f6 = f4 * v2_lower_avg + p2_avg * normal_direction[2] + f5 = f4 * v1_lower_avg + p_lower_avg * normal_direction[1] + f6 = f4 * v2_lower_avg + p_lower_avg * normal_direction[2] return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) end """ - flux_es_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations1D) - -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative -[`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy -variables. - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) -Energy conservative and stable schemes for the two-layer shallow water equations. -[DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. + flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) + +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ -@inline function flux_es_fjordholm_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations2D) +@inline function flux_es_ersing_etal(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) # Compute entropy conservative flux but without the bottom topography - f_ec = flux_fjordholm_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations) + f_ec = flux_wintermeyer_etal(u_ll, u_rr, + orientation_or_normal_direction, + equations) # Get maximum signal velocity λ = max_abs_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations) @@ -926,12 +733,12 @@ end rho_lower = equations.rho_lower v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) - w1 = rho_upper * (equations.gravity * (h_upper + h_lower + b) + - -0.5 * (v1_upper^2 + v2_upper^2)) + w1 = (rho_upper * (equations.gravity * (h_upper + h_lower + b) + + -0.5 * (v1_upper^2 + v2_upper^2))) w2 = rho_upper * v1_upper w3 = rho_upper * v2_upper - w4 = rho_lower * (equations.gravity * (equations.r * h_upper + h_lower + b) + - -0.5 * (v1_lower^2 + v2_lower^2)) + w4 = (rho_lower * (equations.gravity * (equations.r * h_upper + h_lower + b) + + -0.5 * (v1_lower^2 + v2_lower^2))) w5 = rho_lower * v1_lower w6 = rho_lower * v2_lower return SVector(w1, w2, w3, w4, w5, w6, b) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index fd69a0c1d0e..2269e858928 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -96,6 +96,29 @@ end end end +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254838, + 1.6657566141935285e-14, + 0.10416666834254838, + ], + linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), @@ -169,6 +192,33 @@ end end end +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.005774284062933275, + 0.017408601639513584, + 4.43649172561843e-5, + ], + linf=[ + 0.01639116193303547, + 0.05102877460799604, + 9.098379777450205e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl index a504f4f93a6..180fb3ec3b3 100644 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ b/test/test_tree_1d_shallowwater_twolayer.jl @@ -10,95 +10,65 @@ include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @testset "Shallow Water Two layer" begin -#! format: noindent - -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0050681532925156945, 0.002089013899370176, - 0.005105544300292713, 0.002526442122643306, - 0.0004744186597732706], - linf=[0.022256679217306008, 0.005421833004652266, - 0.02233993939574197, 0.008765261497422516, - 0.0008992474511784199], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_convergence.jl"), + l2=[0.005012009872109003, 0.002091035326731071, + 0.005049271397924551, + 0.0024633066562966574, 0.0004744186597732739], + linf=[0.0213772149343594, 0.005385752427290447, + 0.02175023787351349, + 0.008212004668840978, 0.0008992474511784199], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl with flux_es_fjordholm_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0027681377074701345, 0.0018007543202559165, - 0.0028036917433720576, - 0.0013980358596935737, 0.0004744186597732706], - linf=[0.005699303919826093, 0.006432952918256296, - 0.0058507082844360125, 0.002717615543961216, - 0.0008992474511784199], - surface_flux=(flux_es_fjordholm_etal, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[8.949288784402005e-16, 4.0636427176237915e-17, + 0.001002881985401548, + 2.133351105037203e-16, 0.0010028819854016578], + linf=[2.6229018956769323e-15, 1.878451903240623e-16, + 0.005119880996670156, + 8.003199803957679e-16, 0.005119880996670666], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_well_balanced.jl"), - l2=[8.949288784402005e-16, 4.0636427176237915e-17, - 0.001002881985401548, - 2.133351105037203e-16, 0.0010028819854016578], - linf=[2.6229018956769323e-15, 1.878451903240623e-16, - 0.005119880996670156, - 8.003199803957679e-16, 0.005119880996670666], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_dam_break.jl"), + l2=[0.1000774903431289, 0.5670692949571057, 0.08764242501014498, + 0.45412307886094555, 0.013638618139749523], + linf=[0.586718937495144, 2.1215606128311584, 0.5185911311186155, + 1.820382495072612, 0.5], + surface_flux=(flux_lax_friedrichs, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end end -@trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_dam_break.jl"), - l2=[0.10010269243463918, 0.5668733957648654, - 0.08759617327649398, - 0.4538443183566172, 0.013638618139749523], - linf=[ - 0.5854202777756559, - 2.1278930820498934, - 0.5193686074348809, - 1.8071213168086229, - 0.5, - ], - surface_flux=(flux_lax_friedrichs, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end -end - end # module diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index d280e380192..58db7c5f35f 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -116,6 +116,35 @@ end end end +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987146, + 1.0323158914614244e-14, + 1.0276096319430528e-14, + 0.9130579602987147, + ], + linf=[ + 2.11306203761566, + 4.063916419044386e-14, + 3.694484044448245e-14, + 2.1130620376156584, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), @@ -219,6 +248,35 @@ end end end +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.002471853426064005, + 0.05619168608950033, + 0.11844727575152562, + 6.274146767730281e-5, + ], + linf=[ + 0.014332922987500218, + 0.2141204806174546, + 0.5392313755637872, + 0.0001819675955490041, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_conical_island.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_conical_island.jl"), l2=[ diff --git a/test/test_tree_2d_shallowwater_twolayer.jl b/test/test_tree_2d_shallowwater_twolayer.jl index 5959e7ed882..802bf4e021c 100644 --- a/test/test_tree_2d_shallowwater_twolayer.jl +++ b/test/test_tree_2d_shallowwater_twolayer.jl @@ -10,105 +10,79 @@ include("test_trixi.jl") EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @testset "Two-Layer Shallow Water" begin -#! format: noindent - -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0004040147445601598, 0.005466848793475609, - 0.006149138398472166, 0.0002908599437447256, - 0.003011817461911792, 0.0026806180089700674, - 8.873630921431545e-6], - linf=[0.002822006686981293, 0.014859895905040332, - 0.017590546190827894, 0.0016323702636176218, - 0.009361402900653015, 0.008411036357379165, - 3.361991620143279e-5], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_convergence.jl"), + l2=[0.0004016779699408397, 0.005466339651545468, + 0.006148841330156112, + 0.0002882339012602492, 0.0030120142442780313, + 0.002680752838455618, + 8.873630921431545e-6], + linf=[0.002788654460984752, 0.01484602033450666, + 0.017572229756493973, + 0.0016010835493927011, 0.009369847995372549, + 0.008407961775489636, + 3.361991620143279e-5], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl with flux_es_fjordholm_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.00024709443131137236, 0.0019215286339769443, - 0.0023833298173254447, - 0.00021258247976270914, 0.0011299428031136195, - 0.0009191313765262401, - 8.873630921431545e-6], - linf=[0.0016099763244645793, 0.007659242165565017, - 0.009123320235427057, - 0.0013496983982568267, 0.0035573687287770994, - 0.00296823235874899, - 3.361991620143279e-5], - surface_flux=(flux_es_fjordholm_etal, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[3.2935164267930016e-16, 4.6800825611195103e-17, + 4.843057532147818e-17, + 0.0030769233188015013, 1.4809161150389857e-16, + 1.509071695038043e-16, + 0.0030769233188014935], + linf=[2.248201624865942e-15, 2.346382070278936e-16, + 2.208565017494899e-16, + 0.026474051138910493, 9.237568031609006e-16, + 7.520758026187046e-16, + 0.026474051138910267], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_well_balanced.jl"), - l2=[3.2935164267930016e-16, 4.6800825611195103e-17, - 4.843057532147818e-17, - 0.0030769233188015013, 1.4809161150389857e-16, - 1.509071695038043e-16, - 0.0030769233188014935], - linf=[2.248201624865942e-15, 2.346382070278936e-16, - 2.208565017494899e-16, - 0.026474051138910493, 9.237568031609006e-16, - 7.520758026187046e-16, - 0.026474051138910267], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_well_balanced with flux_lax_friedrichs.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[2.0525741072929735e-16, 6.000589392730905e-17, + 6.102759428478984e-17, + 0.0030769233188014905, 1.8421386173122792e-16, + 1.8473184927121752e-16, + 0.0030769233188014935], + linf=[7.355227538141662e-16, 2.960836949170518e-16, + 4.2726562436938764e-16, + 0.02647405113891016, 1.038795478061861e-15, + 1.0401789378532516e-15, + 0.026474051138910267], + surface_flux=(flux_lax_friedrichs, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end end -@trixi_testset "elixir_shallowwater_twolayer_well_balanced with flux_lax_friedrichs.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_well_balanced.jl"), - l2=[2.0525741072929735e-16, 6.000589392730905e-17, - 6.102759428478984e-17, - 0.0030769233188014905, 1.8421386173122792e-16, - 1.8473184927121752e-16, - 0.0030769233188014935], - linf=[7.355227538141662e-16, 2.960836949170518e-16, - 4.2726562436938764e-16, - 0.02647405113891016, 1.038795478061861e-15, - 1.0401789378532516e-15, - 0.026474051138910267], - surface_flux=(flux_lax_friedrichs, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end -end - end # module diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 26483931cf3..d4416ac5b6a 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -349,6 +349,35 @@ end end end +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839083, + 2.590643638636187e-12, + 1.0945471514840143e-12, + 1.2164292510839079, + ], + linf=[ + 1.5138512282315792, + 5.0276441977281156e-11, + 1.9816934589292803e-11, + 1.513851228231574, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ @@ -402,6 +431,35 @@ end end end +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011196687776346434, + 0.044562672453443995, + 0.014306265289763618, + 5.089218476759981e-6, + ], + linf=[ + 0.007825021762002393, + 0.348550815397918, + 0.1115517935018282, + 2.6407324614341476e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ @@ -535,15 +593,15 @@ end @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0007953969898161991, 0.00882074628714633, - 0.0024322572528892934, - 0.0007597425017400447, 0.004501238950166439, - 0.0015784803573661104, + l2=[0.0007935561625451243, 0.008825315509943844, + 0.002429969315645897, + 0.0007580145888686304, 0.004495741879625235, + 0.0015758146898767814, 6.849532064729749e-6], - linf=[0.00592559068081977, 0.08072451118697077, - 0.0344854497419107, 0.005892196680485795, - 0.04262651217675306, 0.014006223513881366, - 2.5829318284764646e-5], + linf=[0.0059205195991136605, 0.08072126590166251, + 0.03463806075399023, + 0.005884818649227186, 0.042658506561995546, + 0.014125956138838602, 2.5829318284764646e-5], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -582,17 +640,17 @@ end @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_dam_break.jl"), - l2=[0.012471300561905669, 0.012363413819726868, - 0.0009541478004413331, - 0.09120260327331643, 0.015269590815749993, - 0.0012064657396853422, - 0.09991983966647647], - linf=[0.04497814714937959, 0.03286959000796511, - 0.010746094385294369, - 0.11138723974511211, 0.03640850605444494, - 0.014368386516056392, 0.10000000000000003], + l2=[0.012447632879122346, 0.012361250464676683, + 0.0009551519536340908, + 0.09119400061322577, 0.015276216721920347, + 0.0012126995108983853, 0.09991983966647647], + linf=[0.044305765721807444, 0.03279620980615845, + 0.010754320388190101, + 0.111309922939555, 0.03663360204931427, + 0.014332822306649284, + 0.10000000000000003], surface_flux=(flux_lax_friedrichs, - flux_nonconservative_fjordholm_etal), + flux_nonconservative_ersing_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 5fd3d73630fd973f9216fab3c40c1c52509cca68 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 10 Nov 2023 16:02:11 +0100 Subject: [PATCH 16/17] HLL MHD Breaking changes for Consistency (#1708) * Breaking changes HLL MHD * format * format examples * hlle * fix * news, tests, example changes * fmt * remove left-right-biased flux from test * Set version to v0.6.0 * Migrate neural network-based indicators to new repository (#1701) * Remove all neural network indicator stuff from `src/` * Migrate neural network tests * Migrate neural network examples * Migrate test dependencies * Update NEWS.md * Fix typo * Remove Requires.jl-based use of Flux.jl * Fix formatting * Add migration of indicators to section with breaking changes --------- Co-authored-by: Hendrik Ranocha * fix hlle noncartesian 2d * remove parantheses * correct test vals * Make parabolic terms nonexperimental (#1714) * Make parabolic terms non-experimental * Make NSE a separate item * Add MPI to supported features * Mention that parabolic terms are now officially supported in NEWS.md Co-authored-by: Hendrik Ranocha * Deprecate some `DGMultiMesh` constructors (#1709) * remove previously deprecated functions * fix typo in NEWS.md about deprecation vs removal * fix literate tutorial * removing other deprecation * format * Revert "fix typo in NEWS.md about deprecation vs removal" This reverts commit 6b03020a8c881a86550484891a0f53bca018e447. * add gradient variable type parameter to `AbstractEquationsParabolic` (#1409) * add gradient variable type parameter * fix parabolic literate test * remove trailing comment * remove unnecessary abstract type * move gradient variable structs * formatting * fix dropped changes * try to fix doc tests * fixing navier stokes 1D * formatting * remove duplicate GradientVariablesPrimitive/Entropy definition * update news * bring downloads back * fix failing test * fmt --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha --- NEWS.md | 5 ++ .../p4est_2d_dgsem/elixir_mhd_alfven_wave.jl | 4 +- .../elixir_mhd_alfven_wave_nonconforming.jl | 4 +- .../elixir_mhd_alfven_wave.jl | 4 +- .../t8code_2d_dgsem/elixir_mhd_alfven_wave.jl | 2 +- .../elixir_mhd_briowu_shock_tube.jl | 2 +- .../elixir_mhd_ryujones_shock_tube.jl | 2 +- .../elixir_mhd_shu_osher_shock_tube.jl | 2 +- .../elixir_mhd_alfven_wave_mortar.jl | 4 +- .../elixir_mhd_alfven_wave_mortar.jl | 4 +- .../elixir_mhd_alfven_wave.jl | 4 +- src/equations/ideal_glm_mhd_1d.jl | 22 +++++- src/equations/ideal_glm_mhd_2d.jl | 60 ++++++++++++++-- src/equations/ideal_glm_mhd_3d.jl | 68 +++++++++++++++++-- test/test_tree_2d_mhd.jl | 3 +- test/test_tree_3d_mhd.jl | 3 +- test/test_unit.jl | 19 +++--- 17 files changed, 177 insertions(+), 35 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5d258fa65bb..54fbb90b8fc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,11 @@ for human readability. #### Changed +- The wave speed estimates for `flux_hll`, `FluxHLL()` are now consistent across equations. + In particular, the functions `min_max_speed_naive`, `min_max_speed_einfeldt` are now + conceptually identical across equations. + Users, who have been using `flux_hll` for MHD have now to use `flux_hlle` in order to use the + Einfeldt wave speed estimate. - Parabolic diffusion terms are now officially supported and not marked as experimental anymore. diff --git a/examples/p4est_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/p4est_2d_dgsem/elixir_mhd_alfven_wave.jl index 93db7bfdbaf..b2b402a25f6 100644 --- a/examples/p4est_2d_dgsem/elixir_mhd_alfven_wave.jl +++ b/examples/p4est_2d_dgsem/elixir_mhd_alfven_wave.jl @@ -12,7 +12,9 @@ initial_condition = initial_condition_convergence_test # Get the DG approximation space volume_flux = (flux_central, flux_nonconservative_powell) -solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_powell), +solver = DGSEM(polydeg = 4, + surface_flux = (flux_hlle, + flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) coordinates_min = (0.0, 0.0) diff --git a/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_nonconforming.jl b/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_nonconforming.jl index 6a62368ef99..12ddf9e4a5f 100644 --- a/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_nonconforming.jl +++ b/examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_nonconforming.jl @@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations3D(5 / 3) initial_condition = initial_condition_convergence_test volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) -solver = DGSEM(polydeg = 3, surface_flux = (flux_hll, flux_nonconservative_powell), +solver = DGSEM(polydeg = 3, + surface_flux = (flux_hlle, + flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) coordinates_min = (-1.0, -1.0, -1.0) diff --git a/examples/structured_3d_dgsem/elixir_mhd_alfven_wave.jl b/examples/structured_3d_dgsem/elixir_mhd_alfven_wave.jl index 6eb35078ef4..03f50ff3e55 100644 --- a/examples/structured_3d_dgsem/elixir_mhd_alfven_wave.jl +++ b/examples/structured_3d_dgsem/elixir_mhd_alfven_wave.jl @@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations3D(5 / 3) initial_condition = initial_condition_convergence_test volume_flux = (flux_central, flux_nonconservative_powell) -solver = DGSEM(polydeg = 5, surface_flux = (flux_hll, flux_nonconservative_powell), +solver = DGSEM(polydeg = 5, + surface_flux = (flux_hlle, + flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # Create the mesh diff --git a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl index 8040f79cafd..1e2362a123c 100644 --- a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl +++ b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl @@ -11,7 +11,7 @@ initial_condition = initial_condition_convergence_test # Get the DG approximation space volume_flux = (flux_central, flux_nonconservative_powell) -solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_powell), +solver = DGSEM(polydeg = 4, surface_flux = (flux_hlle, flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) coordinates_min = (0.0, 0.0) diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl index bb294af68cb..4d537508b47 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl @@ -32,7 +32,7 @@ initial_condition = initial_condition_briowu_shock_tube boundary_conditions = BoundaryConditionDirichlet(initial_condition) -surface_flux = flux_hll +surface_flux = flux_hlle volume_flux = flux_derigs_etal basis = LobattoLegendreBasis(4) diff --git a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl index b6f856fbc64..a7d7689a806 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_ryujones_shock_tube.jl @@ -39,7 +39,7 @@ initial_condition = initial_condition_ryujones_shock_tube boundary_conditions = BoundaryConditionDirichlet(initial_condition) -surface_flux = flux_hll +surface_flux = flux_hlle volume_flux = flux_hindenlang_gassner basis = LobattoLegendreBasis(3) diff --git a/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl index 8c0317277b0..689537ebdd4 100644 --- a/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl +++ b/examples/tree_1d_dgsem/elixir_mhd_shu_osher_shock_tube.jl @@ -61,7 +61,7 @@ initial_condition = initial_condition_shu_osher_shock_tube boundary_conditions = BoundaryConditionDirichlet(initial_condition) -surface_flux = flux_hll +surface_flux = flux_hlle volume_flux = flux_hindenlang_gassner basis = LobattoLegendreBasis(4) diff --git a/examples/tree_2d_dgsem/elixir_mhd_alfven_wave_mortar.jl b/examples/tree_2d_dgsem/elixir_mhd_alfven_wave_mortar.jl index 229360f266e..0200b591844 100644 --- a/examples/tree_2d_dgsem/elixir_mhd_alfven_wave_mortar.jl +++ b/examples/tree_2d_dgsem/elixir_mhd_alfven_wave_mortar.jl @@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations2D(gamma) initial_condition = initial_condition_convergence_test volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) -solver = DGSEM(polydeg = 3, surface_flux = (flux_hll, flux_nonconservative_powell), +solver = DGSEM(polydeg = 3, + surface_flux = (flux_hlle, + flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) coordinates_min = (0.0, 0.0) diff --git a/examples/tree_3d_dgsem/elixir_mhd_alfven_wave_mortar.jl b/examples/tree_3d_dgsem/elixir_mhd_alfven_wave_mortar.jl index 3ce166a7fa7..2fa22d535d7 100644 --- a/examples/tree_3d_dgsem/elixir_mhd_alfven_wave_mortar.jl +++ b/examples/tree_3d_dgsem/elixir_mhd_alfven_wave_mortar.jl @@ -10,7 +10,9 @@ equations = IdealGlmMhdEquations3D(5 / 3) initial_condition = initial_condition_convergence_test volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) -solver = DGSEM(polydeg = 3, surface_flux = (flux_hll, flux_nonconservative_powell), +solver = DGSEM(polydeg = 3, + surface_flux = (flux_hlle, + flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) coordinates_min = (-1.0, -1.0, -1.0) diff --git a/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl index fdafed98c8d..3ed3e828ca8 100644 --- a/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl +++ b/examples/unstructured_2d_dgsem/elixir_mhd_alfven_wave.jl @@ -11,7 +11,9 @@ equations = IdealGlmMhdEquations2D(gamma) initial_condition = initial_condition_convergence_test volume_flux = (flux_central, flux_nonconservative_powell) -solver = DGSEM(polydeg = 7, surface_flux = (flux_hll, flux_nonconservative_powell), +solver = DGSEM(polydeg = 7, + surface_flux = (flux_hlle, + flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # Get the unstructured quad mesh from a file (downloads the file if not available locally) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 7e5c94c7bc3..85030e8a5ad 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -277,6 +277,22 @@ end λ_max = max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr) end +# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations1D) + rho_ll, rho_v1_ll, _ = u_ll + rho_rr, rho_v1_rr, _ = u_rr + + # Calculate primitive variables + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + + λ_min = v1_ll - calc_fast_wavespeed(u_ll, orientation, equations) + λ_max = v1_rr + calc_fast_wavespeed(u_rr, orientation, equations) + + return λ_min, λ_max +end + # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D) @@ -298,15 +314,15 @@ end end """ - min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D) + min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D) Calculate minimum and maximum wave speeds for HLL-type fluxes as in - Li (2005) An HLLC Riemann solver for magneto-hydrodynamics [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020). """ -@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, - equations::IdealGlmMhdEquations1D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations1D) rho_ll, rho_v1_ll, _ = u_ll rho_rr, rho_v1_rr, _ = u_rr diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index e8de0cedde1..43d1991e34b 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -775,6 +775,56 @@ end return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr) end +# Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations2D) + rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr + + # Calculate primitive velocity variables + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + # Approximate the left-most and right-most eigenvalues in the Riemann fan + if orientation == 1 # x-direction + λ_min = v1_ll - calc_fast_wavespeed(u_ll, orientation, equations) + λ_max = v1_rr + calc_fast_wavespeed(u_rr, orientation, equations) + else # y-direction + λ_min = v2_ll - calc_fast_wavespeed(u_ll, orientation, equations) + λ_max = v2_rr + calc_fast_wavespeed(u_rr, orientation, equations) + end + + return λ_min, λ_max +end + +@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::IdealGlmMhdEquations2D) + rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr + + # Calculate primitive velocity variables + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + v_normal_ll = (v1_ll * normal_direction[1] + v2_ll * normal_direction[2]) + v_normal_rr = (v1_rr * normal_direction[1] + v2_rr * normal_direction[2]) + + c_f_ll = calc_fast_wavespeed(u_ll, normal_direction, equations) + c_f_rr = calc_fast_wavespeed(u_rr, normal_direction, equations) + + # Estimate the min/max eigenvalues in the normal direction + λ_min = min(v_normal_ll - c_f_ll, v_normal_rr - c_f_rr) + λ_max = max(v_normal_rr + c_f_rr, v_normal_rr + c_f_rr) + + return λ_min, λ_max +end + # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D) @@ -833,15 +883,15 @@ end end """ - min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D) + min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D) Calculate minimum and maximum wave speeds for HLL-type fluxes as in - Li (2005) An HLLC Riemann solver for magneto-hydrodynamics [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020). """ -@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, - equations::IdealGlmMhdEquations2D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations2D) rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr @@ -870,8 +920,8 @@ Calculate minimum and maximum wave speeds for HLL-type fluxes as in return λ_min, λ_max end -@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, - equations::IdealGlmMhdEquations2D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, normal_direction::AbstractVector, + equations::IdealGlmMhdEquations2D) rho_ll, rho_v1_ll, rho_v2_ll, _ = u_ll rho_rr, rho_v1_rr, rho_v2_rr, _ = u_rr diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 09990837706..321e501b051 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -670,6 +670,64 @@ end return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr) end +# Calculate estimate for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations3D) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr + + # Calculate primitive variables and speed of sound + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + v3_rr = rho_v3_rr / rho_rr + + # Approximate the left-most and right-most eigenvalues in the Riemann fan + if orientation == 1 # x-direction + λ_min = v1_ll - calc_fast_wavespeed(u_ll, orientation, equations) + λ_max = v1_rr + calc_fast_wavespeed(u_rr, orientation, equations) + elseif orientation == 2 # y-direction + λ_min = v2_ll - calc_fast_wavespeed(u_ll, orientation, equations) + λ_max = v2_rr + calc_fast_wavespeed(u_rr, orientation, equations) + else # z-direction + λ_min = v3_ll - calc_fast_wavespeed(u_ll, orientation, equations) + λ_max = v3_rr + calc_fast_wavespeed(u_rr, orientation, equations) + end + + return λ_min, λ_max +end + +@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::IdealGlmMhdEquations3D) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr + + # Calculate primitive velocity variables + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + v3_rr = rho_v3_rr / rho_rr + + v_normal_ll = (v1_ll * normal_direction[1] + + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3]) + v_normal_rr = (v1_rr * normal_direction[1] + + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3]) + + # Estimate the min/max eigenvalues in the normal direction + λ_min = v_normal_ll - calc_fast_wavespeed(u_ll, normal_direction, equations) + λ_max = v_normal_rr + calc_fast_wavespeed(u_rr, normal_direction, equations) + + return λ_min, λ_max +end + # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations3D) @@ -742,15 +800,15 @@ end end """ - min_max_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations::IdealGlmMhdEquations3D) + min_max_speed_einfeldt(u_ll, u_rr, orientation_or_normal_direction, equations::IdealGlmMhdEquations3D) Calculate minimum and maximum wave speeds for HLL-type fluxes as in - Li (2005) An HLLC Riemann solver for magneto-hydrodynamics [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020) """ -@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, - equations::IdealGlmMhdEquations3D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdEquations3D) rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr @@ -787,8 +845,8 @@ Calculate minimum and maximum wave speeds for HLL-type fluxes as in return λ_min, λ_max end -@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, - equations::IdealGlmMhdEquations3D) +@inline function min_max_speed_einfeldt(u_ll, u_rr, normal_direction::AbstractVector, + equations::IdealGlmMhdEquations3D) rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, _ = u_ll rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, _ = u_rr diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index bd6a95bba50..953c077c0a3 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -208,7 +208,8 @@ end 0.04879208429337193, ], tspan=(0.0, 0.06), - surface_flux=(flux_hll, flux_nonconservative_powell)) + surface_flux=(flux_hlle, + flux_nonconservative_powell)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_tree_3d_mhd.jl b/test/test_tree_3d_mhd.jl index 7ce5ef1d18f..e75685f0b43 100644 --- a/test/test_tree_3d_mhd.jl +++ b/test/test_tree_3d_mhd.jl @@ -240,7 +240,8 @@ end 1000 end end, - surface_flux=(flux_hll, flux_nonconservative_powell), + surface_flux=(flux_hlle, + flux_nonconservative_powell), volume_flux=(flux_central, flux_nonconservative_powell), coordinates_min=(0.0, 0.0, 0.0), coordinates_max=(1.0, 1.0, 1.0), diff --git a/test/test_unit.jl b/test/test_unit.jl index 609100793ba..1d768c0df69 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -952,14 +952,12 @@ end end @timed_testset "Consistency check for HLLE flux: MHD" begin - # Note: min_max_speed_naive for MHD is essentially min_max_speed_einfeldt - equations = IdealGlmMhdEquations1D(1.4) u_values = [SVector(1.0, 0.4, -0.5, 0.1, 1.0, 0.1, -0.2, 0.1), SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2)] for u in u_values - @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations) + @test flux_hlle(u, u, 1, equations) ≈ flux(u, 1, equations) end equations = IdealGlmMhdEquations2D(1.4, 5.0) #= c_h =# @@ -973,11 +971,11 @@ end SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)] for u in u_values, orientation in orientations - @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations) + @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations) end for u in u_values, normal_direction in normal_directions - @test flux_hll(u, u, normal_direction, equations) ≈ + @test flux_hlle(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) end @@ -993,11 +991,11 @@ end SVector(1.5, -0.2, 0.1, 0.2, 5.0, -0.1, 0.1, 0.2, 0.2)] for u in u_values, orientation in orientations - @test flux_hll(u, u, orientation, equations) ≈ flux(u, orientation, equations) + @test flux_hlle(u, u, orientation, equations) ≈ flux(u, orientation, equations) end for u in u_values, normal_direction in normal_directions - @test flux_hll(u, u, normal_direction, equations) ≈ + @test flux_hlle(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) end end @@ -1244,8 +1242,8 @@ end fluxes = [ flux_central, flux_hindenlang_gassner, - flux_hll, FluxHLL(min_max_speed_davis), + flux_hlle, ] for f_std in fluxes @@ -1271,8 +1269,8 @@ end fluxes = [ flux_central, flux_hindenlang_gassner, - flux_hll, FluxHLL(min_max_speed_davis), + flux_hlle, ] for f_std in fluxes @@ -1322,7 +1320,8 @@ end dg = DGMulti(polydeg = 1, element_type = Line(), approximation_type = Polynomial(), surface_integral = SurfaceIntegralWeakForm(flux_central), volume_integral = VolumeIntegralFluxDifferencing(flux_central)) - mesh = DGMultiMesh(dg, cells_per_dimension = (1,), periodicity = false) + cells_per_dimension = (1,) + mesh = DGMultiMesh(dg, cells_per_dimension, periodicity = false) @test mesh.boundary_faces[:entire_boundary] == [1, 2] end From 1b75f5ecbd0b21d5cdce25f3e4f1884a06727fbb Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 12 Nov 2023 07:16:33 +0100 Subject: [PATCH 17/17] set development version to v0.6.1-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e95e5347ee3..da618dc78c1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.6.0" +version = "0.6.1-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"