From 61d3766afbc19b27b76c862d9602b047859e0e3e Mon Sep 17 00:00:00 2001 From: Tyagi Date: Wed, 25 Jun 2025 14:56:05 +1000 Subject: [PATCH 1/5] updated checkpoint_xdmf to look for correct vertex field name in h5 file --- src/underworld3/discretisation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/underworld3/discretisation.py b/src/underworld3/discretisation.py index 61d5a93a..c134c6b4 100644 --- a/src/underworld3/discretisation.py +++ b/src/underworld3/discretisation.py @@ -2810,7 +2810,7 @@ def checkpoint_xdmf( DataType="Float" Precision="8" Dimensions="1 {numVertices} {var.num_components}" Format="HDF"> - &{var.clean_name+"_Data"};:/vertex_fields/{var.clean_name+"_P"+str(var.degree)} + &{var.clean_name+"_Data"};:/vertex_fields/{var.clean_name+"_"+var.clean_name} From 27f3da10f40a5bf611d2940e4bffc757898babc4 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Wed, 25 Jun 2025 17:14:30 +1000 Subject: [PATCH 2/5] added test to check vertex fields in xdmf and h5 files --- tests/test_0005_check_xdmf.py | 188 ++++++++++++++++++++++++++++++++++ 1 file changed, 188 insertions(+) create mode 100644 tests/test_0005_check_xdmf.py diff --git a/tests/test_0005_check_xdmf.py b/tests/test_0005_check_xdmf.py new file mode 100644 index 00000000..f11271b3 --- /dev/null +++ b/tests/test_0005_check_xdmf.py @@ -0,0 +1,188 @@ +import pytest +import sympy +import underworld3 as uw + +import re +import h5py +import glob + +import os +import numpy as np + +# + +structured_quad_box = uw.meshing.StructuredQuadBox(elementRes=(3,) * 2) + +unstructured_quad_box_irregular = uw.meshing.UnstructuredSimplexBox( + cellSize=0.2, regular=False, qdegree=2, refinement=1 +) +unstructured_quad_box_regular = uw.meshing.UnstructuredSimplexBox( + cellSize=0.2, regular=True, qdegree=2, refinement=2 +) + +unstructured_quad_box_irregular_3D = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0, 0.0), + maxCoords=(1.0, 1.0, 1.0), + cellSize=0.25, + regular=False, + qdegree=2, +) + + +# - + +def check_xdmf_vertex_fields_exist_in_h5(xdmf_filename, tmp_path=""): + """ + Checks that all vertex field datasets referenced in the XDMF file + exist in their corresponding HDF5 files in tmp_path. + Prints a summary for each check. + """ + # 1. Extract ENTITY->h5 mapping and XDMF content + with open(xdmf_filename, 'r') as f: + content = f.read() + doctype_match = re.search(r'', content, re.DOTALL) + if not doctype_match: + print("No DOCTYPE entity block found.") + return + entity_block = doctype_match.group(1) + entities = dict(re.findall(r'', entity_block)) + + # 2. Extract (&Entity;/vertex_fields/field_name) references + refs = re.findall(r'&(\w+);:(/vertex_fields/[A-Za-z0-9_]+)', content) + print("Checking vertex field dataset references in XDMF:") + for entity_name, dataset_path in refs: + h5_file = entities.get(entity_name) + if not h5_file: + print(f"[ENTITY NOT FOUND] {entity_name}: {dataset_path}") + continue + h5_full_path = os.path.join(tmp_path, h5_file) + try: + with h5py.File(h5_full_path, 'r') as f: + h5_path = dataset_path.lstrip('/') + if h5_path in f: + print(f"[OK] {h5_file}: {dataset_path} found") + else: + print(f"[MISSING] {h5_file}: {dataset_path} not found") + except OSError as e: + print(f"[ERROR] Cannot open {h5_file}: {e}") + + +def remove_test_mesh_files(directory='.'): + """Delete all files starting with test.mesh. in the specified directory.""" + pattern = os.path.join(directory, 'test.mesh.*') + for file_path in glob.glob(pattern): + try: + os.remove(file_path) + print(f"Removed: {file_path}") + except Exception as e: + print(f"Could not remove {file_path}: {e}") + + +@pytest.mark.parametrize( + "mesh", + [ + structured_quad_box, + unstructured_quad_box_irregular, + unstructured_quad_box_regular, + unstructured_quad_box_irregular_3D, + ], +) +def test_stokes_boxmesh(mesh, tmp_path): + print(f"Mesh - Coordinates: {mesh.CoordinateSystem.type}") + mesh.dm.view() + + if mesh.dim == 2: + x, y = mesh.X + else: + x, y, z = mesh.X + + u = uw.discretisation.MeshVariable( + 'u', mesh, mesh.dim, vtype=uw.VarType.VECTOR, degree=2 + ) + p = uw.discretisation.MeshVariable( + 'p', mesh, 1, vtype=uw.VarType.SCALAR, degree=1 + ) + u2 = uw.discretisation.MeshVariable( + 'u2', mesh, mesh.dim, vtype=uw.VarType.VECTOR, degree=2 + ) + p2 = uw.discretisation.MeshVariable( + 'p2', mesh, 1, vtype=uw.VarType.SCALAR, degree=1 + ) + + stokes = uw.systems.Stokes(mesh, velocityField=u, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1 + + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + + stokes.petsc_options["snes_type"] = "newtonls" + stokes.petsc_options["ksp_type"] = "fgmres" + stokes.petsc_options["ksp_monitor"] = None + stokes.petsc_options["snes_monitor"] = None + stokes.tolerance = 1.0e-3 + + # stokes.petsc_options.setValue("fieldsplit_velocity_pc_type", "mg") + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade") + stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w") + + stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" + stokes.petsc_options[f"fieldsplit_velocity_ksp_type"] = "fcg" + stokes.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" + stokes.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_max_it"] = 7 + stokes.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None + + stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "gamg") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative") + stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v") + + if mesh.dim == 2: + stokes.bodyforce = 1.0e6 * sympy.Matrix([0, x]) + + stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((0.0, None), "Top") + + stokes.add_dirichlet_bc((0.0, None), "Left") + stokes.add_condition(conds=(0.0, None), + label="Right", + f_id=0, + c_type='dirichlet' ) + else: + stokes.bodyforce = 1.0e6 * sympy.Matrix([0, x, 0]) + + stokes.add_dirichlet_bc((0.0, 0.0, 0.0), "Bottom") + stokes.add_dirichlet_bc((0.0, 0.0, 0.0), "Top") + + stokes.add_dirichlet_bc((0.0, None, sympy.oo), "Left") + stokes.add_dirichlet_bc((0.0, sympy.oo, None), "Right") + + stokes.add_dirichlet_bc((sympy.oo, 0.0, sympy.oo), "Front") + stokes.add_dirichlet_bc((sympy.oo, 0.0, sympy.oo), "Back") + + stokes.solve() + + print(f"Mesh dimensions {mesh.dim}", flush=True) + # stokes.dm.ds.view() + + assert stokes.snes.getConvergedReason() > 0 + + mesh.write_timestep("test", meshUpdates=False, meshVars=[u, p], outputPath=tmp_path, index=0) + + # Call XDMF/HDF5 checker (assume xdmf file is named "test.mesh.xdmf" and written in tmp_path) + xdmf_filename = os.path.join(tmp_path, "test.mesh.00000.xdmf") + check_xdmf_vertex_fields_exist_in_h5(xdmf_filename, tmp_path=str(tmp_path)) + + u2.read_timestep("test", "u", 0, outputPath=tmp_path) + p2.read_timestep("test", "p", 0, outputPath=tmp_path) + + with mesh.access(): + assert np.allclose(u.data, u2.data) + assert np.allclose(p.data, p2.data) + + remove_test_mesh_files(directory=tmp_path) + + del mesh + del stokes + print('----------------------------------------------------------------------------------------') + return + + From 56e443e0902024dc01889177f6001c6c3d133a39 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Wed, 25 Jun 2025 17:25:19 +1000 Subject: [PATCH 3/5] minor updates to the test --- tests/test_0005_check_xdmf.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/tests/test_0005_check_xdmf.py b/tests/test_0005_check_xdmf.py index f11271b3..19dcd122 100644 --- a/tests/test_0005_check_xdmf.py +++ b/tests/test_0005_check_xdmf.py @@ -31,28 +31,22 @@ # - def check_xdmf_vertex_fields_exist_in_h5(xdmf_filename, tmp_path=""): - """ - Checks that all vertex field datasets referenced in the XDMF file - exist in their corresponding HDF5 files in tmp_path. - Prints a summary for each check. - """ - # 1. Extract ENTITY->h5 mapping and XDMF content + errors = [] with open(xdmf_filename, 'r') as f: content = f.read() doctype_match = re.search(r'', content, re.DOTALL) if not doctype_match: - print("No DOCTYPE entity block found.") - return + raise AssertionError("No DOCTYPE entity block found in XDMF file.") entity_block = doctype_match.group(1) entities = dict(re.findall(r'', entity_block)) - - # 2. Extract (&Entity;/vertex_fields/field_name) references refs = re.findall(r'&(\w+);:(/vertex_fields/[A-Za-z0-9_]+)', content) print("Checking vertex field dataset references in XDMF:") for entity_name, dataset_path in refs: h5_file = entities.get(entity_name) if not h5_file: - print(f"[ENTITY NOT FOUND] {entity_name}: {dataset_path}") + err = f"[ENTITY NOT FOUND] {entity_name}: {dataset_path}" + print(err) + errors.append(err) continue h5_full_path = os.path.join(tmp_path, h5_file) try: @@ -61,9 +55,15 @@ def check_xdmf_vertex_fields_exist_in_h5(xdmf_filename, tmp_path=""): if h5_path in f: print(f"[OK] {h5_file}: {dataset_path} found") else: - print(f"[MISSING] {h5_file}: {dataset_path} not found") + err = f"[MISSING] {h5_file}: {dataset_path} not found" + print(err) + errors.append(err) except OSError as e: - print(f"[ERROR] Cannot open {h5_file}: {e}") + err = f"[ERROR] Cannot open {h5_file}: {e}" + print(err) + errors.append(err) + if errors: + raise AssertionError("Missing or inaccessible vertex fields:\n" + "\n".join(errors)) def remove_test_mesh_files(directory='.'): From b149867ab029c12b57d9ab9d42ee438c2532b816 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Thu, 26 Jun 2025 16:52:40 +1000 Subject: [PATCH 4/5] updated checkpoint_xdmf method save discontinuous fields properly --- src/underworld3/discretisation.py | 34 +++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/src/underworld3/discretisation.py b/src/underworld3/discretisation.py index c134c6b4..cba80e77 100644 --- a/src/underworld3/discretisation.py +++ b/src/underworld3/discretisation.py @@ -2784,39 +2784,55 @@ def checkpoint_xdmf( attributes = "" for var in meshVars: - var_filename = filename + f"mesh.{var.clean_name}.{index:05}.h5" + var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5" + + def get_cell_field_shape(h5_filename): + with h5py.File(h5_filename, 'r') as f: + shape = f[f'cell_fields/{var.clean_name}_{var.clean_name}'].shape[0] + return shape + if var.num_components == 1: variable_type = "Scalar" else: variable_type = "Vector" - # We should add a tensor type here ... + + # Determine if data is stored on nodes (vertex_fields) or cells (cell_fields) + if not getattr(var, "continuous") or getattr(var, "degree")==0: + center = "Cell" + numItems = get_cell_field_shape(var_filename) + field_group = "cell_fields" + else: + center = "Node" + numItems = numVertices + field_group = "vertex_fields" var_attribute = f""" + Center="{center}"> + Dimensions="1 {numItems} {var.num_components}" + Type="HyperSlab"> 0 0 0 1 1 1 - 1 {numVertices} {var.num_components} + 1 {numItems} {var.num_components} - &{var.clean_name+"_Data"};:/vertex_fields/{var.clean_name+"_"+var.clean_name} + &{var.clean_name+"_Data"};:/{field_group}/{var.clean_name+"_"+var.clean_name} - """ + """ attributes += var_attribute + for var in swarmVars: var_filename = filename + f".proxy.{var.clean_name}.{index:05}.h5" if var.num_components == 1: From b6a24080b1df8fcae99d6a85c22b3be7f62eae80 Mon Sep 17 00:00:00 2001 From: Julian Giordani Date: Thu, 26 Jun 2025 17:05:12 +1000 Subject: [PATCH 5/5] Update tests/test_0005_check_xdmf.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- tests/test_0005_check_xdmf.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_0005_check_xdmf.py b/tests/test_0005_check_xdmf.py index 19dcd122..a519b659 100644 --- a/tests/test_0005_check_xdmf.py +++ b/tests/test_0005_check_xdmf.py @@ -115,8 +115,6 @@ def test_stokes_boxmesh(mesh, tmp_path): stokes.petsc_options["snes_type"] = "newtonls" stokes.petsc_options["ksp_type"] = "fgmres" - stokes.petsc_options["snes_type"] = "newtonls" - stokes.petsc_options["ksp_type"] = "fgmres" stokes.petsc_options["ksp_monitor"] = None stokes.petsc_options["snes_monitor"] = None stokes.tolerance = 1.0e-3