From 9284b1111ef3a70b99d9026483c903895a76a740 Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Fri, 26 Apr 2024 22:28:53 +0530
Subject: [PATCH 01/59] Add analysis_surface_2d.jl

---
 Project.toml                                  |   1 +
 ...xir_navierstokes_NACA0012airfoil_mach08.jl |  19 +-
 src/Trixi.jl                                  |   6 +-
 src/callbacks_step/analysis.jl                |   9 +-
 src/callbacks_step/analysis_surface_2d.jl     | 278 ++++++++++++++++++
 5 files changed, 304 insertions(+), 9 deletions(-)
 create mode 100644 src/callbacks_step/analysis_surface_2d.jl

diff --git a/Project.toml b/Project.toml
index 483e9cd8233..8bc75b0cb69 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.7.9-pre"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 1b485913ab2..0659ebeeb5b 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -1,7 +1,6 @@
 using Downloads: download
 using OrdinaryDiffEq
 using Trixi
-
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 
@@ -15,7 +14,7 @@ using Trixi
 #   Structured and Unstructured Grid Methods (2016)
 #   [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623)
 # - Deep Ray, Praveen Chandrashekar (2017)
-#   An entropy stable finite volume scheme for the 
+#   An entropy stable finite volume scheme for the
 #   two dimensional Navier–Stokes equations on triangular grids
 #   [DOI:10.1016/j.amc.2017.07.020](https://doi.org/10.1016/j.amc.2017.07.020)
 
@@ -142,6 +141,18 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
+
+# TODO - Should these two be moved to a different callback to avoid printing "Dummy Value" so often?
+# This will be fixed if the "Dummy Callback" is somehow not printed at all.
+# Can this be done by overloading on `pretty_form_utf`?
+friction_coefficient = AnalysisSurface(semi, force_boundary_names,
+                                       SurfaceFrictionCoefficient(rho_inf(), u_inf(equations), l_inf()))
+
+pressure_coefficient = AnalysisSurface(semi, force_boundary_names,
+                                       SurfacePressureCoefficient(p_inf(), rho_inf(),
+                                       u_inf(equations),
+                                       l_inf()))
+
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",
                                      save_analysis = true,
@@ -149,7 +160,9 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      analysis_integrals = (drag_coefficient,
                                                            lift_coefficient,
                                                            drag_coefficient_shear_force,
-                                                           lift_coefficient_shear_force))
+                                                           lift_coefficient_shear_force,
+                                                           friction_coefficient,
+                                                           pressure_coefficient))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 5511f7e02e2..2c792d114bf 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -42,6 +42,7 @@ import SciMLBase: get_du, get_tmp_cache, u_modified!,
 using Downloads: Downloads
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
+using DelimitedFiles: writedlm
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 @reexport using EllipsisNotation # ..
 using FillArrays: Ones, Zeros
@@ -262,8 +263,9 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
-       DragCoefficientShearStress, LiftCoefficientShearStress
+       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient, SurfaceFrictionCoefficient
+
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 8f89af755a2..1a0e7de14da 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -26,8 +26,8 @@ or `extra_analysis_errors = (:conservation_error,)`.
 If you want to omit the computation (to safe compute-time) of the [`default_analysis_errors`](@ref), specify
 `analysis_errors = Symbol[]`.
 Note: `default_analysis_errors` are `:l2_error` and `:linf_error` for all equations.
-If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e., 
-without `:l2_error, :linf_error` you need to specify 
+If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e.,
+without `:l2_error, :linf_error` you need to specify
 `analysis_errors = [:conservation_error]` instead of `extra_analysis_errors = [:conservation_error]`.
 
 Further scalar functions `func` in `extra_analysis_integrals` are applied to the numerical
@@ -690,15 +690,16 @@ end # @muladd
 include("analysis_dg1d.jl")
 include("analysis_dg2d.jl")
 include("analysis_surface_integral_2d.jl")
+include("analysis_surface_2d.jl")
 include("analysis_dg2d_parallel.jl")
 include("analysis_dg3d.jl")
 include("analysis_dg3d_parallel.jl")
 
 # Special analyze for `SemidiscretizationHyperbolicParabolic` such that
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
-# Note that this needs to be included after `analysis_surface_integral_2d.jl` to 
+# Note that this needs to be included after `analysis_surface_integral_2d.jl` to
 # have `VariableViscous` available.
-function analyze(quantity::AnalysisSurfaceIntegral{Variable},
+function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable}, AnalysisSurface{Variable}},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
                                                                      Variable <:
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
new file mode 100644
index 00000000000..3388091ce84
--- /dev/null
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -0,0 +1,278 @@
+# 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
+
+# This file contains callbacks that are performed on the surface like computation of
+# surface forces
+
+"""
+    AnalysisSurface{Semidiscretization, Variable}(semi,
+                                                          boundary_symbol_or_boundary_symbols,
+                                                          variable)
+
+This struct is used to compute pointwise surface values of a quantity of interest `variable`
+alongside the boundary/boundaries associated with particular name(s) given in `boundary_symbol`
+or `boundary_symbols`.
+For instance, this can be used to compute the surface pressure coefficient [`SurfacePressureCoefficient`](@ref) or
+surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an airfoil with the boundary
+name `:Airfoil` in 2D.
+
+- `semi::Semidiscretization`: Passed in to retrieve boundary condition information
+- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries
+  where the quantity of interest is computed
+- `variable::Variable`: Quantity of interest, like lift or drag
+"""
+struct AnalysisSurface{Variable}
+    indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
+    variable::Variable # Quantity of interest, like lift or drag
+
+    function AnalysisSurface(semi, boundary_symbol, variable)
+        @unpack boundary_symbol_indices = semi.boundary_conditions
+        indices = boundary_symbol_indices[boundary_symbol]
+
+        return new{typeof(variable)}(indices, variable)
+    end
+
+    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable)
+        @unpack boundary_symbol_indices = semi.boundary_conditions
+        indices = Vector{Int}()
+        for name in boundary_symbols
+            append!(indices, boundary_symbol_indices[name])
+        end
+        sort!(indices)
+
+        return new{typeof(variable)}(indices, variable)
+    end
+end
+
+struct SurfacePressureCoefficient{RealT <: Real}
+    pinf::RealT # Free stream pressure
+    force_state::ForceState{RealT}
+end
+
+struct SurfaceFrictionCoefficient{RealT <: Real} <: VariableViscous
+    force_state::ForceState{RealT}
+end
+
+"""
+    SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
+
+Compute the surface pressure coefficient
+```math
+C_p \\coloneqq \\frac{ p - p_{p_\\infty}}
+                        {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+```
+based on the pressure distribution along a boundary.
+Supposed to be used in conjunction with [`AnalysisSurface`](@ref)
+which stores the boundary information and semidiscretization.
+
+- `pinf::Real`: Free-stream pressure
+- `rhoinf::Real`: Free-stream density
+- `uinf::Real`: Free-stream velocity
+- `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
+"""
+function SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
+    psi = (zero(pinf), zero(pinf)) # Dummy value/placeholder, will not be used
+    return SurfacePressureCoefficient(pinf, ForceState(psi, rhoinf, uinf, linf))
+end
+
+"""
+SurfaceFrictionCoefficient(rhoinf, uinf, linf)
+
+Compute the surface skin friction coefficient
+```math
+C_f \\coloneqq \\frac{\\boldsymbol (\\tau_w  \\boldsymbol n) \\cdot \\boldsymbol b^\\perp}
+                        {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+```
+based on the wall shear stress vector ``\\tau_w`` along a boundary.
+Supposed to be used in conjunction with [`AnalysisSurface`](@ref)
+which stores the boundary information and semidiscretization.
+
+- `rhoinf::Real`: Free-stream density
+- `uinf::Real`: Free-stream velocity
+- `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
+"""
+function SurfaceFrictionCoefficient(rhoinf, uinf, linf)
+    psi = (zero(rhoinf), zero(rhoinf)) # Dummy value/placeholder, will not be used
+    return SurfaceFrictionCoefficient(ForceState(psi, rhoinf, uinf, linf))
+end
+
+function (pressure_coefficient::SurfacePressureCoefficient)(u, equations)
+    p = pressure(u, equations)
+    @unpack pinf = pressure_coefficient
+    @unpack rhoinf, uinf, linf = pressure_coefficient.force_state
+    return (p - pinf) / (0.5 * rhoinf * uinf^2 * linf)
+end
+
+function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x, t,
+                                                        equations_parabolic,
+                                                        gradients_1, gradients_2)
+    visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
+                                               gradients_1, gradients_2)
+    @unpack psi, rhoinf, uinf, linf = surface_friction.force_state
+
+    # Normalize as `normal_direction` is not necessarily a unit vector
+    n = normal_direction / norm(normal_direction)
+    n_perp = (-n[2], n[1])
+    return (visc_stress_vector[1] * n_perp[1] + visc_stress_vector[2] * n_perp[2]) /
+           (0.5 * rhoinf * uinf^2 * linf)
+end
+
+function analyze(surface_variable::AnalysisSurface, du, u, t,
+                 mesh::P4estMesh{2},
+                 equations, dg::DGSEM, cache)
+    @unpack boundaries = cache
+    @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
+    @unpack weights = dg.basis
+
+    @unpack indices, variable = surface_variable
+
+    dim = ndims(mesh)
+    n_nodes = nnodes(dg)
+    n_elements = length(indices)
+
+    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
+    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+
+    # TODO - Can variable_array be a vector?
+    # TODO - Decide whether to save element mean values too
+
+    index_range = eachnode(dg)
+
+    global_node_index = 1 # Keeps track of solution point number on the surface
+    for boundary in indices
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+        direction = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in index_range
+            u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
+                                         boundary)
+
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            var = variable(u_node, equations)
+
+            coord_array[global_node_index, 1] = x[1]
+            coord_array[global_node_index, 2] = x[2]
+            variable_array[global_node_index, 1] = var
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_index += 1
+        end
+    end
+    # TODO - Sort coord_array, variable_array increasing x?
+    mkpath("out")
+    t_trunc = @sprintf("%.3f", t)
+    filename = varname(variable) * t_trunc * ".txt"
+    # TODO - Should we start with a bigger array and avoid hcat?
+    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
+    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+end
+
+function analyze(surface_variable::AnalysisSurface{Variable},
+                 du, u, t, mesh::P4estMesh{2},
+                 equations, equations_parabolic,
+                 dg::DGSEM, cache,
+                 cache_parabolic) where {Variable <: VariableViscous}
+    @unpack boundaries = cache
+    @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
+    @unpack weights = dg.basis
+
+    @unpack indices, variable = surface_variable
+
+    dim = ndims(mesh)
+    n_nodes = nnodes(dg)
+    n_elements = length(indices)
+
+    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
+    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+
+    # TODO - Should variable_array be a vector?
+    # TODO - Decide whether to save element mean values too
+
+    # Additions for parabolic
+    @unpack viscous_container = cache_parabolic
+    @unpack gradients = viscous_container
+
+    gradients_x, gradients_y = gradients
+
+    surface_integral = zero(eltype(u))
+    index_range = eachnode(dg)
+    global_node_index = 1 # Keeps track of solution point number on the surface
+    for boundary in indices
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+        direction = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in index_range
+            u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
+                                         boundary)
+
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            # Extract normal direction at nodes which points from the elements outwards,
+            # i.e., *into* the structure.
+            normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                                    i_node, j_node,
+                                                    element)
+
+            gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg, i_node,
+                                        j_node, element)
+            gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg, i_node,
+                                        j_node, element)
+
+            # Integral over whole boundary surface
+            var = variable(u_node, normal_direction, x, t, equations_parabolic,
+                           gradients_1, gradients_2)
+
+            coord_array[global_node_index, 1] = x[1]
+            coord_array[global_node_index, 2] = x[2]
+            variable_array[global_node_index, 1] = var
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_index += 1
+        end
+    end
+    # TODO - Sort coord_array, variable_array increasing x?
+    mkpath("out")
+    t_trunc = @sprintf("%.3f", t)
+    filename = varname(variable) * t_trunc * ".txt"
+    # TODO - Should we start with a bigger array and avoid hcat?
+    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
+    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+end
+
+varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
+varname(pressure_coefficient::SurfacePressureCoefficient) = "Cp"
+varname(friction_coefficient::SurfaceFrictionCoefficient) = "Cf"
+
+
+# AnalysisCallback prints to screen and needs these. TODO - Can it be turned off?
+function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
+    "Dummy value"
+end
+function pretty_form_utf(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
+    "Dummy value"
+end
+
+function pretty_form_ascii(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
+    "Dummy value"
+end
+function pretty_form_utf(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
+    "Dummy value"
+end
+
+
+end # muladd

From 0f6ad41d062dcf22ad0106a35f7fdc1d09ace95d Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Fri, 26 Apr 2024 22:46:06 +0530
Subject: [PATCH 02/59] Format

---
 .../elixir_navierstokes_NACA0012airfoil_mach08.jl        | 9 +++++----
 src/Trixi.jl                                             | 7 ++++---
 src/callbacks_step/analysis.jl                           | 3 ++-
 src/callbacks_step/analysis_surface_2d.jl                | 9 ++++-----
 4 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 0659ebeeb5b..94ee43b3cbd 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -141,17 +141,18 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-
 # TODO - Should these two be moved to a different callback to avoid printing "Dummy Value" so often?
 # This will be fixed if the "Dummy Callback" is somehow not printed at all.
 # Can this be done by overloading on `pretty_form_utf`?
 friction_coefficient = AnalysisSurface(semi, force_boundary_names,
-                                       SurfaceFrictionCoefficient(rho_inf(), u_inf(equations), l_inf()))
+                                       SurfaceFrictionCoefficient(rho_inf(),
+                                                                  u_inf(equations),
+                                                                  l_inf()))
 
 pressure_coefficient = AnalysisSurface(semi, force_boundary_names,
                                        SurfacePressureCoefficient(p_inf(), rho_inf(),
-                                       u_inf(equations),
-                                       l_inf()))
+                                                                  u_inf(equations),
+                                                                  l_inf()))
 
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 2c792d114bf..71873e233d5 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -263,9 +263,10 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
-       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient, SurfaceFrictionCoefficient
-
+       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure,
+       LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient,
+       SurfaceFrictionCoefficient
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 1a0e7de14da..a0fae2b6a85 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -699,7 +699,8 @@ include("analysis_dg3d_parallel.jl")
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
 # Note that this needs to be included after `analysis_surface_integral_2d.jl` to
 # have `VariableViscous` available.
-function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable}, AnalysisSurface{Variable}},
+function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
+                                 AnalysisSurface{Variable}},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
                                                                      Variable <:
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 3388091ce84..5a6b887a596 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -157,7 +157,8 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
+                                element)
             var = variable(u_node, equations)
 
             coord_array[global_node_index, 1] = x[1]
@@ -221,7 +222,8 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
+                                element)
             # Extract normal direction at nodes which points from the elements outwards,
             # i.e., *into* the structure.
             normal_direction = get_normal_direction(direction, contravariant_vectors,
@@ -258,7 +260,6 @@ varname(::Any) = @assert false "Surface variable name not assigned" # This makes
 varname(pressure_coefficient::SurfacePressureCoefficient) = "Cp"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "Cf"
 
-
 # AnalysisCallback prints to screen and needs these. TODO - Can it be turned off?
 function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
     "Dummy value"
@@ -273,6 +274,4 @@ end
 function pretty_form_utf(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
     "Dummy value"
 end
-
-
 end # muladd

From 2e79e4ec40b39743ceca77fc86f0371a139439ad Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Tue, 30 Apr 2024 09:50:27 +0200
Subject: [PATCH 03/59] Apply suggestions from code review

---
 src/callbacks_step/analysis_surface_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 5a6b887a596..6062affeaf4 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 # This file contains callbacks that are performed on the surface like computation of
-# surface forces
+# pointwise surface forces.
 
 """
     AnalysisSurface{Semidiscretization, Variable}(semi,

From 500518a4e01225694b3df168f88244e1849efb7f Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 30 Apr 2024 10:59:50 +0200
Subject: [PATCH 04/59] Pointwise analysis

---
 ...xir_navierstokes_NACA0012airfoil_mach08.jl |  4 +-
 src/callbacks_step/analysis.jl                | 39 ++++++++++-
 src/callbacks_step/analysis_surface_2d.jl     | 69 +++++++++----------
 .../analysis_surface_integral_2d.jl           | 28 ++++----
 4 files changed, 87 insertions(+), 53 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 94ee43b3cbd..477c58cd8a0 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -161,8 +161,8 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      analysis_integrals = (drag_coefficient,
                                                            lift_coefficient,
                                                            drag_coefficient_shear_force,
-                                                           lift_coefficient_shear_force,
-                                                           friction_coefficient,
+                                                           lift_coefficient_shear_force),
+                                     analysis_pointwise = (friction_coefficient,
                                                            pressure_coefficient))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index a0fae2b6a85..5a9cc964c9b 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -43,7 +43,8 @@ evaluating the computational performance, such as the total runtime, the perform
 (time/DOF/rhs!), the time spent in garbage collection (GC), or the current memory usage (alloc'd
 memory).
 """
-mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegrals,
+mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, AnalysisPointwise,
+                                InitialStateIntegrals,
                                 Cache}
     start_time::Float64
     start_time_last_analysis::Float64
@@ -56,6 +57,7 @@ mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegra
     analyzer::Analyzer
     analysis_errors::Vector{Symbol}
     analysis_integrals::AnalysisIntegrals
+    analysis_pointwise::AnalysisPointwise
     initial_state_integrals::InitialStateIntegrals
     cache::Cache
 end
@@ -80,6 +82,9 @@ function Base.show(io::IO, ::MIME"text/plain",
         for (idx, integral) in enumerate(analysis_callback.analysis_integrals)
             push!(setup, "│ integral " * string(idx) => integral)
         end
+        for (idx, quantity) in enumerate(analysis_callback.analysis_pointwise)
+            push!(setup, "│ quantity " * string(idx) => quantity)
+        end
         push!(setup,
               "save analysis to file" => analysis_callback.save_analysis ? "yes" : "no")
         if analysis_callback.save_analysis
@@ -109,6 +114,7 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                           extra_analysis_integrals = (),
                           analysis_integrals = union(default_analysis_integrals(equations),
                                                      extra_analysis_integrals),
+                          analysis_pointwise = (),
                           RealT = real(solver),
                           uEltype = eltype(cache.elements),
                           kwargs...)
@@ -132,6 +138,7 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                                          analysis_filename,
                                          analyzer,
                                          analysis_errors, Tuple(analysis_integrals),
+                                         Tuple(analysis_pointwise),
                                          SVector(ntuple(_ -> zero(uEltype),
                                                         Val(nvariables(equations)))),
                                          cache_analysis)
@@ -158,7 +165,7 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
 
     analysis_callback = cb.affect!
     analysis_callback.initial_state_integrals = initial_state_integrals
