Skip to content

Fix cylindrical angular flux scaling and geometric sources - #1824

Draft
aliencaocao wants to merge 2 commits into
MFlowCode:masterfrom
aliencaocao:fix/cylindrical-angular-flux
Draft

Fix cylindrical angular flux scaling and geometric sources#1824
aliencaocao wants to merge 2 commits into
MFlowCode:masterfrom
aliencaocao:fix/cylindrical-angular-flux

Conversation

@aliencaocao

@aliencaocao aliencaocao commented Sep 5, 2026

Copy link
Copy Markdown
Contributor

Closes #1817. For the context and fix description, please refer to the issue post.

Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

All contributions are expected to demonstrate:

  • A clear understanding of the codebase
  • Alignment with product direction
  • Thoughtful reasoning behind changes
  • Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.


Acknowledgement

  • I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

PR template credit: junegunn

Copilot AI lite review requested due to automatic review settings September 5, 2026 14:13
@aliencaocao
aliencaocao marked this pull request as draft September 5, 2026 14:13
@aliencaocao
aliencaocao force-pushed the fix/cylindrical-angular-flux branch from 5f8763e to 26000e3 Compare September 5, 2026 14:19

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟡 Changes recommended

The new cylindrical z/θ geometric-source flux logic in LF/HLL does not correctly cancel the bubbles_euler pressure term (pres - ptilde), which can produce incorrect geometric sources when bubbles_euler is enabled.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Pull request overview

This PR fixes a cylindrical-coordinate correctness bug in the z/azimuthal sweep so that angular transport and advection-source terms use the proper metric scaling and do not double-count velocity factors, addressing the issue where uniform mixture fractions drift after a step in cylindrical setups (Issue #1817).

Changes:

  • Remove the extra velocity factor from the cylindrical angular (z/θ) flux divergence in m_rhs.fpp.
  • Apply the missing 1/r metric scaling to cylindrical z/θ advection-source terms in m_rhs.fpp.
  • Add cylindrical z/θ geometric-source flux handling in the LF and HLL Riemann solvers.
File summaries
File Description
src/simulation/m_riemann_solver_lf.fpp Adds cylindrical z/θ geometric-source flux population for LF solver.
src/simulation/m_riemann_solver_hll.fpp Adds cylindrical z/θ geometric-source flux population for HLL solver.
src/simulation/m_rhs.fpp Corrects cylindrical z/θ flux divergence scaling and adds missing 1/r scaling to advection-source terms.
Review details
  • Files reviewed: 3/3 changed files
  • Comments generated: 2
  • Review effort level: Lite

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment on lines +671 to +681
#:if (NORM_DIR == 3)
if (grid_geometry == 3) then
$:GPU_LOOP(parallelism='[seq]')
do i = 1, sys_size
flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp
end do
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1) = &
& -flux_rsx_vf(${SF('')}$, eqn_idx%mom%end) + (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%end) = flux_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1)
end if
#:endif
Comment on lines +352 to +362
#:if (NORM_DIR == 3)
if (grid_geometry == 3) then
$:GPU_LOOP(parallelism='[seq]')
do i = 1, sys_size
flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp
end do
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1) = &
& -flux_rsx_vf(${SF('')}$, eqn_idx%mom%end) + (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%end) = flux_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1)
end if
#:endif
@aliencaocao
aliencaocao force-pushed the fix/cylindrical-angular-flux branch from 26000e3 to 654b5cb Compare September 5, 2026 14:23
@sbryngelson

Copy link
Copy Markdown
Member

i don't understand the fix and it's out of date with master

@aliencaocao

aliencaocao commented Sep 5, 2026

Copy link
Copy Markdown
Contributor Author

My understanding is based on code-level rather than from a scientific angle, so please correct me if I have misunderstood the intended discretization.

The reproducer exposed two issues in m_rhs.fpp.
First, flux_n(3) is already the numerical flux returned by the Riemann solver. The cylindrical z sweep was multiplying its face difference by the cell-centered angular velocity again: u_theta * (F_left - F_right) / (r * dtheta)
I changed this to (F_left - F_right) / (r * dtheta)

