DimensionMismatch when taking hessian of vector function using DiffResults

Hessian of vector function can be calculated by composing two jacobians:

f(x) = x.^2
g(x) = ForwardDiff.jacobian(f,x)
x = rand(3)
ForwardDiff.jacobian(g,x)

But if I want to retrieve the first jacobian using DiffResults,

result = DiffResults.JacobianResult(x)
result = ForwardDiff.jacobian!(result,g,x)

This would give

ERROR: DimensionMismatch("new dimensions (9, 3) must be consistent with array size 9")

Any workaround?

Start by reproducing the nested Jacobian example with ForwardDiff.jacobian! and DiffResults.JacobianResult shown in the issue. Inspect those entry points and how the result dimensions are handled. Done means explaining or resolving the DimensionMismatch for this vector-function Hessian case, with a verified workaround or regression test if the project’s existing tests cover it.