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