multiphase_codegen.py 8.06 KB
Newer Older
Markus Holzer's avatar
Markus Holzer committed
1
2
3
4
5
6
7
from pystencils import fields, TypedSymbol
from pystencils.simp import sympy_cse
from pystencils import AssignmentCollection

from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
from lbmpy.stencils import get_stencil

Markus Holzer's avatar
Markus Holzer committed
8
9
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field
from lbmpy_walberla import generate_lb_pack_info
Markus Holzer's avatar
Markus Holzer committed
10
11
12

from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
    initializer_kernel_hydro_lb, interface_tracking_force, \
Markus Holzer's avatar
Markus Holzer committed
13
    hydrodynamic_force, get_collision_assignments_hydro, get_collision_assignments_phase
Markus Holzer's avatar
Markus Holzer committed
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel

import numpy as np

stencil_phase = get_stencil("D3Q15", "walberla")
stencil_hydro = get_stencil("D3Q27", "walberla")
q_phase = len(stencil_phase)
q_hydro = len(stencil_hydro)

maxwellian_moments = False
assert (len(stencil_phase[0]) == len(stencil_hydro[0]))
dimensions = len(stencil_hydro[0])

########################
# PARAMETER DEFINITION #
########################

density_liquid = 1.0
density_gas = 0.001
surface_tension = 0.0001
M = 0.02

# phase-field parameter
drho3 = (density_liquid - density_gas) / 3
# interface thickness
W = 5
# coefficient related to surface tension
beta = 12.0 * (surface_tension / W)
# coefficient related to surface tension
kappa = 1.5 * surface_tension * W
# relaxation rate allen cahn (h)
w_c = 1.0 / (0.5 + (3.0 * M))

########################
# FIELDS #
########################

# velocity field
u = fields(f"vel_field({dimensions}): [{dimensions}D]", layout='fzyx')
# phase-field
C = fields(f"phase_field: [{dimensions}D]", layout='fzyx')
Markus Holzer's avatar
Markus Holzer committed
56
C_tmp = fields(f"phase_field_tmp: [{dimensions}D]", layout='fzyx')
Markus Holzer's avatar
Markus Holzer committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

# phase-field distribution functions
h = fields(f"lb_phase_field({q_phase}): [{dimensions}D]", layout='fzyx')
h_tmp = fields(f"lb_phase_field_tmp({q_phase}): [{dimensions}D]", layout='fzyx')
# hydrodynamic distribution functions
g = fields(f"lb_velocity_field({q_hydro}): [{dimensions}D]", layout='fzyx')
g_tmp = fields(f"lb_velocity_field_tmp({q_hydro}): [{dimensions}D]", layout='fzyx')

########################################
# RELAXATION RATES AND EXTERNAL FORCES #
########################################

# calculate the relaxation rate for the hydro lb as well as the body force
density = density_gas + C.center * (density_liquid - density_gas)

# force acting on all phases of the model
body_force = np.array([0, 1e-6, 0])

relaxation_time = 0.03 + 0.5
relaxation_rate = 1.0 / relaxation_time

###############
# LBM METHODS #
###############

method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible=True)

method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
                                relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1],
                                maxwellian_moments=maxwellian_moments)

# create the kernels for the initialization of the g and h field
h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)


Markus Holzer's avatar
Markus Holzer committed
93
force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W, fd_stencil=get_stencil("D3Q27"))]
Markus Holzer's avatar
Markus Holzer committed
94
95
force_model_h = MultiphaseForceModel(force=force_h)

Markus Holzer's avatar
Markus Holzer committed
96
97
force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta, body_force,
                             fd_stencil=get_stencil("D3Q27"))
Markus Holzer's avatar
Markus Holzer committed
98

Markus Holzer's avatar
Markus Holzer committed
99
force_model_g = MultiphaseForceModel(force=force_g, rho=density)
Markus Holzer's avatar
Markus Holzer committed
100
101
102
103
104

####################
# LBM UPDATE RULES #
####################

Markus Holzer's avatar
Markus Holzer committed
105
106
107
108
109
110
111
phase_field_LB_step = get_collision_assignments_phase(lb_method=method_phase,
                                                      velocity_input=u,
                                                      output={'density': C_tmp},
                                                      force_model=force_model_h,
                                                      symbolic_fields={"symbolic_field": h,
                                                                       "symbolic_temporary_field": h_tmp},
                                                      kernel_type='stream_pull_collide')
Markus Holzer's avatar
Markus Holzer committed
112
113
114
115
116
117
118
119

phase_field_LB_step = sympy_cse(phase_field_LB_step)