-    @unpack save_analysis, output_directory, analysis_filename, analysis_errors, analysis_integrals = analysis_callback
+    @unpack save_analysis, output_directory, analysis_filename, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
 
     if save_analysis && mpi_isroot()
         mkpath(output_directory)
@@ -203,6 +210,10 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
                 @printf(io, "   %-14s", pretty_form_ascii(quantity))
             end
 
+            for quantity in analysis_pointwise
+                @printf(io, "   %-14s", pretty_form_ascii(quantity))
+            end
+
             println(io)
         end
     end
@@ -368,7 +379,7 @@ end
 # This method is just called internally from `(analysis_callback::AnalysisCallback)(integrator)`
 # and serves as a function barrier. Additionally, it makes the code easier to profile and optimize.
 function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
-    @unpack analyzer, analysis_errors, analysis_integrals = analysis_callback
+    @unpack analyzer, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
     cache_analysis = analysis_callback.cache
     _, equations, _, _ = mesh_equations_solver_cache(semi)
 
@@ -486,6 +497,8 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
 
     # additional integrals
     analyze_integrals(analysis_integrals, io, du, u, t, semi)
+    # additional pointwise quantities
+    analyze_pointwise(analysis_pointwise, du, u, t, semi)
 
     return nothing
 end
@@ -606,6 +619,26 @@ function analyze_integrals(analysis_integrals::Tuple{}, io, du, u, t, semi)
     nothing
 end
 
+# Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
+function analyze_pointwise(analysis_quantities::NTuple{N, Any}, du, u, t,
+                           semi) where {N}
+
+    # Extract the first analysis integral and process it; keep the remaining to be processed later
+    quantity = first(analysis_quantities)
+    remaining_quantities = Base.tail(analysis_quantities)
+
+    analyze(quantity, du, u, t, semi)
+
+    # Recursively call this method with the unprocessed integrals
+    analyze_pointwise(remaining_quantities, du, u, t, semi)
+    return nothing
+end
+
+# terminate the type-stable iteration over tuples
+function analyze_pointwise(analysis_quantities::Tuple{}, du, u, t, semi)
+    nothing
+end
+
 # used for error checks and EOC analysis
 function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
                                                                 Affect! <:
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 6062affeaf4..24cf626f8e4 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -48,13 +48,19 @@ struct AnalysisSurface{Variable}
     end
 end
 
+struct FlowState{RealT <: Real}
+    rhoinf::RealT
+    uinf::RealT
+    linf::RealT
+end
+
 struct SurfacePressureCoefficient{RealT <: Real}
     pinf::RealT # Free stream pressure
-    force_state::ForceState{RealT}
+    flow_state::FlowState{RealT}
 end
 
 struct SurfaceFrictionCoefficient{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowState{RealT}
 end
 
 """
@@ -75,8 +81,7 @@ which stores the boundary information and semidiscretization.
 - `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
 """
 function SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
-    psi = (zero(pinf), zero(pinf)) # Dummy value/placeholder, will not be used
-    return SurfacePressureCoefficient(pinf, ForceState(psi, rhoinf, uinf, linf))
+    return SurfacePressureCoefficient(pinf, FlowState(rhoinf, uinf, linf))
 end
 
 """
@@ -96,14 +101,13 @@ which stores the boundary information and semidiscretization.
 - `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
 """
 function SurfaceFrictionCoefficient(rhoinf, uinf, linf)
-    psi = (zero(rhoinf), zero(rhoinf)) # Dummy value/placeholder, will not be used
-    return SurfaceFrictionCoefficient(ForceState(psi, rhoinf, uinf, linf))
+    return SurfaceFrictionCoefficient(FlowState(rhoinf, uinf, linf))
 end
 
 function (pressure_coefficient::SurfacePressureCoefficient)(u, equations)
     p = pressure(u, equations)
     @unpack pinf = pressure_coefficient
-    @unpack rhoinf, uinf, linf = pressure_coefficient.force_state
+    @unpack rhoinf, uinf, linf = pressure_coefficient.flow_state
     return (p - pinf) / (0.5 * rhoinf * uinf^2 * linf)
 end
 
@@ -112,7 +116,7 @@ function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x,
                                                         gradients_1, gradients_2)
     visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = surface_friction.force_state
+    @unpack rhoinf, uinf, linf = surface_friction.flow_state
 
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = normal_direction / norm(normal_direction)
@@ -134,10 +138,9 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     n_nodes = nnodes(dg)
     n_elements = length(indices)
 
-    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
-    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+    coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Can variable_array be a vector?
     # TODO - Decide whether to save element mean values too
 
     index_range = eachnode(dg)
@@ -161,21 +164,20 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
                                 element)
             var = variable(u_node, equations)
 
-            coord_array[global_node_index, 1] = x[1]
-            coord_array[global_node_index, 2] = x[2]
-            variable_array[global_node_index, 1] = var
+            coords[global_node_index, 1] = x[1]
+            coords[global_node_index, 2] = x[2]
+            variables[global_node_index] = var
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coord_array, variable_array increasing x?
+    # TODO - Sort coords, variables increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * t_trunc * ".txt"
+    filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
-    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+    writedlm(joinpath("out", filename), hcat(coords, variables))
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
@@ -193,10 +195,9 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     n_nodes = nnodes(dg)
     n_elements = length(indices)
 
-    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
-    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+    coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Should variable_array be a vector?
     # TODO - Decide whether to save element mean values too
 
     # Additions for parabolic
@@ -239,39 +240,37 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             var = variable(u_node, normal_direction, x, t, equations_parabolic,
                            gradients_1, gradients_2)
 
-            coord_array[global_node_index, 1] = x[1]
-            coord_array[global_node_index, 2] = x[2]
-            variable_array[global_node_index, 1] = var
+            coords[global_node_index, 1] = x[1]
+            coords[global_node_index, 2] = x[2]
+            variables[global_node_index] = var
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coord_array, variable_array increasing x?
+    # TODO - Sort coords, variables increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * t_trunc * ".txt"
+    filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
-    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+    writedlm(joinpath("out", filename), hcat(coords, variables))
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
-varname(pressure_coefficient::SurfacePressureCoefficient) = "Cp"
-varname(friction_coefficient::SurfaceFrictionCoefficient) = "Cf"
+varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
+varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
-# AnalysisCallback prints to screen and needs these. TODO - Can it be turned off?
 function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
-    "Dummy value"
+    "CP(x)"
 end
 function pretty_form_utf(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
-    "Dummy value"
+    "CP(x)"
 end
 
 function pretty_form_ascii(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
-    "Dummy value"
+    "CF(x)"
 end
 function pretty_form_utf(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
-    "Dummy value"
+    "CF(x)"
 end
 end # muladd
diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index 7ae259e5285..e58769f3616 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -48,7 +48,7 @@ struct AnalysisSurfaceIntegral{Variable}
     end
 end
 
-struct ForceState{RealT <: Real}
+struct FlowStateDirectional{RealT <: Real}
     psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream
     rhoinf::RealT
     uinf::RealT
@@ -56,11 +56,11 @@ struct ForceState{RealT <: Real}
 end
 
 struct LiftCoefficientPressure{RealT <: Real}
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 struct DragCoefficientPressure{RealT <: Real}
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 # Abstract base type used for dispatch of `analyze` for quantities
@@ -68,11 +68,11 @@ end
 abstract type VariableViscous end
 
 struct LiftCoefficientShearStress{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 struct DragCoefficientShearStress{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 """
@@ -99,7 +99,7 @@ function LiftCoefficientPressure(aoa, rhoinf, uinf, linf)
     # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same
     # value, but with the opposite sign.
     psi_lift = (-sin(aoa), cos(aoa))
-    return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, linf))
+    return LiftCoefficientPressure(FlowStateDirectional(psi_lift, rhoinf, uinf, linf))
 end
 
 """
@@ -122,7 +122,7 @@ which stores the boundary information and semidiscretization.
 function DragCoefficientPressure(aoa, rhoinf, uinf, linf)
     # `psi_drag` is the unit vector tangent to the freestream direction
     psi_drag = (cos(aoa), sin(aoa))
-    return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf))
+    return DragCoefficientPressure(FlowStateDirectional(psi_drag, rhoinf, uinf, linf))
 end
 
 """
@@ -149,7 +149,8 @@ function LiftCoefficientShearStress(aoa, rhoinf, uinf, linf)
     # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same
     # value, but with the opposite sign.
     psi_lift = (-sin(aoa), cos(aoa))
-    return LiftCoefficientShearStress(ForceState(psi_lift, rhoinf, uinf, linf))
+    return LiftCoefficientShearStress(FlowStateDirectional(psi_lift, rhoinf, uinf,
+                                                           linf))
 end
 
 """
@@ -172,13 +173,14 @@ which stores the boundary information and semidiscretization.
 function DragCoefficientShearStress(aoa, rhoinf, uinf, linf)
     # `psi_drag` is the unit vector tangent to the freestream direction
     psi_drag = (cos(aoa), sin(aoa))
-    return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf))
+    return DragCoefficientShearStress(FlowStateDirectional(psi_drag, rhoinf, uinf,
+                                                           linf))
 end
 
 function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, x, t,
                                                      equations)
     p = pressure(u, equations)
-    @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = dot(normal_direction, psi) / norm(normal_direction)
     return p * n / (0.5 * rhoinf * uinf^2 * linf)
@@ -187,7 +189,7 @@ end
 function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, x, t,
                                                      equations)
     p = pressure(u, equations)
-    @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = dot(normal_direction, psi) / norm(normal_direction)
     return p * n / (0.5 * rhoinf * uinf^2 * linf)
@@ -240,7 +242,7 @@ function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, x,
                                                         gradients_1, gradients_2)
     visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
     return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
@@ -250,7 +252,7 @@ function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x,
                                                         gradients_1, gradients_2)
     visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
     return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end

From 9c1cc7a9649bef6fd8db5e19f2465a89b6972e77 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 3 May 2024 11:17:43 +0200
Subject: [PATCH 05/59] cosmetics

---
 src/callbacks_step/analysis_surface_2d.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 24cf626f8e4..1605149ad94 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -10,8 +10,8 @@
 
 """
     AnalysisSurface{Semidiscretization, Variable}(semi,
-                                                          boundary_symbol_or_boundary_symbols,
-                                                          variable)
+                                                  boundary_symbol_or_boundary_symbols,
+                                                  variable)
 
 This struct is used to compute pointwise surface values of a quantity of interest `variable`
 alongside the boundary/boundaries associated with particular name(s) given in `boundary_symbol`

From 9e4ad1239ed77738fddc84f6430b6f1d7cc753ea Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 3 May 2024 13:52:32 +0200
Subject: [PATCH 06/59] routines similar to timeseries

---
 src/callbacks_step/analysis_surface_2d.jl | 51 ++++++++++++++++-------
 1 file changed, 36 insertions(+), 15 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 1605149ad94..5b407602e49 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -28,15 +28,19 @@ name `:Airfoil` in 2D.
 struct AnalysisSurface{Variable}
     indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
     variable::Variable # Quantity of interest, like lift or drag
+    interval::Int
+    output_directory::String
 
-    function AnalysisSurface(semi, boundary_symbol, variable)
+    function AnalysisSurface(semi, boundary_symbol, variable, interval, 
+                             output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = boundary_symbol_indices[boundary_symbol]
 
-        return new{typeof(variable)}(indices, variable)
+        return new{typeof(variable)}(indices, variable, interval, output_directory)
     end
 
-    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable)
+    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable, interval,
+                             output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = Vector{Int}()
         for name in boundary_symbols
@@ -44,7 +48,7 @@ struct AnalysisSurface{Variable}
         end
         sort!(indices)
 
-        return new{typeof(variable)}(indices, variable)
+        return new{typeof(variable)}(indices, variable, interval, output_directory)
     end
 end
 
@@ -139,7 +143,7 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     n_elements = length(indices)
 
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
-    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
+    values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
     # TODO - Decide whether to save element mean values too
 
@@ -162,22 +166,22 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
 
             x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
                                 element)
-            var = variable(u_node, equations)
+            value = variable(u_node, equations)
 
             coords[global_node_index, 1] = x[1]
             coords[global_node_index, 2] = x[2]
-            variables[global_node_index] = var
+            values[global_node_index] = value
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, variables increasing x?
+    # TODO - Sort coords, values increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
     filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, variables))
+    writedlm(joinpath("out", filename), hcat(coords, values))
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
@@ -196,7 +200,7 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     n_elements = length(indices)
 
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
-    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
+    values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
     # TODO - Decide whether to save element mean values too
 
@@ -237,29 +241,46 @@ function analyze(surface_variable::AnalysisSurface{Variable},
                                         j_node, element)
 
             # Integral over whole boundary surface
-            var = variable(u_node, normal_direction, x, t, equations_parabolic,
-                           gradients_1, gradients_2)
+            value = variable(u_node, normal_direction, x, t, equations_parabolic,
+                             gradients_1, gradients_2)
 
             coords[global_node_index, 1] = x[1]
             coords[global_node_index, 2] = x[2]
-            variables[global_node_index] = var
+            values[global_node_index] = value
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, variables increasing x?
+    # TODO - Sort coords, values increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
     filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, variables))
+    writedlm(joinpath("out", filename), hcat(coords, values))
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
 varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
+function save_pointwise_file(output_directory, varname, coords, values, t)
+    n_points = length(values)
+
+    h5open(joinpath(output_directory, varname) * ".h5", "w") do file
+        # Add context information as attributes
+        attributes(file)["n_points"] = n_points
+        attributes(file)["interval"] = interval
+        attributes(file)["variable_name"] = collect(varname)
+
+        file["time"] = t
+        file["point_coordinates"] = coords
+        for p in 1:n_points
+            file["point_data_$p"] = values[p]
+        end
+    end
+end
+
 function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
     "CP(x)"
 end

From 0fef7b48f6ab2033dbfbfb8c423c14683387618c Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:50:40 +0200
Subject: [PATCH 07/59] write pointwise data to jhdf5

---
 src/callbacks_step/analysis.jl            | 32 ++++++++++++++-----
 src/callbacks_step/analysis_surface_2d.jl | 39 +++++++++--------------
 2 files changed, 39 insertions(+), 32 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 5a9cc964c9b..60f828cf482 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -350,8 +350,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
         @notimeit timer() integrator.f(du_ode, u_ode, semi, t)
         u = wrap_array(u_ode, mesh, equations, solver, cache)
         du = wrap_array(du_ode, mesh, equations, solver, cache)
-        # Compute l2_error, linf_error
-        analysis_callback(io, du, u, u_ode, t, semi)
+        # Compute l2_error, linf_error among other quantities
+        analysis_callback(io, du, u, u_ode, t, semi, iter)
 
         mpi_println("─"^100)
         mpi_println()
@@ -378,7 +378,7 @@ end
 
 # This method is just called internally from `(analysis_callback::AnalysisCallback)(integrator)`
 # and serves as a function barrier. Additionally, it makes the code easier to profile and optimize.
-function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
+function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi, iter)
     @unpack analyzer, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
     cache_analysis = analysis_callback.cache
     _, equations, _, _ = mesh_equations_solver_cache(semi)
@@ -498,7 +498,7 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
     # additional integrals
     analyze_integrals(analysis_integrals, io, du, u, t, semi)
     # additional pointwise quantities
-    analyze_pointwise(analysis_pointwise, du, u, t, semi)
+    analyze_pointwise(analysis_pointwise, du, u, t, semi, iter)
 
     return nothing
 end
@@ -621,21 +621,21 @@ end
 
 # Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
 function analyze_pointwise(analysis_quantities::NTuple{N, Any}, du, u, t,
-                           semi) where {N}
+                           semi, iter) where {N}
 
     # Extract the first analysis integral and process it; keep the remaining to be processed later
     quantity = first(analysis_quantities)
     remaining_quantities = Base.tail(analysis_quantities)
 
-    analyze(quantity, du, u, t, semi)
+    analyze(quantity, du, u, t, semi, iter)
 
     # Recursively call this method with the unprocessed integrals
-    analyze_pointwise(remaining_quantities, du, u, t, semi)
+    analyze_pointwise(remaining_quantities, du, u, t, semi, iter)
     return nothing
 end
 
 # terminate the type-stable iteration over tuples
-function analyze_pointwise(analysis_quantities::Tuple{}, du, u, t, semi)
+function analyze_pointwise(analysis_quantities::Tuple{}, du, u, t, semi, iter)
     nothing
 end
 
@@ -660,6 +660,10 @@ function analyze(quantity, du, u, t, semi::AbstractSemidiscretization)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache)
 end
+function analyze(quantity, du, u, t, semi::AbstractSemidiscretization, iter)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    analyze(quantity, du, u, t, mesh, equations, solver, cache, iter)
+end
 function analyze(quantity, du, u, t, mesh, equations, solver, cache)
     integrate(quantity, u, mesh, equations, solver, cache, normalize = true)
 end
@@ -744,3 +748,15 @@ function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
     analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
             cache_parabolic)
 end
+function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
+                                 AnalysisSurface{Variable}},
+                 du, u, t,
+                 semi::SemidiscretizationHyperbolicParabolic, iter) where {
+                                                                     Variable <:
+                                                                     VariableViscous}
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    equations_parabolic = semi.equations_parabolic
+    cache_parabolic = semi.cache_parabolic
+    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
+            cache_parabolic, iter)
+end
\ No newline at end of file
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 5b407602e49..b5c6e8266f9 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -28,18 +28,17 @@ name `:Airfoil` in 2D.
 struct AnalysisSurface{Variable}
     indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
     variable::Variable # Quantity of interest, like lift or drag
-    interval::Int
     output_directory::String
 
-    function AnalysisSurface(semi, boundary_symbol, variable, interval, 
+    function AnalysisSurface(semi, boundary_symbol, variable, 
                              output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = boundary_symbol_indices[boundary_symbol]
 
-        return new{typeof(variable)}(indices, variable, interval, output_directory)
+        return new{typeof(variable)}(indices, variable, output_directory)
     end
 
-    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable, interval,
+    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable,
                              output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = Vector{Int}()
@@ -48,7 +47,7 @@ struct AnalysisSurface{Variable}
         end
         sort!(indices)
 
-        return new{typeof(variable)}(indices, variable, interval, output_directory)
+        return new{typeof(variable)}(indices, variable, output_directory)
     end
 end
 
@@ -131,7 +130,7 @@ end
 
 function analyze(surface_variable::AnalysisSurface, du, u, t,
                  mesh::P4estMesh{2},
-                 equations, dg::DGSEM, cache)
+                 equations, dg::DGSEM, cache, iter)
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
@@ -177,18 +176,14 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
         end
     end
     # TODO - Sort coords, values increasing x?
-    mkpath("out")
-    t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * "_" * t_trunc * ".txt"
-    # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, values))
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
                  du, u, t, mesh::P4estMesh{2},
                  equations, equations_parabolic,
                  dg::DGSEM, cache,
-                 cache_parabolic) where {Variable <: VariableViscous}
+                 cache_parabolic, iter) where {Variable <: VariableViscous}
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
@@ -253,31 +248,27 @@ function analyze(surface_variable::AnalysisSurface{Variable},
         end
     end
     # TODO - Sort coords, values increasing x?
-    mkpath("out")
-    t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * "_" * t_trunc * ".txt"
-    # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, values))
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
 varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
-function save_pointwise_file(output_directory, varname, coords, values, t)
+function save_pointwise_file(output_directory, varname, coords, values, t, iter)
     n_points = length(values)
 
-    h5open(joinpath(output_directory, varname) * ".h5", "w") do file
+    filename = joinpath(output_directory, varname) * @sprintf("_%06d.h5", iter)
+
+    h5open(filename, "w") do file
         # Add context information as attributes
         attributes(file)["n_points"] = n_points
-        attributes(file)["interval"] = interval
-        attributes(file)["variable_name"] = collect(varname)
+        attributes(file)["variable_name"] = varname
 
         file["time"] = t
+        file["timestep"] = iter
         file["point_coordinates"] = coords
-        for p in 1:n_points
-            file["point_data_$p"] = values[p]
-        end
+        file["values"] = values
     end
 end
 

From dd426b95e2c252a6de11fb4d926b8bb7e7ed0b71 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:51:24 +0200
Subject: [PATCH 08/59] fmt

---
 src/callbacks_step/analysis.jl            | 9 +++++----
 src/callbacks_step/analysis_surface_2d.jl | 8 +++++---
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 60f828cf482..fa1fc0782bd 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -751,12 +751,13 @@ end
 function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
                                  AnalysisSurface{Variable}},
                  du, u, t,
-                 semi::SemidiscretizationHyperbolicParabolic, iter) where {
-                                                                     Variable <:
-                                                                     VariableViscous}
+                 semi::SemidiscretizationHyperbolicParabolic,
+                 iter) where {
+                              Variable <:
+                              VariableViscous}
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     equations_parabolic = semi.equations_parabolic
     cache_parabolic = semi.cache_parabolic
     analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
             cache_parabolic, iter)
-end
\ No newline at end of file
+end
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index b5c6e8266f9..62f9a97c002 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -30,7 +30,7 @@ struct AnalysisSurface{Variable}
     variable::Variable # Quantity of interest, like lift or drag
     output_directory::String
 
-    function AnalysisSurface(semi, boundary_symbol, variable, 
+    function AnalysisSurface(semi, boundary_symbol, variable,
                              output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = boundary_symbol_indices[boundary_symbol]
@@ -176,7 +176,8 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
         end
     end
     # TODO - Sort coords, values increasing x?
-    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
+                        values, t, iter)
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
@@ -248,7 +249,8 @@ function analyze(surface_variable::AnalysisSurface{Variable},
         end
     end
     # TODO - Sort coords, values increasing x?
-    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
+                        values, t, iter)
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting

From f010eabbddc49b627c2efb300dd3702396066180 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:52:21 +0200
Subject: [PATCH 09/59] not needed

---
 Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9db9f6eb00f..dc5fe949f2d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ version = "0.7.10-pre"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

From c10c3c6020b1dd7ce3ff04ec46ee12d6156820d1 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:52:54 +0200
Subject: [PATCH 10/59] rm unnecessary

---
 src/Trixi.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index 71873e233d5..2dce892fea4 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -42,7 +42,6 @@ import SciMLBase: get_du, get_tmp_cache, u_modified!,
 using Downloads: Downloads
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
-using DelimitedFiles: writedlm
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 @reexport using EllipsisNotation # ..
 using FillArrays: Ones, Zeros

From a281eee09d2bedf572f4c8d4fa2e1a53d1e2fe11 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:54:01 +0200
Subject: [PATCH 11/59] re-organize

---
 src/Trixi.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index 2dce892fea4..b4b85d1f754 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -262,10 +262,10 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure,
-       LiftCoefficientPressure,
-       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient,
-       SurfaceFrictionCoefficient
+       AnalysisSurface, AnalysisSurfaceIntegral,
+       DragCoefficientPressure, LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress,
+       SurfacePressureCoefficient, SurfaceFrictionCoefficient
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!

From 178af8c3450a5cf1a84f1514e79643f6dac48d10 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:58:18 +0200
Subject: [PATCH 12/59] remove todo

---
 .../elixir_navierstokes_NACA0012airfoil_mach08.jl              | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 477c58cd8a0..ddbd3dc3652 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -141,9 +141,6 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-# TODO - Should these two be moved to a different callback to avoid printing "Dummy Value" so often?
-# This will be fixed if the "Dummy Callback" is somehow not printed at all.
-# Can this be done by overloading on `pretty_form_utf`?
 friction_coefficient = AnalysisSurface(semi, force_boundary_names,
                                        SurfaceFrictionCoefficient(rho_inf(),
                                                                   u_inf(equations),

From a0abcda9064939ca44c4547b27169bcf0c717fd6 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Tue, 7 May 2024 08:55:43 +0200
Subject: [PATCH 13/59] Apply suggestions from code review

Co-authored-by: Arpit Babbar <arpitbabbar@gmail.com>
---
 src/callbacks_step/analysis.jl            | 6 +++---
 src/callbacks_step/analysis_surface_2d.jl | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index fa1fc0782bd..fcf677dda3a 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -619,17 +619,17 @@ function analyze_integrals(analysis_integrals::Tuple{}, io, du, u, t, semi)
     nothing
 end
 
-# Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
+# Iterate over tuples of pointwise analysis quantities in a type-stable way using "lispy tuple programming".
 function analyze_pointwise(analysis_quantities::NTuple{N, Any}, du, u, t,
                            semi, iter) where {N}
 
-    # Extract the first analysis integral and process it; keep the remaining to be processed later
+    # Extract the first pointwise analysis quantity and process it; keep the remaining to be processed later
     quantity = first(analysis_quantities)
     remaining_quantities = Base.tail(analysis_quantities)
 
     analyze(quantity, du, u, t, semi, iter)
 
-    # Recursively call this method with the unprocessed integrals
+    # Recursively call this method with the unprocessed pointwise analysis quantities
     analyze_pointwise(remaining_quantities, du, u, t, semi, iter)
     return nothing
 end
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 62f9a97c002..fac2c28c734 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -21,7 +21,7 @@ surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an air
 name `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
-- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries
+- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries stored as symbol(s)
   where the quantity of interest is computed
 - `variable::Variable`: Quantity of interest, like lift or drag
 """

