Solution (source code)

= Solution

Direct <differentiation> gives
$$
\nabla a=(y,x,2z),\qquad \nabla\times\mathbf B=(y,-1,-1),\qquad \nabla a\cdot(\nabla\times\mathbf B)=y^2-x-2z.
$$
On the unit disk $D$, <polar coordinates> give $\int_Dy^2\,dA=\pi/4$, $\int_Dx\,dA=0$ and $\int_DdA=\pi$. Consequently the <volume integral> is
$$
\int_0^h\left(\frac\pi4-2\pi z\right)dz=\boxed{\frac{\pi h}{4}-\pi h^2}.
$$
For an independent <surface integral> calculation, the relevant <cross product> is
$$
\mathbf F=\nabla a\times\mathbf B=(x^2+2yz^2,\ 2yz-xy,\ -y^2z-xy).
$$
The bottom face contributes $\int_Dxy\,dA=0$ and the top contributes $\int_D(-hy^2-xy)\,dA=-\pi h/4$. On the cylindrical side, put $x=\cos\theta$, $y=\sin\theta$, $\hat{\mathbf n}=(x,y,0)$ and $dS=d\theta\,dz$. The <flux> density is $x^3+2xyz^2+2y^2z-xy^2$. All terms except $2y^2z$ integrate to zero because they are <odd functions> of $x$ or $y$, giving
$$
\int_{\mathrm{side}}\mathbf F\cdot\hat{\mathbf n}\,dS=\int_0^h2\pi z\,dz=\pi h^2.
$$
Thus minus the total outward <flux> is $-(\pi h^2-\pi h/4)$, in agreement with the <volume integral>.