diff --git a/.gitignore b/.gitignore index e0d6391..9e73f54 100644 --- a/.gitignore +++ b/.gitignore @@ -534,6 +534,7 @@ tests/cases/geomslice_data_fds*/ tests/cases/cface_data_fds*/ tests/cases/slcf_partial_fds*/ tests/cases/slcf_masks_fds*/ +tests/cases/bndf_global_fds*/ # Pickle cache files *.pickle diff --git a/fdsreader/bndf/obstruction.py b/fdsreader/bndf/obstruction.py index d82bd31..cb6ca7c 100644 --- a/fdsreader/bndf/obstruction.py +++ b/fdsreader/bndf/obstruction.py @@ -92,7 +92,8 @@ def get_coordinates(self, ignore_cell_centered: bool = False) -> Dict[Literal["x if self.cell_centered and not ignore_cell_centered: co = co[:-1] + np.diff(co) / 2 - coords[dim] = co[np.where(np.logical_and(co >= self.extent[dim][0], co <= self.extent[dim][1]))] + tol = _tolerance(self.mesh.coordinates[dim]) + coords[dim] = co[(co >= self.extent[dim][0] - tol) & (co <= self.extent[dim][1] + tol)] if coords[dim].size == 0: coords[dim] = np.array([co[np.argmin(np.abs(co - self.extent[dim][0]))]]) @@ -696,163 +697,47 @@ def has_boundary_data(self): def get_global_boundary_data_arrays( self, quantity: Union[str, Quantity] - ) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]]: + ) -> Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]]: """Creates a global numpy ndarray from all subobstruction's boundary data for each orientation. + Each array has the shape (n_t, a, b), where a and b are the two axes of the face in the order x, y, z (e.g. y + and z for the orientations -1 and 1). The grid uses the finest cell size of the meshes the obstruction spans; + values of coarser meshes are repeated, and points not covered by any patch are NaN. The grid is uniform, so + on stretched meshes values are placed at the nearest grid point below their node. Frames a mesh did not write + are NaN as well. + + If the face lies at different positions in different meshes, e.g. because the obstruction spans meshes with + different cell sizes and FDS snaps its bounds to each grid, the orientation maps to a tuple with one array per + position, ordered by position. + :param quantity: The quantity's name or short name for which to load the global arrays. """ - if not self.has_boundary_data: - return dict() + patches_by_orientation: Dict[int, List[Patch]] = dict() for subobst in self._all_subobstructions: - if subobst.has_boundary_data: - cell_centered = next(iter(subobst._boundary_data.values())).cell_centered + try: + boundary = subobst.get_data(quantity) + except KeyError: + continue + for orientation, patch in boundary.data.items(): + patches_by_orientation.setdefault(orientation, list()).append(patch) + if not patches_by_orientation and self.has_boundary_data: + raise KeyError(f"No boundary data found for quantity '{quantity}'") result = dict() - for orientation_int in self.orientations: - subobst_sets = [list(), list()] - dim = ["x", "y", "z"][abs(orientation_int) - 1] - random_subobst = next( - (subobst for subobst in self._all_subobstructions if orientation_int in subobst.get_coordinates()), - None, - ) - if random_subobst is None: - raise KeyError(f"No sub-obstruction found for orientation {orientation_int}") - base_coord = random_subobst.get_coordinates(ignore_cell_centered=False)[orientation_int][dim][0] - - for subobst in self._all_subobstructions: - if not subobst.has_boundary_data: - continue - coords = subobst.get_coordinates(ignore_cell_centered=False) - if orientation_int not in coords: - continue - if coords[orientation_int][dim][0] == base_coord: - subobst_sets[0].append(subobst) + for orientation in sorted(patches_by_orientation): + normal = ("x", "y", "z")[abs(orientation) - 1] + levels: List[Tuple[float, List[Patch]]] = list() + for patch in patches_by_orientation[orientation]: + position = patch.extent[normal][0] + tol = _tolerance(patch.mesh.coordinates[normal]) + level = next((lvl for lvl in levels if abs(lvl[0] - position) <= tol), None) + if level is None: + levels.append((position, [patch])) else: - subobst_sets[1].append(subobst) - if len(subobst_sets[1]) == 0: - subobst_sets.pop(1) - - first_result_grid = None - for subobsts in subobst_sets: - coord_min = {"x": math.inf, "y": math.inf, "z": math.inf} - coord_max = {"x": -math.inf, "y": -math.inf, "z": -math.inf} - for dim in ("x", "y", "z"): - for subobst in subobsts: - patch_extent = subobst.get_data(quantity).data[orientation_int].extent - co = subobst.get_coordinates(ignore_cell_centered=True)[orientation_int][dim] - co = co[ - np.where( - np.logical_and( - co >= patch_extent[dim][0], - co <= patch_extent[dim][1], - ) - ) - ] - coord_min[dim] = min(co[0], coord_min[dim]) - coord_max[dim] = max(co[-1], coord_max[dim]) - - # The global grid will use the finest mesh as base and duplicate values of the coarser meshes - # Therefore we first find the finest mesh and calculate the step size in each dimension - step_sizes_min = { - "x": coord_max["x"] - coord_min["x"], - "y": coord_max["y"] - coord_min["y"], - "z": coord_max["z"] - coord_min["z"], - } - step_sizes_max = {"x": 0, "y": 0, "z": 0} - steps = dict() - global_max = {"x": -math.inf, "y": -math.inf, "z": -math.inf} - - for dim in ("x", "y", "z"): - for subobst in subobsts: - subobst_coords = subobst.get_coordinates(ignore_cell_centered=True)[orientation_int] - if len(subobst_coords[dim]) <= 1: - step_size = 0 - else: - step_size = subobst_coords[dim][1] - subobst_coords[dim][0] - step_sizes_min[dim] = min(step_size, step_sizes_min[dim]) - step_sizes_max[dim] = max(step_size, step_sizes_max[dim]) - global_max[dim] = max(subobst_coords[dim][-1], global_max[dim]) - - for dim in ("x", "y", "z"): - if step_sizes_min[dim] == 0: - step_sizes_min[dim] = math.inf - steps[dim] = 1 - else: - steps[dim] = max(int(round((coord_max[dim] - coord_min[dim]) / step_sizes_min[dim])), 1) + ( - 0 if cell_centered else 1 - ) - - grid = np.full((self.n_t(count_duplicates=False), steps["x"], steps["y"], steps["z"]), np.nan) - - start_coordinates = {"x": coord_min["x"], "y": coord_min["y"], "z": coord_min["z"]} - start_idx = dict() - end_idx = dict() - for subobst in subobsts: - patch_data = np.expand_dims( - subobst.get_data(quantity).data[orientation_int].data, axis=abs(orientation_int) - ) - subobst_coords = subobst.get_coordinates(ignore_cell_centered=True)[orientation_int] - for axis in (0, 1, 2): - dim = ("x", "y", "z")[axis] - if axis == abs(orientation_int) - 1: - start_idx[dim] = 0 - end_idx[dim] = 1 - continue - n_repeat = max( - int(round((subobst_coords[dim][1] - subobst_coords[dim][0]) / step_sizes_min[dim])), 1 - ) - - start_idx[dim] = int( - round((subobst_coords[dim][0] - start_coordinates[dim]) / step_sizes_min[dim]) - ) - end_idx[dim] = int( - round((subobst_coords[dim][-1] - start_coordinates[dim]) / step_sizes_min[dim]) - ) - - # We ignore border points unless they are actually on the border of the simulation space as all - # other border points actually appear twice, as the subobstructions overlap. This only - # applies for face_centered data, as cell_centered data will not overlap. - if not cell_centered: - if axis != abs(orientation_int) - 1: - reduced_shape = list(patch_data.shape) - reduced_shape[axis + 1] -= 1 - reduced_data_slices = tuple(slice(s) for s in reduced_shape) - patch_data = patch_data[reduced_data_slices] - - # Temporarily save border points to add them back to the array again later - if subobst_coords[dim][-1] == global_max[dim]: - end_idx[dim] += 1 - temp_data_slices = [slice(s) for s in patch_data.shape] - temp_data_slices[axis + 1] = slice(patch_data.shape[axis + 1] - 1, None) - temp_data = patch_data[tuple(temp_data_slices)] - - if n_repeat > 1: - patch_data = np.repeat(patch_data, n_repeat, axis=axis + 1) - - # Add border points back again if needed - if not cell_centered and subobst_coords[dim][-1] == global_max[dim]: - patch_data = np.concatenate((patch_data, temp_data), axis=axis + 1) - - grid[ - :, start_idx["x"] : end_idx["x"], start_idx["y"] : end_idx["y"], start_idx["z"] : end_idx["z"] - ] = patch_data.reshape( - ( - self.n_t(count_duplicates=False), - end_idx["x"] - start_idx["x"], - end_idx["y"] - start_idx["y"], - end_idx["z"] - start_idx["z"], - ) - ) - - # Remove empty dimensions, but make sure to note remove the time dimension if there is only a single - # timestep - grid = np.squeeze(grid) - if len(grid.shape) == 2: - grid = grid[np.newaxis, :, :] - if first_result_grid is not None: - result[orientation_int] = (first_result_grid, grid) - else: - result[orientation_int] = grid + level[1].append(patch) + levels.sort(key=lambda lvl: lvl[0]) + grids = tuple(_global_patch_grid(level_patches, normal) for _, level_patches in levels) + result[orientation] = grids[0] if len(grids) == 1 else grids return result def clear_cache(self): @@ -915,3 +800,53 @@ def __repr__(self, *args, **kwargs): + (f", Quantities={[q.short_name for q in self.quantities]}" if self.has_boundary_data else "") + ")" ) + + +def _tolerance(coordinates: np.ndarray) -> float: + """Tolerance for comparing positions on a grid: a thousandth of its smallest cell.""" + return 1e-3 * float(np.min(np.diff(coordinates))) if coordinates.size > 1 else 1e-6 + + +def _global_patch_grid(patches: List[Patch], normal: Literal["x", "y", "z"]) -> np.ndarray: + """Combines patches lying in one plane into one array of shape (n_t, a, b) on the finest grid of the patches.""" + cell_centered = patches[0].cell_centered + axes = [dim for dim in ("x", "y", "z") if dim != normal] + nodes = {dim: [patch.get_coordinates(ignore_cell_centered=True)[dim] for patch in patches] for dim in axes} + + positions = dict() + for dim in axes: + lower = min(n[0] for n in nodes[dim]) + upper = max(n[-1] for n in nodes[dim]) + steps = [np.min(np.diff(n)) for n in nodes[dim] if n.size > 1] + step = min(steps) if steps else 0.0 + n_nodes = int(round((upper - lower) / step)) + 1 if step > 0 else 1 + positions[dim] = lower + step * np.arange(n_nodes) + if cell_centered: + positions[dim] = positions[dim][:-1] + step / 2 + + n_t = max(patch.data.shape[0] for patch in patches) + grid = np.full((n_t,) + tuple(positions[dim].size for dim in axes), np.nan) + # On a node shared by two patches, the patch further up in x, y and z wins, as each is written after those below + order = sorted(range(len(patches)), key=lambda i: tuple(nodes[dim][i][0] for dim in axes)) + for i in order: + patch = patches[i] + global_indices = list() + patch_indices = list() + for dim in axes: + patch_nodes = nodes[dim][i] + tol = _tolerance(patch_nodes) + pos = positions[dim] + if cell_centered: + inside = (pos > patch_nodes[0]) & (pos < patch_nodes[-1]) + n_values = patch_nodes.size - 1 + else: + inside = (pos >= patch_nodes[0] - tol) & (pos <= patch_nodes[-1] + tol) + n_values = patch_nodes.size + # Each point of the global grid takes the value of the patch's cell it lies in, or of the patch's last + # node at or below it, so values of coarser meshes are repeated + idx = np.searchsorted(patch_nodes, pos[inside] + tol, side="right") - 1 + global_indices.append(np.flatnonzero(inside)) + patch_indices.append(np.clip(idx, 0, n_values - 1)) + data = patch.data + grid[(slice(data.shape[0]),) + np.ix_(*global_indices)] = data[(slice(None),) + np.ix_(*patch_indices)] + return grid diff --git a/fdsreader/export/obst_exporter.py b/fdsreader/export/obst_exporter.py index c2b75a5..1b1a508 100644 --- a/fdsreader/export/obst_exporter.py +++ b/fdsreader/export/obst_exporter.py @@ -56,9 +56,7 @@ def export_obst_raw(obst: Obstruction, output_dir: str, ordering: Literal["C", " } ) - def worker( - quantity: str, faces: Dict[int, Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]], vmin: float, vmax: float - ): + def worker(quantity: str, faces: Dict[int, Union[np.ndarray, Tuple[np.ndarray, ...]]], vmin: float, vmax: float): # Abort if no useful data is available if vmax <= 0: return @@ -76,10 +74,12 @@ def worker( with open(os.path.join(output_dir, quantity_name, filename), "wb") as rawfile: for face in faces.values(): # face for each orientation - for d in (face * out["ScaleFactor"]).astype(np.uint8): - if ordering == "F": - d = d.T - d.tofile(rawfile) + # A face lying at several positions (see get_global_boundary_data_arrays) is written part by part + for part in face if isinstance(face, tuple) else (face,): + for d in (part * out["ScaleFactor"]).astype(np.uint8): + if ordering == "F": + d = d.T + d.tofile(rawfile) with lock: meta["Quantities"].append(out) diff --git a/tests/cases/bndf_global_fds611.tgz b/tests/cases/bndf_global_fds611.tgz new file mode 100644 index 0000000..027951a --- /dev/null +++ b/tests/cases/bndf_global_fds611.tgz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6e502ffcf16979705835eb7fe011e7fba41bfae95ddf61390961a1a10fd331d1 +size 25790 diff --git a/tests/cases/fds_inputs/bndf_global.fds b/tests/cases/fds_inputs/bndf_global.fds new file mode 100644 index 0000000..73b9cc7 --- /dev/null +++ b/tests/cases/fds_inputs/bndf_global.fds @@ -0,0 +1,19 @@ +&HEAD CHID='bndf_global', TITLE='Global boundary data arrays for fdsreader' / +&MESH ID='Fine', IJK=10,10,10, XB=0.0,1.0,0.0,1.0,0.0,1.0 / +&MESH ID='Coarse', IJK=5,5,5, XB=1.0,2.0,0.0,1.0,0.0,1.0 / +&TIME T_END=2.0 / +&DUMP DT_BNDF=0.5, SMOKE3D=F / +&REAC FUEL='PROPANE' / +&SURF ID='FIRE', HRRPUA=500.0 / +&VENT XB=0.1,0.5,0.1,0.9,0.0,0.0, SURF_ID='FIRE' / +&VENT MB='XMIN', SURF_ID='OPEN' / &VENT MB='XMAX', SURF_ID='OPEN' / +&VENT MB='ZMAX', SURF_ID='OPEN' / + +Across both meshes: the top at z = 0.53 snaps to 0.5 in 'Fine' and to 0.6 in 'Coarse' +&OBST ID='Table', XB=0.6,1.4,0.2,0.8,0.0,0.53 / +One cell wide in x +&OBST ID='Strip', XB=0.2,0.3,0.2,0.8,0.0,0.3 / + +&BNDF QUANTITY='WALL TEMPERATURE' / +&BNDF QUANTITY='NET HEAT FLUX', CELL_CENTERED=T / +&TAIL / diff --git a/tests/test_boundary_global_arrays.py b/tests/test_boundary_global_arrays.py new file mode 100644 index 0000000..9c102a2 --- /dev/null +++ b/tests/test_boundary_global_arrays.py @@ -0,0 +1,88 @@ +"""Obstruction.get_global_boundary_data_arrays on bndf_global_fds611 (FDS 6.11.1; input in +tests/cases/fds_inputs/bndf_global.fds). + +Two meshes: 'Fine' (x 0..1, 0.1 m cells) and 'Coarse' (x 1..2, 0.2 m cells). 'Table' spans both; its top at +z = 0.53 snaps to 0.5 in 'Fine' and to 0.6 in 'Coarse'. 'Strip' is one cell wide in x. Each face has the node-based +WALL TEMPERATURE and the cell-centred TOTAL HEAT FLUX. +""" + +import numpy as np +import pytest +from fds_cases import case_dirs, case_id, load + +NODES = "WALL TEMPERATURE" +CELLS = "TOTAL HEAT FLUX" + + +@pytest.fixture(scope="module", params=case_dirs("bndf_global_fds611", "bndf_global"), ids=case_id) +def sim(request): + return load(request.param) + + +def _patch(obst, mesh_id, quantity, orientation): + sub = next(sub for sub in obst._all_subobstructions if sub.mesh.id == mesh_id) + return sub.get_data(quantity).data[orientation].data + + +@pytest.mark.parametrize( + "orientation, nodes, cells", + [(-1, (7, 4), (6, 3)), (1, (7, 4), (6, 3)), (-2, (2, 4), (1, 3)), (2, (2, 4), (1, 3)), (3, (2, 7), (1, 6))], +) +def test_time_comes_first_on_faces_one_cell_wide(sim, orientation, nodes, cells): + strip = sim.obstructions.get_by_id("Strip") + n_t = strip.n_t(count_duplicates=False) + assert strip.get_global_boundary_data_arrays(NODES)[orientation].shape == (n_t,) + nodes + assert strip.get_global_boundary_data_arrays(CELLS)[orientation].shape == (n_t,) + cells + + +@pytest.mark.parametrize("quantity", [NODES, CELLS]) +def test_a_face_in_one_mesh_is_its_patch(sim, quantity): + strip = sim.obstructions.get_by_id("Strip") + for orientation, grid in strip.get_global_boundary_data_arrays(quantity).items(): + np.testing.assert_array_equal(grid, _patch(strip, "Fine", quantity, orientation)) + + +@pytest.mark.parametrize("quantity", [NODES, CELLS]) +def test_a_face_at_two_heights_keeps_both(sim, quantity): + table = sim.obstructions.get_by_id("Table") + lower, upper = table.get_global_boundary_data_arrays(quantity)[3] + np.testing.assert_array_equal(lower, _patch(table, "Fine", quantity, 3)) + np.testing.assert_array_equal(upper, _patch(table, "Coarse", quantity, 3)) + + +def test_coarse_values_are_repeated_on_the_finest_grid(sim): + # The face at y = 0.2 spans x 0.6..1.4 and z 0..0.5 in 'Fine', z 0..0.6 in 'Coarse': nodes every 0.1 m + table = sim.obstructions.get_by_id("Table") + grid = table.get_global_boundary_data_arrays(NODES)[-2] + fine = _patch(table, "Fine", NODES, -2) + coarse = _patch(table, "Coarse", NODES, -2) + assert grid.shape == (table.n_t(count_duplicates=False), 9, 7) + np.testing.assert_array_equal(grid[:, :4, :6], fine[:, :4, :]) + np.testing.assert_array_equal(grid[:, 5:7, 2], coarse[:, :2, 1]) # x 1.1 and 1.2 take the nodes 1.0 and 1.2 + np.testing.assert_array_equal(grid[:, 6, :], coarse[:, 1, :].repeat(2, axis=1)[:, :7]) # x 1.2, z 0.0..0.6 + assert np.isnan(grid[:, :4, 6]).all() # 'Fine' ends at z = 0.5 + + grid = table.get_global_boundary_data_arrays(CELLS)[-2] + coarse = _patch(table, "Coarse", CELLS, -2) + assert grid.shape == (table.n_t(count_duplicates=False), 8, 6) + np.testing.assert_array_equal(grid[:, 4:8, :], coarse.repeat(2, axis=1).repeat(2, axis=2)[:, :, :6]) + + +def test_node_based_and_cell_centred_quantities_on_one_obstruction(sim): + table = sim.obstructions.get_by_id("Table") + for quantity in (NODES, CELLS): + assert list(table.get_global_boundary_data_arrays(quantity)) == [-2, -1, 1, 2, 3] # no face on the floor + + +def test_unknown_quantity(sim): + with pytest.raises(KeyError): + sim.obstructions.get_by_id("Table").get_global_boundary_data_arrays("VELOCITY") + + +@pytest.mark.filterwarnings("ignore:invalid value encountered in cast:RuntimeWarning") +def test_export_of_a_face_at_two_heights(sim, tmp_path): + from fdsreader.export import export_obst_raw + + meta = export_obst_raw(sim.obstructions.get_by_id("Table"), str(tmp_path)) + assert (tmp_path / "WALL_TEMPERATURE" / "obst-Table_quantity-WALL_TEMPERATURE.dat").stat().st_size > 0 + assert meta.endswith("obst-Table.yaml")