Error at 256GPU/64Nodes on Lassen

[MTCam]

This is a weak scaling (periodic 3D box, p=3) case to see how far we can push running real simulations on Lassen GPUs. It has 48k elements/rank and exercises most of the components used in our prediction. After running successfully for all mesh sizes (nelem, ngpus)=[(48k, 1), (96k, 2), (192k, 4), (384k, 8), (768k, 16), (1.5M, 32), (3M, 64) , (6M, 128)], it fails at (nelem, ngpus) = (12M, 256).

*Note: To reproduce, make sure to request at least 5 hours on 64 Lassen batch nodes.

Instructions for reproducing on Lassen:

Branch: mirgecom@production
Driver: examples/combozzle-mpi.py
Setup:

• Set the following parameters to control the mesh size, either hard-code or with a yaml file (e.g. setup.yaml)

x_scale: 4
y_scale: 1
z_scale: 1
weak_scale: 20

• Set the following environment:

export PYOPENCL_CTX="port:tesla"
export XDG_CACHE_HOME="/tmp/$USER/xdg-scratch"
export POCL_CACHE_DIR_ROOT="/tmp/$USER/pocl-cache"
export LOOPY_NO_CACHE=1
export CUDA_CACHE_DISABLE=1

• Run the case:

jsrun -a 1 -g 1 -n 256 bash -c 'POCL_CACHE_DIR=$POCL_CACHE_DIR_ROOT/$$ python -u -O -m mpi4py ./combozzle-mpi.py -i setup.yaml --lazy'

