270 efficient hwbb and ibb

examples/advanced_projects/efficient_bounce_back_obstacle/__init__.py

@@ -0,0 +1,11 @@
+from .solid_boundary_data import SolidBoundaryData
+from .fullway_bounce_back_boundary import FullwayBounceBackBoundary
+from .halfway_bounce_back_boundary import HalfwayBounceBackBoundary
+from .linear_interpolated_bounce_back_boundary import LinearInterpolatedBounceBackBoundary
+
+__all__ = [
+    'FullwayBounceBackBoundary',
+    'HalfwayBounceBackBoundary',
+    'LinearInterpolatedBounceBackBoundary',
+    'SolidBoundaryData'
+]

examples/advanced_projects/efficient_bounce_back_obstacle/ebb_simulation.py

@@ -0,0 +1,106 @@
+import torch
+
+from timeit import default_timer as timer
+from typing import List
+
+from lettuce._simulation import *
+from lettuce import Flow
+
+__all__ = ['EbbSimulation']
+
+class EbbSimulation(Simulation):
+    """
+        advanced simulation class, with
+        - (additional) post_streaming_boundaries (substeps: C -> S -> B -> R)
+        - post_collision population information for post_streaming_boundaries
+    """
+    def __init__(self, flow: 'Flow', collision: 'Collision',
+                 reporter: List['Reporter']):
+        flow.context.use_native = False
+        super().__init__(flow, collision, reporter)
+
+        if hasattr(flow.__class__, 'post_streaming_boundaries'):
+            # list of boundaries that are applied AFTER the streaming step
+            self.post_streaming_boundaries = flow.post_streaming_boundaries
+        else:
+            self.post_streaming_boundaries = []
+
+        # adjust No_collision_ and no_streaming_mask for use of post_streaming_boundaries (!)
+        if len(self.post_streaming_boundaries) > 0:
+
+            # create NCM and NSM,
+            # ...if there were no pre_ or post_boundaries that triggered
+            # ...their creation in super-class
+            if self.no_collision_mask is None:
+                # fill masks with value of self.collision_index (= number of pre-boundaries)
+                self.no_collision_mask = self.context.full_tensor(
+                    flow.resolution, self.collision_index, dtype=torch.uint8)
+            if self.no_streaming_mask is None:
+                self.no_streaming_mask = self.context.full_tensor([flow.stencil.q,
+                                                                   *flow.resolution],
+                                                                  self.collision_index,
+                                                                  dtype=torch.uint8)
+
+
+            for i, boundary in enumerate(self.post_streaming_boundaries,
+                                         start=self.collision_index + 1
+                                               + len(self.post_boundaries)):
+# FOR DEBUGGING: print("SIMULATION: creating no_collision_mask entries for: i, boundary = ", i, boundary)
+# FOR DEBUGGING: print("collision index is:", self.collision_index)
+                ncm = boundary.make_no_collision_mask(
+                    [it for it in self.flow.f.shape[1:]], context=self.context)
+                if ncm is not None:
+                    self.no_collision_mask[ncm] = i
+                nsm = boundary.make_no_streaming_mask(
+                    [it for it in self.flow.f.shape], context=self.context)
+                if nsm is not None:
+                    self.no_streaming_mask |= nsm
+
+        # get indices of post_streaming_boundaries, that need f_collided to be stored
+        self.store_f_collided_post_streaming_boundaries_index = []
+        for i, boundary in enumerate(self.post_streaming_boundaries):
+            if hasattr(boundary.__class__, "store_f_collided"):
+                # append boundary to lis of boundaries that need
+                # ...f_collided state between collision and streaming substep
+                self.store_f_collided_post_streaming_boundaries_index.append(i)
+            if hasattr(boundary.__class__, "initialize_f_collided"):
+                # initialize the f_collided storage in each of those boundaries
+                boundary.initialize_f_collided()
+
+        # redefine collide_and_stream() to include the storage of f_collided
+        # ...for use in post_streaming_boundaries
+        def collide_and_stream(*_, **__):
+            self._collide()
+            # pass f to post_streaming_boundaries between collision and streaming substep
+            self._pass_f_collided()
+            self._stream()
+
+        self._collide_and_stream = collide_and_stream
+
+
+    def _post_streaming_boundaries(self):
+        # runs all the post_streaming_boundaries; required for efficient BBBC
+
+        for boundary in self.post_streaming_boundaries:
+            boundary(self.flow)
+
+    def _pass_f_collided(self):
+        # passes f to post_streaming_boundaries, for later use
+
+        for idx in self.store_f_collided_post_streaming_boundaries_index:
+            self.post_streaming_boundaries[idx].store_f_collided(self.flow.f)
+
+    def __call__(self, num_steps):
+        beg = timer()
+
+        if self.flow.i == 0:
+            self._report()
+
+        for _ in range(num_steps):
+            self._collide_and_stream(self)
+            self._post_streaming_boundaries()
+            self.flow.i += 1
+            self._report()
+
+        end = timer()
+        return num_steps * self.flow.rho().numel() / 1e6 / (end - beg)

