Distributed output: zarr and parallel netCDF backends, gather/rank-block write modes - #1403
Distributed output: zarr and parallel netCDF backends, gather/rank-block write modes#1403kotsaloscv wants to merge 21 commits into
Conversation
…k export Adds distributed (MPI) output to the IO layer and the standalone driver: - Output distribution strategies (io/distributed.py): SingleNode, Gather (owned entries collected on the root rank, global order) and RankBlock (every rank writes its own block of a shared store). - ZarrWriter (zarr format 3) alongside the NETCDFWriter, sharing coordinate and CF attribute definitions so both formats stay identical. - Per-field-group `backend` (netcdf|zarr) and `mode` (gather|distributed) config, validated at construction; distributed netCDF is rejected until a parallel netCDF writer exists. - Driver wiring: --output-backend / --output-mode CLI options, IO monitor built with process_props + decomposition_info, output timers and a wall-time share column in the timer report. - icon4py-gridlook CLI: export driver output (zarr or netCDF) to a gridlook-readable store and serve it over HTTP. - Tests: unit tests for distributions, zarr writer and gridlook, MPI tests on a synthetic decomposition and a parallel-vs-single-rank driver check.
msimberg
left a comment
There was a problem hiding this comment.
Looking pretty nice, thanks @kotsaloscv! Some general comments and questions, not everything is blocking.
| backend="zarr", # type: ignore[arg-type] # value strings are coerced on purpose | ||
| mode="distributed", # type: ignore[arg-type] # value strings are coerced on purpose |
There was a problem hiding this comment.
I see the type ignore, but I don't understand what the "on purpose" is meant for. The config should be taking enums, not strings. The yaml config to config conversion should take care of converting strings to enums.
There was a problem hiding this comment.
You are right. Dropping the "on purpose", Claude being peculiar.
There was a problem hiding this comment.
Hmm, dropping the "on purpose" text was not what I meant, sorry for ambiguity. I meant what is the "on purpose" part even referring to. As mentioned in other comments, strings should not be accepted here in my opinion. E.g. this test I think should not be a thing, or should not be allowed:
def test_fieldgroup_config_accepts_backend_and_mode_value_strings() -> None:
config = FieldGroupIOConfig(
filename="a.nc",
variables=["air_density"],
backend="netcdf", # type: ignore[arg-type]
mode="gather", # type: ignore[arg-type]
)
assert config.backend is OutputBackend.NETCDF
assert config.mode is OutputMode.GATHER
There was a problem hiding this comment.
Agreed, and that test is gone. Now that string acceptance is removed from the dataclass (see the __post_init__ thread), the test is inverted: test_fieldgroup_config_rejects_backend_and_mode_strings pins that passing a value string raises InvalidConfigError. String-to-enum conversion is covered where it now lives — the config-file boundary — by the ExperimentConfig yaml roundtrip tests from #1391.
| nc_comment: str = "ICON inspired code in Python and GT4Py" | ||
|
|
||
| def __post_init__(self) -> None: | ||
| # normalize once: value strings ("zarr") are accepted and coerced to the enums |
There was a problem hiding this comment.
Commented elsewhere, but I think we should not be accepting plain strings here. The only place where that should be allowed is in the conversion from yaml dicts to a proper config object.
This is probably a good place to sync with @DropD to decide in which order changes should be made: whether it makes sense to change these later or wait for @DropD's changes so that it can be incorporated here.
There was a problem hiding this comment.
Agree. Let me keep it as is for now, and I will sync with Rico for this.
There was a problem hiding this comment.
I think I still see __post_init__ coercing backend strings to enums... Am I looking at the wrong changes are are there still places where the strings are expected? I'd remove all uses of strings for config options except when reading from a config file.
There was a problem hiding this comment.
You were looking at the right place — the coercion was still there, pending the config-file boundary. With #1391 merged that boundary now exists, so this is done: OutputBackend/OutputMode are registered with common.config.config_io (@config_io.register_enum), and value strings are converted exactly once, when reading a config file. __post_init__ no longer coerces anything — it rejects a non-enum value with an InvalidConfigError pointing to that boundary. The remaining entry points already deliver enums: the driver CLI (typer converts the option strings) and create_io_monitor (enum fields of DriverConfig).
havogt
left a comment
There was a problem hiding this comment.
Just some high-level comments: Mainly, I'd split the unrelated features from this PR (especially gridlook which is completely orthogonal, but huge on its own).
| @@ -0,0 +1,563 @@ | |||
| # ICON4Py - ICON inspired code in Python and GT4Py | |||
There was a problem hiding this comment.
I'd move this out of this PR, because it's an unrelated feature.
There was a problem hiding this comment.
While I fully agree with your comment, I would insist on keeping it here.
My approach with gridlook is that since we are actively developing the driver, it is very essential and useful to have an idea of the generated output. Gridlook does it very efficiently and minimally.
Also, all gridlook additions are extremely self-contained. If we see that it is not as useful, we drop it anytime simply by deleting the corresponding file.
There was a problem hiding this comment.
Not sure I understand why it needs to come in the same PR? Makes history cleaner especially since it's touching dependencies for the different features.
| if output_backend == common_io.OutputBackend.NETCDF: | ||
| return xr.open_dataset(path, decode_times=False) | ||
| return xr.open_zarr(path, decode_times=False) |
There was a problem hiding this comment.
Regarding your comment on the total_wall_time, it is there because I wanted to show the percentage of the output compared to the whole driver execution.
- default output backend/mode: zarr + distributed (config, driver, CLI) - split writers.py into netcdf_writers.py and zarr_writers.py, keeping the shared surface (FieldWriter protocol, dimension names, attributes) in writers.py; writers are context managers now - distributed.py: factor host (numpy) conversion of owner masks/global indices into _host_owner_data with the cupy/numpy rationale documented - make total_wall_time a required argument of show_timer_report - tests: use test_utils.assert_dallclose, match statements for backend dispatch, real types instead of object in the recording monitor, netcdf/gather pinned explicitly where tests assert netCDF artifacts - ci: revert default.yml MODEL_MPI_SUBSETS to datatest
- NETCDFWriter writes distributed (rank-block) output through a shared file opened with parallel=True, mirroring the rank-block zarr layout exactly (padded per-rank blocks, one chunk per block, global_index coordinates); every rank writes its full block so empty ranks still participate in every collective write, and vertical chunks shrink to respect HDF5's 4 GiB chunk limit. - distributed+netcdf is now accepted iff the netCDF4 installation is an MPI-parallel build (PyPI wheels are serial): checked at config validation, writer construction and file open, with the build steps spelled out in the error messages and the io module docs. - The gridlook exporter handles rank-block netCDF sources and keeps the exported store's metadata strict JSON (no bare NaN tokens). - Tests: netCDF rank-block layout through a serial handle, support predicate branches, monitor-to-writer wiring, gridlook rank-block netCDF export, and netcdf+distributed MPI cases that skip with an explicit reason on serial installations.
Chunking and sharding become per-field-group config: - FieldGroupIOConfig.horizontal_chunk_size sets the chunk size of the horizontal axes for both backends (default unchanged: one chunk per rank block in distributed mode); horizontal_shard_size (zarr only) groups whole chunks into one storage file each -- the file-count knob on parallel file systems. - RankBlockDistribution rounds the uniform block size up to the group's chunk/shard granularity (block_alignment), so chunks and shards never cross rank-block boundaries and concurrent rank writes stay in disjoint files; RankBlock.chunk is renamed to size accordingly. IOMonitor shares distributions per (mode, alignment). Review-pass fixes: - both writers validate fields (canonical dims, horizontal dimension, CF attributes) on every rank before any file mutation, so a bad field can no longer leave a phantom time slice or hang non-root ranks in a barrier - the netCDF writer resolves existing variables by name (like the zarr writer) instead of by standard_name, and close() is a no-op when the file was never opened - ZarrWriter rejects shard sizes without a dividing chunk size at construction, on every rank - a UGRID export failure is broadcast so all ranks raise together - time-delta output intervals must be a multiple of the model time step (previously rounded silently) - IOConfig.time_units/calendar are real config fields now, forwarded to the writers Comments and docstrings updated to the implemented behavior and deduplicated (each invariant explained once, at its home).
A rank-block axis never holds fewer chunks or shard files than writing ranks (the default one-chunk-per-block layout is that minimum); dedicated IO ranks (async output) will lower the floor to the IO-rank count.
nfarabullini
left a comment
There was a problem hiding this comment.
very trivial comments on my side. I see that you already have others. Perhaps once everything is tackled, I can do a second round of review
| grid_file_name: pathlib.Path, | ||
| grid_id: uuid.UUID, | ||
| dtime: datetime.timedelta, | ||
| process_props: decomposition_defs.ProcessProperties | None, |
There was a problem hiding this comment.
for mpi-related tests we usually have a separate file, e.g. for the metrics states: model/common/tests/common/metrics/mpi_tests
There was a problem hiding this comment.
These tests don't use MPI — the recorder stub just mirrors IOMonitor's full constructor signature (including the process_props/decomposition_info parameters) to pin the wiring; everything runs single-process on the simple grid. The driver's real MPI tests do live in a separate file, mpi_tests/test_parallel_driver_io.py ;)
| - `nc_title` (optional): Title field of the generated netcdf file. | ||
| - `nc_comment` (optional): Comment to be put to generated netcdf file. | ||
| - `backend` (default="zarr"): File format of the group, `"netcdf"` or `"zarr"`. | ||
| - `mode` (default="distributed"): Write strategy of distributed (MPI) runs: `"gather"` collects all fields on the root rank which writes them in global order; `"distributed"` lets every rank write its owned entries into a rank-contiguous block of a shared store (see `io.distributed`). Single-rank runs write the full state either way, but `"distributed"` with the `"netcdf"` backend requires an MPI-parallel netCDF4 installation and is rejected at configuration time on serial installations regardless of the rank count (see "Parallel netCDF" below). |
There was a problem hiding this comment.
maybe place all modes in indented bullet points?
There was a problem hiding this comment.
Done — the mode entry now lists "gather" and "distributed" as indented bullets (and picked up the new rank-block layout-marker documentation while at it).
| ... | ||
|
|
||
| @property | ||
| def file_horizontal_size(self) -> base.HorizontalGridSize: |
There was a problem hiding this comment.
| def file_horizontal_size(self) -> base.HorizontalGridSize: | |
| def output_horizontal_size(self) -> base.HorizontalGridSize: |
maybe? file is a bit too vague IMO
| every write). ``global_index`` maps the block's entries to their positions in the | ||
| undecomposed global grid of ``global_size`` entries (padding entries carry no | ||
| global index). | ||
| """ |
There was a problem hiding this comment.
can you make the docstring such that it has bullet points which describe what each param indicates?
There was a problem hiding this comment.
Done — the docstring now has an Attributes: section with one bullet per field (start, count, size, padded_size, global_size, global_index).
| #: rounded up to a multiple of this value (see | ||
| #: ``distributed.check_chunks_align_with_blocks``). | ||
| horizontal_chunk_size: int | None = None | ||
| #: Zarr only: entries per shard along the horizontal axes, grouping whole chunks |
There was a problem hiding this comment.
can you reduce the amount of comments here? It seems a bit excessive
There was a problem hiding this comment.
Trimmed — the cross-references and tuning commentary are gone; the details stay in the common.io module docs.
| "(netCDF4.__has_parallel4_support__ is false), as all PyPI wheels are" | ||
| ) | ||
| if importlib.util.find_spec("mpi4py") is None: | ||
| return "the 'mpi4py' package is not installed" |
There was a problem hiding this comment.
should this be an error or at least a warning?
There was a problem hiding this comment.
It effectively is an error: missing_parallel_support is a predicate that only returns the reason, and every caller raises with it — the writer at construction when a multi-rank communicator shares the file, and again at file open. Returning instead of raising keeps it usable for skip conditions too (the MPI test skips with this exact reason).
| self.dataset.createDimension(writers.MODEL_HALF_LEVEL, self.num_interfaces) | ||
| self.dataset.createDimension(writers.CELL, self._horizontal_size.num_cells) | ||
| self.dataset.createDimension(writers.VERTEX, self._horizontal_size.num_vertices) | ||
| self.dataset.createDimension(writers.EDGE, self._horizontal_size.num_edges) |
There was a problem hiding this comment.
what if we have offset dimensions, e.g. E2CDim?
There was a problem hiding this comment.
Fields carrying sparse/offset dimensions are not writable: writers.canonicalize_time_slice accepts only fields with one horizontal (cell/edge/vertex) plus one vertical dimension and raises for anything else, on every rank, before the file is touched — so an offset dimension never reaches createDimension. Output of connectivity-shaped fields would be a new feature (none of the current output variables need it).
| self.dataset.createDimension(writers.EDGE, self._horizontal_size.num_edges) | ||
| log.debug(f"Creating dimensions {self.dataset.dimensions} in {self._file_name}") | ||
| # create time variables | ||
| times = self.dataset.createVariable(writers.TIME, "f8", (writers.TIME,)) |
There was a problem hiding this comment.
why do we not create variables for cell, edge, and vertex as well?
There was a problem hiding this comment.
The horizontal coordinates (lon/lat per cell/edge/vertex) live in the UGRID sidecar file (*_ugrid.nc) the monitor writes once at startup; the data files reference it through the UGRID association instead of duplicating the coordinates into every rollover file. Inlining clon/clat into the data files would be a possible follow-up if a consumer needs fully self-contained files.
| return state | ||
|
|
||
|
|
||
| def create_monitor( |
There was a problem hiding this comment.
maybe move these functions into a common.py file in the test folder
There was a problem hiding this comment.
They're only used by this one file so far, so I'd keep them local for now and move them to a shared module the moment a second file needs them
| if process_props.rank != 0: | ||
| return |
There was a problem hiding this comment.
so is this test not supposed to run for distributed?
There was a problem hiding this comment.
It runs for all modes (gather and distributed are both in the parametrization) and on every rank — the writing part above is fully collective. The early return only limits the verification (reading the finished files back and comparing values) to rank 0, after the barrier; the other ranks have nothing left to do at that point.
jcanton
left a comment
There was a problem hiding this comment.
overall very nice.
I would prefer to separate gridlook into its own PR though.
I'm also not sure it belongs inside the driver: it's a postprocessing visualization tool, not a model running tool
|
|
||
| for dim, dim_name in HORIZONTAL_DIM_NAMES.items(): | ||
| mask, owned_global_index = _host_owner_data(decomposition_info, dim) | ||
| row_counts = self._allgather_counts(owned_global_index.shape[0]) |
There was a problem hiding this comment.
I don't like rows: there are no rows in our unstructured triangular grids,
this is more nr_of_local_cell/edge/vertex right?
horizontal_counts / entry_count / dim_count / local_size (since you're using global_size just below) ?
in any case, these shuold also all be available from DecompositionInfo and not need to be re-computed here
or am I totally misunderstanding?
There was a problem hiding this comment.
again, look at parallel_helpers.py
There was a problem hiding this comment.
Agreed on the naming — there are no rows here. Renamed throughout: entry_counts, _gather_entries, entries. On DecompositionInfo: the local owned count is indeed available from it, but these are the per-rank count vectors (needed for the Gatherv receive counts/displacements and the block sizing), and DecompositionInfo is rank-local and holds no communicator — so an exchange is unavoidable. What was re-computed in two places is now consolidated: both distributions share module-level helpers (_allgather_entry_counts, _gather_entries, _check_partition), and RankBlockDistribution now runs the same partition validation as GatherDistribution through them.
…ation, robust timings - Rank-block variables no longer assert a UGRID mesh association their rank-ordered, padded axes cannot honor: 'mesh'/'location' are replaced by the marker attribute 'icon4py_layout = "rank_block"', documented on RankBlock and in the module docs. - RankBlockDistribution now validates (like GatherDistribution, via a shared helper) that the owner masks partition the global grid, rejects non-positive block alignments, hands out read-only global indices and logs the padding ratio. - The counts/gather/partition-check machinery of both distributions is consolidated into shared module functions; 'rows' terminology replaced by 'entries' (there are no rows in an unstructured grid). - report_timings is communication-free now: a reduction over the ranks could leave surviving ranks blocked when one rank failed earlier. Every rank logs its own totals instead. - The driver assembles diagnostics and samples the output timers only at capture steps; other steps merely advance the schedule counters. - The distributed-netCDF parallel-support check is rank-aware and lives solely in the writer: single-rank runs on serial installations are valid now, and the config-time duplicate is gone. - Review mechanics: output_horizontal_size and _is_distributed renames, generate_name takes an explicit suffix and sits next to its caller, filter_by_standard_name removed (only tests used it), docstring/comment cleanups, .gitignore change reverted.
|
cc: @jcanton @havogt @msimberg |
…mpi_decomposition The output distributions now hold a host-converted DecompositionInfo (new DecompositionInfo.as_host, one device-to-host conversion at construction) and query owner masks and owned global indices from it instead of keeping private copies. The communicator-level helpers (allgather_entry_counts, gather_entries, check_owned_indices_partition) move from io/distributed.py to decomposition/mpi_decomposition.py, where they sit with the rest of the communicator machinery and are reusable beyond IO. io/distributed.py keeps only the IO-specific pieces.
|
remember to check if this needs to be updated with new instructions |
| def allgather_entry_counts(process_props: decomp_defs.ProcessProperties, count: int) -> np.ndarray: | ||
| """Entry counts of all ranks (one entry per rank, in rank order). | ||
|
|
||
| Collective on the communicator; works on a single-rank communicator without an | ||
| MPI installation. | ||
| """ | ||
| if process_props.is_single_rank(): | ||
| return np.asarray([count], dtype=np.int64) | ||
| return np.asarray(process_props.comm.allgather(count), dtype=np.int64) | ||
|
|
||
|
|
||
| def gather_entries( | ||
| process_props: decomp_defs.ProcessProperties, | ||
| local_entries: np.ndarray, | ||
| *, | ||
| entry_counts: np.ndarray, | ||
| ) -> np.ndarray | None: | ||
| """Concatenate the ranks' entries (leading-axis) on the root rank, in rank order. | ||
|
|
||
| ``entry_counts`` holds the per-rank leading-axis sizes (see | ||
| ``allgather_entry_counts``); trailing axes must agree between the ranks. Returns | ||
| the concatenation on the root rank and None on all other ranks. Collective on the | ||
| communicator; works on a single-rank communicator without an MPI installation. | ||
| """ | ||
| if process_props.is_single_rank(): | ||
| return local_entries | ||
| send = np.ascontiguousarray(local_entries) | ||
| entry_elements = int(np.prod(send.shape[1:], dtype=np.int64)) | ||
| if process_props.rank == 0: | ||
| gathered = np.empty((int(entry_counts.sum()), *send.shape[1:]), dtype=send.dtype) | ||
| process_props.comm.Gatherv(send, [gathered, entry_counts * entry_elements], root=0) | ||
| return gathered | ||
| process_props.comm.Gatherv(send, None, root=0) | ||
| return None | ||
|
|
||
|
|
||
| def check_owned_indices_partition( | ||
| process_props: decomp_defs.ProcessProperties, | ||
| dim_name: str, | ||
| owned_global_index: np.ndarray, | ||
| entry_counts: np.ndarray, | ||
| ) -> np.ndarray | None: | ||
| """Check that the ranks' owned global indices partition the global grid. | ||
|
|
||
| Collective: the owned global indices of all ranks are gathered on the root rank | ||
| and verified to be a permutation of ``0..N-1`` -- overlapping or gappy owner | ||
| masks would otherwise reassemble a plausible-looking but wrong global field. The | ||
| verdict is broadcast so all ranks raise together instead of hanging in the next | ||
| collective. | ||
|
|
||
| Returns the gathered indices on the root rank (None on all other ranks), for | ||
| reuse as insertion indices. | ||
|
|
||
| Raises: | ||
| ValueError: if the owned global indices are not a permutation of ``0..N-1``. | ||
| """ | ||
| global_size = int(entry_counts.sum()) | ||
| gathered_index = gather_entries(process_props, owned_global_index, entry_counts=entry_counts) | ||
| is_partition = gathered_index is None or np.array_equal( | ||
| np.sort(gathered_index), np.arange(global_size, dtype=np.int64) | ||
| ) | ||
| if not process_props.is_single_rank(): | ||
| is_partition = process_props.comm.bcast(is_partition, root=0) | ||
| if not is_partition: | ||
| raise ValueError( | ||
| f"Owner masks of dimension '{dim_name}' do not partition the global grid: " | ||
| f"the owned global indices of all ranks are not a permutation of " | ||
| f"0..{global_size - 1}." | ||
| ) | ||
| return gathered_index |
There was a problem hiding this comment.
I think these can be generic over the ndarray type, no? I.e. the functions can receive a data_alloc.NDArray and return an array that is consistent with the input arrays that is determined with array_namespace (see other functions for examples). For allgather_entry_counts there's no input array to determine the array_ns, so array_ns should be explicitly passed to the function instead. This is minor, but means we don't have to worry about these not being generic across cpu/gpu arrays in the future.
There was a problem hiding this comment.
Done, exactly as you describe. Indeed, way more generic ;)
| def check_owned_indices_partition( | ||
| process_props: decomp_defs.ProcessProperties, | ||
| dim_name: str, | ||
| owned_global_index: np.ndarray, | ||
| entry_counts: np.ndarray, | ||
| ) -> np.ndarray | None: | ||
| """Check that the ranks' owned global indices partition the global grid. | ||
|
|
||
| Collective: the owned global indices of all ranks are gathered on the root rank | ||
| and verified to be a permutation of ``0..N-1`` -- overlapping or gappy owner | ||
| masks would otherwise reassemble a plausible-looking but wrong global field. The | ||
| verdict is broadcast so all ranks raise together instead of hanging in the next | ||
| collective. | ||
|
|
||
| Returns the gathered indices on the root rank (None on all other ranks), for | ||
| reuse as insertion indices. | ||
|
|
||
| Raises: | ||
| ValueError: if the owned global indices are not a permutation of ``0..N-1``. | ||
| """ | ||
| global_size = int(entry_counts.sum()) | ||
| gathered_index = gather_entries(process_props, owned_global_index, entry_counts=entry_counts) | ||
| is_partition = gathered_index is None or np.array_equal( | ||
| np.sort(gathered_index), np.arange(global_size, dtype=np.int64) | ||
| ) | ||
| if not process_props.is_single_rank(): | ||
| is_partition = process_props.comm.bcast(is_partition, root=0) | ||
| if not is_partition: | ||
| raise ValueError( | ||
| f"Owner masks of dimension '{dim_name}' do not partition the global grid: " | ||
| f"the owned global indices of all ranks are not a permutation of " | ||
| f"0..{global_size - 1}." | ||
| ) | ||
| return gathered_index |
There was a problem hiding this comment.
This function feels a bit confused. It's both checking and returning values in different use cases. I see one call to it where the return value is ignored and one where the return value is used. The "checking" part is effectively re-checking that our domain decomposition is doing the right thing, which we probably don't need to do every time in production. If there's a gap in check, our testing should rather be expanded (there are already quite a few tests for domain decomposition, are these checks missing from there?). And for the return value ignored or not ignored: if there are two separate use cases, would it instead make sense to separate the function(ality) into two separate functions, one for actually collecting the indices and one for checking, so that they can be used separately.
One level more abstract: we may already have the decomposition on every single node at the beginning of a run, if we've done the decomposition ourselves. Is there no way to reuse that instead?
I realize these might not be straightforward refactorings, but I want to at least ask the questions.
There was a problem hiding this comment.
I did the split as you suggested. Now it is more clear.
| self._append_data(state_to_store, model_time) | ||
| start = timeit.default_timer() | ||
| prepared_state = self._distribution.prepare(state_to_store) | ||
| self._phase_seconds[PHASE_DISTRIBUTE].append(timeit.default_timer() - start) |
There was a problem hiding this comment.
Probably not for now, but this looks like something we'll have to unify sooner or later. We should have a common way of adding timers across modules, not ad-hoc in each...
Are these timers used for anything right now, or more for debugging?
There was a problem hiding this comment.
They are for performance visibility: each rank logs its distribute/write totals at finalize, which is how we compare gather vs. distributed (rank-block) on real runs. Agreed on unifying timers across modules.
Resolves the semantic overlap with the minimal config system (#1391): OutputBackend/OutputMode are registered with config_io so the driver ExperimentConfig yaml roundtrip handles them, and the reference test_config.yml gains the two driver output fields.
…rams, split partition check - FieldGroupIOConfig takes enum members only; value strings are converted at the config-file boundary (common.config.config_io, built on #1391) - process_props/decomposition_info/rank_blocks are explicit (no silent single-node fallback) in the monitors, writers and create_io_monitor - gather helpers are generic over the array namespace; the partition check is split into gather_entries + check_global_index_partition - a filename extension of a different backend is rejected at config time - the vertical coordinate attribute dicts move to states/metadata.py - driver README: distributed output documented (no longer 'disabled in MPI runs'), --output-backend/--output-mode and an mpirun example
|
@jcanton Good catch — it did need updating: the README still claimed output is single-node only and disabled in MPI runs. It now documents distributed output, the |
…nerator With GPU backends the DecompositionInfo arrays live on device, but HaloGenerator.from_gids only accepts host indices, so every multi-rank GPU run failed with a nanobind TypeError during pattern setup. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
msimberg
left a comment
There was a problem hiding this comment.
Looking really good, thanks for the changes @kotsaloscv. I have one minor question/suggestion for the filename handling. Not insisting on it if you feel like it doesn't fit.
… halo generator" This reverts commit 8af717e.
…d from the backend The base name is backend-independent (msimberg): no extension is accepted (a configured '.nc'/'.zarr' would end up inside the stem), generate_name no longer strips anything, and the driver constant is now DEFAULT_OUTPUT_BASENAME.
msimberg
left a comment
There was a problem hiding this comment.
Nice, thanks for the changes @kotsaloscv!
|
cscs-ci run default;SESSIONS=model:model_mpi;MODEL_SUBPACKAGES=common:standalone_driver;MODEL_MPI_SUBPACKAGES=common:standalone_driver;MODEL_MPI_SUBSETS=all;LEVELS=unit:integration |
|
cscs-ci run default;SESSIONS=model:model_mpi;MODEL_SUBPACKAGES=common:standalone_driver;MODEL_MPI_SUBPACKAGES=common:standalone_driver;MODEL_MPI_SUBSETS=all;LEVELS=unit:integration |
|
When developing, you can test your changes on CSCS CI before merge with the You can pass options to override pipeline variables, for example:
Available options are:
For each option, See The Merging Once your PR is approved and ready for merging, add it to the merge queue. The Optional Tests To run benchmarks you can use:
For more detailed information please look at CI in the EXCLAIM universe. |
|
@kotsaloscv #1426 is now merged, so GHEX issues should be gone on main (no need to merge/rebase for the merge queue). |
| ) | ||
| self._group = group | ||
|
|
||
| self._barrier() |
There was a problem hiding this comment.
most failure paths in the IO layer are designed to fail collectively (broadcast verdicts, validation before any file mutation), so a root-only raise cannot strand the other ranks in the next collective. This is the one path that is not: if the root's zarr.open_group or array creation raises (disk full, permissions), the non-root ranks block forever in this barrier. Consider wrapping the root-side initialization in the same broadcast-failure pattern used by _write_ugrid. Nothing blocking.
🤖 Written by an agent on behalf of @jcanton
| dim_name: str, decomposition_info: decomposition.DecompositionInfo, field: xr.DataArray | ||
| ) -> np.ndarray: | ||
| """Drop halo entries from the leading (horizontal) axis of a field.""" | ||
| dim = HORIZONTAL_DIMS_BY_NAME.get(dim_name) |
There was a problem hiding this comment.
prepare expects the horizontal dimension first (field.dims[0]), but the writers canonicalize to horizontal-last (canonicalize_time_slice). A field arriving in canonical (level, cell) order would fail here with a confusing "leading dimension 'level'" error before ever reaching the canonicalizer. It works today because the driver always emits horizontal-first, but the contract is implicit; a comment on OutputDistribution.prepare stating the expected dimension order would help. Nothing blocking.
🤖 Written by an agent on behalf of @jcanton
| pad_value = ( | ||
| np.nan | ||
| if np.issubdtype(data.dtype, np.floating) | ||
| else nc.default_fillvals[f"{data.dtype.kind}{data.dtype.itemsize}"] |
There was a problem hiding this comment.
nc.default_fillvals[f"{data.dtype.kind}{data.dtype.itemsize}"] is unguarded. The table covers the usual u/i/f/c variants, but exotic dtypes (bool, float16) raise a KeyError deep inside append. Nothing blocking; a .get() fallback or a small fill-value table would be more robust.
🤖 Written by an agent on behalf of @jcanton
jcanton
left a comment
There was a problem hiding this comment.
looks good, sorry it took me so long, just a few minor findings from a friend below, nothing blocking
What
Output now works in distributed (MPI) runs. Until now the driver disabled IO with a
warning as soon as
comm_size > 1.Output distributions (
common/io/distributed.py)New strategy layer between the field-group monitor (scheduling) and the writers (format):
GatherDistribution(default): halo entries dropped via the owner masks, ownedentries
Gatherv-ed to rank 0 and placed at their global indices — output identical toa single-rank run, works with both backends, but root-memory-bound.
RankBlockDistribution: no data communication; every rank writes its own ownedentries into a rank-contiguous block of a shared store. Horizontal axes are padded
to a uniform block size;
global_index_<dim>coordinates recover the global order.SingleNodeDistribution: passthrough for single-rank runs.Zarr writer (
common/io/zarr_writers.py)ZarrWriter(zarr format 3) next toNETCDFWriter, both behind aFieldWriterprotocoland sharing all coordinate/CF attribute definitions, so a store reads like the equivalent
netCDF file. Rank-block stores are chunked so that no chunk (or shard) crosses a
rank-block boundary — one chunk per rank block by default — so concurrent writes never
touch the same chunk file; the root rank owns all store-metadata operations, ordered by
one barrier per append.
Parallel netCDF (
common/io/netcdf_writers.py)NETCDFWritersupports the rank-block mode too: on an MPI-parallel netCDF4installation,
distributed+netcdfwrites one shared file (parallel=Trueon themonitor's communicator) with the exact layout of the rank-block zarr store — padded
per-rank blocks, HDF5 chunks never crossing block boundaries,
global_index_<dim>coordinates — so reassembly and any downstream consumer work on either format.
Metadata operations run collectively on every rank; variables touching the unlimited
time dimension use collective access, with every rank writing its full block (padding
included) so ranks owning no entries still participate in every collective write.
Vertical chunk sizes shrink automatically to respect HDF5's 4 GiB chunk limit.
The PyPI wheels of
netCDF4are serial builds (__has_parallel4_support__ == 0), sothe combination is rejected at configuration time on such installations — the error and
the
common.iomodule docs spell out the exact steps to an MPI-parallel build(MPI-enabled netCDF-C/HDF5, then
pip install --no-binary netcdf4 --no-build-isolation --force-reinstall netcdf4withmpi4pypreinstalled, then verify the support flag).The installation is re-verified at writer construction and at file open, and every
parallel open logs the netCDF4/netcdf-c/HDF5 versions in use.
Chunking and sharding (per field group)
horizontal_chunk_size: entries per chunk along the cell/edge/vertex axes, bothbackends. Default: one chunk per rank block (distributed), whole axis / library
default otherwise.
horizontal_shard_size(zarr only): groups whole chunks into one storage file each —the file-count knob on parallel file systems (e.g. one shard file per rank per time
slice with many small chunks inside).
In distributed mode the uniform rank-block size is rounded up to a multiple of the
chunk/shard size (
RankBlockDistributionblock alignment), so chunks and shards nevercross rank-block boundaries and concurrent rank writes stay in disjoint files; the
extra positions are ordinary padding. Misaligned layouts are rejected at writer
construction, on every rank.
Configuration
Per field group:
backend(netcdf|zarr),mode(gather|distributed) andthe chunk/shard sizes above, normalized and validated at construction.
distributed+netcdfrequires an MPI-parallel netCDF4 installation and is rejected otherwise,regardless of rank count. Time-delta output intervals must be a multiple of the model
time step (no silent rounding). Driver exposes
--output-backend/--output-mode(chunk/shard sizes are library-level config for now).
Timings
Output overhead is split into
distributeandwritephases per group and logged as themaximum over the ranks. The driver adds
output_assemble/output_storetimers and a"% of wall time" column to the timer report.
Robustness
Errors raised by the IO layer are raised on all ranks together (overwrite refusal,
partition validation, UGRID export failure, field validation before any file mutation),
so a configuration or data error aborts the job instead of hanging it in the next
collective — and a failed append can never leave a phantom time slice behind.
Tests
rejection), zarr writer (rank-block padding, chunk/shard layouts, misalignment
rejection), netCDF writer rank-block layout and chunking (through a serial handle, so
it runs on any installation), the parallel-support predicate branches, backend/mode
and chunk/shard config validation, monitor-to-writer wiring (rank blocks, chunk/shard
sizes, alignment-aware distribution sharing).
indices, halos, one empty rank) through monitor → distribution → writer with a file
rollover, verified value by value against an analytic field; a sub-chunked + sharded
distributed-zarr case additionally pins the on-disk chunk/shard layout. The
netcdf+distributed case skips with an explicit reason on serial netCDF4 installations.
rank-block zarr must reproduce the single-rank reference; rank-block netCDF joins on
MPI-parallel netCDF4 installations.
MODEL_MPI_SUBSETSset toallso the data-free MPI tests run.Not in this PR (follow-ups)
so the netcdf+distributed MPI cases currently skip there; running them end-to-end
needs an MPI-parallel netCDF4 build in the image.
horizontal_chunk_size/horizontal_shard_sizeare settable onFieldGroupIOConfigbut not yet exposed asdriver CLI options.
storeis still synchronous; the deep copy of the output state isa marked TODO that becomes necessary once writes are off the critical path.
icon4py-gridlookviewer-export CLI has been split out ofthis PR (per review) and lives on the
gridlookbranch; it will come as its own PR.