The case does have an unexplained hang at array context creation time which grows with the number of ranks. For this case, it will hang for about 3 hours, and then compile for about 10 minutes before finally ending with this error:

  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/mirgecom/mirgecom/integrators/lsrk.py", line 66, in euler_step
    return lsrk_step(EulerCoefs, state, t, dt, rhs)
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/mirgecom/mirgecom/integrators/lsrk.py", line 53, in lsrk_step
    k = coefs.A[i]*k + dt*rhs(t + coefs.C[i]*dt, state)
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/arraycontext/arraycontext/impl/pytato/compile.py", line 365, in __call__
    compiled_func = self._dag_to_compiled_func(
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/grudge/grudge/array_context.py", line 286, in _dag_to_compiled_func
    ) = self._dag_to_transformed_pytato_prg(
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/arraycontext/arraycontext/impl/pytato/compile.py", line 441, in _dag_to_transformed_pytato_prg
    pytato_program = (pytato_program
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/pytato/pytato/target/loopy/__init__.py", line 142, in with_transformed_program
    return self.copy(program=f(self.program))
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/meshmode/meshmode/array_context.py", line 1726, in transform_loopy_program
    raise err
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/meshmode/meshmode/array_context.py", line 1723, in transform_loopy_program
    iel_to_idofs = _get_iel_to_idofs(knl)
  File "/p/gpfs1/mtcampbe/CEESD/Experimental/svm/production/meshmode/meshmode/array_context.py", line 1066, in _get_iel_to_idofs
    raise NotImplementedError(f"Cannot fit loop nest '{insn.within_inames}'"

Gist with the entire batch output file:
gist.github.com/MTCam/f48991f584755b6a8530dd9345dc2de4

[inducer]

Do you have the full exception message? (This is just the traceback leading up to the raise AFAICT.)

[MTCam]

Do you have the full exception message? (This is just the traceback leading up to the raise AFAICT.)

I could not find a raise in the batch job output. This appeared to be the only message out-of-the-ordinary when the job failed. Here's a gist with the entire output:
https://gist.github.com/MTCam/f48991f584755b6a8530dd9345dc2de4

[inducer]

What's happening there is that the transform code is running into a pattern of loop nests that it's not familiar with, see this code in @kaushikcfd's branch:

https://github.com/kaushikcfd/meshmode/blob/c44f1b12b1e43d006808877fb5d1d8131bc1e396/meshmode/array_context.py#L1724-L1733

The message I was looking for is

NotImplementedError: Cannot fit loop nest 'frozenset({'idof_ensm10'})' into known set of loop-nest patterns.

which seems very weird; it's a loop over only intra-element DOFs without a loop over elements surrounding it.

As a near-term workaround, you can just hack that code to fall back to the more general single-grid work-balancing strategy also in the compiled-function case, in the code I linked.

Beyond that, I would like to understand how we wound up with a lonely DOF loop.

[kaushikcfd]

Thanks for the report.

For this case, it will hang for about 3 hours,

It's safe to call this an error. However, I would expect cfd_rhs to have the largest compilation time and the log file says:

transform_loopy_program for 'cfd_rhs': completed (112.52s wall 1.00x CPU)

Am I missing the kernel to blame here?

Beyond that, I would like to understand how we wound up with a lonely DOF loop.

I agree. @MTCam: Does the driver work correctly while running with lower #ranks?

[MTCam]

Thanks for the report.

It's safe to call this an error. However, I would expect cfd_rhs to have the largest compilation time and the log file says:

transform_loopy_program for 'cfd_rhs': completed (112.52s wall 1.00x CPU)

Am I missing the kernel to blame here?

This compile exits abnormally after about 10 minutes. The actual, full compile should take more than an hour, as evidenced by the compile times for the lower rank runs:
(nranks, compile time) = [(32, 733s), (64, 1448s), (128, 2455s)]

Beyond that, I would like to understand how we wound up with a lonely DOF loop.

I agree. @MTCam: Does the driver work correctly while running with lower #ranks?

Yes the driver runs successfully for a number of different meshes up to 6M elements. The difference between this run and the previous working run is that it is a factor of 2 more ranks and a factor of 2 more elements.

[kaushikcfd]

(nranks, compile time) = [(32, 733s), (64, 1448s), (128, 2455s)]

Not sure I'm convinced that should be the case. Synchronizations/filesystem issues come to mind. (@inducer any guesses why this might be happening?)

[inducer]

This timing looks suspiciously exponential in the number of ranks. I do not understand why that would be the case.

@MTCam:

• What exactly is being measured?
• Do you happen to know the average/max number of neighbors per rank? The only well-justified cost growth in compilation would be linear in the max number of neighbors.

[kaushikcfd]

The only well-justified cost growth in compilation would be linear in the max number of neighbors

Are you suggesting this might be because of the increased comm pytato nodes? If so, the lines in the log file of the form:

find_partition: Split 136095 nodes into 3 parts, with [41034, 8180, 88338] nodes in each partition.

would be helpful to validate that hypothesis for the 3 mesh sizes with increasing compile times.

[inducer]

Wow, not sure I'd seen that. 136,095 nodes? That's kinda nuts. There's no reason for there to be that many of them, IMO.

[MTCam]

This timing looks suspiciously exponential in the number of ranks. I do not understand why that would be the case.

@MTCam:

• What exactly is being measured?

The timer starts at the first call to logpyle.tick_before. The time will include the compile of the RHS (and only the RHS) and the time of the first step. The step timings here are O(1s). @matthiasdiener can verify about what is being timed for this.

• Do you happen to know the average/max number of neighbors per rank? The only well-justified cost growth in compilation would be linear in the max number of neighbors.

Not offhand - but I do know it isn't going crazy. I will add a snippet to give us the break-down.

The only well-justified cost growth in compilation would be linear in the max number of neighbors

Are you suggesting this might be because of the increased comm pytato nodes? If so, the lines in the log file of the form:

find_partition: Split 136095 nodes into 3 parts, with [41034, 8180, 88338] nodes in each partition.

would be helpful to validate that hypothesis for the 3 mesh sizes with increasing compile times.

FWIW, the compile times are independent of mesh size. For a given RHS, the compile times increase with increasing numbers of ranks. This is not a new "feature" of distributed lazy. It has always been this way - although now it may be a little worse than before. I am not sure because we were never able to run our cases out this far before now.

[kaushikcfd]

This is not a new "feature" of distributed lazy. It has always been this way - although now it may be a little worse than before.

It's safe to call it a bug that ought to be fixed, since, purely based on principles we don't expect it. If you have access to the log files for the 3 mesh sizes that would be helpful to debug this behavior.

[MTCam]

• Do you happen to know the average/max number of neighbors per rank? The only well-justified cost growth in compilation would be linear in the max number of neighbors.

@inducer : I am going to post a running report on the number of rank neighbors for runs. As the data comes in, I'll edit this msg. I hope it can help us suss out whether the number of partition boundaries is playing a significant role.

nranks	nnbrs (min, max, avg)
1	(0, 0, 0)
2	(1 , 1, 1)
4	(3, 3, 3)
8	(7, 7, 7)
16	(10, 11, 10.9)
32	(10, 15, )
64	(11, 19, )
128	(12, 19, ) most with 17, 18

[MTCam]

If you have access to the log files for the 3 mesh sizes that would be helpful to debug this behavior.

1.5M elems, 32 ranks: https://gist.github.com/db1139f6db451f0cd921fab6a6d69223
3M elems, 64 ranks: https://gist.github.com/07237afdfd5a561f250de8f1cb0244a9
6M elems, 128 ranks: https://gist.github.com/b9518b0f26d4a5269d5bf4082c76f9a5

Hope this helps! Thanks a ton for looking into this!

[kaushikcfd]

Thanks, so here is the ranks vs number of pytato nodes:

nranks	max. num pytato nodes across ranks
32	40703
64	76475
128	106285
256	136095

So, at least it's not growing linearly and starts tapering off once our number of neighboring ranks starts to taper off. So @inducer's hunch was correct, but I'm very surprised how come communication nodes outgrow the total numberof nodes on a single rank. Either something in grudge is adding redundant nodes or our materialization strategy in pytato is broken.

Because of these many nodes, we materialize more arrays and the generated loopy kernel has an increasing number of statements.

nranks	max. num lpy insns across ranks and parts	max. time for transform_loopy_program
32	1710	145 secs
64	3204	435 secs
128	4449	807 secs

[inducer]

Either something in grudge is adding redundant nodes or our materialization strategy in pytato is broken.

I completely agree! FWIW, I think that's also consistent with where we left off in terms of the mysteriously-large sum of flux expressions in July. I strongly suspect that something there is causing us to do way more work than necessary. Indirectly, I think this is also what caused the unreasonably-large argument counts. I tried to figure out for a while what that might be, but I ran out of time before my trip. To allow us to run at all, we decided to pursue SVM, but I think the underlying issue remains.