Improve computation of epsilon in jacobian_fd (#2522)

* improve computation of epsilon in jacobian_fd

* add link to comment

* add tests

Author: JoshuaLampert

src/semidiscretization/semidiscretization.jl

@@ -240,7 +240,13 @@ function jacobian_fd(semi::AbstractSemidiscretization;
     # use second order finite difference to estimate Jacobian matrix
     for idx in eachindex(u0_ode)
         # determine size of fluctuation
-        epsilon = sqrt(eps(typeof(u0_ode[idx])))
+        # This is the approach used by FiniteDiff.jl to compute the
+        # step size, which assures that the finite difference is accurate
+        # for very small and very large absolute values `u0_ode[idx]`.
+        # See https://github.com/trixi-framework/Trixi.jl/pull/2514#issuecomment-3190534904.
+        absstep = sqrt(eps(typeof(u0_ode[idx])))
+        relstep = absstep
+        epsilon = max(relstep * abs(u0_ode[idx]), absstep)
 
         # plus fluctuation
         u_ode[idx] = u0_ode[idx] + epsilon

test/test_special_elixirs.jl

@@ -110,6 +110,37 @@ end
     @test A * x ≈ Ax
 end
 
+@testset "Test Jacobian of DG (1D)" begin
+    @timed_testset "Linear advection" begin
+        trixi_include(@__MODULE__,
+                      joinpath(EXAMPLES_DIR, "tree_1d_fdsbp",
+                               "elixir_advection_upwind.jl"),
+                      tspan = (0.0, 0.0))
+
+        A, _ = linear_structure(semi)
+
+        J = jacobian_ad_forward(semi)
+        @test Matrix(A) ≈ J
+        λ = eigvals(J)
+        @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
+
+        J = jacobian_fd(semi)
+        @test Matrix(A) ≈ J
+        λ = eigvals(J)
+        @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
+
+        # See https://github.com/trixi-framework/Trixi.jl/pull/2514
+        @test count(real.(λ) .>= -10) > 5
+        # See https://github.com/trixi-framework/Trixi.jl/pull/2522
+        t0 = zero(real(semi))
+        u0_ode = 1e9 * compute_coefficients(t0, semi)
+        J = jacobian_fd(semi; t0, u0_ode)
+        λ = eigvals(J)
+        @test count((-200 .<= real.(λ) .<= -10) .&& (-100 .<= imag.(λ) .<= 100)) == 0
+        @test count(isapprox.(imag.(λ), 0.0, atol = 10 * sqrt(eps(real(semi))))) == 2
+    end
+end
+
 @testset "Test Jacobian of DG (2D)" begin
     @timed_testset "Linear advection" begin
         trixi_include(@__MODULE__,