From 1310adeb15587f20752cabda0da689d69cc92f58 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 8 May 2024 15:33:54 +0200
Subject: [PATCH 14/59] Sort

---
 src/callbacks_step/analysis_surface_2d.jl | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 62f9a97c002..4f5fbc87b89 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -144,8 +144,6 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
     values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Decide whether to save element mean values too
-
     index_range = eachnode(dg)
 
     global_node_index = 1 # Keeps track of solution point number on the surface
@@ -175,7 +173,11 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, values increasing x?
+    # Sort with ascending x
+    sorting_indices = sortperm(coords[:, 1])
+    coords = coords[sorting_indices, :]
+    values = values[sorting_indices]
+
     save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
                         values, t, iter)
 end
@@ -198,8 +200,6 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
     values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Decide whether to save element mean values too
-
     # Additions for parabolic
     @unpack viscous_container = cache_parabolic
     @unpack gradients = viscous_container
@@ -248,7 +248,14 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, values increasing x?
+    # Sort with ascending x
+    sorting_indices = sortperm(coords[:, 1])
+    coords = coords[sorting_indices, :]
+    values = values[sorting_indices]
+
+    println(size(coords))
+    println(size(values))
+
     save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
                         values, t, iter)
 end

From de991636a76eb3c9ccd80c738fed2e60982f0fc0 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 8 May 2024 15:34:30 +0200
Subject: [PATCH 15/59] remove prnt

---
 src/callbacks_step/analysis_surface_2d.jl | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 4f5fbc87b89..14211c9521b 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -253,9 +253,6 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     coords = coords[sorting_indices, :]
     values = values[sorting_indices]
 
-    println(size(coords))
-    println(size(values))
-
     save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
                         values, t, iter)
 end

From bda1a2b5f4f7be449564a90431f0c2cb4aafb678 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 8 May 2024 15:36:24 +0200
Subject: [PATCH 16/59] name/symbol

---
 src/callbacks_step/analysis_surface_2d.jl          | 2 +-
 src/callbacks_step/analysis_surface_integral_2d.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 35749d010f4..ee09866b323 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -18,7 +18,7 @@ alongside the boundary/boundaries associated with particular name(s) given in `b
 or `boundary_symbols`.
 For instance, this can be used to compute the surface pressure coefficient [`SurfacePressureCoefficient`](@ref) or
 surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an airfoil with the boundary
-name `:Airfoil` in 2D.
+symbol `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
 - `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries stored as symbol(s)
diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index e58769f3616..db1ad0182c1 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -18,7 +18,7 @@ the boundary/boundaries associated with particular name(s) given in `boundary_sy
 or `boundary_symbols`.
 For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or
 drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary
-name `:Airfoil` in 2D.
+symbol `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
 - `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries

From 6f8eb6374cc2c95324230160a115b23d685cb7b7 Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Fri, 26 Apr 2024 22:28:53 +0530
Subject: [PATCH 17/59] Add analysis_surface_2d.jl

---
 Project.toml                                  |   1 +
 ...xir_navierstokes_NACA0012airfoil_mach08.jl |  19 +-
 src/Trixi.jl                                  |   6 +-
 src/callbacks_step/analysis.jl                |   9 +-
 src/callbacks_step/analysis_surface_2d.jl     | 278 ++++++++++++++++++
 5 files changed, 304 insertions(+), 9 deletions(-)
 create mode 100644 src/callbacks_step/analysis_surface_2d.jl

diff --git a/Project.toml b/Project.toml
index 769c1fe7d08..def3691a9e8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.7.13-pre"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 1b485913ab2..0659ebeeb5b 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -1,7 +1,6 @@
 using Downloads: download
 using OrdinaryDiffEq
 using Trixi
-
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 
@@ -15,7 +14,7 @@ using Trixi
 #   Structured and Unstructured Grid Methods (2016)
 #   [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623)
 # - Deep Ray, Praveen Chandrashekar (2017)
-#   An entropy stable finite volume scheme for the 
+#   An entropy stable finite volume scheme for the
 #   two dimensional Navier–Stokes equations on triangular grids
 #   [DOI:10.1016/j.amc.2017.07.020](https://doi.org/10.1016/j.amc.2017.07.020)
 
@@ -142,6 +141,18 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
+
+# TODO - Should these two be moved to a different callback to avoid printing "Dummy Value" so often?
+# This will be fixed if the "Dummy Callback" is somehow not printed at all.
+# Can this be done by overloading on `pretty_form_utf`?
+friction_coefficient = AnalysisSurface(semi, force_boundary_names,
+                                       SurfaceFrictionCoefficient(rho_inf(), u_inf(equations), l_inf()))
+
+pressure_coefficient = AnalysisSurface(semi, force_boundary_names,
+                                       SurfacePressureCoefficient(p_inf(), rho_inf(),
+                                       u_inf(equations),
+                                       l_inf()))
+
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",
                                      save_analysis = true,
@@ -149,7 +160,9 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      analysis_integrals = (drag_coefficient,
                                                            lift_coefficient,
                                                            drag_coefficient_shear_force,
-                                                           lift_coefficient_shear_force))
+                                                           lift_coefficient_shear_force,
+                                                           friction_coefficient,
+                                                           pressure_coefficient))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
diff --git a/src/Trixi.jl b/src/Trixi.jl
index f3977f1f058..80b9c38b850 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -42,6 +42,7 @@ import SciMLBase: get_du, get_tmp_cache, u_modified!,
 using Downloads: Downloads
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
+using DelimitedFiles: writedlm
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 @reexport using EllipsisNotation # ..
 using FillArrays: Ones, Zeros
@@ -263,8 +264,9 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
-       DragCoefficientShearStress, LiftCoefficientShearStress
+       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient, SurfaceFrictionCoefficient
+
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 8f89af755a2..1a0e7de14da 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -26,8 +26,8 @@ or `extra_analysis_errors = (:conservation_error,)`.
 If you want to omit the computation (to safe compute-time) of the [`default_analysis_errors`](@ref), specify
 `analysis_errors = Symbol[]`.
 Note: `default_analysis_errors` are `:l2_error` and `:linf_error` for all equations.
-If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e., 
-without `:l2_error, :linf_error` you need to specify 
+If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e.,
+without `:l2_error, :linf_error` you need to specify
 `analysis_errors = [:conservation_error]` instead of `extra_analysis_errors = [:conservation_error]`.
 
 Further scalar functions `func` in `extra_analysis_integrals` are applied to the numerical
@@ -690,15 +690,16 @@ end # @muladd
 include("analysis_dg1d.jl")
 include("analysis_dg2d.jl")
 include("analysis_surface_integral_2d.jl")
+include("analysis_surface_2d.jl")
 include("analysis_dg2d_parallel.jl")
 include("analysis_dg3d.jl")
 include("analysis_dg3d_parallel.jl")
 
 # Special analyze for `SemidiscretizationHyperbolicParabolic` such that
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
-# Note that this needs to be included after `analysis_surface_integral_2d.jl` to 
+# Note that this needs to be included after `analysis_surface_integral_2d.jl` to
 # have `VariableViscous` available.
-function analyze(quantity::AnalysisSurfaceIntegral{Variable},
+function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable}, AnalysisSurface{Variable}},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
                                                                      Variable <:
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
new file mode 100644
index 00000000000..3388091ce84
--- /dev/null
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -0,0 +1,278 @@
+# 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
+
+# This file contains callbacks that are performed on the surface like computation of
+# surface forces
+
+"""
+    AnalysisSurface{Semidiscretization, Variable}(semi,
+                                                          boundary_symbol_or_boundary_symbols,
+                                                          variable)
+
+This struct is used to compute pointwise surface values of a quantity of interest `variable`
+alongside the boundary/boundaries associated with particular name(s) given in `boundary_symbol`
+or `boundary_symbols`.
+For instance, this can be used to compute the surface pressure coefficient [`SurfacePressureCoefficient`](@ref) or
+surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an airfoil with the boundary
+name `:Airfoil` in 2D.
+
+- `semi::Semidiscretization`: Passed in to retrieve boundary condition information
+- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries
+  where the quantity of interest is computed
+- `variable::Variable`: Quantity of interest, like lift or drag
+"""
+struct AnalysisSurface{Variable}
+    indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
+    variable::Variable # Quantity of interest, like lift or drag
+
+    function AnalysisSurface(semi, boundary_symbol, variable)
+        @unpack boundary_symbol_indices = semi.boundary_conditions
+        indices = boundary_symbol_indices[boundary_symbol]
+
+        return new{typeof(variable)}(indices, variable)
+    end
+
+    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable)
+        @unpack boundary_symbol_indices = semi.boundary_conditions
+        indices = Vector{Int}()
+        for name in boundary_symbols
+            append!(indices, boundary_symbol_indices[name])
+        end
+        sort!(indices)
+
+        return new{typeof(variable)}(indices, variable)
+    end
+end
+
+struct SurfacePressureCoefficient{RealT <: Real}
+    pinf::RealT # Free stream pressure
+    force_state::ForceState{RealT}
+end
+
+struct SurfaceFrictionCoefficient{RealT <: Real} <: VariableViscous
+    force_state::ForceState{RealT}
+end
+
+"""
+    SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
+
+Compute the surface pressure coefficient
+```math
+C_p \\coloneqq \\frac{ p - p_{p_\\infty}}
+                        {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+```
+based on the pressure distribution along a boundary.
+Supposed to be used in conjunction with [`AnalysisSurface`](@ref)
+which stores the boundary information and semidiscretization.
+
+- `pinf::Real`: Free-stream pressure
+- `rhoinf::Real`: Free-stream density
+- `uinf::Real`: Free-stream velocity
+- `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
+"""
+function SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
+    psi = (zero(pinf), zero(pinf)) # Dummy value/placeholder, will not be used
+    return SurfacePressureCoefficient(pinf, ForceState(psi, rhoinf, uinf, linf))
+end
+
+"""
+SurfaceFrictionCoefficient(rhoinf, uinf, linf)
+
+Compute the surface skin friction coefficient
+```math
+C_f \\coloneqq \\frac{\\boldsymbol (\\tau_w  \\boldsymbol n) \\cdot \\boldsymbol b^\\perp}
+                        {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+```
+based on the wall shear stress vector ``\\tau_w`` along a boundary.
+Supposed to be used in conjunction with [`AnalysisSurface`](@ref)
+which stores the boundary information and semidiscretization.
+
+- `rhoinf::Real`: Free-stream density
+- `uinf::Real`: Free-stream velocity
+- `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
+"""
+function SurfaceFrictionCoefficient(rhoinf, uinf, linf)
+    psi = (zero(rhoinf), zero(rhoinf)) # Dummy value/placeholder, will not be used
+    return SurfaceFrictionCoefficient(ForceState(psi, rhoinf, uinf, linf))
+end
+
+function (pressure_coefficient::SurfacePressureCoefficient)(u, equations)
+    p = pressure(u, equations)
+    @unpack pinf = pressure_coefficient
+    @unpack rhoinf, uinf, linf = pressure_coefficient.force_state
+    return (p - pinf) / (0.5 * rhoinf * uinf^2 * linf)
+end
+
+function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x, t,
+                                                        equations_parabolic,
+                                                        gradients_1, gradients_2)
+    visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
+                                               gradients_1, gradients_2)
+    @unpack psi, rhoinf, uinf, linf = surface_friction.force_state
+
+    # Normalize as `normal_direction` is not necessarily a unit vector
+    n = normal_direction / norm(normal_direction)
+    n_perp = (-n[2], n[1])
+    return (visc_stress_vector[1] * n_perp[1] + visc_stress_vector[2] * n_perp[2]) /
+           (0.5 * rhoinf * uinf^2 * linf)
+end
+
+function analyze(surface_variable::AnalysisSurface, du, u, t,
+                 mesh::P4estMesh{2},
+                 equations, dg::DGSEM, cache)
+    @unpack boundaries = cache
+    @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
+    @unpack weights = dg.basis
+
+    @unpack indices, variable = surface_variable
+
+    dim = ndims(mesh)
+    n_nodes = nnodes(dg)
+    n_elements = length(indices)
+
+    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
+    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+
+    # TODO - Can variable_array be a vector?
+    # TODO - Decide whether to save element mean values too
+
+    index_range = eachnode(dg)
+
+    global_node_index = 1 # Keeps track of solution point number on the surface
+    for boundary in indices
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+        direction = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in index_range
+            u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
+                                         boundary)
+
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            var = variable(u_node, equations)
+
+            coord_array[global_node_index, 1] = x[1]
+            coord_array[global_node_index, 2] = x[2]
+            variable_array[global_node_index, 1] = var
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_index += 1
+        end
+    end
+    # TODO - Sort coord_array, variable_array increasing x?
+    mkpath("out")
+    t_trunc = @sprintf("%.3f", t)
+    filename = varname(variable) * t_trunc * ".txt"
+    # TODO - Should we start with a bigger array and avoid hcat?
+    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
+    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+end
+
+function analyze(surface_variable::AnalysisSurface{Variable},
+                 du, u, t, mesh::P4estMesh{2},
+                 equations, equations_parabolic,
+                 dg::DGSEM, cache,
+                 cache_parabolic) where {Variable <: VariableViscous}
+    @unpack boundaries = cache
+    @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
+    @unpack weights = dg.basis
+
+    @unpack indices, variable = surface_variable
+
+    dim = ndims(mesh)
+    n_nodes = nnodes(dg)
+    n_elements = length(indices)
+
+    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
+    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+
+    # TODO - Should variable_array be a vector?
+    # TODO - Decide whether to save element mean values too
+
+    # Additions for parabolic
+    @unpack viscous_container = cache_parabolic
+    @unpack gradients = viscous_container
+
+    gradients_x, gradients_y = gradients
+
+    surface_integral = zero(eltype(u))
+    index_range = eachnode(dg)
+    global_node_index = 1 # Keeps track of solution point number on the surface
+    for boundary in indices
+        element = boundaries.neighbor_ids[boundary]
+        node_indices = boundaries.node_indices[boundary]
+        direction = indices2direction(node_indices)
+
+        i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
+        j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
+
+        i_node = i_node_start
+        j_node = j_node_start
+        for node_index in index_range
+            u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
+                                         boundary)
+
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            # Extract normal direction at nodes which points from the elements outwards,
+            # i.e., *into* the structure.
+            normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                                    i_node, j_node,
+                                                    element)
+
+            gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg, i_node,
+                                        j_node, element)
+            gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg, i_node,
+                                        j_node, element)
+
+            # Integral over whole boundary surface
+            var = variable(u_node, normal_direction, x, t, equations_parabolic,
+                           gradients_1, gradients_2)
+
+            coord_array[global_node_index, 1] = x[1]
+            coord_array[global_node_index, 2] = x[2]
+            variable_array[global_node_index, 1] = var
+            i_node += i_node_step
+            j_node += j_node_step
+            global_node_index += 1
+        end
+    end
+    # TODO - Sort coord_array, variable_array increasing x?
+    mkpath("out")
+    t_trunc = @sprintf("%.3f", t)
+    filename = varname(variable) * t_trunc * ".txt"
+    # TODO - Should we start with a bigger array and avoid hcat?
+    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
+    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+end
+
+varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
+varname(pressure_coefficient::SurfacePressureCoefficient) = "Cp"
+varname(friction_coefficient::SurfaceFrictionCoefficient) = "Cf"
+
+
+# AnalysisCallback prints to screen and needs these. TODO - Can it be turned off?
+function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
+    "Dummy value"
+end
+function pretty_form_utf(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
+    "Dummy value"
+end
+
+function pretty_form_ascii(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
+    "Dummy value"
+end
+function pretty_form_utf(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
+    "Dummy value"
+end
+
+
+end # muladd

From 9c51d7fbfb8e15964da3d91878c0c0418ee79ac6 Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Fri, 26 Apr 2024 22:46:06 +0530
Subject: [PATCH 18/59] Format

---
 .../elixir_navierstokes_NACA0012airfoil_mach08.jl        | 9 +++++----
 src/Trixi.jl                                             | 7 ++++---
 src/callbacks_step/analysis.jl                           | 3 ++-
 src/callbacks_step/analysis_surface_2d.jl                | 9 ++++-----
 4 files changed, 15 insertions(+), 13 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 0659ebeeb5b..94ee43b3cbd 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -141,17 +141,18 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-
 # TODO - Should these two be moved to a different callback to avoid printing "Dummy Value" so often?
 # This will be fixed if the "Dummy Callback" is somehow not printed at all.
 # Can this be done by overloading on `pretty_form_utf`?
 friction_coefficient = AnalysisSurface(semi, force_boundary_names,
-                                       SurfaceFrictionCoefficient(rho_inf(), u_inf(equations), l_inf()))
+                                       SurfaceFrictionCoefficient(rho_inf(),
+                                                                  u_inf(equations),
+                                                                  l_inf()))
 
 pressure_coefficient = AnalysisSurface(semi, force_boundary_names,
                                        SurfacePressureCoefficient(p_inf(), rho_inf(),
-                                       u_inf(equations),
-                                       l_inf()))
+                                                                  u_inf(equations),
+                                                                  l_inf()))
 
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 80b9c38b850..88bb25af004 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -264,9 +264,10 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
-       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient, SurfaceFrictionCoefficient
-
+       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure,
+       LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient,
+       SurfaceFrictionCoefficient
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 1a0e7de14da..a0fae2b6a85 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -699,7 +699,8 @@ include("analysis_dg3d_parallel.jl")
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
 # Note that this needs to be included after `analysis_surface_integral_2d.jl` to
 # have `VariableViscous` available.
