v2.0.0

This is a pre-release version of SOMAR v2.0.0. It is the code that was used to run test cases for an upcoming paper.

Summary of changes

SOMAR v1.0.0 used cell-averaged state variables (for both the velocity and scalars), an approximate projection scheme, and a time-stepping method based on the works of Bell, Colella, and Glaz. Nonlinear advective fluxes were computed using Colella and Woodward's piecewise parabolic method and Saltzman's unsplit 3D assembly scheme. Also, topographies were handled using a non-orthogonal sigma-coordinate system. These choices raised a number of concerns.

1. Incompressible flows are very naturally modeled using staggered meshes. Cell-centered velocity discretizations rely on approximate projections that have poor convergence and stability properties and fail to remove high-wavenumber irrotational flow components, both of which lead to instabilities when the flow is marginally resolved.

2. The time integration scheme of Bell, Collela, and Glaz has an $\mathcal{O}(\Delta t^2)$ local error in the velocity. After proceeding through $1/\Delta t$ timesteps, diffusive errors appear, interfering with the operation of the turbulence model.

3. The piecewise parabolic method is an excellent method for modeling the nonlinear advection terms present in the Navier-Stokes equations, but we found that upwinding methods are not well-trusted by the community of turbulence modelers. Much of this distrust is warranted, as upwinding schemes are typically associated with limiters, artificial viscosity, and implicit LES schemes.

4. Non-orthogonal coordinate systems have a number of disadvantages. With off-diagonal metric terms present, Laplacians have elevated condition numbers and more terms to compute. Poisson solvers are thus made considerably more expensive than their orthogonal counterparts -- about 5-8 times more expensive in our experience. Also, strongly non-orthogonal coordinates have been shown to incur first-order spatial errors. Worse of all, sigma coordinates are only useful when modeling simple, idealized terrain.

To address these issues, SOMARv2 uses Arakawa-C grids for the velocity and pressure with standard finite difference stencils. The advective forcing terms are computed in conservative form with fluxes generated using unbiased, fourth-order averaging stencils to re-center the velocity components where necessary. Time integration is accomplished with a globally third-order Runge-Kutta method with excellent stability properties and error control. (We will soon incorporate an immersed boundary method to model complex bathymetry without the need for a non-orthogonal coordinate system.)

All of these changes led to a more accurate fluid solver, whose AMR capabilities are significantly more robust. Second-order AMR solvers suffer from so-called ''mesh-imprinting errors'' when features that require fine-level discretization slip past the coarse-fine interface (CFI) and into the coarser level. These errors reveal themselves as large, first-order errors in the cells that abut the CFI. In SOMAR v1.0.0, subduing these errors required a great deal of care when choosing refinement criteria. SOMAR v2.0.0's use of Arakawa-C grids and higher-order temporal accuracy eliminates these errors, removing much of the ``art'' of choosing appropriate fine-level grids.