Arakawa & Konor 1996 with mesh distortions

A mountain and a basic terrain-following mesh

I've put a small mountain at the bottom of the domain (schaerCos, peakHeight=250m, halfWidth=125km, wavelength=40km). I've also fixed the fvSchemes to use midPoint interpolation and gradients and, importantly, uncorrected snGrad schemes. Results after two days, Lorenz (left), Charney--Phillips (right):

Again the computational mode is present only with Lorenz. But there are some noticeable numerical errors near the ground with Charney--Phillips that will need investigating. Also, a general cold bias is present in the Charney--Phillips model only.

And the same results but high-resolution (hence no dz-scale perturbation)

I suspect that the warming at the bottom is another reinforcing feedback when grad(theta) is extrapolated to the lower boundary. Here's a zoomed in version showing theta and thetaf:

Yup! Changing top/bottom BCs to be fixedValue rather than calculated gives a better result (although maybe this is a cheat?)

Rectangular domain with mesh distortions

We wonder whether this error is caused by lower boundary condition problems or by the non-orthogonal treatment of theta, or both. To try to hone in on the problem, we'll do the same test with no mountain but use edgeGrading to distort the mesh. In this test, the theta perturbation is not placed on model layers but in Cartesian coordinates, so the perturbation may span more than two model layers.

With Charney--Phillips, a hot blob develops in the middle of the atmosphere. The same blob is present though weaker in amplitude with a less distorted mesh.

Advection-only tests

Perhaps the errors seen in the animation above are due to advection errors? Certainly we've seen bad results from advection-only experiments using the theta advection scheme from exnerFoamCPinterpGrad. I've included the necessary schemes to plug into advectiveFoamF. There is a domain decomposition that causes the tracer to virtually disappear when it is transported from the left-hand block to the right-hand block. So I'm running in serial until that gets fixed.

T_diff errors on orthogonal mesh

Non-orthogonal mesh results are similar

This makes me question whether advection is really the cause of the problem seen in the dynamics test on edge graded meshes.

Trying the same test with div(phi,T) Gauss linear with advectionFoam (Lorenz staggering):

Orthogonal mesh

Non-orthogonal mesh

Errors are much smaller with advectionFoam compared to advectiveFoamF.

Energy time series

Charney--Phillips with edge grading

Lorenz with edge grading

Charney--Phillips orthogonal

Lorenz orthogonal

sum(T_diff) time series (now found to be incorrect!)

This was because of a bug in the cell-to-face interpolation: mag(Tf*g.unitFaceNormal()) (incorrect) is not the same as Tf*mag(g.unitFaceNormal()) (correct) when Tf is negative. The new results are with interpGrad and fvcReconstructCP. The updated discussion is on a separate page.

Charney--Phillips

Lorenz

These (incorrect) plots suggested to me that we do indeed require conservative advection for Charney--Phillips. In this advection test the advective form interpGrad scheme is warming up over time.