Jacobian via ForwardDiff for parabolic part only (#2513)

* Jacobian via ForwardDiff for parabolic part only

* flip

* Update test/test_special_elixirs.jl

Author: DanielDoehring

src/Trixi.jl

@@ -309,7 +309,9 @@ export trixi_include, examples_dir, get_examples, default_example,
 
 export ode_norm, ode_unstable_check
 
-export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure
+export convergence_test,
+       jacobian_fd, jacobian_ad_forward, jacobian_ad_forward_parabolic,
+       linear_structure
 
 export DGMulti, DGMultiBasis, estimate_dt, DGMultiMesh, GaussSBP
 

src/semidiscretization/semidiscretization.jl

@@ -220,7 +220,7 @@ end
 
 Uses the right-hand side operator of the semidiscretization `semi`
 and simple second order finite difference to compute the Jacobian `J`
-of the semidiscretization `semi` at state `u0_ode`.
+of the semidiscretization `semi` at time `t0` and state `u0_ode`.
 """
 function jacobian_fd(semi::AbstractSemidiscretization;
                      t0 = zero(real(semi)),
@@ -267,7 +267,7 @@ end
 
 Uses the right-hand side operator of the semidiscretization `semi`
 and forward mode automatic differentiation to compute the Jacobian `J`
-of the semidiscretization `semi` at state `u0_ode`.
+of the semidiscretization `semi` at time `t0` and state `u0_ode`.
 """
 function jacobian_ad_forward(semi::AbstractSemidiscretization;
                              t0 = zero(real(semi)),

src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl

@@ -375,4 +375,45 @@ function _jacobian_ad_forward(semi::SemidiscretizationHyperbolicParabolic, t0, u
 
     return J
 end
+
+"""
+    jacobian_ad_forward_parabolic(semi::SemidiscretizationHyperbolicParabolic;
+                                  t0=zero(real(semi)),
+                                  u0_ode=compute_coefficients(t0, semi))
+
+Uses the *parabolic part* of the right-hand side operator of the [`SemidiscretizationHyperbolicParabolic`](@ref) `semi`
+and forward mode automatic differentiation to compute the Jacobian `J` of the 
+parabolic/diffusive contribution only at time `t0` and state `u0_ode`.
+
+This might be useful for operator-splitting methods, e.g., the construction of optimized 
+time integrators which optimize different methods for the hyperbolic and parabolic part separately.
+"""
+function jacobian_ad_forward_parabolic(semi::SemidiscretizationHyperbolicParabolic;
+                                       t0 = zero(real(semi)),
+                                       u0_ode = compute_coefficients(t0, semi))
+    jacobian_ad_forward_parabolic(semi, t0, u0_ode)
+end
+
+# The following version is for plain arrays
+function jacobian_ad_forward_parabolic(semi::SemidiscretizationHyperbolicParabolic,
+                                       t0, u0_ode)
+    du_ode = similar(u0_ode)
+    config = ForwardDiff.JacobianConfig(nothing, du_ode, u0_ode)
+
+    # Use a function barrier since the generation of the `config` we use above
+    # is not type-stable
+    _jacobian_ad_forward_parabolic(semi, t0, u0_ode, du_ode, config)
+end
+
+function _jacobian_ad_forward_parabolic(semi, t0, u0_ode, du_ode, config)
+    new_semi = remake(semi, uEltype = eltype(config))
+    # Create anonymous function passed as first argument to `ForwardDiff.jacobian` to match
+    # `ForwardDiff.jacobian(f!, y::AbstractArray, x::AbstractArray, 
+    #                       cfg::JacobianConfig = JacobianConfig(f!, y, x), check=Val{true}())`
+    J = ForwardDiff.jacobian(du_ode, u0_ode, config) do du_ode, u_ode
+        Trixi.rhs_parabolic!(du_ode, u_ode, new_semi, t0)
+    end
+
+    return J
+end
 end # @muladd

test/test_special_elixirs.jl

@@ -138,6 +138,12 @@ end
         J = jacobian_ad_forward(semi)
         λ = eigvals(J)
         @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
+
+        J_parabolic = jacobian_ad_forward_parabolic(semi)
+        λ_parabolic = eigvals(J_parabolic)
+        # Parabolic spectrum is real and negative
+        @test maximum(real, λ_parabolic) < 10^(-14)
+        @test maximum(imag, λ_parabolic) < 10^(-14)
     end
 
     @timed_testset "Compressible Euler equations" begin
@@ -219,6 +225,12 @@ end
         J = jacobian_ad_forward(semi)
         λ = eigvals(J)
         @test maximum(real, λ) < 0.2
+
+        J_parabolic = jacobian_ad_forward_parabolic(semi)
+        λ_parabolic = eigvals(J_parabolic)
+        # Parabolic spectrum is real and negative
+        @test maximum(real, λ_parabolic) < 10^(-16)
+        @test maximum(imag, λ_parabolic) < 10^(-15)
     end
 
     @timed_testset "MHD" begin