diff --git a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md index 929301078..2a07b02fc 100644 --- a/docs/developer/subsystems/checkpoint-output-and-reload-methods.md +++ b/docs/developer/subsystems/checkpoint-output-and-reload-methods.md @@ -11,17 +11,34 @@ but new code should use `write_timestep(..., petsc_reload=True)`. ## Standard API -`write_timestep()` always writes the mesh file and one HDF5 file per mesh -variable. Mesh-variable files always contain raw coordinate/value datasets under -`/fields`, which are the source data used by `MeshVariable.read_timestep()` for -coordinate/KDTree remapping. +`write_timestep()` writes the mesh file and one HDF5 file per mesh variable. +With XDMF enabled, mesh-variable files contain dimensional +coordinates and values under `/fields`. These are the authoritative arrays for +analysis and the source data used by `MeshVariable.read_timestep()` for +coordinate remapping. The two optional payloads are selected with explicit flags: | Flag | Output payload | Reader/use case | | --- | --- | --- | -| `create_xdmf=True` | `/vertex_fields` or `/cell_fields` compatibility datasets plus a companion `.xdmf` file | ParaView/XDMF visualisation | -| `petsc_reload=True` | PETSc DMPlex section/vector metadata under `/topologies/uw_mesh/dms/...` | `MeshVariable.read_checkpoint()` PETSc-native reload | +| `create_xdmf=True` | Dimensional `/fields`, compact high-order reductions when needed, and a companion `.xdmf` file | Analysis and ParaView/XDMF visualisation | +| `petsc_reload=True` | Native PETSc DMPlex section/global-vector data under `/restart/petsc` | `MeshVariable.read_checkpoint()` exact reload | + +The XDMF storage choice follows the finite-element layout: + +| Field layout | XDMF representation | +| --- | --- | +| Continuous P1 | Direct `/fields` node values | +| Continuous P2 triangles/tetrahedra | Direct `/fields` values with `Triangle_6`/`Tetrahedron_10` connectivity | +| DG0 | Direct `/fields` cell values | +| DG1 triangles/tetrahedra | Native `/fields` for reload plus exact disconnected-corner `/visualization` data | +| Continuous P3+ or unsupported P2 layout | Compact P1 dataset under `/visualization` | +| DG2+ or unsupported DG1 | Compact DG0 dataset under `/visualization` | + +DG1 retains native interpolation coordinates and values under `/fields`, so +`read_timestep()` can reload the solver field. XDMF reads an additional exact +basis conversion at disconnected element corners from `/visualization` because +the native interior DG1 points do not describe the full element geometry. ### Visualisation And Remap @@ -48,11 +65,13 @@ output.mesh.Pressure.00000.h5 output.mesh.00000.xdmf ``` -The field files contain coordinate/value datasets such as `/fields/` and -`/fields/coordinates`, plus vertex-field datasets for visualisation. Reloading -uses coordinate-based remapping. In practice this means the target variable is -filled by comparing target coordinates to source coordinates, using a KDTree or -similar nearest-neighbour/remap process. +The field files contain `/fields/` and `/fields/coordinates`. P1, P2 +triangles/tetrahedra, and DG0 are visualized directly from those datasets. DG1 +uses its additional exact disconnected-corner representation. Continuous P3+ +fields use one compact P1 visualization reduction, and DG2+ fields use DG0. +Reloading with `read_timestep()` compares target coordinates to the dimensional +source coordinates and converts the saved values back to the active model's +nondimensional solver frame. ### Unified Visualisation And PETSc Reload @@ -75,9 +94,9 @@ velocity.read_checkpoint( ) ``` -With both `create_xdmf=True` and `petsc_reload=True`, the same variable file can -be used by `read_timestep()` for coordinate/KDTree remapping and by -`read_checkpoint()` for exact PETSc-native reload. +With both flags enabled, the same variable file contains dimensional `/fields` +for analysis, XDMF, and coordinate remapping plus native `/restart/petsc` data +for exact PETSc reload. ### PETSc Reload Without XDMF @@ -94,8 +113,8 @@ mesh.write_timestep( ) ``` -This still writes raw `/fields` datasets, but it does not write -`/vertex_fields`, `/cell_fields`, or a companion `.xdmf` file. +This writes the native checkpoint directly and does not create a companion +`.xdmf` file or duplicate `/fields` coordinate/value arrays. Typical PETSc-reload-only files still use the timestep naming convention: @@ -105,12 +124,17 @@ restart.mesh.Velocity.00000.h5 restart.mesh.Pressure.00000.h5 ``` -The variable files contain raw `/fields` datasets and PETSc reload metadata -under `/topologies/uw_mesh/dms//`. +Each variable file contains only PETSc reload metadata and one native global +vector under `/restart/petsc/topologies/uw_mesh/dms//`. + +`read_checkpoint()` also recognizes the former `/uw_checkpoint` group so +existing restart files remain readable. ### Advantages -- Produces XDMF/HDF5 files suitable for visualisation workflows. +- Produces dimensional HDF5 fields suitable for analysis and visualisation. +- Avoids duplicate P1 field arrays when the finite-element layout can be + represented directly. - Can remap data onto a different mesh or a different node layout. - Useful for postprocessing where exact finite-element section identity is not required. diff --git a/docs/developer/subsystems/checkpointing-system.md b/docs/developer/subsystems/checkpointing-system.md index c0502b894..49b0a6c40 100644 --- a/docs/developer/subsystems/checkpointing-system.md +++ b/docs/developer/subsystems/checkpointing-system.md @@ -18,7 +18,9 @@ including registered meshes, variables, swarms, and Python-side state bearers. `Mesh.write_timestep()` is the standard mesh and mesh-variable output method. It writes one mesh HDF5 file and one HDF5 file per requested mesh variable. -Each variable file always contains `/fields` coordinate/value datasets used by +With `create_xdmf=True`, each variable file contains dimensional +`/fields/coordinates` and `/fields/` datasets. These arrays are the +authoritative analysis output and are also used by `MeshVariable.read_timestep()`. Optional payloads are controlled by explicit flags: @@ -26,17 +28,44 @@ Optional payloads are controlled by explicit flags: | Flag | Payload | Reader / use | | --- | --- | --- | | `create_xdmf=True` | XDMF-compatible visualisation datasets and a companion `.xdmf` file | ParaView and other XDMF tools | -| `petsc_reload=True` | PETSc DMPlex section/vector metadata | `MeshVariable.read_checkpoint()` | +| `petsc_reload=True` | PETSc DMPlex section/global-vector metadata | `MeshVariable.read_checkpoint()` | + +When nondimensional scaling is active, `/fields` is converted during the write +to the mesh and variable units declared in the model. HDF5 attributes and XDMF +`Information` elements record those units. Analysis scripts can therefore read +physical values directly, without maintaining their own conversion table. + +Set `petsc_reload=True` only when an exact restart is needed. It adds the native +nondimensional PETSc payload under `/restart/petsc`; the visualization and +analysis arrays remain dimensional. ### Visualisation and Coordinate Remap +XDMF reads P1 and DG0 values directly from `/fields`. Continuous P2 fields on +triangles and tetrahedra use XDMF `Triangle_6` and `Tetrahedron_10` +connectivity, including their edge nodes, so no P1 projection is stored. DG1 +keeps native interpolation coordinates and values under `/fields`, allowing +`read_timestep()` to recover the solver field. For XDMF, the same element +polynomial is evaluated at disconnected cell corners under `/visualization`, +preserving jumps without averaging traces across shared edges or faces. + +XDMF cannot represent every UW3 finite-element layout directly. Continuous P3+ +fields and unsupported P2 layouts receive one compact P1 dataset under +`/visualization`. DG2+ fields and unsupported DG1 layouts receive one compact +DG0 dataset. Their exact dimensional values remain under `/fields`. Integration +point fields are not supported by this writer. + +DG1 therefore has two representations because their coordinate sets serve +different purposes: native `/fields` for analysis and coordinate reload, and an +exact disconnected-corner `/visualization` representation for XDMF. Add +`/restart/petsc` when PETSc-native restart is also required. + ```python mesh.write_timestep( "output", index=100, outputPath="output", meshVars=[velocity, pressure, temperature], - time=100.0, create_xdmf=True, ) ``` @@ -96,10 +125,14 @@ output/restart.mesh.velocity.00100.h5 output/restart.mesh.pressure.00100.h5 ``` -The variable files contain `/fields` datasets and PETSc reload metadata under -`/topologies/uw_mesh/dms//`. `read_checkpoint()` uses PETSc DMPlex -topology, section, vector, and `PetscSF` metadata. It does not use KDTree -remapping. +The variable files contain PETSc reload metadata and one native global vector +under `/restart/petsc/topologies/uw_mesh/dms//`. `read_checkpoint()` +uses PETSc DMPlex topology, section, vector, and `PetscSF` metadata. It does not +use dimensional `/fields` values or KDTree remapping. Restart-only output does +not write `/fields`, so the native values are stored only once. + +The reader also accepts the former `/uw_checkpoint` group for compatibility +with existing restart files. ### Unified Visualisation and PETSc Reload diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 7f045997e..b1089c084 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -4832,31 +4832,35 @@ def write_timestep( - one mesh HDF5 file, shared across timesteps unless ``meshUpdates=True`` - one HDF5 file per mesh variable - - raw coordinate/value datasets under ``/fields`` for coordinate-based - reload with ``MeshVariable.read_timestep()`` - - The optional payloads are controlled explicitly: - - - ``create_xdmf=True`` writes ParaView/XDMF output. Variable files also - receive ``/vertex_fields`` or ``/cell_fields`` compatibility groups, - and rank 0 writes the companion ``.xdmf`` file. - - ``petsc_reload=True`` writes PETSc DMPlex section/local-vector - metadata and an in-place global-vector payload into the same - per-variable HDF5 files. These files can then be loaded with - ``MeshVariable.read_checkpoint()`` for exact restart. + The variable-file payloads are controlled explicitly: + + - ``create_xdmf=True`` writes a companion XDMF file that reads the + dimensional ``/fields`` datasets directly for P1, P2 triangles and + tetrahedra, DG0, and simplex DG1. All simplex DG1 variables share one + disconnected-corner geometry in the mesh file; each variable stores + one exact corner-nodal value array. Continuous P3+ fields receive one + compact P1 visualization dataset; discontinuous DG2+ fields receive + DG0. + - ``petsc_reload=True`` writes native nondimensional PETSc DMPlex + section/global-vector data under ``/restart/petsc``. + Load that optional payload with ``MeshVariable.read_checkpoint()`` + for an exact solver restart. When XDMF is disabled, this is the only + variable payload; no duplicate ``/fields`` values are written. + - With both flags disabled, the established low-level native + ``/fields`` output is retained for compatibility. Common choices are: - visualisation/remap only: ``create_xdmf=True, petsc_reload=False`` - - PETSc-native reload only: + - PETSc-native reload only (no duplicate ``/fields`` values): ``create_xdmf=False, petsc_reload=True`` - unified visualisation/remap and PETSc reload: ``create_xdmf=True, petsc_reload=True`` - With both flags enabled, the same variable HDF5 file can be used by - ``MeshVariable.read_timestep()`` for coordinate/KDTree remapping and by - ``MeshVariable.read_checkpoint()`` for exact PETSc-native reload. + With both flags enabled, the same variable HDF5 file supports + dimensional analysis, XDMF visualization, coordinate remapping, and + exact PETSc-native restart without duplicate P1 field copies. Parameters ---------- @@ -4876,12 +4880,25 @@ def write_timestep( If ``False``, reuse ``.mesh.00000.h5`` when it already exists. If ``True``, write an indexed mesh file for this timestep. create_xdmf - Write ParaView/XDMF-compatible datasets and companion XDMF file. + Write ParaView/XDMF-compatible dimensional datasets and the + companion XDMF file. DG1 on full-dimensional triangles/tetrahedra + uses independent vertices per cell, preserving jumps without + smoothing. Continuous P2 triangles and tetrahedra use their exact + quadratic connectivity. Higher-order fields use compact P1 or DG0 + visualization reductions while their exact dimensional values + remain in ``/fields``. petsc_reload Write PETSc DMPlex section/vector metadata for reload with ``MeshVariable.read_checkpoint()``. """ + if create_xdmf: + for var in meshVars or []: + if getattr(var, "is_integration_point", False): + raise NotImplementedError( + "Integration-point XDMF output is not supported; use " + "create_xdmf=False for native-only output." + ) options = PETSc.Options() options.setValue("viewer_hdf5_sp_output", True) options.setValue("viewer_hdf5_collective", False) @@ -4910,24 +4927,36 @@ def write_timestep( mesh_file = output_base_name + f".mesh.{index:05}.h5" self.write(mesh_file) - variables = [] + if create_xdmf: + _write_visualisation_geometry(self, mesh_file) + + checkpoint_variables = [] if meshVars is not None: for var in meshVars: save_location = output_base_name + f".mesh.{var.clean_name}.{index:05}.h5" - var.write(save_location) if create_xdmf: - _write_compat_groups(self, var, save_location) - variables.append((var, save_location)) + var.write(save_location) + if petsc_reload: + self._write_petsc_reload_file(save_location, [var], mode="a") + _write_xdmf_field(self, var, save_location, mesh_file) + elif petsc_reload: + # Exact-restart-only output needs the PETSc section/vector + # payload, not a second native /fields copy of the values. + self._write_petsc_reload_file(save_location, [var], mode="w") + else: + # Preserve the established low-level native output when no + # optional visualization or restart payload was requested. + var.write(save_location) if swarmVars is not None: for svar in swarmVars: save_location = output_base_name + f".proxy.{svar.clean_name}.{index:05}.h5" svar.write_proxy(save_location) if petsc_reload: - variables.append((svar._meshVar, save_location)) + checkpoint_variables.append((svar._meshVar, save_location)) if petsc_reload: - for var, save_location in variables: + for var, save_location in checkpoint_variables: self._write_petsc_reload_file(save_location, [var], mode="a") if create_xdmf and uw.mpi.rank == 0: @@ -4952,9 +4981,8 @@ def petsc_save_checkpoint( This is a convenience wrapper around ``write_timestep()`` that provides the simpler interface used by earlier Underworld3 code. - Output uses the same per-variable file layout and XDMF generation - (including vertex/cell compatibility groups, field projection, and - tensor repacking) as ``write_timestep()``. + Output uses the same compact dimensional ``/fields`` layout and XDMF + generation as ``write_timestep()``. Parameters ---------- @@ -5005,22 +5033,26 @@ def _write_petsc_reload_variable(self, viewer, var): if var._lvec is None: var._set_vec(available=True) + var._sync_lvec_to_gvec() iset, subdm = self.dm.createSubDM(var.field_id) subdm.setName(var.clean_name) old_lvec_name = var._lvec.getName() + old_gvec_name = var._gvec.getName() try: var._lvec.setName(var.clean_name) + var._gvec.setName(var.clean_name) self.dm.sectionView(viewer, subdm) - self.dm.localVectorView(viewer, subdm, var._lvec) + self.dm.globalVectorView(viewer, subdm, var._gvec) finally: var._lvec.setName(old_lvec_name) + var._gvec.setName(old_gvec_name) iset.destroy() subdm.destroy() def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): - """Write DMPlex reload metadata and in-place vector payloads.""" + """Write compact DMPlex reload metadata and native global vectors.""" old_dm_name = self.dm.getName() self.dm.setName("uw_mesh") @@ -5029,7 +5061,12 @@ def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): checkpoint_file, mode, comm=PETSc.COMM_WORLD ) viewer.pushFormat(PETSc.Viewer.Format.HDF5_PETSC) + viewer.pushGroup("/restart/petsc") try: + # PETSc needs the complete source section to construct the + # migration SF when the checkpoint is read with a different MPI + # ownership ordering. This is metadata only; owned field values are + # stored once in each variable's global-vector payload. self.dm.sectionView(viewer, self.dm) for var in variables: @@ -5037,28 +5074,12 @@ def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"): uw.mpi.barrier() finally: + viewer.popGroup() viewer.popFormat() viewer.destroy() if old_dm_name is not None: self.dm.setName(old_dm_name) - viewer = PETSc.ViewerHDF5().create( - checkpoint_file, "a", comm=PETSc.COMM_WORLD - ) - try: - viewer.pushGroup("/uw_checkpoint") - for var in variables: - var._sync_lvec_to_gvec() - checkpoint_vec = PETSc.Vec().createWithArray( - var._gvec.array_r, comm=PETSc.COMM_WORLD - ) - checkpoint_vec.setName(var.clean_name) - viewer(checkpoint_vec) - checkpoint_vec.destroy() - viewer.popGroup() - finally: - viewer.destroy() - @timing.routine_timer_decorator def write_checkpoint( self, @@ -5104,8 +5125,8 @@ def write_checkpoint( into one file: ``.checkpoint..h5``. create_xdmf If ``True``, route through ``write_timestep()`` and write XDMF, - vertex/cell compatibility groups, coordinate/KDTree remap data, - and PETSc reload metadata. The output uses the timestep filename + dimensional field/remap data, and PETSc reload metadata. The output + uses the timestep filename convention ``.mesh...h5``. This mode does not support ``unique_id=True`` or ``separate_variable_files=False``. """ @@ -9626,16 +9647,15 @@ def remesh(self, metric_field, verbose=False): return -def _write_compat_groups(mesh, var, var_h5_path): - """Write ``/vertex_fields/`` or ``/cell_fields/`` compatibility groups. +def _write_xdmf_field(mesh, var, var_h5_path, mesh_h5_path): + """Replace native remap arrays with dimensional field output for XDMF. - Uses ``uw.function.write_vertices_to_viewer`` (PETSc interpolation + - ViewerHDF5) for continuous variables, and - ``uw.function.write_cell_field_to_viewer`` for cell/DG-0 variables. - PETSc handles all parallel I/O natively. - - Vertex coordinates are also written to ``/vertex_fields/coordinates`` - for XDMF compatibility. + P1, P2 triangles/tetrahedra, DG0, and simplex DG1 are represented directly + under ``/fields``. DG1 uses exact element-corner values and shares one + disconnected geometry under ``/viz/dg1`` in the mesh file. Unsupported + higher-order layouts retain exact physical values under ``/fields`` and + receive one compact visualization reduction: continuous fields use P1 and + discontinuous fields use DG0. Parameters ---------- @@ -9645,37 +9665,186 @@ def _write_compat_groups(mesh, var, var_h5_path): The variable whose data has already been written to *var_h5_path* by ``var.write()`` (so ``var._gvec`` is up-to-date). var_h5_path : str - Path to the HDF5 file (already contains ``/fields/``). + Path to the HDF5 file. Native restart data, when requested, has already + been written under ``/restart/petsc``. + mesh_h5_path : str + Path to the shared mesh HDF5 file that owns DG1 visualization geometry. """ + import h5py import underworld3 as uw + from underworld3.function.field_projection import ( + _dg1_corner_data, + _physical_visualisation_values, + _write_index_array_to_group, + _write_vec_to_group, + write_field_coordinates_to_viewer, + write_field_to_viewer, + write_p2_simplex_topology_to_viewer, + write_projected_field_to_viewer, + ) - is_cell = (not var.continuous) or (var.degree == 0) - group = "cell_fields" if is_cell else "vertex_fields" - - # Some PETSc versions (3.21+) write /vertex_fields/ or /cell_fields/ - # automatically during var.write(). Remove any pre-existing group so - # that our compat writer can create it afresh (otherwise PETSc error 76 - # on duplicate dataset). - import h5py + direct_p2 = ( + var.continuous + and var.degree == 2 + and mesh.isSimplex + and mesh.dim in (2, 3) + and mesh.cdim == mesh.dim + ) + direct_dg1 = ( + not var.continuous + and var.degree == 1 + and mesh.isSimplex + and mesh.dim in (2, 3) + and mesh.cdim == mesh.dim + ) + needs_projection = (var.continuous and var.degree > 2) or ( + var.continuous and var.degree == 2 and not direct_p2 + ) or (not var.continuous and var.degree > 0 and not direct_dg1) if uw.mpi.rank == 0: with h5py.File(var_h5_path, "a") as f: - if group in f: - del f[group] + for group in ("fields", "vertex_fields", "cell_fields", "dg1", "visualization"): + if group in f: + del f[group] uw.mpi.barrier() viewer = PETSc.ViewerHDF5().create( - var_h5_path, "a", comm=PETSc.COMM_WORLD, + var_h5_path, + "a", + comm=PETSc.COMM_WORLD, ) - if is_cell: - uw.function.write_cell_field_to_viewer(var, viewer) + if direct_dg1: + corner_rows, corner_values, corner_cells = _dg1_corner_data( + var, repack_tensors=False + ) + _write_vec_to_group( + viewer, corner_values, var.clean_name, "/fields", PETSc.COMM_WORLD + ) else: - uw.function.write_vertices_to_viewer(var, viewer) - uw.function.write_coordinates_to_viewer(mesh, viewer) + write_field_to_viewer(var, viewer, "/fields", var.clean_name) + write_field_coordinates_to_viewer(var, viewer, "/fields") + if direct_p2: + write_p2_simplex_topology_to_viewer(var, viewer, group="/fields") + elif needs_projection: + target_degree = 1 if var.continuous else 0 + write_projected_field_to_viewer( + var, + viewer, + target_degree=target_degree, + continuous=var.continuous, + group="/visualization", + name=var.clean_name, + ) + + viewer.destroy() + + if direct_dg1: + if uw.mpi.rank == 0: + with h5py.File(mesh_h5_path, "r") as mesh_handle: + has_dg1_geometry = "viz/dg1/coordinates" in mesh_handle + else: + has_dg1_geometry = None + has_dg1_geometry = uw.mpi.comm.bcast(has_dg1_geometry, root=0) + if not has_dg1_geometry: + mesh_viewer = PETSc.ViewerHDF5().create( + mesh_h5_path, + "a", + comm=PETSc.COMM_WORLD, + ) + _write_vec_to_group( + mesh_viewer, + corner_rows, + "coordinates", + "/viz/dg1", + PETSc.COMM_WORLD, + ) + _write_index_array_to_group( + mesh_viewer, + corner_cells, + "cells", + "/viz/dg1", + PETSc.COMM_WORLD, + ) + mesh_viewer.destroy() + + _, field_units = _physical_visualisation_values(numpy.ones(1), var.units) + _, coordinate_units = _physical_visualisation_values(numpy.ones(1), mesh.units) + with uw.selective_ranks(0) as should_execute: + if should_execute: + with h5py.File(var_h5_path, "a") as handle: + handle.attrs["storage_frame"] = "physical" + field = handle[f"fields/{var.clean_name}"] + field.attrs["units"] = field_units or "dimensionless" + field.attrs["storage_frame"] = "physical" + field.attrs["degree"] = var.degree + field.attrs["continuous"] = var.continuous + field.attrs["representation"] = ( + "dg1_corner_nodal" if direct_dg1 else "exact" + ) + if direct_dg1: + handle.attrs["dg1_mesh_file"] = os.path.basename(mesh_h5_path) + handle.attrs["dg1_geometry_path"] = "/viz/dg1" + else: + coordinates = handle["fields/coordinates"] + coordinates.attrs["units"] = coordinate_units or "dimensionless" + coordinates.attrs["storage_frame"] = "physical" + if needs_projection: + projected = handle[f"visualization/{var.clean_name}"] + projected.attrs["units"] = field_units or "dimensionless" + projected.attrs["source_degree"] = var.degree + projected.attrs["visualization_degree"] = 1 if var.continuous else 0 + projected.attrs["representation"] = "projection" + if direct_dg1: + with h5py.File(mesh_h5_path, "a") as mesh_handle: + coordinates = mesh_handle["viz/dg1/coordinates"] + coordinates.attrs["units"] = coordinate_units or "dimensionless" + coordinates.attrs["storage_frame"] = "physical" + mesh_handle["viz/dg1/cells"].attrs["representation"] = ( + "disconnected_simplex" + ) + uw.mpi.barrier() + +def _write_visualisation_geometry(mesh, mesh_h5_path): + """Write physical mesh coordinates without changing restart geometry.""" + import h5py + import underworld3 as uw + from underworld3.function.field_projection import ( + _physical_visualisation_enabled, + _physical_visualisation_values, + ) + + if not _physical_visualisation_enabled(mesh.units): + return + + with uw.selective_ranks(0) as should_execute: + if should_execute: + with h5py.File(mesh_h5_path, "a") as handle: + if "viz/geometry" in handle: + del handle["viz/geometry"] + uw.mpi.barrier() + + viewer = PETSc.ViewerHDF5().create(mesh_h5_path, "a", comm=PETSc.COMM_WORLD) + uw.function.write_coordinates_to_viewer( + mesh, + viewer, + group="/viz/geometry", + name="vertices", + ) viewer.destroy() + _, unit_label = _physical_visualisation_values(numpy.ones(1), mesh.units) + with uw.selective_ranks(0) as should_execute: + if should_execute: + with h5py.File(mesh_h5_path, "a") as handle: + handle.attrs["checkpoint_units"] = "nondimensional" + handle.attrs["visualisation_units"] = unit_label + handle["geometry/vertices"].attrs["units"] = "nondimensional" + physical = handle["viz/geometry/vertices"] + physical.attrs["units"] = unit_label + uw.mpi.barrier() + def checkpoint_xdmf( filename: str, @@ -9687,6 +9856,7 @@ def checkpoint_xdmf( import h5py import os import warnings + from xml.sax.saxutils import escape """Create xdmf file for checkpoints""" @@ -9722,6 +9892,7 @@ def checkpoint_xdmf( ) vertices = geom["vertices"] + geometry_units = vertices.attrs.get("units") numVertices = vertices.shape[0] spaceDim = vertices.shape[1] cells = topo["cells"] @@ -9763,6 +9934,17 @@ def checkpoint_xdmf( h5.close() + def units_information(units, indent): + """Return an XDMF Information element for an optional unit label.""" + if units is None: + return "" + if isinstance(units, bytes): + units = units.decode() + value = escape(str(units), {'"': """}) + return f'\n{indent}' + + geometry_information = units_information(geometry_units, " ") + # We only use a subset of the possible cell types if spaceDim == 2: if numCorners == 3: @@ -9802,6 +9984,33 @@ def checkpoint_xdmf( header += """ ]>""" + def direct_p2(var): + return ( + var.continuous + and var.degree == 2 + and var.mesh.isSimplex + and var.mesh.dim in (2, 3) + and var.mesh.cdim == var.mesh.dim + ) + + def direct_dg1(var): + return ( + not var.continuous + and var.degree == 1 + and var.mesh.isSimplex + and var.mesh.dim in (2, 3) + and var.mesh.cdim == var.mesh.dim + ) + + p2_vars = [var for var in meshVars if direct_p2(var)] + dg1_vars = [var for var in meshVars if direct_dg1(var)] + special_vars = p2_vars + dg1_vars + collection_start = ( + '' + f'" if special_vars else "" xdmf_start = f""" @@ -9818,58 +10027,61 @@ def checkpoint_xdmf( &MeshData;:/{geomPath}/vertices + {collection_start} - - /Xdmf/Domain/DataItem[@Name="cells"] + + &MeshData;:/{topoPath}/cells - - /Xdmf/Domain/DataItem[@Name="vertices"] - + + &MeshData;:/{geomPath}/vertices + {geometry_information} """ ## The mesh Var attributes - def get_field_info(h5_filename, mesh_var, center): - """ - Return (num_items, num_components, dataset_path) for a mesh variable. - Prefers vertex/cell compatibility groups, falls back to /fields layout. - """ - compat_name = f"{mesh_var.clean_name}_{mesh_var.clean_name}" - candidates = [] - - if center == "Cell": - candidates = [f"cell_fields/{compat_name}", f"fields/{mesh_var.clean_name}"] - else: - candidates = [f"vertex_fields/{compat_name}", f"fields/{mesh_var.clean_name}"] - + def get_field_info(h5_filename, dataset_path): + """Return item/component counts and units for one stored field.""" with h5py.File(h5_filename, "r") as f: - for path in candidates: - if path in f: - shp = f[path].shape - if len(shp) == 1: - return shp[0], 1, path - return shp[0], shp[1], path - - raise RuntimeError( - f"Could not locate data for variable '{mesh_var.clean_name}' in {h5_filename}" - ) + shp = f[dataset_path].shape + units = f[dataset_path].attrs.get("units") + if len(shp) == 1: + return shp[0], 1, units + return shp[0], shp[1], units + + def attribute_kind(var, components): + if var.vtype in (uw.VarType.TENSOR, uw.VarType.SYM_TENSOR) and components == 9: + return "Tensor" + if var.vtype in (uw.VarType.TENSOR, uw.VarType.SYM_TENSOR, uw.VarType.MATRIX): + return "Matrix" + return "Scalar" if components == 1 else "Vector" attributes = "" for var in meshVars: + if direct_p2(var) or direct_dg1(var): + continue var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" - - # Determine if data is stored on nodes (vertex_fields) or cells (cell_fields) - if not getattr(var, "continuous") or getattr(var, "degree") == 0: + projected = (var.continuous and var.degree > 1) or ( + not var.continuous and var.degree > 0 + ) + dataset_path = ( + f"visualization/{var.clean_name}" + if projected + else f"fields/{var.clean_name}" + ) + if not var.continuous or var.degree == 0: center = "Cell" else: center = "Node" - numItems, numComponents, dataset_path = get_field_info(var_filename, var, center) + numItems, numComponents, field_units = get_field_info(var_filename, dataset_path) if center == "Node" and numItems != numVertices: warnings.warn( @@ -9884,17 +10096,7 @@ def get_field_info(h5_filename, mesh_var, center): stacklevel=2, ) - # Use variable type when available, but reflect actual stored component count. - if hasattr(var, "vtype") and var.vtype in ( - uw.VarType.TENSOR, - uw.VarType.SYM_TENSOR, - uw.VarType.MATRIX, - ): - variable_type = "Tensor" - elif numComponents == 1: - variable_type = "Scalar" - else: - variable_type = "Vector" + variable_type = attribute_kind(var, numComponents) data_dimensions = f"{numItems}" if numComponents == 1 else f"{numItems} {numComponents}" var_attribute = f""" @@ -9907,7 +10109,7 @@ def get_field_info(h5_filename, mesh_var, center): Dimensions="{data_dimensions}" Format="HDF"> &{var.clean_name+"_Data"};:/{dataset_path} - + {units_information(field_units, " ")} """ attributes += var_attribute @@ -9946,8 +10148,91 @@ def get_field_info(h5_filename, mesh_var, center): """ attributes += var_attribute + special_grids = "" + for var in p2_vars: + var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" + storage_group = "fields" + with h5py.File(var_filename, "r") as f: + cells_shape = f[f"{storage_group}/cells"].shape + points_shape = f[f"{storage_group}/coordinates"].shape + values_shape = f[f"{storage_group}/{var.clean_name}"].shape + special_geometry_units = f[f"{storage_group}/coordinates"].attrs.get( + "units" + ) + field_units = f[f"{storage_group}/{var.clean_name}"].attrs.get("units") + special_topology = "Triangle_6" if var.mesh.dim == 2 else "Tetrahedron_10" + components = values_shape[1] if len(values_shape) == 2 else 1 + kind = attribute_kind(var, components) + dimensions = " ".join(str(value) for value in values_shape) + special_grids += f""" + + + + &{var.clean_name}_Data;:/{storage_group}/cells + + + + + &{var.clean_name}_Data;:/{storage_group}/coordinates + {units_information(special_geometry_units, " ")} + + + + &{var.clean_name}_Data;:/{storage_group}/{var.clean_name} + {units_information(field_units, " ")} + + """ + + if dg1_vars: + with h5py.File(mesh_filename, "r") as mesh_handle: + dg1_cells = mesh_handle["viz/dg1/cells"] + dg1_coordinates = mesh_handle["viz/dg1/coordinates"] + dg1_cells_shape = dg1_cells.shape + dg1_points_shape = dg1_coordinates.shape + dg1_topology_precision = dg1_cells.dtype.itemsize + dg1_geometry_units = dg1_coordinates.attrs.get("units") + + dg1_attributes = "" + for var in dg1_vars: + var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" + with h5py.File(var_filename, "r") as field_handle: + values = field_handle[f"fields/{var.clean_name}"] + values_shape = values.shape + field_units = values.attrs.get("units") + components = values_shape[1] if len(values_shape) == 2 else 1 + kind = attribute_kind(var, components) + dimensions = " ".join(str(value) for value in values_shape) + dg1_attributes += f""" + + + &{var.clean_name}_Data;:/fields/{var.clean_name} + {units_information(field_units, " ")} + """ + + special_grids += f""" + + + + &MeshData;:/viz/dg1/cells + + + + + &MeshData;:/viz/dg1/coordinates + {units_information(dg1_geometry_units, " ")} + +{dg1_attributes} + """ xdmf_end = f""" + {special_grids} + {collection_end} """ diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index a2622df28..e849cc16c 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/src/underworld3/discretisation/discretisation_mesh_variables.py @@ -1164,10 +1164,12 @@ def read_timestep( verbose=False, ): """ - Read a mesh variable from ``Mesh.write_timestep()`` output using the - coordinate-remap path. The saved mesh and the live mesh may have - different sizes or decompositions; values are matched to the live mesh - nodes by nearest-neighbour KDTree interpolation. + Read a mesh variable from ``Mesh.write_timestep()`` output. Standard + nodal fields use coordinate-remap interpolation. Compact simplex DG1 + fields use their exact element-corner representation: complete cells + are matched by geometry and the local affine polynomial is evaluated + at the live DG1 nodes. The latter supports MPI repartitioning without + averaging the discontinuous traces at shared vertices. This is the flexible remap reader. It is distinct from ``read_checkpoint()``, which loads PETSc DMPlex section/vector metadata @@ -1232,6 +1234,24 @@ def read_timestep( n_components = self.num_components dim = self.mesh.dim + if not is_v1_1: + if uw.mpi.rank == 0: + with h5py.File(data_file, "r") as h5f: + representation = h5f[f"fields/{data_name}"].attrs.get( + "representation" + ) + if isinstance(representation, bytes): + representation = representation.decode() + else: + representation = None + representation = uw.mpi.comm.bcast(representation, root=0) + if representation == "dg1_corner_nodal": + return self._read_dg1_corner_timestep( + data_file, + data_name, + verbose=verbose, + ) + # ---- Phase 1: source swarm carries saved (coord, value) pairs ---- source_swarm = uw.swarm.Swarm(self.mesh) saved = uw.swarm.SwarmVariable( @@ -1257,10 +1277,32 @@ def read_timestep( D_src = D_src.reshape(-1, n_components) else: with h5py.File(data_file, "r") as h5f: - X_src = h5f["fields"]["coordinates"][()].reshape(-1, dim) - D_src = h5f["fields"][data_name][()].reshape( - -1, n_components + coordinate_dataset = h5f["fields"]["coordinates"] + field_dataset = h5f["fields"][data_name] + X_src = coordinate_dataset[()].reshape(-1, dim) + D_src = field_dataset[()].reshape(-1, n_components) + storage_frame = field_dataset.attrs.get( + "storage_frame", h5f.attrs.get("storage_frame") ) + if isinstance(storage_frame, bytes): + storage_frame = storage_frame.decode() + if storage_frame == "physical": + coordinate_units = coordinate_dataset.attrs.get("units") + field_units = field_dataset.attrs.get("units") + if isinstance(coordinate_units, bytes): + coordinate_units = coordinate_units.decode() + if isinstance(field_units, bytes): + field_units = field_units.decode() + if coordinate_units and coordinate_units != "dimensionless": + X_src = np.asarray( + uw.non_dimensionalise( + uw.quantity(X_src, coordinate_units) + ) + ) + if field_units and field_units != "dimensionless": + D_src = np.asarray( + uw.non_dimensionalise(uw.quantity(D_src, field_units)) + ) else: X_src = np.empty((0, dim), dtype=np.float64) D_src = np.empty((0, n_components), dtype=np.float64) @@ -1413,6 +1455,283 @@ def read_timestep( return + def _read_dg1_corner_timestep(self, data_file, data_name, verbose=False): + """Reload one exact simplex DG1 field from disconnected corner values. + + A coordinate-only nearest-neighbour lookup cannot distinguish the + separate traces of cells meeting at one vertex. This reader routes + complete source and destination cells by centroid, matches their full + corner geometry, and evaluates each source cell's affine polynomial at + the live variable's native DG1 nodes. The cell match remains exact when + the same mesh is repartitioned across a different MPI rank count. + """ + import h5py + import numpy as np + from mpi4py import MPI + + mesh = self.mesh + dim = mesh.dim + corner_count = dim + 1 + n_components = self.num_components + if ( + self.continuous + or self.degree != 1 + or not mesh.isSimplex + or dim not in (2, 3) + or mesh.cdim != dim + ): + raise RuntimeError( + "dg1_corner_nodal data require a discontinuous degree-one " + "variable on a full-dimensional triangle or tetrahedron mesh" + ) + + if uw.mpi.rank == 0: + with h5py.File(data_file, "r") as field_handle: + field = field_handle[f"fields/{data_name}"] + values = field[()].reshape(-1, n_components) + field_units = field.attrs.get("units") + mesh_file = field_handle.attrs.get("dg1_mesh_file") + geometry_path = field_handle.attrs.get( + "dg1_geometry_path", "/viz/dg1" + ) + for name, value in ( + ("field units", field_units), + ("mesh file", mesh_file), + ("geometry path", geometry_path), + ): + if isinstance(value, bytes): + if name == "field units": + field_units = value.decode() + elif name == "mesh file": + mesh_file = value.decode() + else: + geometry_path = value.decode() + if not mesh_file: + raise RuntimeError( + f"{data_file} does not identify its shared DG1 mesh geometry" + ) + mesh_path = os.path.join(os.path.dirname(data_file), mesh_file) + with h5py.File(mesh_path, "r") as mesh_handle: + coordinates_dataset = mesh_handle[ + f"{geometry_path.strip('/')}/coordinates" + ] + coordinates = coordinates_dataset[()].reshape(-1, dim) + coordinate_units = coordinates_dataset.attrs.get("units") + cells = mesh_handle[f"{geometry_path.strip('/')}/cells"][()].reshape( + -1, corner_count + ) + if isinstance(coordinate_units, bytes): + coordinate_units = coordinate_units.decode() + if coordinate_units and coordinate_units != "dimensionless": + coordinates = np.asarray( + uw.non_dimensionalise(uw.quantity(coordinates, coordinate_units)) + ) + if field_units and field_units != "dimensionless": + values = np.asarray( + uw.non_dimensionalise(uw.quantity(values, field_units)) + ) + source_corners = coordinates[cells] + source_values = values[cells] + source_centroids = np.round(source_corners.mean(axis=1), decimals=12) + else: + source_corners = np.empty((0, corner_count, dim), dtype=np.float64) + source_values = np.empty( + (0, corner_count, n_components), dtype=np.float64 + ) + source_centroids = np.empty((0, dim), dtype=np.float64) + + source_swarm = uw.swarm.Swarm(mesh) + saved_corners = uw.swarm.SwarmVariable( + "_dg1_source_corners", + source_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * dim), + dtype=float, + _proxy=False, + varsymbol=r"\cal{C}_s", + ) + saved_values = uw.swarm.SwarmVariable( + "_dg1_source_values", + source_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * n_components), + dtype=float, + _proxy=False, + varsymbol=r"\cal{V}_s", + ) + source_before = max(source_swarm.dm.getLocalSize(), 0) + source_swarm.add_particles_with_global_coordinates( + source_centroids, migrate=False + ) + source_swarm._invalidate_canonical_data() + saved_corners.array[source_before:, 0, :] = source_corners.reshape( + -1, corner_count * dim + ) + saved_values.array[source_before:, 0, :] = source_values.reshape( + -1, corner_count * n_components + ) + source_swarm._route_by_nearest_centroid() + + landed_corners = np.asarray(saved_corners.array)[:, 0, :].reshape( + -1, corner_count, dim + ) + landed_values = np.asarray(saved_values.array)[:, 0, :].reshape( + -1, corner_count, n_components + ) + + cell_start, cell_end = mesh.dm.getHeightStratum(0) + owned = np.ones(cell_end - cell_start, dtype=bool) + if mesh.dm.comm.getSize() > 1: + _, leaves, remote = mesh.dm.getPointSF().getGraph() + if leaves is None: + leaves = np.arange(len(remote)) + leaves = np.asarray(leaves) + owned[leaves[(leaves >= cell_start) & (leaves < cell_end)] - cell_start] = ( + False + ) + native_rows = mesh._cell_node_indices(1, False).reshape(-1, corner_count)[owned] + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, corner_count)[owned] + target_nodes = np.asarray(self.coords_nd)[native_rows] + target_corners = mesh._get_coords_for_basis(1, True)[vertex_rows] + target_centroids = np.round(target_corners.mean(axis=1), decimals=12) + n_owned = len(native_rows) + + query_swarm = uw.swarm.Swarm(mesh) + origin_rank = uw.swarm.SwarmVariable( + "_dg1_origin_rank", + query_swarm, + vtype=uw.VarType.SCALAR, + dtype=int, + _proxy=False, + varsymbol=r"\cal{R}_d", + ) + origin_index = uw.swarm.SwarmVariable( + "_dg1_origin_index", + query_swarm, + vtype=uw.VarType.SCALAR, + dtype=int, + _proxy=False, + varsymbol=r"\cal{I}_d", + ) + query_nodes = uw.swarm.SwarmVariable( + "_dg1_query_nodes", + query_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * dim), + dtype=float, + _proxy=False, + varsymbol=r"\cal{X}_d", + ) + query_corners = uw.swarm.SwarmVariable( + "_dg1_query_corners", + query_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * dim), + dtype=float, + _proxy=False, + varsymbol=r"\cal{C}_d", + ) + result = uw.swarm.SwarmVariable( + "_dg1_result", + query_swarm, + vtype=uw.VarType.MATRIX, + size=(1, corner_count * n_components), + dtype=float, + _proxy=False, + varsymbol=r"\cal{D}_d", + ) + query_before = max(query_swarm.dm.getLocalSize(), 0) + query_swarm.add_particles_with_global_coordinates( + target_centroids, migrate=False + ) + query_swarm._invalidate_canonical_data() + origin_rank.array[query_before:, 0, 0] = uw.mpi.rank + origin_index.array[query_before:, 0, 0] = np.arange(n_owned) + query_nodes.array[query_before:, 0, :] = target_nodes.reshape( + -1, corner_count * dim + ) + query_corners.array[query_before:, 0, :] = target_corners.reshape( + -1, corner_count * dim + ) + query_swarm._route_by_nearest_centroid() + + def cell_key(corners): + rounded = np.round(np.asarray(corners), decimals=12) + order = np.lexsort( + tuple(rounded[:, axis] for axis in range(dim - 1, -1, -1)) + ) + return tuple(rounded[order].reshape(-1)) + + source_by_cell = {} + duplicate_source = 0 + for source_index, corners in enumerate(landed_corners): + key = cell_key(corners) + if key in source_by_cell: + duplicate_source += 1 + source_by_cell[key] = source_index + + local_query_nodes = np.asarray(query_nodes.array)[:, 0, :].reshape( + -1, corner_count, dim + ) + local_query_corners = np.asarray(query_corners.array)[:, 0, :].reshape( + -1, corner_count, dim + ) + missing = 0 + for query_index, (query_cell_corners, nodes) in enumerate( + zip(local_query_corners, local_query_nodes, strict=True) + ): + key = cell_key(query_cell_corners) + matched_index = source_by_cell.get(key) + if matched_index is None: + missing += 1 + continue + corners = landed_corners[matched_index] + coefficients = landed_values[matched_index] + matrix = (corners[1:] - corners[:1]).T + local = np.linalg.solve(matrix, (nodes - corners[:1]).T).T + weights = np.column_stack((1.0 - local.sum(axis=1), local)) + result.array[query_index, 0, :] = (weights @ coefficients).reshape(-1) + + failures = uw.mpi.comm.allreduce(missing + duplicate_source, op=MPI.SUM) + if failures: + raise RuntimeError( + "DG1 corner reload could not match every destination cell to " + "one source cell. The cell-aware path supports identical mesh " + "geometry across arbitrary MPI repartitioning." + ) + + query_swarm._rank_var.array[...] = origin_rank.array[...] + query_swarm.dm.migrate(remove_sent_points=True) + uw.mpi.barrier() + query_swarm._invalidate_canonical_data() + + returned_indices = origin_index.array[:, 0, 0].astype(int) + returned = np.asarray(result.array)[:, 0, :].reshape( + -1, corner_count, n_components + ) + cell_values = np.empty((n_owned, corner_count, n_components)) + cell_values[returned_indices] = returned + self.data[...] = 0.0 + self.data[native_rows.reshape(-1), :] = cell_values.reshape( + -1, n_components + ) + + indexset, subdm = mesh.dm.createSubDM(self.field_id) + subdm.localToGlobal(self._lvec, self._gvec, addv=False) + subdm.globalToLocal(self._gvec, self._lvec, addv=False) + indexset.destroy() + subdm.destroy() + if mesh._lvec is not None: + mesh._lvec.destroy() + mesh._lvec = None + mesh._stale_lvec = True + if verbose and uw.mpi.rank == 0: + print( + f"Reloaded DG1 field {data_name!r} from exact corner-nodal data", + flush=True, + ) + return + @timing.routine_timer_decorator def load_from_h5_plex_vector( self, @@ -1453,11 +1772,12 @@ def read_checkpoint( ): """Load this mesh variable from PETSc reload output. - The default path restores DMPlex section/local-vector data through the + By default, DMPlex section/global-vector data are restored through the topology migration SF, so a mesh reconstructed from its checkpoint may - have a different parallel DOF ordering. Set ``same_layout=True`` only - for an in-place restore onto the exact mesh object that wrote the file; - that path reloads the saved global vector directly. + have a different parallel DOF ordering. With ``same_layout=True``, the + existing PETSc variable vector is loaded directly into the original + layout. Both paths read the same stored values; no duplicate checkpoint + vector is required. This method does not use the coordinate/KDTree remapping provided by ``read_timestep()``. New output should be written with @@ -1472,31 +1792,43 @@ def read_checkpoint( if self._lvec is None: self._set_vec(available=True) - if same_layout: - import h5py + import h5py - if uw.mpi.rank == 0: - with h5py.File(filename, "r") as checkpoint_h5: - has_direct_vector = ( - "uw_checkpoint" in checkpoint_h5 - and data_name in checkpoint_h5["uw_checkpoint"] - ) - else: - has_direct_vector = None - has_direct_vector = uw.mpi.comm.bcast( - has_direct_vector, - root=0, + if uw.mpi.rank == 0: + with h5py.File(filename, "r") as checkpoint_h5: + if "restart/petsc/topologies" in checkpoint_h5: + checkpoint_group = "/restart/petsc" + elif "uw_checkpoint/topologies" in checkpoint_h5: + checkpoint_group = "/uw_checkpoint" + else: + checkpoint_group = None + + if f"restart/petsc/{data_name}" in checkpoint_h5: + direct_vector_group = "/restart/petsc" + elif f"uw_checkpoint/{data_name}" in checkpoint_h5: + direct_vector_group = "/uw_checkpoint" + else: + direct_vector_group = None + else: + checkpoint_group = None + direct_vector_group = None + checkpoint_group = uw.mpi.comm.bcast(checkpoint_group, root=0) + direct_vector_group = uw.mpi.comm.bcast(direct_vector_group, root=0) + grouped_checkpoint = checkpoint_group is not None + legacy_direct_vector = direct_vector_group is not None + + if same_layout and not (grouped_checkpoint or legacy_direct_vector): + raise RuntimeError( + f"{filename} has no direct checkpoint vector for {data_name!r}. " + "Reload it with same_layout=False." ) - if not has_direct_vector: - raise RuntimeError( - f"{filename} has no in-place checkpoint vector for " - f"{data_name!r}. Reload it with same_layout=False." - ) indexset, subdm = self.mesh.dm.createSubDM(self.field_id) sectiondm = self.mesh.dm.clone() viewer = PETSc.ViewerHDF5().create(filename, "r", comm=PETSc.COMM_WORLD) viewer.pushFormat(PETSc.Viewer.Format.HDF5_PETSC) + if grouped_checkpoint and not same_layout: + viewer.pushGroup(checkpoint_group) old_mesh_name = self.mesh.dm.getName() old_lvec_name = self._lvec.getName() @@ -1510,25 +1842,37 @@ def read_checkpoint( self._gvec.setName(data_name) if same_layout: + vector_group = ( + f"{checkpoint_group}/topologies/uw_mesh/dms/{data_name}/" + f"vecs/{data_name}" + if grouped_checkpoint + else direct_vector_group + ) + # A DM-associated Vec ignores the viewer group and redirects + # HDF5 I/O through /fields. A plain Vec reads the requested + # restart dataset when /fields is intentionally absent. checkpoint_vec = PETSc.Vec().createMPI( - (self._gvec.getLocalSize(), self._gvec.getSize()), + (self._gvec.getLocalSize(), PETSc.DECIDE), comm=PETSc.COMM_WORLD, ) checkpoint_vec.setName(data_name) - viewer.pushGroup("/uw_checkpoint") + viewer.pushGroup(vector_group) checkpoint_vec.load(viewer) viewer.popGroup() self._gvec.array[...] = checkpoint_vec.array_r - checkpoint_vec.destroy() subdm.globalToLocal(self._gvec, self._lvec, addv=False) + checkpoint_vec.destroy() else: - from underworld3.cython.petsc_discretisation import ( - petsc_dmplex_load_local_vector, + global_sf, local_sf = self.mesh.dm.sectionLoad( + viewer, sectiondm, self.mesh.sf ) - - loaded_lvec = petsc_dmplex_load_local_vector( - self.mesh.dm, viewer, sectiondm, self.mesh.sf, data_name + loaded_gvec = sectiondm.createGlobalVec() + loaded_gvec.setName(data_name) + self.mesh.dm.globalVectorLoad( + viewer, sectiondm, global_sf, loaded_gvec ) + loaded_lvec = sectiondm.createLocalVec() + sectiondm.globalToLocal(loaded_gvec, loaded_lvec, addv=False) source_section = sectiondm.getSection() target_section = subdm.getSection() @@ -1555,12 +1899,19 @@ def read_checkpoint( ) loaded_lvec.destroy() + loaded_gvec.destroy() + if global_sf is not None: + global_sf.destroy() + if local_sf is not None: + local_sf.destroy() self._sync_lvec_to_gvec() finally: self._lvec.setName(old_lvec_name) self._gvec.setName(old_vec_name) if old_mesh_name is not None: self.mesh.dm.setName(old_mesh_name) + if grouped_checkpoint and not same_layout: + viewer.popGroup() viewer.popFormat() viewer.destroy() sectiondm.destroy() diff --git a/src/underworld3/function/field_projection.py b/src/underworld3/function/field_projection.py index 28fc4809d..4ed0c713f 100644 --- a/src/underworld3/function/field_projection.py +++ b/src/underworld3/function/field_projection.py @@ -226,6 +226,59 @@ def _write_vec_to_group(viewer, data_array, name, group, comm): vec.destroy() +def _write_index_array_to_group(viewer, data_array, name, group, comm): + """Write a distributed integer connectivity array to an HDF5 group.""" + indices = PETSc.IS().createGeneral( + np.asarray(data_array, dtype=PETSc.IntType).reshape(-1), comm=comm + ) + if data_array.ndim == 2: + indices.setBlockSize(data_array.shape[1]) + indices.setName(name) + viewer.pushGroup(group) + viewer(indices) + viewer.popGroup() + indices.destroy() + + +def _physical_visualisation_enabled(units): + """Return whether declared units require a physical XDMF copy.""" + import underworld3 as uw + + return ( + units is not None + and uw.get_default_model().has_units_active() + and uw.is_nondimensional_scaling_active() + ) + + +def _physical_visualisation_values(data, units): + """Return output values and their declared unit label. + + Solver and checkpoint vectors use model magnitudes. XDMF arrays are a + user-facing boundary, so an active nondimensional model is converted to + the units declared by the mesh or variable before those arrays are written. + """ + import underworld3 as uw + + if units is None: + return data, None + + target_units = uw.units(units).units if isinstance(units, str) else units + unit_label = str(target_units) + if not _physical_visualisation_enabled(units): + return data, unit_label + + dimensionality = dict(target_units.dimensionality) + if not dimensionality: + return data, unit_label + + physical = uw.dimensionalise( + np.asarray(data), + target_dimensionality=dimensionality, + ).to(target_units) + return np.asarray(physical), unit_label + + def write_vertices_to_viewer( mesh_var: "MeshVariable", viewer: "PETSc.ViewerHDF5", @@ -293,9 +346,140 @@ def write_vertices_to_viewer( if is_tensor: data = _repack_tensor_to_paraview(data, mesh_var.vtype, mesh.dim) + data, _ = _physical_visualisation_values(data, mesh_var.units) + + _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) + + +def write_field_to_viewer( + mesh_var: "MeshVariable", + viewer: "PETSc.ViewerHDF5", + group: str, + name: str, +) -> None: + """Write a variable's owned values in physical units without projection.""" + mesh_var._sync_lvec_to_gvec() + data = mesh_var._gvec.array.reshape(-1, mesh_var.num_components).copy() + data, _ = _physical_visualisation_values(data, mesh_var.units) + _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) + + +def write_field_coordinates_to_viewer( + mesh_var: "MeshVariable", + viewer: "PETSc.ViewerHDF5", + group: str, + name: str = "coordinates", +) -> None: + """Write owned coordinates for a variable's exact finite-element layout.""" + mesh = mesh_var.mesh + coordinate_dm = mesh._basis_coordinate_dm(mesh_var.degree, mesh_var.continuous) + local = coordinate_dm.getLocalVec() + global_vector = coordinate_dm.getGlobalVec() + local.array[...] = np.asarray(mesh_var.coords_nd).reshape(-1) + coordinate_dm.localToGlobal(local, global_vector, addv=False) + coordinates = global_vector.array.reshape(-1, mesh.cdim).copy() + coordinates, _ = _physical_visualisation_values(coordinates, mesh.units) + _write_vec_to_group(viewer, coordinates, name, group, PETSc.COMM_WORLD) + coordinate_dm.restoreGlobalVec(global_vector) + coordinate_dm.restoreLocalVec(local) + coordinate_dm.destroy() + + +def write_projected_field_to_viewer( + mesh_var: "MeshVariable", + viewer: "PETSc.ViewerHDF5", + target_degree: int, + continuous: bool, + group: str, + name: str, +) -> None: + """Write a compact visualization projection for an unsupported layout.""" + data = project_to_degree( + mesh_var, + target_degree=target_degree, + continuous=continuous, + include_ghosts=False, + ) + data, _ = _physical_visualisation_values(data, mesh_var.units) _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) +def write_p2_simplex_topology_to_viewer(mesh_var, viewer, group="/fields"): + """Write VTK-ordered quadratic triangle or tetrahedron connectivity.""" + mesh = mesh_var.mesh + if not ( + mesh_var.continuous + and mesh_var.degree == 2 + and mesh.isSimplex + and mesh.dim in (2, 3) + and mesh.cdim == mesh.dim + ): + raise NotImplementedError( + "direct P2 XDMF requires a full-dimensional triangle or tetrahedron mesh" + ) + + coordinate_dm = mesh._basis_coordinate_dm(2, True) + local_section = coordinate_dm.getLocalSection() + global_section = coordinate_dm.getGlobalSection() + local_to_global = np.full(len(mesh_var.coords_nd), -1, dtype=PETSc.IntType) + point_start, point_end = local_section.getChart() + for point in range(point_start, point_end): + node_count = local_section.getDof(point) // mesh.cdim + if node_count == 0: + continue + local_offset = local_section.getOffset(point) // mesh.cdim + global_offset = global_section.getOffset(point) + if global_offset < 0: + global_offset = -(global_offset + 1) + global_offset //= mesh.cdim + local_to_global[local_offset : local_offset + node_count] = np.arange( + global_offset, global_offset + node_count, dtype=PETSc.IntType + ) + coordinate_dm.destroy() + + cell_start, cell_end = mesh.dm.getHeightStratum(0) + owned = np.ones(cell_end - cell_start, dtype=bool) + if mesh.dm.comm.getSize() > 1: + _, leaves, remote = mesh.dm.getPointSF().getGraph() + if leaves is None: + leaves = np.arange(len(remote)) + leaves = np.asarray(leaves) + owned[leaves[(leaves >= cell_start) & (leaves < cell_end)] - cell_start] = False + + corner_count = mesh.dim + 1 + node_count = 6 if mesh.dim == 2 else 10 + p2_rows = mesh._cell_node_indices(2, True).reshape(-1, node_count)[owned] + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, corner_count)[owned] + p2_coordinates = mesh_var.coords_nd[p2_rows] + corners = mesh._get_coords_for_basis(1, True)[vertex_rows] + negative = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) < 0 + if mesh.dim == 2: + corners[negative] = corners[negative][:, [0, 2, 1]] + edge_pairs = ((0, 1), (1, 2), (2, 0)) + topology_name = "Triangle_6" + else: + corners[negative] = corners[negative][:, [0, 2, 1, 3]] + edge_pairs = ((0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)) + topology_name = "Tetrahedron_10" + + connectivity = np.empty_like(p2_rows, dtype=PETSc.IntType) + for cell_index, (row, nodes, vertices) in enumerate( + zip(p2_rows, p2_coordinates, corners, strict=True) + ): + targets = np.vstack( + (vertices, *(0.5 * (vertices[a] + vertices[b]) for a, b in edge_pairs)) + ) + order = [np.argmin(np.linalg.norm(nodes - target, axis=1)) for target in targets] + if len(set(order)) != node_count: + raise RuntimeError(f"could not map UW3 P2 nodes to {topology_name} ordering") + connectivity[cell_index] = local_to_global[row[order]] + if np.any(connectivity < 0): + raise RuntimeError("P2 XDMF connectivity contains an unmapped global node") + _write_index_array_to_group( + viewer, connectivity, "cells", group, PETSc.COMM_WORLD + ) + + def write_coordinates_to_viewer( mesh, viewer: "PETSc.ViewerHDF5", @@ -316,7 +500,8 @@ def write_coordinates_to_viewer( Dataset name (default ``coordinates``). """ coord_gvec = mesh.dm.getCoordinates() - coords = coord_gvec.array.reshape(-1, mesh.dim).copy() + coords = coord_gvec.array.reshape(-1, mesh.cdim).copy() + coords, _ = _physical_visualisation_values(coords, mesh.units) _write_vec_to_group(viewer, coords, name, group, PETSc.COMM_WORLD) @@ -349,4 +534,75 @@ def write_cell_field_to_viewer( nc = mesh_var.num_components mesh_var._sync_lvec_to_gvec() data = mesh_var._gvec.array.reshape(-1, nc).copy() + data, _ = _physical_visualisation_values(data, mesh_var.units) _write_vec_to_group(viewer, data, name, group, PETSc.COMM_WORLD) + + +def _dg1_corner_data(mesh_var, repack_tensors=False): + """Return owned disconnected-cell geometry and exact DG1 corner values. + + Coordinate-section cell maps preserve element ownership and node ordering; + no point location, coordinate matching, or inter-element averaging is used. + Interior DG interpolation nodes define an affine polynomial, evaluated at + that same cell's vertices. The corner values are an exact alternative basis + for the element-local linear polynomial, not a lower-order projection. + """ + mesh = mesh_var.mesh + if ( + mesh_var.continuous or mesh_var.degree != 1 or not mesh.isSimplex + or mesh.dim not in (2, 3) or mesh.cdim != mesh.dim + ): + raise NotImplementedError("DG1 XDMF requires a full-dimensional triangle/tetrahedron mesh") + cstart, cend = mesh.dm.getHeightStratum(0) + owned = np.ones(cend - cstart, dtype=bool) + # Serial DMPlex may have an unset point SF; no cells are ghosts there. + if mesh.dm.comm.getSize() > 1: + _, leaves, remote = mesh.dm.getPointSF().getGraph() + if leaves is None: + leaves = np.arange(len(remote)) + leaves = np.asarray(leaves) + owned[leaves[(leaves >= cstart) & (leaves < cend)] - cstart] = False + rows = mesh._cell_node_indices(1, False).reshape(-1, mesh.dim + 1)[owned] + vertex_rows = mesh._cell_node_indices(1, True).reshape(-1, mesh.dim + 1)[owned] + corners = mesh._get_coords_for_basis(1, True)[vertex_rows] + # Closure order is arbitrary; give VTK positively oriented simplices. + negative = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) < 0 + corners[negative] = corners[negative][:, [0, 2, 1] if mesh.dim == 2 else [0, 2, 1, 3]] + nodes = mesh_var.coords_nd[rows] + coefficients = mesh_var._lvec.array.reshape(-1, mesh_var.num_components)[rows] + matrix = (nodes[:, 1:] - nodes[:, :1]).transpose(0, 2, 1) + local = np.linalg.solve(matrix, (corners - nodes[:, :1]).transpose(0, 2, 1)) + weights = np.concatenate((1 - local.sum(axis=1, keepdims=True), local), axis=1) + values = np.einsum("cij,cik->cjk", weights, coefficients).reshape(-1, mesh_var.num_components) + if repack_tensors: + values = _repack_tensor_to_paraview(values, mesh_var.vtype, mesh.dim) + corners, _ = _physical_visualisation_values(corners, mesh.units) + values, _ = _physical_visualisation_values(values, mesh_var.units) + corner_rows = corners.reshape(-1, mesh.cdim) + local_count = len(corner_rows) + offset = mesh.dm.comm.tompi4py().exscan(local_count) + if offset is None: + offset = 0 + cells = np.arange(offset, offset + local_count, dtype=PETSc.IntType).reshape( + -1, mesh.dim + 1 + ) + return corner_rows, values, cells + + +def _write_dg1_to_viewer( + mesh_var, + viewer, + group="/dg1", + coordinate_name="vertices", + value_name="values", + repack_tensors=True, +): + """Write owned simplex cells with independent vertices and DG1 traces.""" + corner_rows, values, cells = _dg1_corner_data( + mesh_var, repack_tensors=repack_tensors + ) + _write_vec_to_group( + viewer, corner_rows, coordinate_name, group, PETSc.COMM_WORLD + ) + _write_vec_to_group(viewer, values, value_name, group, PETSc.COMM_WORLD) + _write_index_array_to_group(viewer, cells, "cells", group, PETSc.COMM_WORLD) diff --git a/src/underworld3/meshing/surfaces.py b/src/underworld3/meshing/surfaces.py index d62027e28..030e65ead 100644 --- a/src/underworld3/meshing/surfaces.py +++ b/src/underworld3/meshing/surfaces.py @@ -1318,6 +1318,7 @@ def _compute_distance_field(self) -> None: 1, degree=self.mesh.degree, varsymbol=f"d_{{{self._symbol}}}", + units=self.mesh.units, ) # Get mesh coordinates in model (internal) space. @@ -1351,6 +1352,7 @@ def _compute_distance_field(self) -> None: 1, degree=self.mesh.degree, varsymbol=f"|d_{{{self._symbol}}}|", + units=self.mesh.units, ) with uw.synchronised_array_update(): diff --git a/tests/test_0003_save_load.py b/tests/test_0003_save_load.py index 0b2d1839b..f2f2c98fd 100644 --- a/tests/test_0003_save_load.py +++ b/tests/test_0003_save_load.py @@ -172,9 +172,9 @@ def test_timestep_with_petsc_reload_roundtrip(tmp_path): with h5py.File(var_file, "r") as h5f: assert "fields/x" in h5f assert "fields/coordinates" in h5f - assert "vertex_fields/x_x" in h5f - assert "topologies/uw_mesh/dms/x/section" in h5f - assert "topologies/uw_mesh/dms/x/vecs/x" in h5f + assert "vertex_fields" not in h5f + assert "restart/petsc/topologies/uw_mesh/dms/x/section" in h5f + assert "restart/petsc/topologies/uw_mesh/dms/x/vecs/x" in h5f uw.mpi.barrier() mesh_reloaded = uw.discretisation.Mesh(f"{tmp_path}/unified.mesh.00000.h5") @@ -205,8 +205,8 @@ def test_timestep_with_petsc_reload_roundtrip(tmp_path): assert checkpoint_xdmf_file.is_file() with h5py.File(checkpoint_var_file, "r") as h5f: assert "fields/x" in h5f - assert "vertex_fields/x_x" in h5f - assert "topologies/uw_mesh/dms/x/section" in h5f + assert "vertex_fields" not in h5f + assert "restart/petsc/topologies/uw_mesh/dms/x/section" in h5f uw.mpi.barrier() mesh_reloaded = uw.discretisation.Mesh( diff --git a/tests/test_0005_xdmf_compat.py b/tests/test_0005_xdmf_compat.py index 3c7b34d5c..fae982db6 100644 --- a/tests/test_0005_xdmf_compat.py +++ b/tests/test_0005_xdmf_compat.py @@ -1,11 +1,8 @@ -"""Test XDMF/HDF5 compatibility groups written by write_timestep. +"""Compact XDMF/HDF5 layouts for continuous and cell fields.""" -Validates that /vertex_fields/ and /cell_fields/ groups are created -correctly, with proper component counts and data values. -""" - -import os +from pathlib import Path import re +import xml.etree.ElementTree as ET import h5py import numpy as np @@ -14,425 +11,324 @@ import underworld3 as uw -# --------------------------------------------------------------------------- -# Helpers -# --------------------------------------------------------------------------- - - -def _check_h5_group_exists(h5_path, group_path): - """Return True if group_path exists in the HDF5 file.""" - with h5py.File(h5_path, "r") as f: - return group_path in f - - -def _read_h5_dataset(h5_path, dataset_path): - """Read and return an HDF5 dataset as a numpy array.""" - with h5py.File(h5_path, "r") as f: - return f[dataset_path][:] - - -def _check_xdmf_refs(xdmf_path, tmp_dir): - """Verify all XDMF entity references point to real HDF5 datasets.""" - with open(xdmf_path, "r") as f: - content = f.read() - - doctype_match = re.search(r"", content, re.DOTALL) - assert doctype_match, "No DOCTYPE entity block found in XDMF file" - - entity_block = doctype_match.group(1) - entities = dict(re.findall(r'', entity_block)) - - # Check both vertex_fields and cell_fields references - refs = re.findall(r"&(\w+);:(/(vertex_fields|cell_fields)/[A-Za-z0-9_]+)", content) - errors = [] - for entity_name, dataset_path, _ in refs: - h5_file = entities.get(entity_name) - if not h5_file: - errors.append(f"Entity {entity_name} not found") - continue - h5_full = os.path.join(tmp_dir, h5_file) - if not os.path.exists(h5_full): - errors.append(f"File {h5_file} not found") - continue - with h5py.File(h5_full, "r") as f: - if dataset_path.lstrip("/") not in f: - errors.append(f"{h5_file}: {dataset_path} missing") - - assert not errors, "XDMF reference errors:\n" + "\n".join(errors) - - -# --------------------------------------------------------------------------- -# Test: P1 scalar + P2 vector (2D) -# --------------------------------------------------------------------------- - - -def test_xdmf_compat_2d(tmp_path): - """write_timestep creates correct compat groups for 2D mesh variables.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - - # P1 scalar, P2 vector - p_var = uw.discretisation.MeshVariable("p", mesh, 1, degree=1) - v_var = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2) - - # Initialise with known values - x, y = mesh.X - p_var.data[:, 0] = mesh._coords[:, 0] # p = x coordinate - v_var.data[:, 0] = 1.0 - v_var.data[:, 1] = 2.0 - - mesh.write_timestep( - "test", index=0, outputPath=str(tmp_path), meshVars=[p_var, v_var] - ) - - # Check XDMF file was created - xdmf_file = os.path.join(str(tmp_path), "test.mesh.00000.xdmf") - assert os.path.exists(xdmf_file), "XDMF file not created" - - # Check P1 scalar: /vertex_fields/p_p should exist with correct shape - p_h5 = os.path.join(str(tmp_path), "test.mesh.p.00000.h5") - assert _check_h5_group_exists(p_h5, "vertex_fields/p_p") - assert _check_h5_group_exists(p_h5, "vertex_fields/coordinates") - - p_compat = _read_h5_dataset(p_h5, "vertex_fields/p_p") - # PETSc writes 1-component scalars as 1D (N,); accept both (N,) and (N,1) - effective_ncomp = p_compat.shape[1] if p_compat.ndim == 2 else 1 - assert effective_ncomp == 1, f"P1 scalar should have 1 component, got {effective_ncomp}" - - # P1 scalar: compat values should match var.data exactly - p_original = _read_h5_dataset(p_h5, "fields/p") - np.testing.assert_allclose( - p_compat.ravel(), p_original.ravel(), atol=1e-10, - err_msg="P1 compat data should match original field data exactly" - ) - - # Check P2 vector: /vertex_fields/v_v should exist with dim components - v_h5 = os.path.join(str(tmp_path), "test.mesh.v.00000.h5") - assert _check_h5_group_exists(v_h5, "vertex_fields/v_v") - - v_compat = _read_h5_dataset(v_h5, "vertex_fields/v_v") - # Standalone Vec writes as 1D — infer components from total size / vertex count - n_verts = mesh._coords.shape[0] - v_total = v_compat.size - v_ncomp = v_total // n_verts if n_verts > 0 else 0 - assert v_ncomp == mesh.dim, ( - f"P2 vector should have {mesh.dim} components, got {v_ncomp}" - ) - - # P2 vector: vertex count should match mesh vertices, not P2 DOFs - coords_compat = _read_h5_dataset(v_h5, "vertex_fields/coordinates") - coords_n_verts = coords_compat.size // mesh.dim - assert coords_n_verts == mesh._coords.shape[0], ( - "Vertex count mismatch between compat coords and mesh" - ) - - # Verify XDMF references are valid - _check_xdmf_refs(xdmf_file, str(tmp_path)) - - del mesh +_OLD_GROUPS = ("vertex_fields", "cell_fields", "dg1") -def test_write_timestep_mesh_keeps_restart_and_viz_topology(tmp_path): - """write_timestep mesh output keeps restart data plus XDMF topology.""" +def _shared_path(tmp_path): + return Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - s_var = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s_var.data[:, 0] = mesh._coords[:, 0] - mesh.write_timestep( - "viztopo", index=0, outputPath=str(tmp_path), meshVars=[s_var] - ) - - mesh_h5 = os.path.join(str(tmp_path), "viztopo.mesh.00000.h5") - xdmf_file = os.path.join(str(tmp_path), "viztopo.mesh.00000.xdmf") - - with h5py.File(mesh_h5, "r") as h5f: - assert "labels" in h5f - assert "topology/cells" in h5f - assert "topology/cones" in h5f - assert "viz/topology/cells" in h5f - cells = h5f["viz/topology/cells"] - assert len(cells.shape) == 2 - assert cells.shape[1] > 1 - - with open(xdmf_file, "r") as f: - xdmf_text = f.read() - assert "&MeshData;:/viz/topology/cells" in xdmf_text - assert "&MeshData;:/topology/cells" not in xdmf_text - assert 'ItemType="HyperSlab"' not in xdmf_text - - del mesh - - -def test_write_checkpoint_mesh_uses_petsc_topology(tmp_path): - """write_checkpoint keeps PETSc DMPlex topology for restart output.""" +def _assert_xdmf_references_exist(xdmf_path): + """Check every HDF reference and advertised shape in an XDMF file.""" + text = xdmf_path.read_text() + entities = dict(re.findall(r'', text)) + tree = ET.parse(xdmf_path) + for item in tree.findall(".//DataItem[@Format='HDF']"): + reference, dataset = item.text.strip().split(":", 1) + entity = reference.removeprefix("&").removesuffix(";") + h5_name = entities.get(entity, reference) + with h5py.File(xdmf_path.parent / h5_name, "r") as handle: + assert dataset in handle + dimensions = item.get("Dimensions") + if dimensions: + assert tuple(map(int, dimensions.split())) == handle[dataset].shape - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - s_var = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s_var.data[:, 0] = mesh._coords[:, 0] - with pytest.warns(FutureWarning, match="write_checkpoint\\(\\) is deprecated"): - mesh.write_checkpoint( - "restart", - outputPath=str(tmp_path), - meshUpdates=True, - meshVars=[s_var], - index=0, - ) +def _assert_no_compatibility_copies(handle): + for group in _OLD_GROUPS: + assert group not in handle - mesh_h5 = os.path.join(str(tmp_path), "restart.mesh.00000.h5") - with h5py.File(mesh_h5, "r") as h5f: - assert "topologies/uw_mesh/topology" in h5f - assert "viz/topology/cells" not in h5f - - del mesh - - -# --------------------------------------------------------------------------- -# Test: 3D mesh -# --------------------------------------------------------------------------- - - -def test_xdmf_compat_3d(tmp_path): - """write_timestep creates correct compat groups for 3D mesh.""" - - mesh = uw.meshing.StructuredQuadBox( - elementRes=(3, 3, 3), +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_direct_p1_p2_and_dg0_use_fields_only(tmp_path): + """P1, triangular P2, and DG0 are visualized directly from /fields.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, regular=True, qdegree=3) + p1 = uw.discretisation.MeshVariable("p1", mesh, 1, degree=1) + p2 = uw.discretisation.MeshVariable("p2", mesh, 2, degree=2) + dg0 = uw.discretisation.MeshVariable("dg0", mesh, 1, degree=0, continuous=False) + + p1.array[:, 0, 0] = p1.coords[:, 0] + 2.0 * p1.coords[:, 1] + p2.array[:, 0, 0] = p2.coords[:, 0] + p2.array[:, 0, 1] = p2.coords[:, 1] + dg0.array[:, 0, 0] = 7.0 + mesh.write_timestep("compact", 0, outputPath=str(directory), meshVars=[p1, p2, dg0]) + + if uw.mpi.rank == 0: + p1_file = directory / "compact.mesh.p1.00000.h5" + p2_file = directory / "compact.mesh.p2.00000.h5" + dg0_file = directory / "compact.mesh.dg0.00000.h5" + with h5py.File(p1_file, "r") as handle: + assert set(handle) == {"fields"} + assert set(handle["fields"]) == {"coordinates", "p1"} + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + with h5py.File(dg0_file, "r") as handle: + assert set(handle) == {"fields"} + assert handle["fields/dg0"].shape[0] == handle["fields/coordinates"].shape[0] + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + with h5py.File(p2_file, "r") as handle: + assert set(handle) == {"fields"} + coordinates = handle["fields/coordinates"][:] + cells = handle["fields/cells"][:] + values = handle["fields/p2"][:] + assert cells.shape[1] == 6 + assert values.shape == coordinates.shape + np.testing.assert_allclose(values, coordinates, atol=1.0e-12) + points = coordinates[cells] + np.testing.assert_allclose(points[:, 3], 0.5 * (points[:, 0] + points[:, 1])) + np.testing.assert_allclose(points[:, 4], 0.5 * (points[:, 1] + points[:, 2])) + np.testing.assert_allclose(points[:, 5], 0.5 * (points[:, 2] + points[:, 0])) + _assert_no_compatibility_copies(handle) + + xdmf = directory / "compact.mesh.00000.xdmf" + text = xdmf.read_text() + # ParaView's XDMF3 reader does not accept XML-reference DataItems + # without a Format attribute. Inline the two mesh HDF references so + # both its XDMF3 reader and the legacy XDMF reader can open the file. + assert 'Reference="XML"' not in text + assert 'TopologyType="Triangle_6"' in text + assert "&p1_Data;:/fields/p1" in text + assert "&p2_Data;:/fields/p2" in text + assert "&dg0_Data;:/fields/dg0" in text + _assert_xdmf_references_exist(xdmf) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_direct_p2_tetrahedron_uses_exact_tetrahedron_10(tmp_path): + """A 3D continuous P2 field uses exact VTK-ordered tetrahedral data.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox( minCoords=(0.0, 0.0, 0.0), maxCoords=(1.0, 1.0, 1.0), + cellSize=0.7, + regular=True, + qdegree=3, ) - - s_var = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s_var.data[:, 0] = mesh._coords[:, 2] # s = z - - mesh.write_timestep( - "test3d", index=0, outputPath=str(tmp_path), meshVars=[s_var] - ) - - s_h5 = os.path.join(str(tmp_path), "test3d.mesh.s.00000.h5") - assert _check_h5_group_exists(s_h5, "vertex_fields/s_s") - - s_compat = _read_h5_dataset(s_h5, "vertex_fields/s_s") - s_original = _read_h5_dataset(s_h5, "fields/s") - np.testing.assert_allclose(s_compat.ravel(), s_original.ravel(), atol=1e-10) - - # Verify XDMF - xdmf_file = os.path.join(str(tmp_path), "test3d.mesh.00000.xdmf") - _check_xdmf_refs(xdmf_file, str(tmp_path)) - - mesh_h5 = os.path.join(str(tmp_path), "test3d.mesh.00000.h5") - with h5py.File(mesh_h5, "r") as h5f: - assert "viz/topology/cells" in h5f - cells = h5f["viz/topology/cells"] - assert len(cells.shape) == 2 - assert cells.shape[1] > 1 - - with open(xdmf_file, "r") as f: - xdmf_text = f.read() - assert "&MeshData;:/viz/topology/cells" in xdmf_text - assert "&MeshData;:/topology/cells" not in xdmf_text - assert f'Dimensions="{s_compat.shape[0]}"' in xdmf_text - - del mesh - - -# --------------------------------------------------------------------------- -# Test: Cell (discontinuous) variable -# --------------------------------------------------------------------------- - - -def test_xdmf_compat_cell_variable(tmp_path): - """Cell variables go to /cell_fields/ group.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - - # Discontinuous (cell-centred) variable - c_var = uw.discretisation.MeshVariable( - "c", mesh, 1, degree=0, continuous=False - ) - c_var.data[:, 0] = 42.0 - - mesh.write_timestep( - "testcell", index=0, outputPath=str(tmp_path), meshVars=[c_var] - ) - - c_h5 = os.path.join(str(tmp_path), "testcell.mesh.c.00000.h5") - assert _check_h5_group_exists(c_h5, "cell_fields/c_c"), ( - "/cell_fields/c_c group should exist for discontinuous variable" - ) - - c_compat = _read_h5_dataset(c_h5, "cell_fields/c_c") - c_effective_ncomp = c_compat.shape[1] if c_compat.ndim == 2 else 1 - assert c_effective_ncomp == 1 - - del mesh - - -# --------------------------------------------------------------------------- -# Test: read_timestep round-trip still works (checkpoint integrity) -# --------------------------------------------------------------------------- - - -def test_xdmf_checkpoint_roundtrip(tmp_path): - """Compat groups don't corrupt the /fields/ data used by read_timestep.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - - p1 = uw.discretisation.MeshVariable("p1", mesh, 1, degree=1) - p1.data[:, 0] = mesh._coords[:, 0] + mesh._coords[:, 1] - - mesh.write_timestep( - "roundtrip", index=0, outputPath=str(tmp_path), meshVars=[p1] - ) - - # Read back - p1_check = uw.discretisation.MeshVariable("p1check", mesh, 1, degree=1) - p1_check.read_timestep("roundtrip", "p1", 0, outputPath=str(tmp_path)) - - np.testing.assert_allclose(p1.data, p1_check.data, atol=1e-10) - - del mesh - - -# --------------------------------------------------------------------------- -# Test: create_xdmf=False skips compat groups -# --------------------------------------------------------------------------- - - -def test_no_xdmf_when_disabled(tmp_path): - """create_xdmf=False should not create compat groups or XDMF file.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - s = uw.discretisation.MeshVariable("s", mesh, 1, degree=1) - s.data[:, 0] = 1.0 # initialise so gvec exists - + p2 = uw.discretisation.MeshVariable("p2_tet", mesh, 3, degree=2) + p2.array[:, 0, :] = p2.coords + mesh.write_timestep("tet", 0, outputPath=str(directory), meshVars=[p2]) + + restored = uw.discretisation.MeshVariable("restored", mesh, 3, degree=2) + restored.read_timestep("tet", "p2_tet", 0, outputPath=str(directory)) + np.testing.assert_allclose(restored.array, p2.array, atol=1.0e-12) + + if uw.mpi.rank == 0: + field_file = directory / "tet.mesh.p2_tet.00000.h5" + with h5py.File(field_file, "r") as handle: + coordinates = handle["fields/coordinates"][:] + values = handle["fields/p2_tet"][:] + cells = handle["fields/cells"][:] + assert cells.shape[1] == 10 + assert values.shape == coordinates.shape + np.testing.assert_allclose(values, coordinates, atol=1.0e-12) + points = coordinates[cells] + edge_pairs = ((0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)) + for local_node, (a, b) in enumerate(edge_pairs, start=4): + np.testing.assert_allclose( + points[:, local_node], 0.5 * (points[:, a] + points[:, b]) + ) + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + + xdmf = directory / "tet.mesh.00000.xdmf" + assert 'TopologyType="Tetrahedron_10"' in xdmf.read_text() + assert "&p2_tet_Data;:/fields/p2_tet" in xdmf.read_text() + _assert_xdmf_references_exist(xdmf) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_higher_order_fields_keep_exact_data_and_one_compact_reduction(tmp_path): + """P3+ uses P1 and DG2+ uses DG0 only for the XDMF view.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, regular=True, qdegree=5) + continuous = [ + uw.discretisation.MeshVariable(f"p{degree}", mesh, 1, degree=degree) for degree in (3, 4) + ] + discontinuous = [ + uw.discretisation.MeshVariable(f"dg{degree}", mesh, 1, degree=degree, continuous=False) + for degree in (2, 3, 4) + ] + for field in continuous: + field.array[:, 0, 0] = 2.0 + for field in discontinuous: + field.array[:, 0, 0] = 3.0 + mesh.write_timestep("high", 0, outputPath=str(directory), meshVars=continuous + discontinuous) + + if uw.mpi.rank == 0: + mesh_file = directory / "high.mesh.00000.h5" + with h5py.File(mesh_file, "r") as handle: + vertex_count = handle["geometry/vertices"].shape[0] + cell_count = handle["viz/topology/cells"].shape[0] + for field in continuous: + with h5py.File(directory / f"high.mesh.{field.clean_name}.00000.h5", "r") as handle: + assert handle[f"fields/{field.clean_name}"].shape[0] > vertex_count + assert handle[f"visualization/{field.clean_name}"].shape[0] == vertex_count + np.testing.assert_allclose(handle[f"visualization/{field.clean_name}"][:], 2.0) + assert ( + handle[f"visualization/{field.clean_name}"].attrs["visualization_degree"] == 1 + ) + _assert_no_compatibility_copies(handle) + for field in discontinuous: + with h5py.File(directory / f"high.mesh.{field.clean_name}.00000.h5", "r") as handle: + assert handle[f"fields/{field.clean_name}"].shape[0] > cell_count + assert handle[f"visualization/{field.clean_name}"].shape[0] == cell_count + np.testing.assert_allclose(handle[f"visualization/{field.clean_name}"][:], 3.0) + assert ( + handle[f"visualization/{field.clean_name}"].attrs["visualization_degree"] == 0 + ) + _assert_no_compatibility_copies(handle) + + xdmf = directory / "high.mesh.00000.xdmf" + text = xdmf.read_text() + for field in continuous + discontinuous: + assert f"&{field.clean_name}_Data;:/visualization/{field.clean_name}" in text + _assert_xdmf_references_exist(xdmf) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_compact_tensor_keeps_native_component_count(tmp_path): + """The authoritative field avoids a second nine-component tensor copy.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + tensor = uw.discretisation.MeshVariable("tensor", mesh, degree=1, vtype=uw.VarType.TENSOR) + tensor.array[:, 0, 0] = 1.0 + tensor.array[:, 0, 1] = 2.0 + tensor.array[:, 1, 0] = 3.0 + tensor.array[:, 1, 1] = 4.0 + mesh.write_timestep("tensor", 0, outputPath=str(directory), meshVars=[tensor]) + + if uw.mpi.rank == 0: + with h5py.File(directory / "tensor.mesh.tensor.00000.h5", "r") as handle: + assert handle["fields/tensor"].shape[1] == 4 + assert "visualization" not in handle + _assert_no_compatibility_copies(handle) + text = (directory / "tensor.mesh.00000.xdmf").read_text() + assert 'Type="Matrix"' in text + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_timestep_and_optional_checkpoint_roundtrip(tmp_path): + """Dimensional fields support remap; /restart/petsc supports exact reload.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.4, regular=True) + source = uw.discretisation.MeshVariable("source", mesh, 1, degree=2) + source.array[:, 0, 0] = source.coords[:, 0] + 2.0 * source.coords[:, 1] + same_layout_expected = np.array(source.array) mesh.write_timestep( - "noxdmf", index=0, outputPath=str(tmp_path), - meshVars=[s], create_xdmf=False, + "roundtrip", 0, outputPath=str(directory), meshVars=[source], petsc_reload=True ) - # XDMF file should not be generated. - # Note: we do NOT check HDF5 groups because some PETSc versions (3.21+) - # write /vertex_fields/ automatically during var.write() — that is PETSc - # behaviour, not ours. - xdmf_file = os.path.join(str(tmp_path), "noxdmf.mesh.00000.xdmf") - assert not os.path.exists(xdmf_file), "XDMF file should not exist" - - del mesh - - -# --------------------------------------------------------------------------- -# Test: Tensor variable repacking -# --------------------------------------------------------------------------- - - -def test_tensor_variable_repacking(tmp_path): - """Tensor variables are repacked to 9 components (3x3) in vertex_fields.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) - - # 2D full tensor (4 components internally) - T = uw.discretisation.MeshVariable( - "T", mesh, mesh.dim, degree=1, vtype=uw.VarType.TENSOR - ) - T.data[:, 0] = 1.0 # xx - T.data[:, 1] = 0.1 # xy - T.data[:, 2] = 0.2 # yx - T.data[:, 3] = 2.0 # yy - + field_file = directory / "roundtrip.mesh.source.00000.h5" + source.array[:, 0, 0] = 0.0 + source.read_checkpoint(str(field_file), data_name="source", same_layout=True) + np.testing.assert_allclose(source.array, same_layout_expected, atol=1.0e-12) + + reloaded_mesh = uw.discretisation.Mesh(str(directory / "roundtrip.mesh.00000.h5")) + restored = uw.discretisation.MeshVariable("source", reloaded_mesh, 1, degree=2) + restored.read_checkpoint(str(field_file), data_name="source") + expected = restored.coords[:, 0] + 2.0 * restored.coords[:, 1] + np.testing.assert_allclose(restored.array[:, 0, 0], expected, atol=1.0e-12) + + remapped = uw.discretisation.MeshVariable("remapped", mesh, 1, degree=2) + remapped.read_timestep("roundtrip", "source", 0, outputPath=str(directory)) + expected = remapped.coords[:, 0] + 2.0 * remapped.coords[:, 1] + np.testing.assert_allclose(remapped.array[:, 0, 0], expected, atol=1.0e-12) + + if uw.mpi.rank == 0: + with h5py.File(field_file, "r") as handle: + assert set(handle) == {"fields", "restart"} + assert "restart/petsc/topologies/uw_mesh/dms/source/vecs/source/source" in handle + assert set(handle["restart/petsc/topologies/uw_mesh/dms"]) == { + "source", + "uw_mesh", + } + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_create_xdmf_false_preserves_native_writer(tmp_path): + """Disabling XDMF leaves the established native field file untouched.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + field = uw.discretisation.MeshVariable("field", mesh, 1, degree=1) + field.array[:, 0, 0] = 1.0 + mesh.write_timestep("native", 0, outputPath=str(directory), meshVars=[field], create_xdmf=False) + if uw.mpi.rank == 0: + assert not (directory / "native.mesh.00000.xdmf").exists() + with h5py.File(directory / "native.mesh.field.00000.h5", "r") as handle: + assert "fields/field" in handle + assert "visualization" not in handle + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_petsc_reload_only_omits_duplicate_fields(tmp_path): + """Restart-only output stores the native PETSc payload exactly once.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + field = uw.discretisation.MeshVariable("field", mesh, 1, degree=2) + field.array[:, 0, 0] = field.coords[:, 0] + 2.0 * field.coords[:, 1] + expected = np.asarray(field.array).copy() mesh.write_timestep( - "tensor", index=0, outputPath=str(tmp_path), meshVars=[T] + "restart", + 0, + outputPath=str(directory), + meshVars=[field], + create_xdmf=False, + petsc_reload=True, ) - t_h5 = os.path.join(str(tmp_path), "tensor.mesh.T.00000.h5") - assert _check_h5_group_exists(t_h5, "vertex_fields/T_T") - - t_compat = _read_h5_dataset(t_h5, "vertex_fields/T_T") - n_verts = mesh._coords.shape[0] - t_ncomp = t_compat.size // n_verts - assert t_ncomp == 9, f"Tensor should be repacked to 9 components, got {t_ncomp}" - - # Verify XDMF - xdmf_file = os.path.join(str(tmp_path), "tensor.mesh.00000.xdmf") - _check_xdmf_refs(xdmf_file, str(tmp_path)) - - del mesh - - -# --------------------------------------------------------------------------- -# Test: Valid /viz/topology connectivity generation -# --------------------------------------------------------------------------- - - -def test_xdmf_viz_topology_written_correctly(tmp_path): - """write_timestep uses PETSc-native viz topology and valid geometry.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(3, 3)) - - # Write just the mesh (no vars needed to test mesh topology) - mesh.write_timestep("test_topo", index=0, outputPath=str(tmp_path)) - - mesh_h5 = os.path.join(str(tmp_path), "test_topo.mesh.00000.h5") - assert _check_h5_group_exists(mesh_h5, "viz/topology/cells"), ( - "Mesh HDF5 must contain the /viz/topology/cells dataset" - ) - assert _check_h5_group_exists(mesh_h5, "geometry/vertices"), ( - "Mesh HDF5 must contain the /geometry/vertices dataset" - ) - - # Validate cell-to-vertex connectivity bounds - cells = _read_h5_dataset(mesh_h5, "viz/topology/cells") - vertices = _read_h5_dataset(mesh_h5, "geometry/vertices") - - num_vertices = vertices.shape[0] - assert cells.max() < num_vertices, ( - f"Invalid topology: cells max ({cells.max()}) must be < numVertices ({num_vertices})" - ) - assert cells.min() >= 0, "Invalid topology: cells indices cannot be negative" - - # Validate XDMF actually points to the viz group, not the raw DMPlex group - xdmf_file = os.path.join(str(tmp_path), "test_topo.mesh.00000.xdmf") - assert os.path.exists(xdmf_file), "XDMF file should exist" - - with open(xdmf_file, "r") as f: - xdmf_content = f.read() - - assert "/viz/topology/cells" in xdmf_content, ( - "XDMF file should explicitly point to /viz/topology/cells" + field.array[...] = 0.0 + field.read_checkpoint( + str(directory / "restart.mesh.field.00000.h5"), + data_name="field", + same_layout=True, ) - assert "/geometry/vertices" in xdmf_content, ( - "XDMF file should explicitly point to /geometry/vertices" + np.testing.assert_allclose(field.array, expected) + + if uw.mpi.rank == 0: + assert not (directory / "restart.mesh.00000.xdmf").exists() + with h5py.File(directory / "restart.mesh.field.00000.h5", "r") as handle: + assert set(handle) == {"restart"} + assert set(handle["restart"]) == {"petsc"} + + # Existing files that used the former group name remain readable. + with h5py.File(directory / "restart.mesh.field.00000.h5", "a") as handle: + handle.move("restart/petsc", "uw_checkpoint") + del handle["restart"] + uw.mpi.barrier() + + field.array[...] = 0.0 + field.read_checkpoint( + str(directory / "restart.mesh.field.00000.h5"), + data_name="field", + same_layout=True, ) - - del mesh - - -def test_xdmf_viz_topology_3d_written_correctly(tmp_path): - """write_timestep uses PETSc-native viz topology for 3D meshes.""" - - mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2, 2)) - - # Write just the mesh - mesh.write_timestep("test_topo_3d", index=0, outputPath=str(tmp_path)) - - mesh_h5 = os.path.join(str(tmp_path), "test_topo_3d.mesh.00000.h5") - assert _check_h5_group_exists(mesh_h5, "viz/topology/cells"), ( - "3D Mesh HDF5 must contain the /viz/topology/cells dataset" - ) - - # Validate connectivity - cells = _read_h5_dataset(mesh_h5, "viz/topology/cells") - vertices = _read_h5_dataset(mesh_h5, "geometry/vertices") - - num_vertices = vertices.shape[0] - assert cells.max() < num_vertices, "Invalid 3D topology bounds" - - xdmf_file = os.path.join(str(tmp_path), "test_topo_3d.mesh.00000.xdmf") - assert os.path.exists(xdmf_file) - with open(xdmf_file, "r") as f: - xdmf_content = f.read() - assert "/viz/topology/cells" in xdmf_content - assert "/geometry/vertices" in xdmf_content - - del mesh + np.testing.assert_allclose(field.array, expected) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +@pytest.mark.parametrize("dim", [2, 3]) +def test_mesh_xdmf_topology_is_direct_and_bounded(tmp_path, dim): + """The base mesh XDMF uses valid direct cell-to-vertex connectivity.""" + directory = _shared_path(tmp_path) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2,) * dim) + mesh.write_timestep("topology", 0, outputPath=str(directory)) + if uw.mpi.rank == 0: + with h5py.File(directory / "topology.mesh.00000.h5", "r") as handle: + cells = handle["viz/topology/cells"][:] + vertices = handle["geometry/vertices"][:] + assert cells.ndim == 2 + assert cells.min() >= 0 + assert cells.max() < len(vertices) + text = (directory / "topology.mesh.00000.xdmf").read_text() + assert "&MeshData;:/viz/topology/cells" in text + assert "&MeshData;:/topology/cells" not in text diff --git a/tests/test_0005_xdmf_dg1.py b/tests/test_0005_xdmf_dg1.py new file mode 100644 index 000000000..51708b38b --- /dev/null +++ b/tests/test_0005_xdmf_dg1.py @@ -0,0 +1,140 @@ +"""Direct DG1 XDMF output preserves element-local affine traces.""" + +from pathlib import Path +import xml.etree.ElementTree as ET + +import h5py +import numpy as np +import pytest + +import underworld3 as uw + + +@pytest.mark.level_1 +@pytest.mark.tier_b +@pytest.mark.parametrize("dim", [2, 3]) +def test_dg1_simplex_output(tmp_path, dim): + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0,) * dim, + maxCoords=(1.0,) * dim, + cellSize=0.5, + regular=True, + qdegree=3, + ) + scalar = uw.discretisation.MeshVariable("dg_scalar", mesh, 1, degree=1, continuous=False) + tensor = uw.discretisation.MeshVariable( + "dg_tensor", mesh, degree=1, continuous=False, vtype=uw.VarType.TENSOR + ) + pressure = uw.discretisation.MeshVariable("pressure", mesh, 1, degree=1) + rows = mesh._cell_node_indices(1, False).reshape(-1, dim + 1) + coords = scalar.coords + offsets = np.floor(coords[rows].mean(axis=1)[:, 0] * 7 + 1.0e-8) + scalar.array[rows, 0, 0] = 1 + coords[rows, 0] + 2 * coords[rows, 1] + offsets[:, None] + tensor.array[:] = 0 + tensor.array[:, 0, 0] = scalar.array[:, 0, 0] + tensor.array[:, 0, 1] = 3 + coords[:, 0] + tensor.array[:, 1, 0] = -2 + coords[:, 1] + tensor.array[:, 1, 1] = 5 + pressure.array[:, 0, 0] = pressure.coords[:, 0] + mesh.write_timestep( + "fields", + 0, + outputPath=str(directory), + meshVars=[pressure, scalar, tensor], + petsc_reload=True, + ) + reloaded_mesh = uw.discretisation.Mesh(str(directory / "fields.mesh.00000.h5")) + restored = uw.discretisation.MeshVariable( + "dg_scalar", reloaded_mesh, 1, degree=1, continuous=False + ) + restored.read_checkpoint( + str(directory / "fields.mesh.dg_scalar.00000.h5"), data_name="dg_scalar" + ) + restored_rows = reloaded_mesh._cell_node_indices(1, False).reshape(-1, dim + 1) + restored_coords = restored.coords + restored_offsets = np.floor(restored_coords[restored_rows].mean(axis=1)[:, 0] * 7 + 1.0e-8) + expected_restored = ( + 1 + + restored_coords[restored_rows, 0] + + 2 * restored_coords[restored_rows, 1] + + restored_offsets[:, None] + ) + np.testing.assert_allclose( + restored.array[restored_rows, 0, 0], + expected_restored, + rtol=1.0e-12, + atol=1.0e-12, + ) + + remapped = uw.discretisation.MeshVariable("remapped", mesh, 1, degree=1, continuous=False) + remapped.read_timestep("fields", "dg_scalar", 0, outputPath=str(directory)) + np.testing.assert_allclose(remapped.array, scalar.array, rtol=1.0e-12, atol=1.0e-12) + + if uw.mpi.rank != 0: + return + with h5py.File(directory / "fields.mesh.dg_scalar.00000.h5", "r") as handle: + values = handle["fields/dg_scalar"][:].reshape(-1) + assert handle["fields/dg_scalar"].attrs["representation"] == "dg1_corner_nodal" + assert set(handle["fields"]) == {"dg_scalar"} + assert "visualization" not in handle + assert "dg1" not in handle + with h5py.File(directory / "fields.mesh.00000.h5", "r") as handle: + points = handle["viz/dg1/coordinates"][:] + cells = handle["viz/dg1/cells"][:] + assert len(cells) == len(handle["viz/topology/cells"]) + assert len(points) == len(cells) * (dim + 1) + assert len(values) == len(points) + assert len(np.unique(cells)) == len(points) + + corners = points[cells] + determinants = np.linalg.det((corners[:, 1:] - corners[:, :1]).transpose(0, 2, 1)) + assert np.all(determinants > 0) + centers = corners.mean(axis=1) + offsets = np.repeat(np.floor(centers[:, 0] * 7 + 1.0e-8), dim + 1) + expected = 1 + points[:, 0] + 2 * points[:, 1] + offsets + np.testing.assert_allclose(values, expected, rtol=1.0e-12, atol=1.0e-12) + + _, inverse = np.unique(np.round(points, 10), axis=0, return_inverse=True) + low = np.full(inverse.max() + 1, np.inf) + high = np.full(inverse.max() + 1, -np.inf) + np.minimum.at(low, inverse, values) + np.maximum.at(high, inverse, values) + assert np.max(high - low) >= 1 + + with h5py.File(directory / "fields.mesh.dg_tensor.00000.h5", "r") as handle: + tensor_values = handle["fields/dg_tensor"][:] + assert tensor_values.shape == (len(points), dim * dim) + np.testing.assert_allclose(tensor_values[:, 0], expected) + np.testing.assert_allclose(tensor_values[:, 1], 3 + points[:, 0]) + np.testing.assert_allclose(tensor_values[:, dim], -2 + points[:, 1], atol=1.0e-12) + np.testing.assert_allclose(tensor_values[:, dim + 1], 5) + + tree = ET.parse(directory / "fields.mesh.00000.xdmf") + grids = tree.findall(".//Grid[@GridType='Uniform']") + assert {grid.get("Name") for grid in grids} == {"domain", "DG1"} + grid = next(grid for grid in grids if grid.get("Name") == "DG1") + assert grid.find("Topology/DataItem").text.strip().endswith("/viz/dg1/cells") + assert grid.find("Geometry/DataItem").text.strip().endswith("/viz/dg1/coordinates") + attributes = {attribute.get("Name"): attribute for attribute in grid.findall("Attribute")} + assert set(attributes) == {"dg_scalar", "dg_tensor"} + for name, attribute in attributes.items(): + assert attribute.get("Center") == "Node" + assert f"/fields/{name}" in attribute.find("DataItem").text + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_unsupported_dg_layout_uses_dg0_visualization_reduction(tmp_path): + """Non-simplex DG1 remains exact in /fields and gets one DG0 XDMF array.""" + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2)) + dg = uw.discretisation.MeshVariable("dg", mesh, 1, degree=1, continuous=False) + dg.array[:, 0, 0] = 4.0 + mesh.write_timestep("structured", 0, outputPath=str(directory), meshVars=[dg]) + if uw.mpi.rank == 0: + with h5py.File(directory / "structured.mesh.dg.00000.h5", "r") as handle: + assert "fields/dg" in handle + assert "visualization/dg" in handle + np.testing.assert_allclose(handle["visualization/dg"][:], 4.0) + assert handle["visualization/dg"].attrs["visualization_degree"] == 0 diff --git a/tests/test_0005_xdmf_physical_units.py b/tests/test_0005_xdmf_physical_units.py new file mode 100644 index 000000000..3a3230c6d --- /dev/null +++ b/tests/test_0005_xdmf_physical_units.py @@ -0,0 +1,136 @@ +"""Dimensional /fields output and optional native PETSc checkpoints.""" + +from pathlib import Path + +import h5py +import numpy as np +import pytest + +import underworld3 as uw + + +def _set_reference_scales(): + model = uw.get_default_model() + model.set_scaling_mode("exact") + model.set_reference_quantities( + length=uw.quantity(10, "km"), + velocity=uw.quantity(5, "mm/year"), + pressure=uw.quantity(2, "MPa"), + ) + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_fields_are_dimensional_and_checkpoint_is_native(tmp_path): + """Analysis sees physical values while exact reload restores solver values.""" + _set_reference_scales() + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) + velocity = uw.discretisation.MeshVariable("velocity", mesh, mesh.dim, degree=2, units="mm/year") + pressure = uw.discretisation.MeshVariable( + "pressure", mesh, 1, degree=0, continuous=False, units="MPa" + ) + velocity.array[:, 0, 0] = uw.quantity(10.0, "mm/year") + velocity.array[:, 0, 1] = uw.quantity(15.0, "mm/year") + pressure.array[:, 0, 0] = uw.quantity(8.0, "MPa") + expected_velocity = np.array(velocity.array) + expected_pressure = np.array(pressure.array) + + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh.write_timestep( + "physical", + 0, + outputPath=str(directory), + meshVars=[velocity, pressure], + petsc_reload=True, + ) + + velocity_file = directory / "physical.mesh.velocity.00000.h5" + pressure_file = directory / "physical.mesh.pressure.00000.h5" + velocity_restart = uw.discretisation.MeshVariable( + "velocity_restart", mesh, mesh.dim, degree=2, units="mm/year" + ) + pressure_restart = uw.discretisation.MeshVariable( + "pressure_restart", mesh, 1, degree=0, continuous=False, units="MPa" + ) + velocity_restart.read_checkpoint(str(velocity_file), data_name="velocity") + pressure_restart.read_checkpoint(str(pressure_file), data_name="pressure") + np.testing.assert_allclose(np.array(velocity_restart.array), expected_velocity) + np.testing.assert_allclose(np.array(pressure_restart.array), expected_pressure) + + velocity_remap = uw.discretisation.MeshVariable( + "velocity_remap", mesh, mesh.dim, degree=2, units="mm/year" + ) + velocity_remap.read_timestep("physical", "velocity", 0, outputPath=str(directory)) + np.testing.assert_allclose(np.array(velocity_remap.array), expected_velocity) + + if uw.mpi.rank != 0: + return + with h5py.File(directory / "physical.mesh.00000.h5", "r") as handle: + native_coordinates = handle["geometry/vertices"][:] + physical_coordinates = handle["viz/geometry/vertices"][:] + np.testing.assert_allclose(physical_coordinates, native_coordinates * 10.0) + assert handle["viz/geometry/vertices"].attrs["units"] == "kilometer" + + with h5py.File(velocity_file, "r") as handle: + physical_velocity = handle["fields/velocity"][:] + np.testing.assert_allclose(physical_velocity[:, 0], 10.0) + np.testing.assert_allclose(physical_velocity[:, 1], 15.0) + assert np.isclose(handle["fields/coordinates"][:].max(), 10.0) + assert handle["fields/velocity"].attrs["units"] == "millimeter / year" + assert handle["fields/coordinates"].attrs["units"] == "kilometer" + assert "visualization" not in handle + assert "vertex_fields" not in handle + assert "restart/petsc/topologies/uw_mesh/dms/velocity/vecs/velocity/velocity" in handle + + with h5py.File(pressure_file, "r") as handle: + np.testing.assert_allclose(handle["fields/pressure"][:], 8.0) + assert handle["fields/pressure"].attrs["units"] == "megapascal" + assert "cell_fields" not in handle + + text = (directory / "physical.mesh.00000.xdmf").read_text() + assert "&velocity_Data;:/fields/velocity" in text + assert "&pressure_Data;:/fields/pressure" in text + assert 'Information Name="Units" Value="millimeter / year"' in text + assert 'Information Name="Units" Value="megapascal"' in text + assert 'Information Name="Units" Value="kilometer"' in text + + +@pytest.mark.level_1 +@pytest.mark.tier_b +def test_dg1_fields_are_physical_at_element_corners(tmp_path): + """DG1 uses one physical corner field and shared mesh geometry.""" + _set_reference_scales() + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.5, regular=True) + pressure = uw.discretisation.MeshVariable( + "dg_pressure", mesh, 1, degree=1, continuous=False, units="MPa" + ) + pressure.array[:, 0, 0] = uw.quantity(8.0, "MPa") + directory = Path(uw.mpi.comm.bcast(str(tmp_path), root=0)) + mesh.write_timestep("dg", 0, outputPath=str(directory), meshVars=[pressure]) + + remapped = uw.discretisation.MeshVariable( + "dg_pressure_remapped", + mesh, + 1, + degree=1, + continuous=False, + units="MPa", + ) + remapped.read_timestep("dg", "dg_pressure", 0, outputPath=str(directory)) + np.testing.assert_allclose(np.array(remapped.array), np.array(pressure.array)) + + if uw.mpi.rank == 0: + with h5py.File(directory / "dg.mesh.dg_pressure.00000.h5", "r") as handle: + np.testing.assert_allclose(handle["fields/dg_pressure"][:], 8.0) + assert handle["fields/dg_pressure"].attrs["units"] == "megapascal" + assert ( + handle["fields/dg_pressure"].attrs["representation"] + == "dg1_corner_nodal" + ) + assert set(handle["fields"]) == {"dg_pressure"} + assert "visualization" not in handle + assert "dg1" not in handle + assert "restart" not in handle + with h5py.File(directory / "dg.mesh.00000.h5", "r") as handle: + assert np.isclose(handle["viz/dg1/coordinates"][:].max(), 10.0) + assert handle["viz/dg1/coordinates"].attrs["units"] == "kilometer" diff --git a/tests/test_0010_snapshot_disk_format.py b/tests/test_0010_snapshot_disk_format.py index 0b016e246..518c13431 100644 --- a/tests/test_0010_snapshot_disk_format.py +++ b/tests/test_0010_snapshot_disk_format.py @@ -221,8 +221,9 @@ def test_write_snapshot_produces_wrapper_and_bulk_dir(tmp_path): if filename.endswith(f".{variable_name}.00000.h5") ) with h5py.File(os.path.join(bulk, variable_file), "r") as h5: - assert variable_name in h5["topologies"]["uw_mesh"]["dms"] - assert variable_name in h5["uw_checkpoint"] + dms = h5["restart/petsc/topologies/uw_mesh/dms"] + assert variable_name in dms + assert f"{variable_name}/vecs/{variable_name}/{variable_name}" in dms def test_snapshot_bulk_filenames_do_not_expand_loaded_mesh_name(tmp_path):