Here are two other thoughts:
Maybe you can keep the graph of h on GPU but only move to to disk the part that computes the matrix-Jacobian product.
model = quadratic_fun()
loss_quad = model.forward()
grad_ft = torch.autograd.grad(loss_quad, model.parameters(), create_graph=True)
flat_grad = torch.cat([g.contiguous().view(-1) for g in grad_ft])

# generate the constant matrix V, and compute the matrix-gradient product
torch.manual_seed(0)
V = torch.randn((10, 3))
h = torch.matmul(flat_grad, V)

torch.autograd.graph.set_saved_tensors_default_hooks(pack_hook, unpack_hook)
# compute the matrix-Jacobian product by iterating over the columns of the constant matrix
for i in range(3):
    hvp = torch.autograd.grad(h[i], model.parameters(), retain_graph=True)
    hvp_flat = torch.cat([g.contiguous().view(-1) for g in hvp])
    print(hvp_flat)

torch.autograd.graph.reset_saved_tensors_default_hooks()

or even, instead of moving to disk, moving to CPU:
model = quadratic_fun()
loss_quad = model.forward()
grad_ft = torch.autograd.grad(loss_quad, model.parameters(), create_graph=True)
flat_grad = torch.cat([g.contiguous().view(-1) for g in grad_ft])

# generate the constant matrix V, and compute the matrix-gradient product
torch.manual_seed(0)
V = torch.randn((10, 3))
h = torch.matmul(flat_grad, V)

with torch.autograd.graph.save_on_cpu(pin_memory=True):
    for i in range(3):
        hvp = torch.autograd.grad(h[i], model.parameters(), retain_graph=True)
        hvp_flat = torch.cat([g.contiguous().view(-1) for g in hvp])
        print(hvp_flat)