diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 129c41a3b5c..8f8d674ffa9 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -38,6 +38,7 @@ jobs: with: version: '1.9' show-versioninfo: true + - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@v1 env: PYTHON: "" diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml index 2ea30d6fddb..6aa4809c1c2 100644 --- a/.github/workflows/benchmark.yml +++ b/.github/workflows/benchmark.yml @@ -45,7 +45,7 @@ jobs: run: julia --project=benchmark/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Run benchmarks run: julia --project=benchmark/ --color=yes benchmark/run_benchmarks.jl - - uses: actions/upload-artifact@v3 + - uses: actions/upload-artifact@v4 with: name: my-artifact path: benchmark/results*.md diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index cf8107736e9..f287cc5feb2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -70,6 +70,7 @@ jobs: - p4est_part1 - p4est_part2 - t8code_part1 + - t8code_part2 - unstructured_dgmulti - parabolic - paper_self_gravitating_gas_dynamics @@ -153,7 +154,7 @@ jobs: - shell: bash run: | cp ./lcov.info ./lcov-${{ matrix.trixi_test }}-${{ matrix.os }}-${{ matrix.version }}-${{ matrix.arch }}.info - - uses: actions/upload-artifact@v3 + - uses: actions/upload-artifact@v4 with: name: lcov-${{ matrix.trixi_test }}-${{ matrix.os }}-${{ matrix.version }}-${{ matrix.arch }} path: ./lcov-${{ matrix.trixi_test }}-${{ matrix.os }}-${{ matrix.version }}-${{ matrix.arch }}.info @@ -176,7 +177,7 @@ jobs: # At first, we check out the repository and download all artifacts # (and list files for debugging). - uses: actions/checkout@v4 - - uses: actions/download-artifact@v3 + - uses: actions/download-artifact@v4 - run: ls -R # Next, we merge the individual coverage files and upload # the combined results to Coveralls. @@ -199,7 +200,7 @@ jobs: github-token: ${{ secrets.GITHUB_TOKEN }} path-to-lcov: ./lcov.info # Upload merged coverage data as artifact for debugging - - uses: actions/upload-artifact@v3 + - uses: actions/upload-artifact@v4 with: name: lcov path: ./lcov.info diff --git a/.gitignore b/.gitignore index 3132b9af38b..b4f1cf6bb47 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ *.mesh *.bson *.inp +*.msh **/Manifest.toml out*/ docs/build diff --git a/Project.toml b/Project.toml index faf9f82d335..f246bdfdab4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.6.6-pre" +version = "0.6.7-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/docs/literate/make.jl b/docs/literate/make.jl index a04d8a0b333..84e4fbdced6 100644 --- a/docs/literate/make.jl +++ b/docs/literate/make.jl @@ -10,17 +10,17 @@ function create_files(title, file, repo_src, pages_dir, notebooks_dir; folder="" end binder_logo = "https://mybinder.org/badge_logo.svg" - nbviewer_logo = "https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg" - download_logo = "https://camo.githubusercontent.com/aea75103f6d9f690a19cb0e17c06f984ab0f472d9e6fe4eadaa0cc438ba88ada/68747470733a2f2f696d672e736869656c64732e696f2f62616467652f646f776e6c6f61642d6e6f7465626f6f6b2d627269676874677265656e" + nbviewer_logo = "https://img.shields.io/badge/render-nbviewer-f37726" + raw_notebook_logo = "https://img.shields.io/badge/raw-notebook-4cc61e" notebook_path = "tutorials/notebooks/$notebook_filename" binder_url = "https://mybinder.org/v2/gh/trixi-framework/Trixi.jl/tutorial_notebooks?filepath=$notebook_path" nbviewer_url = "https://nbviewer.jupyter.org/github/trixi-framework/Trixi.jl/blob/tutorial_notebooks/$notebook_path" - download_url = "https://raw.githubusercontent.com/trixi-framework/Trixi.jl/tutorial_notebooks/$notebook_path" + raw_notebook_url = "https://raw.githubusercontent.com/trixi-framework/Trixi.jl/tutorial_notebooks/$notebook_path" binder_badge = "# [![]($binder_logo)]($binder_url)" nbviewer_badge = "# [![]($nbviewer_logo)]($nbviewer_url)" - download_badge = "# [![]($download_logo)]($download_url)" + raw_notebook_badge = "# [![]($raw_notebook_logo)]($raw_notebook_url)" # Generate notebook file function preprocess_notebook(content) @@ -32,7 +32,7 @@ function create_files(title, file, repo_src, pages_dir, notebooks_dir; folder="" # Generate markdown file function preprocess_docs(content) - return string("# # [$title](@id $(splitext(file)[1]))\n $binder_badge\n $nbviewer_badge\n $download_badge\n\n", content) + return string("# # [$title](@id $(splitext(file)[1]))\n $binder_badge\n $nbviewer_badge\n $raw_notebook_badge\n\n", content) end Literate.markdown(joinpath(repo_src, folder, file), joinpath(pages_dir, folder); preprocess=preprocess_docs,) end diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl index d42695611f6..e259d25fb2f 100644 --- a/docs/literate/src/files/index.jl +++ b/docs/literate/src/files/index.jl @@ -5,9 +5,9 @@ # Right now, you are using the classic documentation. The corresponding interactive notebooks can # be opened in [Binder](https://mybinder.org/) and viewed in [nbviewer](https://nbviewer.jupyter.org/) -# via the icons ![](https://mybinder.org/badge_logo.svg) and ![](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg) +# via the icons ![](https://mybinder.org/badge_logo.svg) and ![](https://img.shields.io/badge/render-nbviewer-f37726) # in the respective tutorial. -# You can download the raw notebooks from GitHub via ![](https://camo.githubusercontent.com/aea75103f6d9f690a19cb0e17c06f984ab0f472d9e6fe4eadaa0cc438ba88ada/68747470733a2f2f696d672e736869656c64732e696f2f62616467652f646f776e6c6f61642d6e6f7465626f6f6b2d627269676874677265656e). +# You can also open the raw notebook files via ![](https://img.shields.io/badge/raw-notebook-4cc61e). # **Note:** To improve responsiveness via caching, the notebooks are updated only once a week. They are only # available for the latest stable release of Trixi.jl at the time of caching. diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl new file mode 100644 index 00000000000..05c09d57530 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -0,0 +1,146 @@ +using OrdinaryDiffEq +using Trixi + +# Warm bubble test case from +# - Wicker, L. J., and Skamarock, W. C. (1998) +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + @unpack g, c_p, c_v = setup + + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + +# Source terms +@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack g = setup + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. +surface_flux = FluxLMARS(340.0) + +volume_flux = flux_kennedy_gruber +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +cells_per_dimension = (64, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + source_terms = warm_bubble_setup, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl index 653bab41e2d..0589e76a6a9 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_solution_independent.jl @@ -93,12 +93,10 @@ solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) coordinates_min = (-5.0, -5.0) coordinates_max = (5.0, 5.0) -mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - trees_per_dimension = (1, 1) mesh = T8codeMesh(trees_per_dimension, polydeg = 3, - mapping = mapping, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 1) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index adf1d009a59..f285d24fc6c 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -41,9 +41,9 @@ isfile(mesh_file) || # we can create a t8code mesh. conn = Trixi.read_inp_p4est(mesh_file, Val(2)) -mesh = T8codeMesh{2}(conn, polydeg = 3, - mapping = mapping_flag, - initial_refinement_level = 1) +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 1) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = boundary_conditions) diff --git a/examples/t8code_2d_dgsem/elixir_advection_basic.jl b/examples/t8code_2d_dgsem/elixir_advection_basic.jl index efc51226586..26ced0970fe 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_basic.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_basic.jl @@ -16,12 +16,10 @@ solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) -mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - trees_per_dimension = (8, 8) mesh = T8codeMesh(trees_per_dimension, polydeg = 3, - mapping = mapping, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 1) # A semidiscretization collects data structures and functions for the spatial discretization diff --git a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl index 31a8bc93697..a39f3a7e195 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_nonconforming_flag.jl @@ -20,31 +20,28 @@ f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) faces = (f1, f2, f3, f4) mapping = Trixi.transfinite_mapping(faces) -# Create P4estMesh with 3 x 2 trees and 6 x 4 elements, +# Create T8codeMesh with 3 x 2 trees and 6 x 4 elements, # approximate the geometry with a smaller polydeg for testing. trees_per_dimension = (3, 2) mesh = T8codeMesh(trees_per_dimension, polydeg = 3, mapping = mapping, initial_refinement_level = 1) -function adapt_callback(forest, - forest_from, - which_tree, - lelement_id, - ts, - is_family, - num_elements, - elements_ptr)::Cint - vertex = Vector{Cdouble}(undef, 3) - - elements = unsafe_wrap(Array, elements_ptr, num_elements) - - Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) +# Note: This is actually a `p4est_quadrant_t` which is much bigger than the +# following struct. But we only need the first three fields for our purpose. +struct t8_dquad_t + x::Int32 + y::Int32 + level::Int8 + # [...] # See `p4est.h` in `p4est` for more info. +end - level = Trixi.t8_element_level(ts, elements[1]) +# Refine quadrants of each tree at lower left edge to level 4. +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dquad_t}(elements[1])) - # TODO: Make this condition more general. - if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 4 + if el.x == 0 && el.y == 0 && el.level < 4 # return true (refine) return 1 else @@ -53,26 +50,7 @@ function adapt_callback(forest, end end -Trixi.@T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest)!=0); - -# Init new forest. -new_forest_ref = Ref{Trixi.t8_forest_t}() -Trixi.t8_forest_init(new_forest_ref); -new_forest = new_forest_ref[] - -# Check out `examples/t8_step4_partition_balance_ghost.jl` in -# https://github.com/DLR-AMR/T8code.jl for detailed explanations. -let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0 - Trixi.t8_forest_set_user_data(new_forest, C_NULL) - Trixi.t8_forest_set_adapt(new_forest, mesh.forest, - Trixi.@t8_adapt_callback(adapt_callback), recursive) - Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition) - Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening) - Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES) - Trixi.t8_forest_commit(new_forest) -end - -mesh.forest = new_forest +Trixi.adapt!(mesh, adapt_callback) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, diff --git a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl index df9cbc26f6e..5ba1ab15489 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_unstructured_flag.jl @@ -38,9 +38,9 @@ isfile(mesh_file) || # we can create a t8code mesh. conn = Trixi.read_inp_p4est(mesh_file, Val(2)) -mesh = T8codeMesh{2}(conn, polydeg = 3, - mapping = mapping_flag, - initial_refinement_level = 2) +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 2) # A semidiscretization collects data structures and functions for the spatial discretization. semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, diff --git a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl index 01e0449c67e..37d15f38566 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_free_stream.jl @@ -40,25 +40,17 @@ isfile(mesh_file) || # we can create a t8code mesh. conn = Trixi.read_inp_p4est(mesh_file, Val(2)) -mesh = T8codeMesh{2}(conn, polydeg = 3, - mapping = mapping, - initial_refinement_level = 1) - -function adapt_callback(forest, - forest_from, - which_tree, - lelement_id, - ts, - is_family, - num_elements, - elements_ptr)::Cint - vertex = Vector{Cdouble}(undef, 3) +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping, + initial_refinement_level = 1) - elements = unsafe_wrap(Array, elements_ptr, num_elements) +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + vertex = Vector{Cdouble}(undef, 3) - Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) + Trixi.t8_element_vertex_reference_coords(eclass_scheme, elements[1], 0, vertex) - level = Trixi.t8_element_level(ts, elements[1]) + level = Trixi.t8_element_level(eclass_scheme, elements[1]) # TODO: Make this condition more general. if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 3 @@ -70,26 +62,7 @@ function adapt_callback(forest, end end -Trixi.@T8_ASSERT(Trixi.t8_forest_is_committed(mesh.forest)!=0); - -# Init new forest. -new_forest_ref = Ref{Trixi.t8_forest_t}() -Trixi.t8_forest_init(new_forest_ref); -new_forest = new_forest_ref[] - -# Check out `examples/t8_step4_partition_balance_ghost.jl` in -# https://github.com/DLR-AMR/T8code.jl for detailed explanations. -let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0 - Trixi.t8_forest_set_user_data(new_forest, C_NULL) - Trixi.t8_forest_set_adapt(new_forest, mesh.forest, - Trixi.@t8_adapt_callback(adapt_callback), recursive) - Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition) - Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening) - Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES) - Trixi.t8_forest_commit(new_forest) -end - -mesh.forest = new_forest +Trixi.adapt!(mesh, adapt_callback) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition))) diff --git a/examples/t8code_2d_dgsem/elixir_euler_sedov.jl b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl index 965d794f8dc..82770a4050b 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_sedov.jl @@ -58,12 +58,10 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, coordinates_min = (-1.0, -1.0) coordinates_max = (1.0, 1.0) -mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - trees_per_dimension = (4, 4) mesh = T8codeMesh(trees_per_dimension, polydeg = 4, - mapping = mapping, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 2, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl index 55a9063a001..9ebbd1d28c4 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_shockcapturing_ec.jl @@ -29,12 +29,10 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, coordinates_min = (-1.0, -1.0) coordinates_max = (1.0, 1.0) -mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - trees_per_dimension = (4, 4) mesh = T8codeMesh(trees_per_dimension, polydeg = 4, - mapping = mapping, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 2, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl index 21f26d79ba8..bcc1abc560e 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_flag.jl @@ -40,25 +40,17 @@ isfile(mesh_file) || # we can create a t8code mesh. conn = Trixi.read_inp_p4est(mesh_file, Val(2)) -mesh = T8codeMesh{2}(conn, polydeg = 3, - mapping = mapping_flag, - initial_refinement_level = 1) - -function adapt_callback(forest, - forest_from, - which_tree, - lelement_id, - ts, - is_family, - num_elements, - elements_ptr)::Cint - vertex = Vector{Cdouble}(undef, 3) +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping_flag, + initial_refinement_level = 1) - elements = unsafe_wrap(Array, elements_ptr, num_elements) +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + vertex = Vector{Cdouble}(undef, 3) - Trixi.t8_element_vertex_reference_coords(ts, elements[1], 0, pointer(vertex)) + Trixi.t8_element_vertex_reference_coords(eclass_scheme, elements[1], 0, pointer(vertex)) - level = Trixi.t8_element_level(ts, elements[1]) + level = Trixi.t8_element_level(eclass_scheme, elements[1]) # TODO: Make this condition more general. if vertex[1] < 1e-8 && vertex[2] < 1e-8 && level < 2 @@ -70,26 +62,7 @@ function adapt_callback(forest, end end -@assert(Trixi.t8_forest_is_committed(mesh.forest)!=0); - -# Init new forest. -new_forest_ref = Ref{Trixi.t8_forest_t}() -Trixi.t8_forest_init(new_forest_ref); -new_forest = new_forest_ref[] - -# Check out `examples/t8_step4_partition_balance_ghost.jl` in -# https://github.com/DLR-AMR/T8code.jl for detailed explanations. -let set_from = C_NULL, recursive = 1, set_for_coarsening = 0, no_repartition = 0 - Trixi.t8_forest_set_user_data(new_forest, C_NULL) - Trixi.t8_forest_set_adapt(new_forest, mesh.forest, - Trixi.@t8_adapt_callback(adapt_callback), recursive) - Trixi.t8_forest_set_balance(new_forest, set_from, no_repartition) - Trixi.t8_forest_set_partition(new_forest, set_from, set_for_coarsening) - Trixi.t8_forest_set_ghost(new_forest, 1, Trixi.T8_GHOST_FACES) - Trixi.t8_forest_commit(new_forest) -end - -mesh.forest = new_forest +Trixi.adapt!(mesh, adapt_callback) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms, diff --git a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl index 6d6bb27e0c3..98a9a5521a9 100644 --- a/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl +++ b/examples/t8code_2d_dgsem/elixir_eulergravity_convergence.jl @@ -16,10 +16,8 @@ coordinates_max = (2.0, 2.0) trees_per_dimension = (1, 1) -mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - mesh = T8codeMesh(trees_per_dimension, polydeg = 1, - mapping = mapping, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 2) semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition, diff --git a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl index 1e2362a123c..e184cb3fd05 100644 --- a/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl +++ b/examples/t8code_2d_dgsem/elixir_mhd_alfven_wave.jl @@ -11,18 +11,17 @@ initial_condition = initial_condition_convergence_test # Get the DG approximation space volume_flux = (flux_central, flux_nonconservative_powell) + solver = DGSEM(polydeg = 4, surface_flux = (flux_hlle, flux_nonconservative_powell), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) coordinates_min = (0.0, 0.0) coordinates_max = (sqrt(2.0), sqrt(2.0)) -mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - trees_per_dimension = (8, 8) mesh = T8codeMesh(trees_per_dimension, polydeg = 3, - mapping = mapping, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 0, periodicity = true) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl index 9a4bd99e444..adb154948fb 100644 --- a/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl +++ b/examples/t8code_2d_dgsem/elixir_mhd_rotor.jl @@ -78,9 +78,9 @@ isfile(mesh_file) || # we can create a t8code mesh. conn = Trixi.read_inp_p4est(mesh_file, Val(2)) -mesh = T8codeMesh{2}(conn, polydeg = 4, - mapping = mapping_twist, - initial_refinement_level = 1) +mesh = T8codeMesh(conn, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 1) boundary_condition = BoundaryConditionDirichlet(initial_condition) boundary_conditions = Dict(:all => boundary_condition) diff --git a/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl index b2d5097036f..3610639d554 100644 --- a/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl +++ b/examples/t8code_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -22,12 +22,10 @@ solver = DGSEM(polydeg = 3, coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y)) coordinates_max = (sqrt(2.0), sqrt(2.0)) # maximum coordinates (max(x), max(y)) -mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) - trees_per_dimension = (8, 8) mesh = T8codeMesh(trees_per_dimension, polydeg = 3, - mapping = mapping, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, initial_refinement_level = 1) # A semidiscretization collects data structures and functions for the spatial discretization diff --git a/examples/t8code_3d_dgsem/elixir_advection_amr.jl b/examples/t8code_3d_dgsem/elixir_advection_amr.jl new file mode 100644 index 00000000000..5a4b2218d57 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_amr.jl @@ -0,0 +1,66 @@ +# The same setup as tree_3d_dgsem/elixir_advection_amr.jl +# to verify the T8codeMesh implementation against TreeMesh. + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +initial_condition = initial_condition_gauss +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-5.0, -5.0, -5.0) +coordinates_max = (5.0, 5.0, 5.0) +trees_per_dimension = (1, 1, 1) + +# Note that it is not necessary to use mesh polydeg lower than the solver polydeg +# on a Cartesian mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +mesh = T8codeMesh(trees_per_dimension, polydeg = 1, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 4) + +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, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 4, + med_level = 5, med_threshold = 0.1, + max_level = 6, max_threshold = 0.6) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 1.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl new file mode 100644 index 00000000000..617736afbdd --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_amr_unstructured_curved.jl @@ -0,0 +1,105 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# Solver with polydeg=4 to ensure free stream preservation (FSP) on non-conforming meshes. +# The polydeg of the solver must be at least twice as big as the polydeg of the mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs) + +initial_condition = initial_condition_gauss +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +boundary_conditions = Dict(:all => boundary_condition) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. The original mapping applied to this unstructured mesh +# causes some Jacobians to be negative, which makes the mesh invalid. +function mapping(xi, eta, zeta) + # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries + # xi = 1.5 * xi_ + 1.5 + # eta = 1.5 * eta_ + 1.5 + # zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 4 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 4 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 4 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + # Transform the weird deformed cube to be approximately the size of [-5,5]^3 to match IC + return SVector(5 * x, 5 * y, 5 * z) +end + +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 2, + mapping = mapping, + initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 8.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first), + base_level = 1, + med_level = 2, med_threshold = 0.1, + max_level = 3, max_threshold = 0.6) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 1.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/t8code_3d_dgsem/elixir_advection_basic.jl b/examples/t8code_3d_dgsem/elixir_advection_basic.jl new file mode 100644 index 00000000000..f49462035aa --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_basic.jl @@ -0,0 +1,59 @@ +# The same setup as tree_3d_dgsem/elixir_advection_basic.jl +# to verify the T8codeMesh implementation against TreeMesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# 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, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# Create P4estMesh with 8 x 8 x 8 elements (note `refinement_level=1`) +trees_per_dimension = (4, 4, 4) +mesh = T8codeMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 1) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +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() + +# 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 = 1.2) + +# 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/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl b/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl new file mode 100644 index 00000000000..8d7a48370f5 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_nonconforming.jl @@ -0,0 +1,85 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the linear advection equation. + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# 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, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) +trees_per_dimension = (1, 1, 1) + +# Note that it is not necessary to use mesh polydeg lower than the solver polydeg +# on a Cartesian mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +mesh = T8codeMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 2) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +# Refine bottom left quadrant of each second tree to level 2 +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if iseven(convert(Int, ltreeid)) && el.x == 0 && el.y == 0 && el.z == 0 && + el.level < 3 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +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() + +# 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 = 1.6) + +# 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/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl new file mode 100644 index 00000000000..df358435c9a --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_advection_unstructured_curved.jl @@ -0,0 +1,98 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +# 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) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. The original mapping applied to this unstructured mesh +# causes some Jacobians to be negative, which makes the mesh invalid. +function mapping(xi, eta, zeta) + # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries + # xi = 1.5 * xi_ + 1.5 + # eta = 1.5 * eta_ + 1.5 + # zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 6 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 6 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 6 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + return SVector(x, y, z) +end + +# Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping, + initial_refinement_level = 2) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 0.1 +ode = semidiscretize(semi, (0.0, 0.1)); + +# 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_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +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 = 1.2) + +# 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); + +# Print the timer summary +summary_callback() diff --git a/examples/t8code_3d_dgsem/elixir_euler_ec.jl b/examples/t8code_3d_dgsem/elixir_euler_ec.jl new file mode 100644 index 00000000000..07745c3ac56 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_ec.jl @@ -0,0 +1,92 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(5 / 3) + +initial_condition = initial_condition_weak_blast_wave + +boundary_conditions = Dict(:all => boundary_condition_slip_wall) + +# Get the DG approximation space + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 5, surface_flux = flux_ranocha, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# Get the curved quad mesh from a file + +# Mapping as described in https://arxiv.org/abs/2012.12040 +function mapping(xi_, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 3 / 8 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 3 / 8 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 3 / 8 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + return SVector(x, y, z) +end + +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 5, + mapping = mapping, + initial_refinement_level = 0) + +# Create the semidiscretization object. +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_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/t8code_3d_dgsem/elixir_euler_free_stream.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl new file mode 100644 index 00000000000..e135d464810 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream.jl @@ -0,0 +1,118 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_constant + +boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition)) + +# Solver with polydeg=4 to ensure free stream preservation (FSP) on non-conforming meshes. +# The polydeg of the solver must be at least twice as big as the polydeg of the mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. This can yield problematic geometries if the unrefined mesh +# is not fine enough. +function mapping(xi_, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 6 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 6 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 6 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + return SVector(x, y, z) +end + +# Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connectivity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 2, + mapping = mapping, + initial_refinement_level = 0) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +# Refine bottom left quadrant of each second tree to level 2 +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if iseven(convert(Int, ltreeid)) && el.x == 0 && el.y == 0 && el.z == 0 && + el.level < 2 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_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/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl new file mode 100644 index 00000000000..d129b59826e --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_free_stream_extruded.jl @@ -0,0 +1,106 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_constant + +boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition)) + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D. +# This particular mesh is unstructured in the yz-plane, but extruded in x-direction. +# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded +# in x-direction to ensure free stream preservation on a non-conforming mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +function mapping(xi, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + z = zeta + + 1 / 6 * (cos(1.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + y = eta + 1 / 6 * (cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(2 * pi * (2 * z - 3) / 3)) + + return SVector(xi, y, z) +end + +# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connecvity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +mesh = T8codeMesh(conn, polydeg = 3, + mapping = mapping, + initial_refinement_level = 0) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +# Refine quadrants in y-direction of each tree at one edge to level 2 +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if convert(Int, ltreeid) < 4 && el.x == 0 && el.y == 0 && el.level < 2 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), #maxiters=1, + 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/t8code_3d_dgsem/elixir_euler_sedov.jl b/examples/t8code_3d_dgsem/elixir_euler_sedov.jl new file mode 100644 index 00000000000..618b170b661 --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_sedov.jl @@ -0,0 +1,97 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +""" + initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) + +The Sedov blast wave setup based on Flash +- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 +with smaller strength of the initial discontinuity. +""" +function initial_condition_medium_sedov_blast_wave(x, t, + equations::CompressibleEulerEquations3D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 + r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) + p0_outer = 1.0e-3 + + # Calculate primitive variables + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +initial_condition = initial_condition_medium_sedov_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 5 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +coordinates_min = (-1.0, -1.0, -1.0) +coordinates_max = (1.0, 1.0, 1.0) + +trees_per_dimension = (4, 4, 4) +mesh = T8codeMesh(trees_per_dimension, + polydeg = 4, initial_refinement_level = 0, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = true) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_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/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl new file mode 100644 index 00000000000..d4664522bea --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl @@ -0,0 +1,120 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) + +# Solver with polydeg=4 to ensure free stream preservation (FSP) on non-conforming meshes. +# The polydeg of the solver must be at least twice as big as the polydeg of the mesh. +# See https://doi.org/10.1007/s10915-018-00897-9, Section 6. +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. The original mapping applied to this unstructured mesh +# causes some Jacobians to be negative, which makes the mesh invalid. +function mapping(xi, eta, zeta) + # Don't transform input variables between -1 and 1 onto [0,3] to obtain curved boundaries + # xi = 1.5 * xi_ + 1.5 + # eta = 1.5 * eta_ + 1.5 + # zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 1 / 6 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 1 / 6 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 1 / 6 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + # Transform the weird deformed cube to be approximately the cube [0,2]^3 + return SVector(x + 1, y + 1, z + 1) +end + +# Unstructured mesh with 68 cells of the cube domain [-1, 1]^3 +mesh_file = joinpath(@__DIR__, "cube_unstructured_1.inp") +isfile(mesh_file) || + download("https://gist.githubusercontent.com/efaulhaber/d45c8ac1e248618885fa7cc31a50ab40/raw/37fba24890ab37cfa49c39eae98b44faf4502882/cube_unstructured_1.inp", + mesh_file) + +# INP mesh files are only support by p4est. Hence, we +# create a p4est connecvity object first from which +# we can create a t8code mesh. +conn = Trixi.read_inp_p4est(mesh_file, Val(3)) + +# Mesh polydeg of 2 (half the solver polydeg) to ensure FSP (see above). +mesh = T8codeMesh(conn, polydeg = 2, + mapping = mapping, + initial_refinement_level = 0) + +# Note: This is actually a `p8est_quadrant_t` which is much bigger than the +# following struct. But we only need the first four fields for our purpose. +struct t8_dhex_t + x::Int32 + y::Int32 + z::Int32 + level::Int8 + # [...] # See `p8est.h` in `p4est` for more info. +end + +function adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, + user_data) + el = unsafe_load(Ptr{t8_dhex_t}(elements[1])) + + if el.x == 0 && el.y == 0 && el.z == 0 && el.level < 2 + # return true (refine) + return 1 + else + # return false (don't refine) + return 0 + end +end + +Trixi.adapt!(mesh, adapt_callback) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.045) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.6) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_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/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl new file mode 100644 index 00000000000..7cb03bb312d --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_source_terms_nonperiodic.jl @@ -0,0 +1,62 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:x_neg => boundary_condition, + :x_pos => boundary_condition, + :y_neg => boundary_condition, + :y_pos => boundary_condition, + :z_neg => boundary_condition, + :z_pos => boundary_condition) + +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralWeakForm()) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (2.0, 2.0, 2.0) + +trees_per_dimension = (2, 2, 2) + +mapping = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +mesh = T8codeMesh(trees_per_dimension, polydeg = 1, + mapping = mapping, + periodicity = false, initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.6) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl new file mode 100644 index 00000000000..f2e14273ae7 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl @@ -0,0 +1,150 @@ +using OrdinaryDiffEq +using Trixi + +# Warm bubble test case from +# - Wicker, L. J., and Skamarock, W. C. (1998) +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + @unpack g, c_p, c_v = setup + + # center of perturbation + center_x = 10000.0 + center_z = 2000.0 + # radius of perturbation + radius = 2000.0 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300.0 + potential_temperature_perturbation = 0.0 + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + +# Source terms +@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack g = setup + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. +surface_flux = FluxLMARS(340.0) + +volume_flux = flux_kennedy_gruber +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0) +coordinates_max = (20_000.0, 10_000.0) + +# Same coordinates as in examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +# However TreeMesh will generate a 20_000 x 20_000 square domain instead +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000, + periodicity = (true, false)) + +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + source_terms = warm_bubble_setup, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # 1000 seconds final time + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:entropy_conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + output_directory = "out", + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + maxiters = 1.0e7, + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl new file mode 100644 index 00000000000..b8dbad50680 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -0,0 +1,82 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) + +# An initial condition with a bottom topography and a perturbation in the waterheight to test +# boundary_condition_slip_wall +function initial_condition_perturbation(x, t, equations::ShallowWaterEquations2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + + # Bottom topography + b = 1.5 * exp(-0.5 * ((x[1])^2 + (x[2])^2)) + # Waterheight perturbation + H = H + 0.5 * exp(-10.0 * ((x[1])^2 + (x[2])^2)) + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_perturbation + +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.25) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl index bd781b21c1e..db01476bb86 100644 --- a/src/auxiliary/t8code.jl +++ b/src/auxiliary/t8code.jl @@ -35,7 +35,7 @@ function init_t8code() # production runs this is not mandatory, but is helpful during # development. Hence, this option is only activated when environment # variable TRIXI_T8CODE_SC_FINALIZE exists. - @warn "T8code.jl: sc_finalize will be called during shutdown of Trixi.jl." + @info "T8code.jl: `sc_finalize` will be called during shutdown of Trixi.jl." MPI.add_finalize_hook!(T8code.Libt8.sc_finalize) end else @@ -116,7 +116,6 @@ function trixi_t8_count_interfaces(forest) elseif level < neighbor_level local_num_mortars += 1 end - else local_num_boundary += 1 end @@ -219,38 +218,9 @@ function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundari interfaces.neighbor_ids[1, interface_id] = current_index + 1 interfaces.neighbor_ids[2, interface_id] = neighbor_ielements[1] + 1 - # Iterate over primary and secondary element. - for side in 1:2 - # Align interface in positive coordinate direction of primary element. - # For orientation == 1, the secondary element needs to be indexed backwards - # relative to the interface. - if side == 1 || orientation == 0 - # Forward indexing - indexing = :i_forward - else - # Backward indexing - indexing = :i_backward - end - - if faces[side] == 0 - # Index face in negative x-direction - interfaces.node_indices[side, interface_id] = (:begin, - indexing) - elseif faces[side] == 1 - # Index face in positive x-direction - interfaces.node_indices[side, interface_id] = (:end, - indexing) - elseif faces[side] == 2 - # Index face in negative y-direction - interfaces.node_indices[side, interface_id] = (indexing, - :begin) - else # faces[side] == 3 - # Index face in positive y-direction - interfaces.node_indices[side, interface_id] = (indexing, - :end) - end - end - + # Save interfaces.node_indices dimension specific in containers_3d.jl. + init_interface_node_indices!(interfaces, faces, orientation, + interface_id) # Non-conforming interface. elseif level < neighbor_level local_num_mortars += 1 @@ -262,42 +232,13 @@ function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundari # Last entry is the large element. mortars.neighbor_ids[end, mortar_id] = current_index + 1 - # First `1:end-1` entries are the smaller elements. - mortars.neighbor_ids[1:(end - 1), mortar_id] .= neighbor_ielements .+ - 1 - - for side in 1:2 - # Align mortar in positive coordinate direction of small side. - # For orientation == 1, the large side needs to be indexed backwards - # relative to the mortar. - if side == 1 || orientation == 0 - # Forward indexing for small side or orientation == 0. - indexing = :i_forward - else - # Backward indexing for large side with reversed orientation. - indexing = :i_backward - # Since the orientation is reversed we have to account for this - # when filling the `neighbor_ids` array. - mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + - 1 - mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + - 1 - end - - if faces[side] == 0 - # Index face in negative x-direction - mortars.node_indices[side, mortar_id] = (:begin, indexing) - elseif faces[side] == 1 - # Index face in positive x-direction - mortars.node_indices[side, mortar_id] = (:end, indexing) - elseif faces[side] == 2 - # Index face in negative y-direction - mortars.node_indices[side, mortar_id] = (indexing, :begin) - else # faces[side] == 3 - # Index face in positive y-direction - mortars.node_indices[side, mortar_id] = (indexing, :end) - end - end + # Fill in the `mortars.neighbor_ids` array and reorder if necessary. + init_mortar_neighbor_ids!(mortars, faces[2], faces[1], + orientation, neighbor_ielements, + mortar_id) + + # Fill in the `mortars.node_indices` array. + init_mortar_node_indices!(mortars, faces, orientation, mortar_id) # else: "level > neighbor_level" is skipped since we visit the mortar interface only once. end @@ -309,19 +250,7 @@ function trixi_t8_fill_mesh_info(forest, elements, interfaces, mortars, boundari boundaries.neighbor_ids[boundary_id] = current_index + 1 - if iface == 0 - # Index face in negative x-direction. - boundaries.node_indices[boundary_id] = (:begin, :i_forward) - elseif iface == 1 - # Index face in positive x-direction. - boundaries.node_indices[boundary_id] = (:end, :i_forward) - elseif iface == 2 - # Index face in negative y-direction. - boundaries.node_indices[boundary_id] = (:i_forward, :begin) - else # iface == 3 - # Index face in positive y-direction. - boundaries.node_indices[boundary_id] = (:i_forward, :end) - end + init_boundary_node_indices!(boundaries, iface, boundary_id) # One-based indexing. boundaries.name[boundary_id] = boundary_names[iface + 1, itree + 1] @@ -420,13 +349,15 @@ function trixi_t8_adapt_new(old_forest, indicators) t8_forest_init(new_forest_ref) new_forest = new_forest_ref[] - let set_from = C_NULL, recursive = 0, set_for_coarsening = 0, no_repartition = 0 + let set_from = C_NULL, recursive = 0, set_for_coarsening = 0, no_repartition = 0, + do_ghost = 1 + t8_forest_set_user_data(new_forest, pointer(indicators)) t8_forest_set_adapt(new_forest, old_forest, @t8_adapt_callback(adapt_callback), recursive) t8_forest_set_balance(new_forest, set_from, no_repartition) t8_forest_set_partition(new_forest, set_from, set_for_coarsening) - t8_forest_set_ghost(new_forest, 1, T8_GHOST_FACES) # Note: MPI support not available yet so it is a dummy call. + t8_forest_set_ghost(new_forest, do_ghost, T8_GHOST_FACES) # Note: MPI support not available yet so it is a dummy call. t8_forest_commit(new_forest) end diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 98e531295b7..b816bc06e65 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -396,7 +396,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, old_index = 1 new_index = 1 - # Note: This is true for `quads` only. + # Note: This is true for `quads`. T8_CHILDREN = 4 # Retain current solution data. diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index c8abe6fdb05..392cbba9e28 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -304,9 +304,84 @@ end # this method is called when an `ControllerThreeLevel` is constructed function create_cache(::Type{ControllerThreeLevel}, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, dg::DG, cache) controller_value = Vector{Int}(undef, nelements(dg, cache)) return (; controller_value) end + +# Coarsen and refine elements in the DG solver based on a difference list. +function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, + dg::DGSEM, cache, difference) + + # Return early if there is nothing to do. + if !any(difference .!= 0) + return nothing + end + + # Number of (local) cells/elements. + old_nelems = nelements(dg, cache) + new_nelems = ncells(mesh) + + # Local element indices. + old_index = 1 + new_index = 1 + + # Note: This is only true for `hexs`. + T8_CHILDREN = 8 + + # Retain current solution data. + old_u_ode = copy(u_ode) + + GC.@preserve old_u_ode begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * ndofs(mesh, dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + + while old_index <= old_nelems && new_index <= new_nelems + if difference[old_index] > 0 # Refine. + + # Refine element and store solution directly in new data structure. + refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg, + u_tmp1, u_tmp2) + + old_index += 1 + new_index += T8_CHILDREN + + elseif difference[old_index] < 0 # Coarsen. + + # If an element is to be removed, sanity check if the following elements + # are also marked - otherwise there would be an error in the way the + # cells/elements are sorted. + @assert all(difference[old_index:(old_index + T8_CHILDREN - 1)] .< 0) "bad cell/element order" + + # Coarsen elements and store solution directly in new data structure. + coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, + dg, u_tmp1, u_tmp2) + + old_index += T8_CHILDREN + new_index += 1 + + else # No changes. + + # Copy old element data to new element container. + @views u[:, .., new_index] .= old_u[:, .., old_index] + + old_index += 1 + new_index += 1 + end + end # while + end # GC.@preserve old_u_ode + + return nothing +end end # @muladd diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 81d0795a159..27e8a2b722f 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -35,7 +35,9 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{3}, return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2) end -function create_cache_analysis(analyzer, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, +function create_cache_analysis(analyzer, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DG, cache, RealT, uEltype) @@ -118,7 +120,7 @@ function calc_error_norms(func, u, t, analyzer, end function calc_error_norms(func, u, t, analyzer, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @@ -190,7 +192,8 @@ function integrate_via_indices(func::Func, u, end function integrate_via_indices(func::Func, u, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis @@ -218,7 +221,8 @@ function integrate_via_indices(func::Func, u, end function integrate(func::Func, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DG, cache; normalize = true) where {Func} integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, k, element, equations, dg @@ -248,7 +252,8 @@ function integrate(func::Func, u, end function analyze(::typeof(entropy_timederivative), du, u, t, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, dg::DG, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, @@ -277,7 +282,7 @@ function analyze(::Val{:l2_divb}, du, u, t, end function analyze(::Val{:l2_divb}, du, u, t, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations::IdealGlmMhdEquations3D, dg::DGSEM, cache) @unpack contravariant_vectors = cache.elements @@ -333,7 +338,7 @@ function analyze(::Val{:linf_divb}, du, u, t, end function analyze(::Val{:linf_divb}, du, u, t, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations::IdealGlmMhdEquations3D, dg::DGSEM, cache) @unpack derivative_matrix, weights = dg.basis diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl index c9ab7c478a8..822ab2f87ec 100644 --- a/src/callbacks_step/stepsize_dg3d.jl +++ b/src/callbacks_step/stepsize_dg3d.jl @@ -44,7 +44,7 @@ function max_dt(u, t, mesh::TreeMesh{3}, return 2 / (nnodes(dg) * max_scaled_speed) end -function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, +function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -82,7 +82,7 @@ function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, return 2 / (nnodes(dg) * max_scaled_speed) end -function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}}, +function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, constant_speed::True, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index b0fd5c53f45..3c6f759db2b 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -809,6 +809,98 @@ end return SVector(f1m, f2m, f3m, f4m) end +""" + FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction, + equations::CompressibleEulerEquations2D) + +Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using +an estimate `c` of the speed of sound. + +References: +- Xi Chen et al. (2013) + A Control-Volume Model of the Compressible Euler Equations with a Vertical + Lagrangian Coordinate + [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) +""" +struct FluxLMARS{SpeedOfSound} + # Estimate for the speed of sound + speed_of_sound::SpeedOfSound +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquations2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + if orientation == 1 + v_ll = v1_ll + v_rr = v1_rr + else # orientation == 2 + v_ll = v2_ll + v_rr = v2_rr + end + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4 = v * u_ll + f4 = f4 + p_ll * v + else + f1, f2, f3, f4 = v * u_rr + f4 = f4 + p_rr * v + end + + if orientation == 1 + f2 = f2 + p + else # orientation == 2 + f3 = f3 + p + end + + return SVector(f1, f2, f3, f4) +end + +@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + c = flux_lmars.speed_of_sound + + # Unpack left and right state + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # Note that this is the same as computing v_ll and v_rr with a normalized normal vector + # and then multiplying v by `norm_` again, but this version is slightly faster. + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ + v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p + if v >= 0 + f1, f2, f3, f4 = u_ll * v + f4 = f4 + p_ll * v + else + f1, f2, f3, f4 = u_rr * v + f4 = f4 + p_rr * v + end + + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4) +end + """ splitting_vanleer_haenel(u, orientation::Integer, equations::CompressibleEulerEquations2D) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 82c4a7efa32..292b912f009 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -944,11 +944,6 @@ References: Lagrangian Coordinate [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1) """ -struct FluxLMARS{SpeedOfSound} - # Estimate for the speed of sound - speed_of_sound::SpeedOfSound -end - @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerEquations3D) c = flux_lmars.speed_of_sound @@ -972,10 +967,14 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll + f5 = f5 + p_ll * v else f1, f2, f3, f4, f5 = v * u_rr + f5 = f5 + p_rr * v end if orientation == 1 @@ -985,7 +984,6 @@ end else # orientation == 3 f4 += p end - f5 += p * v return SVector(f1, f2, f3, f4, f5) end @@ -1011,18 +1009,21 @@ end p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_ v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_ + # We treat the energy term analogous to the potential temperature term in the paper by + # Chen et al., i.e. we use p_ll and p_rr, and not p if v >= 0 f1, f2, f3, f4, f5 = v * u_ll + f5 = f5 + p_ll * v else f1, f2, f3, f4, f5 = v * u_rr + f5 = f5 + p_rr * v end - f2 += p * normal_direction[1] - f3 += p * normal_direction[2] - f4 += p * normal_direction[3] - f5 += p * v - - return SVector(f1, f2, f3, f4, f5) + return SVector(f1, + f2 + p * normal_direction[1], + f3 + p * normal_direction[2], + f4 + p * normal_direction[3], + f5) end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index e75c92a27d0..6728d7d5553 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -236,8 +236,12 @@ Should be used together with [`TreeMesh`](@ref). u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4]) end - # compute and return the flux using `boundary_condition_slip_wall` routine above - flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + # 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 diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 13edcc29711..c9665a22af9 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -115,19 +115,49 @@ Non-periodic boundaries will be called ':x_neg', ':x_pos', ':y_neg', ':y_pos', ' - 'polydeg::Integer': polynomial degree used to store the geometry of the mesh. The mapping will be approximated by an interpolation polynomial of the specified degree for each tree. -- 'mapping': a function of 'NDIMS' variables to describe the mapping that transforms - the reference mesh ('[-1, 1]^n') to the physical domain. +- `mapping`: a function of `NDIMS` variables to describe the mapping that transforms + the reference mesh (`[-1, 1]^n`) to the physical domain. + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. +- `faces::NTuple{2*NDIMS}`: a tuple of `2 * NDIMS` functions that describe the faces of the domain. + Each function must take `NDIMS-1` arguments. + `faces[1]` describes the face onto which the face in negative x-direction + of the unit hypercube is mapped. The face in positive x-direction of + the unit hypercube will be mapped onto the face described by `faces[2]`. + `faces[3:4]` describe the faces in positive and negative y-direction respectively + (in 2D and 3D). + `faces[5:6]` describe the faces in positive and negative z-direction respectively (in 3D). + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. +- `coordinates_min`: vector or tuple of the coordinates of the corner in the negative direction of each dimension + to create a rectangular mesh. + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. +- `coordinates_max`: vector or tuple of the coordinates of the corner in the positive direction of each dimension + to create a rectangular mesh. + Use only one of `mapping`, `faces` and `coordinates_min`/`coordinates_max`. - 'RealT::Type': the type that should be used for coordinates. - 'initial_refinement_level::Integer': refine the mesh uniformly to this level before the simulation starts. - 'periodicity': either a 'Bool' deciding if all of the boundaries are periodic or an 'NTuple{NDIMS, Bool}' deciding for each dimension if the boundaries in this dimension are periodic. """ -function T8codeMesh(trees_per_dimension; polydeg, - mapping = coordinates2mapping((-1.0, -1.0), (1.0, 1.0)), - RealT = Float64, initial_refinement_level = 0, periodicity = true) - NDIMS = length(trees_per_dimension) +function T8codeMesh(trees_per_dimension; polydeg = 1, + mapping = nothing, faces = nothing, coordinates_min = nothing, + coordinates_max = nothing, + RealT = Float64, initial_refinement_level = 0, + periodicity = true) + @assert ((coordinates_min === nothing)===(coordinates_max === nothing)) "Either both or none of coordinates_min and coordinates_max must be specified" + + @assert count(i -> i !== nothing, + (mapping, faces, coordinates_min))==1 "Exactly one of mapping, faces and coordinates_min/max must be specified" + + # Extract mapping + if faces !== nothing + validate_faces(faces) + mapping = transfinite_mapping(faces) + elseif coordinates_min !== nothing + mapping = coordinates2mapping(coordinates_min, coordinates_max) + end - @assert NDIMS == 2 # Only support for NDIMS = 2 yet. + NDIMS = length(trees_per_dimension) + @assert (NDIMS == 2||NDIMS == 3) "NDIMS should be 2 or 3." # Convert periodicity to a Tuple of a Bool for every dimension if all(periodicity) @@ -141,10 +171,18 @@ function T8codeMesh(trees_per_dimension; polydeg, periodicity = Tuple(periodicity) end - conn = T8code.Libt8.p4est_connectivity_new_brick(trees_per_dimension..., periodicity...) do_partition = 0 - cmesh = t8_cmesh_new_from_p4est(conn, mpi_comm(), do_partition) - T8code.Libt8.p4est_connectivity_destroy(conn) + if NDIMS == 2 + conn = T8code.Libt8.p4est_connectivity_new_brick(trees_per_dimension..., + periodicity...) + cmesh = t8_cmesh_new_from_p4est(conn, mpi_comm(), do_partition) + T8code.Libt8.p4est_connectivity_destroy(conn) + elseif NDIMS == 3 + conn = T8code.Libt8.p8est_connectivity_new_brick(trees_per_dimension..., + periodicity...) + cmesh = t8_cmesh_new_from_p8est(conn, mpi_comm(), do_partition) + T8code.Libt8.p8est_connectivity_destroy(conn) + end scheme = t8_scheme_new_default_cxx() forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) @@ -156,28 +194,48 @@ function T8codeMesh(trees_per_dimension; polydeg, ntuple(_ -> length(nodes), NDIMS)..., prod(trees_per_dimension)) - # Get cell length in reference mesh: Omega_ref = [-1,1]^2. - dx = 2 / trees_per_dimension[1] - dy = 2 / trees_per_dimension[2] + # Get cell length in reference mesh: Omega_ref = [-1,1]^NDIMS. + dx = [2 / n for n in trees_per_dimension] num_local_trees = t8_cmesh_get_num_local_trees(cmesh) # Non-periodic boundaries. boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) + if mapping === nothing + mapping_ = coordinates2mapping(ntuple(_ -> -1.0, NDIMS), ntuple(_ -> 1.0, NDIMS)) + else + mapping_ = mapping + end + for itree in 1:num_local_trees veptr = t8_cmesh_get_tree_vertices(cmesh, itree - 1) verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS)) # Calculate node coordinates of reference mesh. - cell_x_offset = (verts[1, 1] - 1 / 2 * (trees_per_dimension[1] - 1)) * dx - cell_y_offset = (verts[2, 1] - 1 / 2 * (trees_per_dimension[2] - 1)) * dy - - for j in eachindex(nodes), i in eachindex(nodes) - tree_node_coordinates[:, i, j, itree] .= mapping(cell_x_offset + - dx * nodes[i] / 2, - cell_y_offset + - dy * nodes[j] / 2) + if NDIMS == 2 + cell_x_offset = (verts[1, 1] - 0.5 * (trees_per_dimension[1] - 1)) * dx[1] + cell_y_offset = (verts[2, 1] - 0.5 * (trees_per_dimension[2] - 1)) * dx[2] + + for j in eachindex(nodes), i in eachindex(nodes) + tree_node_coordinates[:, i, j, itree] .= mapping_(cell_x_offset + + dx[1] * nodes[i] / 2, + cell_y_offset + + dx[2] * nodes[j] / 2) + end + elseif NDIMS == 3 + cell_x_offset = (verts[1, 1] - 0.5 * (trees_per_dimension[1] - 1)) * dx[1] + cell_y_offset = (verts[2, 1] - 0.5 * (trees_per_dimension[2] - 1)) * dx[2] + cell_z_offset = (verts[3, 1] - 0.5 * (trees_per_dimension[3] - 1)) * dx[3] + + for k in eachindex(nodes), j in eachindex(nodes), i in eachindex(nodes) + tree_node_coordinates[:, i, j, k, itree] .= mapping_(cell_x_offset + + dx[1] * nodes[i] / 2, + cell_y_offset + + dx[2] * nodes[j] / 2, + cell_z_offset + + dx[3] * nodes[k] / 2) + end end if !periodicity[1] @@ -189,6 +247,13 @@ function T8codeMesh(trees_per_dimension; polydeg, boundary_names[3, itree] = :y_neg boundary_names[4, itree] = :y_pos end + + if NDIMS > 2 + if !periodicity[3] + boundary_names[5, itree] = :z_neg + boundary_names[6, itree] = :z_pos + end + end end return T8codeMesh{NDIMS}(cmesh, scheme, forest, tree_node_coordinates, nodes, @@ -196,9 +261,9 @@ function T8codeMesh(trees_per_dimension; polydeg, end """ - T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}, - mapping=nothing, polydeg=1, RealT=Float64, - initial_refinement_level=0) + T8codeMesh(cmesh::Ptr{t8_cmesh}, + mapping=nothing, polydeg=1, RealT=Float64, + initial_refinement_level=0) Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming mesh from a `t8_cmesh` data structure. @@ -215,10 +280,15 @@ conforming mesh from a `t8_cmesh` data structure. - `RealT::Type`: the type that should be used for coordinates. - `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. """ -function T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}; - mapping = nothing, polydeg = 1, RealT = Float64, - initial_refinement_level = 0) where {NDIMS} - @assert NDIMS == 2 # Only support for NDIMS = 2 yet. +function T8codeMesh(cmesh::Ptr{t8_cmesh}; + mapping = nothing, polydeg = 1, RealT = Float64, + initial_refinement_level = 0) + @assert (t8_cmesh_get_num_trees(cmesh)>0) "Given `cmesh` does not contain any trees." + + # Infer NDIMS from the geometry of the first tree. + NDIMS = Int(t8_geom_get_dimension(t8_cmesh_get_tree_geometry(cmesh, 0))) + + @assert (NDIMS == 2||NDIMS == 3) "NDIMS should be 2 or 3." scheme = t8_scheme_new_default_cxx() forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, 0, mpi_comm()) @@ -234,34 +304,84 @@ function T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}; nodes_in = [-1.0, 1.0] matrix = polynomial_interpolation_matrix(nodes_in, nodes) - data_in = Array{RealT, 3}(undef, 2, 2, 2) - tmp1 = zeros(RealT, 2, length(nodes), length(nodes_in)) - for itree in 0:(num_local_trees - 1) - veptr = t8_cmesh_get_tree_vertices(cmesh, itree) - verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS)) + if NDIMS == 2 + data_in = Array{RealT, 3}(undef, 2, 2, 2) + tmp1 = zeros(RealT, 2, length(nodes), length(nodes_in)) + verts = zeros(3, 4) - u = verts[:, 2] - verts[:, 1] - v = verts[:, 3] - verts[:, 1] - w = [0.0, 0.0, 1.0] + for itree in 0:(num_local_trees - 1) + veptr = t8_cmesh_get_tree_vertices(cmesh, itree) - vol = dot(cross(u, v), w) + # Note, `verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS))` + # sometimes does not work since `veptr` is not necessarily properly + # aligned to 8 bytes. + for icorner in 1:4 + verts[1, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 1) + verts[2, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 2) + end - if vol < 0.0 - @warn "Discovered negative volumes in `cmesh`: vol = $vol" + # Check if tree's node ordering is right-handed or print a warning. + let z = zero(eltype(verts)), o = one(eltype(verts)) + u = verts[:, 2] - verts[:, 1] + v = verts[:, 3] - verts[:, 1] + w = [z, z, o] + + # Triple product gives signed volume of spanned parallelepiped. + vol = dot(cross(u, v), w) + + if vol < z + @warn "Discovered negative volumes in `cmesh`: vol = $vol" + end + end + + # Tree vertices are stored in z-order. + @views data_in[:, 1, 1] .= verts[1:2, 1] + @views data_in[:, 2, 1] .= verts[1:2, 2] + @views data_in[:, 1, 2] .= verts[1:2, 3] + @views data_in[:, 2, 2] .= verts[1:2, 4] + + # Interpolate corner coordinates to specified nodes. + multiply_dimensionwise!(view(tree_node_coordinates, :, :, :, itree + 1), + matrix, matrix, + data_in, + tmp1) end - # Tree vertices are stored in z-order. - @views data_in[:, 1, 1] .= verts[1:2, 1] - @views data_in[:, 2, 1] .= verts[1:2, 2] - @views data_in[:, 1, 2] .= verts[1:2, 3] - @views data_in[:, 2, 2] .= verts[1:2, 4] - - # Interpolate corner coordinates to specified nodes. - multiply_dimensionwise!(view(tree_node_coordinates, :, :, :, itree + 1), - matrix, matrix, - data_in, - tmp1) + elseif NDIMS == 3 + data_in = Array{RealT, 4}(undef, 3, 2, 2, 2) + tmp1 = zeros(RealT, 3, length(nodes), length(nodes_in), length(nodes_in)) + verts = zeros(3, 8) + + for itree in 0:(num_local_trees - 1) + veptr = t8_cmesh_get_tree_vertices(cmesh, itree) + + # Note, `verts = unsafe_wrap(Array, veptr, (3, 1 << NDIMS))` + # sometimes does not work since `veptr` is not necessarily properly + # aligned to 8 bytes. + for icorner in 1:8 + verts[1, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 1) + verts[2, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 2) + verts[3, icorner] = unsafe_load(veptr, (icorner - 1) * 3 + 3) + end + + # Tree vertices are stored in z-order. + @views data_in[:, 1, 1, 1] .= verts[1:3, 1] + @views data_in[:, 2, 1, 1] .= verts[1:3, 2] + @views data_in[:, 1, 2, 1] .= verts[1:3, 3] + @views data_in[:, 2, 2, 1] .= verts[1:3, 4] + + @views data_in[:, 1, 1, 2] .= verts[1:3, 5] + @views data_in[:, 2, 1, 2] .= verts[1:3, 6] + @views data_in[:, 1, 2, 2] .= verts[1:3, 7] + @views data_in[:, 2, 2, 2] .= verts[1:3, 8] + + # Interpolate corner coordinates to specified nodes. + multiply_dimensionwise!(view(tree_node_coordinates, :, :, :, :, itree + 1), + matrix, matrix, matrix, + data_in, + tmp1) + end end map_node_coordinates!(tree_node_coordinates, mapping) @@ -274,9 +394,7 @@ function T8codeMesh{NDIMS}(cmesh::Ptr{t8_cmesh}; end """ - T8codeMesh{NDIMS}(conn::Ptr{p4est_connectivity}, - mapping=nothing, polydeg=1, RealT=Float64, - initial_refinement_level=0) + T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...) Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming mesh from a `p4est_connectivity` data structure. @@ -293,24 +411,45 @@ conforming mesh from a `p4est_connectivity` data structure. - `RealT::Type`: the type that should be used for coordinates. - `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. """ -function T8codeMesh{NDIMS}(conn::Ptr{p4est_connectivity}; kwargs...) where {NDIMS} - @assert NDIMS == 2 # Only support for NDIMS = 2 yet. - +function T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...) cmesh = t8_cmesh_new_from_p4est(conn, mpi_comm(), 0) - return T8codeMesh{NDIMS}(cmesh; kwargs...) + return T8codeMesh(cmesh; kwargs...) +end + +""" + T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) + +Main mesh constructor for the `T8codeMesh` that imports an unstructured, +conforming mesh from a `p4est_connectivity` data structure. + +# Arguments +- `conn::Ptr{p4est_connectivity}`: Pointer to a P4est connectivity object. +- `mapping`: a function of `NDIMS` variables to describe the mapping that transforms + the imported mesh to the physical domain. Use `nothing` for the identity map. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. + The mapping will be approximated by an interpolation polynomial + of the specified degree for each tree. + The default of `1` creates an uncurved geometry. Use a higher value if the mapping + will curve the imported uncurved mesh. +- `RealT::Type`: the type that should be used for coordinates. +- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +""" +function T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) + cmesh = t8_cmesh_new_from_p8est(conn, mpi_comm(), 0) + + return T8codeMesh(cmesh; kwargs...) end """ - T8codeMesh{NDIMS}(meshfile::String; - mapping=nothing, polydeg=1, RealT=Float64, - initial_refinement_level=0) + T8codeMesh{NDIMS}(meshfile::String; kwargs...) Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming mesh from a Gmsh mesh file (`.msh`). # Arguments - `meshfile::String`: path to a Gmsh mesh file. +- `ndims`: Mesh file dimension: `2` or `3`. - `mapping`: a function of `NDIMS` variables to describe the mapping that transforms the imported mesh to the physical domain. Use `nothing` for the identity map. - `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. @@ -321,17 +460,130 @@ mesh from a Gmsh mesh file (`.msh`). - `RealT::Type`: the type that should be used for coordinates. - `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. """ -function T8codeMesh{NDIMS}(meshfile::String; kwargs...) where {NDIMS} - @assert NDIMS == 2 # Only support for NDIMS = 2 yet. +function T8codeMesh(meshfile::String, ndims; kwargs...) # Prevent `t8code` from crashing Julia if the file doesn't exist. @assert isfile(meshfile) meshfile_prefix, meshfile_suffix = splitext(meshfile) - cmesh = t8_cmesh_from_msh_file(meshfile_prefix, 0, mpi_comm(), NDIMS, 0, 0) + cmesh = t8_cmesh_from_msh_file(meshfile_prefix, 0, mpi_comm(), ndims, 0, 0) + + return T8codeMesh(cmesh; kwargs...) +end + +struct adapt_callback_passthrough + adapt_callback::Function + user_data::Any +end - return T8codeMesh{NDIMS}(cmesh; kwargs...) +# Callback function prototype to decide for refining and coarsening. +# If `is_family` equals 1, the first `num_elements` in elements +# form a family and we decide whether this family should be coarsened +# or only the first element should be refined. +# Otherwise `is_family` must equal zero and we consider the first entry +# of the element array for refinement. +# Entries of the element array beyond the first `num_elements` are undefined. +# \param [in] forest the forest to which the new elements belong +# \param [in] forest_from the forest that is adapted. +# \param [in] which_tree the local tree containing `elements` +# \param [in] lelement_id the local element id in `forest_old` in the tree of the current element +# \param [in] ts the eclass scheme of the tree +# \param [in] is_family if 1, the first `num_elements` entries in `elements` form a family. If 0, they do not. +# \param [in] num_elements the number of entries in `elements` that are defined +# \param [in] elements Pointers to a family or, if `is_family` is zero, +# pointer to one element. +# \return greater zero if the first entry in `elements` should be refined, +# smaller zero if the family `elements` shall be coarsened, +# zero else. +function adapt_callback_wrapper(forest, + forest_from, + which_tree, + lelement_id, + ts, + is_family, + num_elements, + elements_ptr)::Cint + passthrough = unsafe_pointer_to_objref(t8_forest_get_user_data(forest))[] + + elements = unsafe_wrap(Array, elements_ptr, num_elements) + + return passthrough.adapt_callback(forest_from, which_tree, ts, lelement_id, elements, + Bool(is_family), passthrough.user_data) +end + +""" + Trixi.adapt!(mesh::T8codeMesh, adapt_callback; kwargs...) + +Adapt a `T8codeMesh` according to a user-defined `adapt_callback`. + +# Arguments +- `mesh::T8codeMesh`: Initialized mesh object. +- `adapt_callback`: A user-defined callback which tells the adaption routines + if an element should be refined, coarsened or stay unchanged. + + The expected callback signature is as follows: + + `adapt_callback(forest, ltreeid, eclass_scheme, lelemntid, elements, is_family, user_data)` + # Arguments + - `forest`: Pointer to the analyzed forest. + - `ltreeid`: Local index of the current tree where the analyzed elements are part of. + - `eclass_scheme`: Element class of `elements`. + - `lelemntid`: Local index of the first element in `elements`. + - `elements`: Array of elements. If consecutive elements form a family + they are passed together, otherwise `elements` consists of just one element. + - `is_family`: Boolean signifying if `elements` represents a family or not. + - `user_data`: Void pointer to some arbitrary user data. Default value is `C_NULL`. + # Returns + -1 : Coarsen family of elements. + 0 : Stay unchanged. + 1 : Refine element. + +- `kwargs`: + - `recursive = true`: Adapt the forest recursively. If true the caller must ensure that the callback + returns 0 for every analyzed element at some point to stop the recursion. + - `balance = true`: Make sure the adapted forest is 2^(NDIMS-1):1 balanced. + - `partition = true`: Partition the forest to redistribute elements evenly among MPI ranks. + - `ghost = true`: Create a ghost layer for MPI data exchange. + - `user_data = C_NULL`: Pointer to some arbitrary user-defined data. +""" +function adapt!(mesh::T8codeMesh, adapt_callback; recursive = true, balance = true, + partition = true, ghost = true, user_data = C_NULL) + # Check that forest is a committed, that is valid and usable, forest. + @assert t8_forest_is_committed(mesh.forest) != 0 + + # Init new forest. + new_forest_ref = Ref{t8_forest_t}() + t8_forest_init(new_forest_ref) + new_forest = new_forest_ref[] + + # Check out `examples/t8_step4_partition_balance_ghost.jl` in + # https://github.com/DLR-AMR/T8code.jl for detailed explanations. + let set_from = C_NULL, set_for_coarsening = 0, no_repartition = !partition + t8_forest_set_user_data(new_forest, + pointer_from_objref(Ref(adapt_callback_passthrough(adapt_callback, + user_data)))) + t8_forest_set_adapt(new_forest, mesh.forest, + @t8_adapt_callback(adapt_callback_wrapper), + recursive) + if balance + t8_forest_set_balance(new_forest, set_from, no_repartition) + end + + if partition + t8_forest_set_partition(new_forest, set_from, set_for_coarsening) + end + + t8_forest_set_ghost(new_forest, ghost, T8_GHOST_FACES) # Note: MPI support not available yet so it is a dummy call. + + # The old forest is destroyed here. + # Call `t8_forest_ref(Ref(mesh.forest))` to keep it. + t8_forest_commit(new_forest) + end + + mesh.forest = new_forest + + return nothing end # TODO: Just a placeholder. Will be implemented later when MPI is supported. @@ -343,3 +595,10 @@ end function partition!(mesh::T8codeMesh; allow_coarsening = true, weight_fn = C_NULL) return nothing end + +#! format: off +@deprecate T8codeMesh{2}(conn::Ptr{p4est_connectivity}; kwargs...) T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...) +@deprecate T8codeMesh{3}(conn::Ptr{p8est_connectivity}; kwargs...) T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) +@deprecate T8codeMesh{2}(meshfile::String; kwargs...) T8codeMesh(meshfile::String, 2; kwargs...) +@deprecate T8codeMesh{3}(meshfile::String; kwargs...) T8codeMesh(meshfile::String, 3; kwargs...) +#! format: on diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index fe7858e31ee..8518cf27fd3 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -253,6 +253,9 @@ end function _jacobian_ad_forward(semi, t0, u0_ode, du_ode, config) new_semi = remake(semi, uEltype = eltype(config)) + # Create anonymous function passed as first argument to `ForwardDiff.jacobian` to match + # `ForwardDiff.jacobian(f!, y::AbstractArray, x::AbstractArray, + # cfg::JacobianConfig = JacobianConfig(f!, y, x), check=Val{true}())` J = ForwardDiff.jacobian(du_ode, u0_ode, config) do du_ode, u_ode Trixi.rhs!(du_ode, u_ode, new_semi, t0) end @@ -279,6 +282,9 @@ end function _jacobian_ad_forward_structarrays(semi, t0, u0_ode_plain, du_ode_plain, config) new_semi = remake(semi, uEltype = eltype(config)) + # Create anonymous function passed as first argument to `ForwardDiff.jacobian` to match + # `ForwardDiff.jacobian(f!, y::AbstractArray, x::AbstractArray, + # cfg::JacobianConfig = JacobianConfig(f!, y, x), check=Val{true}())` J = ForwardDiff.jacobian(du_ode_plain, u0_ode_plain, config) do du_ode_plain, u_ode_plain u_ode = StructArray{SVector{nvariables(semi), eltype(config)}}(ntuple(v -> view(u_ode_plain, diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index e9994fe4569..7e383924ba7 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -6,7 +6,8 @@ #! format: noindent # Initialize data structures in element container -function init_elements!(elements, mesh::P4estMesh{3}, basis::LobattoLegendreBasis) +function init_elements!(elements, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + basis::LobattoLegendreBasis) @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -26,7 +27,7 @@ end # Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis function calc_node_coordinates!(node_coordinates, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, basis::LobattoLegendreBasis) # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index 4c0845ba9af..5b3c5ae5ca8 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.jl @@ -7,8 +7,8 @@ # The methods below are specialized on the mortar type # and called from the basic `create_cache` method at the top. -function create_cache(mesh::P4estMesh{3}, equations, mortar_l2::LobattoLegendreMortarL2, - uEltype) +function create_cache(mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, + mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal compare performance of different types fstar_threaded = [Array{uEltype, 4}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2), 4) @@ -88,7 +88,7 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral, dg::DG) @unpack interfaces = cache index_range = eachnode(dg) @@ -163,7 +163,7 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -244,7 +244,7 @@ end # Inlined function for interface flux computation for conservative flux terms @inline function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -271,7 +271,7 @@ end # Inlined function for interface flux computation for flux + nonconservative terms @inline function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -314,7 +314,7 @@ end end function prolong2boundaries!(cache, u, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral, dg::DG) @unpack boundaries = cache index_range = eachnode(dg) @@ -355,7 +355,7 @@ function prolong2boundaries!(cache, u, end function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral, dg::DG) @unpack boundaries = cache @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements @@ -417,7 +417,7 @@ function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing, end function prolong2mortars!(cache, u, - mesh::P4estMesh{3}, equations, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DGSEM) @unpack fstar_tmp_threaded = cache @@ -521,7 +521,7 @@ function prolong2mortars!(cache, u, end function calc_mortar_flux!(surface_flux_values, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -595,7 +595,7 @@ end # Inlined version of the mortar flux computation on small elements for conservation fluxes @inline function calc_mortar_flux!(fstar, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -616,7 +616,7 @@ end # Inlined version of the mortar flux computation on small elements for conservation fluxes # with nonconservative terms @inline function calc_mortar_flux!(fstar, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, mortar_index, position_index, normal_direction, @@ -643,7 +643,8 @@ end end @inline function mortar_fluxes_to_elements!(surface_flux_values, - mesh::P4estMesh{3}, equations, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, mortar, fstar, u_buffer, fstar_tmp) @@ -727,7 +728,7 @@ end end function calc_surface_integral!(du, u, - mesh::P4estMesh{3}, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index cdb085e9008..1df9f408895 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -58,7 +58,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 =# @inline function weak_form_kernel!(du, u, element, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -115,7 +116,8 @@ end # the physical fluxes in each Cartesian direction @inline function flux_differencing_kernel!(du, u, element, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::False, equations, volume_flux, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -189,7 +191,8 @@ end @inline function flux_differencing_kernel!(du, u, element, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::True, equations, volume_flux, dg::DGSEM, cache, alpha = true) @unpack derivative_split = dg.basis @@ -274,7 +277,8 @@ end # [arXiv: 2008.12044v2](https://arxiv.org/pdf/2008.12044) @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::False, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @@ -369,7 +373,8 @@ end # # Calculate the finite volume fluxes inside curvilinear elements (**with non-conservative terms**). @inline function calcflux_fv!(fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, u, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms::True, equations, volume_flux_fv, dg::DGSEM, element, cache) @unpack contravariant_vectors = cache.elements @@ -783,7 +788,7 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, end function apply_jacobian!(du, - mesh::Union{StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, dg::DG, cache) @threaded for element in eachelement(dg, cache) for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) diff --git a/src/solvers/dgsem_t8code/containers_2d.jl b/src/solvers/dgsem_t8code/containers_2d.jl index 029e6674afb..bf77826a34b 100644 --- a/src/solvers/dgsem_t8code/containers_2d.jl +++ b/src/solvers/dgsem_t8code/containers_2d.jl @@ -1,3 +1,7 @@ +# 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 @@ -27,6 +31,9 @@ function calc_node_coordinates!(node_coordinates, element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) element_level = t8_element_level(eclass_scheme, element) + # Note, `t8_quad_len` is encoded as an integer (Morton encoding) in + # relation to `t8_quad_root_len`. This line transforms the + # "integer" length to a float in relation to the unit interval [0,1]. element_length = t8_quad_len(element_level) / t8_quad_root_len element_coords = Array{Float64}(undef, 3) @@ -55,4 +62,16 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end + +function init_mortar_neighbor_ids!(mortars::P4estMortarContainer{2}, my_face, + other_face, orientation, neighbor_ielements, + mortar_id) + if orientation == 0 + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1 + else + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1 + end +end end # @muladd diff --git a/src/solvers/dgsem_t8code/containers_3d.jl b/src/solvers/dgsem_t8code/containers_3d.jl new file mode 100644 index 00000000000..f2d54ff07da --- /dev/null +++ b/src/solvers/dgsem_t8code/containers_3d.jl @@ -0,0 +1,235 @@ +# 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 + +# Interpolate tree_node_coordinates to each quadrant at the specified nodes +function calc_node_coordinates!(node_coordinates, + mesh::T8codeMesh{3}, + nodes::AbstractVector) + # We use `StrideArray`s here since these buffers are used in performance-critical + # places and the additional information passed to the compiler makes them faster + # than native `Array`s. + tmp1 = StrideArray(undef, real(mesh), + StaticInt(3), static_length(nodes), static_length(mesh.nodes), + static_length(mesh.nodes)) + matrix1 = StrideArray(undef, real(mesh), + static_length(nodes), static_length(mesh.nodes)) + matrix2 = similar(matrix1) + matrix3 = similar(matrix1) + baryweights_in = barycentric_weights(mesh.nodes) + + num_local_trees = t8_forest_get_num_local_trees(mesh.forest) + + current_index = 0 + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(mesh.forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class) + num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree) + + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) + element_level = t8_element_level(eclass_scheme, element) + + # Note, `t8_hex_len` is encoded as an integer (Morton encoding) in + # relation to `t8_hex_root_len`. This line transforms the + # "integer" length to a float in relation to the unit interval [0,1]. + element_length = t8_hex_len(element_level) / t8_hex_root_len + + element_coords = Vector{Float64}(undef, 3) + t8_element_vertex_reference_coords(eclass_scheme, element, 0, + pointer(element_coords)) + + nodes_out_x = (2 * + (element_length * 0.5 * (nodes .+ 1) .+ element_coords[1]) .- + 1) + nodes_out_y = (2 * + (element_length * 0.5 * (nodes .+ 1) .+ element_coords[2]) .- + 1) + nodes_out_z = (2 * + (element_length * 0.5 * (nodes .+ 1) .+ element_coords[3]) .- + 1) + + polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, + baryweights_in) + polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, + baryweights_in) + polynomial_interpolation_matrix!(matrix3, mesh.nodes, nodes_out_z, + baryweights_in) + + multiply_dimensionwise!(view(node_coordinates, :, :, :, :, + current_index += 1), + matrix1, matrix2, matrix3, + view(mesh.tree_node_coordinates, :, :, :, :, + itree + 1), + tmp1) + end + end + + return node_coordinates +end + +# This routine was copied and adapted from `src/dgsem_p4est/containers_3d.jl`: `orientation_to_indices_p4est`. +function init_mortar_neighbor_ids!(mortars::P4estMortarContainer{3}, my_face, + other_face, orientation, neighbor_ielements, + mortar_id) + # my_face and other_face are the face directions (zero-based) + # of "my side" and "other side" respectively. + # Face corner 0 of the face with the lower face direction connects to a corner of the other face. + # The number of this corner is the orientation code in `p4est`. + lower = my_face <= other_face + + # x_pos, y_neg, and z_pos are the directions in which the face has right-handed coordinates + # when looked at from the outside. + my_right_handed = my_face in (1, 2, 5) + other_right_handed = other_face in (1, 2, 5) + + # If both or none are right-handed when looked at from the outside, they will have different + # orientations when looked at from the same side of the interface. + flipped = my_right_handed == other_right_handed + + # In the following illustrations, the face corner numbering of `p4est` is shown. + # ξ and η are the local coordinates of the respective face. + # We're looking at both faces from the same side of the interface, so that "other side" + # (in the illustrations on the left) has right-handed coordinates. + if !flipped + if orientation == 0 + # Corner 0 of other side matches corner 0 of my side + # 2┌──────┐3 2┌──────┐3 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 0└──────┘1 + # η η + # ↑ ↑ + # │ │ + # └───> ξ └───> ξ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[4] + 1 + + elseif ((lower && orientation == 2) # Corner 0 of my side matches corner 2 of other side + || + (!lower && orientation == 1)) # Corner 0 of other side matches corner 1 of my side + # 2┌──────┐3 0┌──────┐2 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 1└──────┘3 + # η ┌───> η + # ↑ │ + # │ ↓ + # └───> ξ ξ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[3] + 1 + + elseif ((lower && orientation == 1) # Corner 0 of my side matches corner 1 of other side + || + (!lower && orientation == 2)) # Corner 0 of other side matches corner 2 of my side + # 2┌──────┐3 3┌──────┐1 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 2└──────┘0 + # η ξ + # ↑ ↑ + # │ │ + # └───> ξ η <───┘ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[2] + 1 + + else # orientation == 3 + # Corner 0 of my side matches corner 3 of other side and + # corner 0 of other side matches corner 3 of my side. + # 2┌──────┐3 1┌──────┐0 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 3└──────┘2 + # η ξ <───┐ + # ↑ │ + # │ ↓ + # └───> ξ η + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[1] + 1 + end + else # flipped + if orientation == 0 + # Corner 0 of other side matches corner 0 of my side + # 2┌──────┐3 1┌──────┐3 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 0└──────┘2 + # η ξ + # ↑ ↑ + # │ │ + # └───> ξ └───> η + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[4] + 1 + + elseif orientation == 2 + # Corner 0 of my side matches corner 2 of other side and + # corner 0 of other side matches corner 2 of my side. + # 2┌──────┐3 0┌──────┐1 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 2└──────┘3 + # η ┌───> ξ + # ↑ │ + # │ ↓ + # └───> ξ η + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[2] + 1 + + elseif orientation == 1 + # Corner 0 of my side matches corner 1 of other side and + # corner 0 of other side matches corner 1 of my side. + # 2┌──────┐3 3┌──────┐2 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 1└──────┘0 + # η η + # ↑ ↑ + # │ │ + # └───> ξ ξ <───┘ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[1] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[3] + 1 + + else # orientation == 3 + # Corner 0 of my side matches corner 3 of other side and + # corner 0 of other side matches corner 3 of my side. + # 2┌──────┐3 2┌──────┐0 + # │ │ │ │ + # │ │ │ │ + # 0└──────┘1 3└──────┘1 + # η η <───┐ + # ↑ │ + # │ ↓ + # └───> ξ ξ + + mortars.neighbor_ids[1, mortar_id] = neighbor_ielements[4] + 1 + mortars.neighbor_ids[2, mortar_id] = neighbor_ielements[2] + 1 + mortars.neighbor_ids[3, mortar_id] = neighbor_ielements[3] + 1 + mortars.neighbor_ids[4, mortar_id] = neighbor_ielements[1] + 1 + end + end +end +end # @muladd diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index 16a9d7d35b1..6e9660c917d 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -28,4 +28,5 @@ end include("containers.jl") include("containers_2d.jl") +include("containers_3d.jl") end # @muladd diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 0955dc38655..02ff338e912 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -36,13 +36,15 @@ end # The methods below are specialized on the volume integral type # and called from the basic `create_cache` method at the top. -function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype) NamedTuple() end -function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] @@ -79,8 +81,8 @@ function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, fstar2_L_threaded, fstar2_R_threaded, fstar3_L_threaded, fstar3_R_threaded) end -function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, - equations, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG, uEltype) A4dp1_x = Array{uEltype, 4} @@ -112,7 +114,8 @@ end # The methods below are specialized on the mortar type # and called from the basic `create_cache` method at the top. -function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, +function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal compare performance of different types A3d = Array{uEltype, 3} @@ -140,7 +143,7 @@ end # TODO: Taal discuss/refactor timer, allowing users to pass a custom timer? function rhs!(du, u, t, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, + mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -209,8 +212,8 @@ function rhs!(du, u, t, end function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, - P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -264,8 +267,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 end function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, - P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) @@ -378,8 +381,8 @@ end # TODO: Taal dimension agnostic function calc_volume_integral!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, - P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DGSEM, cache) @@ -437,7 +440,8 @@ function calc_volume_integral!(du, u, end @inline function fv_kernel!(du, u, - mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}}, + mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, + T8codeMesh{3}}, nonconservative_terms, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded, fstar3_L_threaded, fstar3_R_threaded = cache diff --git a/src/solvers/dgsem_tree/indicators_3d.jl b/src/solvers/dgsem_tree/indicators_3d.jl index 40362889397..a11a8e06e4b 100644 --- a/src/solvers/dgsem_tree/indicators_3d.jl +++ b/src/solvers/dgsem_tree/indicators_3d.jl @@ -101,7 +101,8 @@ end alpha[element] = min(alpha_max, alpha_element) end -function apply_smoothing!(mesh::Union{TreeMesh{3}, P4estMesh{3}}, alpha, alpha_tmp, dg, +function apply_smoothing!(mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, alpha, + alpha_tmp, dg, cache) # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha diff --git a/test/runtests.jl b/test/runtests.jl index 7e195fe7402..49f0977bb70 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -89,6 +89,10 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) include("test_t8code_2d.jl") end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "t8code_part2" + include("test_t8code_3d.jl") + end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "unstructured_dgmulti" include("test_unstructured_2d.jl") include("test_dgmulti_1d.jl") diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 4a2d2112c99..dc5d32b5a04 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -380,18 +380,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_circular_wind_nonconforming.jl"), l2=[ - 1.573832094977477e-7, - 3.863090659429634e-5, - 3.867293305754584e-5, - 3.686550296950078e-5, - 0.05508968493733932, + 1.5737711609657832e-7, + 3.8630261900166194e-5, + 3.8672287531936816e-5, + 3.6865116098660796e-5, + 0.05508620970403884, ], linf=[ - 2.2695202613887133e-6, - 0.0005314968179916946, - 0.0005314969614147458, - 0.0005130280733059617, - 0.7944959432352334, + 2.268845333053271e-6, + 0.000531462302113539, + 0.0005314624461298934, + 0.0005129931254772464, + 0.7942778058932163, ], tspan=(0.0, 2e2), coverage_override=(trees_per_cube_face = (1, 1), polydeg = 3)) # Prevent long compile time in CI @@ -409,18 +409,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_baroclinic_instability.jl"), l2=[ - 6.725065410642336e-7, - 0.00021710117340245454, - 0.000438679759422352, - 0.00020836356588024185, - 0.07602006689579247, + 6.725093801700048e-7, + 0.00021710076010951073, + 0.0004386796338203878, + 0.00020836270267103122, + 0.07601887903440395, ], linf=[ - 1.9101671995258585e-5, - 0.029803626911022396, - 0.04847630924006063, - 0.022001371349740104, - 4.847761006938526, + 1.9107530539574924e-5, + 0.02980358831035801, + 0.048476331898047564, + 0.02200137344113612, + 4.848310144356219, ], tspan=(0.0, 1e2), # Decrease tolerance of adaptive time stepping to get similar results across different systems diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 1addc29e3e6..e5d45ebcc07 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -611,6 +611,33 @@ end end end +@trixi_testset "elixir_euler_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_warm_bubble.jl"), + l2=[ + 0.00019387402388722496, + 0.03086514388623955, + 0.04541427917165, + 43.892826583444716, + ], + linf=[ + 0.0015942305974430138, + 0.17449778969139373, + 0.3729704262394843, + 307.6706958565337, + ], + cells_per_dimension=(32, 16), + tspan=(0.0, 10.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)) < 100 + end +end + @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), l2=[ diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index b3e19471323..ab95e068d02 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -30,7 +30,21 @@ mkdir(outdir) end end +@trixi_testset "test check_for_negative_volumes" begin + @test_warn "Discovered negative volumes" begin + # Unstructured mesh with six cells which have left-handed node ordering. + mesh_file = joinpath(EXAMPLES_DIR, "rectangle_with_negative_volumes.msh") + isfile(mesh_file) || + download("https://gist.githubusercontent.com/jmark/bfe0d45f8e369298d6cc637733819013/raw/cecf86edecc736e8b3e06e354c494b2052d41f7a/rectangle_with_negative_volumes.msh", + mesh_file) + + # This call should throw a warning about negative volumes detected. + mesh = T8codeMesh(mesh_file, 2) + end +end + @trixi_testset "elixir_advection_basic.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), # Expected errors are exactly the same as with TreeMesh! l2=[8.311947673061856e-6], @@ -46,6 +60,7 @@ end end @trixi_testset "elixir_advection_nonconforming_flag.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming_flag.jl"), l2=[3.198940059144588e-5], @@ -61,6 +76,7 @@ end end @trixi_testset "elixir_advection_unstructured_flag.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_flag.jl"), l2=[0.0005379687442422346], linf=[0.007438525029884735]) @@ -91,6 +107,7 @@ end end @trixi_testset "elixir_advection_amr_solution_independent.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_solution_independent.jl"), # Expected errors are exactly the same as with StructuredMesh! @@ -108,6 +125,7 @@ end end @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), l2=[ @@ -133,6 +151,7 @@ end end @trixi_testset "elixir_euler_free_stream.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), l2=[ 2.063350241405049e-15, @@ -153,6 +172,7 @@ end end @trixi_testset "elixir_euler_shockcapturing_ec.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_ec.jl"), l2=[ 9.53984675e-02, @@ -178,6 +198,8 @@ end end @trixi_testset "elixir_euler_sedov.jl" begin + # This test is identical to the one in `test_p4est_2d.jl` besides minor + # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[ 3.76149952e-01, @@ -203,6 +225,7 @@ end end @trixi_testset "elixir_shallowwater_source_terms.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ 9.168126407325352e-5, @@ -228,6 +251,7 @@ end end @trixi_testset "elixir_mhd_alfven_wave.jl" begin + # This test is identical to the one in `test_p4est_2d.jl`. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), l2=[1.0513414461545583e-5, 1.0517900957166411e-6, 1.0517900957304043e-6, 1.511816606372376e-6, @@ -250,6 +274,8 @@ end end @trixi_testset "elixir_mhd_rotor.jl" begin + # This test is identical to the one in `test_p4est_2d.jl` besides minor + # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), l2=[0.44211360369891683, 0.8805178316216257, 0.8262710688468049, 0.0, diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl new file mode 100644 index 00000000000..4232cf04094 --- /dev/null +++ b/test/test_t8code_3d.jl @@ -0,0 +1,279 @@ +module TestExamplesT8codeMesh3D + +using Test +using Trixi + +include("test_trixi.jl") + +EXAMPLES_DIR = joinpath(examples_dir(), "t8code_3d_dgsem") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) +mkdir(outdir) + +@testset "T8codeMesh3D" begin + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_advection_basic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[0.00016263963870641478], + linf=[0.0014537194925779984]) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_advection_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_unstructured_curved.jl"), + l2=[0.0004750004258546538], + linf=[0.026527551737137167]) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_advection_nonconforming.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"), + l2=[0.00253595715323843], + linf=[0.016486952252155795]) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl` besides minor + # deviations from the expected error norms. + @trixi_testset "elixir_advection_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"), + # Expected errors are exactly the same as with TreeMesh! + l2=[1.1302812803902801e-5], + linf=[0.0007889950196294793], + coverage_override=(maxiters = 6, initial_refinement_level = 1, + base_level = 1, med_level = 2, max_level = 3)) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl` besides minor + # deviations from the expected error norms. + @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_amr_unstructured_curved.jl"), + l2=[2.0556575425846923e-5], + linf=[0.00105682693484822], + tspan=(0.0, 1.0), + coverage_override=(maxiters = 6, initial_refinement_level = 0, + base_level = 0, med_level = 1, max_level = 2)) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_curved.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonconforming_unstructured_curved.jl"), + l2=[ + 4.070355207909268e-5, + 4.4993257426833716e-5, + 5.10588457841744e-5, + 5.102840924036687e-5, + 0.00019986264001630542, + ], + linf=[ + 0.0016987332417202072, + 0.003622956808262634, + 0.002029576258317789, + 0.0024206977281964193, + 0.008526972236273522, + ], + tspan=(0.0, 0.01)) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonperiodic.jl"), + l2=[ + 0.0015106060984283647, + 0.0014733349038567685, + 0.00147333490385685, + 0.001473334903856929, + 0.0028149479453087093, + ], + linf=[ + 0.008070806335238156, + 0.009007245083113125, + 0.009007245083121784, + 0.009007245083102688, + 0.01562861968368434, + ], + tspan=(0.0, 1.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 + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_free_stream.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"), + l2=[ + 5.162664597942288e-15, + 1.941857343642486e-14, + 2.0232366394187278e-14, + 2.3381518645408552e-14, + 7.083114561232324e-14, + ], + linf=[ + 7.269740365245525e-13, + 3.289868377720495e-12, + 4.440087186807773e-12, + 3.8686831516088205e-12, + 9.412914891981927e-12, + ], + tspan=(0.0, 0.03)) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_free_stream_extruded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"), + l2=[ + 8.444868392439035e-16, + 4.889826056731442e-15, + 2.2921260987087585e-15, + 4.268460455702414e-15, + 1.1356712092620279e-14, + ], + linf=[ + 7.749356711883593e-14, + 2.8792246364872653e-13, + 1.1121659149182506e-13, + 3.3228975127030935e-13, + 9.592326932761353e-13, + ], + tspan=(0.0, 0.1)) + # 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 + + # This test is identical to the one in `test_p4est_3d.jl`. + @trixi_testset "elixir_euler_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), + l2=[ + 0.010380390326164493, + 0.006192950051354618, + 0.005970674274073704, + 0.005965831290564327, + 0.02628875593094754, + ], + linf=[ + 0.3326911600075694, + 0.2824952141320467, + 0.41401037398065543, + 0.45574161423218573, + 0.8099577682187109, + ], + tspan=(0.0, 0.2), + coverage_override=(polydeg = 3,)) # Prevent long compile time in CI + # 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 + + # This test is identical to the one in `test_p4est_3d.jl` besides minor + # deviations in the expected error norms. + @trixi_testset "elixir_euler_sedov.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), + l2=[ + 7.82070951e-02, + 4.33260474e-02, + 4.33260474e-02, + 4.33260474e-02, + 3.75260911e-01, + ], + linf=[ + 7.45329845e-01, + 3.21754792e-01, + 3.21754792e-01, + 3.21754792e-01, + 4.76151527e+00, + ], + tspan=(0.0, 0.3), + coverage_override=(polydeg = 3,)) # Prevent long compile time in CI + # 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 +@test_nowarn rm(outdir, recursive = true) + +end # module diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 65899cd5263..61b5c54b5e9 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -834,6 +834,32 @@ end end end +@trixi_testset "elixir_euler_warm_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_warm_bubble.jl"), + l2=[ + 0.0001379946769624388, + 0.02078779689715382, + 0.033237241571263176, + 31.36068872331705, + ], + linf=[ + 0.0016286690573188434, + 0.15623770697198225, + 0.3341371832270615, + 334.5373488726036, + ], + tspan=(0.0, 10.0), + initial_refinement_level=4) + # 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)) < 100 + end +end + # Coverage test for all initial conditions @testset "Compressible Euler: Tests for initial conditions" begin @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 58db7c5f35f..1f3dfbf5267 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -327,6 +327,31 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_shallowwater_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), + l2=[ + 0.13517233723296504, + 0.20010876311162215, + 0.20010876311162223, + 2.719538414346464e-7, + ], + linf=[ + 0.5303607982988336, + 0.5080989745682338, + 0.5080989745682352, + 1.1301675764130437e-6, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index 817b4cd550d..e8a8effbe29 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1287,6 +1287,49 @@ end end end +@timed_testset "Consistency check for LMARS flux" begin + equations = CompressibleEulerEquations2D(1.4) + flux_lmars = FluxLMARS(340) + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + orientations = [1, 2] + u_values = [SVector(1.0, 0.5, -0.7, 1.0), + SVector(1.5, -0.2, 0.1, 5.0)] + + for u in u_values, orientation in orientations + @test flux_lmars(u, u, orientation, equations) ≈ + flux(u, orientation, equations) + end + + for u in u_values, normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + equations = CompressibleEulerEquations3D(1.4) + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + orientations = [1, 2, 3] + u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), + SVector(1.5, -0.2, 0.1, 0.2, 5.0)] + + for u in u_values, orientation in orientations + @test flux_lmars(u, u, orientation, equations) ≈ + flux(u, orientation, equations) + end + + for u in u_values, normal_direction in normal_directions + @test flux_lmars(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end +end + @testset "FluxRotated vs. direct implementation" begin @timed_testset "CompressibleEulerMulticomponentEquations2D" begin equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4), @@ -1320,7 +1363,8 @@ end u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc, + ] for f_std in fluxes f_rot = FluxRotated(f_std)