Question about energy in the shallow water equations on a staggered grid

I think this is a question about the energy conservation of a numerical integrator. I'm studying the linearized 1D shallow water equations in python - for reference, here they are:

$$\frac{\partial u}{\partial t} +U\frac{\partial u}{\partial x}+g\frac{\partial h}{\partial x} = 0\\ \frac{\partial h}{\partial t} +U\frac{\partial h}{\partial x}+H\frac{\partial u}{\partial x} = 0\\$$

where $$u$$ is the speed of the water, $$h$$ is the height of a water column, $$U$$ is the average speed, $$H$$ is the average height and $$g$$ is the gravitational constant (9.81 m/s²). I'm trying to solve these equations with a centered space discretization and a leapfrog timestep (with a forward first step) on a staggered grid. The picture below shows an example of a non-staggered and a staggered grid. The integration works and the results are good. Now, I'm investigating the energy of the system. My question is if that seems like a reasonable thing to do. Does it make sense to define the energy of the "waves" as $$E = \frac{1}{2}u^2 + gh$$? For every timestep I calculated the kinetic and potential energy at every position and integrated them over the spatial domain with numpy.trapz() to get an average so that I could plot the total energy over time: but I don't really understand what I am looking at. Shouldn't the energy be sort of conserved? Isn't leapfrog a symplectic integrator? Or does that not mean anything in this context?

• Simply use a FV method for such kind of problems. FD methods are not inherently conservative, nor are they suited to equations which might produce discontinuities. Dec 8 '21 at 22:27
• As long as these are the linear shallow water equations shocks will not be produced unless they were in the initial conditions, but it is a good point. While FV methods would deal with energy conservation there are a number of FD schemes especially designed to maintain things like energy. These are mostly predicated on the idea that the energy is the flux of the originally non-linear equations. The reason why you are seeing $\mathcal{O}(\Delta x)$ (?) variations is well explained by the answer below and the fact that you’re scheme is almost certainly not conserving the discrete version. Dec 10 '21 at 15:25

I'm going to work on the assumption that the disturbances aren't large enough to create shock waves, which are a whole other can of worms.

Your energy functional is close but not quite right -- the Hamiltonian for the shallow water equations is

$$H = \frac{1}{2}\int_\Omega\left(h|u|^2 + gh^2\right)dx,$$

although the Poisson brackets are a little more involved. You're also correct to identify symplectic integrators as the tool to use if you care about the Hamiltonian structure of the problem. But it's important to note that a symplectic integrator conserves the energy of some perturbed Hamiltonian

$$H' = H + \delta t\cdot\delta H + \mathscr{O}(\delta t^2),$$

where very often we can read off what $$\delta H$$ is in terms of a few Poisson brackets using the Baker-Campbell-Hausdorff formula. In other words, some energy drift is to be expected. Some time schemes are guaranteed to exactly conserve all linear and quadratic invariants, but the $$h|u|^2$$ term isn't quadratic. That said, the property of exactly sampling from the trajectory of the modified Hamiltonian $$H'$$ is incredibly useful -- you get a lot of properties for free, like the structure of the system near equilibria, bounded orbits, phase space volume conservation, etc.

I haven't worked through the details of symplectic discretization of the shallow water equations specifically, so take this with a grain of salt. There are several explicit schemes that are only symplectic for simple Hamiltonians that can be written as $$H = K(p) + U(q)$$ with the usual Poisson bracket. The Hamiltonian and the bracket for SWE are neither of those. I think that the implicit midpoint rule

$$\frac{z_{n + 1} - z_n}{\delta t} = f\left(\frac{z_{n + 1} + z_n}{2}\right)$$

for the ODE $$\dot z = f(z)$$ is symplectic even in unusual cases. Now you've got an implicit equation to solve, but nothing's free. So it may be worth trying other timestepping schemes.

Finally, upwinding interacts with the Hamiltonian structure of the problem in ways that are a little bit subtle. Although they took more of a finite element approach, this paper worked through how to do energy-conserving discretization of SWE with upwinding.

Based on the plot you're showing, you should be asking yourself whether the energy drift that you've found indicates an error or failure to preserve the Hamiltonian structure, or merely the expected drift from integrating a slightly perturbed problem. While the Hamiltonian that you wrote down isn't completely correct, the plot you showed indicated an energy drift on the order of 0.035 and a total of around 10, so the drift is less than 1% of the total. That doesn't set off any red flags for me, but if you want to try a few things to debug this, you can:

1. Run the simulation again at resolution $$\delta x / 2$$, $$\delta t / 2$$ and see what the energy drift looks like. Any convergent timestepping scheme (symplectic or not) should halve the drift, so this won't tell you everything you need to know.
2. Run the simulation for a very long time and see what the energy drift looks like. To get an idea of what "long" means, you get a speed scale from the wave speed $$c = \sqrt{gh}$$ and a length scale from the diameter $$L$$ of the domain, so several multiples of $$L / c$$ is a good ballpark. If the solution explodes, your scheme probably isn't symplectic.

I would agree with the suggestion in the comments to just use finite volume methods if the problem you had related to capturing shock waves, but that doesn't appear to be the case.