diff --git a/lbmpy/boundaries/boundaryconditions.py b/lbmpy/boundaries/boundaryconditions.py index e0ba805da2da8f1161e25e7b86e4b66721cb81f8..b82b745e796815d6aa1038a18a34f92c6fee64a7 100644 --- a/lbmpy/boundaries/boundaryconditions.py +++ b/lbmpy/boundaries/boundaryconditions.py @@ -26,8 +26,9 @@ class LbBoundary(abc.ABC): def __init__(self, name=None): self._name = name + self._data_type = None - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): """ This function defines the boundary behavior and must therefore be implemented by all boundaries. The boundary is defined through a list of sympy equations from which a boundary kernel is generated. @@ -107,7 +108,7 @@ class NoSlip(LbBoundary): """Set an optional name here, to mark boundaries, for example for force evaluations""" super(NoSlip, self).__init__(name) - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol)) @@ -212,7 +213,7 @@ class FreeSlip(LbBoundary): else: return [] - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): if self.normal_direction: normal_direction = self.normal_direction mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis) @@ -245,7 +246,7 @@ class UBB(LbBoundary): name: optional name of the boundary. """ - def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'): + def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type=None): self._velocity = velocity self._adaptVelocityToForce = adapt_velocity_to_force if callable(self._velocity) and not dim: @@ -253,7 +254,6 @@ class UBB(LbBoundary): elif not callable(self._velocity): dim = len(velocity) self.dim = dim - self.data_type = data_type super(UBB, self).__init__(name) @@ -262,7 +262,7 @@ class UBB(LbBoundary): """ In case of the UBB boundary additional data is a velocity vector. This vector is added to each cell to realize velocity profiles for the inlet.""" if self.velocity_is_callable: - return [(f'vel_{i}', create_type(self.data_type)) for i in range(self.dim)] + return [(f'vel_{i}', create_type(self._data_type)) for i in range(self.dim)] else: return [] @@ -283,7 +283,7 @@ class UBB(LbBoundary): Returns: list containing LbmWeightInfo and NeighbourOffsetArrays """ - return [LbmWeightInfo(lb_method, data_type=self.data_type), NeighbourOffsetArrays(lb_method.stencil)] + return [LbmWeightInfo(lb_method, data_type=self._data_type), NeighbourOffsetArrays(lb_method.stencil)] @property def velocity_is_callable(self): @@ -291,7 +291,11 @@ class UBB(LbBoundary): This is useful if the inflow velocity should have a certain profile for instance""" return callable(self._velocity) - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): + if data_type is None: + self._data_type = "double" + else: + self._data_type = data_type vel_from_idx_field = callable(self._velocity) vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity @@ -312,7 +316,7 @@ 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_info = LbmWeightInfo(lb_method, data_type=self.data_type) + 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) @@ -376,7 +380,7 @@ class SimpleExtrapolationOutflow(LbBoundary): """ return [NeighbourOffsetArrays(lb_method.stencil)] - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil) tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction)) @@ -416,7 +420,7 @@ class ExtrapolationOutflow(LbBoundary): def __init__(self, normal_direction, lb_method, dt=1, dx=1, name=None, streaming_pattern='pull', zeroth_timestep=Timestep.BOTH, - initial_density=None, initial_velocity=None, data_type='double'): + initial_density=None, initial_velocity=None, data_type=None): self.lb_method = lb_method self.stencil = lb_method.stencil @@ -439,8 +443,6 @@ class ExtrapolationOutflow(LbBoundary): self.initial_velocity = initial_velocity self.equilibrium_calculation = None - self.data_type = data_type - if initial_density and initial_velocity: equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([])) rho = lb_method.zeroth_order_equilibrium_moment_symbol @@ -487,7 +489,7 @@ class ExtrapolationOutflow(LbBoundary): def additional_data(self): """Used internally only. For the ExtrapolationOutflow information of the previous PDF values is needed. This information is stored in the index vector.""" - data = [('pdf', create_type(self.data_type)), ('pdf_nd', create_type(self.data_type))] + data = [('pdf', create_type(self._data_type)), ('pdf_nd', create_type(self._data_type))] return data @property @@ -506,7 +508,11 @@ class ExtrapolationOutflow(LbBoundary): """ return [NeighbourOffsetArrays(lb_method.stencil)] - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): + if data_type is None: + self._data_type = "double" + else: + self._data_type = data_type subexpressions = [] boundary_assignments = [] dtdx = sp.Rational(self.dt, self.dx) @@ -549,7 +555,7 @@ class FixedDensity(LbBoundary): super(FixedDensity, self).__init__(name) - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom): new_main_assignments = [Assignment(a.lhs, get_symmetric_part(a.rhs, degrees_of_freedom)) for a in assignment_collection.main_assignments] @@ -596,11 +602,10 @@ class DiffusionDirichlet(LbBoundary): name: optional name of the boundary. """ - def __init__(self, concentration, name=None, data_type='double'): + def __init__(self, concentration, name=None): if name is None: name = "Diffusion Dirichlet " + str(concentration) self.concentration = concentration - self.data_type = data_type super(DiffusionDirichlet, self).__init__(name) @@ -613,13 +618,17 @@ class DiffusionDirichlet(LbBoundary): Returns: list containing LbmWeightInfo """ - return [LbmWeightInfo(lb_method, self.data_type)] + return [LbmWeightInfo(lb_method, self._data_type)] - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): - weight_info = LbmWeightInfo(lb_method, data_type=self.data_type) + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): + if data_type is None: + self._data_type = "double" + else: + self._data_type = data_type + weight_info = LbmWeightInfo(lb_method, data_type=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))] + + return [Assignment(f_in(inv_dir[dir_symbol]), 2 * w_dir * self.concentration - f_out(dir_symbol))] # end class DiffusionDirichlet @@ -640,7 +649,7 @@ class NeumannByCopy(LbBoundary): """ return [NeighbourOffsetArrays(lb_method.stencil)] - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil) return [Assignment(f_in(inv_dir[dir_symbol]), f_out(inv_dir[dir_symbol])), Assignment(f_out[neighbour_offset](dir_symbol), f_out(dir_symbol))] @@ -673,7 +682,7 @@ class StreamInConstant(LbBoundary): """ return [NeighbourOffsetArrays(lb_method.stencil)] - def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): + def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, data_type=None): neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil) return [Assignment(f_in(inv_dir[dir_symbol]), self.constant), Assignment(f_out[neighbour_offset](dir_symbol), self.constant)] diff --git a/lbmpy/boundaries/boundaryhandling.py b/lbmpy/boundaries/boundaryhandling.py index 42e0638f3318e56a5c9eee96eaabc9ef777110fb..1ca7b972a76e3b704d2ba3d0ed7f5efea8a96913 100644 --- a/lbmpy/boundaries/boundaryhandling.py +++ b/lbmpy/boundaries/boundaryhandling.py @@ -188,7 +188,8 @@ def create_lattice_boltzmann_boundary_kernel(pdf_field, index_field, lb_method, dir_symbol = indexing.dir_symbol inv_dir = indexing.inverse_dir_symbol - boundary_assignments = boundary_functor(f_out, f_in, dir_symbol, inv_dir, lb_method, index_field) + boundary_assignments = boundary_functor(f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, + data_type=pdf_field.dtype.numpy_dtype) boundary_assignments = indexing.substitute_proxies(boundary_assignments) # Code Elements inside the loop