kernel_equations.py 18.8 KB
Newer Older
Markus Holzer's avatar
Markus Holzer committed
1
2
from lbmpy.creationfunctions import update_with_default_parameters
from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAccessor
3
4
5
6
7
8
9
10
11
12
13
14
15
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
16
        stencil: stencil to derive the finite difference for the laplacian (2nd order isotropic)
17
18
19
20
        beta: coefficient related to surface tension and interface thickness
        kappa: coefficient related to surface tension and interface thickness
    """
    dimensions = len(stencil[0])
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
    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)
44
45
46

    # 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
47
        kappa * lap
48
49
50
51
52
53
54
55
    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
56
        stencil: stencil to derive the finite difference for the gradient (2nd order isotropic)
57
58
    """
    dimensions = len(stencil[0])
59
    q = len(stencil)
Markus Holzer's avatar
Markus Holzer committed
60
    deriv = FiniteDifferenceStencilDerivation((0,), stencil)
61
62
63
64
65
66

    deriv.assume_symmetric(0, anti_symmetric=True)
    deriv.assume_symmetric(1, anti_symmetric=False)
    if dimensions == 3:
        deriv.assume_symmetric(2, anti_symmetric=False)

67
68
69
    # 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
70
71
        res = deriv.get_stencil(isotropic=True)
        grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
72
73
74
75
76
77
78
    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))
79
80
        deriv.set_weight((1, 0, 0), sp.Rational(1, 6))

Markus Holzer's avatar
Markus Holzer committed
81
        res = deriv.get_stencil(isotropic=True)
82
        grad = [res.apply(phi_field.center),
83
84
                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
85
    else:
86
87
88
89
90
91
92
        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))]
93
94
95
96

    return grad


97
def normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil=None):
98
99
100
101
102
    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
103
104
        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.
105
    """
106
107
    if fd_stencil is None:
        fd_stencil = stencil
108

109
110
111
    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)]
112
113
114
    return result


115
def pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil=None):
116
117
118
119
120
121
122
    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
123
124
        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.
125
    """
126
127
128
129
    if fd_stencil is None:
        fd_stencil = stencil

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


134
def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light, fd_stencil=None):
135
136
137
138
139
140
141
142
143
    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
144
145
        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.
146
147
148
149
    """
    stencil = mrt_method.stencil
    dimensions = len(stencil[0])

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

    iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
154
155

    moment_matrix = mrt_method.moment_matrix
Markus Holzer's avatar
Markus Holzer committed
156
157
158
159
160
    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)]
161
    m0 = np.dot(moment_matrix.tolist(), g_vals)
162

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

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

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

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

    density_difference = density_heavy - density_light

    # Calculate Viscous Force MRT
Markus Holzer's avatar
Markus Holzer committed
183
184
185
    fmx = (0.5 - tau) * (stress_tensor[0] * iso_grad[0]
                         + stress_tensor[5] * iso_grad[1]
                         + stress_tensor[4] * iso_grad[2]) * density_difference
186

Markus Holzer's avatar
Markus Holzer committed
187
188
189
    fmy = (0.5 - tau) * (stress_tensor[5] * iso_grad[0]
                         + stress_tensor[1] * iso_grad[1]
                         + stress_tensor[3] * iso_grad[2]) * density_difference
190

Markus Holzer's avatar
Markus Holzer committed
191
192
193
    fmz = (0.5 - tau) * (stress_tensor[4] * iso_grad[0]
                         + stress_tensor[3] * iso_grad[1]
                         + stress_tensor[2] * iso_grad[2]) * density_difference
194
195
196
197

    return [fmx, fmy, fmz]


198
def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None):
199
200
201
202
203
204
205
    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
206
207
        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.
208
    """
209
210
211
212
213
    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
214
    return [chemical_potential * x for x in iso_grad]
215
216


Markus Holzer's avatar
Markus Holzer committed
217
def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
218
                       density_heavy, density_light, kappa, beta, body_force, fd_stencil=None):
219
220
221
222
223
    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
224
        lb_method: Lattice boltzmann method used for hydrodynamics
225
226
227
228
229
230
        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
231
232
        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.
233
    """
Markus Holzer's avatar
Markus Holzer committed
234
    stencil = lb_method.stencil
235
    dimensions = len(stencil[0])
236
237
238
239
240
241
242

    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)
243
244
245
246
247
248
249
250

    result = []
    for i in range(dimensions):
        result.append(fs[i] + fp[i] + fm[i] + body_force[i])

    return result


251
def interface_tracking_force(phi_field, stencil, interface_thickness, fd_stencil=None):
252
253
254
255
256
257
    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
258
259
        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.
260
    """
261
262
263
    if fd_stencil is None:
        fd_stencil = stencil

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

    return result


