From 99d0659416ff66de576b73019de5537e8741f509 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 16:26:53 +0530 Subject: [PATCH 01/53] Add AnalysisSurfaceIntegral --- .../elixir_euler_subsonic_cylinder.jl | 122 ++++++++++++++++++ src/callbacks_step/analysis.jl | 1 + src/callbacks_step/analysis_surface_2d.jl | 116 +++++++++++++++++ 3 files changed, 239 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl create mode 100644 src/callbacks_step/analysis_surface_2d.jl 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..e6b438944f3 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -0,0 +1,122 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi +# include("analysis_surface_2d.jl") + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +@inline function initial_condition_mach01_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_mach01_flow + +volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be used +surface_flux = flux_lax_friedrichs + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +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 = 3, + periodicity = (false, true)) + +# The boundary of the outer cylinder is constant but supersonic, 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_mach01_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 state 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 +linf = 1.0 # Diameter of circle + +indices = semi_ -> semi.boundary_conditions.boundary_indices[2] +my_drag_force = Trixi.AnalysisSurfaceIntegral(indices, + Trixi.DragForcePressure(aoa, rho_inf, U_inf, + linf)) + +my_lift_force = Trixi.AnalysisSurfaceIntegral(indices, + Trixi.LiftForcePressure(aoa, rho_inf, U_inf, + linf)) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + analysis_errors = Symbol[], + output_directory = "analysis_results", + save_analysis = true, + analysis_integrals = (my_drag_force, my_lift_force)) + +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, SSPRK43(); + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index ba232032951..2c3351272fa 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_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_2d.jl new file mode 100644 index 00000000000..f984849448b --- /dev/null +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -0,0 +1,116 @@ +using Trixi +using Trixi: integrate_via_indices, norm, apply_jacobian_parabolic!, @threaded, + indices2direction, + index_to_start_step_2d, get_normal_direction, dot, get_node_coords +import Trixi: analyze, pretty_form_ascii, pretty_form_utf + +struct AnalysisSurfaceIntegral{Indices, Variable} + indices::Indices + variable::Variable +end + +struct ForceState{RealT <: Real} + Ψl::Tuple{RealT, RealT} + rhoinf::RealT + uinf::RealT + linf::RealT +end + +# TODO - This should be a struct in ForceState +struct FreeStreamVariables{RealT <: Real} + rhoinf::RealT + uinf::RealT + linf::RealT +end + +struct LiftForcePressure{RealT <: Real} + force_state::ForceState{RealT} +end + +struct DragForcePressure{RealT <: Real} + force_state::ForceState{RealT} +end + +function LiftForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real) + Ψl = (-sin(aoa), cos(aoa)) + force_state = ForceState(Ψl, rhoinf, uinf, linf) + return LiftForcePressure(force_state) +end + +function DragForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real) + Ψd = (cos(aoa), sin(aoa)) + return DragForcePressure(ForceState(Ψd, rhoinf, uinf, linf)) +end + +function (lift_force::LiftForcePressure)(u, normal_direction, equations) + p = pressure(u, equations) + @unpack Ψl, rhoinf, uinf, linf = lift_force.force_state + n = dot(normal_direction, Ψl) / norm(normal_direction) + return p * n / (0.5 * rhoinf * uinf^2 * linf) +end + +function (drag_force::DragForcePressure)(u, normal_direction, equations) + p = pressure(u, equations) + @unpack Ψl, rhoinf, uinf, linf = drag_force.force_state + n = dot(normal_direction, Ψl) / norm(normal_direction) + return p * n / (0.5 * rhoinf * uinf^2 * linf) +end + +function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, + mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, 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 + # TODO - Use initialize callbacks to move boundary_conditions to cache + indices_ = indices(cache) + + surface_integral = zero(eltype(u)) + index_range = eachnode(dg) + for local_index in eachindex(indices_) + # Use the local index to get the global boundary index from the pre-sorted list + boundary = indices_[local_index] + + # Get information on the adjacent element, compute the surface fluxes, + # and store them + 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 eachnode(dg) + u_node = Trixi.get_node_vars(cache.boundaries.u, equations, dg, node_index, + boundary) + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, + element) + + # L2 norm of normal direction is the surface element + # 0.5 factor is NOT needed, the norm(normal_direction) is all the factor needed + dS = weights[node_index] * norm(normal_direction) + 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, <:LiftForcePressure{<:Any}}) + "Pressure_lift" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:LiftForcePressure{<:Any}}) + "Pressure_lift" +end +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:DragForcePressure{<:Any}}) + "Pressure_drag" +end +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:DragForcePressure{<:Any}}) + "Pressure_drag" +end From aee4cfa39933e0e7558f1d9831b0c4a6d964d417 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 16:34:55 +0530 Subject: [PATCH 02/53] Minor change --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index e6b438944f3..e257edbceda 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -1,7 +1,6 @@ using Downloads: download using OrdinaryDiffEq using Trixi -# include("analysis_surface_2d.jl") ############################################################################### # semidiscretization of the compressible Euler equations @@ -91,7 +90,7 @@ rho_inf = 1.4 U_inf = 0.38 linf = 1.0 # Diameter of circle -indices = semi_ -> semi.boundary_conditions.boundary_indices[2] +indices = semi_ -> semi.boundary_conditions.boundary_indices[2] # TODO - Really needs fixing! my_drag_force = Trixi.AnalysisSurfaceIntegral(indices, Trixi.DragForcePressure(aoa, rho_inf, U_inf, linf)) From 9f8daef88900f3adfb6ba79a5e680efe10542200 Mon Sep 17 00:00:00 2001 From: Arpit Babbar Date: Wed, 24 Jan 2024 18:19:58 +0530 Subject: [PATCH 03/53] Update elixir_euler_subsonic_cylinder.jl --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index e257edbceda..678905c3b3c 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -53,7 +53,7 @@ cells_per_dimension = (64, 64) mesh = P4estMesh(cells_per_dimension, mapping = mapping2cylinder, polydeg = 3, periodicity = (false, true)) -# The boundary of the outer cylinder is constant but supersonic, so we cannot compute the +# The boundary of the outer cylinder 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, From 4054ed5fd1033fe3d936f011e0a7294b8f099388 Mon Sep 17 00:00:00 2001 From: Arpit Babbar Date: Wed, 24 Jan 2024 19:30:36 +0530 Subject: [PATCH 04/53] Apply suggestions from code review Co-authored-by: Daniel Doehring --- src/callbacks_step/analysis_surface_2d.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index f984849448b..1f02ba3f0a9 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -1,8 +1,3 @@ -using Trixi -using Trixi: integrate_via_indices, norm, apply_jacobian_parabolic!, @threaded, - indices2direction, - index_to_start_step_2d, get_normal_direction, dot, get_node_coords -import Trixi: analyze, pretty_form_ascii, pretty_form_utf struct AnalysisSurfaceIntegral{Indices, Variable} indices::Indices @@ -32,12 +27,14 @@ struct DragForcePressure{RealT <: Real} end function LiftForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real) + # Ψl is the normal unit vector to the freestream direction Ψl = (-sin(aoa), cos(aoa)) force_state = ForceState(Ψl, rhoinf, uinf, linf) return LiftForcePressure(force_state) end function DragForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real) + # Ψd is the unit vector parallel to the freestream direction Ψd = (cos(aoa), sin(aoa)) return DragForcePressure(ForceState(Ψd, rhoinf, uinf, linf)) end From 6b208c8cebe79d60d1a61cd36207f585c947c6a0 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 21:08:28 +0530 Subject: [PATCH 05/53] Fix labels, add tests --- .../elixir_euler_subsonic_cylinder.jl | 16 +++---- src/callbacks_step/analysis_surface_2d.jl | 42 +++++++++---------- test/test_p4est_2d.jl | 27 ++++++++++++ 3 files changed, 56 insertions(+), 29 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 678905c3b3c..93d4f4eaca6 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -2,6 +2,8 @@ using Downloads: download using OrdinaryDiffEq using Trixi +using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient + ############################################################################### # semidiscretization of the compressible Euler equations @@ -91,19 +93,17 @@ U_inf = 0.38 linf = 1.0 # Diameter of circle indices = semi_ -> semi.boundary_conditions.boundary_indices[2] # TODO - Really needs fixing! -my_drag_force = Trixi.AnalysisSurfaceIntegral(indices, - Trixi.DragForcePressure(aoa, rho_inf, U_inf, - linf)) -my_lift_force = Trixi.AnalysisSurfaceIntegral(indices, - Trixi.LiftForcePressure(aoa, rho_inf, U_inf, - linf)) +drag_coefficient = AnalysisSurfaceIntegral(indices, + DragCoefficient(aoa, rho_inf, U_inf, linf)) +lift_coefficient = AnalysisSurfaceIntegral(indices, + LiftCoefficient(aoa, rho_inf, U_inf, linf)) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - analysis_errors = Symbol[], output_directory = "analysis_results", save_analysis = true, - analysis_integrals = (my_drag_force, my_lift_force)) + analysis_integrals = (drag_coefficient, + lift_coefficient)) alive_callback = AliveCallback(analysis_interval = analysis_interval) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index 1f02ba3f0a9..953b2488b41 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -5,7 +5,7 @@ struct AnalysisSurfaceIntegral{Indices, Variable} end struct ForceState{RealT <: Real} - Ψl::Tuple{RealT, RealT} + Ψ::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream rhoinf::RealT uinf::RealT linf::RealT @@ -18,38 +18,38 @@ struct FreeStreamVariables{RealT <: Real} linf::RealT end -struct LiftForcePressure{RealT <: Real} +struct LiftCoefficient{RealT <: Real} force_state::ForceState{RealT} end -struct DragForcePressure{RealT <: Real} +struct DragCoefficient{RealT <: Real} force_state::ForceState{RealT} end -function LiftForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real) +function LiftCoefficient(aoa, rhoinf, uinf, linf) # Ψl is the normal unit vector to the freestream direction Ψl = (-sin(aoa), cos(aoa)) force_state = ForceState(Ψl, rhoinf, uinf, linf) - return LiftForcePressure(force_state) + return LiftCoefficient(force_state) end -function DragForcePressure(aoa::Real, rhoinf::Real, uinf::Real, linf::Real) +function DragCoefficient(aoa, rhoinf, uinf, linf) # Ψd is the unit vector parallel to the freestream direction Ψd = (cos(aoa), sin(aoa)) - return DragForcePressure(ForceState(Ψd, rhoinf, uinf, linf)) + return DragCoefficient(ForceState(Ψd, rhoinf, uinf, linf)) end -function (lift_force::LiftForcePressure)(u, normal_direction, equations) +function (lift_coefficient::LiftCoefficient)(u, normal_direction, equations) p = pressure(u, equations) - @unpack Ψl, rhoinf, uinf, linf = lift_force.force_state - n = dot(normal_direction, Ψl) / norm(normal_direction) + @unpack Ψ, rhoinf, uinf, linf = lift_coefficient.force_state + n = dot(normal_direction, Ψ) / norm(normal_direction) return p * n / (0.5 * rhoinf * uinf^2 * linf) end -function (drag_force::DragForcePressure)(u, normal_direction, equations) +function (drag_coefficient::DragCoefficient)(u, normal_direction, equations) p = pressure(u, equations) - @unpack Ψl, rhoinf, uinf, linf = drag_force.force_state - n = dot(normal_direction, Ψl) / norm(normal_direction) + @unpack Ψ, rhoinf, uinf, linf = drag_coefficient.force_state + n = dot(normal_direction, Ψ) / norm(normal_direction) return p * n / (0.5 * rhoinf * uinf^2 * linf) end @@ -99,15 +99,15 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, return surface_integral end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:LiftForcePressure{<:Any}}) - "Pressure_lift" +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:LiftCoefficient{<:Any}}) + "CL" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:LiftForcePressure{<:Any}}) - "Pressure_lift" +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:LiftCoefficient{<:Any}}) + "CL" end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:DragForcePressure{<:Any}}) - "Pressure_drag" +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:DragCoefficient{<:Any}}) + "CD" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:DragForcePressure{<:Any}}) - "Pressure_drag" +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:DragCoefficient{<:Any}}) + "CD" end diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 121001b35ff..66931f5c358 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -531,6 +531,33 @@ 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.0001191438553219565, + 0.00010775987550586874, + 6.139942108660038e-5, + 0.0003067690221489065], + linf=[ + 0.16530830127928575, + 0.18684298318465908, + 0.09772804813504012, + 0.43118129113312076], tspan=(0.0, 0.001)) + + 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, 2.6395869624301843e-15, atol = 1e-13) + @test isapprox(drag, 2.5885879058648857, atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory From c89e877980bb23fd3c991e7db180a3bd85e7cfd7 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 22:11:32 +0530 Subject: [PATCH 06/53] Attempt to fix CI --- .../elixir_euler_subsonic_cylinder.jl | 11 +++++++--- test/test_p4est_2d.jl | 22 ++++++++++--------- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 93d4f4eaca6..72e51e28f76 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -112,10 +112,15 @@ save_solution = SaveSolutionCallback(interval = 500, save_final_solution = true, solution_variables = cons2prim) -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) +# Small time step should be used to reach steady state +stepsize_callback = StepsizeCallback(cfl = 0.25) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) ############################################################################### # run the simulation -sol = solve(ode, SSPRK43(); - ode_default_options()..., callback = callbacks); +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + 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/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 66931f5c358..0a37d14d8f0 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -535,15 +535,17 @@ end @trixi_testset "elixir_euler_subsonic_cylinder.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_subsonic_cylinder.jl"), l2=[ - 0.0001191438553219565, - 0.00010775987550586874, - 6.139942108660038e-5, - 0.0003067690221489065], + 0.00011914390523852561, + 0.00010776028621724485, + 6.139954358305467e-5, + 0.0003067693731825959, + ], linf=[ - 0.16530830127928575, - 0.18684298318465908, - 0.09772804813504012, - 0.43118129113312076], tspan=(0.0, 0.001)) + 0.1653075586200805, + 0.1868437275544909, + 0.09772818519679008, + 0.4311796171737692, + ], tspan=(0.0, 0.001)) u_ode = copy(sol.u[end]) du_ode = zero(u_ode) # Just a placeholder in this case @@ -555,8 +557,8 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 2.6395869624301843e-15, atol = 1e-13) - @test isapprox(drag, 2.5885879058648857, atol = 1e-13) + @test isapprox(lift, -6.501138753497174e-15, atol = 1e-13) + @test isapprox(drag, 2.588589856781827, atol = 1e-13) end end From 5b24e2e999ca2656da525e4dca5ee2761c613cb5 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 25 Jan 2024 08:14:52 +0530 Subject: [PATCH 07/53] Obtain indices from user chosen function --- .../elixir_euler_subsonic_cylinder.jl | 9 ++-- src/callbacks_step/analysis_surface_2d.jl | 45 ++++++++++++++----- 2 files changed, 37 insertions(+), 17 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 72e51e28f76..8901d5148b1 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -77,7 +77,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, ############################################################################### # ODE solvers -# Run for a long time to reach state state +# Run for a long time to reach a steady state tspan = (0.0, 100.0) ode = semidiscretize(semi, tspan) @@ -92,11 +92,10 @@ rho_inf = 1.4 U_inf = 0.38 linf = 1.0 # Diameter of circle -indices = semi_ -> semi.boundary_conditions.boundary_indices[2] # TODO - Really needs fixing! - -drag_coefficient = AnalysisSurfaceIntegral(indices, +drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, DragCoefficient(aoa, rho_inf, U_inf, linf)) -lift_coefficient = AnalysisSurfaceIntegral(indices, + +lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, LiftCoefficient(aoa, rho_inf, U_inf, linf)) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index 953b2488b41..c9669e42f9e 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -1,7 +1,28 @@ - -struct AnalysisSurfaceIntegral{Indices, Variable} +# 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 +# TODO - Add explanation for what this file contains + +struct AnalysisSurfaceIntegral{SemiDiscretization, Indices, Variable} + semi::SemiDiscretization indices::Indices - variable::Variable + variable::Variable # Quantity of interest, like lift or drag + function AnalysisSurfaceIntegral(semi, boundary_condition_type, variable) + # The bc list as ordered in digest_boundary_conditions + ordered_bc = semi.boundary_conditions.boundary_condition_types + + # The set of all indices that gives the bc where the surface integral is to be computed + index = sort(findall(x->x==boundary_condition_type, ordered_bc)) + + # Put the bc in function form as they might change under AMR + indices = semi -> Vector([semi.boundary_conditions.boundary_indices[i] + for i in index])[1] # TODO - Should not need Vector and the [1] + + return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, variable) + end end struct ForceState{RealT <: Real} @@ -54,14 +75,13 @@ function (drag_coefficient::DragCoefficient)(u, normal_direction, equations) end function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, + 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 - # TODO - Use initialize callbacks to move boundary_conditions to cache - indices_ = indices(cache) + @unpack semi, indices, variable = surface_variable + indices_ = indices(semi) surface_integral = zero(eltype(u)) index_range = eachnode(dg) @@ -88,7 +108,6 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, element) # L2 norm of normal direction is the surface element - # 0.5 factor is NOT needed, the norm(normal_direction) is all the factor needed dS = weights[node_index] * norm(normal_direction) surface_integral += variable(u_node, normal_direction, equations) * dS @@ -99,15 +118,17 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, return surface_integral end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:LiftCoefficient{<:Any}}) +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, <:LiftCoefficient{<:Any}}) "CL" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:LiftCoefficient{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, <:LiftCoefficient{<:Any}}) "CL" end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:DragCoefficient{<:Any}}) +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, <:DragCoefficient{<:Any}}) "CD" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:DragCoefficient{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, <:DragCoefficient{<:Any}}) "CD" end + +end # muladd \ No newline at end of file From 573ffc3d1a284ee7335c95c47b66670b637e3c13 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 25 Jan 2024 08:18:21 +0530 Subject: [PATCH 08/53] Formatting --- src/callbacks_step/analysis_surface_2d.jl | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index c9669e42f9e..9aa6b303866 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -15,13 +15,14 @@ struct AnalysisSurfaceIntegral{SemiDiscretization, Indices, Variable} ordered_bc = semi.boundary_conditions.boundary_condition_types # The set of all indices that gives the bc where the surface integral is to be computed - index = sort(findall(x->x==boundary_condition_type, ordered_bc)) + index = sort(findall(x -> x == boundary_condition_type, ordered_bc)) # Put the bc in function form as they might change under AMR indices = semi -> Vector([semi.boundary_conditions.boundary_indices[i] - for i in index])[1] # TODO - Should not need Vector and the [1] + for i in index])[1] # TODO - Should not need Vector and the [1] - return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, variable) + return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, + variable) end end @@ -118,17 +119,20 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, return surface_integral end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, <:LiftCoefficient{<:Any}}) +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, + <:LiftCoefficient{<:Any}}) "CL" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, <:LiftCoefficient{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, + <:LiftCoefficient{<:Any}}) "CL" end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, <:DragCoefficient{<:Any}}) +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, + <:DragCoefficient{<:Any}}) "CD" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, <:DragCoefficient{<:Any}}) +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, + <:DragCoefficient{<:Any}}) "CD" end - -end # muladd \ No newline at end of file +end # muladd From f565df6fe44ba14d9bf4847c35537ff662f5520a Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 16:34:55 +0530 Subject: [PATCH 09/53] Minor change --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 8901d5148b1..873d31e7aab 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -92,8 +92,15 @@ rho_inf = 1.4 U_inf = 0.38 linf = 1.0 # Diameter of circle +<<<<<<< HEAD drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, DragCoefficient(aoa, rho_inf, U_inf, linf)) +======= +indices = semi_ -> semi.boundary_conditions.boundary_indices[2] # TODO - Really needs fixing! +my_drag_force = Trixi.AnalysisSurfaceIntegral(indices, + Trixi.DragForcePressure(aoa, rho_inf, U_inf, + linf)) +>>>>>>> d89082be (Minor change) lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, LiftCoefficient(aoa, rho_inf, U_inf, linf)) From 9c0d5e1e0042108e90f4289673f23acee5c19f9c Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 21:08:28 +0530 Subject: [PATCH 10/53] Fix labels, add tests --- .../p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 873d31e7aab..a5074511dc8 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -92,15 +92,8 @@ rho_inf = 1.4 U_inf = 0.38 linf = 1.0 # Diameter of circle -<<<<<<< HEAD drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, DragCoefficient(aoa, rho_inf, U_inf, linf)) -======= -indices = semi_ -> semi.boundary_conditions.boundary_indices[2] # TODO - Really needs fixing! -my_drag_force = Trixi.AnalysisSurfaceIntegral(indices, - Trixi.DragForcePressure(aoa, rho_inf, U_inf, - linf)) ->>>>>>> d89082be (Minor change) lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, LiftCoefficient(aoa, rho_inf, U_inf, linf)) @@ -129,4 +122,4 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), 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 +summary_callback() # print the timer summary \ No newline at end of file From 6de00ddf28d7e384d727a4e91ad99eb60b1890d2 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 22:11:32 +0530 Subject: [PATCH 11/53] Attempt to fix CI --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index a5074511dc8..8901d5148b1 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -122,4 +122,4 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), 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 \ No newline at end of file +summary_callback() # print the timer summary From 8ada85c8c75b0c1bc81edf20dcc34764b0ceb6da Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 25 Jan 2024 08:14:52 +0530 Subject: [PATCH 12/53] Obtain indices from user chosen function --- src/callbacks_step/analysis_surface_2d.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index 9aa6b303866..aa0a00a07b8 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -15,6 +15,7 @@ struct AnalysisSurfaceIntegral{SemiDiscretization, Indices, Variable} ordered_bc = semi.boundary_conditions.boundary_condition_types # The set of all indices that gives the bc where the surface integral is to be computed +<<<<<<< HEAD index = sort(findall(x -> x == boundary_condition_type, ordered_bc)) # Put the bc in function form as they might change under AMR @@ -23,6 +24,15 @@ struct AnalysisSurfaceIntegral{SemiDiscretization, Indices, Variable} return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, variable) +======= + index = sort(findall(x->x==boundary_condition_type, ordered_bc)) + + # Put the bc in function form as they might change under AMR + indices = semi -> Vector([semi.boundary_conditions.boundary_indices[i] + for i in index])[1] # TODO - Should not need Vector and the [1] + + return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, variable) +>>>>>>> 5d22bb6d (Obtain indices from user chosen function) end end From a87ac8f916988b9df2f17d455afe7ecffd02ea5e Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 25 Jan 2024 08:18:21 +0530 Subject: [PATCH 13/53] Add new elixirs --- .../elixir_euler_NACA0012airfoil_mach08.jl | 139 +++++++++++++++++ ...xir_navierstokes_NACA0012airfoil_mach08.jl | 142 ++++++++++++++++++ src/callbacks_step/analysis_surface_2d.jl | 16 +- 3 files changed, 284 insertions(+), 13 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach08.jl create mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach08.jl new file mode 100644 index 00000000000..accd79d1187 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach08.jl @@ -0,0 +1,139 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient + +############################################################################### +# semidiscretization of the compressible Euler equations + +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 +mach_inf() = 0.85 +aoa() = pi / 180.0 +c_inf(equations) = sqrt(equations.gamma * pre_inf() / rho_inf()) +U_inf(equations) = mach_inf() * c_inf(equations) + +@inline function initial_condition_mach085_flow(x, t, + equations::CompressibleEulerEquations2D) + # set the freestream flow parameters + gasGam = equations.gamma + + v1 = U_inf(equations) * cos(aoa()) + v2 = U_inf(equations) * sin(aoa()) + + prim = SVector(rho_inf(), v1, v2, pre_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("", + joinpath(@__DIR__, "NACA0012.inp")) +=# +mesh_file = "NACA0012.inp" + +mesh = P4estMesh{2}(mesh_file) + +# The boundary of the outer cylinder 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 + +linf = 1.0 # Length of airfoil + +drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, + DragCoefficient(aoa(), rho_inf(), + U_inf(equations), linf)) + +lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, + LiftCoefficient(aoa(), rho_inf(), + U_inf(equations), linf)) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "analysis_results", + 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) + +# Small time step should be used to reach steady state +stepsize_callback = StepsizeCallback(cfl = 0.25) + +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_callback = AMRCallback(semi, amr_controller, + interval = 100, + 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, CarpenterKennedy2N54(williamson_condition = false), + 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..178effb2db4 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -0,0 +1,142 @@ +using Downloads: download +using OrdinaryDiffEq +using Trixi + +using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +prandtl_number() = 0.72 +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()) +@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) + + v1 = U_inf * cos(aoa) + v2 = U_inf * sin(aoa) + + prim = SVector(rho_inf, v1, v2, pre_inf) + return prim2cons(prim, equations) +end + +initial_condition = initial_condition_mach08_flow + +surface_flux = flux_lax_friedrichs + +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux) + +#= +mesh_file = Trixi.download("", + joinpath(@__DIR__, "NACA0012.inp")) +=# +mesh_file = "NACA0012.inp" + +mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 1) + +# The boundary of the outer cylinder 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_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) + +velocity_bc_square = NoSlip((x, t, equations) -> initial_condition_mach08_flow(x, t, + equations)[2:3]) +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 steady state +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 2000 + +drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, + DragCoefficient(sw_aoa(), sw_rho_inf(), + sw_U_inf(equations), sw_linf())) + +lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, + LiftCoefficient(sw_aoa(), sw_rho_inf(), + sw_U_inf(equations), sw_linf())) + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + output_directory = "analysis_results", + 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) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +time_int_tol = 1e-11 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index aa0a00a07b8..cb667ee1f9e 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -6,8 +6,8 @@ #! format: noindent # TODO - Add explanation for what this file contains -struct AnalysisSurfaceIntegral{SemiDiscretization, Indices, Variable} - semi::SemiDiscretization +struct AnalysisSurfaceIntegral{Semidiscretization, Indices, Variable} + semi::Semidiscretization indices::Indices variable::Variable # Quantity of interest, like lift or drag function AnalysisSurfaceIntegral(semi, boundary_condition_type, variable) @@ -15,7 +15,6 @@ struct AnalysisSurfaceIntegral{SemiDiscretization, Indices, Variable} ordered_bc = semi.boundary_conditions.boundary_condition_types # The set of all indices that gives the bc where the surface integral is to be computed -<<<<<<< HEAD index = sort(findall(x -> x == boundary_condition_type, ordered_bc)) # Put the bc in function form as they might change under AMR @@ -24,15 +23,6 @@ struct AnalysisSurfaceIntegral{SemiDiscretization, Indices, Variable} return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, variable) -======= - index = sort(findall(x->x==boundary_condition_type, ordered_bc)) - - # Put the bc in function form as they might change under AMR - indices = semi -> Vector([semi.boundary_conditions.boundary_indices[i] - for i in index])[1] # TODO - Should not need Vector and the [1] - - return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, variable) ->>>>>>> 5d22bb6d (Obtain indices from user chosen function) end end @@ -145,4 +135,4 @@ function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, <:DragCoefficient{<:Any}}) "CD" end -end # muladd +end # muladd \ No newline at end of file From 4ac2a21e2da52cab167f3bcddefb99d3af82c2a5 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Sat, 10 Feb 2024 21:49:50 +0530 Subject: [PATCH 14/53] Add tests for NACA0012 --- ....jl => elixir_euler_NACA0012airfoil_mach085.jl} | 10 ++-------- .../elixir_navierstokes_NACA0012airfoil_mach08.jl | 5 +---- test/test_p4est_2d.jl | 13 +++++++++++++ test/test_parabolic_2d.jl | 14 ++++++++++++++ 4 files changed, 30 insertions(+), 12 deletions(-) rename examples/p4est_2d_dgsem/{elixir_euler_NACA0012airfoil_mach08.jl => elixir_euler_NACA0012airfoil_mach085.jl} (95%) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl similarity index 95% rename from examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach08.jl rename to examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index accd79d1187..b7568296602 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -18,9 +18,6 @@ U_inf(equations) = mach_inf() * c_inf(equations) @inline function initial_condition_mach085_flow(x, t, equations::CompressibleEulerEquations2D) - # set the freestream flow parameters - gasGam = equations.gamma - v1 = U_inf(equations) * cos(aoa()) v2 = U_inf(equations) * sin(aoa()) @@ -46,15 +43,12 @@ volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) -#= -mesh_file = Trixi.download("", +mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp", joinpath(@__DIR__, "NACA0012.inp")) -=# -mesh_file = "NACA0012.inp" mesh = P4estMesh{2}(mesh_file) -# The boundary of the outer cylinder is constant but subsonic, so we cannot compute the +# 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, diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 178effb2db4..c15133de4e0 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -44,11 +44,8 @@ polydeg = 3 basis = LobattoLegendreBasis(polydeg) solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux) -#= -mesh_file = Trixi.download("", +mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp", joinpath(@__DIR__, "NACA0012.inp")) -=# -mesh_file = "NACA0012.inp" mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 1) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 0a37d14d8f0..76561fbcb39 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -560,6 +560,19 @@ end @test isapprox(lift, -6.501138753497174e-15, atol = 1e-13) @test isapprox(drag, 2.588589856781827, atol = 1e-13) end + +@trixi_testset "elixir_euler_NACA0012airfoil_mach085.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_NACA0012airfoil_mach085.jl"), + l2=[6.054106141023162e-7, + 7.4947583731596855e-6, + 1.1682455936043817e-5, + 0.0007173721760332341], + linf=[0.002472349864955632, + 0.043496823453937614, + 0.043648671347888836, + 3.046213202855595], tspan=(0.0, 0.0001)) +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 9f1382caa62..423b13995fa 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -689,6 +689,20 @@ 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, + "elixir_navierstokes_NACA0012airfoil_mach08.jl"), + l2=[0.00018648657393597384, + 0.0005076712152849281, + 0.00038074587715240566, + 0.0021281773710793315], + linf=[0.5153387749819276, + 1.1993620992082363, + 0.9077214408394708, + 5.666071686983816], tspan=(0.0, 0.001), + initial_refinement_level=0) +end end # Clean up afterwards: delete Trixi.jl output directory From 4d689b1675e00226ca1291f0a8f468db98c6b8db Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Sat, 10 Feb 2024 23:21:07 +0530 Subject: [PATCH 15/53] Fix tolerance, fix CI, add comments --- .../elixir_euler_NACA0012airfoil_mach085.jl | 3 ++- ...ixir_navierstokes_NACA0012airfoil_mach08.jl | 12 +++++++++++- test/test_parabolic_2d.jl | 18 +++++++++--------- 3 files changed, 22 insertions(+), 11 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index b7568296602..24404904ddd 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -123,7 +123,8 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition_only_refine = true) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, - stepsize_callback, amr_callback) + stepsize_callback, amr_callback + ) ############################################################################### # run the simulation diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index c15133de4e0..56d8ed703cf 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -1,3 +1,13 @@ +# 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. + +# - Roy Charles Swanson, Stefan Langer (2016) +# Structured and Unstructured Grid Methods (2016) +# [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623) + + using Downloads: download using OrdinaryDiffEq using Trixi @@ -133,7 +143,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav ############################################################################### # run the simulation -time_int_tol = 1e-11 +time_int_tol = 1e-8 sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, ode_default_options()..., callback = callbacks) summary_callback() # print the timer summary diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 423b13995fa..c8e33c1084b 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -691,16 +691,16 @@ end end @trixi_testset "elixir_navierstokes_NACA0012airfoil_mach08.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_navierstokes_NACA0012airfoil_mach08.jl"), - l2=[0.00018648657393597384, - 0.0005076712152849281, - 0.00038074587715240566, - 0.0021281773710793315], - linf=[0.5153387749819276, - 1.1993620992082363, - 0.9077214408394708, - 5.666071686983816], tspan=(0.0, 0.001), + 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) end end From 42e4b2d3f7b5c9ed7ce447cc027b34f91cf003c3 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Sat, 10 Feb 2024 23:23:28 +0530 Subject: [PATCH 16/53] Run formatter --- .../p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl | 3 +-- .../elixir_navierstokes_NACA0012airfoil_mach08.jl | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 24404904ddd..b7568296602 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -123,8 +123,7 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition_only_refine = true) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, - stepsize_callback, amr_callback - ) + stepsize_callback, amr_callback) ############################################################################### # run the simulation diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 56d8ed703cf..63483ae8496 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -7,7 +7,6 @@ # Structured and Unstructured Grid Methods (2016) # [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623) - using Downloads: download using OrdinaryDiffEq using Trixi From 79d5eda633a19092773ac64df5293b8ff8c771ac Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 15 Feb 2024 09:47:06 +0530 Subject: [PATCH 17/53] Minor changes --- .../elixir_euler_NACA0012airfoil_mach085.jl | 3 +- .../elixir_euler_subsonic_cylinder.jl | 6 ++-- src/callbacks_step/analysis_surface_2d.jl | 29 +++++++--------- test/test_p4est_2d.jl | 33 ++++++++++++++----- 4 files changed, 41 insertions(+), 30 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index b7568296602..d2e11cce985 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -117,8 +117,9 @@ amr_controller = ControllerThreeLevel(semi, amr_indicator, med_level = 3, med_threshold = 0.05, max_level = 4, max_threshold = 0.1) +amr_interval = 100 amr_callback = AMRCallback(semi, amr_controller, - interval = 100, + interval = amr_interval, adapt_initial_condition = true, adapt_initial_condition_only_refine = true) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 8901d5148b1..a8f5e6c5952 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -9,7 +9,7 @@ using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient equations = CompressibleEulerEquations2D(1.4) -@inline function initial_condition_mach01_flow(x, t, +@inline function initial_condition_mach038_flow(x, t, equations::CompressibleEulerEquations2D) # set the freestream flow parameters rho_freestream = 1.4 @@ -21,7 +21,7 @@ equations = CompressibleEulerEquations2D(1.4) return prim2cons(prim, equations) end -initial_condition = initial_condition_mach01_flow +initial_condition = initial_condition_mach038_flow volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be used surface_flux = flux_lax_friedrichs @@ -63,7 +63,7 @@ mesh = P4estMesh(cells_per_dimension, mapping = mapping2cylinder, polydeg = 3, t, surface_flux_function, equations::CompressibleEulerEquations2D) - u_boundary = initial_condition_mach01_flow(x, t, equations) + u_boundary = initial_condition_mach038_flow(x, t, equations) return surface_flux_function(u_inner, u_boundary, normal_direction, equations) end diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index cb667ee1f9e..f65137a052f 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -27,14 +27,7 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Indices, Variable} end struct ForceState{RealT <: Real} - Ψ::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream - rhoinf::RealT - uinf::RealT - linf::RealT -end - -# TODO - This should be a struct in ForceState -struct FreeStreamVariables{RealT <: Real} + psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream rhoinf::RealT uinf::RealT linf::RealT @@ -49,29 +42,29 @@ struct DragCoefficient{RealT <: Real} end function LiftCoefficient(aoa, rhoinf, uinf, linf) - # Ψl is the normal unit vector to the freestream direction - Ψl = (-sin(aoa), cos(aoa)) - force_state = ForceState(Ψl, rhoinf, uinf, linf) + # psi_lift is the normal unit vector to the freestream direction + psi_lift = (-sin(aoa), cos(aoa)) + force_state = ForceState(psi_lift, rhoinf, uinf, linf) return LiftCoefficient(force_state) end function DragCoefficient(aoa, rhoinf, uinf, linf) - # Ψd is the unit vector parallel to the freestream direction - Ψd = (cos(aoa), sin(aoa)) - return DragCoefficient(ForceState(Ψd, rhoinf, uinf, linf)) + # psi_drag is the unit vector parallel to the freestream direction + psi_drag = (cos(aoa), sin(aoa)) + return DragCoefficient(ForceState(psi_drag, rhoinf, uinf, linf)) end function (lift_coefficient::LiftCoefficient)(u, normal_direction, equations) p = pressure(u, equations) - @unpack Ψ, rhoinf, uinf, linf = lift_coefficient.force_state - n = dot(normal_direction, Ψ) / norm(normal_direction) + @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state + n = dot(normal_direction, psi) / norm(normal_direction) return p * n / (0.5 * rhoinf * uinf^2 * linf) end function (drag_coefficient::DragCoefficient)(u, normal_direction, equations) p = pressure(u, equations) - @unpack Ψ, rhoinf, uinf, linf = drag_coefficient.force_state - n = dot(normal_direction, Ψ) / norm(normal_direction) + @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state + n = dot(normal_direction, psi) / norm(normal_direction) return p * n / (0.5 * rhoinf * uinf^2 * linf) end diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 76561fbcb39..07092c6d531 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -561,17 +561,34 @@ end @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=[6.054106141023162e-7, - 7.4947583731596855e-6, - 1.1682455936043817e-5, - 0.0007173721760332341], - linf=[0.002472349864955632, - 0.043496823453937614, - 0.043648671347888836, - 3.046213202855595], tspan=(0.0, 0.0001)) + l2=[5.756878791758828e-7, + 6.796904726549772e-6, + 1.1127326665302503e-5, + 0.0006841659623491992], + linf=[0.0023954158177901686, + 0.04290641165627283, + 0.04250015332478646, + 2.9640108432161285], + base_level = 0, med_level = 1, max_level = 1, + amr_interval = 1, + tspan=(0.0, 0.0001)) + + 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.029482383224495173, atol = 1e-13) + @test isapprox(drag, 0.14586550052248398, atol = 1e-13) end end From 809c2506292274c0b7e158b909bf7a9452c3ecab Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 15 Feb 2024 09:49:45 +0530 Subject: [PATCH 18/53] Run formatter --- .../p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 2 +- src/callbacks_step/analysis_surface_2d.jl | 2 +- test/test_p4est_2d.jl | 12 ++++++------ test/test_parabolic_2d.jl | 6 +++--- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index a8f5e6c5952..bac64ac5baf 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -10,7 +10,7 @@ using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient equations = CompressibleEulerEquations2D(1.4) @inline function initial_condition_mach038_flow(x, t, - equations::CompressibleEulerEquations2D) + equations::CompressibleEulerEquations2D) # set the freestream flow parameters rho_freestream = 1.4 v1 = 0.38 diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index f65137a052f..971d682fbb3 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -128,4 +128,4 @@ function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, <:DragCoefficient{<:Any}}) "CD" end -end # muladd \ No newline at end of file +end # muladd diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 07092c6d531..8cb90cd2d0e 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -570,12 +570,12 @@ end 1.1127326665302503e-5, 0.0006841659623491992], linf=[0.0023954158177901686, - 0.04290641165627283, - 0.04250015332478646, - 2.9640108432161285], - base_level = 0, med_level = 1, max_level = 1, - amr_interval = 1, - tspan=(0.0, 0.0001)) + 0.04290641165627283, + 0.04250015332478646, + 2.9640108432161285], + base_level=0, med_level=1, max_level=1, + amr_interval=1, + tspan=(0.0, 0.0001)) u_ode = copy(sol.u[end]) du_ode = zero(u_ode) # Just a placeholder in this case diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index c8e33c1084b..3f4f7183a4e 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -698,9 +698,9 @@ end 0.00038074588984354107, 0.002128177239782089], linf=[0.5153387072802718, - 1.199362305026636, - 0.9077214424040279, - 5.666071182328691], tspan=(0.0, 0.001), + 1.199362305026636, + 0.9077214424040279, + 5.666071182328691], tspan=(0.0, 0.001), initial_refinement_level=0) end end From 69b0f67fe0ba73e0a346674907ebf3e8b2a72d8a Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 15 Feb 2024 22:26:54 +0530 Subject: [PATCH 19/53] Fix type instability and need for a vector --- .../elixir_euler_NACA0012airfoil_mach085.jl | 2 + .../elixir_euler_subsonic_cylinder.jl | 4 +- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 4 +- src/callbacks_step/analysis_surface_2d.jl | 47 ++++++++++++------- 4 files changed, 36 insertions(+), 21 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index d2e11cce985..1f9eb90f3e6 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -86,6 +86,8 @@ analysis_interval = 2000 linf = 1.0 # Length of airfoil +# The boundary on which we want to compute the forces is the one where boundary conditions +# are slip wall. drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, DragCoefficient(aoa(), rho_inf(), U_inf(equations), linf)) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index bac64ac5baf..92c188b022e 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -55,8 +55,8 @@ cells_per_dimension = (64, 64) mesh = P4estMesh(cells_per_dimension, mapping = mapping2cylinder, polydeg = 3, periodicity = (false, true)) -# The boundary of the outer cylinder is constant but subsonic, so we cannot compute the -# boundary flux for the external information alone. Thus, we use the numerical flux to distinguish +# The boundary conditions of 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, diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 63483ae8496..d101f6c1a98 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -58,8 +58,8 @@ mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/3396 mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 1) -# The boundary of the outer cylinder is constant but subsonic, so we cannot compute the -# boundary flux for the external information alone. Thus, we use the numerical flux to distinguish +# 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, diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index 971d682fbb3..afcda42c025 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -1,14 +1,26 @@ +# 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 -# TODO - Add explanation for what this file contains -struct AnalysisSurfaceIntegral{Semidiscretization, Indices, Variable} - semi::Semidiscretization - indices::Indices +# The boundary_index is chosen so that +# semi.boundary_conditions.boundary_indices[boundary_index] +# gives the solution point indices on which the function `variable` will compute +# the quantity of interest. The `variable` can, e.g., integrate over all indices +# to compute the coefficient of lift. But it can also be used to loop over all indices +# to save coefficient of pressure versus surface points in a data file +# The struct contains an inner constructor which helps the user choose those indices +# by specifying the `boundary_condition_type` of the indices. The user needs to make +# sure that they choose the `boundary_condition_type` so that it is only applying to the +# parts of boundary that are of interest +struct AnalysisSurfaceIntegral{Semidiscretization, Variable} + semi::Semidiscretization # Semidiscretization of PDE used by the solver + boundary_index::Int # Index in boundary_condition_indices where quantitiy of interest is computed variable::Variable # Quantity of interest, like lift or drag function AnalysisSurfaceIntegral(semi, boundary_condition_type, variable) # The bc list as ordered in digest_boundary_conditions @@ -17,12 +29,12 @@ struct AnalysisSurfaceIntegral{Semidiscretization, Indices, Variable} # The set of all indices that gives the bc where the surface integral is to be computed index = sort(findall(x -> x == boundary_condition_type, ordered_bc)) - # Put the bc in function form as they might change under AMR - indices = semi -> Vector([semi.boundary_conditions.boundary_indices[i] - for i in index])[1] # TODO - Should not need Vector and the [1] + # digest_boundary_conditions clubs all indices with same boundary conditions into + # one. This is just checking that it is indeed the case for the next step. + @assert length(index) == 1 - return new{typeof(semi), typeof(indices), typeof(variable)}(semi, indices, - variable) + return new{typeof(semi), typeof(variable)}(semi, index[1], + variable) end end @@ -74,14 +86,15 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, @unpack boundaries = cache @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements @unpack weights = dg.basis - @unpack semi, indices, variable = surface_variable - indices_ = indices(semi) + @unpack semi, boundary_index, variable = surface_variable + + indices = semi.boundary_conditions.boundary_indices[boundary_index] surface_integral = zero(eltype(u)) index_range = eachnode(dg) - for local_index in eachindex(indices_) + for local_index in eachindex(indices) # Use the local index to get the global boundary index from the pre-sorted list - boundary = indices_[local_index] + boundary = indices[local_index] # Get information on the adjacent element, compute the surface fluxes, # and store them @@ -112,19 +125,19 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, return surface_integral end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:LiftCoefficient{<:Any}}) "CL" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:LiftCoefficient{<:Any}}) "CL" end -function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:Any, +function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, <:DragCoefficient{<:Any}}) "CD" end -function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:Any, +function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, <:DragCoefficient{<:Any}}) "CD" end From feab5b34a35b0eeffcc2e47a6fcc798e2fc879ff Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 15 Feb 2024 22:43:52 +0530 Subject: [PATCH 20/53] Fix typo! --- 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 afcda42c025..8de396fd05c 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -20,7 +20,7 @@ # parts of boundary that are of interest struct AnalysisSurfaceIntegral{Semidiscretization, Variable} semi::Semidiscretization # Semidiscretization of PDE used by the solver - boundary_index::Int # Index in boundary_condition_indices where quantitiy of interest is computed + boundary_index::Int # Index in boundary_condition_indices where quantity of interest is computed variable::Variable # Quantity of interest, like lift or drag function AnalysisSurfaceIntegral(semi, boundary_condition_type, variable) # The bc list as ordered in digest_boundary_conditions From 007bee60b4d456ee99b812c5893614e5f3364f30 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Fri, 16 Feb 2024 08:10:08 +0530 Subject: [PATCH 21/53] Lower CFL to pass CI? --- test/test_p4est_2d.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 8cb90cd2d0e..892cbeefb53 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -565,16 +565,16 @@ end @trixi_testset "elixir_euler_NACA0012airfoil_mach085.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), - l2=[5.756878791758828e-7, - 6.796904726549772e-6, - 1.1127326665302503e-5, - 0.0006841659623491992], - linf=[0.0023954158177901686, - 0.04290641165627283, - 0.04250015332478646, - 2.9640108432161285], + l2=[5.763969153227137e-7, + 6.7983731359493415e-6, + 1.1136096550405798e-5, + 0.000684989462539266], + linf=[0.002411375295350659, + 0.04305328756849549, + 0.04272712024785543, + 2.981598654532575], base_level=0, med_level=1, max_level=1, - amr_interval=1, + amr_interval=1, cfl=0.1, tspan=(0.0, 0.0001)) u_ode = copy(sol.u[end]) @@ -587,8 +587,8 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 0.029482383224495173, atol = 1e-13) - @test isapprox(drag, 0.14586550052248398, atol = 1e-13) + @test isapprox(lift, 0.029723886821464478, atol = 1e-13) + @test isapprox(drag, 0.1463432877615894, atol = 1e-13) end end From 1b85765350fdaa9682a79e3a161046ede10d1673 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Fri, 16 Feb 2024 11:17:32 +0530 Subject: [PATCH 22/53] Reduce amr_interval to pass CI? --- test/test_p4est_2d.jl | 34 ++++++++++++++++++++++------------ 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 892cbeefb53..3e31de66f33 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -565,17 +565,23 @@ end @trixi_testset "elixir_euler_NACA0012airfoil_mach085.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), - l2=[5.763969153227137e-7, - 6.7983731359493415e-6, - 1.1136096550405798e-5, - 0.000684989462539266], - linf=[0.002411375295350659, - 0.04305328756849549, - 0.04272712024785543, - 2.981598654532575], + l2=[ + 5.59290009961641e-7, + 6.709050529319647e-6, + 1.0832036354579161e-5, + 0.00066312990043519, + ], + linf=[ + 0.0020566846007074946, + 0.038004977184127715, + 0.03680677713072339, + 2.5456859522627036, + ], base_level=0, med_level=1, max_level=1, - amr_interval=1, cfl=0.1, - tspan=(0.0, 0.0001)) + amr_interval=5, cfl=0.25, + tspan=(0.0, 0.0001), + adapt_initial_condition=false, + adapt_initial_condition_only_refine=false) u_ode = copy(sol.u[end]) du_ode = zero(u_ode) # Just a placeholder in this case @@ -587,8 +593,12 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 0.029723886821464478, atol = 1e-13) - @test isapprox(drag, 0.1463432877615894, atol = 1e-13) + @show drag, lift + @show analysis_callback(sol).l2 + @show analysis_callback(sol).linf + + @test isapprox(lift, 0.030280394867287057, atol = 1e-13) + @test isapprox(drag, 0.13085519947527666, atol = 1e-13) end end From dc8365ce154695b1af15577aeba3e21084c50976 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Fri, 16 Feb 2024 12:10:18 +0530 Subject: [PATCH 23/53] Maybe positivity limiter will fix CI? --- .../p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl | 5 ++++- test/test_p4est_2d.jl | 4 ---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 1f9eb90f3e6..fc298378fc6 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -128,9 +128,12 @@ amr_callback = AMRCallback(semi, amr_controller, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback, amr_callback) +stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6), + variables = (Trixi.density, pressure)) + ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false), 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/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 3e31de66f33..9257fdd2200 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -593,10 +593,6 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @show drag, lift - @show analysis_callback(sol).l2 - @show analysis_callback(sol).linf - @test isapprox(lift, 0.030280394867287057, atol = 1e-13) @test isapprox(drag, 0.13085519947527666, atol = 1e-13) end From 62a57fa15347d29f796db73f57e4039e44be46bd Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Fri, 16 Feb 2024 13:16:24 +0530 Subject: [PATCH 24/53] Remove AMR from CI --- .../elixir_euler_NACA0012airfoil_mach085.jl | 5 +---- test/test_p4est_2d.jl | 17 ++++++----------- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index fc298378fc6..1f9eb90f3e6 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -128,12 +128,9 @@ amr_callback = AMRCallback(semi, amr_controller, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback, amr_callback) -stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-6), - variables = (Trixi.density, pressure)) - ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition = false), +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), 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/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 9257fdd2200..82d8c22c262 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -566,19 +566,14 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), l2=[ - 5.59290009961641e-7, - 6.709050529319647e-6, - 1.0832036354579161e-5, - 0.00066312990043519, + 5.376797508531838e-7, 6.419368395509845e-6, + 1.0333204566846605e-5, 0.0006354281448373737, ], linf=[ - 0.0020566846007074946, - 0.038004977184127715, - 0.03680677713072339, - 2.5456859522627036, + 0.00162885239831281, 0.028512705456618496, + 0.029917405696784197, 1.951032056018755, ], base_level=0, med_level=1, max_level=1, - amr_interval=5, cfl=0.25, tspan=(0.0, 0.0001), adapt_initial_condition=false, adapt_initial_condition_only_refine=false) @@ -593,8 +588,8 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 0.030280394867287057, atol = 1e-13) - @test isapprox(drag, 0.13085519947527666, atol = 1e-13) + @test isapprox(lift, 0.10909749380301453, atol = 1e-13) + @test isapprox(drag, 0.026361875495671695, atol = 1e-13) end end From 2afbc38aff3d3b21b28b178660b15fc2206f7362 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Fri, 16 Feb 2024 13:51:37 +0530 Subject: [PATCH 25/53] That CI fail was on me! --- test/test_p4est_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 82d8c22c262..ef203a395b3 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -588,8 +588,8 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 0.10909749380301453, atol = 1e-13) - @test isapprox(drag, 0.026361875495671695, atol = 1e-13) + @test isapprox(lift, 0.026361875495671695, atol = 1e-13) + @test isapprox(drag, 0.10909749380301453, atol = 1e-13) end end From 881810c73a367988ed10adf97e15ea10c54d57ea Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 19 Mar 2024 18:28:18 +0100 Subject: [PATCH 26/53] Test for write permit --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 1 - .../p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 92c188b022e..3a340555dbe 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -27,7 +27,6 @@ volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be u surface_flux = flux_lax_friedrichs polydeg = 3 -basis = LobattoLegendreBasis(polydeg) solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index d101f6c1a98..aeb311fbf6b 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -50,7 +50,6 @@ initial_condition = initial_condition_mach08_flow surface_flux = flux_lax_friedrichs polydeg = 3 -basis = LobattoLegendreBasis(polydeg) solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux) mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp", From d914526786fd5f74333cc8318359580dd42fe80f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 19 Mar 2024 19:02:27 +0100 Subject: [PATCH 27/53] forces via names --- .../elixir_euler_NACA0012airfoil_mach085.jl | 9 +++--- src/callbacks_step/analysis_surface_2d.jl | 30 +++++++++++-------- .../sort_boundary_conditions.jl | 14 ++++++++- 3 files changed, 34 insertions(+), 19 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 1f9eb90f3e6..d6a38e3827a 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -86,18 +86,17 @@ analysis_interval = 2000 linf = 1.0 # Length of airfoil -# The boundary on which we want to compute the forces is the one where boundary conditions -# are slip wall. -drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, +force_boundary_names = [:AirfoilBottom, :AirfoilTop] +drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, DragCoefficient(aoa(), rho_inf(), U_inf(equations), linf)) -lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, +lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, LiftCoefficient(aoa(), rho_inf(), U_inf(equations), linf)) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - output_directory = "analysis_results", + output_directory = "out", save_analysis = true, analysis_integrals = (drag_coefficient, lift_coefficient)) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index 8de396fd05c..715f7cb95f3 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -20,20 +20,26 @@ # parts of boundary that are of interest struct AnalysisSurfaceIntegral{Semidiscretization, Variable} semi::Semidiscretization # Semidiscretization of PDE used by the solver - boundary_index::Int # Index in boundary_condition_indices where quantity of interest is computed + 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_condition_type, variable) - # The bc list as ordered in digest_boundary_conditions - ordered_bc = semi.boundary_conditions.boundary_condition_types - # The set of all indices that gives the bc where the surface integral is to be computed - index = sort(findall(x -> x == boundary_condition_type, ordered_bc)) + function AnalysisSurfaceIntegral(semi, boundary_symbol, variable) + @unpack boundary_symbol_indices = semi.boundary_conditions + indices = boundary_symbol_indices[boundary_symbol] - # digest_boundary_conditions clubs all indices with same boundary conditions into - # one. This is just checking that it is indeed the case for the next step. - @assert length(index) == 1 + return new{typeof(semi), typeof(variable)}(semi, indices, + variable) + end - return new{typeof(semi), typeof(variable)}(semi, index[1], + 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 @@ -86,9 +92,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, @unpack boundaries = cache @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements @unpack weights = dg.basis - @unpack semi, boundary_index, variable = surface_variable - - indices = semi.boundary_conditions.boundary_indices[boundary_index] + @unpack semi, indices, variable = surface_variable surface_integral = zero(eltype(u)) index_range = eachnode(dg) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index b5388cadc8b..5c153af2640 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -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 + for (symbol, _) in boundary_dictionary + indices = findall(x -> x === symbol, cache.boundaries.name) + # Store the indices in `boundary_symbol_indices` + boundary_types_container.boundary_symbol_indices[symbol] = sort!(indices) + end + return boundary_types_container end end # @muladd From a198ee0256e59f164e87fa24a112a0fae4b7b38b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 19 Mar 2024 19:10:29 +0100 Subject: [PATCH 28/53] Boundary Names for other examples --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 6 +++--- .../elixir_navierstokes_NACA0012airfoil_mach08.jl | 7 ++++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 3a340555dbe..4300bc76d5e 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -91,14 +91,14 @@ rho_inf = 1.4 U_inf = 0.38 linf = 1.0 # Diameter of circle -drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, +drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, DragCoefficient(aoa, rho_inf, U_inf, linf)) -lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, +lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, LiftCoefficient(aoa, rho_inf, U_inf, linf)) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - output_directory = "analysis_results", + output_directory = "out", save_analysis = true, analysis_integrals = (drag_coefficient, lift_coefficient)) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index aeb311fbf6b..39a19b0c4e6 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -115,16 +115,17 @@ summary_callback = SummaryCallback() analysis_interval = 2000 -drag_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, +force_boundary_names = [:AirfoilBottom, :AirfoilTop] +drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, DragCoefficient(sw_aoa(), sw_rho_inf(), sw_U_inf(equations), sw_linf())) -lift_coefficient = AnalysisSurfaceIntegral(semi, boundary_condition_slip_wall, +lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, LiftCoefficient(sw_aoa(), sw_rho_inf(), sw_U_inf(equations), sw_linf())) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - output_directory = "analysis_results", + output_directory = "out", save_analysis = true, analysis_integrals = (drag_coefficient, lift_coefficient)) From ca9c0e8426ff014bf5c54f11ef81b6683cfaf863 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 20 Mar 2024 01:06:06 +0530 Subject: [PATCH 29/53] Shift to SSPRK54 to pass CI? --- .../elixir_euler_NACA0012airfoil_mach085.jl | 2 +- test/test_p4est_2d.jl | 14 ++++++++------ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index d6a38e3827a..d33f49b06e3 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -129,7 +129,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +sol = solve(ode, SSPRK54(), 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/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index ef203a395b3..bf65a2490e8 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -566,12 +566,12 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), l2=[ - 5.376797508531838e-7, 6.419368395509845e-6, - 1.0333204566846605e-5, 0.0006354281448373737, + 5.37640294704164e-7, 6.418954697992308e-6, + 1.0332595763357066e-5, 0.0006353813399881882, ], linf=[ - 0.00162885239831281, 0.028512705456618496, - 0.029917405696784197, 1.951032056018755, + 0.0016284523634016628, 0.028505821812697323, + 0.029918806073518, 1.9505653217814127, ], base_level=0, med_level=1, max_level=1, tspan=(0.0, 0.0001), @@ -588,8 +588,10 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 0.026361875495671695, atol = 1e-13) - @test isapprox(drag, 0.10909749380301453, atol = 1e-13) + @show analysis_callback(sol) + + @test isapprox(lift, 0.026397498239816828, atol = 1e-13) + @test isapprox(drag, 0.10908833968266139, atol = 1e-13) end end From b2bec6f77bc5b4af69902ea9610409a5c0d6e560 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 19 Mar 2024 23:36:48 +0100 Subject: [PATCH 30/53] docstrings --- src/callbacks_step/analysis_surface_2d.jl | 32 ++++++++----------- .../sort_boundary_conditions.jl | 4 +-- 2 files changed, 16 insertions(+), 20 deletions(-) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index 715f7cb95f3..f0a3bf6b35a 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -8,18 +8,17 @@ @muladd begin #! format: noindent -# The boundary_index is chosen so that -# semi.boundary_conditions.boundary_indices[boundary_index] -# gives the solution point indices on which the function `variable` will compute -# the quantity of interest. The `variable` can, e.g., integrate over all indices -# to compute the coefficient of lift. But it can also be used to loop over all indices -# to save coefficient of pressure versus surface points in a data file -# The struct contains an inner constructor which helps the user choose those indices -# by specifying the `boundary_condition_type` of the indices. The user needs to make -# sure that they choose the `boundary_condition_type` so that it is only applying to the -# parts of boundary that are of interest +""" + 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`. + For instance, this can be used to compute the lift or drag coefficient of e.g. an airfoil in 2D. +""" struct AnalysisSurfaceIntegral{Semidiscretization, Variable} - semi::Semidiscretization # Semidiscretization of PDE used by the solver + 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 @@ -92,16 +91,12 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, @unpack boundaries = cache @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements @unpack weights = dg.basis - @unpack semi, indices, variable = surface_variable + + @unpack indices, variable = surface_variable surface_integral = zero(eltype(u)) index_range = eachnode(dg) - for local_index in eachindex(indices) - # Use the local index to get the global boundary index from the pre-sorted list - boundary = indices[local_index] - - # Get information on the adjacent element, compute the surface fluxes, - # and store them + for boundary in indices element = boundaries.neighbor_ids[boundary] node_indices = boundaries.node_indices[boundary] direction = indices2direction(node_indices) @@ -120,6 +115,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, # L2 norm of normal direction is the surface element dS = weights[node_index] * norm(normal_direction) + # Integral over whole boundary surface surface_integral += variable(u_node, normal_direction, equations) * dS i_node += i_node_step diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 5c153af2640..b648cf844e3 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. """ From 9fcbb04a903cea36bdd9a2c43824f9420976e0ea Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 19 Mar 2024 23:49:27 +0100 Subject: [PATCH 31/53] alloc tests --- test/test_p4est_2d.jl | 19 ++++++++++++++++++- test/test_parabolic_2d.jl | 8 ++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index bf65a2490e8..f5ccf2ffd38 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -546,7 +546,15 @@ end 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 @@ -578,6 +586,15 @@ end adapt_initial_condition=false, adapt_initial_condition_only_refine=false) + # 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 diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 3f4f7183a4e..9ba09a9268e 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -702,6 +702,14 @@ end 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 From 90ef9a810f15a67ae891f56c268dfec8de891396 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 19 Mar 2024 23:53:02 +0100 Subject: [PATCH 32/53] fmt --- test/test_p4est_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index f5ccf2ffd38..17a45870d68 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -554,7 +554,7 @@ 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 From e3f6cbc21906b992b6c4136338c897163120526c Mon Sep 17 00:00:00 2001 From: Arpit Babbar Date: Wed, 20 Mar 2024 07:05:45 +0530 Subject: [PATCH 33/53] Apply suggestions from code review Co-authored-by: Daniel Doehring Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- .../p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl | 2 +- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 2 +- test/test_p4est_2d.jl | 6 ++++-- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index d33f49b06e3..668d8dad60c 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -12,7 +12,7 @@ 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 mach_inf() = 0.85 -aoa() = pi / 180.0 +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) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 4300bc76d5e..052ef4173b9 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -51,7 +51,7 @@ cells_per_dimension = (64, 64) # 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 = 3, +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 diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 17a45870d68..bf3d4268a4d 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -584,7 +584,10 @@ end base_level=0, med_level=1, max_level=1, tspan=(0.0, 0.0001), adapt_initial_condition=false, - adapt_initial_condition_only_refine=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) @@ -605,7 +608,6 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @show analysis_callback(sol) @test isapprox(lift, 0.026397498239816828, atol = 1e-13) @test isapprox(drag, 0.10908833968266139, atol = 1e-13) From d1a3507774b3ee419771a5a2b4ffcec4a9775a2f Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 20 Mar 2024 07:09:29 +0530 Subject: [PATCH 34/53] Format --- test/test_p4est_2d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index bf3d4268a4d..cd4e6f7759b 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -608,7 +608,6 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 0.026397498239816828, atol = 1e-13) @test isapprox(drag, 0.10908833968266139, atol = 1e-13) end From 0fef9ad795726fcb780105473ed46e6916cb415e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 24 Mar 2024 15:50:23 +0100 Subject: [PATCH 35/53] non-allocating BC --- .../elixir_navierstokes_NACA0012airfoil_mach08.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 39a19b0c4e6..8428429aafa 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -84,8 +84,12 @@ heat_airfoil = Adiabatic((x, t, equations) -> 0.0) boundary_conditions_airfoil = BoundaryConditionNavierStokesWall(velocity_airfoil, heat_airfoil) -velocity_bc_square = NoSlip((x, t, equations) -> initial_condition_mach08_flow(x, t, - equations)[2:3]) +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) @@ -127,6 +131,7 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", save_analysis = true, + analysis_errors = Symbol[], analysis_integrals = (drag_coefficient, lift_coefficient)) @@ -143,6 +148,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav # run the simulation time_int_tol = 1e-8 -sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, +sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = time_int_tol, reltol = time_int_tol, ode_default_options()..., callback = callbacks) -summary_callback() # print the timer summary +summary_callback() # print the timer summary \ No newline at end of file From 89f055945c03635d5a104522952d369dc2ccae24 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 27 Mar 2024 09:41:21 +0100 Subject: [PATCH 36/53] comments --- src/callbacks_step/analysis_surface_2d.jl | 11 +++++++---- src/solvers/dgsem_p4est/dg.jl | 3 ++- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index f0a3bf6b35a..52d12cfaad6 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -61,12 +61,11 @@ end function LiftCoefficient(aoa, rhoinf, uinf, linf) # psi_lift is the normal unit vector to the freestream direction psi_lift = (-sin(aoa), cos(aoa)) - force_state = ForceState(psi_lift, rhoinf, uinf, linf) - return LiftCoefficient(force_state) + return LiftCoefficient(ForceState(psi_lift, rhoinf, uinf, linf)) end function DragCoefficient(aoa, rhoinf, uinf, linf) - # psi_drag is the unit vector parallel to the freestream direction + # `psi_drag` is the unit vector in direction of the freestream psi_drag = (cos(aoa), sin(aoa)) return DragCoefficient(ForceState(psi_drag, rhoinf, uinf, linf)) end @@ -74,6 +73,7 @@ end function (lift_coefficient::LiftCoefficient)(u, normal_direction, equations) p = pressure(u, equations) @unpack psi, rhoinf, uinf, linf = 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 * linf) end @@ -81,6 +81,7 @@ end function (drag_coefficient::DragCoefficient)(u, normal_direction, equations) p = pressure(u, equations) @unpack psi, rhoinf, uinf, linf = 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 * linf) end @@ -106,9 +107,11 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, i_node = i_node_start j_node = j_node_start - for node_index in eachnode(dg) + 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) 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 From e877dbc2c12627f31e4ff53663d125ecf1cd812e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 27 Mar 2024 10:15:22 +0100 Subject: [PATCH 37/53] comments --- src/callbacks_step/analysis_surface_2d.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/analysis_surface_2d.jl b/src/callbacks_step/analysis_surface_2d.jl index 52d12cfaad6..d346b7e351f 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_2d.jl @@ -59,13 +59,17 @@ struct DragCoefficient{RealT <: Real} end function LiftCoefficient(aoa, rhoinf, uinf, linf) - # psi_lift is the normal unit vector to the freestream direction + # 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 + # value, but with the opposite sign. psi_lift = (-sin(aoa), cos(aoa)) return LiftCoefficient(ForceState(psi_lift, rhoinf, uinf, linf)) end function DragCoefficient(aoa, rhoinf, uinf, linf) - # `psi_drag` is the unit vector in direction of the freestream + # `psi_drag` is the unit vector in direction of the freestream. psi_drag = (cos(aoa), sin(aoa)) return DragCoefficient(ForceState(psi_drag, rhoinf, uinf, linf)) end @@ -116,7 +120,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, i_node, j_node, element) - # L2 norm of normal direction is the surface element + # L2 norm of normal direction (contravariant_vector) is the surface element dS = weights[node_index] * norm(normal_direction) # Integral over whole boundary surface surface_integral += variable(u_node, normal_direction, equations) * dS From 71a5509c174d597bf85793cbc9653f639dd1c2ec Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 27 Mar 2024 10:16:00 +0100 Subject: [PATCH 38/53] fmt --- .../elixir_navierstokes_NACA0012airfoil_mach08.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 8428429aafa..1c0a38ae15d 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -85,10 +85,11 @@ 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]) + 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)) +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, @@ -148,6 +149,7 @@ 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 = time_int_tol, + reltol = time_int_tol, ode_default_options()..., callback = callbacks) -summary_callback() # print the timer summary \ No newline at end of file +summary_callback() # print the timer summary From fd50740bbf3b400647c7e683da5e8522f7653950 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 27 Mar 2024 10:16:10 +0100 Subject: [PATCH 39/53] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index b648cf844e3..f4ad93adbd4 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -102,7 +102,7 @@ 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 + # 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` From 23306b65230c318a3ca08891317edee335377b75 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 27 Mar 2024 10:16:22 +0100 Subject: [PATCH 40/53] Update src/solvers/dgsem_unstructured/sort_boundary_conditions.jl --- src/solvers/dgsem_unstructured/sort_boundary_conditions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index f4ad93adbd4..2c2c6876d70 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -105,7 +105,7 @@ function initialize!(boundary_types_container::UnstructuredSortedBoundaryTypes{N # 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` + # Store the indices in `boundary_symbol_indices` dictionary boundary_types_container.boundary_symbol_indices[symbol] = sort!(indices) end From e0977a924a32e3831a409487769e899ce6e584ad Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 27 Mar 2024 10:17:45 +0100 Subject: [PATCH 41/53] comments --- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 24 ++++++++++++------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 1c0a38ae15d..8679ebac2d3 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -1,12 +1,3 @@ -# 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. - -# - Roy Charles Swanson, Stefan Langer (2016) -# Structured and Unstructured Grid Methods (2016) -# [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623) - using Downloads: download using OrdinaryDiffEq using Trixi @@ -16,6 +7,21 @@ using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient ############################################################################### # semidiscretization of the compressible Euler equations +# 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. + +# - 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 +# [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 From d0a1e3ecbc1335f58dcf3f1446f98c4b2abce1fa Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 27 Mar 2024 16:42:38 +0100 Subject: [PATCH 42/53] make ready for review --- NEWS.md | 2 +- .../elixir_euler_NACA0012airfoil_mach085.jl | 10 ++- .../elixir_euler_subsonic_cylinder.jl | 8 +-- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 12 ++-- src/Trixi.jl | 3 +- src/callbacks_step/analysis.jl | 2 +- ..._2d.jl => analysis_surface_integral_2d.jl} | 70 ++++++++++++++----- 7 files changed, 71 insertions(+), 36 deletions(-) rename src/callbacks_step/{analysis_surface_2d.jl => analysis_surface_integral_2d.jl} (71%) diff --git a/NEWS.md b/NEWS.md index 5b08d51ab89..5bd93e437a4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,7 @@ for human readability. #### Added - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. - +- Pressure lift and drag coefficients can now be computed for 2D `P4estMesh`. ## 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 index 668d8dad60c..59c46977267 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -2,8 +2,6 @@ using Downloads: download using OrdinaryDiffEq using Trixi -using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient - ############################################################################### # semidiscretization of the compressible Euler equations @@ -88,12 +86,12 @@ linf = 1.0 # Length of airfoil force_boundary_names = [:AirfoilBottom, :AirfoilTop] drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, - DragCoefficient(aoa(), rho_inf(), - U_inf(equations), linf)) + DragCoefficientPressure(aoa(), rho_inf(), + U_inf(equations), linf)) lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, - LiftCoefficient(aoa(), rho_inf(), - U_inf(equations), linf)) + LiftCoefficientPressure(aoa(), rho_inf(), + U_inf(equations), linf)) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index 052ef4173b9..df8782cef82 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -2,8 +2,6 @@ using Downloads: download using OrdinaryDiffEq using Trixi -using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient - ############################################################################### # semidiscretization of the compressible Euler equations @@ -92,10 +90,12 @@ U_inf = 0.38 linf = 1.0 # Diameter of circle drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, - DragCoefficient(aoa, rho_inf, U_inf, linf)) + DragCoefficientPressure(aoa, rho_inf, U_inf, + linf)) lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, - LiftCoefficient(aoa, rho_inf, U_inf, linf)) + LiftCoefficientPressure(aoa, rho_inf, U_inf, + linf)) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 8679ebac2d3..80b5a840c32 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -2,8 +2,6 @@ using Downloads: download using OrdinaryDiffEq using Trixi -using Trixi: AnalysisSurfaceIntegral, DragCoefficient, LiftCoefficient - ############################################################################### # semidiscretization of the compressible Euler equations @@ -128,12 +126,14 @@ analysis_interval = 2000 force_boundary_names = [:AirfoilBottom, :AirfoilTop] drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, - DragCoefficient(sw_aoa(), sw_rho_inf(), - sw_U_inf(equations), sw_linf())) + DragCoefficientPressure(sw_aoa(), sw_rho_inf(), + sw_U_inf(equations), + sw_linf())) lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, - LiftCoefficient(sw_aoa(), sw_rho_inf(), - sw_U_inf(equations), sw_linf())) + LiftCoefficientPressure(sw_aoa(), sw_rho_inf(), + sw_U_inf(equations), + sw_linf())) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", diff --git a/src/Trixi.jl b/src/Trixi.jl index da7359999c5..38fa99504a1 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -262,7 +262,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 2c3351272fa..62048853b53 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -691,7 +691,7 @@ end # @muladd # specialized implementations specific to some solvers include("analysis_dg1d.jl") include("analysis_dg2d.jl") -include("analysis_surface_2d.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_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl similarity index 71% rename from src/callbacks_step/analysis_surface_2d.jl rename to src/callbacks_step/analysis_surface_integral_2d.jl index d346b7e351f..441413ad04a 100644 --- a/src/callbacks_step/analysis_surface_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -15,7 +15,8 @@ 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`. - For instance, this can be used to compute the lift or drag coefficient of e.g. an airfoil in 2D. + For instance, this can be used to compute the lift [`LiftCoefficientPressure`](@ref) or + drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil in 2D. """ struct AnalysisSurfaceIntegral{Semidiscretization, Variable} semi::Semidiscretization # passed in to retrieve boundary condition information @@ -50,31 +51,65 @@ struct ForceState{RealT <: Real} linf::RealT end -struct LiftCoefficient{RealT <: Real} +struct LiftCoefficientPressure{RealT <: Real} force_state::ForceState{RealT} end -struct DragCoefficient{RealT <: Real} +struct DragCoefficientPressure{RealT <: Real} force_state::ForceState{RealT} end -function LiftCoefficient(aoa, rhoinf, uinf, linf) +""" + LiftCoefficientPressure(aoa, rhoinf, uinf, linf) + +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 conjuction 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 +- `linf::Real`: Reference length of geoemtry (e.g. airfoil chord length) +""" +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 # value, but with the opposite sign. psi_lift = (-sin(aoa), cos(aoa)) - return LiftCoefficient(ForceState(psi_lift, rhoinf, uinf, linf)) + return LiftCoefficientPressure(ForceState(psi_lift, rhoinf, uinf, linf)) end -function DragCoefficient(aoa, rhoinf, uinf, linf) +""" + DragCoefficientPressure(aoa, rhoinf, uinf, linf) + +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 conjuction 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 +- `linf::Real`: Reference length of geoemtry (e.g. airfoil chord length) +""" +function DragCoefficientPressure(aoa, rhoinf, uinf, linf) # `psi_drag` is the unit vector in direction of the freestream. psi_drag = (cos(aoa), sin(aoa)) - return DragCoefficient(ForceState(psi_drag, rhoinf, uinf, linf)) + return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf)) end -function (lift_coefficient::LiftCoefficient)(u, normal_direction, equations) +function (lift_coefficient::LiftCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) @unpack psi, rhoinf, uinf, linf = lift_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector @@ -82,7 +117,7 @@ function (lift_coefficient::LiftCoefficient)(u, normal_direction, equations) return p * n / (0.5 * rhoinf * uinf^2 * linf) end -function (drag_coefficient::DragCoefficient)(u, normal_direction, equations) +function (drag_coefficient::DragCoefficientPressure)(u, normal_direction, equations) p = pressure(u, equations) @unpack psi, rhoinf, uinf, linf = drag_coefficient.force_state # Normalize as `normal_direction` is not necessarily a unit vector @@ -133,19 +168,20 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, end function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficient{<:Any}}) - "CL" + <:LiftCoefficientPressure{<:Any}}) + "CL_p" end function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:LiftCoefficient{<:Any}}) - "CL" + <:LiftCoefficientPressure{<:Any}}) + "CL_p" end + function pretty_form_ascii(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficient{<:Any}}) - "CD" + <:DragCoefficientPressure{<:Any}}) + "CD_p" end function pretty_form_utf(::AnalysisSurfaceIntegral{<:Any, - <:DragCoefficient{<:Any}}) - "CD" + <:DragCoefficientPressure{<:Any}}) + "CD_p" end end # muladd From ad4e9956ad76747b365b9b9b634272e6dbb2a256 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 27 Mar 2024 16:47:06 +0100 Subject: [PATCH 43/53] typos --- src/callbacks_step/analysis_surface_integral_2d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 441413ad04a..beb2306e820 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -68,13 +68,13 @@ C_{L,p} \coloneqq \frac{\oint_{\partial \Omega} p \boldsymbol n \cdot \psi_L \, {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 conjuction with [`AnalysisSurfaceIntegral`](@ref) +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 -- `linf::Real`: Reference length of geoemtry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ function LiftCoefficientPressure(aoa, rhoinf, uinf, linf) # psi_lift is the normal unit vector to the freestream direction. @@ -95,13 +95,13 @@ C_{D,p} \coloneqq \frac{\oint_{\partial \Omega} p \boldsymbol n \cdot \psi_D \, {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 conjuction with [`AnalysisSurfaceIntegral`](@ref) +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 -- `linf::Real`: Reference length of geoemtry (e.g. airfoil chord length) +- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) """ function DragCoefficientPressure(aoa, rhoinf, uinf, linf) # `psi_drag` is the unit vector in direction of the freestream. From 7779eb1c2a84035211b828166341629ea2d0d060 Mon Sep 17 00:00:00 2001 From: Arpit Babbar Date: Thu, 28 Mar 2024 09:45:57 +0530 Subject: [PATCH 44/53] 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 --- .../elixir_euler_NACA0012airfoil_mach085.jl | 18 +++---- .../elixir_euler_subsonic_cylinder.jl | 16 +++--- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 53 +++++++++---------- .../analysis_surface_integral_2d.jl | 16 +++--- 4 files changed, 51 insertions(+), 52 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 59c46977267..0a7d5d026ec 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/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", diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index df8782cef82..f2b8892c948 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/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", diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 80b5a840c32..7c2584fb432 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/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 diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index beb2306e820..e59a668a9c1 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -1,6 +1,3 @@ -# 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. @@ -8,15 +5,20 @@ @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 From 16fc47f620a75c9d22e637645141da5609f8eec3 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 28 Mar 2024 09:59:51 +0530 Subject: [PATCH 45/53] Some fixes --- .../elixir_euler_NACA0012airfoil_mach085.jl | 4 +- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 4 +- .../analysis_surface_integral_2d.jl | 52 +++++++++---------- 3 files changed, 30 insertions(+), 30 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 0a7d5d026ec..625512158a7 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -87,11 +87,11 @@ 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), l_inf)) lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, LiftCoefficientPressure(aoa(), rho_inf(), - u_inf(equations), linf)) + u_inf(equations), l_inf)) analysis_callback = AnalysisCallback(semi, interval = analysis_interval, output_directory = "out", diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 7c2584fb432..01adfea4343 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -5,7 +5,7 @@ using Trixi ############################################################################### # semidiscretization of the compressible Euler equations -# Laminar transonic flow around an airfoil. +# 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. @@ -113,7 +113,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers -# Run for a long time to reach a steady state +# Run for a long time to reach a state where forces stabilize up to 3 digits tspan = (0.0, 10.0) ode = semidiscretize(semi, tspan) diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index e59a668a9c1..47a32a2ce88 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -9,14 +9,14 @@ # surface forces """ - AnalysisSurfaceIntegral{Semidiscretization, Variable}(semi, - boundary_symbol_or_boundary_symbols, + 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` + 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 + 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. """ @@ -50,7 +50,7 @@ struct ForceState{RealT <: Real} psi::Tuple{RealT, RealT} # Unit vector normal or parallel to freestream rhoinf::RealT uinf::RealT - linf::RealT + l_inf::RealT end struct LiftCoefficientPressure{RealT <: Real} @@ -61,70 +61,70 @@ struct DragCoefficientPressure{RealT <: Real} force_state::ForceState{RealT} end -""" - LiftCoefficientPressure(aoa, rhoinf, uinf, linf) +""" + LiftCoefficientPressure(aoa, rhoinf, uinf, l_inf) -Compute the lift coefficient +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) +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 -- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function LiftCoefficientPressure(aoa, rhoinf, uinf, linf) +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 + # 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(ForceState(psi_lift, rhoinf, uinf, l_inf)) end -""" - DragCoefficientPressure(aoa, rhoinf, uinf, linf) +""" + DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) -Compute the drag coefficient +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) +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 -- `linf::Real`: Reference length of geometry (e.g. airfoil chord length) +- `l_inf::Real`: Reference length of geometry (e.g. airfoil chord length) """ -function DragCoefficientPressure(aoa, rhoinf, uinf, linf) +function DragCoefficientPressure(aoa, rhoinf, uinf, l_inf) # `psi_drag` is the unit vector in direction of the freestream. psi_drag = (cos(aoa), sin(aoa)) - return DragCoefficientPressure(ForceState(psi_drag, rhoinf, uinf, linf)) + 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, linf = lift_coefficient.force_state + @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 * linf) + 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, linf = drag_coefficient.force_state + @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 * linf) + return p * n / (0.5 * rhoinf * uinf^2 * l_inf) end function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, @@ -152,7 +152,7 @@ function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, 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. + # i.e., *into* the structure. normal_direction = get_normal_direction(direction, contravariant_vectors, i_node, j_node, element) From 545e42a13ea0ec731fe2879f3141f5309e941e71 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Thu, 28 Mar 2024 10:02:28 +0530 Subject: [PATCH 46/53] Format --- .../elixir_navierstokes_NACA0012airfoil_mach08.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 01adfea4343..2c61c988570 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -153,6 +153,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav ############################################################################### # run the simulation -sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1e-8, reltol = 1e-8, +sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1e-8, + reltol = 1e-8, ode_default_options()..., callback = callbacks) summary_callback() # print the timer summary From d4981c3037d14531ed6cd059b9d9df7f13cce587 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 29 Mar 2024 10:35:33 +0100 Subject: [PATCH 47/53] Update NEWS.md --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 2ade59c8039..5516fc9add5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,7 +10,7 @@ for human readability. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D`. - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension to 1D and 3D on `TreeMesh`. -- Pressure lift and drag coefficients can now be computed for 2D `P4estMesh`. +- 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 From 7376dd43993d869e1fd00eee6eab7aad2744cfbe Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 29 Mar 2024 10:36:47 +0100 Subject: [PATCH 48/53] Update src/callbacks_step/analysis_surface_integral_2d.jl --- 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 47a32a2ce88..59fc64992b7 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -106,7 +106,7 @@ which stores the boundary information and semidiscretization. - `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 in direction of the freestream. + # `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 From c156876608580237691aef49fbf043799116e406 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 29 Mar 2024 10:57:04 +0100 Subject: [PATCH 49/53] Polish --- .../elixir_euler_NACA0012airfoil_mach085.jl | 6 +++--- src/callbacks_step/analysis_surface_integral_2d.jl | 8 ++++---- test/test_p4est_2d.jl | 12 ++++++------ 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 625512158a7..e4a9d2b4c76 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -19,7 +19,7 @@ u_inf(equations) = mach_inf() * c_inf(equations) v1 = u_inf(equations) * cos(aoa()) v2 = u_inf(equations) * sin(aoa()) - prim = SVector(rho_inf(), v1, v2, pre_inf()) + prim = SVector(rho_inf(), v1, v2, p_inf()) return prim2cons(prim, equations) end @@ -107,7 +107,7 @@ save_solution = SaveSolutionCallback(interval = 500, solution_variables = cons2prim) # Small time step should be used to reach steady state -stepsize_callback = StepsizeCallback(cfl = 0.25) +stepsize_callback = StepsizeCallback(cfl = 1.0) amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) @@ -127,7 +127,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav ############################################################################### # run the simulation -sol = solve(ode, SSPRK54(), +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/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 59fc64992b7..902fb0b4e85 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -66,8 +66,8 @@ end 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}} +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) @@ -93,8 +93,8 @@ end 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}} +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) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index cd4e6f7759b..1bacccbf812 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -574,12 +574,12 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), l2=[ - 5.37640294704164e-7, 6.418954697992308e-6, - 1.0332595763357066e-5, 0.0006353813399881882, + 5.371568111383228e-7, 6.4158131303956445e-6, + 1.0324346542348325e-5, 0.0006348064933187732, ], linf=[ - 0.0016284523634016628, 0.028505821812697323, - 0.029918806073518, 1.9505653217814127, + 0.0016263400091978443, 0.028471072159724428, + 0.02986133204785877, 1.9481060511014872, ], base_level=0, med_level=1, max_level=1, tspan=(0.0, 0.0001), @@ -608,8 +608,8 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache) - @test isapprox(lift, 0.026397498239816828, atol = 1e-13) - @test isapprox(drag, 0.10908833968266139, atol = 1e-13) + @test isapprox(lift, 0.0262382560809345, atol = 1e-13) + @test isapprox(drag, 0.10898248971932244, atol = 1e-13) end end From e76c617b787944665199537bd6b362342ae607e8 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Fri, 29 Mar 2024 10:59:31 +0100 Subject: [PATCH 50/53] Update examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl --- examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index e4a9d2b4c76..fb5f29bd038 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -106,7 +106,6 @@ save_solution = SaveSolutionCallback(interval = 500, save_final_solution = true, solution_variables = cons2prim) -# Small time step should be used to reach steady state stepsize_callback = StepsizeCallback(cfl = 1.0) amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) From 62cc4c8e5cf778d9da521ae10f9f145c47057d70 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 29 Mar 2024 11:06:54 +0100 Subject: [PATCH 51/53] increase cfl --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index f2b8892c948..fd0e1bd2c11 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -110,15 +110,14 @@ save_solution = SaveSolutionCallback(interval = 500, save_final_solution = true, solution_variables = cons2prim) -# Small time step should be used to reach steady state -stepsize_callback = StepsizeCallback(cfl = 0.25) +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), +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 From 815a7244ccf96e4d0ed113a22400f0a3d23a0566 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 29 Mar 2024 11:08:47 +0100 Subject: [PATCH 52/53] fmt --- examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index fd0e1bd2c11..dc23e192de8 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -117,7 +117,9 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav ############################################################################### # run the simulation -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false; thread = OrdinaryDiffEq.True()), +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 From 28104996e917d12696b9b8e7bdcf9ce551b36e4a Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 29 Mar 2024 18:32:26 +0100 Subject: [PATCH 53/53] fix --- .../elixir_navierstokes_NACA0012airfoil_mach08.jl | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 2c61c988570..5bfa21c5322 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -35,16 +35,12 @@ 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 - mach_inf = mach_inf() - aoa = aoa() - rho_inf = rho_inf() - p_inf = p_inf() - u_inf = mach_inf * sqrt(gamma * p_inf / rho_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, p_inf) + prim = SVector(rho_inf(), v1, v2, p_inf()) return prim2cons(prim, equations) end @@ -114,7 +110,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol # ODE solvers # Run for a long time to reach a state where forces stabilize up to 3 digits -tspan = (0.0, 10.0) +tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) # Callbacks