-function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable}, AnalysisSurface{Variable}},
+function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
+                                 AnalysisSurface{Variable}},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
                                                                      Variable <:
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 3388091ce84..5a6b887a596 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -157,7 +157,8 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
+                                element)
             var = variable(u_node, equations)
 
             coord_array[global_node_index, 1] = x[1]
@@ -221,7 +222,8 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
+                                element)
             # Extract normal direction at nodes which points from the elements outwards,
             # i.e., *into* the structure.
             normal_direction = get_normal_direction(direction, contravariant_vectors,
@@ -258,7 +260,6 @@ varname(::Any) = @assert false "Surface variable name not assigned" # This makes
 varname(pressure_coefficient::SurfacePressureCoefficient) = "Cp"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "Cf"
 
-
 # AnalysisCallback prints to screen and needs these. TODO - Can it be turned off?
 function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
     "Dummy value"
@@ -273,6 +274,4 @@ end
 function pretty_form_utf(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
     "Dummy value"
 end
-
-
 end # muladd

From 022dc2829d01f95594787783790d32e1872af356 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Tue, 30 Apr 2024 09:50:27 +0200
Subject: [PATCH 19/59] Apply suggestions from code review

---
 src/callbacks_step/analysis_surface_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 5a6b887a596..6062affeaf4 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -6,7 +6,7 @@
 #! format: noindent
 
 # This file contains callbacks that are performed on the surface like computation of
-# surface forces
+# pointwise surface forces.
 
 """
     AnalysisSurface{Semidiscretization, Variable}(semi,

From 84190880d211053f8df4e301663c4f9be0f930f6 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 30 Apr 2024 10:59:50 +0200
Subject: [PATCH 20/59] Pointwise analysis

---
 ...xir_navierstokes_NACA0012airfoil_mach08.jl |  4 +-
 src/callbacks_step/analysis.jl                | 39 ++++++++++-
 src/callbacks_step/analysis_surface_2d.jl     | 69 +++++++++----------
 .../analysis_surface_integral_2d.jl           | 28 ++++----
 4 files changed, 87 insertions(+), 53 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 94ee43b3cbd..477c58cd8a0 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -161,8 +161,8 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      analysis_integrals = (drag_coefficient,
                                                            lift_coefficient,
                                                            drag_coefficient_shear_force,
-                                                           lift_coefficient_shear_force,
-                                                           friction_coefficient,
+                                                           lift_coefficient_shear_force),
+                                     analysis_pointwise = (friction_coefficient,
                                                            pressure_coefficient))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index a0fae2b6a85..5a9cc964c9b 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -43,7 +43,8 @@ evaluating the computational performance, such as the total runtime, the perform
 (time/DOF/rhs!), the time spent in garbage collection (GC), or the current memory usage (alloc'd
 memory).
 """
-mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegrals,
+mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, AnalysisPointwise,
+                                InitialStateIntegrals,
                                 Cache}
     start_time::Float64
     start_time_last_analysis::Float64
@@ -56,6 +57,7 @@ mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegra
     analyzer::Analyzer
     analysis_errors::Vector{Symbol}
     analysis_integrals::AnalysisIntegrals
+    analysis_pointwise::AnalysisPointwise
     initial_state_integrals::InitialStateIntegrals
     cache::Cache
 end
@@ -80,6 +82,9 @@ function Base.show(io::IO, ::MIME"text/plain",
         for (idx, integral) in enumerate(analysis_callback.analysis_integrals)
             push!(setup, "│ integral " * string(idx) => integral)
         end
+        for (idx, quantity) in enumerate(analysis_callback.analysis_pointwise)
+            push!(setup, "│ quantity " * string(idx) => quantity)
+        end
         push!(setup,
               "save analysis to file" => analysis_callback.save_analysis ? "yes" : "no")
         if analysis_callback.save_analysis
@@ -109,6 +114,7 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                           extra_analysis_integrals = (),
                           analysis_integrals = union(default_analysis_integrals(equations),
                                                      extra_analysis_integrals),
+                          analysis_pointwise = (),
                           RealT = real(solver),
                           uEltype = eltype(cache.elements),
                           kwargs...)
@@ -132,6 +138,7 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                                          analysis_filename,
                                          analyzer,
                                          analysis_errors, Tuple(analysis_integrals),
+                                         Tuple(analysis_pointwise),
                                          SVector(ntuple(_ -> zero(uEltype),
                                                         Val(nvariables(equations)))),
                                          cache_analysis)
@@ -158,7 +165,7 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
 
     analysis_callback = cb.affect!
     analysis_callback.initial_state_integrals = initial_state_integrals
-    @unpack save_analysis, output_directory, analysis_filename, analysis_errors, analysis_integrals = analysis_callback
+    @unpack save_analysis, output_directory, analysis_filename, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
 
     if save_analysis && mpi_isroot()
         mkpath(output_directory)
@@ -203,6 +210,10 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
                 @printf(io, "   %-14s", pretty_form_ascii(quantity))
             end
 
+            for quantity in analysis_pointwise
+                @printf(io, "   %-14s", pretty_form_ascii(quantity))
+            end
+
             println(io)
         end
     end
@@ -368,7 +379,7 @@ end
 # This method is just called internally from `(analysis_callback::AnalysisCallback)(integrator)`
 # and serves as a function barrier. Additionally, it makes the code easier to profile and optimize.
 function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
-    @unpack analyzer, analysis_errors, analysis_integrals = analysis_callback
+    @unpack analyzer, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
     cache_analysis = analysis_callback.cache
     _, equations, _, _ = mesh_equations_solver_cache(semi)
 
@@ -486,6 +497,8 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
 
     # additional integrals
     analyze_integrals(analysis_integrals, io, du, u, t, semi)
+    # additional pointwise quantities
+    analyze_pointwise(analysis_pointwise, du, u, t, semi)
 
     return nothing
 end
@@ -606,6 +619,26 @@ function analyze_integrals(analysis_integrals::Tuple{}, io, du, u, t, semi)
     nothing
 end
 
+# Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
+function analyze_pointwise(analysis_quantities::NTuple{N, Any}, du, u, t,
+                           semi) where {N}
+
+    # Extract the first analysis integral and process it; keep the remaining to be processed later
+    quantity = first(analysis_quantities)
+    remaining_quantities = Base.tail(analysis_quantities)
+
+    analyze(quantity, du, u, t, semi)
+
+    # Recursively call this method with the unprocessed integrals
+    analyze_pointwise(remaining_quantities, du, u, t, semi)
+    return nothing
+end
+
+# terminate the type-stable iteration over tuples
+function analyze_pointwise(analysis_quantities::Tuple{}, du, u, t, semi)
+    nothing
+end
+
 # used for error checks and EOC analysis
 function (cb::DiscreteCallback{Condition, Affect!})(sol) where {Condition,
                                                                 Affect! <:
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 6062affeaf4..24cf626f8e4 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -48,13 +48,19 @@ struct AnalysisSurface{Variable}
     end
 end
 
+struct FlowState{RealT <: Real}
+    rhoinf::RealT
+    uinf::RealT
+    linf::RealT
+end
+
 struct SurfacePressureCoefficient{RealT <: Real}
     pinf::RealT # Free stream pressure
-    force_state::ForceState{RealT}
+    flow_state::FlowState{RealT}
 end
 
 struct SurfaceFrictionCoefficient{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowState{RealT}
 end
 
 """
@@ -75,8 +81,7 @@ which stores the boundary information and semidiscretization.
 - `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
 """
 function SurfacePressureCoefficient(pinf, rhoinf, uinf, linf)
-    psi = (zero(pinf), zero(pinf)) # Dummy value/placeholder, will not be used
-    return SurfacePressureCoefficient(pinf, ForceState(psi, rhoinf, uinf, linf))
+    return SurfacePressureCoefficient(pinf, FlowState(rhoinf, uinf, linf))
 end
 
 """
@@ -96,14 +101,13 @@ which stores the boundary information and semidiscretization.
 - `linf::Real`: Reference length of geometry (e.g. airfoil chord length)
 """
 function SurfaceFrictionCoefficient(rhoinf, uinf, linf)
-    psi = (zero(rhoinf), zero(rhoinf)) # Dummy value/placeholder, will not be used
-    return SurfaceFrictionCoefficient(ForceState(psi, rhoinf, uinf, linf))
+    return SurfaceFrictionCoefficient(FlowState(rhoinf, uinf, linf))
 end
 
 function (pressure_coefficient::SurfacePressureCoefficient)(u, equations)
     p = pressure(u, equations)
     @unpack pinf = pressure_coefficient
-    @unpack rhoinf, uinf, linf = pressure_coefficient.force_state
+    @unpack rhoinf, uinf, linf = pressure_coefficient.flow_state
     return (p - pinf) / (0.5 * rhoinf * uinf^2 * linf)
 end
 
@@ -112,7 +116,7 @@ function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x,
                                                         gradients_1, gradients_2)
     visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = surface_friction.force_state
+    @unpack rhoinf, uinf, linf = surface_friction.flow_state
 
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = normal_direction / norm(normal_direction)
@@ -134,10 +138,9 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     n_nodes = nnodes(dg)
     n_elements = length(indices)
 
-    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
-    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+    coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Can variable_array be a vector?
     # TODO - Decide whether to save element mean values too
 
     index_range = eachnode(dg)
@@ -161,21 +164,20 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
                                 element)
             var = variable(u_node, equations)
 
-            coord_array[global_node_index, 1] = x[1]
-            coord_array[global_node_index, 2] = x[2]
-            variable_array[global_node_index, 1] = var
+            coords[global_node_index, 1] = x[1]
+            coords[global_node_index, 2] = x[2]
+            variables[global_node_index] = var
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coord_array, variable_array increasing x?
+    # TODO - Sort coords, variables increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * t_trunc * ".txt"
+    filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
-    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+    writedlm(joinpath("out", filename), hcat(coords, variables))
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
@@ -193,10 +195,9 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     n_nodes = nnodes(dg)
     n_elements = length(indices)
 
-    coord_array = Array{real(dg)}(undef, n_elements * n_nodes, dim)   # physical coordinates of indices
-    variable_array = Array{real(dg)}(undef, n_elements * n_nodes, 1) # variable values at indices
+    coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Should variable_array be a vector?
     # TODO - Decide whether to save element mean values too
 
     # Additions for parabolic
@@ -239,39 +240,37 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             var = variable(u_node, normal_direction, x, t, equations_parabolic,
                            gradients_1, gradients_2)
 
-            coord_array[global_node_index, 1] = x[1]
-            coord_array[global_node_index, 2] = x[2]
-            variable_array[global_node_index, 1] = var
+            coords[global_node_index, 1] = x[1]
+            coords[global_node_index, 2] = x[2]
+            variables[global_node_index] = var
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coord_array, variable_array increasing x?
+    # TODO - Sort coords, variables increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * t_trunc * ".txt"
+    filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coord_array, variable_array))
-    return 0.0 # A dummy value is needed by an analysis callback. TODO - Can this be avoided?
+    writedlm(joinpath("out", filename), hcat(coords, variables))
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
-varname(pressure_coefficient::SurfacePressureCoefficient) = "Cp"
-varname(friction_coefficient::SurfaceFrictionCoefficient) = "Cf"
+varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
+varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
-# AnalysisCallback prints to screen and needs these. TODO - Can it be turned off?
 function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
-    "Dummy value"
+    "CP(x)"
 end
 function pretty_form_utf(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
-    "Dummy value"
+    "CP(x)"
 end
 
 function pretty_form_ascii(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
-    "Dummy value"
+    "CF(x)"
 end
 function pretty_form_utf(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
-    "Dummy value"
+    "CF(x)"
 end
 end # muladd
diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index 7ae259e5285..e58769f3616 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -48,7 +48,7 @@ struct AnalysisSurfaceIntegral{Variable}
     end
 end
 
-struct ForceState{RealT <: Real}
+struct FlowStateDirectional{RealT <: Real}
     psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream
     rhoinf::RealT
     uinf::RealT
@@ -56,11 +56,11 @@ struct ForceState{RealT <: Real}
 end
 
 struct LiftCoefficientPressure{RealT <: Real}
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 struct DragCoefficientPressure{RealT <: Real}
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 # Abstract base type used for dispatch of `analyze` for quantities
@@ -68,11 +68,11 @@ end
 abstract type VariableViscous end
 
 struct LiftCoefficientShearStress{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 struct DragCoefficientShearStress{RealT <: Real} <: VariableViscous
-    force_state::ForceState{RealT}
+    flow_state::FlowStateDirectional{RealT}
 end
 
 """
@@ -99,7 +99,7 @@ function LiftCoefficientPressure(aoa, rhoinf, uinf, linf)
     # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same
     # value, but with the opposite sign.
     psi_lift = (-sin(aoa), cos(aoa))
-    return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, linf))
+    return LiftCoefficientPressure(FlowStateDirectional(psi_lift, rhoinf, uinf, linf))
 end
 
 """
@@ -122,7 +122,7 @@ which stores the boundary information and semidiscretization.
 function DragCoefficientPressure(aoa, rhoinf, uinf, linf)
     # `psi_drag` is the unit vector tangent to the freestream direction
     psi_drag = (cos(aoa), sin(aoa))
-    return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf))
+    return DragCoefficientPressure(FlowStateDirectional(psi_drag, rhoinf, uinf, linf))
 end
 
 """
@@ -149,7 +149,8 @@ function LiftCoefficientShearStress(aoa, rhoinf, uinf, linf)
     # One could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same
     # value, but with the opposite sign.
     psi_lift = (-sin(aoa), cos(aoa))
-    return LiftCoefficientShearStress(ForceState(psi_lift, rhoinf, uinf, linf))
+    return LiftCoefficientShearStress(FlowStateDirectional(psi_lift, rhoinf, uinf,
+                                                           linf))
 end
 
 """
@@ -172,13 +173,14 @@ which stores the boundary information and semidiscretization.
 function DragCoefficientShearStress(aoa, rhoinf, uinf, linf)
     # `psi_drag` is the unit vector tangent to the freestream direction
     psi_drag = (cos(aoa), sin(aoa))
-    return DragCoefficientShearStress(ForceState(psi_drag, rhoinf, uinf, linf))
+    return DragCoefficientShearStress(FlowStateDirectional(psi_drag, rhoinf, uinf,
+                                                           linf))
 end
 
 function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, x, t,
                                                      equations)
     p = pressure(u, equations)
-    @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = dot(normal_direction, psi) / norm(normal_direction)
     return p * n / (0.5 * rhoinf * uinf^2 * linf)
@@ -187,7 +189,7 @@ end
 function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, x, t,
                                                      equations)
     p = pressure(u, equations)
