Sure, though I am not sure whether I can convince you or not:

I try not to use np.array as much as possible and define them as  torch.tensor(dtype=torch.float64), such as x5 and omega5 in my example.

eval_g returns [g0], not np.array([g0]).
I do not use jacobian script you proposed above. Instead I define gradient by myself:

x_tensor = torch.from_numpy(x).requires_grad_(True)
jac = []
    for i in range(ncon):
        grad, = torch.autograd.grad(eval_g(x_tensor)[i], x_tensor)
        jac.append(grad)
    return torch.stack(jac).flatten().numpy()


Finally I use torch.autograd.gradcheck function with the following syntax

k_tensor = torch.from_numpy(np.array([k])).requires_grad_(True)
chk = torch.autograd.gradcheck(eval_g, k_tensor)
print('Is the gradient check passed? {}'.format(chk))

chk returns True