Skip to content

Introduce MPI parallelization - #197

Draft
BenjaminRodenberg wants to merge 13 commits into
mainfrom
mpi-parallelization
Draft

BenjaminRodenberg wants to merge 13 commits into
mainfrom
mpi-parallelization

Conversation

@BenjaminRodenberg

@BenjaminRodenberg BenjaminRodenberg commented Oct 20, 2026 •

Copy link
Copy Markdown
Collaborator

This PR adds MPI parallelization + domain decomposition.

Marking this as draft because the scope of this PR is currently huge and there are a lot of restrictions and loose ends. I will try to reduce the size of this PR by merging smaller PRs first.

Author's Todos

@BenjaminRodenberg BenjaminRodenberg self-assigned this Oct 20, 2026
@BenjaminRodenberg BenjaminRodenberg added the enhancement New feature or request label Oct 20, 2026
BenjaminRodenberg added a commit that referenced this pull request Oct 20, 2026
Assisted-by: Claude Code:claude-opus-5-5
@BenjaminRodenberg BenjaminRodenberg added this to the v0.9.0 milestone Oct 20, 2026
BenjaminRodenberg added a commit that referenced this pull request Oct 20, 2026
Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
Every column gets an id that is unique over the whole domain,
grid["column_id"]: the Elmer element id, or the flat index into the
2-D fields of a MATLAB or NetCDF grid. ColumnGather collects
per-column arrays of all ranks on rank 0 in id order, and the output
file and --dump-reference are written from it on rank 0, as a single
rank would write them. Only structured grids are gathered for the
output file; a partitioned mesh still skips it.

timelastsnow_netCDF is now always float64, since date2num returns
integers when all times are whole days, which ranks could disagree
on.

Taken from 8b5a9eb, def15dc, c62c9ff and 72e8d91 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The shading look-up table is built from the terrain around each
column and the grid spacing of the 2-D DEM, which a rank holding only
part of the columns does not have. Shading on several ranks is now a
usage error asking for --no-shading, and _precompute_shading_matlab
asserts that it is not called in such a run.

No grid that supports shading runs on several ranks yet; the next
commits split MATLAB grids over the ranks.

Taken from 43a676e of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
On several ranks, each rank now computes a contiguous block of the
columns of a MATLAB grid, the first block on rank 0. The 2-D fields
stay whole on every rank: the slope needs the neighbours of each cell,
and the output writer places the columns of all ranks into them.

So that the results equal those of a single rank, the mean latitude
that selects the water vapour coefficient of SWin and the lapse rate
of potential temperature are computed for the whole domain, the latter
fitted on rank 0 in column order. This also applies to a partitioned
mesh, whose ranks previously used the values of their own partition.
Without --random-seed, rank 0 draws the seed and shares it with all
ranks.

Taken from 46f47a2, dad8dc4, f9ad61a and 071f0c3 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The MATLAB example runs for one day on a single rank and on 2, 3 and
4 ranks under mpirun; each parallel run must give a byte-identical
output file and an equal --dump-reference snapshot. CI also runs the
MATLAB example on 4 ranks.

Taken from db2b582 and e26d3e8 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
Every rank printed every message, so a run on N ranks repeated each
line N times. Ranks other than 0 now only log warnings and errors to
the console; log files per rank keep the full level. print_time logs
instead of printing, and --print-diagnostics gathers the columns of
all ranks on rank 0 and reports the whole domain.

Taken from ade4fba of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
A hanging mpirun ran until GitHub's default limit of 6 hours. The MPI
steps now stop after 5 minutes and the whole job after 20; the job
normally takes a few minutes, mostly for installing and the serial
Greenland examples.

Taken from 47b4c8f of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
setup-mpi makes Open MPI use TCP between ranks. On GitHub's runners it
then also uses the docker0 interface, over which messages never arrive,
so large messages hung on 4 ranks, e.g. the gathers of column ids in
the parallel tests and the build-and-test MPI example.

Taken from 5e28f4d of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
On several ranks, each rank now computes a contiguous block of the
columns of a MATLAB grid, the first block on rank 0. The 2-D fields
stay whole on every rank: the slope needs the neighbours of each cell,
and the output writer places the columns of all ranks into them.

