Skip to content

bug: 1D flow perturbation applies random factor twice due to index aliasing #1369

@sbryngelson

Description

@sbryngelson

Summary

In src/pre_process/m_perturbation.fpp, s_perturb_surrounding_flow applies a random perturbation to momentum:

q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
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)

In 1D (n == 0, p == 0), num_vels == 1, so eqn_idx%mom%end == eqn_idx%mom%beg. Line 1 writes rand * v0 into the only velocity slot, then line 2 reads that same (now-modified) slot and writes (1 + rand) * rand * v0. The perturbation is applied twice and compounded, giving a result of rand*(1+rand)*v0 instead of the intended (1+rand)*v0.

Fix

Swap the two lines so mom%beg is scaled first, then mom%end is set from the original value:

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)
if (num_vels > 1) q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)

Discovered

Found during review of PR #1365. Pre-existing in master.

Metadata

Metadata

Assignees

Labels

bugSomething isn't working or doesn't seem right

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions