From ec68bfc095091802f8774f27fe244363f85479b7 Mon Sep 17 00:00:00 2001 From: Hesam Salehipour Date: Fri, 18 Oct 2024 16:27:39 -0400 Subject: [PATCH] moved np.unique to examples and a helper function --- examples/cfd/flow_past_sphere_3d.py | 23 ++++++++----------- examples/cfd/lid_driven_cavity_2d.py | 17 ++++++++------ examples/cfd/turbulent_channel_3d.py | 4 ++-- examples/cfd/windtunnel_3d.py | 22 ++++++++---------- examples/performance/mlups_3d.py | 13 ++++------- xlb/grid/grid.py | 1 - xlb/helper/__init__.py | 1 + xlb/helper/check_boundary_overlaps.py | 22 ++++++++++++++++++ .../indices_boundary_masker.py | 9 ++------ 9 files changed, 61 insertions(+), 51 deletions(-) create mode 100644 xlb/helper/check_boundary_overlaps.py diff --git a/examples/cfd/flow_past_sphere_3d.py b/examples/cfd/flow_past_sphere_3d.py index 994d18b..8b958fc 100644 --- a/examples/cfd/flow_past_sphere_3d.py +++ b/examples/cfd/flow_past_sphere_3d.py @@ -1,7 +1,7 @@ import xlb from xlb.compute_backend import ComputeBackend from xlb.precision_policy import PrecisionPolicy -from xlb.helper import create_nse_fields, initialize_eq +from xlb.helper import create_nse_fields, initialize_eq, check_bc_overlaps from xlb.operator.stepper import IncompressibleNavierStokesStepper from xlb.operator.boundary_condition import ( FullwayBounceBackBC, @@ -48,15 +48,12 @@ def _setup(self, omega): self.setup_stepper(omega) def define_boundary_indices(self): - inlet = self.grid.boundingBoxIndices["left"] - outlet = self.grid.boundingBoxIndices["right"] - walls = [ - self.grid.boundingBoxIndices["bottom"][i] - + self.grid.boundingBoxIndices["top"][i] - + self.grid.boundingBoxIndices["front"][i] - + self.grid.boundingBoxIndices["back"][i] - for i in range(self.velocity_set.d) - ] + box = self.grid.bounding_box_indices() + box_noedge = self.grid.bounding_box_indices(remove_edges=True) + inlet = box_noedge["left"] + outlet = box_noedge["right"] + walls = [box["bottom"][i] + box["top"][i] + box["front"][i] + box["back"][i] for i in range(self.velocity_set.d)] + walls = np.unique(np.array(walls), axis=-1).tolist() sphere_radius = self.grid_shape[1] // 12 x = np.arange(self.grid_shape[0]) @@ -79,12 +76,12 @@ def setup_boundary_conditions(self): # bc_outlet = DoNothingBC(indices=outlet) bc_outlet = ExtrapolationOutflowBC(indices=outlet) bc_sphere = HalfwayBounceBackBC(indices=sphere) - self.boundary_conditions = [bc_walls, bc_left, bc_outlet, bc_sphere] - # Note: it is important to add bc_walls before bc_outlet/bc_inlet because - # of the corner nodes. This way the corners are treated as wall and not inlet/outlet. def setup_boundary_masker(self): + # check boundary condition list for duplicate indices before creating bc mask + check_bc_overlaps(self.boundary_conditions, self.velocity_set.d, self.backend) + indices_boundary_masker = IndicesBoundaryMasker( velocity_set=self.velocity_set, precision_policy=self.precision_policy, diff --git a/examples/cfd/lid_driven_cavity_2d.py b/examples/cfd/lid_driven_cavity_2d.py index 300614d..a77cc62 100644 --- a/examples/cfd/lid_driven_cavity_2d.py +++ b/examples/cfd/lid_driven_cavity_2d.py @@ -1,15 +1,16 @@ import xlb from xlb.compute_backend import ComputeBackend from xlb.precision_policy import PrecisionPolicy -from xlb.helper import create_nse_fields, initialize_eq +from xlb.helper import create_nse_fields, initialize_eq, check_bc_overlaps from xlb.operator.boundary_masker import IndicesBoundaryMasker from xlb.operator.stepper import IncompressibleNavierStokesStepper from xlb.operator.boundary_condition import HalfwayBounceBackBC, EquilibriumBC from xlb.operator.macroscopic import Macroscopic from xlb.utils import save_fields_vtk, save_image +import xlb.velocity_set import warp as wp import jax.numpy as jnp -import xlb.velocity_set +import numpy as np class LidDrivenCavity2D: @@ -39,11 +40,11 @@ def _setup(self, omega): self.setup_stepper(omega) def define_boundary_indices(self): - lid = self.grid.boundingBoxIndices["top"] - walls = [ - self.grid.boundingBoxIndices["bottom"][i] + self.grid.boundingBoxIndices["left"][i] + self.grid.boundingBoxIndices["right"][i] - for i in range(self.velocity_set.d) - ] + box = self.grid.bounding_box_indices() + box_noedge = self.grid.bounding_box_indices(remove_edges=True) + lid = box_noedge["top"] + walls = [box["bottom"][i] + box["left"][i] + box["right"][i] for i in range(self.velocity_set.d)] + walls = np.unique(np.array(walls), axis=-1).tolist() return lid, walls def setup_boundary_conditions(self): @@ -53,6 +54,8 @@ def setup_boundary_conditions(self): self.boundary_conditions = [bc_walls, bc_top] def setup_boundary_masker(self): + # check boundary condition list for duplicate indices before creating bc mask + check_bc_overlaps(self.boundary_conditions, self.velocity_set.d, self.backend) indices_boundary_masker = IndicesBoundaryMasker( velocity_set=self.velocity_set, precision_policy=self.precision_policy, diff --git a/examples/cfd/turbulent_channel_3d.py b/examples/cfd/turbulent_channel_3d.py index 2ec5560..eb73fdc 100644 --- a/examples/cfd/turbulent_channel_3d.py +++ b/examples/cfd/turbulent_channel_3d.py @@ -77,8 +77,8 @@ def _setup(self): def define_boundary_indices(self): # top and bottom sides of the channel are no-slip and the other directions are periodic - boundingBoxIndices = self.grid.bounding_box_indices(remove_edges=True) - walls = [boundingBoxIndices["bottom"][i] + boundingBoxIndices["top"][i] for i in range(self.velocity_set.d)] + box = self.grid.bounding_box_indices(remove_edges=True) + walls = [box["bottom"][i] + box["top"][i] for i in range(self.velocity_set.d)] return walls def setup_boundary_conditions(self): diff --git a/examples/cfd/windtunnel_3d.py b/examples/cfd/windtunnel_3d.py index b0ee5b9..a7567d0 100644 --- a/examples/cfd/windtunnel_3d.py +++ b/examples/cfd/windtunnel_3d.py @@ -3,7 +3,7 @@ import time from xlb.compute_backend import ComputeBackend from xlb.precision_policy import PrecisionPolicy -from xlb.helper import create_nse_fields, initialize_eq +from xlb.helper import create_nse_fields, initialize_eq, check_bc_overlaps from xlb.operator.stepper import IncompressibleNavierStokesStepper from xlb.operator.boundary_condition import ( FullwayBounceBackBC, @@ -67,15 +67,12 @@ def voxelize_stl(self, stl_filename, length_lbm_unit): return mesh_matrix, pitch def define_boundary_indices(self): - inlet = self.grid.boundingBoxIndices["left"] - outlet = self.grid.boundingBoxIndices["right"] - walls = [ - self.grid.boundingBoxIndices["bottom"][i] - + self.grid.boundingBoxIndices["top"][i] - + self.grid.boundingBoxIndices["front"][i] - + self.grid.boundingBoxIndices["back"][i] - for i in range(self.velocity_set.d) - ] + box = self.grid.bounding_box_indices() + box_noedge = self.grid.bounding_box_indices(remove_edges=True) + inlet = box_noedge["left"] + outlet = box_noedge["right"] + walls = [box["bottom"][i] + box["top"][i] + box["front"][i] + box["back"][i] for i in range(self.velocity_set.d)] + walls = np.unique(np.array(walls), axis=-1).tolist() # Load the mesh stl_filename = "examples/cfd/stl-files/DrivAer-Notchback.stl" @@ -105,10 +102,11 @@ def setup_boundary_conditions(self): bc_car = GradsApproximationBC(mesh_vertices=car) # bc_car = FullwayBounceBackBC(mesh_vertices=car) self.boundary_conditions = [bc_walls, bc_left, bc_do_nothing, bc_car] - # Note: it is important to add bc_walls before bc_outlet/bc_inlet because - # of the corner nodes. This way the corners are treated as wall and not inlet/outlet. def setup_boundary_masker(self): + # check boundary condition list for duplicate indices before creating bc mask + check_bc_overlaps(self.boundary_conditions, self.velocity_set.d, self.backend) + indices_boundary_masker = IndicesBoundaryMasker( velocity_set=self.velocity_set, precision_policy=self.precision_policy, diff --git a/examples/performance/mlups_3d.py b/examples/performance/mlups_3d.py index 1812d95..4c337e8 100644 --- a/examples/performance/mlups_3d.py +++ b/examples/performance/mlups_3d.py @@ -48,15 +48,10 @@ def create_grid_and_fields(cube_edge): def define_boundary_indices(grid): - lid = grid.boundingBoxIndices["top"] - walls = [ - grid.boundingBoxIndices["bottom"][i] - + grid.boundingBoxIndices["left"][i] - + grid.boundingBoxIndices["right"][i] - + grid.boundingBoxIndices["front"][i] - + grid.boundingBoxIndices["back"][i] - for i in range(len(grid.shape)) - ] + box = grid.bounding_box_indices() + box_noedge = grid.bounding_box_indices(remove_edges=True) + lid = box_noedge["top"] + walls = [box["bottom"][i] + box["left"][i] + box["right"][i] + box["front"][i] + box["back"][i] for i in range(len(grid.shape))] return lid, walls diff --git a/xlb/grid/grid.py b/xlb/grid/grid.py index 7494c3e..53139fc 100644 --- a/xlb/grid/grid.py +++ b/xlb/grid/grid.py @@ -25,7 +25,6 @@ def __init__(self, shape: Tuple[int, ...], compute_backend: ComputeBackend): self.shape = shape self.dim = len(shape) self.compute_backend = compute_backend - self.boundingBoxIndices = self.bounding_box_indices() self._initialize_backend() @abstractmethod diff --git a/xlb/helper/__init__.py b/xlb/helper/__init__.py index 92d3583..4c63aa6 100644 --- a/xlb/helper/__init__.py +++ b/xlb/helper/__init__.py @@ -1,2 +1,3 @@ from xlb.helper.nse_solver import create_nse_fields as create_nse_fields from xlb.helper.initializers import initialize_eq as initialize_eq +from xlb.helper.check_boundary_overlaps import check_bc_overlaps as check_bc_overlaps diff --git a/xlb/helper/check_boundary_overlaps.py b/xlb/helper/check_boundary_overlaps.py new file mode 100644 index 0000000..18e17a1 --- /dev/null +++ b/xlb/helper/check_boundary_overlaps.py @@ -0,0 +1,22 @@ +import numpy as np +from xlb.compute_backend import ComputeBackend + + +def check_bc_overlaps(bclist, dim, backend): + index_list = [[] for _ in range(dim)] + for bc in bclist: + if bc.indices is None: + continue + # Detect duplicates within bc.indices + index_arr = np.unique(bc.indices, axis=-1) + if index_arr.shape[-1] != len(bc.indices[0]): + if backend == ComputeBackend.WARP: + raise ValueError(f"Boundary condition {bc.__class__.__name__} has duplicate indices!") + for d in range(dim): + index_list[d] += bc.indices[d] + + # Detect duplicates within bclist + index_arr = np.unique(index_list, axis=-1) + if index_arr.shape[-1] != len(index_list[0]): + if backend == ComputeBackend.WARP: + raise ValueError("Boundary condition list containes duplicate indices!") diff --git a/xlb/operator/boundary_masker/indices_boundary_masker.py b/xlb/operator/boundary_masker/indices_boundary_masker.py index 10d72c4..848d74e 100644 --- a/xlb/operator/boundary_masker/indices_boundary_masker.py +++ b/xlb/operator/boundary_masker/indices_boundary_masker.py @@ -225,14 +225,9 @@ def warp_implementation(self, bclist, bc_mask, missing_mask, start_index=None): # We are done with bc.indices. Remove them from BC objects bc.__dict__.pop("indices", None) - # Remove duplicates indices to avoid race conditioning - index_arr, unique_loc = np.unique(index_list, axis=-1, return_index=True) - id_arr = np.array(id_list)[unique_loc] - is_interior = np.array(is_interior)[unique_loc] - # convert to warp arrays - indices = wp.array2d(index_arr, dtype=wp.int32) - id_number = wp.array1d(id_arr, dtype=wp.uint8) + indices = wp.array2d(index_list, dtype=wp.int32) + id_number = wp.array1d(id_list, dtype=wp.uint8) is_interior = wp.array1d(is_interior, dtype=wp.bool) if start_index is None: