diff --git a/NEWS.md b/NEWS.md index 022252e61a9..5516fc9add5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,9 +7,10 @@ for human readability. ## Changes in the v0.7 lifecycle #### Added +- Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh`. - +- New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients. ## Changes when updating to v0.7 from v0.6.x diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl new file mode 100644 index 00000000000..fb5f29bd038 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -0,0 +1,132 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +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 * 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()) + + prim = SVector(rho_inf(), v1, v2, p_inf()) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach085_flow + +volume_flux = flux_ranocha_turbo +surface_flux = flux_lax_friedrichs + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +shock_indicator = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp", + joinpath(@__DIR__, "NACA0012.inp")) + +mesh = P4estMesh{2}(mesh_file) + +# The outer boundary is constant but subsonic, so we cannot compute the +# boundary flux for 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, + normal_direction::AbstractVector, x, + t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach085_flow(x, t, equations) + + return Trixi.flux_hll(u_inner, u_boundary, normal_direction, equations) +end + +boundary_conditions = Dict(:Left => boundary_condition_subsonic_constant, + :Right => boundary_condition_subsonic_constant, + :Top => boundary_condition_subsonic_constant, + :Bottom => boundary_condition_subsonic_constant, + :AirfoilBottom => boundary_condition_slip_wall, + :AirfoilTop => boundary_condition_slip_wall) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers + +# Run for a long time to reach a steady state +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 2000 + +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), l_inf)) + +lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, + LiftCoefficientPressure(aoa(), rho_inf(), + u_inf(equations), l_inf)) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "out", + save_analysis = true, + analysis_integrals = (drag_coefficient, + lift_coefficient)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 1, + med_level = 3, med_threshold = 0.05, + max_level = 4, max_threshold = 0.1) + +amr_interval = 100 +amr_callback = AMRCallback(semi, amr_controller, + interval = amr_interval, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback, amr_callback) + +############################################################################### +# run the simulation +sol = solve(ode, SSPRK54(thread = OrdinaryDiffEq.True()), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl new file mode 100644 index 00000000000..dc23e192de8 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -0,0 +1,125 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +@inline function initial_condition_mach038_flow(x, t, + equations::CompressibleEulerEquations2D) + # set the freestream flow parameters + rho_freestream = 1.4 + v1 = 0.38 + v2 = 0.0 + p_freestream = 1.0 + + prim = SVector(rho_freestream, v1, v2, p_freestream) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach038_flow + +volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be used +surface_flux = flux_lax_friedrichs + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +function mapping2cylinder(xi, eta) + xi_, eta_ = 0.5 * (xi + 1), 0.5 * (eta + 1.0) # Map from [-1,1] to [0,1] for simplicity + + R2 = 50.0 # Bigger circle + R1 = 0.5 # Smaller circle + + # Ensure an isotropic mesh by using elements with smaller radial length near the inner circle + + r = R1 * exp(xi_ * log(R2 / R1)) + theta = 2.0 * pi * eta_ + + x = r * cos(theta) + y = r * sin(theta) + return (x, y) +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 +# 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 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, + normal_direction::AbstractVector, x, + t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach038_flow(x, t, equations) + + return surface_flux_function(u_inner, u_boundary, normal_direction, equations) +end + +boundary_conditions = Dict(:x_neg => boundary_condition_slip_wall, + :x_pos => boundary_condition_subsonic_constant) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers + +# Run for a long time to reach a steady state +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 2000 + +aoa = 0.0 +rho_inf = 1.4 +u_inf = 0.38 +l_inf = 1.0 # Diameter of circle + +drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, + DragCoefficientPressure(aoa, rho_inf, u_inf, + l_inf)) + +lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, + LiftCoefficientPressure(aoa, rho_inf, u_inf, + l_inf)) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "out", + save_analysis = true, + analysis_integrals = (drag_coefficient, + lift_coefficient)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 2.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, + CarpenterKennedy2N54(williamson_condition = false; + thread = OrdinaryDiffEq.True()), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl new file mode 100644 index 00000000000..5bfa21c5322 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -0,0 +1,155 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +# Laminar transonic flow around a NACA0012 airfoil. + +# 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) +# - Deep Ray, Praveen Chandrashekar (2017) +# 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) + +equations = CompressibleEulerEquations2D(1.4) + +prandtl_number() = 0.72 +mu() = 0.0031959974968701088 +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), + Prandtl = prandtl_number()) + +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 + gamma = equations.gamma + u_inf = mach_inf() * sqrt(gamma * p_inf() / rho_inf()) + + v1 = u_inf * cos(aoa()) + v2 = u_inf * sin(aoa()) + + prim = SVector(rho_inf(), v1, v2, p_inf()) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach08_flow + +surface_flux = flux_lax_friedrichs + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux) + +mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp", + joinpath(@__DIR__, "NACA0012.inp")) + +mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 1) + +# The boundary values across outer boundary are 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, + normal_direction::AbstractVector, x, + t, + surface_flux_function, + equations::CompressibleEulerEquations2D) + u_boundary = initial_condition_mach08_flow(x, t, equations) + + return Trixi.flux_hll(u_inner, u_boundary, normal_direction, equations) +end + +boundary_conditions = Dict(:Left => boundary_condition_subsonic_constant, + :Right => boundary_condition_subsonic_constant, + :Top => boundary_condition_subsonic_constant, + :Bottom => boundary_condition_subsonic_constant, + :AirfoilBottom => boundary_condition_slip_wall, + :AirfoilTop => boundary_condition_slip_wall) + +velocity_airfoil = NoSlip((x, t, equations) -> SVector(0.0, 0.0)) + +heat_airfoil = Adiabatic((x, t, equations) -> 0.0) + +boundary_conditions_airfoil = BoundaryConditionNavierStokesWall(velocity_airfoil, + heat_airfoil) + +function momenta_initial_condition_mach08_flow(x, t, equations) + u = initial_condition_mach08_flow(x, t, equations) + momenta = SVector(u[2], u[3]) +end +velocity_bc_square = NoSlip((x, t, equations) -> momenta_initial_condition_mach08_flow(x, t, + equations)) + +heat_bc_square = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_square = BoundaryConditionNavierStokesWall(velocity_bc_square, + heat_bc_square) + +boundary_conditions_parabolic = Dict(:Left => boundary_condition_square, + :Right => boundary_condition_square, + :Top => boundary_condition_square, + :Bottom => boundary_condition_square, + :AirfoilBottom => boundary_conditions_airfoil, + :AirfoilTop => boundary_conditions_airfoil) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +############################################################################### +# ODE solvers + +# Run for a long time to reach a state where forces stabilize up to 3 digits +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 2000 + +force_boundary_names = [:AirfoilBottom, :AirfoilTop] +drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, + DragCoefficientPressure(aoa(), rho_inf(), + u_inf(equations), + l_inf())) + +lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, + LiftCoefficientPressure(aoa(), 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)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1e-8, + reltol = 1e-8, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index 9375c80d77e..168441513ad 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -261,7 +261,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, AveragingCallback, AMRCallback, StepsizeCallback, GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback, - TrivialCallback, AnalysisCallbackCoupled + TrivialCallback, AnalysisCallbackCoupled, + AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure 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 ba232032951..62048853b53 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -691,6 +691,7 @@ end # @muladd # specialized implementations specific to some solvers include("analysis_dg1d.jl") include("analysis_dg2d.jl") +include("analysis_surface_integral_2d.jl") include("analysis_dg2d_parallel.jl") include("analysis_dg3d.jl") include("analysis_dg3d_parallel.jl") diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl new file mode 100644 index 00000000000..902fb0b4e85 --- /dev/null +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -0,0 +1,189 @@ +# 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 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 with the boundary + name `:Airfoil` in 2D. +""" +struct AnalysisSurfaceIntegral{Semidiscretization, Variable} + semi::Semidiscretization # passed in to retrieve boundary condition information + indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed + variable::Variable # Quantity of interest, like lift or drag + + function AnalysisSurfaceIntegral(semi, boundary_symbol, variable) + @unpack boundary_symbol_indices = semi.boundary_conditions + indices = boundary_symbol_indices[boundary_symbol] + + return new{typeof(semi), typeof(variable)}(semi, indices, + variable) + end + + function AnalysisSurfaceIntegral(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(semi), typeof(variable)}(semi, indices, + variable) + end +end + +struct ForceState{RealT <: Real} + psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream + rhoinf::RealT + uinf::RealT + l_inf::RealT +end + +struct LiftCoefficientPressure{RealT <: Real} + force_state::ForceState{RealT} +end + +struct DragCoefficientPressure{RealT <: Real} + force_state::ForceState{RealT} +end + +""" + LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) + +Compute the lift coefficient +```math +C_{L,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_L \\, \\mathrm{d} S} + {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} +``` +based on the pressure distribution along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) + # 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. + # 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, l_inf)) +end + +""" + DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) + +Compute the drag coefficient +```math +C_{D,p} \\coloneqq \\frac{\\oint_{\\partial \\Omega} p \\boldsymbol n \\cdot \\psi_D \\, \\mathrm{d} S} + {0.5 \\cdot \\rho_{\\infty} \\cdot U_{\\infty}^2 \\cdot L_{\\infty}} +``` +based on the pressure distribution along a boundary. +Supposed to be used in conjunction with [`AnalysisSurfaceIntegral`](@ref) +which stores the boundary information and semidiscretization. + +- `aoa::Real`: Angle of attack in radians (for airfoils etc.) +- `rhoinf::Real`: Free-stream density +- `uinf::Real`: Free-stream velocity +- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) +""" +function DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) + # `psi_drag` is the unit vector tangent to the freestream direction + psi_drag = (cos(aoa), sin(aoa)) + return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, l_inf)) +end + +function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations) + p = pressure(u, equations) + @unpack psi, rhoinf, uinf, l_inf = lift_coefficient.force_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 * l_inf) +end + +function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations) + p = pressure(u, equations) + @unpack psi, rhoinf, uinf, l_inf = drag_coefficient.force_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 * l_inf) +end + +function analyze(surface_variable::AnalysisSurfaceIntegral, 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 + + surface_integral = zero(eltype(u)) + index_range = eachnode(dg) + 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) + # 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) + + # L2 norm of normal direction (contravariant_vector) is the surface element + dS = weights[node_index] * norm(normal_direction) + # Integral over entire boundary surface + surface_integral += variable(u_node, normal_direction, equations) * dS + + i_node += i_node_step + j_node += j_node_step + end + end + return surface_integral +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, + <:LiftCoefficientPressure{<:Any}}) + "CL_p" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, + <:LiftCoefficientPressure{<:Any}}) + "CL_p" +end + +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, + <:DragCoefficientPressure{<:Any}}) + "CD_p" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, + <:DragCoefficientPressure{<:Any}}) + "CD_p" +end +end # muladd diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index ec50627d3ef..10cc075089c 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -36,7 +36,8 @@ end orientation = (direction + 1) >> 1 normal = get_contravariant_vector(orientation, contravariant_vectors, indices...) - # Contravariant vectors at interfaces in negative coordinate direction are pointing inwards + # Contravariant vectors at interfaces in negative coordinate direction are pointing inwards, + # flip sign to make them point outwards if isodd(direction) return -normal else diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index b5388cadc8b..2c2c6876d70 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -8,8 +8,8 @@ """ UnstructuredSortedBoundaryTypes -General container to sort the boundary conditions by type for some unstructured meshes/solvers. -It stores a set of global indices for each boundary condition type to expedite computation +General container to sort the boundary conditions by type and name for some unstructured meshes/solvers. +It stores a set of global indices for each boundary condition type and name to expedite computation during the call to `calc_boundary_flux!`. The original dictionary form of the boundary conditions set by the user in the elixir file is also stored for printing. """ @@ -17,6 +17,7 @@ mutable struct UnstructuredSortedBoundaryTypes{N, BCs <: NTuple{N, Any}} boundary_condition_types::BCs # specific boundary condition type(s), e.g. BoundaryConditionDirichlet boundary_indices::NTuple{N, Vector{Int}} # integer vectors containing global boundary indices boundary_dictionary::Dict{Symbol, Any} # boundary conditions as set by the user in the elixir file + boundary_symbol_indices::Dict{Symbol, Vector{Int}} # integer vectors containing global boundary indices per boundary identifier end # constructor that "eats" the original boundary condition dictionary and sorts the information @@ -28,10 +29,14 @@ function UnstructuredSortedBoundaryTypes(boundary_conditions::Dict, cache) n_boundary_types = length(boundary_condition_types) boundary_indices = ntuple(_ -> [], n_boundary_types) + # Initialize `boundary_symbol_indices` as an empty dictionary, filled later in `initialize!` + boundary_symbol_indices = Dict{Symbol, Vector{Int}}() + container = UnstructuredSortedBoundaryTypes{n_boundary_types, typeof(boundary_condition_types)}(boundary_condition_types, boundary_indices, - boundary_conditions) + boundary_conditions, + boundary_symbol_indices) initialize!(container, cache) end @@ -97,6 +102,13 @@ function initialize!(boundary_types_container::UnstructuredSortedBoundaryTypes{N # convert the work array with the boundary indices into a tuple boundary_types_container.boundary_indices = Tuple(_boundary_indices) + # Store boundary indices per symbol (required for force computations, for instance) + for (symbol, _) in boundary_dictionary + indices = findall(x -> x === symbol, cache.boundaries.name) + # Store the indices in `boundary_symbol_indices` dictionary + boundary_types_container.boundary_symbol_indices[symbol] = sort!(indices) + end + return boundary_types_container end end # @muladd diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 121001b35ff..1bacccbf812 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -531,6 +531,86 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_subsonic_cylinder.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_subsonic_cylinder.jl"), + l2=[ + 0.00011914390523852561, + 0.00010776028621724485, + 6.139954358305467e-5, + 0.0003067693731825959, + ], + linf=[ + 0.1653075586200805, + 0.1868437275544909, + 0.09772818519679008, + 0.4311796171737692, + ], tspan=(0.0, 0.001)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + 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 + end + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + + @test isapprox(lift, -6.501138753497174e-15, atol = 1e-13) + @test isapprox(drag, 2.588589856781827, atol = 1e-13) +end + +# Forces computation test in an AMR code +@trixi_testset "elixir_euler_NACA0012airfoil_mach085.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_NACA0012airfoil_mach085.jl"), + l2=[ + 5.371568111383228e-7, 6.4158131303956445e-6, + 1.0324346542348325e-5, 0.0006348064933187732, + ], + linf=[ + 0.0016263400091978443, 0.028471072159724428, + 0.02986133204785877, 1.9481060511014872, + ], + base_level=0, med_level=1, max_level=1, + tspan=(0.0, 0.0001), + adapt_initial_condition=false, + adapt_initial_condition_only_refine=false, + # With the default `maxiters = 1` in coverage tests, + # the values for `drag` and `lift` below would differ. + coverage_override=(maxiters = 100_000,)) + + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + 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 + end + + u_ode = copy(sol.u[end]) + du_ode = zero(u_ode) # Just a placeholder in this case + + u = Trixi.wrap_array(u_ode, semi) + du = Trixi.wrap_array(du_ode, semi) + drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, + semi.cache) + + @test isapprox(lift, 0.0262382560809345, atol = 1e-13) + @test isapprox(drag, 0.10898248971932244, atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index d47c34f9e75..edbdfe5a84f 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -707,6 +707,28 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_navierstokes_NACA0012airfoil_mach08.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_navierstokes_NACA0012airfoil_mach08.jl"), + l2=[0.000186486564226516, + 0.0005076712323400374, + 0.00038074588984354107, + 0.002128177239782089], + linf=[0.5153387072802718, + 1.199362305026636, + 0.9077214424040279, + 5.666071182328691], tspan=(0.0, 0.001), + initial_refinement_level=0) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + 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 + end +end end # Clean up afterwards: delete Trixi.jl output directory