add diffusion test

lbmpy_tests/test_diffusion.py

+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
+from lbmpy.session                   import *
+from lbmpy.stencils                  import get_stencil
+import math
+
+def test_diffusion():
+  """
+  Runs the "Diffusion from Plate in Uniform Flow" benchmark as it is described
+  in [ch. 8.6.3, The Lattice Boltzmann Method, Krüger et al.].
+
+            dC/dy = 0
+        ┌───────────────┐
+        │     → → →     │
+  C = 0 │     → u →     │ dC/dx = 0
+        │     → → →     │
+        └───────────────┘
+              C = 1
+
+  The analytical solution is given by:
+    C(x,y) = 1 * erfc(y / sqrt(4Dx/u))
+
+  The hydrodynamic field is not simulated, instead a constant velocity is assumed.
+  """
+
+  # Parameters
+  domain_size = (1600, 160)
+  omega       = 1.38
+  D           = (1/omega - 0.5) / 3
+  velocity    = 0.05
+  time_steps  = 50000
+  stencil     = get_stencil('D2Q9')
+
+  # Data Handling
+  dh = ps.create_data_handling(domain_size=domain_size)
+
+  vel_field = dh.add_array('vel_field', values_per_cell=dh.dim)
+  dh.fill('vel_field', velocity, 0, ghost_layers=True)
+  dh.fill('vel_field', 0.0     , 1, ghost_layers=True)
+
+  con_field = dh.add_array('con_field', values_per_cell=1)
+  dh.fill('con_field', 0.0, ghost_layers=True)
+
+  pdfs     = dh.add_array('pdfs'    , values_per_cell=len(stencil))
+  pdfs_tmp = dh.add_array('pdfs_tmp', values_per_cell=len(stencil))
+
+  # Lattice Boltzmann method
+  params = { 'stencil'          : stencil,
+             'method'           : 'mrt',
+             'relaxation_rates' : [1, 1.5, 1, 1.5, 1],
+             'velocity_input'   : vel_field,
+             'output'           : {'density' : con_field},
+             'compressible'     : True,
+             'weighted'         : True,
+             'optimization'     : {'symbolic_field'           : pdfs,
+                                   'symbolic_temporary_field' : pdfs_tmp},
+             'kernel_type'      : 'stream_pull_collide'
+           }
+
+  method = create_lb_method(**params)
+  method.set_conserved_moments_relaxation_rate(omega)
+
+  update_rule = create_lb_update_rule(lb_method=method, **params)
+  kernel = ps.create_kernel(update_rule, target=dh.default_target).compile()
+
+  # PDF initalization
+  init = pdf_initialization_assignments(method, con_field.center,
+                                        vel_field.center_vector, pdfs.center_vector)
+  dh.run_kernel(ps.create_kernel(init, target=dh.default_target).compile())
+
+  # Boundary Handling
+  bh = LatticeBoltzmannBoundaryHandling(update_rule.method, dh, 'pdfs', name="bh")
+  add_box_boundary(bh, boundary=NeumannByCopy())
+  bh.set_boundary(DiffusionDirichlet(0), slice_from_direction('W', dh.dim))
+  bh.set_boundary(DiffusionDirichlet(1), slice_from_direction('S', dh.dim))
+
+  # Timeloop
+  for i in range(time_steps):
+      bh()
+      dh.run_kernel(kernel)
+      dh.swap("pdfs", "pdfs_tmp")
+
+  # Verification
+  x = np.arange(1, domain_size[0], 1)
+  y = np.arange(0, domain_size[1], 1)
+  X, Y = np.meshgrid(x, y)
+  analytical       = np.zeros(domain_size)
+  analytical[1:,:] = np.vectorize(math.erfc)(Y / np.vectorize(math.sqrt)(4 * D * X / velocity)).transpose()
+  simulated = dh.gather_array('con_field', ghost_layers=False)
+
+  residual = 0
+  for i in x:
+    for j in y:
+      residual += (simulated[i,j] - analytical[i,j])**2
+  residual = math.sqrt(residual / (domain_size[0] * domain_size[1]))
+
+  assert residual < 1e-2