Boundary integration of viscous forces

examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl

@@ -110,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, 1.0)
+tspan = (0.0, 10.0)
 ode = semidiscretize(semi, tspan)
 
 # Callbacks
@@ -130,12 +130,26 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
                                                                    u_inf(equations),
                                                                    l_inf()))
 
+drag_coefficient_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+                                                 DragCoefficientShearStress(aoa(),
+                                                                            rho_inf(),
+                                                                            u_inf(equations),
+                                                                            l_inf()))
+
+lift_coefficient_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+                                                 LiftCoefficientShearStress(aoa(),
+                                                                            rho_inf(),
+                                                                            u_inf(equations),
+                                                                            l_inf()))
+
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      output_directory = "out",
                                      save_analysis = true,
                                      analysis_errors = Symbol[],
                                      analysis_integrals = (drag_coefficient,
-                                                           lift_coefficient))
+                                                           lift_coefficient,
+                                                           drag_coefficient_force,
+                                                           lift_coefficient_force))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 

src/Trixi.jl

@@ -262,7 +262,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        AMRCallback, StepsizeCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled,
-       AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure
+       AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure,
+       DragCoefficientShearStress, LiftCoefficientShearStress
 
 export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
        load_adaptive_time_integrator!

src/callbacks_step/analysis.jl

@@ -9,11 +9,11 @@
 # - analysis_interval part as PeriodicCallback called after a certain amount of simulation time
 """
     AnalysisCallback(semi; interval=0,
-                           save_analysis=false,
-                           output_directory="out",
-                           analysis_filename="analysis.dat",
-                           extra_analysis_errors=Symbol[],
-                           extra_analysis_integrals=())
+                            save_analysis=false,
+                            output_directory="out",
+                            analysis_filename="analysis.dat",
+                            extra_analysis_errors=Symbol[],
+                            extra_analysis_integrals=())
 
 Analyze a numerical solution every `interval` time steps and print the
 results to the screen. If `save_analysis`, the results are also saved in
@@ -633,24 +633,6 @@ end
 pretty_form_utf(quantity) = get_name(quantity)
 pretty_form_ascii(quantity) = get_name(quantity)
 
-# Special analyze for `SemidiscretizationHyperbolicParabolic` such that
-# precomputed gradients are available. For now only implemented for the `enstrophy`
-#!!! warning "Experimental code"
-#    This code is experimental and may be changed or removed in any future release.
-function analyze(quantity::typeof(enstrophy), du, u, t,
-                 semi::SemidiscretizationHyperbolicParabolic)
-    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
-    equations_parabolic = semi.equations_parabolic
-    cache_parabolic = semi.cache_parabolic
-    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
-            cache_parabolic)
-end
-function analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver,
-                 cache, cache_parabolic)
-    integrate(quantity, u, mesh, equations, equations_parabolic, solver, cache,
-              cache_parabolic, normalize = true)
-end
-
 function entropy_timederivative end
 pretty_form_utf(::typeof(entropy_timederivative)) = "∑∂S/∂U ⋅ Uₜ"
 pretty_form_ascii(::typeof(entropy_timederivative)) = "dsdu_ut"
@@ -695,3 +677,25 @@ include("analysis_surface_integral_2d.jl")
 include("analysis_dg2d_parallel.jl")
 include("analysis_dg3d.jl")
 include("analysis_dg3d_parallel.jl")
+
+# Special analyze for `SemidiscretizationHyperbolicParabolic` such that
+# precomputed gradients are available. Required for `enstrophy` and viscous forces.
+# Note that this needs to be included after `analysis_surface_integral_2d.jl` to 
+# have `VariableViscous` available.
+function analyze(quantity::Union{typeof(enstrophy),
+                                 AnalysisSurfaceIntegral{Semidiscretization, Variable}},
+                 du, u, t,
+                 semi::SemidiscretizationHyperbolicParabolic) where {Semidiscretization,
+                                                                     Variable <:
+                                                                     VariableViscous}
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    equations_parabolic = semi.equations_parabolic
+    cache_parabolic = semi.cache_parabolic
+    analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache,
+            cache_parabolic)
+end
+function analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver,
+                 cache, cache_parabolic)
+    integrate(quantity, u, mesh, equations, equations_parabolic, solver, cache,
+              cache_parabolic, normalize = true)
+end