kernel_equations.py 18.4 KB
Newer Older
Markus Holzer's avatar
Markus Holzer committed
1
2
3
from pystencils.field import Field
from lbmpy.creationfunctions import update_with_default_parameters
from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAccessor
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
17
        stencil: stencil to derive the finite difference for the laplacian (2nd order isotropic)
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])
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)
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's avatar
Markus Holzer committed
48
        kappa * lap
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
57
        stencil: stencil to derive the finite difference for the gradient (2nd order isotropic)
58
59
    """
    dimensions = len(stencil[0])
60
    q = len(stencil)
Markus Holzer's avatar
Markus Holzer committed
61
    deriv = FiniteDifferenceStencilDerivation((0,), stencil)
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)

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's avatar
Markus Holzer committed
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]
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))
80
81
        deriv.set_weight((1, 0, 0), sp.Rational(1, 6))

Markus Holzer's avatar
Markus Holzer committed
82
        res = deriv.get_stencil(isotropic=True)
83
        grad = [res.apply(phi_field.center),
84
85
                res.rotate_weights_and_apply(phi_field.center, (0, 1)),
                res.rotate_weights_and_apply(phi_field.center, (1, 2))]
Markus Holzer's avatar
Markus Holzer committed
86
    else:
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))]
94
95
96
97

    return grad


98
def normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil=None):
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
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.
106
    """
107
108
    if fd_stencil is None:
        fd_stencil = stencil
109

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)]
113
114
115
    return result


116
def pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil=None):
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
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.
126
    """
127
128
129
130
    if fd_stencil is None:
        fd_stencil = stencil

    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
Markus Holzer's avatar
Markus Holzer committed
131
    result = list(map(lambda x: sp.Rational(-1, 3) * sp.symbols("rho") * (density_heavy - density_light) * x, iso_grad))
132
133
134
    return result


135
def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light, fd_stencil=None):
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
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.
147
148
149
150
    """
    stencil = mrt_method.stencil
    dimensions = len(stencil[0])

151
152
153
154
    if fd_stencil is None:
        fd_stencil = stencil

    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
155
156

    moment_matrix = mrt_method.moment_matrix
Markus Holzer's avatar
Markus Holzer committed
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)
163

Markus Holzer's avatar
Markus Holzer committed
164
165
166
    m = m0 - eq
    m = m * rel
    non_equilibrium = np.dot(moment_matrix.inv(), m)
167
168
169
170

    stress_tensor = [0] * 6
    # Calculate Stress Tensor MRT
    for i, d in enumerate(stencil):
Markus Holzer's avatar
Markus Holzer committed
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]))
173
174

        if dimensions == 3:
Markus Holzer's avatar
Markus Holzer committed
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]))
178

Markus Holzer's avatar
Markus Holzer committed
179
        stress_tensor[5] = sp.Add(stress_tensor[5], non_equilibrium[i] * (d[0] * d[1]))
180
181
182
183

    density_difference = density_heavy - density_light

    # Calculate Viscous Force MRT
Markus Holzer's avatar
Markus Holzer committed
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
187

Markus Holzer's avatar
Markus Holzer committed
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
191

Markus Holzer's avatar
Markus Holzer committed
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
195
196
197
198

    return [fmx, fmy, fmz]


199
def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None):
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
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.
209
    """
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's avatar
Markus Holzer committed
215
    return [chemical_potential * x for x in iso_grad]
216
217


Markus Holzer's avatar
Markus Holzer committed
218
def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
219
                       density_heavy, density_light, kappa, beta, body_force, fd_stencil=None):
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's avatar
Markus Holzer committed
225
        lb_method: Lattice boltzmann method used for hydrodynamics
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
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.
234
    """
Markus Holzer's avatar
Markus Holzer committed
235
    stencil = lb_method.stencil
236
    dimensions = len(stencil[0])
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)
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


252
def interface_tracking_force(phi_field, stencil, interface_thickness, fd_stencil=None):
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
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.
261
    """
262
263
264
    if fd_stencil is None:
        fd_stencil = stencil

265
    dimensions = len(stencil[0])
266
    normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
267
268
    result = []
    for i in range(dimensions):
269
        result.append(((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness) * normal_fd[i])
270
271
272
273

    return result


Markus Holzer's avatar
Markus Holzer committed
274
def get_update_rules_velocity(src_field, u_in, lb_method, force, density):
275
    r"""
