Fix labels, add tests

trixi-framework · Jan 24, 2024 · 2badb17 · 2badb17
1 parent 05d5bcf
commit 2badb17
Show file tree

Hide file tree

Showing 3 changed files with 52 additions and 29 deletions.
diff --git 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,14 @@ 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
@@ -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
@@ -510,6 +510,32 @@ 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