-    @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = dot(normal_direction, psi) / norm(normal_direction)
     return p * n / (0.5 * rhoinf * uinf^2 * linf)
@@ -240,7 +242,7 @@ function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, x,
                                                         gradients_1, gradients_2)
     visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
     return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
@@ -250,7 +252,7 @@ function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x,
                                                         gradients_1, gradients_2)
     visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                gradients_1, gradients_2)
-    @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state
+    @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
     return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end

From 10cfb003df2eab30d85059160c28326fa88f7c5a Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 3 May 2024 11:17:43 +0200
Subject: [PATCH 21/59] cosmetics

---
 src/callbacks_step/analysis_surface_2d.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 24cf626f8e4..1605149ad94 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -10,8 +10,8 @@
 
 """
     AnalysisSurface{Semidiscretization, Variable}(semi,
-                                                          boundary_symbol_or_boundary_symbols,
-                                                          variable)
+                                                  boundary_symbol_or_boundary_symbols,
+                                                  variable)
 
 This struct is used to compute pointwise surface values of a quantity of interest `variable`
 alongside the boundary/boundaries associated with particular name(s) given in `boundary_symbol`

From 57d2b380a1f7b8f77f9ebf6f3903a0f018dca259 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 3 May 2024 13:52:32 +0200
Subject: [PATCH 22/59] routines similar to timeseries

---
 src/callbacks_step/analysis_surface_2d.jl | 51 ++++++++++++++++-------
 1 file changed, 36 insertions(+), 15 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 1605149ad94..5b407602e49 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -28,15 +28,19 @@ name `:Airfoil` in 2D.
 struct AnalysisSurface{Variable}
     indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
     variable::Variable # Quantity of interest, like lift or drag
+    interval::Int
+    output_directory::String
 
-    function AnalysisSurface(semi, boundary_symbol, variable)
+    function AnalysisSurface(semi, boundary_symbol, variable, interval, 
+                             output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = boundary_symbol_indices[boundary_symbol]
 
-        return new{typeof(variable)}(indices, variable)
+        return new{typeof(variable)}(indices, variable, interval, output_directory)
     end
 
-    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable)
+    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable, interval,
+                             output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = Vector{Int}()
         for name in boundary_symbols
@@ -44,7 +48,7 @@ struct AnalysisSurface{Variable}
         end
         sort!(indices)
 
-        return new{typeof(variable)}(indices, variable)
+        return new{typeof(variable)}(indices, variable, interval, output_directory)
     end
 end
 
@@ -139,7 +143,7 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     n_elements = length(indices)
 
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
-    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
+    values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
     # TODO - Decide whether to save element mean values too
 
@@ -162,22 +166,22 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
 
             x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
                                 element)
-            var = variable(u_node, equations)
+            value = variable(u_node, equations)
 
             coords[global_node_index, 1] = x[1]
             coords[global_node_index, 2] = x[2]
-            variables[global_node_index] = var
+            values[global_node_index] = value
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, variables increasing x?
+    # TODO - Sort coords, values increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
     filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, variables))
+    writedlm(joinpath("out", filename), hcat(coords, values))
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
@@ -196,7 +200,7 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     n_elements = length(indices)
 
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
-    variables = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
+    values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
     # TODO - Decide whether to save element mean values too
 
@@ -237,29 +241,46 @@ function analyze(surface_variable::AnalysisSurface{Variable},
                                         j_node, element)
 
             # Integral over whole boundary surface
-            var = variable(u_node, normal_direction, x, t, equations_parabolic,
-                           gradients_1, gradients_2)
+            value = variable(u_node, normal_direction, x, t, equations_parabolic,
+                             gradients_1, gradients_2)
 
             coords[global_node_index, 1] = x[1]
             coords[global_node_index, 2] = x[2]
-            variables[global_node_index] = var
+            values[global_node_index] = value
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, variables increasing x?
+    # TODO - Sort coords, values increasing x?
     mkpath("out")
     t_trunc = @sprintf("%.3f", t)
     filename = varname(variable) * "_" * t_trunc * ".txt"
     # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, variables))
+    writedlm(joinpath("out", filename), hcat(coords, values))
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
 varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
+function save_pointwise_file(output_directory, varname, coords, values, t)
+    n_points = length(values)
+
+    h5open(joinpath(output_directory, varname) * ".h5", "w") do file
+        # Add context information as attributes
+        attributes(file)["n_points"] = n_points
+        attributes(file)["interval"] = interval
+        attributes(file)["variable_name"] = collect(varname)
+
+        file["time"] = t
+        file["point_coordinates"] = coords
+        for p in 1:n_points
+            file["point_data_$p"] = values[p]
+        end
+    end
+end
+
 function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
     "CP(x)"
 end

From b995e3fc4d448aeda9a12b7ff25790e9c015ff5c Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:50:40 +0200
Subject: [PATCH 23/59] write pointwise data to jhdf5

---
 src/callbacks_step/analysis.jl            | 32 ++++++++++++++-----
 src/callbacks_step/analysis_surface_2d.jl | 39 +++++++++--------------
 2 files changed, 39 insertions(+), 32 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 5a9cc964c9b..60f828cf482 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -350,8 +350,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
         @notimeit timer() integrator.f(du_ode, u_ode, semi, t)
         u = wrap_array(u_ode, mesh, equations, solver, cache)
         du = wrap_array(du_ode, mesh, equations, solver, cache)
-        # Compute l2_error, linf_error
-        analysis_callback(io, du, u, u_ode, t, semi)
+        # Compute l2_error, linf_error among other quantities
+        analysis_callback(io, du, u, u_ode, t, semi, iter)
 
         mpi_println("─"^100)
         mpi_println()
@@ -378,7 +378,7 @@ end
 
 # This method is just called internally from `(analysis_callback::AnalysisCallback)(integrator)`
 # and serves as a function barrier. Additionally, it makes the code easier to profile and optimize.
-function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
+function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi, iter)
     @unpack analyzer, analysis_errors, analysis_integrals, analysis_pointwise = analysis_callback
     cache_analysis = analysis_callback.cache
     _, equations, _, _ = mesh_equations_solver_cache(semi)
@@ -498,7 +498,7 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
     # additional integrals
     analyze_integrals(analysis_integrals, io, du, u, t, semi)
     # additional pointwise quantities
-    analyze_pointwise(analysis_pointwise, du, u, t, semi)
+    analyze_pointwise(analysis_pointwise, du, u, t, semi, iter)
 
     return nothing
 end
@@ -621,21 +621,21 @@ end
 
 # Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
 function analyze_pointwise(analysis_quantities::NTuple{N, Any}, du, u, t,
-                           semi) where {N}
+                           semi, iter) where {N}
 
     # Extract the first analysis integral and process it; keep the remaining to be processed later
     quantity = first(analysis_quantities)
     remaining_quantities = Base.tail(analysis_quantities)
 
-    analyze(quantity, du, u, t, semi)
+    analyze(quantity, du, u, t, semi, iter)
 
     # Recursively call this method with the unprocessed integrals
-    analyze_pointwise(remaining_quantities, du, u, t, semi)
+    analyze_pointwise(remaining_quantities, du, u, t, semi, iter)
     return nothing
 end
 
 # terminate the type-stable iteration over tuples
-function analyze_pointwise(analysis_quantities::Tuple{}, du, u, t, semi)
+function analyze_pointwise(analysis_quantities::Tuple{}, du, u, t, semi, iter)
     nothing
 end
 
@@ -660,6 +660,10 @@ function analyze(quantity, du, u, t, semi::AbstractSemidiscretization)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache)
 end
+function analyze(quantity, du, u, t, semi::AbstractSemidiscretization, iter)
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    analyze(quantity, du, u, t, mesh, equations, solver, cache, iter)
+end
 function analyze(quantity, du, u, t, mesh, equations, solver, cache)
     integrate(quantity, u, mesh, equations, solver, cache, normalize = true)
 end
@@ -744,3 +748,15 @@ function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
     analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
             cache_parabolic)
 end
+function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
+                                 AnalysisSurface{Variable}},
+                 du, u, t,
+                 semi::SemidiscretizationHyperbolicParabolic, iter) where {
+                                                                     Variable <:
+                                                                     VariableViscous}
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    equations_parabolic = semi.equations_parabolic
+    cache_parabolic = semi.cache_parabolic
+    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
+            cache_parabolic, iter)
+end
\ No newline at end of file
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 5b407602e49..b5c6e8266f9 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -28,18 +28,17 @@ name `:Airfoil` in 2D.
 struct AnalysisSurface{Variable}
     indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
     variable::Variable # Quantity of interest, like lift or drag
-    interval::Int
     output_directory::String
 
-    function AnalysisSurface(semi, boundary_symbol, variable, interval, 
+    function AnalysisSurface(semi, boundary_symbol, variable, 
                              output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = boundary_symbol_indices[boundary_symbol]
 
-        return new{typeof(variable)}(indices, variable, interval, output_directory)
+        return new{typeof(variable)}(indices, variable, output_directory)
     end
 
-    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable, interval,
+    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable,
                              output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = Vector{Int}()
@@ -48,7 +47,7 @@ struct AnalysisSurface{Variable}
         end
         sort!(indices)
 
-        return new{typeof(variable)}(indices, variable, interval, output_directory)
+        return new{typeof(variable)}(indices, variable, output_directory)
     end
 end
 
@@ -131,7 +130,7 @@ end
 
 function analyze(surface_variable::AnalysisSurface, du, u, t,
                  mesh::P4estMesh{2},
-                 equations, dg::DGSEM, cache)
+                 equations, dg::DGSEM, cache, iter)
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
@@ -177,18 +176,14 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
         end
     end
     # TODO - Sort coords, values increasing x?
-    mkpath("out")
-    t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * "_" * t_trunc * ".txt"
-    # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, values))
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
                  du, u, t, mesh::P4estMesh{2},
                  equations, equations_parabolic,
                  dg::DGSEM, cache,
-                 cache_parabolic) where {Variable <: VariableViscous}
+                 cache_parabolic, iter) where {Variable <: VariableViscous}
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
@@ -253,31 +248,27 @@ function analyze(surface_variable::AnalysisSurface{Variable},
         end
     end
     # TODO - Sort coords, values increasing x?
-    mkpath("out")
-    t_trunc = @sprintf("%.3f", t)
-    filename = varname(variable) * "_" * t_trunc * ".txt"
-    # TODO - Should we start with a bigger array and avoid hcat?
-    writedlm(joinpath("out", filename), hcat(coords, values))
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
 varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
-function save_pointwise_file(output_directory, varname, coords, values, t)
+function save_pointwise_file(output_directory, varname, coords, values, t, iter)
     n_points = length(values)
 
-    h5open(joinpath(output_directory, varname) * ".h5", "w") do file
+    filename = joinpath(output_directory, varname) * @sprintf("_%06d.h5", iter)
+
+    h5open(filename, "w") do file
         # Add context information as attributes
         attributes(file)["n_points"] = n_points
-        attributes(file)["interval"] = interval
-        attributes(file)["variable_name"] = collect(varname)
+        attributes(file)["variable_name"] = varname
 
         file["time"] = t
+        file["timestep"] = iter
         file["point_coordinates"] = coords
-        for p in 1:n_points
-            file["point_data_$p"] = values[p]
-        end
+        file["values"] = values
     end
 end
 

From cb826db39d2cb81a29f21684802a4b6681b565c2 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:51:24 +0200
Subject: [PATCH 24/59] fmt

---
 src/callbacks_step/analysis.jl            | 9 +++++----
 src/callbacks_step/analysis_surface_2d.jl | 8 +++++---
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 60f828cf482..fa1fc0782bd 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -751,12 +751,13 @@ end
 function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
                                  AnalysisSurface{Variable}},
                  du, u, t,
-                 semi::SemidiscretizationHyperbolicParabolic, iter) where {
-                                                                     Variable <:
-                                                                     VariableViscous}
+                 semi::SemidiscretizationHyperbolicParabolic,
+                 iter) where {
+                              Variable <:
+                              VariableViscous}
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     equations_parabolic = semi.equations_parabolic
     cache_parabolic = semi.cache_parabolic
     analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
             cache_parabolic, iter)
-end
\ No newline at end of file
+end
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index b5c6e8266f9..62f9a97c002 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -30,7 +30,7 @@ struct AnalysisSurface{Variable}
     variable::Variable # Quantity of interest, like lift or drag
     output_directory::String
 
-    function AnalysisSurface(semi, boundary_symbol, variable, 
+    function AnalysisSurface(semi, boundary_symbol, variable,
                              output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = boundary_symbol_indices[boundary_symbol]
@@ -176,7 +176,8 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
         end
     end
     # TODO - Sort coords, values increasing x?
-    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
+                        values, t, iter)
 end
 
 function analyze(surface_variable::AnalysisSurface{Variable},
@@ -248,7 +249,8 @@ function analyze(surface_variable::AnalysisSurface{Variable},
         end
     end
     # TODO - Sort coords, values increasing x?
-    save_pointwise_file(surface_variable.output_directory, varname(variable), coords, values, t, iter)
+    save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
+                        values, t, iter)
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting

From f39d7f372371041bb8b8848feeb742ffbaa15d2b Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:52:21 +0200
Subject: [PATCH 25/59] not needed

---
 Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index def3691a9e8..769c1fe7d08 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ version = "0.7.13-pre"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

From cf2fb11e89cdca11c76d452e8fa8c3767c547ada Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:52:54 +0200
Subject: [PATCH 26/59] rm unnecessary

---
 src/Trixi.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index 88bb25af004..cc5eaac7510 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -42,7 +42,6 @@ import SciMLBase: get_du, get_tmp_cache, u_modified!,
 using Downloads: Downloads
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
-using DelimitedFiles: writedlm
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
 @reexport using EllipsisNotation # ..
 using FillArrays: Ones, Zeros

From 2bc600184cfc497841fb9510ffd3aa9b8ef50939 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:54:01 +0200
Subject: [PATCH 27/59] re-organize

---
 src/Trixi.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index cc5eaac7510..802c1c46e30 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -263,10 +263,10 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurface, AnalysisSurfaceIntegral, DragCoefficientPressure,
-       LiftCoefficientPressure,
-       DragCoefficientShearStress, LiftCoefficientShearStress, SurfacePressureCoefficient,
-       SurfaceFrictionCoefficient
+       AnalysisSurface, AnalysisSurfaceIntegral,
+       DragCoefficientPressure, LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress,
+       SurfacePressureCoefficient, SurfaceFrictionCoefficient
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!

From f8b247345f8319469b4dd3261ae70fd8ed6a6365 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Sun, 5 May 2024 12:58:18 +0200
Subject: [PATCH 28/59] remove todo

---
 .../elixir_navierstokes_NACA0012airfoil_mach08.jl              | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 477c58cd8a0..ddbd3dc3652 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -141,9 +141,6 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-# TODO - Should these two be moved to a different callback to avoid printing "Dummy Value" so often?
-# This will be fixed if the "Dummy Callback" is somehow not printed at all.
-# Can this be done by overloading on `pretty_form_utf`?
 friction_coefficient = AnalysisSurface(semi, force_boundary_names,
                                        SurfaceFrictionCoefficient(rho_inf(),
                                                                   u_inf(equations),

From cd18b8367ea688212d318e21d20fed640bc31b0e Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 8 May 2024 15:33:54 +0200
Subject: [PATCH 29/59] Sort

---
 src/callbacks_step/analysis_surface_2d.jl | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 62f9a97c002..4f5fbc87b89 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -144,8 +144,6 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
     values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Decide whether to save element mean values too
-
     index_range = eachnode(dg)
 
     global_node_index = 1 # Keeps track of solution point number on the surface
@@ -175,7 +173,11 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, values increasing x?
+    # Sort with ascending x
+    sorting_indices = sortperm(coords[:, 1])
+    coords = coords[sorting_indices, :]
+    values = values[sorting_indices]
+
     save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
                         values, t, iter)
 end
@@ -198,8 +200,6 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
     values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
-    # TODO - Decide whether to save element mean values too
-
     # Additions for parabolic
     @unpack viscous_container = cache_parabolic
     @unpack gradients = viscous_container
@@ -248,7 +248,14 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             global_node_index += 1
         end
     end
-    # TODO - Sort coords, values increasing x?
+    # Sort with ascending x
+    sorting_indices = sortperm(coords[:, 1])
+    coords = coords[sorting_indices, :]
+    values = values[sorting_indices]
+
+    println(size(coords))
+    println(size(values))
+
     save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
                         values, t, iter)
 end

From bdb82f1f54a91b9ea275e15b5c0aa1336df9f7bb Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 8 May 2024 15:34:30 +0200
Subject: [PATCH 30/59] remove prnt

---
 src/callbacks_step/analysis_surface_2d.jl | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 4f5fbc87b89..14211c9521b 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -253,9 +253,6 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     coords = coords[sorting_indices, :]
     values = values[sorting_indices]
 
-    println(size(coords))
-    println(size(values))
-
     save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
                         values, t, iter)
 end

From 5f5f4644bc7d647e809a6f45ca369f96adc352bf Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Tue, 7 May 2024 08:55:43 +0200
Subject: [PATCH 31/59] Apply suggestions from code review

Co-authored-by: Arpit Babbar <arpitbabbar@gmail.com>
---
 src/callbacks_step/analysis.jl            | 6 +++---
 src/callbacks_step/analysis_surface_2d.jl | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index fa1fc0782bd..fcf677dda3a 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -619,17 +619,17 @@ function analyze_integrals(analysis_integrals::Tuple{}, io, du, u, t, semi)
     nothing
 end
 
-# Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
+# Iterate over tuples of pointwise analysis quantities in a type-stable way using "lispy tuple programming".
 function analyze_pointwise(analysis_quantities::NTuple{N, Any}, du, u, t,
                            semi, iter) where {N}
 
-    # Extract the first analysis integral and process it; keep the remaining to be processed later
+    # Extract the first pointwise analysis quantity and process it; keep the remaining to be processed later
     quantity = first(analysis_quantities)
     remaining_quantities = Base.tail(analysis_quantities)
 
     analyze(quantity, du, u, t, semi, iter)
 
-    # Recursively call this method with the unprocessed integrals
+    # Recursively call this method with the unprocessed pointwise analysis quantities
     analyze_pointwise(remaining_quantities, du, u, t, semi, iter)
     return nothing
 end
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 14211c9521b..35749d010f4 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -21,7 +21,7 @@ surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an air
 name `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
-- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries
+- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries stored as symbol(s)
   where the quantity of interest is computed
 - `variable::Variable`: Quantity of interest, like lift or drag
 """

