We were having a discussion about the orthogonality of ROMS grids and the effect of loss of such orthogonality on conservation properties of ROMS.
To be more specific, we could not immediately agree on whether grids need to be discretely orthogonal or merely up to discretization errors.
Does anyone have a discrete expression for the required orthogonality of a ROMS grid (in spherical coordinates) and does the violation of such expression break any discrete conservation property?
Jeroen Molemaker, UCLA
Orthogonal grids and conservation laws
-
arango
- Site Admin
- Posts: 1364
- Joined: Wed Feb 26, 2003 4:41 pm
- Location: DMCS, Rutgers University
- Contact:
I pretty much avoid curvilinear grids that follow a particular coastline because of the orthogonality problem and its implication on the conservation properties of the model. I prefer uniform grids with land-sea masking.
Now, it is possible to build orthogonal curvilinear grids with analytical expressions. In ROMS we can have Polar Coordinates (see ana_grid.h for the LAB_CANYON test case) or Spherical Coordinates (see ana_grid.h for the BENCHMARCK test cases). Both applications are fully orthogonal but discrete. I have not evaluated the conservation properties for these test cases.
Is the variable grid-spacing affecting the conservation properties? This is a good question for which I don't have an answer from the discrete point of view. Are the additional curvilinear terms d(1/n)/d(xi) and d(1/m)/d(eta) used in the horizontal advection terms preserving the discrete conservation properties? These terms are also known as dndx(i,j) and dmde(i,j). They are computed in a discrete (spatial derivatives) stencil where m and n are the grid (xi,eta) distances in meters. In Polar and Spherical coordinates we have continuous equations for these terms but they are still discretized.
In other Spherical (geographical) applications, for example, we can have a Cartesian representation of a Mercator projection grid. This grid is orthogonal. However, the curvilinear metrics n and m are discretized and saved as double-precision quantities in the grid NetCDF to avoid roundoff and achieve a discrete conservation.
Things can get more complicated because of the terrain-following coordinates nature of ROMS and the nonorthogonality of the omega vertical velocity. The conservation of total energy is very tricky to derive in the continuous equations and unlikely to achieve in the discretized equations. Applications with open boundary conditions aggravate the conservation problem further.
Anyway, these are excellent questions that we would like to explore. We need to design a good test case to explore the orthogonality and conservation properties. There are several nonorthogonal models out there that make me uneasy. The conservation properties are essential in long simulations and climate research. Models like MOM have discrete conservation integral constraints which are needed in climate simulations over hundreds of years.
Now, it is possible to build orthogonal curvilinear grids with analytical expressions. In ROMS we can have Polar Coordinates (see ana_grid.h for the LAB_CANYON test case) or Spherical Coordinates (see ana_grid.h for the BENCHMARCK test cases). Both applications are fully orthogonal but discrete. I have not evaluated the conservation properties for these test cases.
Is the variable grid-spacing affecting the conservation properties? This is a good question for which I don't have an answer from the discrete point of view. Are the additional curvilinear terms d(1/n)/d(xi) and d(1/m)/d(eta) used in the horizontal advection terms preserving the discrete conservation properties? These terms are also known as dndx(i,j) and dmde(i,j). They are computed in a discrete (spatial derivatives) stencil where m and n are the grid (xi,eta) distances in meters. In Polar and Spherical coordinates we have continuous equations for these terms but they are still discretized.
In other Spherical (geographical) applications, for example, we can have a Cartesian representation of a Mercator projection grid. This grid is orthogonal. However, the curvilinear metrics n and m are discretized and saved as double-precision quantities in the grid NetCDF to avoid roundoff and achieve a discrete conservation.
Things can get more complicated because of the terrain-following coordinates nature of ROMS and the nonorthogonality of the omega vertical velocity. The conservation of total energy is very tricky to derive in the continuous equations and unlikely to achieve in the discretized equations. Applications with open boundary conditions aggravate the conservation problem further.
Anyway, these are excellent questions that we would like to explore. We need to design a good test case to explore the orthogonality and conservation properties. There are several nonorthogonal models out there that make me uneasy. The conservation properties are essential in long simulations and climate research. Models like MOM have discrete conservation integral constraints which are needed in climate simulations over hundreds of years.
Last edited by arango on Wed Feb 13, 2008 6:47 am, edited 6 times in total.
just throwing in my 2 cents:
dont confuse conservation properties with accuracy.
I think the model will still completely conserve tracer mass, with or without the dmde and dndx curv terms. However, the accuracy could suffer if the grid curvs and the curvilinear terms are not a 'good' representation of the true structure of the grid.
dont confuse conservation properties with accuracy.
I think the model will still completely conserve tracer mass, with or without the dmde and dndx curv terms. However, the accuracy could suffer if the grid curvs and the curvilinear terms are not a 'good' representation of the true structure of the grid.
I agree with John to some degree. If you are careful when implementing conservation equations such that you use exactly the same formulation for both the flux leaving a cell and entering a neighboring cell (or even better: compute the flux once and use it twice, once with + and once with - sign) then you would have conservation of the considered quantity no matter the grid characteristics. For mass conservation (of water, sediment or tracer) this is relatively straightforward. Whether this is true for ROMS I don't know (still too limited knowledge of the code) but I would expect it.
However, for momentum conservation (and associated energy conservation) this is more difficult due to the staggering and the importance of curvature terms.
Additionally there is the accuracy issue that John already mentioned: forcing terms (e.g. water level slopes, concentration/density gradients, Coriolis) are generally assumed to be perpendicular for the two grid directions. If the grid is not orthogonal you must be careful not to mess up the forcing; even though the conservation equations may then be satisfied your answer could be completely wrong. In a vertical plane you have similar problems distinguishing horizontal gradients and vertical gradients in a terrain following coordinate system as Arango referred to.
However, for momentum conservation (and associated energy conservation) this is more difficult due to the staggering and the importance of curvature terms.
Additionally there is the accuracy issue that John already mentioned: forcing terms (e.g. water level slopes, concentration/density gradients, Coriolis) are generally assumed to be perpendicular for the two grid directions. If the grid is not orthogonal you must be careful not to mess up the forcing; even though the conservation equations may then be satisfied your answer could be completely wrong. In a vertical plane you have similar problems distinguishing horizontal gradients and vertical gradients in a terrain following coordinate system as Arango referred to.
-
sak007
This is an interesting question. When generating a grid for D'Entrecasteaux Channel (http://www.marine.csiro.au/~sak007/dent.pdf) for CSIRO's in-house model MECO (now SHOC), I have initially generated it in lat-lon rather than using an orthogonal projection like Mercator. In regions other than close to the equator this results in systematic deviations from orthogonality -- small squares in lat/lon become diamonds in physical space. Some time later, after correcting the error, I compared the output of the two otherwise identical models (different only in projection used for the grid generation) and could not notice any visual difference in output fields.
In my view, changes in horizontal grid scale within a grid to get a "focus" on certain area may be of more concern in regard to the model precision than small deviations from orthogonality caused by the finite size of the grid cells.
Take care,
Pavel Sakov
In my view, changes in horizontal grid scale within a grid to get a "focus" on certain area may be of more concern in regard to the model precision than small deviations from orthogonality caused by the finite size of the grid cells.
Take care,
Pavel Sakov
The circle test problem can be used to get some understanding of this. I ran it with grids that are Cartesian with masking and grids that are circles with orthogonality errors in the "corners". You get errors from the staircase in the first case, errors from the funky resolution in the second case. The nature of orthogonal grids is that the resolution focuses on the peninsulas jutting into the domain and away from the bays. For the circle, you get essentially four bays and four peninsula-like corners.
At all resolutions, the staircases generated noise on the grid scale. Convergence to the right answer was linear in dx, as I recall. For the curvy grids, the errors are smaller at coarse resolution, but adding gridpoints wasn't very good at reducing the errors since it is very hard to get that resolution into the "bays" where it is needed, especially in the cross-shelf direction. Doing so would also cause the corners to get an incredibly small dx and therefore dt.
At all resolutions, the staircases generated noise on the grid scale. Convergence to the right answer was linear in dx, as I recall. For the curvy grids, the errors are smaller at coarse resolution, but adding gridpoints wasn't very good at reducing the errors since it is very hard to get that resolution into the "bays" where it is needed, especially in the cross-shelf direction. Doing so would also cause the corners to get an incredibly small dx and therefore dt.