Skip to content

Commit 381ef1e

Browse files
committed
Fix 1D flow perturbation applying random factor twice due to index aliasing
In 1D, eqn_idx%mom%end == eqn_idx%mom%beg (same slot), so the old code wrote rand*v0 into the slot then read that modified value to compute (1+rand)*rand*v0, doubling the perturbation. Fix: save the original mom%beg value before any writes, scale mom%beg first, and only assign mom%end when num_vels > 1 (i.e., multi-D). Fixes #1369
1 parent d513442 commit 381ef1e

1 file changed

Lines changed: 4 additions & 3 deletions

File tree

src/pre_process/m_perturbation.fpp

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ contains
6363

6464
type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
6565
integer :: i, j, k
66-
real(wp) :: rand_real
66+
real(wp) :: rand_real, v_beg
6767

6868
call random_seed()
6969

@@ -72,8 +72,9 @@ contains
7272
do i = 0, m
7373
call random_number(rand_real)
7474
rand_real = rand_real*perturb_flow_mag
75-
q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
76-
q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
75+
v_beg = q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
76+
q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k) = (1._wp + rand_real)*v_beg
77+
if (num_vels > 1) q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*v_beg
7778
if (bubbles_euler) then
7879
q_prim_vf(eqn_idx%alf)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(eqn_idx%alf)%sf(i, j, k)
7980
end if

0 commit comments

Comments
 (0)