Merge branch 'main' into subcell-limiting-outsource-saving-errors

NEWS.md

@@ -4,13 +4,19 @@ 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`.
+
+
 ## 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 +32,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 +46,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.2-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/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

src/callbacks_step/time_series.jl

@@ -23,8 +23,8 @@ After the last time step, the results are stored in an HDF5 file `filename` in d
 The real data type `RealT` and data type for solution variables `uEltype` default to the respective
 types used in the solver and the cache.
 
-!!! warning "Experimental implementation"
-    This is an experimental feature and may change in future releases.
+Currently this callback is only implemented for [`TreeMesh`](@ref) in 2D
+and [`UnstructuredMesh2D`](@ref).
 """
 mutable struct TimeSeriesCallback{RealT <: Real, uEltype <: Real, SolutionVariables,
                                   VariableNames, Cache}
@@ -96,6 +96,11 @@ function TimeSeriesCallback(mesh, equations, solver, cache, point_coordinates;
         throw(ArgumentError("`point_coordinates` must be a matrix of size n_points × ndims"))
     end
 
+    # create the output folder if it does not exist already
+    if mpi_isroot() && !isdir(output_directory)
+        mkpath(output_directory)
+    end
+
     # Transpose point_coordinates to our usual format [ndims, n_points]
     # Note: They are accepted in a different format to allow direct input from `readdlm`
     point_coordinates_ = permutedims(point_coordinates)

src/callbacks_step/time_series_dg.jl

@@ -5,8 +5,10 @@
 @muladd begin
 #! format: noindent
 
-# Store time series file for a TreeMesh with a DG solver
-function save_time_series_file(time_series_callback, mesh::TreeMesh, equations, dg::DG)
+# Store time series file for a DG solver
+function save_time_series_file(time_series_callback,
+                               mesh::Union{TreeMesh, UnstructuredMesh2D},
+                               equations, dg::DG)
     @unpack (interval, solution_variables, variable_names,
     output_directory, filename, point_coordinates,
     point_data, time, step, time_series_cache) = time_series_callback