Merge branch 'main' into t8codemesh-FV-1order

trixi-framework · Jan 19, 2024 · ee713ee · ee713ee
2 parents 51c5d25 + 1946f9d
commit ee713ee
Show file tree

Hide file tree

Showing 60 changed files with 2,753 additions and 382 deletions.
diff --git 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
@@ -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
@@ -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
@@ -10,6 +10,7 @@
 *.mesh
 *.bson
 *.inp
+*.msh
 **/Manifest.toml
 out*/
 docs/build

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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,
-Original file line number
+Diff line change
@@ Expand Up / @@ -10,6 +10,7 @@ @@
     *.mesh
     *.bson
     *.inp
+    *.msh
     **/Manifest.toml
     out*/
     docs/build
@@ Expand Down @@