Skip to content

Commit 98d0887

Browse files
wilfonbawilfonba
committed
simplified transfer lists and new lagrange advection example
1 parent 151b511 commit 98d0887

2 files changed

Lines changed: 241 additions & 116 deletions

File tree

‎examples/3D_lagrange_advection/case.py‎

Lines changed: 224 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,224 @@
1+
#!/usr/bin/env python3
2+
# 3D tri-periodic advection of a cubic lattice of Lagrangian bubbles
3+
# with log-normally distributed radii. Background flow is uniform along
4+
# the (1,1,1) diagonal.
5+
import json
6+
import math
7+
import os
8+
9+
import numpy as np
10+
11+
# Reference values for nondimensionalization
12+
L0 = 1e-3 # length - m
13+
rho0 = 950 # density - kg/m3
14+
c0 = 1449.0 # speed of sound - m/s
15+
p0 = rho0 * c0 * c0 # pressure - Pa
16+
T0 = 298 # temperature - K
17+
18+
# Domain: cube [-1, 1] mm, nondimensionalized to [-1, 1]
19+
xb, xe = -1.0e-3 / L0, 1.0e-3 / L0
20+
yb, ye = -1.0e-3 / L0, 1.0e-3 / L0
21+
zb, ze = -1.0e-3 / L0, 1.0e-3 / L0
22+
23+
# Host (water)
24+
gamma_host = 6.12
25+
pi_inf_host = 3.43e8
26+
mu_host = 0.001
27+
rho_host = 950
28+
29+
# Lagrangian bubbles' properties (matched to the reference 3D lag case)
30+
R_uni = 8314
31+
MW_g = 28.0
32+
MW_v = 18.0
33+
gamma_g = 1.4
34+
gamma_v = 1.333
35+
pv = 2350
36+
cp_g = 1.0e3
37+
cp_v = 2.1e3
38+
k_g = 0.025
39+
k_v = 0.02
40+
diffVapor = 2.5e-5
41+
sigBubble = 0.020
42+
mu_g = 1.48e-5
43+
rho_g = 1
44+
45+
patm = 1e5
46+
rho_host_nd = rho_host / rho0
47+
rho_g_nd = rho_g / rho0
48+
pres = patm / p0
49+
50+
# Diagonal advection: equal (nondim) velocity components. ~Mach 0.01.
51+
u_adv = 0.05
52+
v_adv = 0.05
53+
w_adv = 0.05
54+
55+
# Time stepping. With m=n=p=49 (dx=0.04) and unit nondim sound speed,
56+
# CFL=0.5 at dt=0.02. One full traversal of the domain along x at
57+
# u_adv=0.01 takes (xe-xb)/u_adv = 200 nondim time units.
58+
Ncells = 50
59+
m = n = p = Ncells - 1
60+
dx = (xe - xb) / Ncells
61+
dt = 0.02
62+
t_traverse = (xe - xb) / u_adv
63+
t_step_stop = int(round(t_traverse / dt))
64+
t_step_save = t_step_stop // 100
65+
66+
# Lattice and log-normal radius distribution
67+
N_side = 2
68+
nBubs = N_side**3
69+
R_median_nd = 0.02 # 20 micron median radius
70+
sigma_ln = 0.3
71+
72+
data = {
73+
"run_time_info": "T",
74+
"x_domain%beg": xb,
75+
"x_domain%end": xe,
76+
"y_domain%beg": yb,
77+
"y_domain%end": ye,
78+
"z_domain%beg": zb,
79+
"z_domain%end": ze,
80+
"m": m,
81+
"n": n,
82+
"p": p,
83+
"dt": dt,
84+
"t_step_start": 0,
85+
"t_step_stop": t_step_stop,
86+
"t_step_save": t_step_save,
87+
"model_eqns": 2,
88+
"alt_soundspeed": "F",
89+
"mixture_err": "F",
90+
"mpp_lim": "T",
91+
"time_stepper": 3,
92+
"weno_order": 5,
93+
"mapped_weno": "T",
94+
"mp_weno": "F",
95+
"avg_state": 2,
96+
"weno_eps": 1e-16,
97+
"riemann_solver": 2,
98+
"wave_speeds": 1,
99+
# Tri-periodic boundaries
100+
"bc_x%beg": -1,
101+
"bc_x%end": -1,
102+
"bc_y%beg": -1,
103+
"bc_y%end": -1,
104+
"bc_z%beg": -1,
105+
"bc_z%end": -1,
106+
"num_patches": 1,
107+
"num_fluids": 2,
108+
"format": 1,
109+
"precision": 2,
110+
"prim_vars_wrt": "T",
111+
"parallel_io": "T",
112+
"lag_db_wrt": "T",
113+
"lag_rad_wrt": "T",
114+
# Fluid 1: water host
115+
"fluid_pp(1)%gamma": 1.0 / (gamma_host - 1.0),
116+
"fluid_pp(1)%pi_inf": gamma_host * (pi_inf_host / p0) / (gamma_host - 1.0),
117+
# Fluid 2: gas (inside bubbles)
118+
"fluid_pp(2)%gamma": 1.0 / (gamma_g - 1.0),
119+
"fluid_pp(2)%pi_inf": 0.0,
120+
# Bubble physical parameters
121+
"bub_pp%R0ref": 1.0,
122+
"bub_pp%p0ref": 1.0,
123+
"bub_pp%rho0ref": 1.0,
124+
"bub_pp%T0ref": 1.0,
125+
"bub_pp%ss": sigBubble / (rho0 * L0 * c0 * c0),
126+
"bub_pp%pv": pv / p0,
127+
"bub_pp%vd": diffVapor / (L0 * c0),
128+
"bub_pp%mu_l": mu_host / (rho0 * L0 * c0),
129+
"bub_pp%gam_v": gamma_v,
130+
"bub_pp%gam_g": gamma_g,
131+
"bub_pp%M_v": MW_v,
132+
"bub_pp%M_g": MW_g,
133+
"bub_pp%k_v": k_v * (T0 / (L0 * rho0 * c0 * c0 * c0)),
134+
"bub_pp%k_g": k_g * (T0 / (L0 * rho0 * c0 * c0 * c0)),
135+
"bub_pp%cp_v": cp_v * (T0 / (c0 * c0)),
136+
"bub_pp%cp_g": cp_g * (T0 / (c0 * c0)),
137+
"bub_pp%R_v": (R_uni / MW_v) * (T0 / (c0 * c0)),
138+
"bub_pp%R_g": (R_uni / MW_g) * (T0 / (c0 * c0)),
139+
# Viscosity
140+
"viscous": "T",
141+
"fluid_pp(1)%Re(1)": 1.0 / (mu_host / (rho0 * c0 * L0)),
142+
"fluid_pp(2)%Re(1)": 1.0 / (mu_g / (rho0 * c0 * L0)),
143+
# Single patch: uniform background flow advecting diagonally
144+
"patch_icpp(1)%geometry": 9,
145+
"patch_icpp(1)%x_centroid": (xb + xe) / 2,
146+
"patch_icpp(1)%y_centroid": (yb + ye) / 2,
147+
"patch_icpp(1)%z_centroid": (zb + ze) / 2,
148+
"patch_icpp(1)%length_x": (xe - xb),
149+
"patch_icpp(1)%length_y": (ye - yb),
150+
"patch_icpp(1)%length_z": (ze - zb),
151+
"patch_icpp(1)%vel(1)": u_adv,
152+
"patch_icpp(1)%vel(2)": v_adv,
153+
"patch_icpp(1)%vel(3)": w_adv,
154+
"patch_icpp(1)%pres": pres,
155+
"patch_icpp(1)%alpha_rho(1)": rho_host_nd,
156+
"patch_icpp(1)%alpha_rho(2)": 0,
157+
"patch_icpp(1)%alpha(1)": 1,
158+
"patch_icpp(1)%alpha(2)": 0,
159+
# Lagrangian bubbles
160+
"bubbles_lagrange": "T",
161+
"fd_order": 4,
162+
"bubble_model": 3,
163+
"thermal": 3,
164+
"polytropic": "F",
165+
"adap_dt": "T",
166+
"lag_params%nBubs_glb": nBubs,
167+
"lag_params%vel_model": 1,
168+
"lag_params%drag_model": 1,
169+
"lag_params%solver_approach": 1,
170+
"lag_params%cluster_type": 2,
171+
"lag_params%pressure_corrector": "F",
172+
"lag_params%smooth_type": 1,
173+
"lag_params%epsilonb": 1.0,
174+
"lag_params%valmaxvoid": 0.9,
175+
"lag_params%write_bubbles": "T",
176+
"lag_params%write_bubbles_stats": "T",
177+
"lag_params%write_void_evol": "T",
178+
}
179+
180+
181+
def write_lag_bubbles_file():
182+
# Cubic lattice spanning the periodic domain. Lattice spacing equals
183+
# the domain length divided by N_side; positions are shifted by
184+
# spacing/2 so that no bubble sits exactly on a periodic face.
185+
Lx = xe - xb
186+
Ly = ye - yb
187+
Lz = ze - zb
188+
sx = Lx / N_side
189+
sy = Ly / N_side
190+
sz = Lz / N_side
191+
192+
xs = xb + 0.5 * sx + np.arange(N_side) * sx
193+
ys = yb + 0.5 * sy + np.arange(N_side) * sy
194+
zs = zb + 0.5 * sz + np.arange(N_side) * sz
195+
X, Y, Z = np.meshgrid(xs, ys, zs, indexing="ij")
196+
x_pos = X.flatten()
197+
y_pos = Y.flatten()
198+
z_pos = Z.flatten()
199+
200+
# Log-normal radii (nondim by L0). Clip the tails so a freak draw
201+
# cannot overlap with neighbors on the lattice.
202+
rng = np.random.default_rng(seed=42)
203+
R = rng.lognormal(mean=math.log(R_median_nd), sigma=sigma_ln, size=nBubs)
204+
R_max_allowed = 0.45 * min(sx, sy, sz)
205+
R = np.minimum(R, R_max_allowed)
206+
207+
vx = 0.0 * np.full(nBubs, u_adv)
208+
vy = 0.0 * np.full(nBubs, v_adv)
209+
vz = 0.0 * np.full(nBubs, w_adv)
210+
# Small initial radial velocity to kick off bubble dynamics.
211+
# For R~0.02 at atmospheric pres, Minnaert ω_n ≈ 0.73 (nondim),
212+
# so rdot=1e-3 → ε/R₀ ≈ rdot/(R₀·ω_n) ≈ 7% amplitude oscillations.
213+
rdot = 0.0 * np.full(nBubs, 1.0e-3)
214+
215+
input_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), "input")
216+
os.makedirs(input_dir, exist_ok=True)
217+
with open(os.path.join(input_dir, "lag_bubbles.dat"), "w") as f:
218+
for i in range(nBubs):
219+
f.write(f"{x_pos[i]:12.6f}\t{y_pos[i]:12.6f}\t{z_pos[i]:12.6f}\t{vx[i]:12.6f}\t{vy[i]:12.6f}\t{vz[i]:12.6f}\t{R[i]:12.6f}\t{rdot[i]:12.6f}\n")
220+
221+
222+
write_lag_bubbles_file()
223+
224+
print(json.dumps(data, indent=4))

