Poor scaling for runs with stochastic perturbations (in RRFSv1 system)

[MatthewPyle-NOAA]

While testing the RRFS ensemble system using more nodes, we noticed that the timing improvements with more resources wasn't very good, and seemed to level off pretty significantly for higher node counts. As background, we've been running our experiments with 52 nodes, but plan to run in production with about 92 nodes. Out of curiosity we increased the node count first from 92 to 110 nodes (a 19.5% increase) but the timing only decreased by about 3.5%. Our deterministic runs not using perturbations scale quite a bit better. Wasn't sure if this is a known feature - we are running on a 3950 x 2700, 3 km grid spacing grid, and due to stability needs have to regenerate the stochastic pattern generation quite frequently (lndpint=180 - 5 time steps; sppint=36 - every time step)

WCOSS2 nodes	wall clock (s)
52	11525
92	8463 (from 52 to 92: ~77% more nodes, about 26.6% faster)
102	8251
110	8163 (from 92 to 110: ~19% more nodes, about 3.5% faster)

Tagging @JiliDong-NOAA as well.

[MatthewPyle-NOAA]

@JiliDong-NOAA Could you describe the threading work you did on the stochastic physics?

[JiliDong-NOAA]

Two things I have tried to improve stochastically perturbed parameterization SPP efficiency are as follows:

1. add OMP directives to subroutines getnoise and patterngenerator_advance
2. limit soil moisture and temperature (smc/stc) random pattern generation only to initial time

The changes are here:

JiliDong-NOAA@8f9abed

[MatthewPyle-NOAA]

@JeffBeck-NOAA Your name came up in a meeting I had yesterday where this item was discussed briefly. Not sure if you have any thoughts or suggestions on improving the speed of perturbed REFS members.

One item that never made a lot of sense to me is why our runs are much more unstable when creating the stochastic patterns less frequently. Once you'd added that nice generality to the code we had hopes of generating them quite infrequently (hourly?) with a huge improvement in model speed, but we've gradually had to increase the frequency of the generation almost back to every timestep for RRFS/REFS.

Thanks for any thoughts you might have!

[JeffBeck-NOAA]

@MatthewPyle-NOAA, thanks for tagging me. Unfortunately, we have never tested scaling characteristics of the stochastic physics pattern generator in the FV3-LAM with SPP or the ad-hoc schemes at CAM scales. Since perturbations are applied each time step regardless of the SPPINT setting in the namelist, and assuming you kept SPPINT constant across your test runs, the lack of effective scaling should be related to the calculation of the pattern, not the application of the perturbations. I'm tagging @pjpegion to see if he has any insight on possible reasons why.

Aside from mitigating the scaling issue through optimization, my next recommendation would be to increase the value of SPPINT, which I see you've done, but this has led to instability. Freezing the pattern through use of SPPINT will result in identical perturbations being applied at each time step for every point in the domain until the pattern is updated. If SPPINT is large, this means that areas of the random pattern with large positive or negative perturbations will be applied continuously for a long time over a given area. I could definitely see how this could result in instability in the model. At least at CAM scales, ideally, SPPINT shouldn't be used at all, since it interrupts the normal evolution of the stochastic physics random pattern, which is based on scale-specific magnitudes and decorrelation lengths and times set in the namelist. This evolution keeps large perturbations from occurring in the same locations for long periods of time. However, for computational purposes, the *INT options can be useful.

In order to make this work (aside from refactoring the pattern generator), you will likely need to reduce the magnitudes or decorrelation lengths and times that are being used in order to avoid large perturbations persistently being applied at the same location with a large SPPINT. Again, that's my best guess as to why the runs are becoming unstable. That said, smaller perturbations will reduce stochasticity, and reduced decorrelation lengths/times may change the result in unknown ways (although our previous tests with WRF didn't show much sensitivity to decorrelation changes when we ran tests at 3-km).

Finally, one last thought: stochastic physics are only as stable as the underlying physics. If the baseline physics are pushing the envelope for some reason, which I believe is what happened previously in REFS with SPPT, it will be very difficult to keep things stable. Is there (or have there been) any last minute physics changes that may be causing instabilities that could be targeted for mitigation?