Markus Holzer's avatar
Markus Holzer committed
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
285
286
    dimensions = len(stencil[0])

Markus Holzer's avatar
Markus Holzer committed
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]))
300
301

    u_symp = sp.symbols("u_:{}".format(dimensions))
Markus Holzer's avatar
Markus Holzer committed
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

312

313
def get_collision_assignments_hydro(density=1, optimization=None, **kwargs):
Markus Holzer's avatar
Markus Holzer committed
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
     """
321
322
    if optimization is None:
        optimization = {}
Markus Holzer's avatar
Markus Holzer committed
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)
332

Markus Holzer's avatar
Markus Holzer committed
333
334
    u_in = params['velocity_input']
    force = params['force']
335

Markus Holzer's avatar
Markus Holzer committed
336
337
338
    if opt_params['symbolic_field'] is not None:
        src_field = opt_params['symbolic_field']
    else:
339
        src_field = Field.create_generic(params['field_name'], spatial_dimensions=lb_method.dim,
Markus Holzer's avatar
Markus Holzer committed
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()
353
    for i in range(dimensions):
Markus Holzer's avatar
Markus Holzer committed
354
355
356
        indices.append(eq.index(first_eqs[i]))

    eq = np.array(eq)
357

Markus Holzer's avatar
Markus Holzer committed
358
359
    g_vals = [src_field.center(i) for i, _ in enumerate(stencil)]
    m0 = np.dot(moment_matrix, g_vals)
360

Markus Holzer's avatar
Markus Holzer committed
361
    mf = np.zeros(len(stencil), dtype=object)
362
    for i in range(dimensions):
Markus Holzer's avatar
Markus Holzer committed
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)
384

Markus Holzer's avatar
Markus Holzer committed
385
386
    for i in range(0, len(stencil)):
        update_g.append(Assignment(post_collision_accesses[i], var[i]))
387

Markus Holzer's avatar
Markus Holzer committed
388
389
    hydro_lb_update_rule = AssignmentCollection(main_assignments=update_g,
                                                subexpressions=update_m)
390

Markus Holzer's avatar
Markus Holzer committed
391
    return hydro_lb_update_rule
392
393


394
395
def initializer_kernel_phase_field_lb(lb_phase_field, phi_field, velocity_field, mrt_method, interface_thickness,
                                      fd_stencil=None):
396
397
398
    r"""
    Returns an assignment list for initializing the phase-field distribution functions
    Args:
399
        lb_phase_field: source field of phase-field distribution function
400
        phi_field: phase-field
401
        velocity_field: velocity field
402
403
        mrt_method: lattice Boltzmann method of the phase-field lattice Boltzmann step
        interface_thickness: interface thickness
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.
406
407
408
    """
    stencil = mrt_method.stencil
    dimensions = len(stencil[0])
409
410
411
412

    if fd_stencil is None:
        fd_stencil = stencil

413
414
415
    weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
    u_symp = sp.symbols("u_:{}".format(dimensions))

416
    normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
417
418
419

    gamma = mrt_method.get_equilibrium_terms()
    gamma = gamma.subs({sp.symbols("rho"): 1})
420
    gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)})
Markus Holzer's avatar
Markus Holzer committed
421
    # create the kernels for the initialization of the h field
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):
429
        f.append(weights[i] * ((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness)
430
431
432
                 * scalar_product(d, normal_fd[0:dimensions]))

    for i, _ in enumerate(stencil):
433
        h_updates.append(Assignment(lb_phase_field.center(i), phi_field.center * gamma_init[i] - 0.5 * f[i]))
434
435
436
437

    return h_updates


438
def initializer_kernel_hydro_lb(lb_velocity_field, velocity_field, mrt_method):
439
440
441
    r"""
    Returns an assignment list for initializing the velocity distribution functions
    Args:
442
443
        lb_velocity_field: source field of velocity distribution function
        velocity_field: velocity field
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})
453
    gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)})
454
455
456

    g_updates = list()
    for i, _ in enumerate(stencil):
457
        g_updates.append(Assignment(lb_velocity_field.center(i), gamma_init[i] - weights[i]))
458
459

    return g_updates