2D field in Julia

2d.jl

using FFTW
using DiffEqSensitivity, OrdinaryDiffEq, Zygote

nlon = 64
nlat = 64
K = (1:(div(nlat,2)+1)) .* (1:nlon)' / (nlat*nlon)

function f(x,p,t)
    k = p[1]
    lc = k * irfft(K .* rfft(x), nlat)
    @. x - x^3/3 + lc
end

p = [1.5];
u0 = randn(nlat, nlon);
prob = ODEProblem(f, u0, (0.0, 10.0), p);
loss(u0,p) = sum(solve(prob,Tsit5(),u0=u0,p=p,saveat=0.1))
du01,dp1 = Zygote.gradient(loss,u0,p)

pkg status

(@v1.6) pkg> status
      Status `~/.julia/environments/v1.6/Project.toml`
  [6e4b80f9] BenchmarkTools v1.2.0
  [aae7a2af] DiffEqFlux v1.43.0
  [41bf760c] DiffEqSensitivity v6.59.0
  [7a1cc6ca] FFTW v1.4.5
  [587475ba] Flux v0.12.7
  [28b8d3ca] GR v0.61.0
  [15e1cf62] NPZ v0.4.1
  [1dea7af3] OrdinaryDiffEq v5.64.1
  [91a5bcdd] Plots v1.22.6
  [c3e4b0f8] Pluto v0.16.1
  [29a0d76e] StochasticDelayDiffEq v1.1.1
  [e88e6eb3] Zygote v0.6.27

traceback.jl

ERROR: LoadError: type ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}} not supported
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] _fftfloat(#unused#::Type{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:22
  [3] _fftfloat(x::ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:23
  [4] fftfloat(x::ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:18
  [5] realfloat(x::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:33
  [6] plan_rfft(x::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}, region::UnitRange{Int64}; kws::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:205
  [7] plan_rfft(x::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}, region::UnitRange{Int64})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:205
  [8] plan_rfft(x::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}; kws::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:52
  [9] plan_rfft(x::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:52
 [10] rfft(x::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}})
    @ AbstractFFTs ~/.julia/packages/AbstractFFTs/JebmH/src/definitions.jl:50
 [11] f(dx::Matrix{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}}, x::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}, p::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, t::ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}})
    @ Main ~/src/nfs/2d.jl:10
 [12] (::ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing})(::Matrix{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}}, ::Vararg{Any, N} where N)
    @ SciMLBase ~/.julia/packages/SciMLBase/h4Gxc/src/scimlfunctions.jl:334
 [13] (::DiffEqSensitivity.var"#110#125"{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}})(u::ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}, p::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, t::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/adjoint_common.jl:146
 [14] ReverseDiff.GradientTape(f::Function, input::Tuple{Matrix{Float64}, Vector{Float64}, Vector{Float64}}, cfg::ReverseDiff.GradientConfig{Tuple{ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/E4Tzn/src/api/tape.jl:207
 [15] ReverseDiff.GradientTape(f::Function, input::Tuple{Matrix{Float64}, Vector{Float64}, Vector{Float64}})
    @ ReverseDiff ~/.julia/packages/ReverseDiff/E4Tzn/src/api/tape.jl:204
 [16] adjointdiffcache(g::DiffEqSensitivity.var"#df#233"{FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, discrete::Bool, sol::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, dg::Nothing, f::ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}; quad::Bool, noiseterm::Bool)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/adjoint_common.jl:144
 [17] DiffEqSensitivity.ODEInterpolatingAdjointSensitivityFunction(g::Function, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, discrete::Bool, sol::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, dg::Nothing, f::Function, checkpoints::Vector{Float64}, tols::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, tstops::Nothing; noiseterm::Bool)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/interpolating_adjoint.jl:72
 [18] DiffEqSensitivity.ODEInterpolatingAdjointSensitivityFunction(g::Function, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, discrete::Bool, sol::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, dg::Nothing, f::Function, checkpoints::Vector{Float64}, tols::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, tstops::Nothing)
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/interpolating_adjoint.jl:23
 [19] ODEAdjointProblem(sol::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, g::DiffEqSensitivity.var"#df#233"{FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing; checkpoints::Vector{Float64}, callback::Nothing, reltol::Float64, abstol::Float64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/interpolating_adjoint.jl:255
 [20] _adjoint_sensitivities(sol::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, g::DiffEqSensitivity.var"#df#233"{FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}}, Colon}, t::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, dg::Nothing; abstol::Float64, reltol::Float64, checkpoints::Vector{Float64}, corfunc_analytical::Nothing, callback::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/sensitivity_interface.jl:17
 [21] adjoint_sensitivities(::ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}}, DiffEqBase.DEStats}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, ::Vararg{Any, N} where N; sensealg::InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, kwargs::Base.Iterators.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/sensitivity_interface.jl:6
 [22] (::DiffEqSensitivity.var"#adjoint_sensitivity_backpass#232"{Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, Matrix{Float64}, Vector{Float64}, Tuple{}, Colon, NamedTuple{(), Tuple{}}})(Δ::FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ DiffEqSensitivity ~/.julia/packages/DiffEqSensitivity/eq2wv/src/concrete_solve.jl:231
 [23] ZBack
    @ ~/.julia/packages/Zygote/Y1xbK/src/compiler/chainrules.jl:168 [inlined]
 [24] (::Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#232"{Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, Matrix{Float64}, Vector{Float64}, Tuple{}, Colon, NamedTuple{(), Tuple{}}}})(dy::FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/chainrules.jl:194
 [25] #203
    @ ~/.julia/packages/Zygote/Y1xbK/src/lib/lib.jl:203 [inlined]
 [26] (::Zygote.var"#1733#back#205"{Zygote.var"#203#204"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, Zygote.var"#kw_zpullback#40"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#232"{Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, InterpolatingAdjoint{0, true, Val{:central}, ReverseDiffVJP{false}, Bool}, Matrix{Float64}, Vector{Float64}, Tuple{}, Colon, NamedTuple{(), Tuple{}}}}}})(Δ::FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [27] Pullback
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:73 [inlined]
 [28] (::typeof(∂(#solve#43)))(Δ::FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface2.jl:0
 [29] (::Zygote.var"#203#204"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))})(Δ::FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/lib/lib.jl:203
 [30] (::Zygote.var"#1733#back#205"{Zygote.var"#203#204"{Tuple{NTuple{6, Nothing}, Tuple{Nothing}}, typeof(∂(#solve#43))}})(Δ::FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [31] Pullback
    @ ~/.julia/packages/DiffEqBase/OPDgm/src/solve.jl:68 [inlined]
 [32] (::typeof(∂(solve##kw)))(Δ::FillArrays.Fill{Float64, 3, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}})
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface2.jl:0
 [33] Pullback
    @ ~/src/nfs/2d.jl:17 [inlined]
 [34] (::typeof(∂(loss)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface2.jl:0
 [35] (::Zygote.var"#50#51"{typeof(∂(loss))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface.jl:41
 [36] gradient(::Function, ::Matrix{Float64}, ::Vararg{Any, N} where N)
    @ Zygote ~/.julia/packages/Zygote/Y1xbK/src/compiler/interface.jl:76
 [37] top-level scope
    @ ~/src/nfs/2d.jl:18
in expression starting at ~/src/nfs/2d.jl:18