From 9e95c37719af25da6990bc6b592d014b0d246b38 Mon Sep 17 00:00:00 2001 From: Arpit-Babbar Date: Wed, 24 Jan 2024 21:08:28 +0530 Subject: [PATCH] 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 cebc2917d52..a92b1c8aa5e 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -510,6 +510,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