I hope this helps!

[MatthewPyle-NOAA]

Thanks @JeffBeck-NOAA for all of your thoughts! The other sticky aspect is that our science is frozen at this point, and our ensemble spread was on the low side of desired during the evaluation, so a significant tuning down the magnitude of perturbations likely isn't possible at this stage. But I've been pondering if a slightly larger SPPINT could be combined with some sort of decaying amplitude perturbations between pattern regenerations. But this likely wouldn't be a game changer given the other constraints.

[dkokron]

I working this issue on the WCOSS2 test cluster.

There is a loop in subroutine get_random_pattern_sfc (file ufs-weather-model/stochastic_physics/get_stochy_pattern.F90) that has either 5 or 7 iterations.

There is a call to mp_reduce_sum(workg,lonf,latg) inside that loop which just does a MPI_Allreduce with the MPI_SUM operation. The size of the workg array appears to be fixed and represent a global field.
allocate(workg(lonf,latg))

I found the following in the output listing. That output comes from either compns_stochy() or compns_stochy_ocn().

gaussian grid not set, defining here 1608 800

MPI_Allreduce() with a fixed size global array argument will slow down with increasing process count.

What is the science reason for combining workg from all ranks?

[JeffBeck-NOAA]

@MatthewPyle-NOAA, I had wondered about how much is truly possible to modify given that the science evaluation has already taken place. A slightly larger SPPINT should have a minimal impact on results, but decaying perturbations may not. It's tough when there is an issue of instability involved.

@dkokron, the spherical harmonics used to create the perturbation pattern are global, regardless of whether a regional domain is being used, but I can't speak to the MPI question. @pjpegion may be able to help.

[pjpegion]

I do not see a way around the MPI calls. The spherical harmonics are indeed defined globally, and the spectral transforms are distributed across all of the processors, which each processor getting latitude bands. All of this was inherited from the spectral GFS. The next step is to gather all of the individual latitude bands on each processor so each processor can interpolate from the global gaussian grid to its sets of model grid points.
Perhaps only the root processor can gather all of the global gaussian grid points, the interpolate to the full model grid, then distribute each processor its set of points, but that would require the stochastic physics package to know the mpi decomposition, and I don't know if this would be much faster.

[dkokron]

Is there anything special about the random numbers that come out of that spherical harmonics code? Could we generate a random sequence with similar statistics by calling mersenne twister on the native grid?

[pjpegion]

The random patterns need to be correlated in space and time. A few years ago, I tried computing random numbers only correlated in time at each grid point, then used the model's diffusion operators to smooth the noise what is white in space to create spatial correlation, but it turned out to be difficult to control the amplitude of the perturbations, and wound up being slow for a high resolution model, when you needed large scale spatial patterns.

[dkokron]

A deeper analysis of the glolal() and workg arrays.
glolal is declared (lonf,gis_stochy%lats_node_a) with gis_stochy%lats_node_a being 0, 1 or 2 depending on rank. of the 1792 ranks in this job, 994 have gis_stochy%lats_node_a set to zero, 796 have it set to 1 and 2 of the 1792 have it set to 2.

Presumably, this means that less than half of the ranks contribute non-zeros to the final workg array. Furthermore, even when a rank does contribute non-zeros, it's only one or two latitude bands. The vast majority of the local workg is zeros! This means the MPI_Allreduce is mostly just summing zeros. What a horrible waste!

What you want is an MPI_AllGatherv() not MPI_Allreduce(). I'll see if I can code that up and test.

[dkokron]

I provided an approach using MPI_AllGatherv() to Matt for testing.

[dkokron]

Does anyone know if the pattern generation code is thread safe? Can the lndp_loop patterns be generated in parallel?

[JiliDong-NOAA]

Does anyone know if the pattern generation code is thread safe? Can the lndp_loop patterns be generated in parallel?

I believe so. The same as npatterns in subroutine get_random_pattern_spp. This could be potentially optimized with OMP.

The pattern generation is also independent from the forecast itself. Potentially, we could use dedicated nodes for pattern generation ahead of the integration?