Gradient calculation in PINN

I may have found the information I was looking for on how to find the first derivative in the discussion shown in the link below.

Gradient of Flux model wrt to weights

Based on this information, the part of the PINN programme I would like to create would look like this.

using Flux
model = Flux.Dense(24,32)
baseline = Flux.params(model)
input = [-1 1 -1 1 -1 1 -1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 -1 -1 1 1 1 1]
output = model(input)
grad = Flux.jacobian(x -> model(x), input)

In the program, input is location of 8 points in 3-dimensional domain (X[-1, 1], Y[-1, 1], Z[-1, 1]), and output is velocity in X, Y, Z direction and pressure at the each points.

I am looking for a way to calculate the second order derivative of the output to the input of this neural network, and when I use the Hessian function I get the following error message

grad2 = Flux.hessian(x -> model(x), some_input)
ERROR: output an array, so the gradient is not defined. Perhaps you wanted jacobian.

Could you tell me how to solve this error?