So that the results equal those of a single rank, the mean latitude
that selects the water vapour coefficient of SWin and the lapse rate
of potential temperature are computed for the whole domain, the latter
fitted on rank 0 in column order. This also applies to a partitioned
mesh, whose ranks previously used the values of their own partition.
Without --random-seed, rank 0 draws the seed and shares it with all
ranks.

Taken from 46f47a2, dad8dc4, f9ad61a and 071f0c3 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The MATLAB example runs for one day on a single rank and on 2, 3 and
4 ranks under mpirun; each parallel run must give a byte-identical
output file and an equal --dump-reference snapshot. CI also runs the
MATLAB example on 4 ranks.

Taken from db2b582 and e26d3e8 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
Every rank printed every message, so a run on N ranks repeated each
line N times. Ranks other than 0 now only log warnings and errors to
the console; log files per rank keep the full level. print_time logs
instead of printing, and --print-diagnostics gathers the columns of
all ranks on rank 0 and reports the whole domain.

Taken from ade4fba of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
A hanging mpirun ran until GitHub's default limit of 6 hours. The MPI
steps now stop after 5 minutes and the whole job after 20; the job
normally takes a few minutes, mostly for installing and the serial
Greenland examples.

Taken from 47b4c8f of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
setup-mpi makes Open MPI use TCP between ranks. On GitHub's runners it
then also uses the docker0 interface, over which messages never arrive,
so large messages hung on 4 ranks, e.g. the gathers of column ids in
the parallel tests and the build-and-test MPI example.

Taken from 5e28f4d of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
On a partitioned Elmer mesh, each rank used the mean latitude of its
own partition to select the water vapour coefficient of SWin, fitted
the lapse rate of potential temperature over its own columns, and,
without --random-seed, drew its own random seed. The results thus
depended on the partitioning.

The mean latitude is now reduced over all ranks in init_grid
(grid["lat_mean_whole_domain"]), the lapse rate is fitted on rank 0
in column id order and broadcast, and rank 0 draws the seed and
shares it. On a single rank each of these is a no-op.

A new test splits a small generated mesh into 3 partitions and
compares --dump-reference against a serial run on the whole mesh.

Taken from 3cc2b1e of mpi-parallelization-mvp-matlab, from #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
Every column gets an id that is unique over the whole domain,
grid["column_id"]: the Elmer element id, or the flat index into the
2-D fields of a MATLAB or NetCDF grid. ColumnGather collects
per-column arrays of all ranks on rank 0 in id order.
--dump-reference and --print-diagnostics now use it, so on a
partitioned mesh rank 0 writes and reports the whole domain instead
of every rank its own partition.

timelastsnow_netCDF is now always float64: date2num returns integers
when all times are whole days, ranks can disagree on that, and
ColumnGather requires equal dtypes.

Taken from 261e83f and 88e576c of mpi-parallelization-mvp-matlab,
from #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
On a partitioned Elmer mesh, each rank used the mean latitude of its
own partition to select the water vapour coefficient of SWin, fitted
the lapse rate of potential temperature over its own columns, and,
without --random-seed, drew its own random seed. The results thus
depended on the partitioning.

The mean latitude is now reduced over all ranks in init_grid
(grid["lat_mean_whole_domain"]), the lapse rate is fitted on rank 0
in column id order and broadcast, and rank 0 draws the seed and
shares it. On a single rank each of these is a no-op.

A new test splits a small generated mesh into 3 partitions and
compares --dump-reference against a serial run on the whole mesh.

Taken from 3cc2b1e of mpi-parallelization-mvp-matlab, from #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
setup-mpi makes Open MPI use TCP between ranks, and on GitHub's
runners Open MPI's TCP transport also picks the docker0 interface,
over which messages never arrive. Small messages still get through,
but large ones hang: in #197 a gather on 4 ranks hung until GitHub's
limit of 6 hours, while the same run passed locally, where the ranks
use shared memory. The next commits run the MATLAB example on up to
4 ranks, with gathers large enough to hang.

Open MPI now only uses the loopback interface in both workflows. It
is set in mca-params.conf, since tox does not pass OMPI_MCA_*
variables to the tests. In addition, the build-and-test job now stops
after 20 minutes and its MPI steps after 5, so a hang fails quickly
instead of blocking a runner for 6 hours.

