Add functionality for TimeSeries callback on UnstructuredMesh2D (#…

…1855) * add functionality for TimeSeries callback on UnstructuredMesh2D * Update src/callbacks_step/time_series_dg.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Apply suggestions from code review Co-authored-by: Daniel Doehring <[email protected]> * add strategy to correctly locate a gauge point within a curvilinear element * add sanity check that the Newton solution is correct * run formatter * implement a more general approach that also works on curved element without issue * run formatter * forgot to format the examples * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> Co-authored-by: Daniel Doehring <[email protected]> * working version of the element finding routine * run formatter * add new elixir for the time series callback * add additional test for the time series callback on an unstructured mesh * add appropriate test * update docstring * add comment about the barycenter computation * add simplifications and comments from code review * adjust variable name to avoid ugly formatting * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * fix variable name * remove Experimental status from the TimeSeriesCallback * move new TimeSeries test into the unit testing * add output_directory creation if not already done. Necessary if this callback is used without the SaveSolution callback * formatting * update test mesh to have one straight-sided element to trigger inverse bi-linear interpolation * update test values * add news item * forgot to update all new test values on the new mesh * update tests and use coverage override to avoid redundancy --------- Co-authored-by: Hendrik Ranocha <[email protected]> Co-authored-by: Daniel Doehring <[email protected]>
trixi-framework · Mar 14, 2024 · bd060b6 · bd060b6
1 parent f235619
commit bd060b6
Show file tree

Hide file tree

Showing 7 changed files with 443 additions and 12 deletions.
diff --git a/NEWS.md b/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.

diff --git a/examples/unstructured_2d_dgsem/elixir_euler_time_series.jl b/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
diff --git a/src/callbacks_step/time_series.jl b/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)

diff --git a/src/callbacks_step/time_series_dg.jl b/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