From db33b6437d51c7a5b976c9067703a9a147118bd2 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 8 May 2024 15:36:24 +0200
Subject: [PATCH 32/59] name/symbol

---
 src/callbacks_step/analysis_surface_2d.jl          | 2 +-
 src/callbacks_step/analysis_surface_integral_2d.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 35749d010f4..ee09866b323 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -18,7 +18,7 @@ alongside the boundary/boundaries associated with particular name(s) given in `b
 or `boundary_symbols`.
 For instance, this can be used to compute the surface pressure coefficient [`SurfacePressureCoefficient`](@ref) or
 surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an airfoil with the boundary
-name `:Airfoil` in 2D.
+symbol `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
 - `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries stored as symbol(s)
diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index e58769f3616..db1ad0182c1 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -18,7 +18,7 @@ the boundary/boundaries associated with particular name(s) given in `boundary_sy
 or `boundary_symbols`.
 For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or
 drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary
-name `:Airfoil` in 2D.
+symbol `:Airfoil` in 2D.
 
 - `semi::Semidiscretization`: Passed in to retrieve boundary condition information
 - `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries

From b9021a43d6a42dc9796b2f7e18260bb5cdcb0116 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Thu, 9 May 2024 09:36:31 +0200
Subject: [PATCH 33/59] Update src/callbacks_step/analysis_surface_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/analysis_surface_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index ee09866b323..928012f68fa 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -274,7 +274,7 @@ function save_pointwise_file(output_directory, varname, coords, values, t, iter)
         file["time"] = t
         file["timestep"] = iter
         file["point_coordinates"] = coords
-        file["values"] = values
+        file["point_data"] = values
     end
 end
 

From bac4e4d320d5e2cdb5f4449b27b40685ecdc5748 Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Mon, 13 May 2024 14:20:10 +0530
Subject: [PATCH 34/59] Add test

---
 .../elixir_advection_diffusion_nonperiodic.jl |  8 ++---
 test/test_parabolic_2d.jl                     | 31 ++++++++++++++++++-
 2 files changed, 34 insertions(+), 5 deletions(-)

diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl
index 8da542b0b5d..fbd95c1565f 100644
--- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl
+++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl
@@ -7,7 +7,7 @@ using Trixi
 diffusivity() = 5.0e-2
 advection_velocity = (1.0, 0.0)
 equations = LinearScalarAdvectionEquation2D(advection_velocity)
-equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+equations_parabolic = LaplaceDiffusion2D(5.0e-2, equations)
 
 # 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)
@@ -17,7 +17,7 @@ coordinates_max = (0.0, 0.5) # maximum coordinates (max(x), max(y))
 
 # Create a uniformly refined mesh with periodic boundaries
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 4,
+                initial_refinement_level = 3,
                 periodicity = false,
                 n_cells_max = 30_000) # set maximum capacity of tree data structure
 
@@ -29,7 +29,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 #   [DOI](https://doi.org/10.1007/978-3-319-41640-3_6).
 function initial_condition_eriksson_johnson(x, t, equations)
     l = 4
-    epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt
+    epsilon = 5.0e-2 # TODO: this requires epsilon < .6 due to sqrt
     lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
     lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
     r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
@@ -79,7 +79,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
 # run the simulation
 
 # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
-time_int_tol = 1.0e-11
+time_int_tol = 1.0e-13
 sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)
 
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 7749a0c4780..dada007c36a 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -729,7 +729,7 @@ end
         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
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 3000
     end
 
     u_ode = copy(sol.u[end])
@@ -755,6 +755,35 @@ end
 
     @test isapprox(drag_f, 1.5427441885921553, atol = 1e-13)
     @test isapprox(lift_f, 0.005621910087395724, atol = 1e-13)
+
+    # Check cp, cf written to outfiles
+
+    cp_vals = read(Trixi.h5open("out/CP_x_000007.h5"))["point_data"]
+    cf_vals = read(Trixi.h5open("out/CF_x_000007.h5"))["point_data"]
+    @test cp_vals[1:10] ≈ [
+        2.0956920203526193,
+        2.098639791879171,
+        2.0354491762572806,
+        2.1591959225932777,
+        1.7417814818119175,
+        2.0400855926847936,
+        1.874247697759559,
+        1.871513931283643,
+        1.5182860925755697,
+        1.5152278762705806,
+    ]
+    @show cf_vals[1:10] ≈ [
+        0.2038354965561312,
+        0.2634810993266428,
+        0.42919006623618894,
+        0.045526882648679115,
+        0.6611394672752093,
+        -0.23515678100990495,
+        -0.37053727238603884,
+        -0.4537078569805583,
+        0.8113349266282898,
+        0.7458959144321556,
+    ]
 end
 end
 

From 29d2663ed3101845fcbf2e0ae1e5da55000c8108 Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Mon, 13 May 2024 14:23:11 +0530
Subject: [PATCH 35/59] Bug fix

---
 test/test_parabolic_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index dada007c36a..8818979eea2 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -772,7 +772,7 @@ end
         1.5182860925755697,
         1.5152278762705806,
     ]
-    @show cf_vals[1:10] ≈ [
+    @test cf_vals[1:10] ≈ [
         0.2038354965561312,
         0.2634810993266428,
         0.42919006623618894,

From a7b9ec19adb4ed11ea034fe91028ac3b8d5f784d Mon Sep 17 00:00:00 2001
From: Arpit-Babbar <arpitbabbar@gmail.com>
Date: Mon, 13 May 2024 14:25:25 +0530
Subject: [PATCH 36/59] Accidental add fix

---
 .../elixir_advection_diffusion_nonperiodic.jl             | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl
index fbd95c1565f..8da542b0b5d 100644
--- a/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl
+++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_nonperiodic.jl
@@ -7,7 +7,7 @@ using Trixi
 diffusivity() = 5.0e-2
 advection_velocity = (1.0, 0.0)
 equations = LinearScalarAdvectionEquation2D(advection_velocity)
-equations_parabolic = LaplaceDiffusion2D(5.0e-2, equations)
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
 
 # 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)
@@ -17,7 +17,7 @@ coordinates_max = (0.0, 0.5) # maximum coordinates (max(x), max(y))
 
 # Create a uniformly refined mesh with periodic boundaries
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 3,
+                initial_refinement_level = 4,
                 periodicity = false,
                 n_cells_max = 30_000) # set maximum capacity of tree data structure
 
@@ -29,7 +29,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 #   [DOI](https://doi.org/10.1007/978-3-319-41640-3_6).
 function initial_condition_eriksson_johnson(x, t, equations)
     l = 4
-    epsilon = 5.0e-2 # TODO: this requires epsilon < .6 due to sqrt
+    epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt
     lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
     lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon)
     r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon)
@@ -79,7 +79,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
 # run the simulation
 
 # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
-time_int_tol = 1.0e-13
+time_int_tol = 1.0e-11
 sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
             ode_default_options()..., callback = callbacks)
 

From efe8ddf1766e1df47d53a60d5046116612a39685 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Tue, 14 May 2024 09:02:16 +0200
Subject: [PATCH 37/59] Apply suggestions from code review

---
 src/callbacks_step/analysis_surface_2d.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 928012f68fa..5148efae7ea 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -71,8 +71,8 @@ end
 
 Compute the surface pressure coefficient
 ```math
