Fix cylindrical angular flux scaling and geometric sources - #1824
aliencaocao wants to merge 4 commits into
Conversation
5f8763e to
26000e3
Compare
There was a problem hiding this comment.
🟡 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/rmetric scaling to cylindrical z/θ advection-source terms inm_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.
| #: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 |
| #: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 |
26000e3 to
654b5cb
Compare
|
i don't understand the fix and it's out of date with master |
|
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. Second, the corresponding nonconservative volume-fraction terms used 1/dtheta, while the conservative angular flux used The LF/HLL additions are separate consistency changes. I based them on the existing HLLC NORM_DIR == 3 handling: MFC/src/simulation/m_riemann_solver_hllc.fpp Lines 1529 to 1543 in 22c2aae 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. |
|
I have no idea if this is correct (and it sounds like you do not either). Maybe someone else does. |
There was a problem hiding this comment.
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.
| $:GPU_LOOP(parallelism='[seq]') | ||
| do i = 1, sys_size | ||
| flux_gsrc_rsx_vf(${SF('')}$, i) = 0._wp | ||
| end do |
| #: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 |
| 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) |
Conflict in m_rhs.fpp was the local declarations: master added pres_K, alpha_K, alpha_rho_K and blkmod_K for per-phase EOS evaluation, this branch dropped velocity_val. Kept master's additions, which the merged body uses, and the removal of velocity_val, which it no longer references.
Lines of Code
|
Fixes #1817.
In MFC’s cylindrical coordinates, the third sweep represents the angular direction. Its numerical flux difference must be divided by r times the angular cell width. The current RHS additionally multiplies this difference by angular velocity, changing transport rates and suppressing angular pressure forces when that velocity is zero.
The volume fraction source terms also need the same radial metric. The reproducer uses model_eqns=2 with alt_soundspeed=F and no phase relaxation, so its initially uniform fractions must remain constant. Applying the same metric to their flux and source terms restores the required cancellation.
This PR removes the extra velocity factor and adds the missing radial scaling. It also supplies the angular geometric momentum sources in HLL and LF: the centrifugal contribution to radial momentum and the coupling between radial and angular momentum. Removing the matching pressure contribution from the numerical momentum flux gives the required geometric term.
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:
If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.
Acknowledgement
PR template credit: junegunn