From d86a0f47bebe0945b3a70143015adb0fe95e6174 Mon Sep 17 00:00:00 2001 From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> Date: Tue, 26 Mar 2024 08:22:48 +0100 Subject: [PATCH 01/13] Use @batch reduction functionality for subcell bounds check (#1888) * Use @batch reduction functionality for bounds check * fmt * Adapt compat bound for Polyester.jl * Add comments --- Project.toml | 2 +- docs/src/performance.md | 11 ---- .../subcell_bounds_check_2d.jl | 66 ++++++++----------- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 11 +--- 4 files changed, 32 insertions(+), 58 deletions(-) diff --git a/Project.toml b/Project.toml index 27df49ed4fa..3db90c9928e 100644 --- a/Project.toml +++ b/Project.toml @@ -75,7 +75,7 @@ MuladdMacro = "0.2.2" Octavian = "0.3.21" OffsetArrays = "1.12" P4est = "0.4.9" -Polyester = "0.7.5" +Polyester = "0.7.10" PrecompileTools = "1.1" Preferences = "1.3" Printf = "1" diff --git a/docs/src/performance.md b/docs/src/performance.md index 40970e58c5c..9f81d3c3d8e 100644 --- a/docs/src/performance.md +++ b/docs/src/performance.md @@ -282,14 +282,3 @@ requires. It can thus be seen as a proxy for "energy used" and, as an extension, timing result, you need to set the analysis interval such that the `AnalysisCallback` is invoked at least once during the course of the simulation and discard the first PID value. - -## Performance issues with multi-threaded reductions -[False sharing](https://en.wikipedia.org/wiki/False_sharing) is a known performance issue -for systems with distributed caches. It also occurred for the implementation of a thread -parallel bounds checking routine for the subcell IDP limiting -in [PR #1736](https://github.com/trixi-framework/Trixi.jl/pull/1736). -After some [testing and discussion](https://github.com/trixi-framework/Trixi.jl/pull/1736#discussion_r1423881895), -it turned out that initializing a vector of length `n * Threads.nthreads()` and only using every -n-th entry instead of a vector of length `Threads.nthreads()` fixes the problem. -Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)`. -Now, the bounds checking routine of the IDP limiting scales as hoped. diff --git a/src/callbacks_stage/subcell_bounds_check_2d.jl b/src/callbacks_stage/subcell_bounds_check_2d.jl index 19d73968c9a..3a56ea71f62 100644 --- a/src/callbacks_stage/subcell_bounds_check_2d.jl +++ b/src/callbacks_stage/subcell_bounds_check_2d.jl @@ -12,35 +12,37 @@ (; variable_bounds) = limiter.cache.subcell_limiter_coefficients (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache - # Note: Accessing the threaded memory vector `idp_bounds_delta_local` with - # `deviation = idp_bounds_delta_local[key][Threads.threadid()]` causes critical performance - # issues due to False Sharing. - # Initializing a vector with n times the length and using every n-th entry fixes this - # problem and allows proper scaling: - # `deviation = idp_bounds_delta_local[key][n * Threads.threadid()]` - # Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)` - stride_size = div(128, sizeof(eltype(u))) # = n + # Note: In order to get the maximum deviation from the target bounds, this bounds check + # requires a reduction in every RK stage and for every enabled limiting option. To make + # this Thread-parallel we are using Polyester.jl's (at least v0.7.10) `@batch reduction` + # functionality. + # Although `@threaded` and `@batch` are currently used equivalently in Trixi.jl, we use + # `@batch` here to allow a possible redefinition of `@threaded` without creating errors here. + # See also https://github.com/trixi-framework/Trixi.jl/pull/1888#discussion_r1537785293. if local_minmax for v in limiter.local_minmax_variables_cons v_string = string(v) key_min = Symbol(v_string, "_min") key_max = Symbol(v_string, "_max") - deviation_min_threaded = idp_bounds_delta_local[key_min] - deviation_max_threaded = idp_bounds_delta_local[key_max] - @threaded for element in eachelement(solver, cache) - deviation_min = deviation_min_threaded[stride_size * Threads.threadid()] - deviation_max = deviation_max_threaded[stride_size * Threads.threadid()] + deviation_min = idp_bounds_delta_local[key_min] + deviation_max = idp_bounds_delta_local[key_max] + @batch reduction=((max, deviation_min), (max, deviation_max)) for element in eachelement(solver, + cache) for j in eachnode(solver), i in eachnode(solver) var = u[v, i, j, element] + # Note: We always save the absolute deviations >= 0 and therefore use the + # `max` operator for the lower and upper bound. The different directions of + # upper and lower bound are considered in their calculations with a + # different sign. deviation_min = max(deviation_min, variable_bounds[key_min][i, j, element] - var) deviation_max = max(deviation_max, var - variable_bounds[key_max][i, j, element]) end - deviation_min_threaded[stride_size * Threads.threadid()] = deviation_min - deviation_max_threaded[stride_size * Threads.threadid()] = deviation_max end + idp_bounds_delta_local[key_min] = deviation_min + idp_bounds_delta_local[key_max] = deviation_max end end if positivity @@ -49,40 +51,35 @@ continue end key = Symbol(string(v), "_min") - deviation_threaded = idp_bounds_delta_local[key] - @threaded for element in eachelement(solver, cache) - deviation = deviation_threaded[stride_size * Threads.threadid()] + deviation = idp_bounds_delta_local[key] + @batch reduction=(max, deviation) for element in eachelement(solver, cache) for j in eachnode(solver), i in eachnode(solver) var = u[v, i, j, element] deviation = max(deviation, variable_bounds[key][i, j, element] - var) end - deviation_threaded[stride_size * Threads.threadid()] = deviation end + idp_bounds_delta_local[key] = deviation end for variable in limiter.positivity_variables_nonlinear key = Symbol(string(variable), "_min") - deviation_threaded = idp_bounds_delta_local[key] - @threaded for element in eachelement(solver, cache) - deviation = deviation_threaded[stride_size * Threads.threadid()] + deviation = idp_bounds_delta_local[key] + @batch reduction=(max, deviation) for element in eachelement(solver, cache) for j in eachnode(solver), i in eachnode(solver) var = variable(get_node_vars(u, equations, solver, i, j, element), equations) deviation = max(deviation, variable_bounds[key][i, j, element] - var) end - deviation_threaded[stride_size * Threads.threadid()] = deviation end + idp_bounds_delta_local[key] = deviation end end for (key, _) in idp_bounds_delta_local - # Calculate maximum deviations of all threads - idp_bounds_delta_local[key][stride_size] = maximum(idp_bounds_delta_local[key][stride_size * i] - for i in 1:Threads.nthreads()) # Update global maximum deviations idp_bounds_delta_global[key] = max(idp_bounds_delta_global[key], - idp_bounds_delta_local[key][stride_size]) + idp_bounds_delta_local[key]) end if save_errors @@ -92,10 +89,8 @@ if local_minmax for v in limiter.local_minmax_variables_cons v_string = string(v) - print(f, ", ", - idp_bounds_delta_local[Symbol(v_string, "_min")][stride_size], - ", ", - idp_bounds_delta_local[Symbol(v_string, "_max")][stride_size]) + print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")], + ", ", idp_bounds_delta_local[Symbol(v_string, "_max")]) end end if positivity @@ -103,21 +98,18 @@ if v in limiter.local_minmax_variables_cons continue end - print(f, ", ", - idp_bounds_delta_local[Symbol(string(v), "_min")][stride_size]) + print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")]) end for variable in limiter.positivity_variables_nonlinear print(f, ", ", - idp_bounds_delta_local[Symbol(string(variable), "_min")][stride_size]) + idp_bounds_delta_local[Symbol(string(variable), "_min")]) end end println(f) end # Reset local maximum deviations for (key, _) in idp_bounds_delta_local - for i in 1:Threads.nthreads() - idp_bounds_delta_local[key][stride_size * i] = zero(eltype(idp_bounds_delta_local[key][stride_size])) - end + idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key])) end end diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 3f7954c8958..9343cee4397 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -18,18 +18,11 @@ function create_cache(limiter::Type{SubcellLimiterIDP}, equations::AbstractEquat # Memory for bounds checking routine with `BoundsCheckCallback`. # Local variable contains the maximum deviation since the last export. - # Using a threaded vector to parallelize bounds check. - idp_bounds_delta_local = Dict{Symbol, Vector{real(basis)}}() + idp_bounds_delta_local = Dict{Symbol, real(basis)}() # Global variable contains the total maximum deviation. idp_bounds_delta_global = Dict{Symbol, real(basis)}() - # Note: False sharing causes critical performance issues on multiple threads when using a vector - # of length `Threads.nthreads()`. Initializing a vector of length `n * Threads.nthreads()` - # and then only using every n-th entry, fixes the problem and allows proper scaling. - # Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)` - stride_size = div(128, sizeof(eltype(basis.nodes))) # = n for key in bound_keys - idp_bounds_delta_local[key] = [zero(real(basis)) - for _ in 1:(stride_size * Threads.nthreads())] + idp_bounds_delta_local[key] = zero(real(basis)) idp_bounds_delta_global[key] = zero(real(basis)) end From cd6fdaaa1443376ad97bc408cafc27d82a647175 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Tue, 26 Mar 2024 08:23:15 +0100 Subject: [PATCH 02/13] abort initial AMR loop after 10 iterations (#1890) Co-authored-by: Hendrik Ranocha --- src/callbacks_step/amr.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 6f57d6647fc..1ab65a3553e 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -138,11 +138,18 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u, t, # iterate until mesh does not change anymore has_changed = amr_callback(integrator, only_refine = amr_callback.adapt_initial_condition_only_refine) + iterations = 1 while has_changed compute_coefficients!(integrator.u, t, semi) u_modified!(integrator, true) has_changed = amr_callback(integrator, only_refine = amr_callback.adapt_initial_condition_only_refine) + iterations = iterations + 1 + if iterations > 10 + @warn "AMR for initial condition did not settle within 10 iterations!\n" * + "Consider adjusting thresholds or setting `adapt_initial_condition_only_refine`." + break + end end end From 27026de6e2be4f0f6655efb68c082515f1383e84 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 26 Mar 2024 10:07:11 +0100 Subject: [PATCH 03/13] set version to v0.7.5 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3db90c9928e..9ca3086306f 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.7.5-pre" +version = "0.7.5" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 909abb4473c95042a87e93e46d443dc381c2b523 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 26 Mar 2024 10:07:23 +0100 Subject: [PATCH 04/13] set development version to v0.7.6-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9ca3086306f..6ff7f29686d 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.7.5" +version = "0.7.6-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 73da92c446c6f60fc09ad1cf5e579bbf5eccb34d Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 27 Mar 2024 15:55:37 +0100 Subject: [PATCH 05/13] Sutherlands Law for temperature dependent viscosity (#1808) * sample 2D * bug fixes * typo * sutherlands 1d * sutherlands 3d * fmt * main merge * test var mu * variable mu * fmt * example using sutherland * comment * example * reset toml * Shorten * Unit test * fmt * Update examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl * remove todo * Apply suggestions from code review Co-authored-by: Andrew Winters * test vals * non-functionalized mu * comments & dosctrings * docstrings * fix tests * avoid if * docstring and comments * Update src/equations/compressible_navier_stokes.jl Co-authored-by: Michael Schlottke-Lakemper * fmt --------- Co-authored-by: Andrew Winters Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> --- .../elixir_navierstokes_lid_driven_cavity.jl | 5 +- ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- .../elixir_navierstokes_lid_driven_cavity.jl | 5 +- ...ixir_navierstokes_lid_driven_cavity_amr.jl | 5 +- .../elixir_navierstokes_blast_wave_amr.jl | 5 +- ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- ...ir_navierstokes_taylor_green_vortex_amr.jl | 5 +- ...lixir_navierstokes_convergence_periodic.jl | 1 - .../elixir_navierstokes_lid_driven_cavity.jl | 5 +- .../elixir_navierstokes_shearlayer_amr.jl | 5 +- ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- ...erstokes_taylor_green_vortex_sutherland.jl | 94 +++++++++++++++++++ ...elixir_navierstokes_taylor_green_vortex.jl | 5 +- src/equations/compressible_navier_stokes.jl | 14 +++ .../compressible_navier_stokes_1d.jl | 26 +++-- .../compressible_navier_stokes_2d.jl | 26 +++-- .../compressible_navier_stokes_3d.jl | 26 +++-- test/test_parabolic_2d.jl | 34 +++++-- test/test_unit.jl | 39 ++++++++ 19 files changed, 243 insertions(+), 72 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl index 7c55cbf0ccf..a612dd0e0d5 100644 --- a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl index dedd8267a3b..9ae90ac47b6 100644 --- a/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index bc28ae6ffb3..728736fe49e 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl index 898366969a4..d9eaf4fdb72 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_lid_driven_cavity_amr.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl index 5df89fbcdf2..d556d0ab70d 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) function initial_condition_3d_blast_wave(x, t, equations::CompressibleEulerEquations3D) diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index d785013f5a9..9c90e4d3218 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl index c15227a1c29..2741f0df174 100644 --- a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl index 33ad0d5271c..eab0840f385 100644 --- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_periodic.jl @@ -5,7 +5,6 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 mu() = 6.25e-4 # equivalent to Re = 1600 diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl index 70d76fc9075..b8e20e27f66 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_lid_driven_cavity.jl @@ -4,12 +4,11 @@ using Trixi ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 0.001 +mu = 0.001 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 4d92ea261e9..a7492bafb47 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 1.0 / 3.0 * 10^(-5) # equivalent to Re = 30,000 +mu = 1.0 / 3.0 * 10^(-4) # equivalent to Re = 30,000 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index a3e38bf6d92..c6e5f0bc40a 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl new file mode 100644 index 00000000000..9598ae184cd --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex_sutherland.jl @@ -0,0 +1,94 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +prandtl_number() = 0.72 + +# Use Sutherland's law for a temperature-dependent viscosity. +# For details, see e.g. +# Frank M. White: Viscous Fluid Flow, 2nd Edition. +# 1991, McGraw-Hill, ISBN, 0-07-069712-4 +# Pages 28 and 29. +@inline function mu(u, equations) + T_ref = 291.15 + + R_specific_air = 287.052874 + T = R_specific_air * Trixi.temperature(u, equations) + + C_air = 120.0 + mu_ref_air = 1.827e-5 + + return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5 +end + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, + Prandtl = prandtl_number()) + +""" + initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D) + +The classical viscous Taylor-Green vortex in 2D. +This forms the basis behind the 3D case found for instance in + - Jonathan R. Bull and Antony Jameson + Simulation of the Compressible Taylor Green Vortex using High-Order Flux Reconstruction Schemes + [DOI: 10.2514/6.2014-3210](https://doi.org/10.2514/6.2014-3210) +""" +function initial_condition_taylor_green_vortex(x, t, + equations::CompressibleEulerEquations2D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) + v2 = -A * cos(x[1]) * sin(x[2]) + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + 1.0 / 4.0 * A^2 * rho * (cos(2 * x[1]) + cos(2 * x[2])) + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0) .* pi +coordinates_max = (1.0, 1.0) .* pi +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 100_000) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-9 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl index 65bd9aa133d..3e54c791ec6 100644 --- a/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl +++ b/examples/tree_3d_dgsem/elixir_navierstokes_taylor_green_vortex.jl @@ -5,12 +5,11 @@ using Trixi ############################################################################### # semidiscretization of the compressible Navier-Stokes equations -# TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 6.25e-4 # equivalent to Re = 1600 +mu = 6.25e-4 # equivalent to Re = 1600 equations = CompressibleEulerEquations3D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu, Prandtl = prandtl_number()) """ diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl index 3059771197c..c530625f226 100644 --- a/src/equations/compressible_navier_stokes.jl +++ b/src/equations/compressible_navier_stokes.jl @@ -62,3 +62,17 @@ Under `GradientVariablesEntropy`, the Navier-Stokes discretization is provably e """ struct GradientVariablesPrimitive end struct GradientVariablesEntropy end + +""" + dynamic_viscosity(u, equations) + +Wrapper for the dynamic viscosity that calls +`dynamic_viscosity(u, equations.mu, equations)`, which dispatches on the type of +`equations.mu`. +For constant `equations.mu`, i.e., `equations.mu` is of `Real`-type it is returned directly. +In all other cases, `equations.mu` is assumed to be a function with arguments +`u` and `equations` and is called with these arguments. +""" +dynamic_viscosity(u, equations) = dynamic_viscosity(u, equations.mu, equations) +dynamic_viscosity(u, mu::Real, equations) = mu +dynamic_viscosity(u, mu::T, equations) where {T} = mu(u, equations) diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index d2c46ecc7d8..3dbdf11c56b 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -21,6 +21,9 @@ the [`CompressibleEulerEquations1D`](@ref). Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. +The viscosity ``\mu`` may be a constant or a function of the current state, e.g., +depending on temperature (Sutherland's law): ``\mu = \mu(T)``. +In the latter case, the function `mu` needs to have the signature `mu(u, equations)`. The particular form of the compressible Navier-Stokes implemented is ```math @@ -80,7 +83,7 @@ where w_2 = \frac{\rho v1}{p},\, w_3 = -\frac{\rho}{p} ``` """ -struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, +struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, Mu, E <: AbstractCompressibleEulerEquations{1}} <: AbstractCompressibleNavierStokesDiffusion{1, 3, GradientVariables} # TODO: parabolic @@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion1D{GradientVariables, RealT <: Real, gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - mu::RealT # viscosity + mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law @@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion1D(equations::CompressibleEulerEquatio gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - μ, Pr = promote(mu, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity - # constant is kappa = gamma μ / ((gamma-1) Pr). + # constant is kappa = gamma μ / ((gamma-1) Prandtl). # Important note! Factor of μ is accounted for later in `flux`. - kappa = gamma * inv_gamma_minus_one / Pr + # This avoids recomputation of kappa for non-constant μ. + kappa = gamma * inv_gamma_minus_one / Prandtl CompressibleNavierStokesDiffusion1D{typeof(gradient_variables), typeof(gamma), + typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, - μ, Pr, kappa, + mu, Prandtl, kappa, equations, gradient_variables) end @@ -159,10 +163,12 @@ function flux(u, gradients, orientation::Integer, # in the implementation q1 = equations.kappa * dTdx - # Constant dynamic viscosity is copied to a variable for readability. - # Offers flexibility for dynamic viscosity via Sutherland's law where it depends - # on temperature and reference values, Ts and Tref such that mu(T) - mu = equations.mu + # In the simplest cases, the user passed in `mu` or `mu()` + # (which returns just a constant) but + # more complex functions like Sutherland's law are possible. + # `dynamic_viscosity` is a helper function that handles both cases + # by dispatching on the type of `equations.mu`. + mu = dynamic_viscosity(u, equations) # viscous flux components in the x-direction f1 = zero(rho) diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 5df7c01ca5c..3256343703a 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -21,6 +21,9 @@ the [`CompressibleEulerEquations2D`](@ref). Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. +The viscosity ``\mu`` may be a constant or a function of the current state, e.g., +depending on temperature (Sutherland's law): ``\mu = \mu(T)``. +In the latter case, the function `mu` needs to have the signature `mu(u, equations)`. The particular form of the compressible Navier-Stokes implemented is ```math @@ -80,7 +83,7 @@ where w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = -\frac{\rho}{p} ``` """ -struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, +struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, Mu, E <: AbstractCompressibleEulerEquations{2}} <: AbstractCompressibleNavierStokesDiffusion{2, 4, GradientVariables} # TODO: parabolic @@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion2D{GradientVariables, RealT <: Real, gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - mu::RealT # viscosity + mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law @@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion2D(equations::CompressibleEulerEquatio gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - μ, Pr = promote(mu, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity - # constant is kappa = gamma μ / ((gamma-1) Pr). + # constant is kappa = gamma μ / ((gamma-1) Prandtl). # Important note! Factor of μ is accounted for later in `flux`. - kappa = gamma * inv_gamma_minus_one / Pr + # This avoids recomputation of kappa for non-constant μ. + kappa = gamma * inv_gamma_minus_one / Prandtl CompressibleNavierStokesDiffusion2D{typeof(gradient_variables), typeof(gamma), + typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, - μ, Pr, kappa, + mu, Prandtl, kappa, equations, gradient_variables) end @@ -168,10 +172,12 @@ function flux(u, gradients, orientation::Integer, q1 = equations.kappa * dTdx q2 = equations.kappa * dTdy - # Constant dynamic viscosity is copied to a variable for readability. - # Offers flexibility for dynamic viscosity via Sutherland's law where it depends - # on temperature and reference values, Ts and Tref such that mu(T) - mu = equations.mu + # In the simplest cases, the user passed in `mu` or `mu()` + # (which returns just a constant) but + # more complex functions like Sutherland's law are possible. + # `dynamic_viscosity` is a helper function that handles both cases + # by dispatching on the type of `equations.mu`. + mu = dynamic_viscosity(u, equations) if orientation == 1 # viscous flux components in the x-direction diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index e5567ae5789..9833122eb32 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -21,6 +21,9 @@ the [`CompressibleEulerEquations3D`](@ref). Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g., [``\mu``] = kg m⁻¹ s⁻¹. +The viscosity ``\mu`` may be a constant or a function of the current state, e.g., +depending on temperature (Sutherland's law): ``\mu = \mu(T)``. +In the latter case, the function `mu` needs to have the signature `mu(u, equations)`. The particular form of the compressible Navier-Stokes implemented is ```math @@ -80,7 +83,7 @@ where 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} ``` """ -struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, +struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, Mu, E <: AbstractCompressibleEulerEquations{3}} <: AbstractCompressibleNavierStokesDiffusion{3, 5, GradientVariables} # TODO: parabolic @@ -89,7 +92,7 @@ struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, gamma::RealT # ratio of specific heats inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - mu::RealT # viscosity + mu::Mu # viscosity Pr::RealT # Prandtl number kappa::RealT # thermal diffusivity for Fick's law @@ -103,16 +106,17 @@ function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquatio gradient_variables = GradientVariablesPrimitive()) gamma = equations.gamma inv_gamma_minus_one = equations.inv_gamma_minus_one - μ, Pr = promote(mu, Prandtl) # Under the assumption of constant Prandtl number the thermal conductivity - # constant is kappa = gamma μ / ((gamma-1) Pr). + # constant is kappa = gamma μ / ((gamma-1) Prandtl). # Important note! Factor of μ is accounted for later in `flux`. - kappa = gamma * inv_gamma_minus_one / Pr + # This avoids recomputation of kappa for non-constant μ. + kappa = gamma * inv_gamma_minus_one / Prandtl CompressibleNavierStokesDiffusion3D{typeof(gradient_variables), typeof(gamma), + typeof(mu), typeof(equations)}(gamma, inv_gamma_minus_one, - μ, Pr, kappa, + mu, Prandtl, kappa, equations, gradient_variables) end @@ -181,10 +185,12 @@ function flux(u, gradients, orientation::Integer, q2 = equations.kappa * dTdy q3 = equations.kappa * dTdz - # Constant dynamic viscosity is copied to a variable for readability. - # Offers flexibility for dynamic viscosity via Sutherland's law where it depends - # on temperature and reference values, Ts and Tref such that mu(T) - mu = equations.mu + # In the simplest cases, the user passed in `mu` or `mu()` + # (which returns just a constant) but + # more complex functions like Sutherland's law are possible. + # `dynamic_viscosity` is a helper function that handles both cases + # by dispatching on the type of `equations.mu`. + mu = dynamic_viscosity(u, equations) if orientation == 1 # viscous flux components in the x-direction diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 9f1382caa62..d47c34f9e75 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -476,20 +476,38 @@ end @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_shearlayer_amr.jl"), l2=[ - 0.00526017743452336, - 0.4130430692895672, - 0.4310996183791349, - 1.1544344171604635, + 0.005155557460409018, + 0.4048446934219344, + 0.43040068852937047, + 1.1255130552079322, ], linf=[ - 0.03492185879198495, - 1.392635891671335, - 1.357551616406459, - 8.713760873018146, + 0.03287305649809613, + 1.1656793717431393, + 1.3917196016246969, + 8.146587380114653, ], tspan=(0.0, 0.7)) end +@trixi_testset "TreeMesh2D: elixir_navierstokes_taylor_green_vortex_sutherland.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", + "elixir_navierstokes_taylor_green_vortex_sutherland.jl"), + l2=[ + 0.001452856280034929, + 0.0007538775539989481, + 0.0007538775539988681, + 0.011035506549989587, + ], + linf=[ + 0.003291912841311362, + 0.002986462478096974, + 0.0029864624780958637, + 0.0231954665514138, + ], + tspan=(0.0, 1.0)) +end + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"), diff --git a/test/test_unit.jl b/test/test_unit.jl index 03a78f6918a..79950f58d59 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1576,6 +1576,45 @@ end @test mesh.boundary_faces[:entire_boundary] == [1, 2] end + +@testset "Sutherlands Law" begin + function mu(u, equations) + T_ref = 291.15 + + R_specific_air = 287.052874 + T = R_specific_air * Trixi.temperature(u, equations) + + C_air = 120.0 + mu_ref_air = 1.827e-5 + + return mu_ref_air * (T_ref + C_air) / (T + C_air) * (T / T_ref)^1.5 + end + + function mu_control(u, equations, T_ref, R_specific, C, mu_ref) + T = R_specific * Trixi.temperature(u, equations) + + return mu_ref * (T_ref + C) / (T + C) * (T / T_ref)^1.5 + end + + # Dry air (values from Wikipedia: https://de.wikipedia.org/wiki/Sutherland-Modell) + T_ref = 291.15 + C = 120.0 # Sutherland's constant + R_specific = 287.052874 + mu_ref = 1.827e-5 + prandtl_number() = 0.72 + gamma = 1.4 + + equations = CompressibleEulerEquations2D(gamma) + equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu, + Prandtl = prandtl_number()) + + # Flow at rest + u = prim2cons(SVector(1.0, 0.0, 0.0, 1.0), equations_parabolic) + + # Comparison value from https://www.engineeringtoolbox.com/air-absolute-kinematic-viscosity-d_601.html at 18°C + @test isapprox(mu_control(u, equations_parabolic, T_ref, R_specific, C, mu_ref), + 1.803e-5, atol = 5e-8) +end end end #module From 57f4b2042208ca8a7c123b1719ad3a7944b59164 Mon Sep 17 00:00:00 2001 From: Arpit Babbar Date: Sat, 30 Mar 2024 13:32:37 +0530 Subject: [PATCH 06/13] Add AnalysisSurfaceIntegral (#1812) * Add AnalysisSurfaceIntegral * Minor change * Update elixir_euler_subsonic_cylinder.jl * Apply suggestions from code review Co-authored-by: Daniel Doehring * Fix labels, add tests * Attempt to fix CI * Obtain indices from user chosen function * Formatting * Minor change * Fix labels, add tests * Attempt to fix CI * Obtain indices from user chosen function * Add new elixirs * Add tests for NACA0012 * Fix tolerance, fix CI, add comments * Run formatter * Minor changes * Run formatter * Fix type instability and need for a vector * Fix typo! * Lower CFL to pass CI? * Reduce amr_interval to pass CI? * Maybe positivity limiter will fix CI? * Remove AMR from CI * That CI fail was on me! * Test for write permit * forces via names * Boundary Names for other examples * Shift to SSPRK54 to pass CI? * docstrings * alloc tests * fmt * Apply suggestions from code review Co-authored-by: Daniel Doehring Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Format * non-allocating BC * comments * comments * fmt * Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl * Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl * comments * make ready for review * typos * Apply suggestions from code review Some key points 1) Change pre to p 2) Don't capitalize U 3) Don't append Swanson quantities with 'sw' 4) linf -> l_inf (may be pending in some other places) 5) Other formatting improvements Co-authored-by: Andrew Winters * Some fixes * Format * Update NEWS.md * Update src/callbacks_step/analysis_surface_integral_2d.jl * Polish * Update examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl * increase cfl * fmt * fix --------- Co-authored-by: Daniel Doehring Co-authored-by: Daniel_Doehring Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Co-authored-by: Andrew Winters --- NEWS.md | 3 +- .../elixir_euler_NACA0012airfoil_mach085.jl | 132 ++++++++++++ .../elixir_euler_subsonic_cylinder.jl | 125 ++++++++++++ ...xir_navierstokes_NACA0012airfoil_mach08.jl | 155 ++++++++++++++ src/Trixi.jl | 3 +- src/callbacks_step/analysis.jl | 1 + .../analysis_surface_integral_2d.jl | 189 ++++++++++++++++++ src/solvers/dgsem_p4est/dg.jl | 3 +- .../sort_boundary_conditions.jl | 18 +- test/test_p4est_2d.jl | 80 ++++++++ test/test_parabolic_2d.jl | 22 ++ 11 files changed, 725 insertions(+), 6 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl create mode 100644 examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl create mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl create mode 100644 src/callbacks_step/analysis_surface_integral_2d.jl diff --git a/NEWS.md b/NEWS.md index 022252e61a9..5516fc9add5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,9 +7,10 @@ for human readability. ## Changes in the v0.7 lifecycle #### Added +- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh`. - +- New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients. ## Changes when updating to v0.7 from v0.6.x diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl new file mode 100644 index 00000000000..fb5f29bd038 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -0,0 +1,132 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +p_inf() = 1.0 +rho_inf() = p_inf() / (1.0 * 287.87) # p_inf = 1.0, T = 1, R = 287.87 +mach_inf() = 0.85 +aoa() = pi / 180.0 # 1 Degree angle of attack +c_inf(equations) = sqrt(equations.gamma * p_inf() / rho_inf()) +u_inf(equations) = mach_inf() * c_inf(equations) + +@inline function initial_condition_mach085_flow(x, t, + equations::CompressibleEulerEquations2D) + v1 = u_inf(equations) * cos(aoa()) + v2 = u_inf(equations) * sin(aoa()) + + prim = SVector(rho_inf(), v1, v2, p_inf()) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach085_flow + +volume_flux = flux_ranocha_turbo +surface_flux = flux_lax_friedrichs + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +shock_indicator = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp", + joinpath(@__DIR__, "NACA0012.inp")) + +mesh = P4estMesh{2}(mesh_file) + +# The outer boundary is constant but subsonic, so we cannot compute the +# boundary flux for the external information alone. Thus, we use the numerical flux to distinguish +# between inflow and outflow characteristics +@inline function boundary_condition_subsonic_constant(u_inner, + normal_direction::AbstractVector, x, + t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach085_flow(x, t, equations) + + return Trixi.flux_hll(u_inner, u_boundary, normal_direction, equations) +end + +boundary_conditions = Dict(:Left => boundary_condition_subsonic_constant, + :Right => boundary_condition_subsonic_constant, + :Top => boundary_condition_subsonic_constant, + :Bottom => boundary_condition_subsonic_constant, + :AirfoilBottom => boundary_condition_slip_wall, + :AirfoilTop => boundary_condition_slip_wall) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers + +# Run for a long time to reach a steady state +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 2000 + +l_inf = 1.0 # Length of airfoil + +force_boundary_names = [:AirfoilBottom, :AirfoilTop] +drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, + DragCoefficientPressure(aoa(), rho_inf(), + u_inf(equations), l_inf)) + +lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, + LiftCoefficientPressure(aoa(), rho_inf(), + u_inf(equations), l_inf)) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "out", + save_analysis = true, + analysis_integrals = (drag_coefficient, + lift_coefficient)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 1, + med_level = 3, med_threshold = 0.05, + max_level = 4, max_threshold = 0.1) + +amr_interval = 100 +amr_callback = AMRCallback(semi, amr_controller, + interval = amr_interval, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback, amr_callback) + +############################################################################### +# run the simulation +sol = solve(ode, SSPRK54(thread = OrdinaryDiffEq.True()), + 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/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl new file mode 100644 index 00000000000..dc23e192de8 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -0,0 +1,125 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +@inline function initial_condition_mach038_flow(x, t, + equations::CompressibleEulerEquations2D) + # set the freestream flow parameters + rho_freestream = 1.4 + v1 = 0.38 + v2 = 0.0 + p_freestream = 1.0 + + prim = SVector(rho_freestream, v1, v2, p_freestream) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach038_flow + +volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be used +surface_flux = flux_lax_friedrichs + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +function mapping2cylinder(xi, eta) + xi_, eta_ = 0.5 * (xi + 1), 0.5 * (eta + 1.0) # Map from [-1,1] to [0,1] for simplicity + + R2 = 50.0 # Bigger circle + R1 = 0.5 # Smaller circle + + # Ensure an isotropic mesh by using elements with smaller radial length near the inner circle + + r = R1 * exp(xi_ * log(R2 / R1)) + theta = 2.0 * pi * eta_ + + x = r * cos(theta) + y = r * sin(theta) + return (x, y) +end + +cells_per_dimension = (64, 64) +# xi = -1 maps to the inner circle and xi = +1 maps to the outer circle and we can specify boundary +# conditions there. However, the image of eta = -1, +1 coincides at the line y = 0. There is no +# physical boundary there so we specify `periodicity = true` there and the solver treats the +# element across eta = -1, +1 as neighbours which is what we want +mesh = P4estMesh(cells_per_dimension, mapping = mapping2cylinder, polydeg = polydeg, + periodicity = (false, true)) + +# The boundary condition on the outer cylinder is constant but subsonic, so we cannot compute the +# boundary flux from the external information alone. Thus, we use the numerical flux to distinguish +# between inflow and outflow characteristics +@inline function boundary_condition_subsonic_constant(u_inner, + normal_direction::AbstractVector, x, + t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach038_flow(x, t, equations) + + return surface_flux_function(u_inner, u_boundary, normal_direction, equations) +end + +boundary_conditions = Dict(:x_neg => boundary_condition_slip_wall, + :x_pos => boundary_condition_subsonic_constant) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers + +# Run for a long time to reach a steady state +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 2000 + +aoa = 0.0 +rho_inf = 1.4 +u_inf = 0.38 +l_inf = 1.0 # Diameter of circle + +drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, + DragCoefficientPressure(aoa, rho_inf, u_inf, + l_inf)) + +lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, + LiftCoefficientPressure(aoa, rho_inf, u_inf, + l_inf)) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "out", + save_analysis = true, + analysis_integrals = (drag_coefficient, + lift_coefficient)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 2.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, + CarpenterKennedy2N54(williamson_condition = false; + thread = OrdinaryDiffEq.True()), + 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/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl new file mode 100644 index 00000000000..5bfa21c5322 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -0,0 +1,155 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +# Laminar transonic flow around a NACA0012 airfoil. + +# This test is taken from the paper of Swanson and Langer. The values for the drag and lift coefficients +# from Case 5 in Table 3 are used to validate the scheme and computation of surface forces. + +# References: +# - Roy Charles Swanson, Stefan Langer (2016) +# Structured and Unstructured Grid Methods (2016) +# [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623) +# - Deep Ray, Praveen Chandrashekar (2017) +# An entropy stable finite volume scheme for the +# two dimensional Navier–Stokes equations on triangular grids +# [DOI:10.1016/j.amc.2017.07.020](https://doi.org/10.1016/j.amc.2017.07.020) + +equations = CompressibleEulerEquations2D(1.4) + +prandtl_number() = 0.72 +mu() = 0.0031959974968701088 +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), + Prandtl = prandtl_number()) + +rho_inf() = 1.0 +p_inf() = 2.85 +aoa() = 10.0 * pi / 180.0 # 10 degree angle of attack +l_inf() = 1.0 +mach_inf() = 0.8 +u_inf(equations) = mach_inf() * sqrt(equations.gamma * p_inf() / rho_inf()) +@inline function initial_condition_mach08_flow(x, t, equations) + # set the freestream flow parameters + gamma = equations.gamma + u_inf = mach_inf() * sqrt(gamma * p_inf() / rho_inf()) + + v1 = u_inf * cos(aoa()) + v2 = u_inf * sin(aoa()) + + prim = SVector(rho_inf(), v1, v2, p_inf()) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach08_flow + +surface_flux = flux_lax_friedrichs + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux) + +mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp", + joinpath(@__DIR__, "NACA0012.inp")) + +mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 1) + +# The boundary values across outer boundary are constant but subsonic, so we cannot compute the +# boundary flux from the external information alone. Thus, we use the numerical flux to distinguish +# between inflow and outflow characteristics +@inline function boundary_condition_subsonic_constant(u_inner, + normal_direction::AbstractVector, x, + t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach08_flow(x, t, equations) + + return Trixi.flux_hll(u_inner, u_boundary, normal_direction, equations) +end + +boundary_conditions = Dict(:Left => boundary_condition_subsonic_constant, + :Right => boundary_condition_subsonic_constant, + :Top => boundary_condition_subsonic_constant, + :Bottom => boundary_condition_subsonic_constant, + :AirfoilBottom => boundary_condition_slip_wall, + :AirfoilTop => boundary_condition_slip_wall) + +velocity_airfoil = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) + +heat_airfoil = Adiabatic((x, t, equations) -> 0.0) + +boundary_conditions_airfoil = BoundaryConditionNavierStokesWall(velocity_airfoil, + heat_airfoil) + +function momenta_initial_condition_mach08_flow(x, t, equations) + u = initial_condition_mach08_flow(x, t, equations) + momenta = SVector(u[2], u[3]) +end +velocity_bc_square = NoSlip((x, t, equations) -> momenta_initial_condition_mach08_flow(x, t, + equations)) + +heat_bc_square = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_square = BoundaryConditionNavierStokesWall(velocity_bc_square, + heat_bc_square) + +boundary_conditions_parabolic = Dict(:Left => boundary_condition_square, + :Right => boundary_condition_square, + :Top => boundary_condition_square, + :Bottom => boundary_condition_square, + :AirfoilBottom => boundary_conditions_airfoil, + :AirfoilTop => boundary_conditions_airfoil) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers + +# Run for a long time to reach a state where forces stabilize up to 3 digits +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 2000 + +force_boundary_names = [:AirfoilBottom, :AirfoilTop] +drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, + DragCoefficientPressure(aoa(), rho_inf(), + u_inf(equations), + l_inf())) + +lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, + LiftCoefficientPressure(aoa(), rho_inf(), + u_inf(equations), + l_inf())) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "out", + save_analysis = true, + analysis_errors = Symbol[], + analysis_integrals = (drag_coefficient, + lift_coefficient)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1e-8, + reltol = 1e-8, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index 9375c80d77e..168441513ad 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -261,7 +261,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, AveragingCallback, AMRCallback, StepsizeCallback, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, - TrivialCallback, AnalysisCallbackCoupled + TrivialCallback, AnalysisCallbackCoupled, + AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure export load_mesh, load_time, load_timestep, load_timestep!, load_dt, load_adaptive_time_integrator! diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index ba232032951..62048853b53 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -691,6 +691,7 @@ end # @muladd # specialized implementations specific to some solvers include("analysis_dg1d.jl") include("analysis_dg2d.jl") +include("analysis_surface_integral_2d.jl") include("analysis_dg2d_parallel.jl") include("analysis_dg3d.jl") include("analysis_dg3d_parallel.jl") diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl new file mode 100644 index 00000000000..902fb0b4e85 --- /dev/null +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -0,0 +1,189 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# This file contains callbacks that are performed on the surface like computation of +# surface forces + +""" + AnalysisSurfaceIntegral{Semidiscretization, Variable}(semi, + boundary_symbol_or_boundary_symbols, + variable) + + This struct is used to compute the surface integral of a quantity of interest `variable` alongside + the boundary/boundaries associated with particular name(s) given in `boundary_symbol` + or `boundary_symbols`. + For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or + drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary + name `:Airfoil` in 2D. +""" +struct AnalysisSurfaceIntegral{Semidiscretization, Variable} + semi::Semidiscretization # passed in to retrieve boundary condition information + indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed + variable::Variable # Quantity of interest, like lift or drag + + function AnalysisSurfaceIntegral(semi, boundary_symbol, variable) + @unpack boundary_symbol_indices = semi.boundary_conditions + indices = boundary_symbol_indices[boundary_symbol] + + return new{typeof(semi), typeof(variable)}(semi, indices, + variable) + end + + function AnalysisSurfaceIntegral(semi, boundary_symbols::Vector{Symbol}, variable) + @unpack boundary_symbol_indices = semi.boundary_conditions + indices = Vector{Int}() + for name in boundary_symbols + append!(indices, boundary_symbol_indices[name]) + end + sort!(indices) + + return new{typeof(semi), typeof(variable)}(semi, indices, + variable) + end +end + +struct ForceState{RealT <: Real} + psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream + rhoinf::RealT + uinf::RealT + l_inf::RealT +end + +struct LiftCoefficientPressure{RealT <: Real} + force_state::ForceState{RealT} +end + +struct DragCoefficientPressure{RealT <: Real} + force_state::ForceState{RealT} +end + +""" + LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) + +Compute the lift coefficient +```math +C_{L,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_L \\, \\mathrm{d} S} + {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} +``` +based on the pressure distribution along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) + # psi_lift is the normal unit vector to the freestream direction. + # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa)) + # leads to positive lift coefficients for positive angles of attack for airfoils. + # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same + # value, but with the opposite sign. + psi_lift = (-sin(aoa), cos(aoa)) + return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, l_inf)) +end + +""" + DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) + +Compute the drag coefficient +```math +C_{D,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_D \\, \\mathrm{d} S} + {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} +``` +based on the pressure distribution along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) + # `psi_drag` is the unit vector tangent to the freestream direction + psi_drag = (cos(aoa), sin(aoa)) + return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, l_inf)) +end + +function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations) + p = pressure(u, equations) + @unpack psi, rhoinf, uinf, l_inf = lift_coefficient.force_state + # Normalize as `normal_direction` is not necessarily a unit vector + n = dot(normal_direction, psi) / norm(normal_direction) + return p * n / (0.5 * rhoinf * uinf^2 * l_inf) +end + +function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations) + p = pressure(u, equations) + @unpack psi, rhoinf, uinf, l_inf = drag_coefficient.force_state + # Normalize as `normal_direction` is not necessarily a unit vector + n = dot(normal_direction, psi) / norm(normal_direction) + return p * n / (0.5 * rhoinf * uinf^2 * l_inf) +end + +function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, + mesh::P4estMesh{2}, + equations, dg::DGSEM, cache) + @unpack boundaries = cache + @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements + @unpack weights = dg.basis + + @unpack indices, variable = surface_variable + + surface_integral = zero(eltype(u)) + index_range = eachnode(dg) + for boundary in indices + element = boundaries.neighbor_ids[boundary] + node_indices = boundaries.node_indices[boundary] + direction = indices2direction(node_indices) + + i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) + j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + + i_node = i_node_start + j_node = j_node_start + for node_index in index_range + u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index, + boundary) + # Extract normal direction at nodes which points from the elements outwards, + # i.e., *into* the structure. + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, + element) + + # L2 norm of normal direction (contravariant_vector) is the surface element + dS = weights[node_index] * norm(normal_direction) + # Integral over entire boundary surface + surface_integral += variable(u_node, normal_direction, equations) * dS + + i_node += i_node_step + j_node += j_node_step + end + end + return surface_integral +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, + <:LiftCoefficientPressure{<:Any}}) + "CL_p" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, + <:LiftCoefficientPressure{<:Any}}) + "CL_p" +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, + <:DragCoefficientPressure{<:Any}}) + "CD_p" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, + <:DragCoefficientPressure{<:Any}}) + "CD_p" +end +end # muladd diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index ec50627d3ef..10cc075089c 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -36,7 +36,8 @@ end orientation = (direction + 1) >> 1 normal = get_contravariant_vector(orientation, contravariant_vectors, indices...) - # Contravariant vectors at interfaces in negative coordinate direction are pointing inwards + # Contravariant vectors at interfaces in negative coordinate direction are pointing inwards, + # flip sign to make them point outwards if isodd(direction) return -normal else diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index b5388cadc8b..2c2c6876d70 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -8,8 +8,8 @@ """ UnstructuredSortedBoundaryTypes -General container to sort the boundary conditions by type for some unstructured meshes/solvers. -It stores a set of global indices for each boundary condition type to expedite computation +General container to sort the boundary conditions by type and name for some unstructured meshes/solvers. +It stores a set of global indices for each boundary condition type and name to expedite computation during the call to `calc_boundary_flux!`. The original dictionary form of the boundary conditions set by the user in the elixir file is also stored for printing. """ @@ -17,6 +17,7 @@ mutable struct UnstructuredSortedBoundaryTypes{N, BCs <: NTuple{N, Any}} boundary_condition_types::BCs # specific boundary condition type(s), e.g. BoundaryConditionDirichlet boundary_indices::NTuple{N, Vector{Int}} # integer vectors containing global boundary indices boundary_dictionary::Dict{Symbol, Any} # boundary conditions as set by the user in the elixir file + boundary_symbol_indices::Dict{Symbol, Vector{Int}} # integer vectors containing global boundary indices per boundary identifier end # constructor that "eats" the original boundary condition dictionary and sorts the information @@ -28,10 +29,14 @@ function UnstructuredSortedBoundaryTypes(boundary_conditions::Dict, cache) n_boundary_types = length(boundary_condition_types) boundary_indices = ntuple(_ -> [], n_boundary_types) + # Initialize `boundary_symbol_indices` as an empty dictionary, filled later in `initialize!` + boundary_symbol_indices = Dict{Symbol, Vector{Int}}() + container = UnstructuredSortedBoundaryTypes{n_boundary_types, typeof(boundary_condition_types)}(boundary_condition_types, boundary_indices, - boundary_conditions) + boundary_conditions, + boundary_symbol_indices) initialize!(container, cache) end @@ -97,6 +102,13 @@ function initialize!(boundary_types_container::UnstructuredSortedBoundaryTypes{N # convert the work array with the boundary indices into a tuple boundary_types_container.boundary_indices = Tuple(_boundary_indices) + # Store boundary indices per symbol (required for force computations, for instance) + for (symbol, _) in boundary_dictionary + indices = findall(x -> x === symbol, cache.boundaries.name) + # Store the indices in `boundary_symbol_indices` dictionary + boundary_types_container.boundary_symbol_indices[symbol] = sort!(indices) + end + return boundary_types_container end end # @muladd diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 121001b35ff..1bacccbf812 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -531,6 +531,86 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_subsonic_cylinder.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_subsonic_cylinder.jl"), + l2=[ + 0.00011914390523852561, + 0.00010776028621724485, + 6.139954358305467e-5, + 0.0003067693731825959, + ], + linf=[ + 0.1653075586200805, + 0.1868437275544909, + 0.09772818519679008, + 0.4311796171737692, + ], tspan=(0.0, 0.001)) + # 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 + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + + @test isapprox(lift, -6.501138753497174e-15, atol = 1e-13) + @test isapprox(drag, 2.588589856781827, atol = 1e-13) +end + +# Forces computation test in an AMR code +@trixi_testset "elixir_euler_NACA0012airfoil_mach085.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_NACA0012airfoil_mach085.jl"), + l2=[ + 5.371568111383228e-7, 6.4158131303956445e-6, + 1.0324346542348325e-5, 0.0006348064933187732, + ], + linf=[ + 0.0016263400091978443, 0.028471072159724428, + 0.02986133204785877, 1.9481060511014872, + ], + base_level=0, med_level=1, max_level=1, + tspan=(0.0, 0.0001), + adapt_initial_condition=false, + adapt_initial_condition_only_refine=false, + # With the default `maxiters = 1` in coverage tests, + # the values for `drag` and `lift` below would differ. + coverage_override=(maxiters = 100_000,)) + + # 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 + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + + @test isapprox(lift, 0.0262382560809345, atol = 1e-13) + @test isapprox(drag, 0.10898248971932244, atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index d47c34f9e75..edbdfe5a84f 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -707,6 +707,28 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_navierstokes_NACA0012airfoil_mach08.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_navierstokes_NACA0012airfoil_mach08.jl"), + l2=[0.000186486564226516, + 0.0005076712323400374, + 0.00038074588984354107, + 0.002128177239782089], + linf=[0.5153387072802718, + 1.199362305026636, + 0.9077214424040279, + 5.666071182328691], tspan=(0.0, 0.001), + initial_refinement_level=0) + # 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 # Clean up afterwards: delete Trixi.jl output directory From 641fde8c40e5b182193320fb805def4266d894d8 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 2 Apr 2024 09:31:47 +0200 Subject: [PATCH 07/13] Bump crate-ci/typos from 1.18.2 to 1.19.0 (#1894) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.18.2 to 1.19.0. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.18.2...v1.19.0) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 87e34cb50f3..e8f2271f831 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.18.2 + uses: crate-ci/typos@v1.19.0 From 09600b56915296e28e5792f87c91a820e9c4eee4 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 3 Apr 2024 08:15:22 +0200 Subject: [PATCH 08/13] Bump julia-actions/setup-julia from 1 to 2 (#1895) Bumps [julia-actions/setup-julia](https://github.com/julia-actions/setup-julia) from 1 to 2. - [Release notes](https://github.com/julia-actions/setup-julia/releases) - [Commits](https://github.com/julia-actions/setup-julia/compare/v1...v2) --- updated-dependencies: - dependency-name: julia-actions/setup-julia dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Documenter.yml | 2 +- .github/workflows/Downgrade.yml | 2 +- .github/workflows/Invalidations.yml | 2 +- .github/workflows/benchmark.yml | 2 +- .github/workflows/ci.yml | 2 +- .github/workflows/downstream.yml | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 8f8d674ffa9..1ac862d2c85 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -34,7 +34,7 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: '1.9' show-versioninfo: true diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index dd5d8ee7e32..234939a8017 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -64,7 +64,7 @@ jobs: - threaded steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml index 18048d26be8..b2d34cbc856 100644 --- a/.github/workflows/Invalidations.yml +++ b/.github/workflows/Invalidations.yml @@ -16,7 +16,7 @@ jobs: if: github.base_ref == github.event.repository.default_branch runs-on: ubuntu-latest steps: - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: '1' - uses: actions/checkout@v4 diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 4531c3aee0a..15575717f11 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -21,7 +21,7 @@ jobs: - run: | git fetch --tags git branch --create-reflog main origin/main - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a7dfe033a90..27f7af7c528 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -107,7 +107,7 @@ jobs: trixi_test: threaded steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml index b40d5d365c0..83ed6cc79c7 100644 --- a/.github/workflows/downstream.yml +++ b/.github/workflows/downstream.yml @@ -64,7 +64,7 @@ jobs: - TrixiShallowWater.jl steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} From b3bb45745db5803c65ea151a17690b3b4b42dc57 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 3 Apr 2024 08:33:04 +0200 Subject: [PATCH 09/13] 1D Linearized Euler (#1867) * 1D Linearized Euler * remove redundant comment * Fix docu * Update src/equations/linearized_euler_1d.jl * Update examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl * code review * fmt * Apply suggestions from code review Co-authored-by: Andrew Winters * Update src/equations/linearized_euler_1d.jl * Comments * comments * news --------- Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha --- NEWS.md | 1 + ...r_linearizedeuler_characteristic_system.jl | 112 ++++++++++++++ .../elixir_linearizedeuler_convergence.jl | 59 +++++++ .../elixir_linearizedeuler_gauss_wall.jl | 65 ++++++++ src/Trixi.jl | 2 +- src/equations/equations.jl | 1 + src/equations/linearized_euler_1d.jl | 144 ++++++++++++++++++ test/test_structured_1d.jl | 15 ++ test/test_tree_1d.jl | 3 + test/test_tree_1d_linearizedeuler.jl | 51 +++++++ 10 files changed, 452 insertions(+), 1 deletion(-) create mode 100644 examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl create mode 100644 examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl create mode 100644 examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl create mode 100644 src/equations/linearized_euler_1d.jl create mode 100644 test/test_tree_1d_linearizedeuler.jl diff --git a/NEWS.md b/NEWS.md index 5516fc9add5..e1238562fc9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,7 @@ for human readability. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh`. +- Implementation of 1D Linearized Euler Equations. - New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients. ## Changes when updating to v0.7 from v0.6.x diff --git a/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl new file mode 100644 index 00000000000..663b25b18c0 --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_linearizedeuler_characteristic_system.jl @@ -0,0 +1,112 @@ + +using OrdinaryDiffEq +using LinearAlgebra: dot +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +rho_0 = 1.0 +v_0 = 1.0 +c_0 = 1.0 +equations = LinearizedEulerEquations1D(rho_0, v_0, c_0) + +solver = DGSEM(polydeg = 3, surface_flux = flux_hll) + +coordinates_min = (0.0,) # minimum coordinate +coordinates_max = (1.0,) # maximum coordinate +cells_per_dimension = (64,) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +# For eigensystem of the linearized Euler equations see e.g. +# https://www.nas.nasa.gov/assets/nas/pdf/ams/2018/introtocfd/Intro2CFD_Lecture1_Pulliam_Euler_WaveEQ.pdf +# Linearized Euler: Eigensystem +lin_euler_eigvals = [v_0 - c_0; v_0; v_0 + c_0] +lin_euler_eigvecs = [-rho_0/c_0 1 rho_0/c_0; + 1 0 1; + -rho_0*c_0 0 rho_0*c_0] +lin_euler_eigvecs_inv = inv(lin_euler_eigvecs) + +# Trace back characteristics. +# See https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf, p.95 +function compute_char_initial_pos(x, t) + return SVector(x[1], x[1], x[1]) .- t * lin_euler_eigvals +end + +function compute_primal_sol(char_vars) + return lin_euler_eigvecs * char_vars +end + +# Initial condition is in principle arbitrary, only periodicity is required +function initial_condition_entropy_wave(x, t, equations::LinearizedEulerEquations1D) + # Parameters + alpha = 1.0 + beta = 150.0 + center = 0.5 + + rho_prime = alpha * exp(-beta * (x[1] - center)^2) + v_prime = 0.0 + p_prime = 0.0 + + return SVector(rho_prime, v_prime, p_prime) +end + +function initial_condition_char_vars(x, t, equations::LinearizedEulerEquations1D) + # Trace back characteristics + x_char = compute_char_initial_pos(x, t) + + # Employ periodicity + for p in 1:3 + while x_char[p] < coordinates_min[1] + x_char[p] += coordinates_max[1] - coordinates_min[1] + end + while x_char[p] > coordinates_max[1] + x_char[p] -= coordinates_max[1] - coordinates_min[1] + end + end + + # Set up characteristic variables + w = zeros(3) + t_0 = 0 # Assumes t_0 = 0 + for p in 1:3 + u_char = initial_condition_entropy_wave(x_char[p], t_0, equations) + w[p] = dot(lin_euler_eigvecs_inv[p, :], u_char) + end + + return compute_primal_sol(w) +end + +initial_condition = initial_condition_char_vars + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.3) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_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_1d_dgsem/elixir_linearizedeuler_convergence.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl new file mode 100644 index 00000000000..5b17ab4f3dc --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_convergence.jl @@ -0,0 +1,59 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations1D(v_mean_global = 0.0, c_mean_global = 1.0, + rho_mean_global = 1.0) + +initial_condition = initial_condition_convergence_test + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-1.0) +coordinates_max = (1.0) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +analysis_interval = 100 + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.8) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +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_linearizedeuler_gauss_wall.jl b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl new file mode 100644 index 00000000000..0884249559a --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_linearizedeuler_gauss_wall.jl @@ -0,0 +1,65 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linearized Euler equations + +equations = LinearizedEulerEquations1D(v_mean_global = 0.5, c_mean_global = 1.0, + rho_mean_global = 1.0) + +solver = DGSEM(polydeg = 5, surface_flux = flux_hll) + +coordinates_min = (0.0,) +coordinates_max = (90.0,) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 100_000, + periodicity = false) + +# Initialize density and pressure perturbation with a Gaussian bump +# that is advected to left with v - c and to the right with v + c. +# Correspondigly, the bump splits in half. +function initial_condition_gauss_wall(x, t, equations::LinearizedEulerEquations1D) + v1_prime = 0.0 + rho_prime = p_prime = 2 * exp(-(x[1] - 45)^2 / 25) + return SVector(rho_prime, v1_prime, p_prime) +end +initial_condition = initial_condition_gauss_wall + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition_wall) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 30.0 +tspan = (0.0, 30.0) +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.7) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +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) + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index 168441513ad..883f8d66f07 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -160,7 +160,7 @@ export AcousticPerturbationEquations2D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D, ShallowWaterEquationsQuasi1D, - LinearizedEulerEquations2D, + LinearizedEulerEquations1D, LinearizedEulerEquations2D, PolytropicEulerEquations2D, TrafficFlowLWREquations1D diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8f476cf6f16..ef4dac8b1a5 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -501,6 +501,7 @@ include("acoustic_perturbation_2d.jl") # Linearized Euler equations abstract type AbstractLinearizedEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end +include("linearized_euler_1d.jl") include("linearized_euler_2d.jl") abstract type AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} <: diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl new file mode 100644 index 00000000000..22f452219b5 --- /dev/null +++ b/src/equations/linearized_euler_1d.jl @@ -0,0 +1,144 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + LinearizedEulerEquations1D(v_mean_global, c_mean_global, rho_mean_global) + +Linearized euler equations in one space dimension. The equations are given by +```math +\partial_t +\begin{pmatrix} + \rho' \\ v_1' \\ p' +\end{pmatrix} ++ +\partial_x +\begin{pmatrix} + \bar{\rho} v_1' + \bar{v_1} \rho ' \\ \bar{v_1} v_1' + \frac{p'}{\bar{\rho}} \\ \bar{v_1} p' + c^2 \bar{\rho} v_1' +\end{pmatrix} += +\begin{pmatrix} + 0 \\ 0 \\ 0 +\end{pmatrix} +``` +The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. +The unknowns are the perturbation quantities of the acoustic velocity ``v_1'``, the pressure ``p'`` +and the density ``\rho'``. +""" +struct LinearizedEulerEquations1D{RealT <: Real} <: + AbstractLinearizedEulerEquations{1, 3} + v_mean_global::RealT + c_mean_global::RealT + rho_mean_global::RealT +end + +function LinearizedEulerEquations1D(v_mean_global::Real, + c_mean_global::Real, rho_mean_global::Real) + if rho_mean_global < 0 + throw(ArgumentError("rho_mean_global must be non-negative")) + elseif c_mean_global < 0 + throw(ArgumentError("c_mean_global must be non-negative")) + end + + return LinearizedEulerEquations1D(v_mean_global, c_mean_global, + rho_mean_global) +end + +# Constructor with keywords +function LinearizedEulerEquations1D(; v_mean_global::Real, + c_mean_global::Real, rho_mean_global::Real) + return LinearizedEulerEquations1D(v_mean_global, c_mean_global, + rho_mean_global) +end + +function varnames(::typeof(cons2cons), ::LinearizedEulerEquations1D) + ("rho_prime", "v1_prime", "p_prime") +end +function varnames(::typeof(cons2prim), ::LinearizedEulerEquations1D) + ("rho_prime", "v1_prime", "p_prime") +end + +""" + initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations1D) + +A smooth initial condition used for convergence tests. +""" +function initial_condition_convergence_test(x, t, equations::LinearizedEulerEquations1D) + rho_prime = -cospi(2 * t) * sinpi(2 * x[1]) + v1_prime = sinpi(2 * t) * cospi(2 * x[1]) + p_prime = rho_prime + + return SVector(rho_prime, v1_prime, p_prime) +end + +""" + boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function, + equations::LinearizedEulerEquations1D) + +Boundary conditions for a solid wall. +""" +function boundary_condition_wall(u_inner, orientation, direction, x, t, + surface_flux_function, + equations::LinearizedEulerEquations1D) + # Boundary state is equal to the inner state except for the velocity. For boundaries + # in the -x/+x direction, we multiply the velocity (in the x direction by) -1. + u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3]) + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + + return flux +end + +# Calculate 1D flux for a single point +@inline function flux(u, orientation::Integer, equations::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global, rho_mean_global = equations + rho_prime, v1_prime, p_prime = u + f1 = v_mean_global * rho_prime + rho_mean_global * v1_prime + f2 = v_mean_global * v1_prime + p_prime / rho_mean_global + f3 = v_mean_global * p_prime + c_mean_global^2 * rho_mean_global * v1_prime + + return SVector(f1, f2, f3) +end + +@inline have_constant_speed(::LinearizedEulerEquations1D) = True() + +@inline function max_abs_speeds(equations::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global = equations + return abs(v_mean_global) + c_mean_global +end + +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global = equations + return abs(v_mean_global) + c_mean_global +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::LinearizedEulerEquations1D) + min_max_speed_davis(u_ll, u_rr, orientation, equations) +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::LinearizedEulerEquations1D) + @unpack v_mean_global, c_mean_global = equations + + λ_min = v_mean_global - c_mean_global + λ_max = v_mean_global + c_mean_global + + return λ_min, λ_max +end + +# Convert conservative variables to primitive +@inline cons2prim(u, equations::LinearizedEulerEquations1D) = u +@inline cons2entropy(u, ::LinearizedEulerEquations1D) = u +end # muladd diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index fea06554c57..f97696d089a 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -139,6 +139,21 @@ end end end +@trixi_testset "elixir_linearizedeuler_characteristic_system.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_linearizedeuler_characteristic_system.jl"), + l2=[2.9318078842789714e-6, 0.0, 0.0], + linf=[4.291208715723194e-5, 0.0, 0.0]) + # 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_traffic_flow_lwr_greenlight.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_traffic_flow_lwr_greenlight.jl"), diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 4a25a51a45e..76129f15e07 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -48,6 +48,9 @@ isdir(outdir) && rm(outdir, recursive = true) # Traffic flow LWR include("test_tree_1d_traffic_flow_lwr.jl") + + # Linearized Euler + include("test_tree_1d_linearizedeuler.jl") end # Coverage test for all initial conditions diff --git a/test/test_tree_1d_linearizedeuler.jl b/test/test_tree_1d_linearizedeuler.jl new file mode 100644 index 00000000000..c7cffee3f66 --- /dev/null +++ b/test/test_tree_1d_linearizedeuler.jl @@ -0,0 +1,51 @@ + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") + +@testset "Linearized Euler Equations 1D" begin +#! format: noindent + +@trixi_testset "elixir_linearizedeuler_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_convergence.jl"), + l2=[ + 0.00010894927270421941, + 0.00014295255695912358, + 0.00010894927270421941, + ], + linf=[ + 0.0005154647164193893, + 0.00048457837684242266, + 0.0005154647164193893, + ]) + # 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_linearizedeuler_gauss_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_gauss_wall.jl"), + l2=[0.650082087850354, 0.2913911415488769, 0.650082087850354], + linf=[ + 1.9999505145390108, + 0.9999720404625275, + 1.9999505145390108, + ]) + # 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 From 876bdeb130099d4d218c817805432f24af1c2285 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 3 Apr 2024 08:50:13 +0200 Subject: [PATCH 10/13] remove duplicate entry from NEWS (#1896) --- NEWS.md | 1 - 1 file changed, 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index e1238562fc9..66dd1e58dc6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,7 +7,6 @@ for human readability. ## Changes in the v0.7 lifecycle #### Added -- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh`. - Implementation of 1D Linearized Euler Equations. From 01153c503ebc85dd0f5685ed139844af057e7cb5 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 3 Apr 2024 08:51:19 +0200 Subject: [PATCH 11/13] set version to v0.7.6 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6ff7f29686d..60f340d1306 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.7.6-pre" +version = "0.7.6" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 90d6d6588cca6a968125c0d1e91a758f566b98a7 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 3 Apr 2024 08:51:32 +0200 Subject: [PATCH 12/13] set development version to v0.7.7-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 60f340d1306..89052859e42 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.7.6" +version = "0.7.7-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 0f60ebc55189260ec108a07fdbf85d000a7c92d1 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 5 Apr 2024 06:27:08 +0200 Subject: [PATCH 13/13] Math rendering speed of sound docstring linearized Euler Eqs (#1897) * Update linearized_euler_1d.jl * Update linearized_euler_2d.jl --- src/equations/linearized_euler_1d.jl | 2 +- src/equations/linearized_euler_2d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl index 22f452219b5..856a53a3d5d 100644 --- a/src/equations/linearized_euler_1d.jl +++ b/src/equations/linearized_euler_1d.jl @@ -24,7 +24,7 @@ Linearized euler equations in one space dimension. The equations are given by 0 \\ 0 \\ 0 \end{pmatrix} ``` -The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. +The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and ``c`` is the speed of sound. The unknowns are the perturbation quantities of the acoustic velocity ``v_1'``, the pressure ``p'`` and the density ``\rho'``. """ diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl index d497762bf62..d9d0d44e553 100644 --- a/src/equations/linearized_euler_2d.jl +++ b/src/equations/linearized_euler_2d.jl @@ -29,7 +29,7 @@ Linearized euler equations in two space dimensions. The equations are given by 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} ``` -The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and c is the speed of sound. +The bar ``\bar{(\cdot)}`` indicates uniform mean flow variables and ``c`` is the speed of sound. The unknowns are the acoustic velocities ``v' = (v_1', v_2')``, the pressure ``p'`` and the density ``\rho'``. """ struct LinearizedEulerEquations2D{RealT <: Real} <: