Pointwise boundary surface values

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"

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,14 +141,29 @@ 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,
                                      analysis_errors = Symbol[],
                                      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)
 

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,10 @@ 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!

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,17 @@ 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 <:

src/callbacks_step/analysis_surface_2d.jl

@@ -0,0 +1,277 @@
+# 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