Apply suggestions from code review

Some key points

1) Change pre to p
2) Don't capitalize U
3) Don't append Swanson quantities with 'sw' 
4) linf -> l_inf (may be pending in some other places)
5) Other formatting improvements

Co-authored-by: Andrew Winters <[email protected]>

examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl

@@ -7,17 +7,17 @@ using Trixi
 
 equations = CompressibleEulerEquations2D(1.4)
 
-pre_inf() = 1.0
-rho_inf() = pre_inf() / (1.0 * 287.87) # pre_inf = 1.0,  T = 1, R = 287.87
+p_inf() = 1.0
+rho_inf() = p_inf() / (1.0 * 287.87) # p_inf = 1.0,  T = 1, R = 287.87
 mach_inf() = 0.85
 aoa() = pi / 180.0 # 1 Degree angle of attack
-c_inf(equations) = sqrt(equations.gamma * pre_inf() / rho_inf())
-U_inf(equations) = mach_inf() * c_inf(equations)
+c_inf(equations) = sqrt(equations.gamma * p_inf() / rho_inf())
+u_inf(equations) = mach_inf() * c_inf(equations)
 
 @inline function initial_condition_mach085_flow(x, t,
                                                 equations::CompressibleEulerEquations2D)
-    v1 = U_inf(equations) * cos(aoa())
-    v2 = U_inf(equations) * sin(aoa())
+    v1 = u_inf(equations) * cos(aoa())
+    v2 = u_inf(equations) * sin(aoa())
 
     prim = SVector(rho_inf(), v1, v2, pre_inf())
     return prim2cons(prim, equations)
@@ -82,16 +82,16 @@ summary_callback = SummaryCallback()
 
 analysis_interval = 2000
 
-linf = 1.0 # Length of airfoil
+l_inf = 1.0 # Length of airfoil
 
 force_boundary_names = [:AirfoilBottom, :AirfoilTop]
 drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
-                                                                   U_inf(equations), linf))
+                                                                   u_inf(equations), linf))
 
 lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
                                            LiftCoefficientPressure(aoa(), rho_inf(),
-                                                                   U_inf(equations), linf))
+                                                                   u_inf(equations), linf))
 
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",

examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl

@@ -47,12 +47,12 @@ end
 cells_per_dimension = (64, 64)
 # xi = -1 maps to the inner circle and xi = +1 maps to the outer circle and we can specify boundary
 # conditions there. However, the image of eta = -1, +1 coincides at the line y = 0. There is no
-# physical boundary there so we specify periodicity = true there and the solver treats the
+# physical boundary there so we specify `periodicity = true` there and the solver treats the
 # element across eta = -1, +1 as neighbours which is what we want
 mesh = P4estMesh(cells_per_dimension, mapping = mapping2cylinder, polydeg = polydeg,
                  periodicity = (false, true))
 
-# The boundary conditions of the outer cylinder is constant but subsonic, so we cannot compute the
+# The boundary condition on the outer cylinder is constant but subsonic, so we cannot compute the
 # boundary flux from the external information alone. Thus, we use the numerical flux to distinguish
 # between inflow and outflow characteristics
 @inline function boundary_condition_subsonic_constant(u_inner,
@@ -86,16 +86,16 @@ analysis_interval = 2000
 
 aoa = 0.0
 rho_inf = 1.4
-U_inf = 0.38
-linf = 1.0 # Diameter of circle
+u_inf = 0.38
+l_inf = 1.0 # Diameter of circle
 
 drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg,
-                                           DragCoefficientPressure(aoa, rho_inf, U_inf,
-                                                                   linf))
+                                           DragCoefficientPressure(aoa, rho_inf, u_inf,
+                                                                   l_inf))
 
 lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg,
-                                           LiftCoefficientPressure(aoa, rho_inf, U_inf,
-                                                                   linf))
+                                           LiftCoefficientPressure(aoa, rho_inf, u_inf,
+                                                                   l_inf))
 
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",

examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl

@@ -7,14 +7,13 @@ using Trixi
 
 # Laminar transonic flow around an airfoil.
 
-# This test is taken from the paper below. The values from Case 5 in Table 3 are used to validate
-# the scheme and computation of surface forces.
+# This test is taken from the paper of Swanson and Langer. The values for the drag and lift coefficients
+# from Case 5 in Table 3 are used to validate the scheme and computation of surface forces.
 
+# References:
 # - Roy Charles Swanson, Stefan Langer (2016)
 #   Structured and Unstructured Grid Methods (2016)
 #   [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623)