‎src/simulation/m_mpi_proxy.fpp‎

Lines changed: 17 additions & 116 deletions
Original file line numberDiff line numberDiff line change
@@ -304,6 +304,7 @@ contains
304304
real(wp), dimension(:,:), intent(in) :: pos, posPrev
305305
integer :: bubID
306306
integer :: i, j, k
307+
integer :: dx, dy, dz
307308

308309
do k = nidx(3)%beg, nidx(3)%end
309310
do j = nidx(2)%beg, nidx(2)%end
@@ -314,129 +315,29 @@ contains
314315
end do
315316

316317
do k = 1, nbub
318+
dx = 0; dy = 0; dz = 0
317319
if (f_crosses_boundary(k, 1, -1, pos, posPrev)) then
318-
call s_add_particle_to_direction(k, -1, 0, 0)
319-
if (n > 0) then
320-
if (f_crosses_boundary(k, 2, -1, pos, posPrev)) then
321-
call s_add_particle_to_direction(k, -1, -1, 0)
322-
call s_add_particle_to_direction(k, 0, -1, 0)
323-
if (p > 0) then
324-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
325-
call s_add_particle_to_direction(k, -1, -1, -1)
326-
call s_add_particle_to_direction(k, 0, -1, -1)
327-
call s_add_particle_to_direction(k, -1, 0, -1)
328-
call s_add_particle_to_direction(k, 0, 0, -1)
329-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
330-
call s_add_particle_to_direction(k, -1, -1, 1)
331-
call s_add_particle_to_direction(k, 0, -1, 1)
332-
call s_add_particle_to_direction(k, -1, 0, 1)
333-
call s_add_particle_to_direction(k, 0, 0, 1)
334-
end if
335-
end if
336-
else if (f_crosses_boundary(k, 2, 1, pos, posPrev)) then
337-
call s_add_particle_to_direction(k, -1, 1, 0)
338-
call s_add_particle_to_direction(k, 0, 1, 0)
339-
if (p > 0) then
340-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
341-
call s_add_particle_to_direction(k, -1, 1, -1)
342-
call s_add_particle_to_direction(k, 0, 1, -1)
343-
call s_add_particle_to_direction(k, -1, 0, -1)
344-
call s_add_particle_to_direction(k, 0, 0, -1)
345-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
346-
call s_add_particle_to_direction(k, -1, 1, 1)
347-
call s_add_particle_to_direction(k, 0, 1, 1)
348-
call s_add_particle_to_direction(k, -1, 0, 1)
349-
call s_add_particle_to_direction(k, 0, 0, 1)
350-
end if
351-
end if
352-
else
353-
if (p > 0) then
354-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
355-
call s_add_particle_to_direction(k, -1, 0, -1)
356-
call s_add_particle_to_direction(k, 0, 0, -1)
357-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
358-
call s_add_particle_to_direction(k, -1, 0, 1)
359-
call s_add_particle_to_direction(k, 0, 0, 1)
360-
end if
361-
end if
362-
end if
363-
end if
320+
dx = -1
364321
else if (f_crosses_boundary(k, 1, 1, pos, posPrev)) then
365-
call s_add_particle_to_direction(k, 1, 0, 0)
366-
if (n > 0) then
367-
if (f_crosses_boundary(k, 2, -1, pos, posPrev)) then
368-
call s_add_particle_to_direction(k, 1, -1, 0)
369-
call s_add_particle_to_direction(k, 0, -1, 0)
370-
if (p > 0) then
371-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
372-
call s_add_particle_to_direction(k, 1, -1, -1)
373-
call s_add_particle_to_direction(k, 0, -1, -1)
374-
call s_add_particle_to_direction(k, 1, 0, -1)
375-
call s_add_particle_to_direction(k, 0, 0, -1)
376-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
377-
call s_add_particle_to_direction(k, 1, -1, 1)
378-
call s_add_particle_to_direction(k, 0, -1, 1)
379-
call s_add_particle_to_direction(k, 1, 0, 1)
380-
call s_add_particle_to_direction(k, 0, 0, 1)
381-
end if
382-
end if
383-
else if (f_crosses_boundary(k, 2, 1, pos, posPrev)) then
384-
call s_add_particle_to_direction(k, 1, 1, 0)
385-
call s_add_particle_to_direction(k, 0, 1, 0)
386-
if (p > 0) then
387-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
388-
call s_add_particle_to_direction(k, 1, 1, -1)
389-
call s_add_particle_to_direction(k, 0, 1, -1)
390-
call s_add_particle_to_direction(k, 1, 0, -1)
391-
call s_add_particle_to_direction(k, 0, 0, -1)
392-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
393-
call s_add_particle_to_direction(k, 1, 1, 1)
394-
call s_add_particle_to_direction(k, 0, 1, 1)
395-
call s_add_particle_to_direction(k, 1, 0, 1)
396-
call s_add_particle_to_direction(k, 0, 0, 1)
397-
end if
398-
end if
399-
else
400-
if (p > 0) then
401-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
402-
call s_add_particle_to_direction(k, 1, 0, -1)
403-
call s_add_particle_to_direction(k, 0, 0, -1)
404-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
405-
call s_add_particle_to_direction(k, 1, 0, 1)
406-
call s_add_particle_to_direction(k, 0, 0, 1)
407-
end if
408-
end if
409-
end if
410-
end if
411-
else if (f_crosses_boundary(k, 2, -1, pos, posPrev)) then
412-
call s_add_particle_to_direction(k, 0, -1, 0)
413-
if (p > 0) then
414-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
415-
call s_add_particle_to_direction(k, 0, -1, -1)
416-
call s_add_particle_to_direction(k, 0, 0, -1)
417-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
418-
call s_add_particle_to_direction(k, 0, -1, 1)
419-
call s_add_particle_to_direction(k, 0, 0, 1)
420-
end if
421-
end if
422-
else if (f_crosses_boundary(k, 2, 1, pos, posPrev)) then
423-
call s_add_particle_to_direction(k, 0, 1, 0)
424-
if (p > 0) then
425-
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
426-
call s_add_particle_to_direction(k, 0, 1, -1)
427-
call s_add_particle_to_direction(k, 0, 0, -1)
428-
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
429-
call s_add_particle_to_direction(k, 0, 1, 1)
430-
call s_add_particle_to_direction(k, 0, 0, 1)
431-
end if
322+
dx = 1
323+
end if
324+
if (n > 0) then
325+
if (f_crosses_boundary(k, 2, -1, pos, posPrev)) then
326+
dy = -1
327+
else if (f_crosses_boundary(k, 2, 1, pos, posPrev)) then
328+
dy = 1
432329
end if
433-
else if (p > 0) then
330+
end if
331+
if (p > 0) then
434332
if (f_crosses_boundary(k, 3, -1, pos, posPrev)) then
435-
call s_add_particle_to_direction(k, 0, 0, -1)
333+
dz = -1
436334
else if (f_crosses_boundary(k, 3, 1, pos, posPrev)) then
437-
call s_add_particle_to_direction(k, 0, 0, 1)
335+
dz = 1
438336
end if
439337
end if
338+
if (abs(dx) + abs(dy) + abs(dz) /= 0) then
339+
call s_add_particle_to_direction(k, dx, dy, dz)
340+
end if
440341
end do
441342

442343
contains

0 commit comments

Comments
 (0)