test_codegen_3d.py 8.11 KB
Newer Older
1
2
3
from collections import OrderedDict

import numpy as np
Markus Holzer's avatar
Markus Holzer committed
4
5
6

from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
7
from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
Markus Holzer's avatar
Markus Holzer committed
8
from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
9
10
11
12
13
14
15
from lbmpy.phasefield_allen_cahn.kernel_equations import (
    get_collision_assignments_hydro, hydrodynamic_force, initializer_kernel_hydro_lb, initializer_kernel_phase_field_lb,
    interface_tracking_force)
from lbmpy.phasefield_allen_cahn.parameter_calculation import (
    calculate_dimensionless_rising_bubble, calculate_parameters_rti)
from lbmpy.stencils import get_stencil
from pystencils import AssignmentCollection, fields
Markus Holzer's avatar
Markus Holzer committed
16
17


18
def test_codegen_3d():
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
56
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
93
94
95
96
97
98
99
100
    stencil_phase = get_stencil("D3Q15")
    stencil_hydro = get_stencil("D3Q27")
    assert (len(stencil_phase[0]) == len(stencil_hydro[0]))
    dimensions = len(stencil_hydro[0])

    parameters = calculate_dimensionless_rising_bubble(reference_time=18000,
                                                       density_heavy=1.0,
                                                       bubble_radius=16,
                                                       bond_number=30,
                                                       reynolds_number=420,
                                                       density_ratio=1000,
                                                       viscosity_ratio=100)

    np.isclose(parameters["density_light"], 0.001, rtol=1e-05, atol=1e-08, equal_nan=False)
    np.isclose(parameters["gravitational_acceleration"], -9.876543209876543e-08, rtol=1e-05, atol=1e-08, equal_nan=False)

    parameters = calculate_parameters_rti(reference_length=128,
                                          reference_time=18000,
                                          density_heavy=1.0,
                                          capillary_number=9.1,
                                          reynolds_number=128,
                                          atwood_number=1.0,
                                          peclet_number=744,
                                          density_ratio=3,
                                          viscosity_ratio=3)

    np.isclose(parameters["density_light"], 1/3, rtol=1e-05, atol=1e-08, equal_nan=False)
    np.isclose(parameters["gravitational_acceleration"], -3.9506172839506174e-07, rtol=1e-05, atol=1e-08, equal_nan=False)
    np.isclose(parameters["mobility"], 0.0012234169653524492, rtol=1e-05, atol=1e-08, equal_nan=False)

    rs = analytic_rising_speed(1-6, 20, 0.01)
    np.isclose(rs, 16666.666666666668, rtol=1e-05, atol=1e-08, equal_nan=False)

    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
    u = fields("vel_field(" + str(dimensions) + "): [" + str(dimensions) + "D]", layout='fzyx')
    C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')

    h = fields("lb_phase_field(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')
    h_tmp = fields("lb_phase_field_tmp(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')

    g = fields("lb_velocity_field(" + str(len(stencil_hydro)) + "): [" + str(dimensions) + "D]", layout='fzyx')
    g_tmp = fields("lb_velocity_field_tmp(" + str(len(stencil_hydro)) + "): [" + str(dimensions) + "D]", layout='fzyx')

    # 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

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

    mrt = create_lb_method(method="mrt", weighted=False, stencil=stencil_hydro,
                           relaxation_rates=[1, 1, relaxation_rate, 1, 1, 1, 1])
    rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))

    method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=False)

    # 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)

    force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
    force_model_h = MultiphaseForceModel(force=force_h)

101
102
    force_g = hydrodynamic_force(g, C, method_hydro,
                                 relaxation_time, density_liquid, density_gas, kappa, beta, body_force)
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    force_model_g = MultiphaseForceModel(force=force_g, rho=density)

    h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
    sum_h = np.sum(h_tmp_symbol_list[:])

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

    allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,
                                          velocity_input=u,
                                          compressible=True,
                                          optimization={"symbolic_field": h,
                                                        "symbolic_temporary_field": h_tmp},
                                          kernel_type='stream_pull_collide')

    allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})
    allen_cahn_update_rule = AssignmentCollection(main_assignments=allen_cahn_lb.main_assignments,
                                                  subexpressions=allen_cahn_lb.subexpressions)
    # ---------------------------------------------------------------------------------------------------------

    method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
    hydro_lb_update_rule_normal = get_collision_assignments_hydro(lb_method=method_hydro,
                                                                  density=density,
                                                                  velocity_input=u,
                                                                  force=force_g,
                                                                  optimization={"symbolic_field": g,
                                                                                "symbolic_temporary_field": g_tmp},
                                                                  kernel_type='collide_only')

    hydro_lb_update_rule_push = get_collision_assignments_hydro(lb_method=method_hydro,
                                                                density=density,
                                                                velocity_input=u,
                                                                force=force_g,
                                                                optimization={"symbolic_field": g,
                                                                              "symbolic_temporary_field": g_tmp},
                                                                kernel_type='collide_stream_push')

    hydro_lb_update_rule_generic_fields = get_collision_assignments_hydro(lb_method=method_hydro,
                                                                          density=density,
                                                                          velocity_input=u,
                                                                          force=force_g,
                                                                          kernel_type='collide_only')