-
-# See also:
 # - Deep Ray, Praveen Chandrashekar (2017)
 #   An entropy stable finite volume scheme for the 
 #   two dimensional Navier–Stokes equations on triangular grids
@@ -27,25 +26,25 @@ mu() = 0.0031959974968701088
 equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
                                                           Prandtl = prandtl_number())
 
-sw_rho_inf() = 1.0
-sw_pre_inf() = 2.85
-sw_aoa() = 10.0 * pi / 180.0
-sw_linf() = 1.0
-sw_mach_inf() = 0.8
-sw_U_inf(equations) = sw_mach_inf() * sqrt(equations.gamma * sw_pre_inf() / sw_rho_inf())
+rho_inf() = 1.0
+p_inf() = 2.85
+aoa() = 10.0 * pi / 180.0 # 10 degree angle of attack
+l_inf() = 1.0
+mach_inf() = 0.8
+u_inf(equations) = mach_inf() * sqrt(equations.gamma * p_inf() / rho_inf())
 @inline function initial_condition_mach08_flow(x, t, equations)
     # set the freestream flow parameters
-    gasGam = equations.gamma
-    mach_inf = sw_mach_inf()
-    aoa = sw_aoa()
-    rho_inf = sw_rho_inf()
-    pre_inf = sw_pre_inf()
-    U_inf = mach_inf * sqrt(gasGam * pre_inf / rho_inf)
+    gamma = equations.gamma
+    mach_inf = mach_inf()
+    aoa = aoa()
+    rho_inf = rho_inf()
+    p_inf = p_inf()
+    u_inf = mach_inf * sqrt(gamma * p_inf / rho_inf)
 
-    v1 = U_inf * cos(aoa)
-    v2 = U_inf * sin(aoa)
+    v1 = u_inf * cos(aoa)
+    v2 = u_inf * sin(aoa)
 
-    prim = SVector(rho_inf, v1, v2, pre_inf)
+    prim = SVector(rho_inf, v1, v2, p_inf)
     return prim2cons(prim, equations)
 end
 
@@ -126,14 +125,14 @@ analysis_interval = 2000
 
 force_boundary_names = [:AirfoilBottom, :AirfoilTop]
 drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
-                                           DragCoefficientPressure(sw_aoa(), sw_rho_inf(),
-                                                                   sw_U_inf(equations),
-                                                                   sw_linf()))
+                                           DragCoefficientPressure(aoa(), rho_inf(),
+                                                                   u_inf(equations),
+                                                                   l_inf()))
 
 lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
-                                           LiftCoefficientPressure(sw_aoa(), sw_rho_inf(),
-                                                                   sw_U_inf(equations),
-                                                                   sw_linf()))
+                                           LiftCoefficientPressure(aoa(), rho_inf(),
+                                                                   u_inf(equations),
+                                                                   l_inf()))
 
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",
@@ -154,8 +153,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav
 ###############################################################################
 # run the simulation
 
-time_int_tol = 1e-8
-sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = time_int_tol,
-            reltol = time_int_tol,
+sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1e-8, reltol = 1e-8,
             ode_default_options()..., callback = callbacks)
 summary_callback() # print the timer summary

src/callbacks_step/analysis_surface_integral_2d.jl

@@ -1,22 +1,24 @@
-# This file contains callbacks that are performed on the surface like computation of
-# surface forces
-
 # 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
+
 """
     AnalysisSurfaceIntegral{Semidiscretization, Variable}(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 `boundary_symbol` or `boundary_symbols`.
+    the boundary/boundaries associated with particular name(s) given in `boundary_symbol` 
+    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 in 2D.
+    drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the boundary
+    name `:Airfoil` in 2D.
 """
 struct AnalysisSurfaceIntegral{Semidiscretization, Variable}
     semi::Semidiscretization # passed in to retrieve boundary condition information
@@ -80,7 +82,7 @@ function LiftCoefficientPressure(aoa, rhoinf, uinf, linf)
     # psi_lift is the normal unit vector to the freestream direction.
     # Note: The choice of the normal vector psi_lift = (-sin(aoa), cos(aoa))
     # leads to positive lift coefficients for positive angles of attack for airfoils.
-    # Note that one could also use psi_lift = (sin(aoa), -cos(aoa)) which results in the same 
+    # 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))
@@ -157,7 +159,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t,
 
             # L2 norm of normal direction (contravariant_vector) is the surface element
             dS = weights[node_index] * norm(normal_direction)
-            # Integral over whole boundary surface
+            # Integral over entire boundary surface
             surface_integral += variable(u_node, normal_direction, equations) * dS
 
             i_node += i_node_step