diff --git a/lbmpy/boundaries/boundaryconditions.py b/lbmpy/boundaries/boundaryconditions.py index 96251e0e858ddd09c1398ba9e0a02878dd7e6783..29e20bc4c0ff3969da4f4ef47c1678396dcc0019 100644 --- a/lbmpy/boundaries/boundaryconditions.py +++ b/lbmpy/boundaries/boundaryconditions.py @@ -283,7 +283,7 @@ class UBB(LbBoundary): Returns: list containing LbmWeightInfo and NeighbourOffsetArrays """ - return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)] + return [LbmWeightInfo(lb_method, self.data_type), NeighbourOffsetArrays(lb_method.stencil)] @property def velocity_is_callable(self): @@ -312,7 +312,8 @@ class UBB(LbBoundary): velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.first_order_moment_symbols).main_assignments] c_s_sq = sp.Rational(1, 3) - weight_of_direction = LbmWeightInfo.weight_of_direction + weight_info = LbmWeightInfo(lb_method, data_type=self.data_type) + weight_of_direction = weight_info.weight_of_direction vel_term = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction( dir_symbol, lb_method) @@ -595,10 +596,11 @@ class DiffusionDirichlet(LbBoundary): name: optional name of the boundary. """ - def __init__(self, concentration, name=None): + def __init__(self, concentration, name=None, data_type='double'): if name is None: name = "Diffusion Dirichlet " + str(concentration) self.concentration = concentration + self._data_type = data_type super(DiffusionDirichlet, self).__init__(name) @@ -611,10 +613,11 @@ class DiffusionDirichlet(LbBoundary): Returns: list containing LbmWeightInfo """ - return [LbmWeightInfo(lb_method)] + return [LbmWeightInfo(lb_method, self._data_type)] def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): - w_dir = LbmWeightInfo.weight_of_direction(dir_symbol, lb_method) + weight_info = LbmWeightInfo(lb_method, self._data_type) + w_dir = weight_info.weight_of_direction(dir_symbol, lb_method) return [Assignment(f_in(inv_dir[dir_symbol]), 2 * w_dir * self.concentration - f_out(dir_symbol))] diff --git a/lbmpy/boundaries/boundaryhandling.py b/lbmpy/boundaries/boundaryhandling.py index 71bded3b415ac7874d110ec568173337cd9693cd..fc1fae4a5e04d09ebdb0a671f00a1294e713a839 100644 --- a/lbmpy/boundaries/boundaryhandling.py +++ b/lbmpy/boundaries/boundaryhandling.py @@ -38,7 +38,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling): self._prev_timestep = None def add_fixed_steps(self, fixed_loop, **kwargs): - if self._inplace: # Fixed Loop can't do timestep selection + if self._inplace: # Fixed Loop can't do timestep selection raise NotImplementedError("Adding to fixed loop is currently not supported for inplace kernels") super(LatticeBoltzmannBoundaryHandling, self).add_fixed_steps(fixed_loop, **kwargs) @@ -52,10 +52,12 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling): if boundary_obj not in self._boundary_object_to_boundary_info: sym_index_field = Field.create_generic('indexField', spatial_dimensions=1, dtype=numpy_data_type_for_boundary_object(boundary_obj, self.dim)) - kernels = [self._create_boundary_kernel( - self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.EVEN).compile(), - self._create_boundary_kernel( - self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.ODD).compile()] + + ast_even = self._create_boundary_kernel(self._data_handling.fields[self._field_name], sym_index_field, + boundary_obj, Timestep.EVEN) + ast_odd = self._create_boundary_kernel(self._data_handling.fields[self._field_name], sym_index_field, + boundary_obj, Timestep.ODD) + kernels = [ast_even.compile(), ast_odd.compile()] if flag is None: flag = self.flag_interface.reserve_next_flag() boundary_info = self.InplaceStreamingBoundaryInfo(self, boundary_obj, flag, kernels) @@ -84,6 +86,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling): self.boundary_object = boundary_obj self.flag = flag self._kernels = kernels + # end class InplaceStreamingBoundaryInfo # ------------------------------ Force On Boundary ------------------------------------------------------------ @@ -148,29 +151,32 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling): return dh.reduce_float_sequence(list(result), 'sum') + # end class LatticeBoltzmannBoundaryHandling class LbmWeightInfo(CustomCodeNode): + def __init__(self, lb_method, data_type='double'): + self.weights_symbol = TypedSymbol("weights", data_type) + data_type_string = "double" if self.weights_symbol.dtype.numpy_dtype == np.float64 else "float" - # --------------------------- Functions to be used by boundaries -------------------------- + weights = [str(w.evalf()) for w in lb_method.weights] + if data_type_string == "float": + weights = "f, ".join(weights) + weights += "f" # suffix for the last element + else: + weights = ", ".join(weights) + w_sym = self.weights_symbol + code = f"const {data_type_string} {w_sym.name} [] = {{{ weights }}};\n" + super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym}) - @staticmethod - def weight_of_direction(dir_idx, lb_method=None): + def weight_of_direction(self, dir_idx, lb_method=None): if isinstance(sp.sympify(dir_idx), sp.Integer): return lb_method.weights[dir_idx].evalf() else: - return sp.IndexedBase(LbmWeightInfo.WEIGHTS_SYMBOL, shape=(1,))[dir_idx] - - # ---------------------------------- Internal --------------------------------------------- + return sp.IndexedBase(self.weights_symbol, shape=(1,))[dir_idx] - WEIGHTS_SYMBOL = TypedSymbol("weights", "double") - def __init__(self, lb_method): - weights = [str(w.evalf()) for w in lb_method.weights] - w_sym = LbmWeightInfo.WEIGHTS_SYMBOL - code = "const double %s [] = { %s };\n" % (w_sym.name, ",".join(weights)) - super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym}) # end class LbmWeightInfo