examples/advanced_projects/efficient_bounce_back_obstacle/fullway_bounce_back_boundary.py

@@ -0,0 +1,185 @@
+import torch
+import numpy as np
+
+from typing import List, Optional
+from lettuce import Boundary, Flow, Context
+
+__all__ = ["FullwayBounceBackBoundary"]
+
+class FullwayBounceBackBoundary(Boundary):
+    """
+        FULLWAY BOUNCE BACK BOUNDARY CONDITION (FWBB):
+        - inverts populations on solid boundary nodes,
+            instead of colliding them
+        - (!) mind observable-artifacts on solid regions. They do not affect
+            the flow in the fluid region, but can mess with visualization and
+            non-local statistics (mean velocity etc.)
+        - able to calculate fluid force on boundary by momentum exchange
+    """
+
+    def __init__(self, context, flow, mask, global_solid_mask=None,
+                 periodicity: tuple[bool,bool,bool|None] = None, calc_force=False):
+        self.context = context
+        self.flow = flow
+
+        self.mask = mask
+
+        if periodicity is None:
+            periodicity = (False, False, False if self.flow.stencil.d == 3 else None)
+
+        # global_solid_mask (optional) to filter out all "fake" fluid neighbors,
+        # ...which are outside the FWBB but not in the fluid region
+        if global_solid_mask is None:
+            global_solid_mask = np.zeros_like(self.mask, dtype=bool)
+
+        # exclude self.mask from global_solid_mask for "other" solid mask
+        other_solid_bc_mask = np.where(~self.mask, global_solid_mask, False)
+
+        if calc_force:
+            # summed force vector on all boundary nodes, in D dimensions (x,y,(z))
+            self.force_sum = torch.zeros_like(self.context.convert_to_tensor(
+                self.flow.torch_stencil.e[0]))
+            self.calc_force = True
+        else:
+            self.calc_force = False
+
+        ### create f_index_fwbb, needed for force-calculation
+        # ...(marks all fs which streamed into the boundary in prior streaming step)
+        # ... in other words: marks all fs that need to be bounced
+        self.f_index_fwbb = []
+
+        if self.flow.stencil.d == 2:
+            nx, ny = mask.shape  # domain size in x and y
+            ix_sp, iy_sp = np.where(mask)
+                # np.arrays: list of (ix_sp) x-indizes and (iy_sp) y-indizes in the boundary.mask
+                # ...to enable iteration over all boundary/wall/object-nodes
+            for sp_index in range(0, len(ix_sp)):  # for all TRUE-nodes in boundary.mask
+                for q_index in range(0, self.flow.stencil.q):
+                    # for all stencil-directions c_i (lattice.stencil.e in lettuce)
+
+                    # check for boundary-nodes neighboring the domain-border.
+                    # ...they have to take the periodicity into account...
+                    border = np.zeros(self.flow.stencil.d, dtype=int)
+                    if (ix_sp[sp_index] == 0 and self.flow.stencil.e[q_index][ 0] == -1
+                            and periodicity[0]):  # searching border on left
+                        border[0] = -1
+                    elif (ix_sp[sp_index] == nx - 1 and self.flow.stencil.e[q_index][ 0] == 1
+                          and periodicity[0]):  # searching border on right
+                        border[0] = 1
+                    if (iy_sp[sp_index] == 0 and self.flow.stencil.e[q_index][ 1] == -1
+                            and periodicity[1]):  # searching border on left
+                        border[1] = -1
+                    elif (iy_sp[sp_index] == ny - 1 and self.flow.stencil.e[q_index][ 1] == 1
+                          and periodicity[1]):  # searching border on right
+                        border[1] = 1
+                    try:  # try in case the neighboring cell does not exist
+                            # (= an f pointing out of the simulation domain)
+                        if (not mask[ix_sp[sp_index] + self.flow.stencil.e[q_index][ 0] - border[0]*nx,
+                                    iy_sp[sp_index] + self.flow.stencil.e[q_index][ 1] - border[1]*ny]
+                                and not other_solid_bc_mask[ix_sp[sp_index] + self.flow.stencil.e[q_index][ 0] - border[0]*nx,
+                                    iy_sp[sp_index] + self.flow.stencil.e[q_index][ 1] - border[1]*ny]):
+                            # if the neighbour of sp_index is False in the boundary.mask,
+                            # ...sp_index is ix_sp solid node, neighbouring ix_sp fluid node:
+                            # the direction pointing from the fluid neighbour to solid
+                            # ...sp_index is marked on the solid sp_index
+
+                            # list of [q, nx, ny], marks all fs on the boundary-border, which point into the boundary/solid
+                            self.f_index_fwbb.append([self.flow.stencil.opposite[q_index], ix_sp[sp_index], iy_sp[sp_index]])
+                    except IndexError:
+                        pass  # just ignore this iteration since there is no neighbor there
+        if self.flow.stencil.d == 3:  # like 2D, but in 3D...guess what...
+            nx, ny, nz = mask.shape
+            ix_sp, iy_sp, iz_sp = np.where(mask)
+            for sp_index in range(0, len(ix_sp)):
+                for q_index in range(0, self.flow.stencil.q):
+                    border = np.zeros(self.flow.stencil.d, dtype=int)
+                    if (ix_sp[sp_index] == 0 and self.flow.stencil.e[q_index][ 0] == -1
+                            and periodicity[0]):  # searching border on left
+                        border[0] = -1
+                    elif (ix_sp[sp_index] == nx - 1 and self.flow.stencil.e[q_index][ 0] == 1
+                          and periodicity[0]):  # searching border on right
+                        border[0] = 1
+                    if (iy_sp[sp_index] == 0 and self.flow.stencil.e[q_index][ 1] == -1
+                            and periodicity[1]):  # searching border on left
+                        border[1] = -1
+                    elif (iy_sp[sp_index] == ny - 1 and self.flow.stencil.e[q_index][ 1] == 1
+                          and periodicity[1]):  # searching border on right
+                        border[1] = 1
+                    if (iz_sp[sp_index] == 0 and self.flow.stencil.e[q_index][ 2] == -1
+                            and periodicity[2]):  # searching border on left
+                        border[2] = -1
+                    elif (iz_sp[sp_index] == nz - 1 and self.flow.stencil.e[q_index][ 2] == 1
+                          and periodicity[2]):  # searching border on right
+                        border[2] = 1
+                    try:  # try in case the neighboring cell does not exist
+                            # (and f pointing out of simulation domain)
+                        if (not mask[ix_sp[sp_index] + self.flow.stencil.e[q_index][ 0] - border[0] * nx,
+                                    iy_sp[sp_index] + self.flow.stencil.e[q_index][ 1] - border[1] * ny,
+                                    iz_sp[sp_index] + self.flow.stencil.e[q_index][ 2] - border[2] * nz]
+                                and not other_solid_bc_mask[ix_sp[sp_index] + self.flow.stencil.e[q_index][ 0] - border[0] * nx,
+                                                            iy_sp[sp_index] + self.flow.stencil.e[q_index][ 1] - border[1] * ny,
+                                                            iz_sp[sp_index] + self.flow.stencil.e[q_index][ 2] - border[2] * nz]):
+                            self.f_index_fwbb.append([self.flow.stencil.opposite[q_index], ix_sp[sp_index], iy_sp[sp_index], iz_sp[sp_index]])
+                    except IndexError:
+                        pass  # just ignore this iteration since there is no neighbor there
+
+        self.f_index_fwbb = torch.tensor(np.array(self.f_index_fwbb),
+                                         device=flow.context.device,
+                                         dtype=torch.int64) # the batch-index has to be integer
+        self.opposite_tensor = torch.tensor(self.flow.stencil.opposite,
+                                            device=flow.context.device,
+                                            dtype=torch.int64)  # batch-index has to be ix_sp tensor
+
+
+    def __call__(self, flow):
+        # FULLWAY-BBBC: inverts populations on all boundary nodes
+
+        # calc force on boundary:#
+        if self.calc_force:
+            self.calc_force_on_boundary(flow.f)
+
+        # bounce (invert populations on boundary nodes)
+        # -> basically an efficient version of:
+        #    f = torch.where(self.mask, f[self.flow.stencil.opposite], f)
+
+        if self.flow.torch_stencil.d == 2:
+            flow.f[self.opposite_tensor[self.f_index_fwbb[:, 0]],
+            self.f_index_fwbb[:, 1],
+            self.f_index_fwbb[:, 2]] = flow.f[self.f_index_fwbb[:, 0],
+            self.f_index_fwbb[:, 1],
+            self.f_index_fwbb[:, 2]]
+        if self.flow.torch_stencil.d == 3:
+            flow.f[self.opposite_tensor[self.f_index_fwbb[:, 0]],
+            self.f_index_fwbb[:, 1],
+            self.f_index_fwbb[:, 2],
+            self.f_index_fwbb[:, 3]] = flow.f[self.f_index_fwbb[:, 0],
+            self.f_index_fwbb[:, 1],
+            self.f_index_fwbb[:, 2],
+            self.f_index_fwbb[:, 3]]
+
+    def calc_force_on_boundary(self, f):
+        """calculate fluid force on all relevant boundary nodes by
+            momentum exchange"""
+        if self.flow.torch_stencil.d == 2:
+            self.force_sum = 2 * torch.einsum('i..., id -> d', f[self.f_index_fwbb[:, 0],
+            self.f_index_fwbb[:, 1],
+            self.f_index_fwbb[:, 2]], self.flow.torch_stencil.e[self.f_index_fwbb[:, 0]])
+        if self.flow.torch_stencil.d == 3:
+            self.force_sum = 2 * torch.einsum('i..., id -> d', f[self.f_index_fwbb[:, 0],
+            self.f_index_fwbb[:, 1],
+            self.f_index_fwbb[:, 2],
+            self.f_index_fwbb[:, 3]], self.flow.torch_stencil.e[self.f_index_fwbb[:, 0]])
+
+    def make_no_collision_mask(self, f_shape: List[int], context: 'Context') -> Optional[torch.Tensor]:
+        return self.context.convert_to_tensor(self.mask, dtype=bool)
+
+    def make_no_streaming_mask(self, shape: List[int], context: 'Context') -> Optional[torch.Tensor]:
+        """FWBB needs streaming to invert the populations on the solid itself!"""
+        pass
+
+    def native_available(self) -> bool:
+        return False
+
+    def native_generator(self, index: int):
+        # not implemented yet
+        pass