boundaryconditions.py 23 KB
Newer Older
 Frederik Hennig committed Jan 06, 2021 1 2 from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep from pystencils.simp.assignment_collection import AssignmentCollection  Martin Bauer committed Jul 11, 2019 3 from pystencils import Assignment, Field  Frederik Hennig committed Nov 09, 2020 4 from lbmpy.boundaries.boundaryhandling import LbmWeightInfo  Martin Bauer committed Apr 10, 2018 5 from pystencils.data_types import create_type  Martin Bauer committed Jul 11, 2019 6 from pystencils.sympyextensions import get_symmetric_part  Frederik Hennig committed Nov 09, 2020 7 8 from lbmpy.simplificationfactory import create_simplification_strategy from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays  Frederik Hennig committed Jan 20, 2021 9 from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction  Martin Bauer committed Feb 09, 2017 10   Frederik Hennig committed Feb 02, 2021 11 12 import sympy as sp  Martin Bauer committed Feb 09, 2017 13   Frederik Hennig committed Nov 09, 2020 14 class LbBoundary:  Frederik Hennig committed Jan 06, 2021 15 16 17 18 19  """Base class that all boundaries should derive from. Args: name: optional name of the boundary. """  Martin Bauer committed Feb 09, 2017 20   Martin Bauer committed Jan 30, 2019 21 22 23  inner_or_boundary = True single_link = False  Martin Bauer committed Dec 11, 2017 24 25 26  def __init__(self, name=None): self._name = name  Frederik Hennig committed Nov 09, 2020 27  def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):  Martin Bauer committed May 16, 2017 28 29  """ This function defines the boundary behavior and must therefore be implemented by all boundaries.  Frederik Hennig committed Nov 09, 2020 30 31  The boundary is defined through a list of sympy equations from which a boundary kernel is generated.  Frederik Hennig committed Jan 20, 2021 32 33 34 35 36 37 38 39 40 41 42 43  Args: f_out: a pystencils field acting as a proxy to access the populations streaming out of the current cell, i.e. the post-collision PDFs of the previous LBM step f_in: a pystencils field acting as a proxy to access the populations streaming into the current cell, i.e. the pre-collision PDFs for the next LBM step dir_symbol: a sympy symbol that can be used as an index to f_out and f_in. It describes the direction pointing from the fluid to the boundary cell. inv_dir: an indexed sympy symbol which describes the inversion of a direction index. It can be used in the indices of f_out and f_in for retrieving the PDF of the inverse direction. lb_method: an instance of the LB method used. Use this to adapt the boundary to the method (e.g. compressibility) index_field: the boundary index field that can be used to retrieve and update boundary data  Frederik Hennig committed Nov 09, 2020 44   Frederik Hennig committed Jan 20, 2021 45 46  Returns: list of pystencils assignments, or pystencils.AssignmentCollection  Martin Bauer committed May 16, 2017 47 48 49 50  """ raise NotImplementedError("Boundary class has to overwrite __call__") @property  Martin Bauer committed Apr 10, 2018 51  def additional_data(self):  Martin Bauer committed Oct 21, 2017 52  """Return a list of (name, type) tuples for additional data items required in this boundary  Martin Bauer committed Apr 10, 2018 53 54  These data items can either be initialized in separate kernel see additional_data_kernel_init or by Python callbacks - see additional_data_callback """  Martin Bauer committed May 16, 2017 55 56  return []  Martin Bauer committed Oct 21, 2017 57  @property  Martin Bauer committed Apr 10, 2018 58  def additional_data_init_callback(self):  Martin Bauer committed Feb 06, 2018 59 60  """Return a callback function called with a boundary data setter object and returning a dict of data-name to data for each element that should be initialized"""  Martin Bauer committed Oct 21, 2017 61 62  return None  Frederik Hennig committed Nov 09, 2020 63  def get_additional_code_nodes(self, lb_method):  Frederik Hennig committed Jan 06, 2021 64 65 66 67 68  """Return a list of code nodes that will be added in the generated code before the index field loop. Args: lb_method: lattice Boltzmann method. See :func:lbmpy.creationfunctions.create_lb_method """  Frederik Hennig committed Nov 09, 2020 69 70  return []  Martin Bauer committed May 16, 2017 71 72  @property def name(self):  Martin Bauer committed Dec 11, 2017 73 74 75 76 77 78  if self._name: return self._name else: return type(self).__name__ @name.setter  Martin Bauer committed Apr 10, 2018 79 80  def name(self, new_value): self._name = new_value  Martin Bauer committed May 16, 2017 81   Frederik Hennig committed Jan 20, 2021 82   Frederik Hennig committed Nov 09, 2020 83 84 # end class Boundary  Martin Bauer committed May 16, 2017 85   Frederik Hennig committed Nov 09, 2020 86 87 class NoSlip(LbBoundary): """  Markus Holzer committed May 12, 2021 88 89  No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle. Extended for use with any streaming pattern.  Frederik Hennig committed Jan 06, 2021 90 91 92  Args: name: optional name of the boundary.  Frederik Hennig committed Nov 09, 2020 93 94  """  Frederik Hennig committed Jan 06, 2021 95 96 97 98  def __init__(self, name=None): """Set an optional name here, to mark boundaries, for example for force evaluations""" super(NoSlip, self).__init__(name)  Frederik Hennig committed Nov 09, 2020 99 100  def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))  Martin Bauer committed May 16, 2017 101 102  def __hash__(self):  Martin Bauer committed Nov 02, 2017 103  return hash(self.name)  Martin Bauer committed May 16, 2017 104 105  def __eq__(self, other):  Martin Bauer committed Nov 02, 2017 106 107 108 109  if not isinstance(other, NoSlip): return False return self.name == other.name  Frederik Hennig committed Jan 20, 2021 110   Frederik Hennig committed Nov 09, 2020 111 # end class NoSlip  Martin Bauer committed May 16, 2017 112   Frederik Hennig committed Nov 09, 2020 113 114  class UBB(LbBoundary):  Frederik Hennig committed Jan 06, 2021 115 116 117 118 119 120 121 122 123 124 125 126  """Velocity bounce back boundary condition, enforcing specified velocity at obstacle Args: velocity: can either be a constant, an access into a field, or a callback function. The callback functions gets a numpy record array with members, 'x','y','z', 'dir' (direction) and 'velocity' which has to be set to the desired velocity of the corresponding link adapt_velocity_to_force: adapts the velocity to the correct equilibrium when the lattice Boltzmann method holds a forcing term. If no forcing term is set and adapt_velocity_to_force is set to True it has no effect. dim: number of spatial dimensions name: optional name of the boundary. """  Martin Bauer committed Feb 09, 2017 127   Helen Schottenhamml committed Mar 24, 2021 128  def __init__(self, velocity, adapt_velocity_to_force=False, dim=None, name=None, data_type='double'):  Martin Bauer committed Dec 11, 2017 129  super(UBB, self).__init__(name)  Martin Bauer committed May 16, 2017 130  self._velocity = velocity  Martin Bauer committed Apr 10, 2018 131  self._adaptVelocityToForce = adapt_velocity_to_force  Martin Bauer committed Oct 21, 2017 132 133 134 135 136  if callable(self._velocity) and not dim: raise ValueError("When using a velocity callback the dimension has to be specified with the dim parameter") elif not callable(self._velocity): dim = len(velocity) self.dim = dim  Helen Schottenhamml committed Mar 24, 2021 137  self.data_type = data_type  Martin Bauer committed Feb 09, 2017 138   Martin Bauer committed Oct 21, 2017 139  @property  Martin Bauer committed Apr 10, 2018 140  def additional_data(self):  Frederik Hennig committed Jan 06, 2021 141 142  """ 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."""  Frederik Hennig committed Feb 02, 2021 143  if self.velocity_is_callable:  Helen Schottenhamml committed Mar 24, 2021 144  return [('vel_%d' % (i,), create_type(self.data_type)) for i in range(self.dim)]  Martin Bauer committed Oct 21, 2017 145 146 147 148  else: return [] @property  Martin Bauer committed Apr 10, 2018 149  def additional_data_init_callback(self):  Frederik Hennig committed Jan 06, 2021 150 151 152  """Initialise additional data of the boundary. For an example see tutorial 02 _ or lbmpy.geometry.add_pipe_inflow_boundary"""  Martin Bauer committed Oct 21, 2017 153 154 155  if callable(self._velocity): return self._velocity  Frederik Hennig committed Nov 09, 2020 156  def get_additional_code_nodes(self, lb_method):  Frederik Hennig committed Jan 06, 2021 157 158 159 160 161 162 163 164  """Return a list of code nodes that will be added in the generated code before the index field loop. Args: lb_method: Lattice Boltzmann method. See :func:lbmpy.creationfunctions.create_lb_method Returns: list containing LbmWeightInfo and NeighbourOffsetArrays """  Frederik Hennig committed Nov 09, 2020 165 166  return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]  Frederik Hennig committed Feb 02, 2021 167 168 169 170 171 172  @property def velocity_is_callable(self): """Returns True is velocity is callable. This means the velocity should be initialised via a callback function. This is useful if the inflow velocity should have a certain profile for instance""" return callable(self._velocity)  Frederik Hennig committed Nov 09, 2020 173  def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):  Martin Bauer committed Apr 10, 2018 174  vel_from_idx_field = callable(self._velocity)  Frederik Hennig committed Nov 09, 2020 175 176  vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity direction = dir_symbol  Martin Bauer committed Apr 11, 2017 177   Frederik Hennig committed Nov 09, 2020 178 179  assert self.dim == lb_method.dim, \ f"Dimension of UBB ({self.dim}) does not match dimension of method ({lb_method.dim})"  Martin Bauer committed Oct 21, 2017 180   Frederik Hennig committed Nov 09, 2020 181  neighbor_offset = NeighbourOffsetArrays.neighbour_offset(direction, lb_method.stencil)  Martin Bauer committed May 16, 2017 182   Frederik Hennig committed Nov 09, 2020 183 184 185  velocity = tuple(v_i.get_shifted(*neighbor_offset) if isinstance(v_i, Field.Access) and not vel_from_idx_field else v_i  Martin Bauer committed Oct 21, 2017 186  for v_i in vel)  Martin Bauer committed May 16, 2017 187 188  if self._adaptVelocityToForce:  Martin Bauer committed Apr 10, 2018 189 190 191  cqc = lb_method.conserved_quantity_computation shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity) velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.first_order_moment_symbols).main_assignments]  Martin Bauer committed May 16, 2017 192 193  c_s_sq = sp.Rational(1, 3)  Martin Bauer committed Apr 10, 2018 194  weight_of_direction = LbmWeightInfo.weight_of_direction  Frederik Hennig committed Jan 20, 2021 195 196  vel_term = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction( direction, lb_method)  Martin Bauer committed May 16, 2017 197 198  # Better alternative: in conserved value computation  Martin Bauer committed Apr 10, 2018 199  # rename what is currently called density to "virtual_density"  Martin Bauer committed May 16, 2017 200  # provide a new quantity density, which is constant in case of incompressible models  Martin Bauer committed Apr 10, 2018 201 202 203  if not lb_method.conserved_quantity_computation.zero_centered_pdfs: cqc = lb_method.conserved_quantity_computation density_symbol = sp.Symbol("rho")  Frederik Hennig committed Nov 09, 2020 204  pdf_field_accesses = [f_out(i) for i in range(len(lb_method.stencil))]  Martin Bauer committed Apr 10, 2018 205 206 207  density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol}) density_symbol = lb_method.conserved_quantity_computation.defined_symbols()['density'] result = density_equations.all_assignments  Frederik Hennig committed Nov 09, 2020 208 209  result += [Assignment(f_in(inv_dir[direction]), f_out(direction) - vel_term * density_symbol)]  Martin Bauer committed May 16, 2017 210 211  return result else:  Frederik Hennig committed Nov 09, 2020 212 213  return [Assignment(f_in(inv_dir[direction]), f_out(direction) - vel_term)]  Frederik Hennig committed Jan 20, 2021 214 215   Frederik Hennig committed Nov 09, 2020 216 # end class UBB  Martin Bauer committed May 16, 2017 217 218   Frederik Hennig committed Jan 06, 2021 219 220 class SimpleExtrapolationOutflow(LbBoundary): r"""  Markus Holzer committed May 12, 2021 221 222 223  Simple Outflow boundary condition :cite:geier2015, equation F.1 (listed below). This boundary condition extrapolates missing populations from the last layer of fluid cells onto the boundary by copying them in the normal direction.  Frederik Hennig committed Jan 06, 2021 224 225 226 227 228 229 230  .. math :: f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yzt} Args: normal_direction: direction vector normal to the outflow  Frederik Hennig committed Jan 20, 2021 231  stencil: stencil used by the lattice boltzmann method  Frederik Hennig committed Jan 06, 2021 232 233 234 235 236 237 238 239 240 241 242 243 244  name: optional name of the boundary. """ def __init__(self, normal_direction, stencil, name=None): if isinstance(normal_direction, str): normal_direction = direction_string_to_offset(normal_direction, dim=len(stencil[0])) if name is None: name = f"Simple Outflow: {offset_to_direction_string(normal_direction)}" self.normal_direction = normal_direction super(SimpleExtrapolationOutflow, self).__init__(name)  Frederik Hennig committed Jan 20, 2021 245 246 247 248 249 250 251 252 253 254 255 256  def get_additional_code_nodes(self, lb_method): """Return a list of code nodes that will be added in the generated code before the index field loop. Args: lb_method: Lattice Boltzmann method. See :func:lbmpy.creationfunctions.create_lb_method Returns: list containing NeighbourOffsetArrays """ return [NeighbourOffsetArrays(lb_method.stencil)]  Frederik Hennig committed Jan 06, 2021 257  def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):  Frederik Hennig committed Jan 20, 2021 258 259 260 261  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)) return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](inv_dir[dir_symbol]))  Frederik Hennig committed Jan 06, 2021 262 263 264 265 266 # end class SimpleExtrapolationOutflow class ExtrapolationOutflow(LbBoundary): r"""  Markus Holzer committed May 12, 2021 267 268 269 270 271  Outflow boundary condition :cite:geier2015, equation F.2, with u neglected (listed below). This boundary condition interpolates populations missing on the boundary in normal direction. For this interpolation, the PDF values of the last time step are used. They are interpolated between fluid cell and boundary cell. To get the PDF values from the last time step an index array is used which stores them.  Frederik Hennig committed Jan 06, 2021 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294  .. math :: f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yz(t - \Delta t)} c \theta^{\frac{1}{2}} \frac{\Delta t}{\Delta x} + \left(1 - c \theta^{\frac{1}{2}} \frac{\Delta t}{\Delta x} \right) f_{\overline{1}jk(x - \Delta x)yzt} Args: normal_direction: direction vector normal to the outflow lb_method: the lattice boltzman method to be used in the simulation dt: lattice time step size dx: lattice spacing distance name: optional name of the boundary. streaming_pattern: Streaming pattern to be used in the simulation zeroth_timestep: for in-place patterns, whether the initial setup corresponds to an even or odd time step initial_density: floating point constant or callback taking spatial coordinates (x, y [,z]) as positional arguments, specifying the initial density on boundary nodes initial_velocity: tuple of floating point constants or callback taking spatial coordinates (x, y [,z]) as positional arguments, specifying the initial velocity on boundary nodes """ def __init__(self, normal_direction, lb_method, dt=1, dx=1, name=None, streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,  Helen Schottenhamml committed Mar 24, 2021 295 296 297  initial_density=None, initial_velocity=None, data_type='double'): self.data_type = data_type  Frederik Hennig committed Jan 06, 2021 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338  self.lb_method = lb_method self.stencil = lb_method.stencil self.dim = len(self.stencil[0]) if isinstance(normal_direction, str): normal_direction = direction_string_to_offset(normal_direction, dim=self.dim) if name is None: name = f"Outflow: {offset_to_direction_string(normal_direction)}" self.normal_direction = normal_direction self.streaming_pattern = streaming_pattern self.zeroth_timestep = zeroth_timestep self.dx = sp.Number(dx) self.dt = sp.Number(dt) self.c = sp.sqrt(sp.Rational(1, 3)) * (self.dx / self.dt) self.initial_density = initial_density self.initial_velocity = initial_velocity self.equilibrium_calculation = None if initial_density and initial_velocity: equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([])) rho = lb_method.zeroth_order_equilibrium_moment_symbol u_vec = lb_method.first_order_equilibrium_moment_symbols eq_lambda = equilibrium.lambdify((rho,) + u_vec) post_pdf_symbols = lb_method.post_collision_pdf_symbols def calc_eq_pdfs(density, velocity, j): return eq_lambda(density, *velocity)[post_pdf_symbols[j]] self.equilibrium_calculation = calc_eq_pdfs super(ExtrapolationOutflow, self).__init__(name) def init_callback(self, boundary_data, **_): dim = boundary_data.dim coord_names = ['x', 'y', 'z'][:dim] pdf_acc = AccessPdfValues(self.stencil, streaming_pattern=self.streaming_pattern, timestep=self.zeroth_timestep, streaming_dir='out')  Frederik Hennig committed Feb 02, 2021 339  def get_boundary_cell_pdfs(f_cell, b_cell, direction):  Frederik Hennig committed Jan 06, 2021 340 341  if self.equilibrium_calculation is not None: density = self.initial_density(  Frederik Hennig committed Feb 02, 2021 342  *b_cell) if callable(self.initial_density) else self.initial_density  Frederik Hennig committed Jan 06, 2021 343  velocity = self.initial_velocity(  Frederik Hennig committed Feb 02, 2021 344 345  *b_cell) if callable(self.initial_velocity) else self.initial_velocity return self.equilibrium_calculation(density, velocity, direction)  Frederik Hennig committed Jan 06, 2021 346  else:  Frederik Hennig committed Feb 02, 2021 347  return pdf_acc.read_pdf(boundary_data.pdf_array, f_cell, direction)  Frederik Hennig committed Jan 06, 2021 348 349  for entry in boundary_data.index_array:  Frederik Hennig committed Jan 20, 2021 350 351 352 353 354 355  center = tuple(entry[c] for c in coord_names) direction = self.stencil[entry["dir"]] inv_dir = self.stencil.index(inverse_direction(direction)) tangential_offset = tuple(offset - normal for offset, normal in zip(direction, self.normal_direction)) domain_cell = tuple(f + o for f, o in zip(center, tangential_offset)) outflow_cell = tuple(f + o for f, o in zip(center, direction))  Frederik Hennig committed Jan 06, 2021 356 357  # Initial fluid cell PDF values  Frederik Hennig committed Jan 20, 2021 358 359  entry['pdf'] = pdf_acc.read_pdf(boundary_data.pdf_array, domain_cell, inv_dir) entry['pdf_nd'] = get_boundary_cell_pdfs(domain_cell, outflow_cell, inv_dir)  Frederik Hennig committed Jan 06, 2021 360 361 362  @property def additional_data(self):  Frederik Hennig committed Jan 20, 2021 363 364  """Used internally only. For the ExtrapolationOutflow information of the previous PDF values is needed. This information is stored in the index vector."""  Helen Schottenhamml committed Mar 24, 2021 365  data = [('pdf', create_type(self.data_type)), ('pdf_nd', create_type(self.data_type))]  Frederik Hennig committed Jan 06, 2021 366 367 368 369  return data @property def additional_data_init_callback(self):  Frederik Hennig committed Jan 20, 2021 370 371 372 373 374 375 376 377 378 379 380 381 382  return self.init_callback def get_additional_code_nodes(self, lb_method): """Return a list of code nodes that will be added in the generated code before the index field loop. Args: lb_method: Lattice Boltzmann method. See :func:lbmpy.creationfunctions.create_lb_method Returns: list containing NeighbourOffsetArrays """ return [NeighbourOffsetArrays(lb_method.stencil)]  Frederik Hennig committed Jan 06, 2021 383 384 385 386 387 388  def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): subexpressions = [] boundary_assignments = [] dtdx = sp.Rational(self.dt, self.dx)  Frederik Hennig committed Jan 20, 2021 389 390  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))  Frederik Hennig committed Jan 06, 2021 391   Frederik Hennig committed Jan 20, 2021 392 393 394 395  interpolated_pdf_sym = sp.Symbol('pdf_inter') interpolated_pdf_asm = Assignment(interpolated_pdf_sym, (index_field[0]('pdf') * (self.c * dtdx)) + ((sp.Number(1) - self.c * dtdx) * index_field[0]('pdf_nd'))) subexpressions.append(interpolated_pdf_asm)  Frederik Hennig committed Jan 06, 2021 396   Frederik Hennig committed Jan 20, 2021 397 398  asm = Assignment(f_in.center(inv_dir[dir_symbol]), interpolated_pdf_sym) boundary_assignments.append(asm)  Frederik Hennig committed Jan 06, 2021 399   Frederik Hennig committed Jan 20, 2021 400 401  asm = Assignment(index_field[0]('pdf'), f_out[tangential_offset](inv_dir[dir_symbol])) boundary_assignments.append(asm)  Frederik Hennig committed Jan 06, 2021 402   Frederik Hennig committed Jan 20, 2021 403 404  asm = Assignment(index_field[0]('pdf_nd'), interpolated_pdf_sym) boundary_assignments.append(asm)  Frederik Hennig committed Jan 06, 2021 405 406  return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)  Frederik Hennig committed Jan 20, 2021 407 408   Frederik Hennig committed Jan 06, 2021 409 410 411 # end class ExtrapolationOutflow  Frederik Hennig committed Nov 09, 2020 412 class FixedDensity(LbBoundary):  Frederik Hennig committed Jan 06, 2021 413 414 415 416 417 418  """Boundary condition that fixes the density/pressure at the obstacle. Args: density: value of the density which should be set. name: optional name of the boundary. """  Martin Bauer committed May 16, 2017 419   Martin Bauer committed Dec 11, 2017 420  def __init__(self, density, name=None):  Martin Bauer committed Apr 18, 2018 421 422  if name is None: name = "Fixed Density " + str(density)  Martin Bauer committed Dec 11, 2017 423  super(FixedDensity, self).__init__(name)  Martin Bauer committed May 16, 2017 424  self._density = density  Martin Bauer committed Feb 13, 2017 425   Frederik Hennig committed Nov 09, 2020 426  def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):  Martin Bauer committed Apr 10, 2018 427 428 429 430  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] return assignment_collection.copy(new_main_assignments)  Martin Bauer committed May 16, 2017 431   Martin Bauer committed Apr 10, 2018 432  cqc = lb_method.conserved_quantity_computation  Martin Bauer committed Apr 10, 2018 433  velocity = cqc.defined_symbols()['velocity']  Martin Bauer committed Apr 10, 2018 434 435  symmetric_eq = remove_asymmetric_part_of_main_assignments(lb_method.get_equilibrium(), degrees_of_freedom=velocity)  Frederik Hennig committed Nov 09, 2020 436  substitutions = {sym: f_out(i) for i, sym in enumerate(lb_method.pre_collision_pdf_symbols)}  Martin Bauer committed Apr 10, 2018 437  symmetric_eq = symmetric_eq.new_with_substitutions(substitutions)  Martin Bauer committed May 16, 2017 438   Martin Bauer committed Apr 10, 2018 439 440  simplification = create_simplification_strategy(lb_method) symmetric_eq = simplification(symmetric_eq)  Martin Bauer committed May 16, 2017 441   Martin Bauer committed Apr 10, 2018 442  density_symbol = cqc.defined_symbols()['density']  Martin Bauer committed May 16, 2017 443 444  density = self._density  Martin Bauer committed Apr 13, 2018 445 446  equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density) equilibrium_input = equilibrium_input.new_without_subexpressions()  Martin Bauer committed Apr 10, 2018 447 448 449  density_eq = equilibrium_input.main_assignments[0] assert density_eq.lhs == density_symbol transformed_density = density_eq.rhs  Martin Bauer committed May 16, 2017 450   Frederik Hennig committed Nov 09, 2020 451  conditions = [(eq_i.rhs, sp.Equality(dir_symbol, i))  Martin Bauer committed Apr 10, 2018 452  for i, eq_i in enumerate(symmetric_eq.main_assignments)] + [(0, True)]  Martin Bauer committed May 16, 2017 453 454  eq_component = sp.Piecewise(*conditions)  Martin Bauer committed Apr 10, 2018 455 456  subexpressions = [Assignment(eq.lhs, transformed_density if eq.lhs == density_symbol else eq.rhs) for eq in symmetric_eq.subexpressions]  Martin Bauer committed Jul 26, 2017 457   Frederik Hennig committed Nov 09, 2020 458 459  return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]), 2 * eq_component - f_out(dir_symbol))]  Frederik Hennig committed Jan 20, 2021 460 461   Frederik Hennig committed Nov 09, 2020 462 # end class FixedDensity  Martin Bauer committed Feb 09, 2017 463   Martin Bauer committed Jul 26, 2017 464   Frederik Hennig committed Nov 09, 2020 465 class NeumannByCopy(LbBoundary):  Frederik Hennig committed Jan 06, 2021 466 467  """Neumann boundary condition which is implemented by coping the PDF values to achieve similar values at the fluid and the boundary node"""  Frederik Hennig committed Nov 09, 2020 468 469  def get_additional_code_nodes(self, lb_method):  Frederik Hennig committed Jan 06, 2021 470 471 472 473 474 475 476 477  """Return a list of code nodes that will be added in the generated code before the index field loop. Args: lb_method: Lattice Boltzmann method. See :func:lbmpy.creationfunctions.create_lb_method Returns: list containing NeighbourOffsetArrays """  Frederik Hennig committed Nov 09, 2020 478 479 480 481 482 483  return [NeighbourOffsetArrays(lb_method.stencil)] def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): 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))]  Martin Bauer committed Jul 26, 2017 484 485 486 487 488 489 490  def __hash__(self): # All boundaries of these class behave equal -> should also be equal return hash("NeumannByCopy") def __eq__(self, other): return type(other) == NeumannByCopy  Frederik Hennig committed Jan 20, 2021 491 492   Frederik Hennig committed Nov 09, 2020 493 # end class NeumannByCopy  Martin Bauer committed Jul 26, 2017 494   Martin Bauer committed Mar 05, 2018 495   Frederik Hennig committed Nov 09, 2020 496 class StreamInConstant(LbBoundary):  Frederik Hennig committed Jan 06, 2021 497 498 499 500 501 502 503  """Boundary condition that takes a constant and overrides the boundary PDFs with this value. This is used for debugging mainly. Args: constant: value which should be set for the PDFs at the boundary cell. name: optional name of the boundary. """  Frederik Hennig committed Jan 20, 2021 504   Martin Bauer committed Apr 30, 2018 505 506 507  def __init__(self, constant, name=None): super(StreamInConstant, self).__init__(name) self._constant = constant  Martin Bauer committed Mar 05, 2018 508   Frederik Hennig committed Nov 09, 2020 509  def get_additional_code_nodes(self, lb_method):  Frederik Hennig committed Jan 06, 2021 510 511 512 513 514 515 516 517  """Return a list of code nodes that will be added in the generated code before the index field loop. Args: lb_method: Lattice Boltzmann method. See :func:lbmpy.creationfunctions.create_lb_method Returns: list containing NeighbourOffsetArrays """  Frederik Hennig committed Nov 09, 2020 518 519 520 521 522 523  return [NeighbourOffsetArrays(lb_method.stencil)] def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field): 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)]  Martin Bauer committed Mar 05, 2018 524 525  def __hash__(self):  Martin Bauer committed Apr 30, 2018 526 527  # All boundaries of these class behave equal -> should also be equal return hash("StreamInConstant")  Martin Bauer committed Mar 05, 2018 528 529  def __eq__(self, other):  Martin Bauer committed Apr 30, 2018 530  return type(other) == StreamInConstant  Frederik Hennig committed Nov 09, 2020 531 # end class StreamInConstant