Merge branch 'main' into subcell-limiting

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/[email protected].0
+        uses: crate-ci/[email protected].2

NEWS.md

@@ -4,13 +4,20 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes in the v0.7 lifecycle
+
+#### Added
+- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension
+  to 1D and 3D on `TreeMesh`.
+
+
 ## Changes when updating to v0.7 from v0.6.x
 
 #### Added
 
 #### Changed
 
-- The default wave speed estimate used within `flux_hll` is now `min_max_speed_davis` 
+- The default wave speed estimate used within `flux_hll` is now `min_max_speed_davis`
   instead of `min_max_speed_naive`.
 
 #### Deprecated
@@ -26,7 +33,7 @@ for human readability.
 #### Added
 - AMR for hyperbolic-parabolic equations on 3D `P4estMesh`
 - `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D`
-- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh, 
+- Different boundary conditions for quad/hex meshes in Abaqus format, even if not generated by HOHQMesh,
   can now be digested by Trixi in 2D and 3D.
 - Subcell (positivity) limiting support for nonlinear variables in 2D for `TreeMesh`
 - Added Lighthill-Whitham-Richards (LWR) traffic model
@@ -40,7 +47,7 @@ for human readability.
 #### Changed
 
 - The wave speed estimates for `flux_hll`, `FluxHLL()` are now consistent across equations.
-  In particular, the functions `min_max_speed_naive`, `min_max_speed_einfeldt` are now 
+  In particular, the functions `min_max_speed_naive`, `min_max_speed_einfeldt` are now
   conceptually identical across equations.
   Users, who have been using `flux_hll` for MHD have now to use `flux_hlle` in order to use the
   Einfeldt wave speed estimate.

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.7.1-pre"
+version = "0.7.4-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

docs/literate/src/files/scalar_linear_advection_1d.jl

@@ -115,7 +115,7 @@ integral = sum(nodes.^3 .* weights)
 # To approximate the solution, we need to get the polynomial coefficients $\{u_j^{Q_l}\}_{j=0}^N$
 # for every element $Q_l$.
 
-# After defining all nodes, we can implement the spatial coordinate $x$ and its initial value $u0$
+# After defining all nodes, we can implement the spatial coordinate $x$ and its initial value $u0 = u(t_0)$
 # for every node.
 x = Matrix{Float64}(undef, length(nodes), n_elements)
 for element in 1:n_elements

docs/src/meshes/p4est_mesh.md

@@ -256,6 +256,8 @@ By doing so, only nodesets with a label present in `boundary_symbols` are treate
 Other nodesets that could be used for diagnostics are not treated as external boundaries.
 Note that there is a leading colon `:` compared to the label in the `.inp` mesh file.
 This is required to turn the label into a [`Symbol`](https://docs.julialang.org/en/v1/manual/metaprogramming/#Symbols).
+**Important**: In Julia, a symbol _cannot_ contain a hyphen/dash `-`, i.e., `:BC-1` is _not_ a valid symbol.
+Keep this in mind when importing boundaries, you might have to convert hyphens/dashes `-` to underscores `_` in the `.inp` mesh file, i.e., `BC_1` instead of `BC-1`.
 
 A 2D example for this mesh, which is read-in for an unstructured mesh file created with `gmsh`, is presented in 
 `examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl`.

docs/src/parallelization.md

@@ -69,8 +69,7 @@ installations. Follow the steps described [here](https://github.com/DLR-AMR/T8co
 [here](https://github.com/trixi-framework/P4est.jl/blob/main/README.md#installation) for the configuration.
 The paths that point to `libp4est.so` (and potentially to `libsc.so`) need to be
 the same for P4est.jl and T8code.jl. This could, e.g., be `libp4est.so` that usually can be found
-in `lib/` or `local/lib/` in the installation directory of `t8code`. Note that the `T8codeMesh`, however,
-does not support MPI yet.
+in `lib/` or `local/lib/` in the installation directory of `t8code`.
 The preferences for [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) always need to be set, even if you
 do not want to use `HDF5` from Trixi.jl, see also [issue #1079 in HDF5.jl](https://github.com/JuliaIO/HDF5.jl/issues/1079).
 To set the preferences for HDF5.jl, follow the instructions described

examples/structured_2d_dgsem/elixir_advection_coupled.jl

@@ -53,7 +53,8 @@ cells_per_dimension = (8, 8)
 coordinates_min1 = (-1.0, 0.0) # minimum coordinates (min(x), min(y))
 coordinates_max1 = (0.0, 1.0) # maximum coordinates (max(x), max(y))
 
-mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1)
+mesh1 = StructuredMesh(cells_per_dimension, coordinates_min1, coordinates_max1,
+                       periodicity = false)
 
 # Define the coupling functions
 coupling_function12 = (x, u, equations_other, equations_own) -> u
@@ -84,7 +85,8 @@ semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_converg
 coordinates_min2 = (0.0, 0.0) # minimum coordinates (min(x), min(y))
 coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y))
 
-mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2)
+mesh2 = StructuredMesh(cells_per_dimension, coordinates_min2, coordinates_max2,
+                       periodicity = false)
 
 # Define the coupling functions
 coupling_function21 = (x, u, equations_other, equations_own) -> u
