Skip to content

Commit 3d7c11d

Browse files
sbryngelsonsbryngelson
committed
simulation: deterministic RNG + solenoidal polarization + mixture density in synthetic forcing
1 parent 5cfd0a9 commit 3d7c11d

1 file changed

Lines changed: 38 additions & 42 deletions

File tree

‎src/simulation/m_body_forces.fpp‎

Lines changed: 38 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ module m_body_forces
1111
use m_global_parameters
1212
use m_variables_conversion
1313
use m_mpi_proxy
14+
use m_helper, only: s_prng, f_unit_vector, f_cross
1415
use m_nvtx
1516

1617
! $:USE_GPU_MODULE()
@@ -44,10 +45,10 @@ contains
4445
!> broadcasts to all MPI ranks and copies to GPU.
4546
impure subroutine s_initialize_body_forces_module
4647

47-
integer :: s, m_wave, m_global
48-
integer :: n_seed
49-
integer, allocatable :: seed_arr(:)
50-
real(wp) :: rn1, rn2, kz, r_xy, k_mag
48+
integer :: s, m_wave, m_global
49+
integer :: seed
50+
real(wp) :: rn1, rn2, k_mag
51+
real(wp), dimension(3) :: khat, xi, sig, sig_tmp
5152

5253
if (n > 0) then
5354
if (p > 0) then
@@ -81,16 +82,12 @@ contains
8182
@:ALLOCATE(synthetic_k_z(1:num_synthetic_wave_numbers))
8283
end if
8384

84-
! Generate random wave vectors and phases on rank 0, then broadcast
85+
! Generate random wave vectors and phases on rank 0, then broadcast. Uses the
86+
! compiler-independent LCG (s_prng) so the forcing is reproducible across
87+
! compilers; the 3-D polarization is built perpendicular to k via a double
88+
! cross product, guaranteeing a divergence-free (solenoidal) mode.
8589
if (proc_rank == 0) then
86-
call random_seed(size=n_seed)
87-
allocate (seed_arr(n_seed))
88-
! Vary each element so the seed state is not degenerate
89-
do s = 1, n_seed
90-
seed_arr(s) = synth_seed + (s - 1)*1000003
91-
end do
92-
call random_seed(put=seed_arr)
93-
deallocate (seed_arr)
90+
seed = synth_seed
9491

