Earth Science Stack Exchange is a question and answer site for those interested in the geology, meteorology, oceanography, and environmental sciences. It only takes a minute to sign up.
Sign up to join this community
I am working with a coastal hydrodynamic model in the MIKE3 modelling package. MIKE3 uses the finite volume method to solve the shallow water equations with explicit timestepping, over a triangular mesh.
My model frequently "blows up" (numerically diverges) in a few cells around a headland, at times of peak flow. What steps can I take to improve the numerical stability, either through modifications to the mesh or other measures?
Almost all finite methods that use forward time models adhere to the Courant-Friedrichs-Lewy law which calculates a courant number and compares it to a $C_{max}$, which is what determines stability, so for 2-D:
$$C = \frac {V_xdt}{dx} + \frac{V_ydt}{dy} \ge C_{max}$$
Where $C$ is the courant number, $V_i$ is the velocity in subscripted direction, $dx$ or $dy$ is the length interval in specified direction and $dt$ is the specified time step. If you are getting numerical instability, increasing the mesh spacing will lower the courant number, but then you might be missing the dynamics you want to capture. If you decrease the time step, it will take longer for the dynamics to develop computationally.
One way to solve this computational problem is to code in a dynamic time step: Keep calculating the $C$ and $C_{max}$. If $C > C_{max}$, then have the model lower the step by half. Check again, then reset for the next step so you can have larger timesteps for "less dynamic" times. I have implemented this many times in other models and it has worked well.
Another possibility might be that you have high aspect ratios in your mesh: if your $x$ direction is 10 km and your $y$ is 1 km, a $10:1$ aspect ratio is very high. In your case, since the model package you use are triangles, you want the angles in the triangle grid to be as acute as possible. Obtuse angles would mean high aspect ration, in general.
Like llmari, I am not familiar with hydrodynamic models, so I do not do specifics. My expertise are in both mantle convection and plate flexure models. I hope this helps.
The usual panacea for numeric stability problems is to decrease the timestep. It doesn't always help (and, indeed, in some cases can even make the instability worse), but it's generally the first thing to try.
Of course, a shorter timestep costs more computer time. Another approach to try, if reducing the timestep doesn't help, or leads to unacceptably poor performance, is to either refine or, somewhat paradoxically, coarsen the mesh in the problematic regions. (Which one of these is more likely to help depends on the exact cause of the instability.) If this alone doesn't help either, try adjusting the mesh and decreasing the timestep.
Ultimately, the best way to fix such problems is usually to switch to a more stable timestepping method. As a general rule of thumb, higher-order methods are more stable than lower-order methods (up to some point, anyway), and implicit methods are more stable than explicit. In both cases, stability comes at the cost of extra computational effort, but it's still generally cheaper than just decreasing the timestep enough to achieve comparable stability.
Unfortunately, I'm not personally familiar with the MIKE3 software, and have no idea what timestepping algorithm(s) it uses, or whether there's any way for the end user to adjust the algorithm.
