kernel_equations.py 18.4 KB
 Markus Holzer committed Nov 03, 2019 1 2 3 ``````from pystencils.field import Field from lbmpy.creationfunctions import update_with_default_parameters from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAccessor `````` Markus Holzer committed Sep 20, 2019 4 5 6 7 8 9 10 11 12 13 14 15 16 ``````from pystencils.fd.derivation import FiniteDifferenceStencilDerivation from lbmpy.maxwellian_equilibrium import get_weights from pystencils import Assignment, AssignmentCollection import sympy as sp import numpy as np def chemical_potential_symbolic(phi_field, stencil, beta, kappa): r""" Get a symbolic expression for the chemical potential according to equation (5) in PhysRevE.96.053301. Args: phi_field: the phase-field on which the chemical potential is applied `````` Markus Holzer committed Dec 08, 2019 17 `````` stencil: stencil to derive the finite difference for the laplacian (2nd order isotropic) `````` Markus Holzer committed Sep 20, 2019 18 19 20 21 `````` beta: coefficient related to surface tension and interface thickness kappa: coefficient related to surface tension and interface thickness """ dimensions = len(stencil[0]) `````` Markus Holzer committed Dec 08, 2019 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 `````` q = len(stencil) lap = sp.simplify(0) for i in range(dimensions): deriv = FiniteDifferenceStencilDerivation((i, i), stencil) for j in range(dimensions): # assume the stencil is symmetric deriv.assume_symmetric(dim=j, anti_symmetric=False) # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic if q == 9: res = deriv.get_stencil(isotropic=True) lap += res.apply(phi_field.center) elif q == 15: deriv.set_weight((0, 0, 0), sp.Rational(-32, 27)) res = deriv.get_stencil(isotropic=True) lap += res.apply(phi_field.center) elif q == 19: res = deriv.get_stencil(isotropic=True) lap += res.apply(phi_field.center) else: deriv.set_weight((0, 0, 0), sp.Rational(-38, 27)) res = deriv.get_stencil(isotropic=True) lap += res.apply(phi_field.center) `````` Markus Holzer committed Sep 20, 2019 45 46 47 `````` # get the chemical potential mu = 4.0 * beta * phi_field.center * (phi_field.center - 1.0) * (phi_field.center - 0.5) - \ `````` Markus Holzer committed Nov 03, 2019 48 `````` kappa * lap `````` Markus Holzer committed Sep 20, 2019 49 50 51 52 53 54 55 56 `````` return mu def isotropic_gradient_symbolic(phi_field, stencil): r""" Get a symbolic expression for the isotropic gradient of the phase-field Args: phi_field: the phase-field on which the isotropic gradient is applied `````` Markus Holzer committed Dec 08, 2019 57 `````` stencil: stencil to derive the finite difference for the gradient (2nd order isotropic) `````` Markus Holzer committed Sep 20, 2019 58 59 `````` """ dimensions = len(stencil[0]) `````` Markus Holzer committed Dec 08, 2019 60 `````` q = len(stencil) `````` Markus Holzer committed Nov 03, 2019 61 `````` deriv = FiniteDifferenceStencilDerivation((0,), stencil) `````` Markus Holzer committed Sep 20, 2019 62 63 64 65 66 67 `````` deriv.assume_symmetric(0, anti_symmetric=True) deriv.assume_symmetric(1, anti_symmetric=False) if dimensions == 3: deriv.assume_symmetric(2, anti_symmetric=False) `````` Markus Holzer committed Dec 08, 2019 68 69 70 `````` # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic # furthermore the stencils gets rotated to get the y and z components if q == 9: `````` Markus Holzer committed Nov 03, 2019 71 72 `````` res = deriv.get_stencil(isotropic=True) grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0] `````` Markus Holzer committed Dec 08, 2019 73 74 75 76 77 78 79 `````` elif q == 15: res = deriv.get_stencil(isotropic=True) grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), res.rotate_weights_and_apply(phi_field.center, (1, 2))] elif q == 19: deriv.set_weight((0, 0, 0), sp.sympify(0)) `````` Markus Holzer committed Sep 20, 2019 80 81 `````` deriv.set_weight((1, 0, 0), sp.Rational(1, 6)) `````` Markus Holzer committed Nov 03, 2019 82 `````` res = deriv.get_stencil(isotropic=True) `````` Markus Holzer committed Sep 20, 2019 83 `````` grad = [res.apply(phi_field.center), `````` Markus Holzer committed Dec 08, 2019 84 85 `````` res.rotate_weights_and_apply(phi_field.center, (0, 1)), res.rotate_weights_and_apply(phi_field.center, (1, 2))] `````` Markus Holzer committed Nov 03, 2019 86 `````` else: `````` Markus Holzer committed Dec 08, 2019 87 88 89 90 91 92 93 `````` deriv.set_weight((0, 0, 0), sp.sympify(0)) deriv.set_weight((1, 0, 0), sp.Rational(2, 9)) res = deriv.get_stencil(isotropic=True) grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), res.rotate_weights_and_apply(phi_field.center, (1, 2))] `````` Markus Holzer committed Sep 20, 2019 94 95 96 97 `````` return grad `````` Markus Holzer committed Dec 08, 2019 98 ``````def normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil=None): `````` Markus Holzer committed Sep 20, 2019 99 100 101 102 103 `````` r""" Get a symbolic expression for the normalized isotropic gradient of the phase-field Args: phi_field: the phase-field on which the normalized isotropic gradient is applied stencil: stencil of the lattice Boltzmann step `````` Markus Holzer committed Dec 08, 2019 104 105 `````` fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase field. If it is not given the stencil of the LB method will be applied. `````` Markus Holzer committed Sep 20, 2019 106 `````` """ `````` Markus Holzer committed Dec 08, 2019 107 108 `````` if fd_stencil is None: fd_stencil = stencil `````` Markus Holzer committed Sep 20, 2019 109 `````` `````` Markus Holzer committed Dec 08, 2019 110 111 112 `````` tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, fd_stencil))) + 1.e-32) ** 0.5 result = [x / tmp for x in isotropic_gradient_symbolic(phi_field, fd_stencil)] `````` Markus Holzer committed Sep 20, 2019 113 114 115 `````` return result `````` Markus Holzer committed Dec 08, 2019 116 ``````def pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil=None): `````` Markus Holzer committed Sep 20, 2019 117 118 119 120 121 122 123 `````` r""" Get a symbolic expression for the pressure force Args: phi_field: phase-field stencil: stencil of the lattice Boltzmann step density_heavy: density of the heavier fluid density_light: density of the lighter fluid `````` Markus Holzer committed Dec 08, 2019 124 125 `````` fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase field. If it is not given the stencil of the LB method will be applied. `````` Markus Holzer committed Sep 20, 2019 126 `````` """ `````` Markus Holzer committed Dec 08, 2019 127 128 129 130 `````` if fd_stencil is None: fd_stencil = stencil iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil) `````` Markus Holzer committed Nov 03, 2019 131 `````` result = list(map(lambda x: sp.Rational(-1, 3) * sp.symbols("rho") * (density_heavy - density_light) * x, iso_grad)) `````` Markus Holzer committed Sep 20, 2019 132 133 134 `````` return result `````` Markus Holzer committed Dec 08, 2019 135 ``````def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light, fd_stencil=None): `````` Markus Holzer committed Sep 20, 2019 136 137 138 139 140 141 142 143 144 `````` r""" Get a symbolic expression for the viscous force Args: lb_velocity_field: hydrodynamic distribution function phi_field: phase-field mrt_method: mrt lattice boltzmann method used for hydrodynamics tau: relaxation time of the hydrodynamic lattice boltzmann step density_heavy: density of the heavier fluid density_light: density of the lighter fluid `````` Markus Holzer committed Dec 08, 2019 145 146 `````` fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase field. If it is not given the stencil of the LB method will be applied. `````` Markus Holzer committed Sep 20, 2019 147 148 149 150 `````` """ stencil = mrt_method.stencil dimensions = len(stencil[0]) `````` Markus Holzer committed Dec 08, 2019 151 152 153 154 `````` if fd_stencil is None: fd_stencil = stencil iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil) `````` Markus Holzer committed Sep 20, 2019 155 156 `````` moment_matrix = mrt_method.moment_matrix `````` Markus Holzer committed Nov 03, 2019 157 158 159 160 161 162 `````` rel = mrt_method.relaxation_rates eq = mrt_method.moment_equilibrium_values eq = np.array(eq) g_vals = [lb_velocity_field.center(i) for i, _ in enumerate(stencil)] m0 = np.dot(moment_matrix, g_vals) `````` Markus Holzer committed Sep 20, 2019 163 `````` `````` Markus Holzer committed Nov 03, 2019 164 165 166 `````` m = m0 - eq m = m * rel non_equilibrium = np.dot(moment_matrix.inv(), m) `````` Markus Holzer committed Sep 20, 2019 167 168 169 170 `````` stress_tensor = [0] * 6 # Calculate Stress Tensor MRT for i, d in enumerate(stencil): `````` Markus Holzer committed Nov 03, 2019 171 172 `````` stress_tensor[0] = sp.Add(stress_tensor[0], non_equilibrium[i] * (d[0] * d[0])) stress_tensor[1] = sp.Add(stress_tensor[1], non_equilibrium[i] * (d[1] * d[1])) `````` Markus Holzer committed Sep 20, 2019 173 174 `````` if dimensions == 3: `````` Markus Holzer committed Nov 03, 2019 175 176 177 `````` stress_tensor[2] = sp.Add(stress_tensor[2], non_equilibrium[i] * (d[2] * d[2])) stress_tensor[3] = sp.Add(stress_tensor[3], non_equilibrium[i] * (d[1] * d[2])) stress_tensor[4] = sp.Add(stress_tensor[4], non_equilibrium[i] * (d[0] * d[2])) `````` Markus Holzer committed Sep 20, 2019 178 `````` `````` Markus Holzer committed Nov 03, 2019 179 `````` stress_tensor[5] = sp.Add(stress_tensor[5], non_equilibrium[i] * (d[0] * d[1])) `````` Markus Holzer committed Sep 20, 2019 180 181 182 183 `````` density_difference = density_heavy - density_light # Calculate Viscous Force MRT `````` Markus Holzer committed Nov 03, 2019 184 185 186 `````` fmx = (0.5 - tau) * (stress_tensor[0] * iso_grad[0] + stress_tensor[5] * iso_grad[1] + stress_tensor[4] * iso_grad[2]) * density_difference `````` Markus Holzer committed Sep 20, 2019 187 `````` `````` Markus Holzer committed Nov 03, 2019 188 189 190 `````` fmy = (0.5 - tau) * (stress_tensor[5] * iso_grad[0] + stress_tensor[1] * iso_grad[1] + stress_tensor[3] * iso_grad[2]) * density_difference `````` Markus Holzer committed Sep 20, 2019 191 `````` `````` Markus Holzer committed Nov 03, 2019 192 193 194 `````` fmz = (0.5 - tau) * (stress_tensor[4] * iso_grad[0] + stress_tensor[3] * iso_grad[1] + stress_tensor[2] * iso_grad[2]) * density_difference `````` Markus Holzer committed Sep 20, 2019 195 196 197 198 `````` return [fmx, fmy, fmz] `````` Markus Holzer committed Dec 08, 2019 199 ``````def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None): `````` Markus Holzer committed Sep 20, 2019 200 201 202 203 204 205 206 `````` r""" Get a symbolic expression for the surface tension force Args: phi_field: the phase-field on which the chemical potential is applied stencil: stencil of the lattice Boltzmann step beta: coefficient related to surface tension and interface thickness kappa: coefficient related to surface tension and interface thickness `````` Markus Holzer committed Dec 08, 2019 207 208 `````` fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase field. If it is not given the stencil of the LB method will be applied. `````` Markus Holzer committed Sep 20, 2019 209 `````` """ `````` Markus Holzer committed Dec 08, 2019 210 211 212 213 214 `````` if fd_stencil is None: fd_stencil = stencil chemical_potential = chemical_potential_symbolic(phi_field, fd_stencil, beta, kappa) iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil) `````` Markus Holzer committed Nov 03, 2019 215 `````` return [chemical_potential * x for x in iso_grad] `````` Markus Holzer committed Sep 20, 2019 216 217 `````` `````` Markus Holzer committed Nov 03, 2019 218 ``````def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau, `````` Markus Holzer committed Dec 08, 2019 219 `````` density_heavy, density_light, kappa, beta, body_force, fd_stencil=None): `````` Markus Holzer committed Sep 20, 2019 220 221 222 223 224 `````` r""" Get a symbolic expression for the hydrodynamic force Args: lb_velocity_field: hydrodynamic distribution function phi_field: phase-field `````` Markus Holzer committed Nov 03, 2019 225 `````` lb_method: Lattice boltzmann method used for hydrodynamics `````` Markus Holzer committed Sep 20, 2019 226 227 228 229 230 231 `````` tau: relaxation time of the hydrodynamic lattice boltzmann step density_heavy: density of the heavier fluid density_light: density of the lighter fluid beta: coefficient related to surface tension and interface thickness kappa: coefficient related to surface tension and interface thickness body_force: force acting on the fluids. Usually the gravity `````` Markus Holzer committed Dec 08, 2019 232 233 `````` fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase field. If it is not given the stencil of the LB method will be applied. `````` Markus Holzer committed Sep 20, 2019 234 `````` """ `````` Markus Holzer committed Nov 03, 2019 235 `````` stencil = lb_method.stencil `````` Markus Holzer committed Sep 20, 2019 236 `````` dimensions = len(stencil[0]) `````` Markus Holzer committed Dec 08, 2019 237 238 239 240 241 242 243 `````` if fd_stencil is None: fd_stencil = stencil fp = pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil) fm = viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light, fd_stencil) fs = surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil) `````` Markus Holzer committed Sep 20, 2019 244 245 246 247 248 249 250 251 `````` result = [] for i in range(dimensions): result.append(fs[i] + fp[i] + fm[i] + body_force[i]) return result `````` Markus Holzer committed Dec 08, 2019 252 ``````def interface_tracking_force(phi_field, stencil, interface_thickness, fd_stencil=None): `````` Markus Holzer committed Sep 20, 2019 253 254 255 256 257 258 `````` r""" Get a symbolic expression for the hydrodynamic force Args: phi_field: phase-field stencil: stencil of the phase-field distribution lattice Boltzmann step interface_thickness: interface thickness `````` Markus Holzer committed Dec 08, 2019 259 260 `````` fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase field. If it is not given the stencil of the LB method will be applied. `````` Markus Holzer committed Sep 20, 2019 261 `````` """ `````` Markus Holzer committed Dec 08, 2019 262 263 264 `````` if fd_stencil is None: fd_stencil = stencil `````` Markus Holzer committed Sep 20, 2019 265 `````` dimensions = len(stencil[0]) `````` Markus Holzer committed Dec 08, 2019 266 `````` normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil) `````` Markus Holzer committed Sep 20, 2019 267 268 `````` result = [] for i in range(dimensions): `````` Martin Bauer committed Nov 28, 2019 269 `````` result.append(((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness) * normal_fd[i]) `````` Markus Holzer committed Sep 20, 2019 270 271 272 273 `````` return result `````` Markus Holzer committed Nov 03, 2019 274 ``````def get_update_rules_velocity(src_field, u_in, lb_method, force, density): `````` Markus Holzer committed Sep 20, 2019 275 `````` r""" `````` Markus Holzer committed Nov 03, 2019 276 277 278 279 280 281 282 283 284 `````` Get assignments to update the velocity with a force shift Args: src_field: the source field of the hydrodynamic distribution function u_in: velocity field lb_method: mrt lattice boltzmann method used for hydrodynamics force: force acting on the hydrodynamic lb step density: the interpolated density of the simulation """ stencil = lb_method.stencil `````` Markus Holzer committed Sep 20, 2019 285 286 `````` dimensions = len(stencil[0]) `````` Markus Holzer committed Nov 03, 2019 287 288 289 290 291 292 293 294 295 296 297 298 299 `````` moment_matrix = lb_method.moment_matrix eq = lb_method.moment_equilibrium_values first_eqs = lb_method.first_order_equilibrium_moment_symbols indices = list() for i in range(dimensions): indices.append(eq.index(first_eqs[i])) src = [src_field.center(i) for i, _ in enumerate(stencil)] m0 = np.dot(moment_matrix, src) update_u = list() update_u.append(Assignment(sp.symbols("rho"), m0[0])) `````` Markus Holzer committed Sep 20, 2019 300 301 `````` u_symp = sp.symbols("u_:{}".format(dimensions)) `````` Markus Holzer committed Nov 03, 2019 302 303 304 305 306 307 308 309 310 311 `````` zw = sp.symbols("zw_:{}".format(dimensions)) for i in range(dimensions): update_u.append(Assignment(zw[i], u_in.center_vector[i])) subs_dict = dict(zip(u_symp, zw)) for i in range(dimensions): update_u.append(Assignment(u_in.center_vector[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2)) return update_u `````` Markus Holzer committed Sep 20, 2019 312 `````` `````` Markus Holzer committed Dec 08, 2019 313 ``````def get_collision_assignments_hydro(density=1, optimization=None, **kwargs): `````` Markus Holzer committed Nov 03, 2019 314 315 316 317 318 319 320 `````` r""" Get collision assignments for the hydrodynamic lattice Boltzmann step. Here the force gets applied in the moment space. Afterwards the transformation back to the pdf space happens. Args: density: the interpolated density of the simulation optimization: for details see createfunctions.py """ `````` Markus Holzer committed Dec 08, 2019 321 322 `````` if optimization is None: optimization = {} `````` Markus Holzer committed Nov 03, 2019 323 324 325 326 327 328 329 330 331 `````` params, opt_params = update_with_default_parameters(kwargs, optimization) lb_method = params['lb_method'] stencil = lb_method.stencil dimensions = len(stencil[0]) field_data_type = 'float64' if opt_params['double_precision'] else 'float32' q = len(stencil) `````` Markus Holzer committed Sep 20, 2019 332 `````` `````` Markus Holzer committed Nov 03, 2019 333 334 `````` u_in = params['velocity_input'] force = params['force'] `````` Markus Holzer committed Sep 20, 2019 335 `````` `````` Markus Holzer committed Nov 03, 2019 336 337 338 `````` if opt_params['symbolic_field'] is not None: src_field = opt_params['symbolic_field'] else: `````` Markus Holzer committed Dec 08, 2019 339 `````` src_field = Field.create_generic(params['field_name'], spatial_dimensions=lb_method.dim, `````` Markus Holzer committed Nov 03, 2019 340 341 342 343 344 345 346 347 348 349 350 351 352 `````` index_shape=(q,), layout=opt_params['field_layout'], dtype=field_data_type) if opt_params['symbolic_temporary_field'] is not None: dst_field = opt_params['symbolic_temporary_field'] else: dst_field = src_field.new_field_with_different_name(params['temporary_field_name']) moment_matrix = lb_method.moment_matrix rel = lb_method.relaxation_rates eq = lb_method.moment_equilibrium_values first_eqs = lb_method.first_order_equilibrium_moment_symbols indices = list() `````` Markus Holzer committed Sep 20, 2019 353 `````` for i in range(dimensions): `````` Markus Holzer committed Nov 03, 2019 354 355 356 `````` indices.append(eq.index(first_eqs[i])) eq = np.array(eq) `````` Markus Holzer committed Sep 20, 2019 357 `````` `````` Markus Holzer committed Nov 03, 2019 358 359 `````` g_vals = [src_field.center(i) for i, _ in enumerate(stencil)] m0 = np.dot(moment_matrix, g_vals) `````` Markus Holzer committed Sep 20, 2019 360 `````` `````` Markus Holzer committed Nov 03, 2019 361 `````` mf = np.zeros(len(stencil), dtype=object) `````` Markus Holzer committed Sep 20, 2019 362 `````` for i in range(dimensions): `````` Markus Holzer committed Nov 03, 2019 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 `````` mf[indices[i]] = force[i] / density m = sp.symbols("m_:{}".format(len(stencil))) update_m = get_update_rules_velocity(src_field, u_in, lb_method, force, density) u_symp = sp.symbols("u_:{}".format(dimensions)) for i in range(dimensions): update_m.append(Assignment(u_symp[i], u_in.center_vector[i])) for i in range(0, len(stencil)): update_m.append(Assignment(m[i], m0[i] - (m0[i] - eq[i] + mf[i] / 2) * rel[i] + mf[i])) update_g = list() var = np.dot(moment_matrix.inv(), m) if params['kernel_type'] == 'collide_stream_push': push_accessor = StreamPushTwoFieldsAccessor() post_collision_accesses = push_accessor.write(dst_field, stencil) else: collide_accessor = CollideOnlyInplaceAccessor() post_collision_accesses = collide_accessor.write(src_field, stencil) `````` Markus Holzer committed Sep 20, 2019 384 `````` `````` Markus Holzer committed Nov 03, 2019 385 386 `````` for i in range(0, len(stencil)): update_g.append(Assignment(post_collision_accesses[i], var[i])) `````` Markus Holzer committed Sep 20, 2019 387 `````` `````` Markus Holzer committed Nov 03, 2019 388 389 `````` hydro_lb_update_rule = AssignmentCollection(main_assignments=update_g, subexpressions=update_m) `````` Markus Holzer committed Sep 20, 2019 390 `````` `````` Markus Holzer committed Nov 03, 2019 391 `````` return hydro_lb_update_rule `````` Markus Holzer committed Sep 20, 2019 392 393 `````` `````` Markus Holzer committed Dec 08, 2019 394 395 ``````def initializer_kernel_phase_field_lb(lb_phase_field, phi_field, velocity_field, mrt_method, interface_thickness, fd_stencil=None): `````` Markus Holzer committed Sep 20, 2019 396 397 398 `````` r""" Returns an assignment list for initializing the phase-field distribution functions Args: `````` Markus Holzer committed Dec 08, 2019 399 `````` lb_phase_field: source field of phase-field distribution function `````` Markus Holzer committed Sep 20, 2019 400 `````` phi_field: phase-field `````` Markus Holzer committed Dec 08, 2019 401 `````` velocity_field: velocity field `````` Markus Holzer committed Sep 20, 2019 402 403 `````` mrt_method: lattice Boltzmann method of the phase-field lattice Boltzmann step interface_thickness: interface thickness `````` Markus Holzer committed Dec 08, 2019 404 405 `````` fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase field. If it is not given the stencil of the LB method will be applied. `````` Markus Holzer committed Sep 20, 2019 406 407 408 `````` """ stencil = mrt_method.stencil dimensions = len(stencil[0]) `````` Markus Holzer committed Dec 08, 2019 409 410 411 412 `````` if fd_stencil is None: fd_stencil = stencil `````` Markus Holzer committed Sep 20, 2019 413 414 415 `````` weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3)) u_symp = sp.symbols("u_:{}".format(dimensions)) `````` Markus Holzer committed Dec 08, 2019 416 `````` normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil) `````` Markus Holzer committed Sep 20, 2019 417 418 419 `````` gamma = mrt_method.get_equilibrium_terms() gamma = gamma.subs({sp.symbols("rho"): 1}) `````` Markus Holzer committed Dec 08, 2019 420 `````` gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)}) `````` Markus Holzer committed Nov 03, 2019 421 `````` # create the kernels for the initialization of the h field `````` Markus Holzer committed Sep 20, 2019 422 423 424 425 426 427 428 `````` h_updates = list() def scalar_product(a, b): return sum(a_i * b_i for a_i, b_i in zip(a, b)) f = [] for i, d in enumerate(stencil): `````` Martin Bauer committed Nov 28, 2019 429 `````` f.append(weights[i] * ((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness) `````` Markus Holzer committed Sep 20, 2019 430 431 432 `````` * scalar_product(d, normal_fd[0:dimensions])) for i, _ in enumerate(stencil): `````` Markus Holzer committed Dec 08, 2019 433 `````` h_updates.append(Assignment(lb_phase_field.center(i), phi_field.center * gamma_init[i] - 0.5 * f[i])) `````` Markus Holzer committed Sep 20, 2019 434 435 436 437 `````` return h_updates `````` Markus Holzer committed Dec 08, 2019 438 ``````def initializer_kernel_hydro_lb(lb_velocity_field, velocity_field, mrt_method): `````` Markus Holzer committed Sep 20, 2019 439 440 441 `````` r""" Returns an assignment list for initializing the velocity distribution functions Args: `````` Markus Holzer committed Dec 08, 2019 442 443 `````` lb_velocity_field: source field of velocity distribution function velocity_field: velocity field `````` Markus Holzer committed Sep 20, 2019 444 445 446 447 448 449 450 451 452 `````` mrt_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step """ stencil = mrt_method.stencil dimensions = len(stencil[0]) weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3)) u_symp = sp.symbols("u_:{}".format(dimensions)) gamma = mrt_method.get_equilibrium_terms() gamma = gamma.subs({sp.symbols("rho"): 1}) `````` Markus Holzer committed Dec 08, 2019 453 `````` gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)}) `````` Markus Holzer committed Sep 20, 2019 454 455 456 `````` g_updates = list() for i, _ in enumerate(stencil): `````` Markus Holzer committed Dec 08, 2019 457 `````` g_updates.append(Assignment(lb_velocity_field.center(i), gamma_init[i] - weights[i])) `````` Markus Holzer committed Sep 20, 2019 458 459 `````` return g_updates``````