See the following example. I do calculate the derivative of f w.r.t. to x, called fx. Then I calculate the derivative of fx w.r.t. to y, called fxy. The function fxyz is computed similarly.
The output of fx is correct, but fxy and fxyz give incorrect results. What is the proper way to do this with ReverseDiff?
using ReverseDiff: gradient
# the function of interest
f(x, y, z) = x.*y.*z
fx(x, y, z) = gradient(x->f(x, y, z), x)
fxy(x, y, z) = gradient(y->fx(x, y, z), y)
fxyz(x, y, z) = gradient(z->fxy(x, y, z), z)
inputs = ([1.f0], [2.f0], [3.f0])
fx(inputs...) # output is [6.f0], CORRECT
fxy(inputs...) # output is [0.f0], WRONG, shall be [3.f0]
fxyz(inputs...) # output is [0.f0], WRONG, shall be [1.f0]