Taken from 5e28f4d and 47b4c8f of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The output writer now gathers the columns of all ranks with
ColumnGather and rank 0 writes them to the file a single rank would
write, one variable at a time to limit its memory. Only structured
grids are gathered; a partitioned mesh still skips the output file.

Taken from 8b5a9eb of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The shading look-up table is built from the terrain around each
column and the grid spacing of the 2-D DEM, which a rank holding only
part of the columns does not have. Shading on several ranks is now a
usage error asking for --no-shading, and _precompute_shading_matlab
asserts that it is not called in such a run.

No grid that supports shading runs on several ranks yet; the next
commits split MATLAB grids over the ranks.

Taken from 43a676e of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The shading look-up table is computed by walking from every glacier
cell of the 2-D grid, not only from the columns of one rank, so a
rank holding part of the columns cannot build its rows of the table.
Shading on several ranks is now a usage error asking for --no-shading,
and _precompute_shading_matlab asserts that it is not called in such
a run.

No grid that supports shading runs on several ranks yet; the next
commits split MATLAB grids over the ranks.

Taken from 43a676e of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
On several ranks, each rank now computes a contiguous block of the
columns of a MATLAB grid, the first block on rank 0. The 2-D fields
stay whole on every rank: the slope needs the neighbours of each cell,
and the output writer places the columns of all ranks into them.

The domain-wide latitude mean, lapse rate and random seed of #203 make
the results equal those of a single rank.

Taken from 46f47a2 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The MATLAB example runs for one day on a single rank and on 2, 3 and
4 ranks under mpirun; each parallel run must give a byte-identical
output file and an equal --dump-reference snapshot. CI also runs the
MATLAB example on 4 ranks.

Taken from db2b582 and e26d3e8 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
setup-mpi makes Open MPI use TCP between ranks, and on GitHub's
runners Open MPI's TCP transport also picks the docker0 interface,
over which messages never arrive. Small messages still get through,
but large ones hang: in #197 a gather on 4 ranks hung until GitHub's
limit of 6 hours, while the same run passed locally, where the ranks
use shared memory. The next commits run the MATLAB example on up to
4 ranks, with gathers large enough to hang.

Open MPI now only uses the loopback interface in both workflows. It
is set in mca-params.conf, since tox does not pass OMPI_MCA_*
variables to the tests. In addition, the build-and-test job now stops
after 20 minutes and its MPI steps after 5, so a hang fails quickly
instead of blocking a runner for 6 hours.

Taken from 5e28f4d and 47b4c8f of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The output writer now gathers the columns of all ranks with
ColumnGather and rank 0 writes them to the file a single rank would
write, one variable at a time to limit its memory. Only structured
grids are gathered; a partitioned mesh still skips the output file.

Taken from 8b5a9eb of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The shading look-up table is computed by walking from every glacier
cell of the 2-D grid, not only from the columns of one rank, so a
rank holding part of the columns cannot build its rows of the table.
Shading on several ranks is now a usage error asking for --no-shading,
and _precompute_shading_matlab asserts that it is not called in such
a run.

No grid that supports shading runs on several ranks yet; the next
commits split MATLAB grids over the ranks.

Taken from 43a676e of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
On several ranks, each rank now computes a contiguous block of the
columns of a MATLAB grid, the first block on rank 0. The 2-D fields
stay whole on every rank: the slope needs the neighbours of each cell,
and the output writer places the columns of all ranks into them.

The domain-wide latitude mean, lapse rate and random seed of #203 make
the results equal those of a single rank.

Taken from 46f47a2 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
The MATLAB example runs for one day on a single rank and on 2, 3 and
4 ranks under mpirun; each parallel run must give a byte-identical
output file and an equal --dump-reference snapshot. CI also runs the
MATLAB example on 4 ranks.

Taken from db2b582 and e26d3e8 of #197.

Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 21, 2026
Assisted-by: Claude Code:claude-opus-5-5
BenjaminRodenberg added a commit that referenced this pull request Oct 22, 2026
Assisted-by: Claude Code:claude-opus-5-5
On several ranks, a column held by more than one rank would be computed
and written twice. Before the time loop, rank 0 gathers the column ids
of all ranks, and every rank stops with an error if an id occurs
twice. The check sends the ids of all columns to rank 0, so
--no-check-grid-distribution skips it; rank 0 says so when it passes.

