Adapt jacobian_ad_forward for hyperbolic-parabolic semidiscretizations

examples/tree_2d_dgsem/elixir_navierstokes_taylor_green_vortex.jl

@@ -0,0 +1,74 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 6.25e-4 # equivalent to Re = 1600
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(),
+                                                          Prandtl=prandtl_number())
+
+"""
+    initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D)
+
+The classical inviscid Taylor-Green vortex.
+"""
+function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations2D)
+  A  = 1.0 # magnitude of speed
+  Ms = 0.1 # maximum Mach number
+
+  rho = 1.0
+  v1  =  A * sin(x[1]) * cos(x[2])
+  v2  = -A * cos(x[1]) * sin(x[2])
+  p   = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
+  p   = p + 1.0/4.0 * A^2 * rho * (cos(2*x[1]) + cos(2*x[2]))
+
+  return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_taylor_green_vortex
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg=3, surface_flux=flux_hllc,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0, -1.0) .* pi
+coordinates_max = ( 1.0,  1.0) .* pi
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                n_cells_max=100_000)
+
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 20.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))
+
+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/semidiscretization/semidiscretization.jl

@@ -251,15 +251,6 @@ function jacobian_ad_forward(semi::AbstractSemidiscretization, t0, u0_ode)
     _jacobian_ad_forward(semi, t0, u0_ode, du_ode, config)
 end
 
-function _jacobian_ad_forward(semi, t0, u0_ode, du_ode, config)
-    new_semi = remake(semi, uEltype = eltype(config))
-    J = ForwardDiff.jacobian(du_ode, u0_ode, config) do du_ode, u_ode
-        Trixi.rhs!(du_ode, u_ode, new_semi, t0)
-    end
-
-    return J
-end
-
 # This version is specialized to `StructArray`s used by some `DGMulti` solvers.
 # We need to convert the numerical solution vectors since ForwardDiff cannot
 # handle arrays of `SVector`s.

src/semidiscretization/semidiscretization_hyperbolic.jl

@@ -353,4 +353,14 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolic, t)
 
     return nothing
 end
+
+function _jacobian_ad_forward(semi::SemidiscretizationHyperbolic, t0, u0_ode, du_ode,
+                              config)
+    new_semi = remake(semi, uEltype = eltype(config))
+    J = ForwardDiff.jacobian(du_ode, u0_ode, config) do du_ode, u_ode
+        Trixi.rhs!(du_ode, u_ode, new_semi, t0)
+    end
+
+    return J
+end
 end # @muladd

src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl

@@ -330,4 +330,25 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol
 
     return nothing
 end
+
+function rhs_hyperbolic_parabolic!(du_ode, u_ode,
+                                   semi::SemidiscretizationHyperbolicParabolic, t)
+    @trixi_timeit timer() "hyperbolic-parabolic rhs!" begin
+        # Implementation of split ODE problem in OrdinaryDiffEq
+        du_ode_hyp = similar(du_ode)
+        rhs!(du_ode_hyp, u_ode, semi, t)
+        rhs_parabolic!(du_ode, u_ode, semi, t)
+        du_ode .+= du_ode_hyp
+    end
+end
+
+function _jacobian_ad_forward(semi::SemidiscretizationHyperbolicParabolic, t0, u0_ode,
+                              du_ode, config)
+    new_semi = remake(semi, uEltype = eltype(config))
+    J = ForwardDiff.jacobian(du_ode, u0_ode, config) do du_ode, u_ode
+        Trixi.rhs_hyperbolic_parabolic!(du_ode, u_ode, new_semi, t0)
+    end
+
+    return J
+end
 end # @muladd

test/test_special_elixirs.jl

@@ -107,6 +107,15 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
       @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
     end
 
+    @timed_testset "Linear advection-diffusion" begin
+      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion.jl"),
+                    tspan=(0.0, 0.0), initial_refinement_level=2)
+
+      J = jacobian_ad_forward(semi)
+      λ = eigvals(J)
+      @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
+    end
+
     @timed_testset "Compressible Euler equations" begin
       trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_euler_density_wave.jl"),
                     tspan=(0.0, 0.0), initial_refinement_level=1)
@@ -165,6 +174,15 @@ coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none",
       end
     end
 
+    @timed_testset "Navier-Stokes" begin
+      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_navierstokes_taylor_green_vortex.jl"),
+                    tspan=(0.0, 0.0), initial_refinement_level=2)
+
+      J = jacobian_ad_forward(semi)
+      λ = eigvals(J)
+      @test maximum(real, λ) < 0.2
+    end
+
     @timed_testset "MHD" begin
       trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_mhd_alfven_wave.jl"),
                     tspan=(0.0, 0.0), initial_refinement_level=0)