-C_p \\coloneqq \\frac{ p - p_{p_\\infty}}
-                        {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+C_p \\coloneqq \\frac{p - p_{p_\\infty}}
+                     {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
 ```
 based on the pressure distribution along a boundary.
 Supposed to be used in conjunction with [`AnalysisSurface`](@ref)
@@ -92,8 +92,8 @@ SurfaceFrictionCoefficient(rhoinf, uinf, linf)
 
 Compute the surface skin friction coefficient
 ```math
-C_f \\coloneqq \\frac{\\boldsymbol (\\tau_w  \\boldsymbol n) \\cdot \\boldsymbol b^\\perp}
-                        {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
+C_f \\coloneqq \\frac{\\boldsymbol \\tau_w  \\boldsymbol n^\\perp}
+                     {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
 ```
 based on the wall shear stress vector ``\\tau_w`` along a boundary.
 Supposed to be used in conjunction with [`AnalysisSurface`](@ref)

From 4180757d4cdb2f8a19b21b69cc2eef749cea2fe4 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Tue, 21 May 2024 08:20:41 +0200
Subject: [PATCH 38/59] Update
 examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl

---
 .../p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index ddbd3dc3652..a1c74048671 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -1,6 +1,7 @@
 using Downloads: download
 using OrdinaryDiffEq
 using Trixi
+
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 

From 9785a83898943cb34c098f1aeb06548cda5e9126 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 11:08:59 +0200
Subject: [PATCH 39/59] remove sorting

---
 src/callbacks_step/analysis_surface_2d.jl | 39 ++++++++++-------------
 1 file changed, 16 insertions(+), 23 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 5148efae7ea..7f119fbaf00 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -141,7 +141,7 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     n_nodes = nnodes(dg)
     n_elements = length(indices)
 
-    coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    coordinates = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
     values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
     index_range = eachnode(dg)
@@ -150,7 +150,6 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
     for boundary in indices
         element = boundaries.neighbor_ids[boundary]
         node_indices = boundaries.node_indices[boundary]
-        direction = indices2direction(node_indices)
 
         i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
         j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)
@@ -161,24 +160,21 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
-                                element)
+            x = get_node_coordinates(node_coordinates, equations, dg, i_node, j_node,
+                                     element)
             value = variable(u_node, equations)
 
-            coords[global_node_index, 1] = x[1]
-            coords[global_node_index, 2] = x[2]
+            coordinates[global_node_index, 1] = x[1]
+            coordinates[global_node_index, 2] = x[2]
             values[global_node_index] = value
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # Sort with ascending x
-    sorting_indices = sortperm(coords[:, 1])
-    coords = coords[sorting_indices, :]
-    values = values[sorting_indices]
 
-    save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
+    save_pointwise_file(surface_variable.output_directory, varname(variable),
+                        coordinates,
                         values, t, iter)
 end
 
@@ -197,7 +193,7 @@ function analyze(surface_variable::AnalysisSurface{Variable},
     n_nodes = nnodes(dg)
     n_elements = length(indices)
 
-    coords = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
+    coordinates = Matrix{real(dg)}(undef, n_elements * n_nodes, dim) # physical coordinates of indices
     values = Vector{real(dg)}(undef, n_elements * n_nodes) # variable values at indices
 
     # Additions for parabolic
@@ -223,8 +219,8 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
-                                element)
+            x = get_node_coordinates(node_coordinates, equations, dg, i_node, j_node,
+                                     element)
             # Extract normal direction at nodes which points from the elements outwards,
             # i.e., *into* the structure.
             normal_direction = get_normal_direction(direction, contravariant_vectors,
@@ -240,20 +236,17 @@ function analyze(surface_variable::AnalysisSurface{Variable},
             value = variable(u_node, normal_direction, x, t, equations_parabolic,
                              gradients_1, gradients_2)
 
-            coords[global_node_index, 1] = x[1]
-            coords[global_node_index, 2] = x[2]
+            coordinates[global_node_index, 1] = x[1]
+            coordinates[global_node_index, 2] = x[2]
             values[global_node_index] = value
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
         end
     end
-    # Sort with ascending x
-    sorting_indices = sortperm(coords[:, 1])
-    coords = coords[sorting_indices, :]
-    values = values[sorting_indices]
 
-    save_pointwise_file(surface_variable.output_directory, varname(variable), coords,
+    save_pointwise_file(surface_variable.output_directory, varname(variable),
+                        coordinates,
                         values, t, iter)
 end
 
@@ -261,7 +254,7 @@ varname(::Any) = @assert false "Surface variable name not assigned" # This makes
 varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
-function save_pointwise_file(output_directory, varname, coords, values, t, iter)
+function save_pointwise_file(output_directory, varname, coordinates, values, t, iter)
     n_points = length(values)
 
     filename = joinpath(output_directory, varname) * @sprintf("_%06d.h5", iter)
@@ -273,7 +266,7 @@ function save_pointwise_file(output_directory, varname, coords, values, t, iter)
 
         file["time"] = t
         file["timestep"] = iter
-        file["point_coordinates"] = coords
+        file["point_coordinates"] = coordinates
         file["point_data"] = values
     end
 end

From e25e4626d6235127c46df514174cf6178eb7259c Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Fri, 24 May 2024 11:09:22 +0200
Subject: [PATCH 40/59] Update
 examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 .../p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index ddbd3dc3652..a1c74048671 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -1,6 +1,7 @@
 using Downloads: download
 using OrdinaryDiffEq
 using Trixi
+
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 

From f0e7047e51d041699d05ddf9c142ce9fe36ef78b Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Fri, 24 May 2024 11:09:36 +0200
Subject: [PATCH 41/59] Update src/callbacks_step/analysis.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
---
 src/callbacks_step/analysis.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index fcf677dda3a..215d11b281f 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -83,7 +83,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             push!(setup, "│ integral " * string(idx) => integral)
         end
         for (idx, quantity) in enumerate(analysis_callback.analysis_pointwise)
-            push!(setup, "│ quantity " * string(idx) => quantity)
+            push!(setup, "│ pointwise " * string(idx) => quantity)
         end
         push!(setup,
               "save analysis to file" => analysis_callback.save_analysis ? "yes" : "no")

From 012680c9111fa413bf7d1f93c3d3aa37be6d9065 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 11:14:28 +0200
Subject: [PATCH 42/59] pointwise

---
 ...xir_navierstokes_NACA0012airfoil_mach08.jl |  4 +--
 src/Trixi.jl                                  |  2 +-
 src/callbacks_step/analysis.jl                |  4 +--
 src/callbacks_step/analysis_surface_2d.jl     | 28 +++++++++----------
 4 files changed, 19 insertions(+), 19 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index a1c74048671..c76e2a24d0d 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -142,12 +142,12 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-friction_coefficient = AnalysisSurface(semi, force_boundary_names,
+friction_coefficient = AnalysisSurfacePointwise(semi, force_boundary_names,
                                        SurfaceFrictionCoefficient(rho_inf(),
                                                                   u_inf(equations),
                                                                   l_inf()))
 
-pressure_coefficient = AnalysisSurface(semi, force_boundary_names,
+pressure_coefficient = AnalysisSurfacePointwise(semi, force_boundary_names,
                                        SurfacePressureCoefficient(p_inf(), rho_inf(),
                                                                   u_inf(equations),
                                                                   l_inf()))
diff --git a/src/Trixi.jl b/src/Trixi.jl
index dca3934307c..96d4b86142a 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -263,7 +263,7 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurface, AnalysisSurfaceIntegral,
+       AnalysisSurfacePointwise, AnalysisSurfaceIntegral,
        DragCoefficientPressure, LiftCoefficientPressure,
        DragCoefficientShearStress, LiftCoefficientShearStress,
        SurfacePressureCoefficient, SurfaceFrictionCoefficient
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 215d11b281f..046f6327da5 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -737,7 +737,7 @@ include("analysis_dg3d_parallel.jl")
 # Note that this needs to be included after `analysis_surface_integral_2d.jl` to
 # have `VariableViscous` available.
 function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
-                                 AnalysisSurface{Variable}},
+                                 AnalysisSurfacePointwise{Variable}},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
                                                                      Variable <:
@@ -749,7 +749,7 @@ function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
             cache_parabolic)
 end
 function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
-                                 AnalysisSurface{Variable}},
+                                 AnalysisSurfacePointwise{Variable}},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic,
                  iter) where {
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl
index 7f119fbaf00..5e9503051d6 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_2d.jl
@@ -9,7 +9,7 @@
 # pointwise surface forces.
 
 """
-    AnalysisSurface{Semidiscretization, Variable}(semi,
+    AnalysisSurfacePointwise{Semidiscretization, Variable}(semi,
                                                   boundary_symbol_or_boundary_symbols,
                                                   variable)
 
@@ -25,21 +25,21 @@ symbol `:Airfoil` in 2D.
   where the quantity of interest is computed
 - `variable::Variable`: Quantity of interest, like lift or drag
 """
-struct AnalysisSurface{Variable}
+struct AnalysisSurfacePointwise{Variable}
     indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
     variable::Variable # Quantity of interest, like lift or drag
     output_directory::String
 
-    function AnalysisSurface(semi, boundary_symbol, variable,
-                             output_directory = "out")
+    function AnalysisSurfacePointwise(semi, boundary_symbol, variable,
+                                      output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = boundary_symbol_indices[boundary_symbol]
 
         return new{typeof(variable)}(indices, variable, output_directory)
     end
 
-    function AnalysisSurface(semi, boundary_symbols::Vector{Symbol}, variable,
-                             output_directory = "out")
+    function AnalysisSurfacePointwise(semi, boundary_symbols::Vector{Symbol}, variable,
+                                      output_directory = "out")
         @unpack boundary_symbol_indices = semi.boundary_conditions
         indices = Vector{Int}()
         for name in boundary_symbols
@@ -75,7 +75,7 @@ C_p \\coloneqq \\frac{p - p_{p_\\infty}}
                      {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
 ```
 based on the pressure distribution along a boundary.
-Supposed to be used in conjunction with [`AnalysisSurface`](@ref)
+Supposed to be used in conjunction with [`AnalysisSurfacePointwise`](@ref)
 which stores the boundary information and semidiscretization.
 
 - `pinf::Real`: Free-stream pressure
@@ -96,7 +96,7 @@ C_f \\coloneqq \\frac{\\boldsymbol \\tau_w  \\boldsymbol n^\\perp}
                      {0.5 \\rho_{\\infty} U_{\\infty}^2 L_{\\infty}}
 ```
 based on the wall shear stress vector ``\\tau_w`` along a boundary.
-Supposed to be used in conjunction with [`AnalysisSurface`](@ref)
+Supposed to be used in conjunction with [`AnalysisSurfacePointwise`](@ref)
 which stores the boundary information and semidiscretization.
 
 - `rhoinf::Real`: Free-stream density
@@ -128,7 +128,7 @@ function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x,
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
-function analyze(surface_variable::AnalysisSurface, du, u, t,
+function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
                  mesh::P4estMesh{2},
                  equations, dg::DGSEM, cache, iter)
     @unpack boundaries = cache
@@ -178,7 +178,7 @@ function analyze(surface_variable::AnalysisSurface, du, u, t,
                         values, t, iter)
 end
 
-function analyze(surface_variable::AnalysisSurface{Variable},
+function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
                  du, u, t, mesh::P4estMesh{2},
                  equations, equations_parabolic,
                  dg::DGSEM, cache,
@@ -271,17 +271,17 @@ function save_pointwise_file(output_directory, varname, coordinates, values, t,
     end
 end
 
-function pretty_form_ascii(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
+function pretty_form_ascii(::AnalysisSurfacePointwise{<:SurfacePressureCoefficient{<:Any}})
     "CP(x)"
 end
-function pretty_form_utf(::AnalysisSurface{<:SurfacePressureCoefficient{<:Any}})
+function pretty_form_utf(::AnalysisSurfacePointwise{<:SurfacePressureCoefficient{<:Any}})
     "CP(x)"
 end
 
-function pretty_form_ascii(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
+function pretty_form_ascii(::AnalysisSurfacePointwise{<:SurfaceFrictionCoefficient{<:Any}})
     "CF(x)"
 end
-function pretty_form_utf(::AnalysisSurface{<:SurfaceFrictionCoefficient{<:Any}})
+function pretty_form_utf(::AnalysisSurfacePointwise{<:SurfaceFrictionCoefficient{<:Any}})
     "CF(x)"
 end
 end # muladd

From 880674e89466b0eee658a25fab1da2a9bf9378ff Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 11:25:50 +0200
Subject: [PATCH 43/59] rename

---
 src/callbacks_step/analysis.jl                                 | 2 +-
 ...analysis_surface_2d.jl => analysis_surface_pointwise_2d.jl} | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)
 rename src/callbacks_step/{analysis_surface_2d.jl => analysis_surface_pointwise_2d.jl} (99%)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 046f6327da5..54702295c17 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -727,7 +727,7 @@ end # @muladd
 include("analysis_dg1d.jl")
 include("analysis_dg2d.jl")
 include("analysis_surface_integral_2d.jl")
-include("analysis_surface_2d.jl")
+include("analysis_surface_pointwise_2d.jl")
 include("analysis_dg2d_parallel.jl")
 include("analysis_dg3d.jl")
 include("analysis_dg3d_parallel.jl")
diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
similarity index 99%
rename from src/callbacks_step/analysis_surface_2d.jl
rename to src/callbacks_step/analysis_surface_pointwise_2d.jl
index 5e9503051d6..05edc658bc5 100644
--- a/src/callbacks_step/analysis_surface_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -167,6 +167,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
             coordinates[global_node_index, 1] = x[1]
             coordinates[global_node_index, 2] = x[2]
             values[global_node_index] = value
+
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1
@@ -202,7 +203,6 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
 
     gradients_x, gradients_y = gradients
 
-    surface_integral = zero(eltype(u))
     index_range = eachnode(dg)
     global_node_index = 1 # Keeps track of solution point number on the surface
     for boundary in indices
@@ -239,6 +239,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
             coordinates[global_node_index, 1] = x[1]
             coordinates[global_node_index, 2] = x[2]
             values[global_node_index] = value
+
             i_node += i_node_step
             j_node += j_node_step
             global_node_index += 1

From 1b52899548d6e37e64cb963e857b631e9cb40b5a Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 11:35:47 +0200
Subject: [PATCH 44/59] viscous spelled out

---
 .../analysis_surface_integral_2d.jl           | 20 ++++++++++---------
 .../analysis_surface_pointwise_2d.jl          |  7 ++++---
 2 files changed, 15 insertions(+), 12 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index db1ad0182c1..335101c8705 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -231,29 +231,31 @@ function viscous_stress_vector(u, normal_direction, equations_parabolic,
                                                    gradients_1, gradients_2)
 
     # Viscous stress vector: Stress tensor * normal vector
-    visc_stress_vector_1 = tau_11 * n_normal[1] + tau_12 * n_normal[2]
-    visc_stress_vector_2 = tau_12 * n_normal[1] + tau_22 * n_normal[2]
+    viscous_stress_vector_1 = tau_11 * n_normal[1] + tau_12 * n_normal[2]
+    viscous_stress_vector_2 = tau_12 * n_normal[1] + tau_22 * n_normal[2]
 
-    return (visc_stress_vector_1, visc_stress_vector_2)
+    return (viscous_stress_vector_1, viscous_stress_vector_2)
 end
 
 function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
-                                               gradients_1, gradients_2)
+    viscous_stress_vector = viscous_stress_vector(u, normal_direction,
+                                                  equations_parabolic,
+                                                  gradients_1, gradients_2)
     @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
-    return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
+    return (viscous_stress_vector[1] * psi[1] + viscous_stress_vector[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
 function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
-                                               gradients_1, gradients_2)
+    viscous_stress_vector = viscous_stress_vector(u, normal_direction,
+                                                  equations_parabolic,
+                                                  gradients_1, gradients_2)
     @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
-    return (visc_stress_vector[1] * psi[1] + visc_stress_vector[2] * psi[2]) /
+    return (viscous_stress_vector[1] * psi[1] + viscous_stress_vector[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
index 05edc658bc5..f870d9e4010 100644
--- a/src/callbacks_step/analysis_surface_pointwise_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -117,14 +117,15 @@ end
 function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    visc_stress_vector = viscous_stress_vector(u, normal_direction, equations_parabolic,
-                                               gradients_1, gradients_2)
+    viscous_stress_vector = viscous_stress_vector(u, normal_direction,
+                                                  equations_parabolic,
+                                                  gradients_1, gradients_2)
     @unpack rhoinf, uinf, linf = surface_friction.flow_state
 
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = normal_direction / norm(normal_direction)
     n_perp = (-n[2], n[1])
-    return (visc_stress_vector[1] * n_perp[1] + visc_stress_vector[2] * n_perp[2]) /
+    return (viscous_stress_vector[1] * n_perp[1] + viscous_stress_vector[2] * n_perp[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 

From c37b20311a39059135cf4de2f13b224c6a189ce1 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 11:36:09 +0200
Subject: [PATCH 45/59] fmt

---
 .../elixir_navierstokes_NACA0012airfoil_mach08.jl   | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index c76e2a24d0d..221ff3f9b3a 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -143,14 +143,15 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   l_inf()))
 
 friction_coefficient = AnalysisSurfacePointwise(semi, force_boundary_names,
-                                       SurfaceFrictionCoefficient(rho_inf(),
-                                                                  u_inf(equations),
-                                                                  l_inf()))
+                                                SurfaceFrictionCoefficient(rho_inf(),
+                                                                           u_inf(equations),
+                                                                           l_inf()))
 
 pressure_coefficient = AnalysisSurfacePointwise(semi, force_boundary_names,
-                                       SurfacePressureCoefficient(p_inf(), rho_inf(),
-                                                                  u_inf(equations),
-                                                                  l_inf()))
+                                                SurfacePressureCoefficient(p_inf(),
+                                                                           rho_inf(),
+                                                                           u_inf(equations),
+                                                                           l_inf()))
 
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",

From a7b5f8f7e66235c1a2f46aee5af35ff23e70e57f Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 12:00:39 +0200
Subject: [PATCH 46/59] distinguish func and var

---
 .../analysis_surface_integral_2d.jl              | 16 ++++++++--------
 .../analysis_surface_pointwise_2d.jl             |  9 +++++----
 2 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index 335101c8705..6c4c17304b4 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -240,22 +240,22 @@ end
 function (lift_coefficient::LiftCoefficientShearStress)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    viscous_stress_vector = viscous_stress_vector(u, normal_direction,
-                                                  equations_parabolic,
-                                                  gradients_1, gradients_2)
+    viscous_stress_vector_ = viscous_stress_vector(u, normal_direction,
+                                                   equations_parabolic,
+                                                   gradients_1, gradients_2)
     @unpack psi, rhoinf, uinf, linf = lift_coefficient.flow_state
-    return (viscous_stress_vector[1] * psi[1] + viscous_stress_vector[2] * psi[2]) /
+    return (viscous_stress_vector_[1] * psi[1] + viscous_stress_vector_[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
 function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    viscous_stress_vector = viscous_stress_vector(u, normal_direction,
-                                                  equations_parabolic,
-                                                  gradients_1, gradients_2)
+    viscous_stress_vector_ = viscous_stress_vector(u, normal_direction,
+                                                   equations_parabolic,
+                                                   gradients_1, gradients_2)
     @unpack psi, rhoinf, uinf, linf = drag_coefficient.flow_state
-    return (viscous_stress_vector[1] * psi[1] + viscous_stress_vector[2] * psi[2]) /
+    return (viscous_stress_vector_[1] * psi[1] + viscous_stress_vector_[2] * psi[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
index f870d9e4010..a367ff9c4a2 100644
--- a/src/callbacks_step/analysis_surface_pointwise_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -117,15 +117,16 @@ end
 function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
-    viscous_stress_vector = viscous_stress_vector(u, normal_direction,
-                                                  equations_parabolic,
-                                                  gradients_1, gradients_2)
+    viscous_stress_vector_ = viscous_stress_vector(u, normal_direction,
+                                                   equations_parabolic,
+                                                   gradients_1, gradients_2)
     @unpack rhoinf, uinf, linf = surface_friction.flow_state
 
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = normal_direction / norm(normal_direction)
     n_perp = (-n[2], n[1])
-    return (viscous_stress_vector[1] * n_perp[1] + viscous_stress_vector[2] * n_perp[2]) /
+    return (viscous_stress_vector_[1] * n_perp[1] +
+            viscous_stress_vector_[2] * n_perp[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 

From 842eaebe896579e6ad05d5f1357aa123be7813db Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 12:05:38 +0200
Subject: [PATCH 47/59] fix coords

---
 src/callbacks_step/analysis_surface_pointwise_2d.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
index a367ff9c4a2..e0f3b5c9708 100644
--- a/src/callbacks_step/analysis_surface_pointwise_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -162,7 +162,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coordinates(node_coordinates, equations, dg, i_node, j_node,
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
                                      element)
             value = variable(u_node, equations)
 
@@ -221,7 +221,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
             u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index,
                                          boundary)
 
-            x = get_node_coordinates(node_coordinates, equations, dg, i_node, j_node,
+            x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
                                      element)
             # Extract normal direction at nodes which points from the elements outwards,
             # i.e., *into* the structure.

From 0374850c2a5d066f36ad6650ca72cdcf370d133d Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 24 May 2024 12:05:49 +0200
Subject: [PATCH 48/59] fmt

---
 src/callbacks_step/analysis_surface_pointwise_2d.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
index e0f3b5c9708..b82955dc854 100644
--- a/src/callbacks_step/analysis_surface_pointwise_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -163,7 +163,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
                                          boundary)
 
             x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
-                                     element)
+                                element)
             value = variable(u_node, equations)
 
             coordinates[global_node_index, 1] = x[1]
@@ -222,7 +222,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
                                          boundary)
 
             x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
-                                     element)
+                                element)
             # Extract normal direction at nodes which points from the elements outwards,
             # i.e., *into* the structure.
             normal_direction = get_normal_direction(direction, contravariant_vectors,

From f04c1ee8885956bff2e4ef2cf8aacf3ed4799f97 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 5 Jul 2024 14:20:59 +0200
Subject: [PATCH 49/59] update

---
 ...xir_navierstokes_NACA0012airfoil_mach08.jl |  4 +-
 src/callbacks_step/analysis.jl                | 20 +++++---
 .../analysis_surface_integral_2d.jl           |  6 +--
 .../analysis_surface_pointwise_2d.jl          | 48 ++++++++-----------
 4 files changed, 37 insertions(+), 41 deletions(-)

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index bc9f21fc903..e000c8f5566 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -142,12 +142,12 @@ lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_name
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-friction_coefficient = AnalysisSurfacePointwise(semi, force_boundary_names,
+friction_coefficient = AnalysisSurfacePointwise(force_boundary_names,
                                                 SurfaceFrictionCoefficient(rho_inf(),
                                                                            u_inf(equations),
                                                                            l_inf()))
 
-pressure_coefficient = AnalysisSurfacePointwise(semi, force_boundary_names,
+pressure_coefficient = AnalysisSurfacePointwise(force_boundary_names,
                                                 SurfacePressureCoefficient(p_inf(),
                                                                            rho_inf(),
                                                                            u_inf(equations),
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 1a29f198620..dfb7efde893 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -698,8 +698,8 @@ include("analysis_dg3d_parallel.jl")
 # This version of `analyze` is used for [`AnalysisSurfaceIntegral`](@ref) which requires
 # `semi` to be passed along to retrieve the current boundary indices, which are non-static 
 # in the case of AMR.
-function analyze(quantity::AnalysisSurfaceIntegral, du, u, t,
-                 semi::AbstractSemidiscretization)
+function analyze(quantity::AnalysisSurfaceIntegral{Variable}, du, u, t,
+                 semi::AbstractSemidiscretization) where {Variable}
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache, semi)
 end
@@ -708,8 +708,7 @@ end
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
 # Note that this needs to be included after `analysis_surface_integral_2d.jl` to
 # have `VariableViscous` available.
-function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
-                                 AnalysisSurfacePointwise{Variable}},
+function analyze(quantity::AnalysisSurfaceIntegral{Variable},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
                                                                      Variable <:
@@ -720,8 +719,15 @@ function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
     analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, semi,
             cache_parabolic)
 end
-function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
-                                 AnalysisSurfacePointwise{Variable}},
+
+function analyze(quantity::AnalysisSurfacePointwise{Variable},
+                 du, u, t,
+                 semi::AbstractSemidiscretization,
+                 iter) where {Variable}
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    analyze(quantity, du, u, t, mesh, equations, solver, cache, semi, iter)
+end
+function analyze(quantity::AnalysisSurfacePointwise{Variable},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic,
                  iter) where {
@@ -730,6 +736,6 @@ function analyze(quantity::Union{AnalysisSurfaceIntegral{Variable},
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     equations_parabolic = semi.equations_parabolic
     cache_parabolic = semi.cache_parabolic
-    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
+    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, semi,
             cache_parabolic, iter)
 end
diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index 95330e885c8..383739c08e6 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -9,9 +9,9 @@
 # surface forces
 
 """
-    AnalysisSurfaceIntegral{Semidiscretization, Variable}(semi,
-                                                          boundary_symbol_or_boundary_symbols,
-                                                          variable)
+    AnalysisSurfaceIntegral{Variable, NBoundaries}(semi,
+                                                   boundary_symbol_or_boundary_symbols,
+                                                   variable)
 
 This struct is used to compute the surface integral of a quantity of interest `variable` alongside
 the boundary/boundaries associated with particular name(s) given in `boundary_symbol`
diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
index b82955dc854..7e7de2ef0e9 100644
--- a/src/callbacks_step/analysis_surface_pointwise_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -9,9 +9,8 @@
 # pointwise surface forces.
 
 """
-    AnalysisSurfacePointwise{Semidiscretization, Variable}(semi,
-                                                  boundary_symbol_or_boundary_symbols,
-                                                  variable)
+    AnalysisSurfacePointwise{Variable, NBoundaries}(boundary_symbol_or_boundary_symbols,
+                                                    variable, output_directory = "out")
 
 This struct is used to compute pointwise surface values of a quantity of interest `variable`
 alongside the boundary/boundaries associated with particular name(s) given in `boundary_symbol`
@@ -20,34 +19,21 @@ For instance, this can be used to compute the surface pressure coefficient [`Sur
 surface friction coefficient [`SurfaceFrictionCoefficient`](@ref) of e.g. an airfoil with the boundary
 symbol `:Airfoil` in 2D.
 
-- `semi::Semidiscretization`: Passed in to retrieve boundary condition information
-- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries stored as symbol(s)
+- `boundary_symbols::NTuple{NBoundaries, Symbol}`: Name(s) of the boundary/boundaries
   where the quantity of interest is computed
 - `variable::Variable`: Quantity of interest, like lift or drag
+- `output_directory = "out"`: Directory where the pointwise value files are stored.
 """
-struct AnalysisSurfacePointwise{Variable}
-    indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed
+struct AnalysisSurfacePointwise{Variable, NBoundaries}
     variable::Variable # Quantity of interest, like lift or drag
+    boundary_symbols::NTuple{NBoundaries, Symbol} # Name(s) of the boundary/boundaries
     output_directory::String
 
-    function AnalysisSurfacePointwise(semi, boundary_symbol, variable,
-                                      output_directory = "out")
-        @unpack boundary_symbol_indices = semi.boundary_conditions
-        indices = boundary_symbol_indices[boundary_symbol]
-
-        return new{typeof(variable)}(indices, variable, output_directory)
-    end
-
-    function AnalysisSurfacePointwise(semi, boundary_symbols::Vector{Symbol}, variable,
-                                      output_directory = "out")
-        @unpack boundary_symbol_indices = semi.boundary_conditions
-        indices = Vector{Int}()
-        for name in boundary_symbols
-            append!(indices, boundary_symbol_indices[name])
-        end
-        sort!(indices)
-
-        return new{typeof(variable)}(indices, variable, output_directory)
+    function AnalysisSurfacePointwise(boundary_symbols::NTuple{NBoundaries, Symbol},
+                                      variable,
+                                      output_directory = "out") where {NBoundaries}
+        return new{typeof(variable), NBoundaries}(variable, boundary_symbols,
+                                                  output_directory)
     end
 end
 
@@ -132,12 +118,14 @@ end
 
 function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
                  mesh::P4estMesh{2},
-                 equations, dg::DGSEM, cache, iter)
+                 equations, dg::DGSEM, cache, semi, iter)
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
 
-    @unpack indices, variable = surface_variable
+    @unpack variable, boundary_symbols = surface_variable
+    @unpack boundary_symbol_indices = semi.boundary_conditions
+    indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices)
 
     dim = ndims(mesh)
     n_nodes = nnodes(dg)
@@ -184,13 +172,15 @@ end
 function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
                  du, u, t, mesh::P4estMesh{2},
                  equations, equations_parabolic,
-                 dg::DGSEM, cache,
+                 dg::DGSEM, cache, semi,
                  cache_parabolic, iter) where {Variable <: VariableViscous}
     @unpack boundaries = cache
     @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
     @unpack weights = dg.basis
 
-    @unpack indices, variable = surface_variable
+    @unpack variable, boundary_symbols = surface_variable
+    @unpack boundary_symbol_indices = semi.boundary_conditions
+    indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices)
 
     dim = ndims(mesh)
     n_nodes = nnodes(dg)

From 42afa387e6196c7576b57beca5e6c20bc12d486c Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 5 Jul 2024 14:24:20 +0200
Subject: [PATCH 50/59] correct docstring

---
 src/callbacks_step/analysis_surface_integral_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl
index 383739c08e6..1a95db66d81 100644
--- a/src/callbacks_step/analysis_surface_integral_2d.jl
+++ b/src/callbacks_step/analysis_surface_integral_2d.jl
@@ -10,7 +10,7 @@
 
 """
     AnalysisSurfaceIntegral{Variable, NBoundaries}(semi,
-                                                   boundary_symbol_or_boundary_symbols,
+                                                   boundary_symbols::NTuple{NBoundaries, Symbol},
                                                    variable)
 
 This struct is used to compute the surface integral of a quantity of interest `variable` alongside

From b9d2fcba45926c6c3a0a8099b29b42b31cae1634 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 5 Jul 2024 17:01:23 +0200
Subject: [PATCH 51/59] sort arrays

---
 test/test_parabolic_2d.jl | 44 ++++++++++++++++++---------------------
 1 file changed, 20 insertions(+), 24 deletions(-)

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 96e216517f8..38bacf5d41d 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -760,30 +760,26 @@ end
 
     cp_vals = read(Trixi.h5open("out/CP_x_000007.h5"))["point_data"]
     cf_vals = read(Trixi.h5open("out/CF_x_000007.h5"))["point_data"]
-    @test cp_vals[1:10] ≈ [
-        2.0956920203526193,
-        2.098639791879171,
-        2.0354491762572806,
-        2.1591959225932777,
-        1.7417814818119175,
-        2.0400855926847936,
-        1.874247697759559,
-        1.871513931283643,
-        1.5182860925755697,
-        1.5152278762705806,
-    ]
-    @test cf_vals[1:10] ≈ [
-        0.2038354965561312,
-        0.2634810993266428,
-        0.42919006623618894,
-        0.045526882648679115,
-        0.6611394672752093,
-        -0.23515678100990495,
-        -0.37053727238603884,
-        -0.4537078569805583,
-        0.8113349266282898,
-        0.7458959144321556,
-    ]
+    @test sort(cp_vals[1:10]) ≈ [1.5152278762705806
+           1.5182860925755697
+           1.7417814818119175
+           1.871513931283643
+           1.874247697759559
+           2.0354491762572806
+           2.0400855926847936
+           2.0956920203526193
+           2.098639791879171
+           2.1591959225932777]
+    @test sort!(cf_vals[1:10]) ≈ [-0.4537078569805583
+           -0.37053727238603884
+           -0.23515678100990495
+           0.045526882648679115
+           0.2038354965561312
+           0.2634810993266428
+           0.42919006623618894
+           0.6611394672752093
+           0.7458959144321556
+           0.8113349266282898]
 end
 end
 

From a289eb1b4ac2e6c259fa51ece0b5a9da1b805739 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 5 Jul 2024 17:01:37 +0200
Subject: [PATCH 52/59] mno in place

---
 test/test_parabolic_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 38bacf5d41d..ed84322b7c0 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -770,7 +770,7 @@ end
            2.0956920203526193
            2.098639791879171
            2.1591959225932777]
-    @test sort!(cf_vals[1:10]) ≈ [-0.4537078569805583
+    @test sort(cf_vals[1:10]) ≈ [-0.4537078569805583
            -0.37053727238603884
            -0.23515678100990495
            0.045526882648679115

From 634f60cf37c8beb76a3cbfc920bc635b82bd59b0 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 5 Sep 2024 12:23:26 +0200
Subject: [PATCH 53/59] Use counter instead of index to avoid confusion

---
 .../analysis_surface_pointwise_2d.jl          | 20 +++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
index 7e7de2ef0e9..4fb05db0d6d 100644
--- a/src/callbacks_step/analysis_surface_pointwise_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -136,7 +136,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
 
     index_range = eachnode(dg)
 
-    global_node_index = 1 # Keeps track of solution point number on the surface
+    global_node_counter = 1 # Keeps track of solution point number on the surface
     for boundary in indices
         element = boundaries.neighbor_ids[boundary]
         node_indices = boundaries.node_indices[boundary]
@@ -154,13 +154,13 @@ function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
                                 element)
             value = variable(u_node, equations)
 
-            coordinates[global_node_index, 1] = x[1]
-            coordinates[global_node_index, 2] = x[2]
-            values[global_node_index] = value
+            coordinates[global_node_counter, 1] = x[1]
+            coordinates[global_node_counter, 2] = x[2]
+            values[global_node_counter] = value
 
             i_node += i_node_step
             j_node += j_node_step
-            global_node_index += 1
+            global_node_counter += 1
         end
     end
 
@@ -196,7 +196,7 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
     gradients_x, gradients_y = gradients
 
     index_range = eachnode(dg)
-    global_node_index = 1 # Keeps track of solution point number on the surface
+    global_node_counter = 1 # Keeps track of solution point number on the surface
     for boundary in indices
         element = boundaries.neighbor_ids[boundary]
         node_indices = boundaries.node_indices[boundary]
@@ -228,13 +228,13 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
             value = variable(u_node, normal_direction, x, t, equations_parabolic,
                              gradients_1, gradients_2)
 
-            coordinates[global_node_index, 1] = x[1]
-            coordinates[global_node_index, 2] = x[2]
-            values[global_node_index] = value
+            coordinates[global_node_counter, 1] = x[1]
+            coordinates[global_node_counter, 2] = x[2]
+            values[global_node_counter] = value
 
             i_node += i_node_step
             j_node += j_node_step
-            global_node_index += 1
+            global_node_counter += 1
         end
     end
 

From d3a57edf950f488a4a5f97b3a45bf9f4b10960bc Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 5 Sep 2024 14:36:40 +0200
Subject: [PATCH 54/59] Comments

---
 .../analysis_surface_pointwise_2d.jl          | 55 +++++++++++++------
 1 file changed, 39 insertions(+), 16 deletions(-)

diff --git a/src/callbacks_step/analysis_surface_pointwise_2d.jl b/src/callbacks_step/analysis_surface_pointwise_2d.jl
index 4fb05db0d6d..145bba56374 100644
--- a/src/callbacks_step/analysis_surface_pointwise_2d.jl
+++ b/src/callbacks_step/analysis_surface_pointwise_2d.jl
@@ -93,6 +93,9 @@ function SurfaceFrictionCoefficient(rhoinf, uinf, linf)
     return SurfaceFrictionCoefficient(FlowState(rhoinf, uinf, linf))
 end
 
+# Compute local pressure coefficient.
+# Works for both purely hyperbolic and hyperbolic-parabolic systems.
+# C_p(x) = (p(x) - p_inf) / (0.5 * rho_inf * u_inf^2 * l_inf)
 function (pressure_coefficient::SurfacePressureCoefficient)(u, equations)
     p = pressure(u, equations)
     @unpack pinf = pressure_coefficient
@@ -100,6 +103,9 @@ function (pressure_coefficient::SurfacePressureCoefficient)(u, equations)
     return (p - pinf) / (0.5 * rhoinf * uinf^2 * linf)
 end
 
+# Compute local friction coefficient.
+# Works only in conjunction with a hyperbolic-parabolic system.
+# C_f(x) = (tau_w(x) * n_perp(x)) / (0.5 * rho_inf * u_inf^2 * l_inf)
 function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x, t,
                                                         equations_parabolic,
                                                         gradients_1, gradients_2)
@@ -110,12 +116,18 @@ function (surface_friction::SurfaceFrictionCoefficient)(u, normal_direction, x,
 
     # Normalize as `normal_direction` is not necessarily a unit vector
     n = normal_direction / norm(normal_direction)
-    n_perp = (-n[2], n[1])
-    return (viscous_stress_vector_[1] * n_perp[1] +
-            viscous_stress_vector_[2] * n_perp[2]) /
+    # Tangent vector = perpendicular vector to normal vector
+    t = (-n[2], n[1])
+    return (viscous_stress_vector_[1] * t[1] +
+            viscous_stress_vector_[2] * t[2]) /
            (0.5 * rhoinf * uinf^2 * linf)
 end
 
+# Compute and save to disk a space-dependent `surface_variable`.
+# For the purely hyperbolic, i.e., non-parabolic case, this is for instance 
+# the pressure coefficient `SurfacePressureCoefficient`.
+# The boundary/boundaries along which this quantity is to be integrated is determined by
+# `boundary_symbols`, which is retrieved from `surface_variable`.
 function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
                  mesh::P4estMesh{2},
                  equations, dg::DGSEM, cache, semi, iter)
@@ -164,11 +176,16 @@ function analyze(surface_variable::AnalysisSurfacePointwise, du, u, t,
         end
     end
 
+    # Save to disk
     save_pointwise_file(surface_variable.output_directory, varname(variable),
-                        coordinates,
-                        values, t, iter)
+                        coordinates, values, t, iter)
 end
 
+# Compute and save to disk a space-dependent `surface_variable`.
+# For the purely hyperbolic-parabolic case, this is may be for instance 
+# the surface skin fricition coefficient `SurfaceFrictionCoefficient`.
+# The boundary/boundaries along which this quantity is to be integrated is determined by
+# `boundary_symbols`, which is retrieved from `surface_variable`.
 function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
                  du, u, t, mesh::P4estMesh{2},
                  equations, equations_parabolic,
@@ -213,16 +230,16 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
 
             x = get_node_coords(node_coordinates, equations, dg, i_node, j_node,
                                 element)
-            # Extract normal direction at nodes which points from the elements outwards,
-            # i.e., *into* the structure.
-            normal_direction = get_normal_direction(direction, contravariant_vectors,
-                                                    i_node, j_node,
-                                                    element)
 
-            gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg, i_node,
-                                        j_node, element)
-            gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg, i_node,
-                                        j_node, element)
+            gradients_1 = get_node_vars(gradients_x, equations_parabolic, dg,
+                                        i_node, j_node, element)
+            gradients_2 = get_node_vars(gradients_y, equations_parabolic, dg,
+                                        i_node, j_node, element)
+
+            # Extract normal direction at nodes which points from the 
+            # fluid cells *outwards*, i.e., *into* the structure.
+            normal_direction = get_normal_direction(direction, contravariant_vectors,
+                                                    i_node, j_node, element)
 
             # Integral over whole boundary surface
             value = variable(u_node, normal_direction, x, t, equations_parabolic,
@@ -238,15 +255,21 @@ function analyze(surface_variable::AnalysisSurfacePointwise{Variable},
         end
     end
 
+    # Save to disk
     save_pointwise_file(surface_variable.output_directory, varname(variable),
-                        coordinates,
-                        values, t, iter)
+                        coordinates, values, t, iter)
 end
 
 varname(::Any) = @assert false "Surface variable name not assigned" # This makes sure default behaviour is not overwriting
 varname(pressure_coefficient::SurfacePressureCoefficient) = "CP_x"
 varname(friction_coefficient::SurfaceFrictionCoefficient) = "CF_x"
 
+# Helper function that saves a space-dependent quantity `values`
+# at every solution/quadrature point `coordinates` at 
+# time `t` and iteration `iter` to disk.
+# The file is written to the `output_directory` with name `varname` in HDF5 (.h5) format.
+# The latter two are retrieved from the `surface_variable`,
+# an instantiation of `AnalysisSurfacePointwise`.
 function save_pointwise_file(output_directory, varname, coordinates, values, t, iter)
     n_points = length(values)
 

From 84f1352f07fd78271c6f9a26c4e2e397e66cfb4c Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 5 Sep 2024 14:58:24 +0200
Subject: [PATCH 55/59] Update docstring

---
 src/callbacks_step/analysis.jl | 24 +++++++++++++++++++-----
 1 file changed, 19 insertions(+), 5 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index dfb7efde893..707a284ab5c 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -9,11 +9,15 @@
 # - analysis_interval part as PeriodicCallback called after a certain amount of simulation time
 """
     AnalysisCallback(semi; interval=0,
-                            save_analysis=false,
-                            output_directory="out",
-                            analysis_filename="analysis.dat",
-                            extra_analysis_errors=Symbol[],
-                            extra_analysis_integrals=())
+                           save_analysis=false,
+                           output_directory="out",
+                           analysis_filename="analysis.dat",
+                           analysis_errors = union(default_analysis_errors(equations),
+                                                   extra_analysis_errors),
+                           extra_analysis_integrals = (),
+                           analysis_integrals = union(default_analysis_integrals(equations),
+                                                      extra_analysis_integrals),
+                           analysis_pointwise = ())
 
 Analyze a numerical solution every `interval` time steps and print the
 results to the screen. If `save_analysis`, the results are also saved in
@@ -35,6 +39,16 @@ solution and integrated over the computational domain. Some examples for this ar
 [`entropy`](@ref), [`energy_kinetic`](@ref), [`energy_internal`](@ref), and [`energy_total`](@ref).
 You can also write your own function with the same signature as the examples listed above and
 pass it via `extra_analysis_integrals`.
+The default `analysis_integrals` is `(entropy_timederivative,)`.
+You can also request `extra_analysis_integrals` such as [`LiftCoefficientPressure`](@ref) or
+[`DragCoefficientPressure`](@ref) by constructing an [`AnalysisSurfaceIntegral`](@ref) with one of 
+the previously mentioned functions.
+
+Similarly, pointwise, i.e., per quadrature/interpolation point, quantities such at
+[`SurfacePressureCoefficient`](@ref) or [`SurfaceFrictionCoefficient`](@ref) can be computed.
+Instances of these need to be passed into [`AnalysisSurfacePointwise`](@ref) which is then in turn
+passed to `analysis_pointwise`.
+
 See the developer comments about `Trixi.analyze`, `Trixi.pretty_form_utf`, and
 `Trixi.pretty_form_ascii` for further information on how to create custom analysis quantities.
 

From 3c952f77ac6edd6aba4a640ed4b28dd0c356aaa9 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 5 Sep 2024 15:12:55 +0200
Subject: [PATCH 56/59] More comments

---
 src/callbacks_step/analysis.jl | 20 +++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 707a284ab5c..e6649e1d471 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -637,6 +637,9 @@ function analyze(quantity, du, u, t, semi::AbstractSemidiscretization)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache)
 end
+# `analyze` function that passes also the iteration number `iter`along.
+# Required for callbacks that handle the write-off of results to disk themselves,
+# such as `AnalysisSurfacePointwise`.
 function analyze(quantity, du, u, t, semi::AbstractSemidiscretization, iter)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache, iter)
@@ -712,7 +715,8 @@ include("analysis_dg3d_parallel.jl")
 # This version of `analyze` is used for [`AnalysisSurfaceIntegral`](@ref) which requires
 # `semi` to be passed along to retrieve the current boundary indices, which are non-static 
 # in the case of AMR.
-function analyze(quantity::AnalysisSurfaceIntegral{Variable}, du, u, t,
+function analyze(quantity::AnalysisSurfaceIntegral{Variable}, 
+                 du, u, t,
                  semi::AbstractSemidiscretization) where {Variable}
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache, semi)
@@ -720,8 +724,10 @@ end
 
 # Special analyze for `SemidiscretizationHyperbolicParabolic` such that
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
-# Note that this needs to be included after `analysis_surface_integral_2d.jl` to
-# have `VariableViscous` available.
+# Note that this needs to be defined after `analysis_surface_integral_2d.jl` is included 
+# to have `VariableViscous` available.
+# We require `semi` to be passed along to retrieve the current boundary indices, which are non-static 
+# in the case of AMR.
 function analyze(quantity::AnalysisSurfaceIntegral{Variable},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic) where {
@@ -734,6 +740,9 @@ function analyze(quantity::AnalysisSurfaceIntegral{Variable},
             cache_parabolic)
 end
 
+# This version of `analyze` is used for `AnalysisSurfacePointwise` such as `SurfacePressureCoefficient`.
+# We need the iteration number `iter` to be passed in here 
+# as for `AnalysisSurfacePointwise` the writing to disk is handled by the callback itself.
 function analyze(quantity::AnalysisSurfacePointwise{Variable},
                  du, u, t,
                  semi::AbstractSemidiscretization,
@@ -741,6 +750,11 @@ function analyze(quantity::AnalysisSurfacePointwise{Variable},
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
     analyze(quantity, du, u, t, mesh, equations, solver, cache, semi, iter)
 end
+# Special analyze for `SemidiscretizationHyperbolicParabolic` such that
+# precomputed gradients are available. Required for `AnalysisSurfacePointwise` equipped 
+# with `VariableViscous` such as `SurfaceFrictionCoefficient`.
+# As for the inviscid version, we need to pass in the iteration number `iter` as 
+# for `AnalysisSurfacePointwise` the writing to disk is handled by the callback itself.
 function analyze(quantity::AnalysisSurfacePointwise{Variable},
                  du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic,

From ed971acba5b44e2096d5ebabd017764ad6405410 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 5 Sep 2024 15:15:13 +0200
Subject: [PATCH 57/59] fmt

---
 src/callbacks_step/analysis.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index e6649e1d471..a8c83d88b44 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -715,7 +715,7 @@ include("analysis_dg3d_parallel.jl")
 # This version of `analyze` is used for [`AnalysisSurfaceIntegral`](@ref) which requires
 # `semi` to be passed along to retrieve the current boundary indices, which are non-static 
 # in the case of AMR.
-function analyze(quantity::AnalysisSurfaceIntegral{Variable}, 
+function analyze(quantity::AnalysisSurfaceIntegral{Variable},
                  du, u, t,
                  semi::AbstractSemidiscretization) where {Variable}
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)

From 2f79a036f71f6c45d402b7b011402aaed1a4914d Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 5 Sep 2024 16:17:30 +0200
Subject: [PATCH 58/59] update test vals

---
 test/test_parabolic_2d.jl | 44 +++++++++++++++++++++------------------
 1 file changed, 24 insertions(+), 20 deletions(-)

diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index ed84322b7c0..39511a6d4eb 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -760,26 +760,30 @@ end
 
     cp_vals = read(Trixi.h5open("out/CP_x_000007.h5"))["point_data"]
     cf_vals = read(Trixi.h5open("out/CF_x_000007.h5"))["point_data"]
-    @test sort(cp_vals[1:10]) ≈ [1.5152278762705806
-           1.5182860925755697
-           1.7417814818119175
-           1.871513931283643
-           1.874247697759559
-           2.0354491762572806
-           2.0400855926847936
-           2.0956920203526193
-           2.098639791879171
-           2.1591959225932777]
-    @test sort(cf_vals[1:10]) ≈ [-0.4537078569805583
-           -0.37053727238603884
-           -0.23515678100990495
-           0.045526882648679115
-           0.2038354965561312
-           0.2634810993266428
-           0.42919006623618894
-           0.6611394672752093
-           0.7458959144321556
-           0.8113349266282898]
+    @test sort(cp_vals[1:10]) ≈ [
+        1.7642563337066723,
+        1.8309926967660561,
+        1.8676797245949768,
+        1.8715139312836404,
+        1.8742476977595561,
+        1.9596647512001575,
+        1.9749658702816537,
+        2.0400855926847923,
+        2.095692020352621,
+        2.1591959225932738,
+    ]
+    @test sort(cf_vals[1:10]) ≈ [
+        -1.2838932978862392,
+        -1.1230998392551437,
+        -1.1167735453217076,
+        -0.9947974567347279,
+        -0.6878311452704875,
+        -0.45370785698055793,
+        -0.37053727238603945,
+        -0.23515678100990403,
+        0.045526882648679205,
+        0.20383549655613004,
+    ]
 end
 end
 

From 9d7231d583c7d2962d5cbece75f5321c142b679c Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 5 Sep 2024 16:24:45 +0200
Subject: [PATCH 59/59] do not pollute `analysis.dat`

---
 src/callbacks_step/analysis.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index a8c83d88b44..f1b9b6c3d28 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -221,10 +221,8 @@ function initialize!(cb::DiscreteCallback{Condition, Affect!}, u_ode, du_ode, t,
             for quantity in analysis_integrals
                 @printf(io, "   %-14s", pretty_form_ascii(quantity))
             end
-
-            for quantity in analysis_pointwise
-                @printf(io, "   %-14s", pretty_form_ascii(quantity))
-            end
+            # Pointwise quantities are not saved in `analysis_filename`,
+            # i.e., `analysis.dat` but handle their own output.
 
             println(io)
         end