Second, the corresponding nonconservative volume-fraction terms used 1/dtheta, while the conservative angular flux used 1/(r*dtheta). This prevented the two terms from cancelling for a uniform volume fraction. The reproducer starts with alpha values of 0.4 and 0.6 everywhere, but they drift after one step. Applying the same 1/(r*dtheta) metric to both terms removes that drift.

The LF/HLL additions are separate consistency changes. I based them on the existing HLLC NORM_DIR == 3 handling:

#:if (NORM_DIR == 3)
if (grid_geometry == 3) then
$:GPU_LOOP(parallelism='[seq]')
do i = 1, sys_size
flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp
end do
flux_gsrc_rsx_vf(${SF('')}$, &
& eqn_idx%mom%beg + 1) = -f_compute_hllc_star_momentum_flux(rho_L, &
& rho_R, vel_L(dir_idx(1)), vel_R(dir_idx(1)), s_M, s_P, s_S, xi_L, &
& xi_R, xi_M, xi_P, dir_flg(dir_idx(1)))
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%end) = flux_rsx_vf(${SF('')}$, &
& eqn_idx%mom%beg + 1)
end if
#:endif

HLLC explicitly fills the angular geometric-source entries for radial and angular momentum. LF and HLL did not. Since LF/HLL include pressure in their momentum flux directly, the new code removes the matching face-pressure contribution to obtain the pressure-free geometric term.

If you prefer, I can separate the LF/HLL consistency change from the two m_rhs.fpp changes that directly address the reproducer.

@sbryngelson

Copy link
Copy Markdown
Member

I have no idea if this is correct (and it sounds like you do not either). Maybe someone else does.

@github-actions

github-actions Bot commented Sep 5, 2026

Copy link
Copy Markdown

Lines of Code

File Lines Diff
src/simulation/m_riemann_solver_hll.fpp 612 +17
src/simulation/m_riemann_solver_lf.fpp 534 +17
src/simulation/m_rhs.fpp 1950 -1
Directory Lines Diff
simulation 27701 +33
total 45471 +33

@sbryngelson
sbryngelson requested a balanced review from Copilot September 5, 2026 20:58

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Warning

Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.

Pull request overview

Copilot reviewed 3 out of 3 changed files in this pull request and generated 3 comments.

Comment on lines +686 to +689
$:GPU_LOOP(parallelism='[seq]')
do i = 1, sys_size
flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp
end do
Comment on lines +328 to +344
#:if (NORM_DIR == 3)
if (grid_geometry == 3) then
$:GPU_LOOP(parallelism='[seq]')
do i = 1, sys_size
flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp
end do
if (bubbles_euler) then
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1) = -flux_rsx_vf(${SF('')}$, &
& eqn_idx%mom%end) + (s_M*(pres_R - ptilde_R) - s_P*(pres_L - ptilde_L)) &
& /(s_M - s_P)
else
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1) = -flux_rsx_vf(${SF('')}$, &
& eqn_idx%mom%end) + (s_M*pres_R - s_P*pres_L)/(s_M - s_P)
end if
flux_gsrc_rsx_vf(${SF('')}$, eqn_idx%mom%end) = flux_rsx_vf(${SF('')}$, eqn_idx%mom%beg + 1)
end if
#:endif
Comment thread src/simulation/m_rhs.fpp
Comment on lines 1318 to +1321
inv_ds = 1._wp/(dz(k)*y_cc(q))
velocity_val = q_prim_vf%vf(eqn_idx%cont%end + idir)%sf(l, q, k)
flux_face1 = flux_n(3)%vf(j)%sf(l, q, k - 1)
flux_face2 = flux_n(3)%vf(j)%sf(l, q, k)
rhs_vf(j)%sf(l, q, k) = rhs_vf(j)%sf(l, q, k) + inv_ds*velocity_val*(flux_face1 - flux_face2)
rhs_vf(j)%sf(l, q, k) = rhs_vf(j)%sf(l, q, k) + inv_ds*(flux_face1 - flux_face2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

Bug: cylindrical angular flux has an extra velocity factor and missing metric scaling

3 participants