# ---------------------------------------------------------------------------------------------------------

hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
                                                density=density,
                                                velocity_input=u,
Markus Holzer's avatar
Markus Holzer committed
120
121
                                                force_model=force_model_g,
                                                sub_iterations=2,
122
123
                                                symbolic_fields={"symbolic_field": g,
                                                                 "symbolic_temporary_field": g_tmp},
Markus Holzer's avatar
Markus Holzer committed
124
125
126
127
128
129
                                                kernel_type='collide_stream_push')

###################
# GENERATE SWEEPS #
###################

Michael Kuron's avatar
Michael Kuron committed
130
cpu_vec = {'assume_inner_stride_one': True, 'nontemporal': True}
Markus Holzer's avatar
Markus Holzer committed
131
132
133
134
135
136
137
138
139
140

vp = [('int32_t', 'cudaBlockSize0'),
      ('int32_t', 'cudaBlockSize1')]

sweep_block_size = (TypedSymbol("cudaBlockSize0", np.int32),
                    TypedSymbol("cudaBlockSize1", np.int32),
                    1)
sweep_params = {'block_size': sweep_block_size}

info_header = f"""
141
#include "stencil/D3Q{q_phase}.h"\nusing Stencil_phase_T = walberla::stencil::D3Q{q_phase};
Markus Holzer's avatar
Markus Holzer committed
142
143
144
145
#include "stencil/D3Q{q_hydro}.h"\nusing Stencil_hydro_T = walberla::stencil::D3Q{q_hydro};
"""

with CodeGeneration() as ctx:
146
    if not ctx.cuda:
Markus Holzer's avatar
Markus Holzer committed
147
148
149
150
151
152
153
        if not ctx.optimize_for_localhost:
            cpu_vec = {'instruction_set': None}

        generate_sweep(ctx, 'initialize_phase_field_distributions', h_updates)
        generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates)

        generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
Markus Holzer's avatar
Markus Holzer committed
154
                       field_swaps=[(h, h_tmp), (C, C_tmp)],
Markus Holzer's avatar
Markus Holzer committed
155
156
157
158
159
160
161
162
163
                       inner_outer_split=True,
                       cpu_vectorize_info=cpu_vec)

        generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
                       field_swaps=[(g, g_tmp)],
                       inner_outer_split=True,
                       cpu_vectorize_info=cpu_vec)

        # communication
Markus Holzer's avatar
Markus Holzer committed
164
165
166
167
168
169
170
        generate_lb_pack_info(ctx, 'PackInfo_phase_field_distributions', stencil_phase, h,
                              streaming_pattern='pull', target='cpu')

        generate_lb_pack_info(ctx, 'PackInfo_velocity_based_distributions', stencil_hydro, g,
                              streaming_pattern='push', target='cpu')

        generate_pack_info_for_field(ctx, 'PackInfo_phase_field', C, target='cpu')
Markus Holzer's avatar
Markus Holzer committed
171
172
173
174
175
176
177
178
179
180

        ctx.write_file("GenDefines.h", info_header)

    if ctx.cuda:
        generate_sweep(ctx, 'initialize_phase_field_distributions',
                       h_updates, target='gpu')
        generate_sweep(ctx, 'initialize_velocity_based_distributions',
                       g_updates, target='gpu')

        generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
Markus Holzer's avatar
Markus Holzer committed
181
                       field_swaps=[(h, h_tmp), (C, C_tmp)],
Markus Holzer's avatar
Markus Holzer committed
182
183
184
185
186
187
188
189
190
191
192
193
                       inner_outer_split=True,
                       target='gpu',
                       gpu_indexing_params=sweep_params,
                       varying_parameters=vp)

        generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
                       field_swaps=[(g, g_tmp)],
                       inner_outer_split=True,
                       target='gpu',
                       gpu_indexing_params=sweep_params,
                       varying_parameters=vp)
        # communication
Markus Holzer's avatar
Markus Holzer committed
194
195
196
197
198
199
200
        generate_lb_pack_info(ctx, 'PackInfo_phase_field_distributions', stencil_phase, h,
                              streaming_pattern='pull', target='gpu')

        generate_lb_pack_info(ctx, 'PackInfo_velocity_based_distributions', stencil_hydro, g,
                              streaming_pattern='push', target='gpu')

        generate_pack_info_for_field(ctx, 'PackInfo_phase_field', C, target='gpu')
Markus Holzer's avatar
Markus Holzer committed
201
202
203
204

        ctx.write_file("GenDefines.h", info_header)

print("finished code generation successfully")