@@ -115,7 +117,8 @@ semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_converg
 coordinates_min3 = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
 coordinates_max3 = (0.0, 0.0) # maximum coordinates (max(x), max(y))
 
-mesh3 = StructuredMesh(cells_per_dimension, coordinates_min3, coordinates_max3)
+mesh3 = StructuredMesh(cells_per_dimension, coordinates_min3, coordinates_max3,
+                       periodicity = false)
 
 # Define the coupling functions
 coupling_function34 = (x, u, equations_other, equations_own) -> u
@@ -146,7 +149,8 @@ semi3 = SemidiscretizationHyperbolic(mesh3, equations, initial_condition_converg
 coordinates_min4 = (0.0, -1.0) # minimum coordinates (min(x), min(y))
 coordinates_max4 = (1.0, 0.0) # maximum coordinates (max(x), max(y))
 
-mesh4 = StructuredMesh(cells_per_dimension, coordinates_min4, coordinates_max4)
+mesh4 = StructuredMesh(cells_per_dimension, coordinates_min4, coordinates_max4,
+                       periodicity = false)
 
 # Define the coupling functions
 coupling_function43 = (x, u, equations_other, equations_own) -> u

examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl

@@ -69,7 +69,7 @@ mapping(xi, eta) = SVector(0.25 * 0.5 * (1.0 + xi), 0.5 * (1.0 + eta))
 
 num_elements_per_dimension = 32
 cells_per_dimension = (num_elements_per_dimension, num_elements_per_dimension * 4)
-mesh = StructuredMesh(cells_per_dimension, mapping)
+mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false)
 
 initial_condition = initial_condition_rayleigh_taylor_instability
 boundary_conditions = (x_neg = boundary_condition_slip_wall,

examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl

@@ -100,7 +100,8 @@ 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)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = (true, false))
 
 semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver,
                                     source_terms = warm_bubble_setup,

examples/tree_1d_dgsem/elixir_euler_source_terms.jl

@@ -44,9 +44,12 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 0.8)
 
+time_series = TimeSeriesCallback(semi, [0.0, 0.33, 1.0], interval = 10)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
                         save_solution,
+                        time_series,
                         stepsize_callback)
 
 ###############################################################################

examples/tree_3d_dgsem/elixir_euler_source_terms.jl

@@ -41,9 +41,14 @@ save_solution = SaveSolutionCallback(interval = 100,
 
 stepsize_callback = StepsizeCallback(cfl = 0.6)
 
+time_series = TimeSeriesCallback(semi,
+                                 [(0.0, 0.0, 0.0), (0.33, 0.33, 0.33), (1.0, 1.0, 1.0)],
+                                 interval = 10)
+
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
                         save_solution,
+                        time_series,
                         stepsize_callback)
 
 ###############################################################################

examples/unstructured_2d_dgsem/elixir_euler_time_series.jl

