Revert HLL/LF adv gsrc change; update comment to explain PR #800 intent

The previous commit incorrectly set flux_gsrc(adv) = 0 for HLL and LF,
which would have reverted PR #800. Per the analysis in issue #1402 comment
by ChrisZYJ: flux_rs(adv) in HLL/LF is the numerical diffusion component
(not transport), and any face-flux in cylindrical coordinates requires the
full geometric divergence (1/r)*d(rF)/dr = dF/dr + F/r. The flux_gsrc term
provides the F/r correction for the diffusion. The transport half is handled
separately via flux_src_rs (non-conservative source). PR #800 was correct.

The comment is updated to document this so future readers do not repeat the
mistake. The convergence test is unaffected: for num_fluids=1, alpha=1
everywhere so flux_rs(adv) = 0 at all faces regardless.

Author: sbryngelson

src/simulation/m_riemann_solvers.fpp

@@ -772,10 +772,10 @@ contains
                                     ! Recalculating the radial momentum geometric source flux
                                     flux_gsrc_rs${XYZ}$_vf(j, k, l, eqn_idx%cont%end + 2) = flux_rs${XYZ}$_vf(j, k, l, &
                                                            & eqn_idx%cont%end + 2) - (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
-                                    ! Geometrical source of the void fraction(s) is zero
+                                    ! Cylindrical correction for HLL/LF numerical diffusion of void fraction(s)
                                     $:GPU_LOOP(parallelism='[seq]')
                                     do i = eqn_idx%adv%beg, eqn_idx%adv%end
-                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
+                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
                                     end do
                                 end if
 
@@ -1382,10 +1382,10 @@ contains
                                     ! Recalculating the radial momentum geometric source flux
                                     flux_gsrc_rs${XYZ}$_vf(j, k, l, eqn_idx%cont%end + 2) = flux_rs${XYZ}$_vf(j, k, l, &
                                                            & eqn_idx%cont%end + 2) - (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
-                                    ! Geometrical source of the void fraction(s) is zero
+                                    ! Cylindrical correction for HLL/LF numerical diffusion of void fraction(s)
                                     $:GPU_LOOP(parallelism='[seq]')
                                     do i = eqn_idx%adv%beg, eqn_idx%adv%end
-                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
+                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
                                     end do
                                 end if