Enstrophy for 2D Navier-Stokes (#1591)

* Doubly periodic shear layer

* test if prject toml shows up in git diff

* remove chnages

* Enstrophy for 2D Navier-Stokes

Project.toml

@@ -60,8 +60,8 @@ HDF5 = "0.14, 0.15, 0.16"
 IfElse = "0.1"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.118"
-Makie = "0.19"
 MPI = "0.20"
+Makie = "0.19"
 MuladdMacro = "0.2.2"
 Octavian = "0.3.5"
 OffsetArrays = "1.3"

examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl

@@ -0,0 +1,71 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 1.0/3.0 * 10^(-3) # equivalent to Re = 3000
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
+                                                          Prandtl=prandtl_number())
+
+function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D)
+  k = 80
+  delta = 0.05
+  u0 = 1.0
+  Ms = 0.1 # maximum Mach number
+
+  rho = 1.0
+  v1  = x[2] <= 0.5 ? u0*tanh(k*(x[2]*0.5 - 0.25)) : tanh(k*(0.75 -x[2]*0.5))
+  v2  = u0*delta * sin(2*pi*(x[1]*0.5 + 0.25))
+  p   = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_shear_layer
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg=3, surface_flux=flux_hllc,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=100_000)
+
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 50
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
+                                     extra_analysis_integrals=(energy_kinetic,
+                                                               energy_internal,
+                                                               enstrophy))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval,)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,
+            ode_default_options()..., callback=callbacks)
+summary_callback() # print the timer summary

src/callbacks_step/analysis_dg2d.jl

@@ -213,6 +213,24 @@ function integrate(func::Func, u,
     end
 end
 
+function integrate(func::Func, u,
+                   mesh::Union{TreeMesh{2}, P4estMesh{2}},
+                   equations, equations_parabolic,
+                   dg::DGSEM,
+                   cache, cache_parabolic; normalize = true) where {Func}
+    gradients_x, gradients_y = cache_parabolic.gradients
+    integrate_via_indices(u, mesh, equations, dg, cache;
+                          normalize = normalize) do u, i, j, element, equations, dg
+        u_local = get_node_vars(u, equations, dg, i, j, element)
+        gradients_1_local = get_node_vars(gradients_x, equations_parabolic, dg, i, j,
+                                          element)
+        gradients_2_local = get_node_vars(gradients_y, equations_parabolic, dg, i, j,
+                                          element)
+        return func(u_local, (gradients_1_local, gradients_2_local),
+                    equations_parabolic)
+    end
+end
+
 function analyze(::typeof(entropy_timederivative), du, u, t,
                  mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
                              P4estMesh{2}, T8codeMesh{2}},

src/equations/compressible_navier_stokes_2d.jl

@@ -300,6 +300,21 @@ end
     return T
 end
 
+@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion2D)
+    # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v
+
+    omega = vorticity(u, gradients, equations)
+    return 0.5 * u[1] * omega^2
+end
+
+@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion2D)
+    # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function.
+    _, dv1dx, dv2dx, _ = convert_derivative_to_primitive(u, gradients[1], equations)
+    _, dv1dy, dv2dy, _ = convert_derivative_to_primitive(u, gradients[2], equations)
+
+    return dv2dx - dv1dy
+end
+
 # TODO: can we generalize this to MHD?
 """
     struct BoundaryConditionNavierStokesWall

test/test_parabolic_2d.jl

@@ -136,6 +136,10 @@ isdir(outdir) && rm(outdir, recursive=true)
   @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"),
       initial_refinement_level = 2, tspan=(0.0, 0.1),
+      analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+      extra_analysis_integrals=(energy_kinetic,
+                                energy_internal,
+                                enstrophy)),
       l2 = [0.002111672530658797, 0.0034322351490857846, 0.0038742528195910416, 0.012469246082568561],
       linf = [0.012006418939223495, 0.035520871209746126, 0.024512747492231427, 0.11191122588756564]
     )