diff --git a/.gitignore b/.gitignore index 81c0b038..304ae631 100755 --- a/.gitignore +++ b/.gitignore @@ -153,7 +153,9 @@ docs/doxyhtml/ docs/doxyxml/ ### Output #### -## Plot files -**/plt* +## Plot files (e.g., plt00000 and test outputs like plt_read_hdr, +## but not source files such as pltfile_to_openpmd.py) +**/plt[0-9]* +**/plt_* ## Checkpoints **/chk* diff --git a/docs/source/usage/workflows.rst b/docs/source/usage/workflows.rst index fd60f341..b3e2de86 100644 --- a/docs/source/usage/workflows.rst +++ b/docs/source/usage/workflows.rst @@ -10,6 +10,7 @@ This section collects typical user workflows and best practices for pyAMReX. :maxdepth: 2 workflows/read_plotfiles + workflows/pltfile_to_openpmd .. workflows/parallelization .. workflows/profiling .. workflows/debugging diff --git a/docs/source/usage/workflows/pltfile_to_openpmd.rst b/docs/source/usage/workflows/pltfile_to_openpmd.rst new file mode 100644 index 00000000..fdfa4cfb --- /dev/null +++ b/docs/source/usage/workflows/pltfile_to_openpmd.rst @@ -0,0 +1,46 @@ +.. _usage-how-to-pltfile-to-openpmd: + +Convert Plotfiles to openPMD +============================ + +The ``pltfile-to-openpmd`` tool converts native AMReX plotfiles (mesh fields and particles) to `openPMD `__ series, readable by openPMD-viewer, ParaView, VisIt and the wider openPMD ecosystem. +The openPMD backend - HDF5 (``.h5``), ADIOS2 (``.bp``) or JSON (``.json``) - is selected by the output file extension. +It requires the `openpmd-api `__ Python package (``pip install openpmd-api``). + +Each plotfile becomes one openPMD iteration, indexed by its level-0 step number: + +.. code-block:: bash + + pltfile-to-openpmd -o sim_%T.h5 diags/plt00000 diags/plt00100 + + # equivalent, e.g., if the entry point is not on PATH: + python -m amrex.tools.pltfile_to_openpmd -o sim_%T.h5 diags/plt????? + +or, from Python: + +.. code-block:: python + + from amrex.tools.pltfile_to_openpmd import convert + + convert(["diags/plt00000", "diags/plt00100"], "sim_%T.h5") + +Run ``pltfile-to-openpmd --help`` for all options (field/species selection, skipping particles, recording a time step, quiet mode). +The plotfile's dimensionality is detected automatically and the matching ``amrex.space{1,2,3}d`` module is used. + +Data Mapping +------------ + +The conversion is information-preserving: + +* Field data is copied at its on-disk precision, per AMR level, with every AMReX grid stored as one chunk of the level's dataset. + AMReX's Fortran axis order is reversed into openPMD's C order (``axisLabels`` e.g. ``["z", "y", "x"]``). +* Mesh refinement levels follow the openPMD `PatchBasedMeshRefinement `__ extension proposal: the coarsest level keeps the plain record name (readable by every openPMD tool), finer levels are suffixed ``_lvl`` and carry a ``refinementRatio`` attribute. +* Particles are converted per species (discovered via :py:func:`~amrex.space3d.list_particle_species` and read via :py:func:`~amrex.space3d.read_particles`, see :ref:`Read Back Plotfiles `), with their component names verbatim, unpacked ``id`` and ``amrex_cpu`` records, and a constant ``positionOffset`` of zero. +* AMReX metadata without an openPMD equivalent is stored in ``amrex_``-prefixed attributes: per-level steps, box arrays, ghost-cell widths, the coordinate system, and per-species file metadata. + Together with the chunk layout, this suffices to reconstruct the plotfile structure. + +Limitations, by design of the source format and this tool: + +* Plotfiles carry no unit metadata, so ``unitSI`` is 1 and ``unitDimension`` is dimensionless; record a time step with ``--dt`` if needed. +* The on-disk *ordering* of particles is not preserved (identities are, via ``id``/``amrex_cpu``); their mesh-refinement level assignment is recoverable from positions and the stored per-level box arrays - the same rule AMReX applies in ``Redistribute()``. +* Ghost cell *values* are not written (the valid region is); the ghost width is recorded in ``amrex_n_grow``. diff --git a/setup.py b/setup.py index f3d9b0f8..7dce8f61 100644 --- a/setup.py +++ b/setup.py @@ -278,6 +278,14 @@ def build_extension(self, ext): python_requires=">=3.11", tests_require=["pytest"], install_requires=install_requires, + extras_require={ + "openpmd": ["openpmd-api"], + }, + entry_points={ + "console_scripts": [ + "pltfile-to-openpmd = amrex.tools.pltfile_to_openpmd:main", + ], + }, # cmdclass={'test': PyTest}, # platforms='any', classifiers=[ diff --git a/src/amrex/tools/__init__.py b/src/amrex/tools/__init__.py new file mode 100644 index 00000000..319a9391 --- /dev/null +++ b/src/amrex/tools/__init__.py @@ -0,0 +1,2 @@ +# -*- coding: utf-8 -*- +"""Standalone tools built on pyAMReX, usable as modules and CLIs.""" diff --git a/src/amrex/tools/pltfile_to_openpmd.py b/src/amrex/tools/pltfile_to_openpmd.py new file mode 100644 index 00000000..8f8aec7e --- /dev/null +++ b/src/amrex/tools/pltfile_to_openpmd.py @@ -0,0 +1,463 @@ +# -*- coding: utf-8 -*- +""" +Convert AMReX plotfiles to openPMD series. + +This tool reads native AMReX plotfiles (mesh fields and particles) and writes +them as an openPMD series (HDF5, ADIOS2 or JSON, chosen by the output file +extension). It aims to be information-preserving: + +* every field value is copied at its on-disk precision, per AMR level, +* mesh refinement levels are encoded following the openPMD + ``PatchBasedMeshRefinement`` extension proposal + (https://github.com/openPMD/openPMD-standard/pull/252): the coarsest level + keeps the plain record name, finer levels are suffixed ``_lvl`` and carry + a ``refinementRatio`` attribute; every AMReX grid is stored as one chunk, +* particles are converted per species with their component names, ids and + cpus, +* AMReX metadata that has no openPMD equivalent (level steps, box arrays, + ghost cell widths, ...) is stored in ``amrex_``-prefixed attributes. + +Not preserved: the on-disk *ordering* of particles (their identity is kept via +``id``/``cpu``) and their mesh-refinement level assignment, which is +recoverable from the positions and the stored per-level box arrays - the same +rule AMReX itself applies in ``Redistribute()``. Ghost cell *values* are not +written (the valid region is); the ghost width is recorded. + +Usage (CLI):: + + pltfile-to-openpmd -o sim_%T.h5 plt00000 plt00100 ... + python -m amrex.tools.pltfile_to_openpmd -o sim_%T.bp plt????? + +Usage (module):: + + from amrex.tools.pltfile_to_openpmd import convert + + convert(["plt00000", "plt00100"], "sim_%T.h5") + +Requires the ``openpmd_api`` Python package (``pip install openpmd-api``). +""" + +import argparse +import glob +import os +import re + +import numpy as np + + +def _peek_spacedim(plotfile): + """Read the spatial dimension from a plotfile ``Header`` without pyAMReX. + + The header starts with: version string, number of components, that many + component names, then the spatial dimension. + """ + header = os.path.join(plotfile, "Header") + if not os.path.isfile(header): + raise FileNotFoundError( + f"pltfile-to-openpmd: no plotfile Header found at '{header}'. " + f"Is '{plotfile}' an AMReX plotfile directory?" + ) + with open(header) as f: + f.readline() # version + ncomp = int(f.readline()) + for _ in range(ncomp): + f.readline() # component names + return int(f.readline()) + + +def _import_amr(spacedim): + """Import the pyAMReX module matching the plotfile's dimensionality.""" + import importlib + + try: + return importlib.import_module(f"amrex.space{spacedim}d") + except ImportError as e: + raise ImportError( + f"pltfile-to-openpmd: the plotfile is {spacedim}D, but pyAMReX was " + f"built without amrex.space{spacedim}d." + ) from e + + +def _fab_dtype(plotfile, lev): + """Detect the on-disk floating point precision of a level's field data. + + Every FAB in the level's binary data files starts with an ASCII header + ``FAB ((, ...`` - 8 for double, 4 for single precision. + Returns None if the layout is non-standard and detection fails. + """ + for vismf_header in sorted( + glob.glob(os.path.join(plotfile, f"Level_{lev}", "*_H")) + ): + with open(vismf_header) as f: + m = re.search(r"^FabOnDisk: (\S+) (\d+)", f.read(), re.MULTILINE) + if m is None: + continue + fab_file = os.path.join(os.path.dirname(vismf_header), m.group(1)) + try: + with open(fab_file, "rb") as f: + f.seek(int(m.group(2))) + fab = f.read(16) + except OSError: + continue + m = re.match(rb"FAB \(\((\d+),", fab) + if m is not None: + return {4: np.float32, 8: np.float64}.get(int(m.group(1))) + return None + + +def _reversed_list(seq): + """AMReX orders axes Fortran-style (x fastest); openPMD datasets are + written in C order, so all per-axis vectors are reversed.""" + return list(seq)[::-1] + + +def _axis_labels(spacedim, coord_sys): + """Axis labels in AMReX (Fortran) order, before reversal.""" + if int(coord_sys) == 1 and spacedim == 2: # RZ / cylindrical + return ["r", "z"] + return ["x", "y", "z"][:spacedim] + + +def _mesh_name(varname, lev): + return varname if lev == 0 else f"{varname}_lvl{lev}" + + +def _convert_mesh(amr, io, plt, iteration, plotfile, fields=None, verbose=False): + """Write all (selected) field components of all levels of one plotfile.""" + spacedim = plt.spaceDim() + varnames = [str(v) for v in plt.varNames()] + if fields is not None: + unknown = set(fields) - set(varnames) + if unknown: + raise ValueError( + f"pltfile-to-openpmd: unknown field(s) {sorted(unknown)}; " + f"the plotfile contains {varnames}" + ) + varnames = [v for v in varnames if v in fields] + + axis_labels_c = _reversed_list(_axis_labels(spacedim, plt.coordSys())) + geometry = { + 0: io.Geometry.cartesian, + 1: io.Geometry.cylindrical, + 2: io.Geometry.spherical, + }.get(int(plt.coordSys()), io.Geometry.other) + + for lev in range(plt.finestLevel() + 1): + domain = plt.probDomain(lev) + extent_c = _reversed_list( + [int(b) - int(s) + 1 for s, b in zip(domain.small_end, domain.big_end)] + ) + dtype = _fab_dtype(plotfile, lev) + n_grow = plt.nGrowVect(lev) + if verbose and max(n_grow) > 0: + print( + f" note: level {lev} stores {list(n_grow)} ghost cells; " + "their values are not written (the valid region is)" + ) + + for varname in varnames: + mf = plt.get(lev, varname) + mesh = iteration.meshes[_mesh_name(varname, lev)] + mrc = mesh[io.Mesh_Record_Component.SCALAR] + + mesh.geometry = geometry + mesh.axis_labels = axis_labels_c + mesh.data_order = "C" + mesh.grid_spacing = _reversed_list([float(x) for x in plt.cellSize(lev)]) + mesh.grid_global_offset = _reversed_list([float(x) for x in plt.probLo()]) + if lev > 0: + # refinement ratio towards the previous (coarser) level, + # ordered like axis_labels (openPMD-standard PR #252) + ratio = amr.IntVect(plt.refRatio(lev - 1)) + mesh.set_attribute( + "refinementRatio", _reversed_list([int(r) for r in ratio]) + ) + + # in-cell position: 0.5 for cell centers, 0.0 on nodes + ix_type = domain.ix_type + mrc.position = _reversed_list( + [0.0 if ix_type.node_centered(d) else 0.5 for d in range(spacedim)] + ) + + data_dtype = ( + dtype if dtype is not None else mf.array(next(iter(mf))).to_xp().dtype + ) + mrc.reset_dataset(io.Dataset(np.dtype(data_dtype), extent_c)) + + for mfi in mf: + arr = mf.array(mfi).to_xp() # (nx, ny, nz[, ...], ncomp) w/ ghosts + box = mfi.validbox() + lo = [int(s) for s in box.small_end] + hi = [int(b) for b in box.big_end] + # strip ghost cells: array indices start at the grown lower end + sl = tuple( + slice(int(n_grow[d]), int(n_grow[d]) + (hi[d] - lo[d] + 1)) + for d in range(spacedim) + ) + valid = arr[sl + (0,)] + # to C order: reverse the axes, match the dataset precision + chunk = np.ascontiguousarray(valid.transpose()).astype( + data_dtype, copy=False + ) + offset_c = _reversed_list( + [lo[d] - int(domain.small_end[d]) for d in range(spacedim)] + ) + mrc.store_chunk(chunk, offset_c, list(chunk.shape)) + + +def _soa_to_numpy(podvector): + """Copy a PODVector view to a host numpy array.""" + return np.array(podvector.to_numpy(copy=True)) + + +def _convert_particles(amr, io, iteration, plotfile, species=None, verbose=False): + """Write all (selected) particle species of one plotfile.""" + found = amr.list_particle_species(plotfile) + if species is not None: + unknown = set(species) - set(found) + if unknown: + raise ValueError( + f"pltfile-to-openpmd: unknown species {sorted(unknown)}; " + f"the plotfile contains {found}" + ) + found = [s for s in found if s in species] + + spacedim = amr.Config.spacedim + axes = _axis_labels(spacedim, 0) + + for name in found: + if verbose: + print(f" particles: {name}") + header = amr.ParticleHeader.read(plotfile, name) + pc = amr.read_particles(plotfile, name) + np_total = header.num_particles + + sp = iteration.particles[name] + + # gather rank-local tiles (this tool is serial); positions are the + # first AMREX_SPACEDIM SoA real components in a pure-SoA container + idcpu_parts, real_parts, int_parts = [], [], [] + n_real = pc.num_real_comps + n_int = pc.num_int_comps + for lvl in range(pc.finest_level + 1): + for pti in pc.iterator(level=lvl): + soa = pti.soa() + idcpu_parts.append(_soa_to_numpy(soa.get_idcpu_data())) + real_parts.append( + [_soa_to_numpy(soa.get_real_data(j)) for j in range(n_real)] + ) + int_parts.append( + [_soa_to_numpy(soa.get_int_data(j)) for j in range(n_int)] + ) + + def concat(parts, j=None): + arrs = [p if j is None else p[j] for p in parts] + return ( + np.concatenate(arrs) + if arrs + else np.array([], dtype=np.uint64 if j is None else np.float64) + ) + + idcpu = concat(idcpu_parts) + assert idcpu.size == np_total, ( + f"read {idcpu.size} particles, header announces {np_total}" + ) + + def store(record_component, data): + record_component.reset_dataset(io.Dataset(data.dtype, [np_total])) + record_component.store_chunk(np.ascontiguousarray(data), [0], [data.size]) + + # position + constant positionOffset (openPMD base records) + for d, ax in enumerate(axes): + store(sp["position"][ax], concat(real_parts, d)) + poff = sp["positionOffset"][ax] + poff.reset_dataset(io.Dataset(np.dtype(np.float64), [np_total])) + poff.make_constant(0.0) + + # identity: unpacked AMReX id and cpu + store(sp["id"][io.Record_Component.SCALAR], amr.unpack_ids(idcpu)) + store(sp["amrex_cpu"][io.Record_Component.SCALAR], amr.unpack_cpus(idcpu)) + + # named runtime components, verbatim as scalar records + for j, comp in enumerate(header.real_comp_names): + store( + sp[str(comp)][io.Record_Component.SCALAR], + concat(real_parts, spacedim + j), + ) + for j, comp in enumerate(header.int_comp_names): + store(sp[str(comp)][io.Record_Component.SCALAR], concat(int_parts, j)) + + # AMReX metadata: file layout details with no openPMD equivalent + sp.set_attribute("amrex_version", header.version) + sp.set_attribute("amrex_is_checkpoint", int(header.is_checkpoint)) + sp.set_attribute("amrex_next_id", int(header.next_id)) + sp.set_attribute( + "amrex_num_particles_per_level", + [int(sum(e.count for e in entries)) for entries in header.grids], + ) + + +def convert( + plotfiles, + output, + fields=None, + species=None, + no_particles=False, + dt=0.0, + author=None, + verbose=True, +): + """Convert AMReX plotfiles into one openPMD series. + + Parameters + ---------- + plotfiles : list of str + AMReX plotfile directories; each becomes one iteration, indexed by its + level-0 step number. + output : str + openPMD series path; the extension selects the backend (``.h5``, + ``.bp``, ``.json``) and a ``%T`` placeholder selects file-based + iteration encoding. + fields : list of str, optional + Only convert these field components (default: all). + species : list of str, optional + Only convert these particle species (default: all). + no_particles : bool, optional + Skip particle data entirely. + dt : float, optional + Time step to record per iteration; plotfiles do not store one. + author : str, optional + openPMD author attribute, e.g. ``"Jane Doe "``. + verbose : bool, optional + Print progress. + """ + try: + import openpmd_api as io + except ImportError as e: + raise ImportError( + "pltfile-to-openpmd requires the openpmd_api package: " + "https://openpmd-api.readthedocs.io - e.g. 'pip install openpmd-api'" + ) from e + + if not plotfiles: + raise ValueError("pltfile-to-openpmd: no input plotfiles given") + + spacedim = _peek_spacedim(plotfiles[0]) + for p in plotfiles[1:]: + if _peek_spacedim(p) != spacedim: + raise ValueError( + f"pltfile-to-openpmd: '{p}' is not {spacedim}D like " + f"'{plotfiles[0]}'; convert equal-dimension files together" + ) + amr = _import_amr(spacedim) + + initialized_here = False + if not amr.initialized(): + amr.initialize([]) + initialized_here = True + + try: + # open each plotfile, order iterations by ascending step + plts = {} + for p in plotfiles: + plt = amr.PlotFileData(p.rstrip("/")) + step = plt.levelStep(0) + if step in plts: + raise ValueError( + f"pltfile-to-openpmd: '{p}' and '{plts[step][0]}' both have " + f"step {step}; cannot write both into one series" + ) + plts[step] = (p, plt) + + series = io.Series(output, io.Access.create) + series.set_software("pyAMReX", amr.__version__) + if author: + series.author = author + + for step in sorted(plts): + p, plt = plts[step] + if verbose: + print(f"converting {p} -> iteration {step}") + it = series.write_iterations()[step] + it.time = float(plt.time()) + it.dt = float(dt) + it.time_unit_SI = 1.0 + + it.set_attribute("amrex_plotfile_version", "unknown") + it.set_attribute("amrex_finest_level", plt.finestLevel()) + it.set_attribute( + "amrex_level_steps", + [int(plt.levelStep(lev)) for lev in range(plt.finestLevel() + 1)], + ) + it.set_attribute("amrex_coord_sys", int(plt.coordSys())) + it.set_attribute( + "amrex_n_grow", + [ + int(g) + for lev in range(plt.finestLevel() + 1) + for g in plt.nGrowVect(lev) + ], + ) + for lev in range(plt.finestLevel() + 1): + ba = plt.boxArray(lev) + flat = [] + for i in range(ba.size): + b = ba[i] + flat += [int(x) for x in b.small_end] + [int(x) for x in b.big_end] + it.set_attribute(f"amrex_box_array_lvl{lev}", flat) + + _convert_mesh(amr, io, plt, it, p, fields=fields, verbose=verbose) + if not no_particles: + _convert_particles(amr, io, it, p, species=species, verbose=verbose) + + it.close() + + series.close() + finally: + if initialized_here: + amr.finalize() + + if verbose: + print(f"wrote {output}") + + +def main(argv=None): + parser = argparse.ArgumentParser( + prog="pltfile-to-openpmd", + description="Convert AMReX plotfiles to an openPMD series " + "(HDF5/ADIOS2/JSON by output extension).", + ) + parser.add_argument("plotfiles", nargs="+", help="AMReX plotfile directories") + parser.add_argument( + "-o", + "--output", + default="openpmd_%T.h5", + help="output series (default: %(default)s); '%%T' selects file-based " + "iteration encoding", + ) + parser.add_argument("--fields", nargs="+", help="only convert these fields") + parser.add_argument("--species", nargs="+", help="only convert these species") + parser.add_argument( + "--no-particles", action="store_true", help="skip particle data" + ) + parser.add_argument( + "--dt", type=float, default=0.0, help="time step to record (not in plotfiles)" + ) + parser.add_argument("--author", help="openPMD author attribute") + parser.add_argument("-q", "--quiet", action="store_true", help="no progress output") + args = parser.parse_args(argv) + + convert( + args.plotfiles, + args.output, + fields=args.fields, + species=args.species, + no_particles=args.no_particles, + dt=args.dt, + author=args.author, + verbose=not args.quiet, + ) + + +if __name__ == "__main__": + main() diff --git a/tests/test_pltfile_to_openpmd.py b/tests/test_pltfile_to_openpmd.py new file mode 100644 index 00000000..cc6660b9 --- /dev/null +++ b/tests/test_pltfile_to_openpmd.py @@ -0,0 +1,206 @@ +# -*- coding: utf-8 -*- + +import importlib.util +import shutil + +import numpy as np +import pytest + +import amrex.space3d as amr + +has_openpmd = importlib.util.find_spec("openpmd_api") is not None +if has_openpmd: + import openpmd_api as io + + from amrex.tools.pltfile_to_openpmd import convert + +pytestmark = [ + pytest.mark.skipif(amr.Config.spacedim != 3, reason="Requires AMREX_SPACEDIM = 3"), + pytest.mark.skipif(not has_openpmd, reason="Requires openpmd_api"), +] + + +def write_single_level_plotfile(filename): + """32^3 cells in 16^3 boxes, one linear-ramp component.""" + domain_box = amr.Box([0, 0, 0], [31, 31, 31]) + real_box = amr.RealBox([-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]) + geom = amr.Geometry(domain_box, real_box, amr.CoordSys.cartesian, [0, 0, 0]) + + ba = amr.BoxArray(domain_box) + ba.max_size(16) + dm = amr.DistributionMapping(ba, 1) + mf = amr.MultiFab(ba, dm, 1, 0) + # a position-dependent value, so chunk placement errors are caught + for mfi in mf: + bx = mfi.tilebox() + marr = mf.array(mfi).to_xp() + i_s, j_s, k_s = tuple(bx.small_end) + nx, ny, nz, _ = marr.shape + i = np.arange(i_s, i_s + nx)[:, None, None] + j = np.arange(j_s, j_s + ny)[None, :, None] + k = np.arange(k_s, k_s + nz)[None, None, :] + marr[..., 0] = i + 100 * j + 10000 * k + + amr.write_single_level_plotfile( + filename, mf, amr.Vector_string(["ramp"]), geom, 1.5, 300 + ) + return mf + + +def test_convert_single_level_mesh(tmp_path): + plt_name = str(tmp_path / "plt00300") + write_single_level_plotfile(plt_name) + + out = str(tmp_path / "series_%T.h5") + convert([plt_name], out, verbose=False) + + series = io.Series(out, io.Access.read_only) + it = series.iterations[300] + assert np.isclose(it.time, 1.5) + + mesh = it.meshes["ramp"] + assert mesh.geometry == io.Geometry.cartesian + assert mesh.axis_labels == ["z", "y", "x"] + np.testing.assert_allclose(mesh.grid_spacing, [1.0 / 32.0] * 3) + np.testing.assert_allclose(mesh.grid_global_offset, [-0.5] * 3) + + mrc = mesh[io.Mesh_Record_Component.SCALAR] + assert list(mrc.shape) == [32, 32, 32] + data = mrc.load_chunk() + series.flush() + + # data is (z, y, x); rebuild the expected ramp + k, j, i = np.meshgrid(np.arange(32), np.arange(32), np.arange(32), indexing="ij") + np.testing.assert_array_equal(data, i + 100 * j + 10000 * k) + + # lossless-reconstruction metadata (1-element arrays may read as scalars) + assert list(np.atleast_1d(it.get_attribute("amrex_level_steps"))) == [300] + assert it.get_attribute("amrex_coord_sys") == 0 + ba_flat = it.get_attribute("amrex_box_array_lvl0") + assert len(ba_flat) == 8 * 6 # 8 boxes, small_end+big_end each + + series.close() + + +def test_convert_multi_level_mesh(tmp_path): + """Two-level hierarchy: level 1 refines the upper-x half of the domain.""" + plt_name = str(tmp_path / "plt00007") + + real_box = amr.RealBox([-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]) + domain0 = amr.Box([0, 0, 0], [15, 15, 15]) + geom0 = amr.Geometry(domain0, real_box, amr.CoordSys.cartesian, [0, 0, 0]) + domain1 = amr.Box([0, 0, 0], [31, 31, 31]) + geom1 = amr.Geometry(domain1, real_box, amr.CoordSys.cartesian, [0, 0, 0]) + + ba0 = amr.BoxArray(domain0) + ba0.max_size(8) + ba1 = amr.BoxArray(amr.Box([16, 0, 0], [31, 31, 31])) # refined upper-x half + ba1.max_size(16) + + mfs = [] + for ba, val in ((ba0, 1.0), (ba1, 2.0)): + dm = amr.DistributionMapping(ba, 1) + mf = amr.MultiFab(ba, dm, 1, 0) + mf.set_val(val) + mfs.append(mf) + + amr.write_multi_level_plotfile( + plt_name, mfs, ["density"], [geom0, geom1], 0.25, [7, 7], [amr.IntVect(2)] + ) + + out = str(tmp_path / "ml_%T.h5") + convert([plt_name], out, verbose=False) + + series = io.Series(out, io.Access.read_only) + it = series.iterations[7] + + # level 0: plain name, full domain + m0 = it.meshes["density"][io.Mesh_Record_Component.SCALAR] + assert list(m0.shape) == [16, 16, 16] + d0 = m0.load_chunk() + series.flush() + np.testing.assert_array_equal(d0, 1.0) + + # level 1: _lvl1 suffix, refinementRatio, sized to the refined index space + m1 = it.meshes["density_lvl1"] + assert m1.get_attribute("refinementRatio") == [2, 2, 2] + mrc1 = m1[io.Mesh_Record_Component.SCALAR] + assert list(mrc1.shape) == [32, 32, 32] + np.testing.assert_allclose(m1.grid_spacing, [1.0 / 32.0] * 3) # refined spacing + # only the covered half is defined: read it back chunk-wise + d1 = mrc1.load_chunk([0, 0, 16], [32, 32, 16]) + series.flush() + np.testing.assert_array_equal(d1, 2.0) + + series.close() + + +def test_convert_particles(tmp_path): + plt_name = str(tmp_path / "plt00042") + n_part = 21 + + domain_box = amr.Box([0, 0, 0], [31, 31, 31]) + real_box = amr.RealBox([-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]) + geom = amr.Geometry(domain_box, real_box, amr.CoordSys.cartesian, [0, 0, 0]) + ba = amr.BoxArray(domain_box) + dm = amr.DistributionMapping(ba, 1) + + mf = amr.MultiFab(ba, dm, 1, 0) + mf.set_val(0.0) + amr.write_single_level_plotfile( + plt_name, mf, amr.Vector_string(["dummy"]), geom, 0.0, 42 + ) + + pc = amr.ParticleContainer_pureSoA_3_0_polymorphic(geom, dm, ba) + pc.arena = amr.The_Arena() + myt = amr.ParticleInitType_pureSoA_3_0() + myt.real_array_data = [0.0, 0.0, 0.0] + myt.int_array_data = [] + pc.init_random(n_part, 7, myt, False, real_box) + pc.add_real_comp("w", True) + pc.add_int_comp("tag", True) + for pti in pc.iterator(level=0): + soa = pti.soa() + soa.get_real_data(3).assign(3.25) + soa.get_int_data(0).assign(9) + pc.redistribute() + pc.write_plotfile( + plt_name, "electrons", amr.Vector_string(["w"]), amr.Vector_string(["tag"]) + ) + + out = str(tmp_path / "parts_%T.h5") + convert([plt_name], out, verbose=False) + + series = io.Series(out, io.Access.read_only) + sp = series.iterations[42].particles["electrons"] + + x = sp["position"]["x"].load_chunk() + w = sp["w"][io.Record_Component.SCALAR].load_chunk() + tag = sp["tag"][io.Record_Component.SCALAR].load_chunk() + pid = sp["id"][io.Record_Component.SCALAR].load_chunk() + series.flush() + + assert x.size == n_part + assert np.all((x >= -0.5) & (x <= 0.5)) + np.testing.assert_allclose(w, 3.25) + np.testing.assert_array_equal(tag, 9) + assert np.unique(pid).size == n_part # ids are unique particle identities + + assert list(np.atleast_1d(sp.get_attribute("amrex_num_particles_per_level"))) == [ + n_part + ] + + series.close() + + +def test_convert_field_selection_errors(tmp_path): + plt_name = str(tmp_path / "plt00300") + write_single_level_plotfile(plt_name) + + with pytest.raises(ValueError, match="unknown field"): + convert([plt_name], str(tmp_path / "e_%T.h5"), fields=["nope"], verbose=False) + + with pytest.raises(FileNotFoundError, match="plotfile Header"): + convert([str(tmp_path / "not_a_plotfile")], str(tmp_path / "f_%T.h5")) + + shutil.rmtree(plt_name)