9592
m_global = 0
9693
do s = 1, synth_n_shells
@@ -105,41 +102,36 @@ contains
105102
synthetic_ey(m_global) = 0._wp
106103
synthetic_ez(m_global) = 0._wp
107104
else if (num_dims == 2) then
108-
! Azimuthal angle theta uniform on [0, 2*pi)
109-
call random_number(rn1)
105+
! In-plane wavevector at azimuth theta; solenoidal dir perpendicular to k
106+
call s_prng(rn1, seed)
110107
rn1 = rn1*2._wp*pi
111108
synthetic_k_x(m_global) = k_mag*cos(rn1)
112109
synthetic_k_y(m_global) = k_mag*sin(rn1)
113-
! Solenoidal direction perpendicular to k: (-sin theta, cos theta)
114110
synthetic_ex(m_global) = -sin(rn1)
115111
synthetic_ey(m_global) = cos(rn1)
116112
synthetic_ez(m_global) = 0._wp
117113
else
118-
! Uniform on sphere: azimuthal phi, cos(polar) uniform in [-1,1]
119-
call random_number(rn1)
120-
call random_number(rn2)
121-
rn1 = rn1*2._wp*pi
122-
kz = 2._wp*rn2 - 1._wp
123-
r_xy = sqrt(max(0._wp, 1._wp - kz*kz))
124-
synthetic_k_x(m_global) = k_mag*r_xy*cos(rn1)
125-
synthetic_k_y(m_global) = k_mag*r_xy*sin(rn1)
126-
synthetic_k_z(m_global) = k_mag*kz
127-
! Solenoidal direction: random unit vector perpendicular to k.
128-
! Build an orthonormal frame: k_hat x (0,0,1) if not degenerate.
129-
if (abs(kz) < 0.9_wp) then
130-
! Cross k_hat with z-hat -> (-k_y, k_x, 0) / r_xy
131-
synthetic_ex(m_global) = -sin(rn1)
132-
synthetic_ey(m_global) = cos(rn1)
133-
synthetic_ez(m_global) = 0._wp
134-
else
135-
! k nearly parallel to z; cross with x-hat instead k_hat x x_hat = (0, -kz, k_y)/norm (normalised)
136-
synthetic_ex(m_global) = 0._wp
137-
synthetic_ey(m_global) = -kz/max(sqrt(kz*kz + (k_mag*sin(rn1))**2), 1e-10_wp)
138-
synthetic_ez(m_global) = (k_mag*sin(rn1))/max(sqrt(kz*kz + (k_mag*sin(rn1))**2), 1e-10_wp)
139-
end if
114+
! Random unit wavevector uniform on the sphere
115+
call s_prng(rn1, seed)
116+
call s_prng(rn2, seed)
117+
khat = f_unit_vector(rn1, rn2)
118+
! Random reference vector projected perpendicular to k by a double
119+
! cross product: sig = khat x (xi x khat) is a unit vector with k.sig = 0
120+
call s_prng(rn1, seed)
121+
call s_prng(rn2, seed)
122+
xi = f_unit_vector(rn1, rn2)
123+
sig_tmp = f_cross(xi, khat)
124+
sig_tmp = sig_tmp/max(sqrt(sum(sig_tmp**2._wp)), 1.e-10_wp)
125+
sig = f_cross(khat, sig_tmp)
126+
synthetic_k_x(m_global) = k_mag*khat(1)
127+
synthetic_k_y(m_global) = k_mag*khat(2)
128+
synthetic_k_z(m_global) = k_mag*khat(3)
129+
synthetic_ex(m_global) = sig(1)
130+
synthetic_ey(m_global) = sig(2)
131+
synthetic_ez(m_global) = sig(3)
140132
end if
141133

142-
call random_number(rn1)
134+
call s_prng(rn1, seed)
143135
synthetic_phase(m_global) = rn1*2._wp*pi
144136

145137
synthetic_amp(m_global) = synth_amp_shell(s)
@@ -293,7 +285,7 @@ contains
293285
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
294286
type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
295287
type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf
296-
integer :: i, j, k, l, n_mode, turb_idx
288+
integer :: i, j, k, l, n_mode, turb_idx, iq
297289
real(wp) :: pos_x, pos_y, pos_z
298290
real(wp) :: gauss_env, G_norm, f_scale
299291
real(wp) :: a_m, phase_arg
@@ -323,7 +315,7 @@ contains
323315
adv_offset = synth_U_inf*mytime
324316

325317
do turb_idx = 1, num_turbulent_sources
326-
$:GPU_PARALLEL_LOOP(private='[j, k, l, n_mode, pos_x, pos_y, pos_z, gauss_env, G_norm, f_scale, a_m, phase_arg, &
318+
$:GPU_PARALLEL_LOOP(private='[j, k, l, n_mode, iq, pos_x, pos_y, pos_z, gauss_env, G_norm, f_scale, a_m, phase_arg, &
327319
& force_x, force_y, force_z, rho_local, u_local_x, u_local_y, u_local_z, in_box]', collapse=3)
328320
do l = 0, p
329321
do k = 0, n
@@ -364,7 +356,11 @@ contains
364356
end do
365357

366358
! Scale: rho * (U_inf/L_x) * G_norm
367-
rho_local = q_prim_vf(eqn_idx%cont%beg)%sf(j, k, l)
359+
! Mixture density = sum of all partial densities (continuity components)
360+
rho_local = 0._wp
361+
do iq = eqn_idx%cont%beg, eqn_idx%cont%end
362+
rho_local = rho_local + q_prim_vf(iq)%sf(j, k, l)
363+
end do
368364
f_scale = rho_local*(synth_U_inf/synth_L(turb_idx, 1))*G_norm
369365

370366
force_x = force_x*f_scale

0 commit comments

Comments
 (0)