Markus Holzer's avatar
Markus Holzer committed
273
def get_update_rules_velocity(src_field, u_in, lb_method, force, density, sub_iterations=2):
274
    r"""
Markus Holzer's avatar
Markus Holzer committed
275
276
277
278
279
280
281
     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
Markus Holzer's avatar
Markus Holzer committed
282
         sub_iterations: number of updates of the velocity field
Markus Holzer's avatar
Markus Holzer committed
283
284
     """
    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
    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)]
296
    m0 = np.dot(moment_matrix.tolist(), src)
Markus Holzer's avatar
Markus Holzer committed
297
298
299

    update_u = list()
    update_u.append(Assignment(sp.symbols("rho"), m0[0]))
300

Markus Holzer's avatar
Markus Holzer committed
301
    index = 0
302
    u_symp = sp.symbols("u_:{}".format(dimensions))
Markus Holzer's avatar
Markus Holzer committed
303
304
    aleph = sp.symbols("aleph_:{}".format(dimensions * sub_iterations))

Markus Holzer's avatar
Markus Holzer committed
305
    for i in range(dimensions):
Markus Holzer's avatar
Markus Holzer committed
306
307
308
309
310
311
312
313
        update_u.append(Assignment(aleph[i], u_in.center_vector[i]))
        index += 1

    for k in range(sub_iterations - 1):
        subs_dict = dict(zip(u_symp, aleph[k * dimensions:index]))
        for i in range(dimensions):
            update_u.append(Assignment(aleph[index], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
            index += 1
Markus Holzer's avatar
Markus Holzer committed
314

Markus Holzer's avatar
Markus Holzer committed
315
    subs_dict = dict(zip(u_symp, aleph[index - dimensions:index]))
Markus Holzer's avatar
Markus Holzer committed
316
    for i in range(dimensions):
Markus Holzer's avatar
Markus Holzer committed
317
318
        update_u.append(Assignment(u_symp[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
        # update_u.append(Assignment(u_in.center_vector[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
Markus Holzer's avatar
Markus Holzer committed
319
320
321

    return update_u

322

323
324
def get_collision_assignments_hydro(lb_method, density, velocity_input, force, sub_iterations, symbolic_fields,
                                    kernel_type):
Markus Holzer's avatar
Markus Holzer committed
325
326
327
328
    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:
329
         lb_method: moment based lattice Boltzmann method
Markus Holzer's avatar
Markus Holzer committed
330
         density: the interpolated density of the simulation
331
332
         velocity_input: velocity field for the hydrodynamic and Allen-Chan LB step
         force: force vector containing a summation of the surface tension-, pressure-, viscous- and bodyforce vector
Markus Holzer's avatar
Markus Holzer committed
333
         sub_iterations: number of updates of the velocity field
334
335
         symbolic_fields: PDF fields for source and destination
         kernel_type: collide_stream_push or collide_only
Markus Holzer's avatar
Markus Holzer committed
336
337
338
339
340
     """

    stencil = lb_method.stencil
    dimensions = len(stencil[0])

341
342
    src_field = symbolic_fields['symbolic_field']
    dst_field = symbolic_fields['symbolic_temporary_field']
Markus Holzer's avatar
Markus Holzer committed
343
344
345
346
347
348
349

    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()
350
    for i in range(dimensions):
Markus Holzer's avatar
Markus Holzer committed
351
352
353
        indices.append(eq.index(first_eqs[i]))

    eq = np.array(eq)
354

Markus Holzer's avatar
Markus Holzer committed
355
    g_vals = [src_field.center(i) for i, _ in enumerate(stencil)]
356
    m0 = np.dot(moment_matrix.tolist(), g_vals)
357

Markus Holzer's avatar
Markus Holzer committed
358
    mf = np.zeros(len(stencil), dtype=object)
359
    for i in range(dimensions):
Markus Holzer's avatar
Markus Holzer committed
360
361
362
363
        mf[indices[i]] = force[i] / density

    m = sp.symbols("m_:{}".format(len(stencil)))

364
365
    update_m = get_update_rules_velocity(src_field, velocity_input, lb_method, force,
                                         density, sub_iterations=sub_iterations)
Markus Holzer's avatar
Markus Holzer committed
366
367
368
369
370
371
    u_symp = sp.symbols("u_:{}".format(dimensions))

    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()
372
    var = np.dot(moment_matrix.inv().tolist(), m)
373
    if kernel_type == 'collide_stream_push':
Markus Holzer's avatar
Markus Holzer committed
374
375
376
377
378
        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)
379

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

Markus Holzer's avatar
Markus Holzer committed
383
    for i in range(dimensions):
384
        update_g.append(Assignment(velocity_input.center_vector[i], u_symp[i]))
Markus Holzer's avatar
Markus Holzer committed
385

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

Markus Holzer's avatar
Markus Holzer committed
389
    return hydro_lb_update_rule
390
391


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

    if fd_stencil is None:
        fd_stencil = stencil

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

414
    normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
415
416
417

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

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

    return h_updates


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

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

    return g_updates