Optimize `niter`

[jbusecke]

One of the question that came up in our sprint planning was what value we should choose for the iteration parameter niter in our work.

The values used in aerobulk vary: I found niter=10 in the examples while the default value is apparently 4.

I will start to investigate how the performance of the algorithm scales with niter, but we will also have to discuss how we judge whether the values have 'converged'. I suspect that this will to some degree depend on the location and combination.

Do people have suggestions how to do this?

[TomNicholas]

how we judge whether the values have 'converged'

I mean I would just start with plotting the globally integrated flux as a function of niter. Could also plot some measure of the difference between niter=n and niter=n+1, against against niter.

[jbusecke]

TLDR; niter as 6 to 7 (depending on the algo) is a very solid choice, giving us converged values for 95% of the values represented in our CM2.6 example.

Nice! That is very similar to what I cooked up. I did however started with a single grid point first:

Synthetic example for single grid point

I am using our [canonical test values]( sst = np.full(shape, 290.0, order=order) t_zt = np.full(shape, 280.0, order=order) hum_zt = np.full(shape, 0.001, order=order) u_zu = np.full(shape, 1.0, order=order) v_zu = np.full(shape, -1.0, order=order) slp = np.full(shape, 101000.0, order=order) flux_noskin(sst, t_zt, hum_zt, u_zu, v_zu, slp, algo=algo)) and compute them for each algorithm and for niter in [1,...,8]. I then define the 'relative iteration error':

$\epsilon(niter) = |(aero(niter) - aero(niter-1))/aero(8)|$ for niter in [2, ..., 8]

We would like this quantity to be small, which means the output for niter is very similar to the output of (niter-1).

Each panel shows a different output variable as function of niter. The algorithms are colorcoded). The vertical lines indicate the value where $episilon(niter) < 0.2 $%.

Depending on the algo, niter values from 3-4 seem to work well for this example. $episilon(niter) < 0.2 $% seems to work well as criterion for 'convergence' but is somewhat arbitrary.

It is interesting to see how much these algos differ in the converged values! Might be relevant to #26.

This conclusion could however be heavily biased due to the values chosen. So next I tried to replicate this pointwise analysis for a realistic model dataset (CM2.6).

Full model output

I am currently restricting the domain to the Northern and parts of the South Pacific and two daily snapshot half a year apart, to reduce the amount of memory used. I believe this provides a wide range of scenarios at each grid cell, and would be surprised if the results change drastically for a global domain and more time steps, but I should still recompute these results, once we have a working xarray wrapper.

I basically calculated the relative iteration error like above for every point in space and time. To summarize the relative iteration error, I opted for a high quartile along x,y, and time. Basically what this is showing is when 'most' of the grid cells have converged to a low relative iteration error. This is a relatively strict criterion, but as you can see in the notebook the relative iteration error is actually very spatially dependent.

Each panel shows a different output variable. Lines are color coded for algo, and the values shown indicate the 95th percentile of the relative iteration error. Dashed line shows same value used in the synthetic example

The results suggest that we should chose a slightly higher niter value than the synthetic example suggested above. We can probably settle on 6 (or 7/8 if we are using 'andreas').

Very keen on comments on this methodology. I think ultimately this would be a section of the Supp Materials for any paper resulting from this work.

[rabernat]

This is great work! Let's set niter=6 as the default in the code.

[paigem]

Closing since this seems to be solved for now!