Boundary integration of viscous forces (trixi-framework#1893)

* Viscous force computation * Modularize stress tensor comp for skin friction * typo * try aqua * remove semi from AnalysisSurfaceIntegral * bugfix * two methods for passing aqua tests * rename * test viscous forces * Update src/callbacks_step/analysis_surface_integral_2d.jl Co-authored-by: Andrew Winters <[email protected]> * docstring * try to get test to run * Try nicer formatting * revert fmt --------- Co-authored-by: Andrew Winters <[email protected]>
DanielDoehring · Apr 8, 2024 · c025873 · c025873
1 parent b4fb1db
commit c025873
Show file tree

Hide file tree

Showing 5 changed files with 296 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,26 @@ lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
                                                                    u_inf(equations),
                                                                    l_inf()))
 
+drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+                                                       DragCoefficientShearStress(aoa(),
+                                                                                  rho_inf(),
+                                                                                  u_inf(equations),
+                                                                                  l_inf()))
+
+lift_coefficient_shear_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_shear_force,
+                                                           lift_coefficient_shear_force))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 

diff --git a/src/Trixi.jl b/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!

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
@@ -634,9 +634,7 @@ 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.
+# precomputed gradients are available.
 function analyze(quantity::typeof(enstrophy), du, u, t,
                  semi::SemidiscretizationHyperbolicParabolic)
     mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
@@ -695,3 +693,19 @@ 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` (see above) and viscous forces.
+# Note that this needs to be included after `analysis_surface_integral_2d.jl` to 
+# have `VariableViscous` available.
+function analyze(quantity::AnalysisSurfaceIntegral{Variable},
+                 du, u, t,
+                 semi::SemidiscretizationHyperbolicParabolic) where {
+                                                                     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