I have a numerical scheme to solve an hyperbolic system of equations of this type (this a more simple version just for clarity):
$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{\partial H}{\partial x}$$ $$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial x}=\frac{\partial H}{\partial y}$$ $$\frac{\partial H}{\partial t}+ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial x}=0$$
I have a steady state analytical solution, and I am checking the order of accuracy of the scheme. It is a direction splitting scheme, and all the terms should be second order accurate. I used 4 grids, with successive refinements $i=1\ldots 4$ (cartesian uniform grids with grid size $\Delta x_i=\Delta x_{i-1}/2$ and $\Delta x=\Delta y$). I saw that the order of accuracy is somehow around 1.5, therefore I decided to replace each derivative with the correspondent analytical term until I found the term responsible for the loss of accuracy. The culprits were: $\frac{\partial H}{\partial x}$ and $\frac{\partial H}{\partial y}$. Now my questions are:
1) Am I doing right by using a fixed CFL number, in the sense that every time I half the grid size I also half the time step? What if I use a fixed dt for all the grids? Note that the scheme is implicit and unconditionally stable.
2) Im using a staggered grid, i.e. H is cell centered and u and v are on the center of the cell edges. The terms $\frac{\partial H}{\partial x}$ and $\frac{\partial H}{\partial y}$ are discretized by second order central differences. When I print the error of $\frac{\partial H}{\partial x}$ after a single time step starting from the exact solution I see that it converges perfectly with second order accuracy (it is 4 times smaller at each refinement). However when I start from a second-order accurate approximation of H, $\frac{\partial H}{\partial x}$ has first order rate of convergence. I am trying to figure this out and it might make sense since $\frac{\partial H}{\partial x}= \frac{H_{m+1}-H_m}{\Delta x} + O(\Delta x^2)$ that is a second order accurate expression at the edge $m+1/2$ if $H_{m+1}$ and $H_m$ are exact. However, if I start from a second order accurate approximation of $H_{m+1}$ and $H_m$: $$H_{m+1} = H^{exact}_{m+1}+ O(\Delta x^2)$$ $$H_{m} = H^{exact}_{m}+ O(\Delta x^2)$$ (as it would be the second order accurate solution of the system) I get $$\frac{\partial H}{\partial x}= \frac{H^{exact}_{m+1}-H^{exact}_m}{\Delta x} + \frac{O(\Delta x^2)}{\Delta x} + O(\Delta x^2)$$ that is first order, and thats what I get (first order) if I print the numeric value of $\frac{\partial H}{\partial x}$ for the 4 different grids. However, if I print $H$ for the first half time step it is second order (it is the second half of the splitting that cause the order of accuracy). Am I doing something wrong here?