@@ -0,0 +1,115 @@
+# An elixir that has an alternative convergence test that uses
+# the `TimeSeriesCallback` on several gauge points. Many of the
+# gauge points are selected as "stress tests" for the element
+# identification, e.g., a gauge point that lies on an
+# element corner of a curvilinear mesh
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+# Modify the manufactured solution test to use `L = sqrt(2)`
+# in the initial condition and source terms
+function initial_condition_convergence_shifted(x, t,
+                                               equations::CompressibleEulerEquations2D)
+    c = 2
+    A = 0.1
+    L = sqrt(2)
+    f = 1 / L
+    ω = 2 * pi * f
+    ini = c + A * sin(ω * (x[1] + x[2] - t))
+
+    rho = ini
+    rho_v1 = ini
+    rho_v2 = ini
+    rho_e = ini^2
+
+    return SVector(rho, rho_v1, rho_v2, rho_e)
+end
+
+@inline function source_terms_convergence_shifted(u, x, t,
+                                                  equations::CompressibleEulerEquations2D)
+    # Same settings as in `initial_condition`
+    c = 2
+    A = 0.1
+    L = sqrt(2)
+    f = 1 / L
+    ω = 2 * pi * f
+    γ = equations.gamma
+
+    x1, x2 = x
+    si, co = sincos(ω * (x1 + x2 - t))
+    rho = c + A * si
+    rho_x = ω * A * co
+    # Note that d/dt rho = -d/dx rho = -d/dy rho.
+
+    tmp = (2 * rho - 1) * (γ - 1)
+
+    du1 = rho_x
+    du2 = rho_x * (1 + tmp)
+    du3 = du2
+    du4 = 2 * rho_x * (rho + tmp)
+
+    return SVector(du1, du2, du3, du4)
+end
+
+initial_condition = initial_condition_convergence_shifted
+
+source_term = source_terms_convergence_shifted
+
+###############################################################################
+# Get the DG approximation space
+
+solver = DGSEM(polydeg = 6, surface_flux = flux_lax_friedrichs)
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/b434e724e3972a9c4ee48d58c80cdcdb/raw/55c916cd8c0294a2d4a836e960dac7247b7c8ccf/mesh_multiple_flips.mesh",
+                           joinpath(@__DIR__, "mesh_multiple_flips.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = true)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_term)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+time_series = TimeSeriesCallback(semi,
+                                 [(0.75, 0.7), (1.23, 0.302), (0.8, 1.0),
+                                     (0.353553390593274, 0.353553390593274),
+                                     (0.505, 1.125), (1.37, 0.89), (0.349, 0.7153),
+                                     (0.883883476483184, 0.406586401289607),
+                                     (sqrt(2), sqrt(2))];
+                                 interval = 10)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        time_series,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-6, reltol = 1.0e-6,
+            ode_default_options()..., callback = callbacks);
+
+summary_callback() # print the timer summary

examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl

@@ -0,0 +1,86 @@
+# !!! warning "Experimental implementation (upwind SBP)"
+#     This is an experimental feature and may change in future releases.
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_constant
+
+# Boundary conditions for free-stream preservation test
+boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
+
+boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream,
+                           :cone1 => boundary_condition_free_stream,
+                           :cone2 => boundary_condition_free_stream,
+                           :iceCream => boundary_condition_free_stream)
+
+###############################################################################
+# Get the Upwind FDSBP approximation space
+
+# TODO: FDSBP
+# Note, one must set `xmin=-1` and `xmax=1` due to the reuse
+# of interpolation routines from `calc_node_coordinates!` to create
+# the physical coordinates in the mappings.
+D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017,
+                         derivative_order = 1,
+                         accuracy_order = 8,
+                         xmin = -1.0, xmax = 1.0,
+                         N = 17)
+
+flux_splitting = splitting_vanleer_haenel
+solver = FDSBP(D_upw,
+               surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)),
+               volume_integral = VolumeIntegralUpwind(flux_splitting))
+
+###############################################################################
+# Get the curved quad mesh from a file (downloads the file if not available locally)
+
+# Mesh with second-order boundary polynomials requires an upwind SBP operator
+# with (at least) 4th order boundary closure to guarantee the approximation is
+# free-stream preserving
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/ec9a345f09199ebe471d35d5c1e4e08f/raw/15975943d8642e42f8292235314b6f1b30aa860d/mesh_inner_outer_boundaries.mesh",
+                           joinpath(@__DIR__, "mesh_inner_outer_boundaries.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+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)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        save_solution,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+# set small tolerances for the free-stream preservation test
+sol = solve(ode, SSPRK43(), abstol = 1.0e-12, reltol = 1.0e-12,
+            save_everystep = false, callback = callbacks)
+
+summary_callback() # print the timer summary