Viscous force computation

trixi-framework · Mar 31, 2024 · 4fd1593 · 4fd1593
1 parent 57f4b20
commit 4fd1593
Show file tree

Hide file tree

Showing 3 changed files with 248 additions and 47 deletions.
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/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,24 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
                                                                    u_inf(equations),
                                                                    l_inf()))
 
+drag_coefficient_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+                                                 Trixi.DragCoefficientShearStress(aoa(), rho_inf(),
+                                                                      u_inf(equations),
+                                                                      l_inf()))
+
+lift_coefficient_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+                                                 Trixi.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)
 
@@ -144,7 +156,7 @@ save_solution = SaveSolutionCallback(interval = 500,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)#, save_solution)
 
 ###############################################################################
 # run the simulation

diff --git a/src/callbacks_step/analysis.jl b/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,26 @@ 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