Assisted-by: Claude Code:claude-opus-5-5
A partitioned mesh skipped writing output, and the output of an
unstructured grid could not be gathered over the ranks. Rank 0 now
also gathers the cells of an unstructured grid, ordered by Elmer
element id, together with the meshes of all ranks; vertices that
several ranks share are written once.

Assisted-by: Claude Code:claude-opus-5-5
On several ranks, each rank now keeps a contiguous block of the cells
of the mesh, the first block on rank 0, and only the vertices those
cells use. The projection of the vertices and the cell centres moves
from the Mesh constructor to Mesh.compute_geometry(), which INIT calls
after the split, so each rank transforms only its own cells. The test
elevations of the ELMER grid type are drawn for the whole mesh before
the split, so every column keeps its serial value.

The startup check accepts ELMER and CUSTOM grids on a non-partitioned
mesh now. Rank 0 warns that fitting the lapse rate over all columns
synchronizes the ranks every time step.

Assisted-by: Claude Code:claude-opus-5-5
Each rank computes exactly one partition, so a run on fewer ranks than
the mesh has partitions skipped the others without notice; the grid
distribution check only finds columns held twice, not missing ones.
GridConfig.num_partitions now also counts the partitions of a
partitioned mesh, by its part.N.header files, and a mismatch with the
number of ranks is a usage error. --allow-partition-mismatch turns it
into a warning, e.g. to compute a single partition with --use-part.

Assisted-by: Claude Code:claude-opus-5-5
A NetCDF grid is now split into contiguous blocks of its columns like a
MATLAB grid, and the startup check accepts it on several ranks.
GridConfig now splits every non-partitioned grid.

Assisted-by: Claude Code:claude-opus-5-5
A partitioned mesh skipped writing the restart file with only a
warning at the end of the run, also on a single rank with --use-part.
--restart-dir and --restart-init are now a usage error for partitioned
meshes, as on several ranks, and the restart writer runs in every
run: it writes nothing without --restart-dir.

That makes --dump-reference contain timelastsnow_netCDF on partitioned
meshes too, so the partitioned-mesh test compares it as well.

Assisted-by: Claude Code:claude-opus-5-5
The serial-vs-parallel comparison gets Elmer cases: a small synthetic
Elmer mesh, committed under tests/data, on 3 ranks, and the same mesh
with fake coupling to Elmer/Ice, ICON and ICON-Land on 2 ranks.
tests/data also holds an XIOS DEM for the mesh, kept out of Git LFS.

CI also runs the pseudo-coupled Greenland example on 4 ranks.

Assisted-by: Claude Code:claude-opus-5-5
Assisted-by: Claude Code:claude-opus-5-5
Assisted-by: Claude Code:claude-opus-5-5
read_dem_xios looked up zs and h by the mesh's local vertex indices.
After the split over the ranks each rank renumbers its vertices, so
--netcdf-dem-mesh-unstructured was rejected on several ranks. It now
looks the values up by Elmer node id, node id i at position i - 1,
which holds for every part of the mesh and for any order of
mesh.nodes. A file without values for all node ids of the mesh is
rejected.

Elmer/Ice output for greenland_mesh_v0 stores its nodes in this order;
in that mesh it is also the order of mesh.nodes, so serial results are
unchanged.

With that, the startup check of the grid type on several ranks has
nothing left to reject and is removed.

Assisted-by: Claude Code:claude-opus-5-5
A non-partitioned Elmer mesh is split over the ranks only after every
rank has read all of it. Rank 0 now warns about this for the ELMER,
CUSTOM and ELMERXIOS grid types and points to a mesh with one
partition per rank, which requires --netcdf-dem-mesh.

Assisted-by: Claude Code:claude-opus-5-5
Restart files were rejected on several ranks and with partitioned
meshes. Rank 0 now gathers the columns of all ranks one variable at a
time and writes them ordered by column id, together with the ids in a
new variable column_id. Reading, rank 0 matches the rows of the file
to the columns by id and scatters them with the new
ColumnGather.scatter, so a file continues on any number of ranks and
on a partitioned mesh.

A file without column_id holds the columns in the order of a serial
run, which for a split grid is the ranks' columns in rank order; on
several ranks its x and y are checked against the grid.

Assisted-by: Claude Code:claude-opus-5-5

This branch was successfully deployed

1 active (outdated) deployment
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant