refactor: unify equation/range indices into sys_idx; clean BC globals#1365
Open
sbryngelson wants to merge 12 commits intoMFlowCode:masterfrom
Open
refactor: unify equation/range indices into sys_idx; clean BC globals#1365sbryngelson wants to merge 12 commits intoMFlowCode:masterfrom
sbryngelson wants to merge 12 commits intoMFlowCode:masterfrom
Conversation
…remove dead bcxb/bcxe module vars - Add bc_xyz_info type to m_derived_types.fpp grouping bc x/y/z info - s_populate_capillary_buffers now takes explicit bc argument instead of implicitly depending on bc_x/bc_y/bc_z module globals - Access pattern is now bc%x%beg, bc%x%end etc. (was implicit global) - Remove unused bcxb/bcxe/bcyb/bcye/bczb/bcze scalars from m_cbc.fpp (declared and GPU-declared but never read)
- Add eqn_idx_info type to m_derived_types.fpp with members E, n, alf, gamma, pi_inf, c - Replace 6 standalone integer globals (E_idx, alf_idx, c_idx, n_idx, gamma_idx, pi_inf_idx) with type(eqn_idx_info) :: eqn_idx in all three targets' m_global_parameters.fpp - Update GPU_DECLARE and GPU_UPDATE in simulation to use eqn_idx struct - ~672 use-sites updated to eqn_idx%E, eqn_idx%alf, etc.
… e_idx typo in m_rhs
Claude Code ReviewHead SHA: 28b88d0 Files changed:
FindingsGPU correctness:
|
sbryngelson
changed the title
…ences In m_rhs.fpp, replace 'sys_size - 1' with 'eqn_idx%c - 1' when surface_tension is active and remove the comment noting the fragile assumption (eqn_idx%c makes the intent explicit). In m_start_up.fpp, remove misleading '! eqn_idx%adv%end' comments on restart loops that correctly iterate over all sys_size variables.
Contributor
Code Review by Qodo
|
![qodo-code-review[bot]](https://avatars.githubusercontent.com/in/484649?s=60&v=4){kind=small}
Comment on lines
+922
to
929
| @:ALLOCATE(qbmm_idx%rs(nb), qbmm_idx%vs(nb)) | ||
| @:ALLOCATE(qbmm_idx%ps(nb), qbmm_idx%ms(nb)) | ||
|
|
||
| gam = bub_pp%gam_g | ||
|
|
||
| if (qbmm) then | ||
| @:ALLOCATE(bub_idx%moms(nb, nmom)) | ||
| @:ALLOCATE(qbmm_idx%moms(nb, nmom)) | ||
| do i = 1, nb |
Contributor
There was a problem hiding this comment.
1. qbmm_idx allocations not freed 📘 Rule violation ☼ Reliability
qbmm_idx%rs/vs/ps/ms (and sometimes qbmm_idx%moms) are allocated via @:ALLOCATE but are never deallocated in s_finalize_global_parameters_module, risking memory leaks and non-deterministic resource cleanup. This violates the required @:ALLOCATE/@:DEALLOCATE pairing policy.
Agent Prompt
## Issue description
`qbmm_idx` members are allocated with `@:ALLOCATE` in `s_initialize_global_parameters_module` but are never freed in `s_finalize_global_parameters_module`.
## Issue Context
This PR replaced `bub_idx` allocatables with `qbmm_idx` allocatables; the finalization path must be updated to deallocate the new allocations (including `qbmm_idx%moms` when `qbmm` is enabled).
## Fix Focus Areas
- src/simulation/m_global_parameters.fpp[922-929]
- src/simulation/m_global_parameters.fpp[995-996]
- src/simulation/m_global_parameters.fpp[1299-1340]
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
Comment on lines
+727
to
736
| if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_SLIP_WALL) & | ||
| & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_SLIP_WALL)) then | ||
| call s_compute_slip_wall_L(lambda, L, rho, c, dpres_ds, dvel_ds) | ||
| else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_BUFFER) & | ||
| & .or. (cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_BUFFER)) then | ||
| else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_BUFFER) & | ||
| & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_BUFFER)) then | ||
| call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, & | ||
| & dvel_ds, dadv_ds, dYs_ds) | ||
| else if ((cbc_loc == -1 .and. bc${XYZ}$b == BC_CHAR_NR_SUB_INFLOW) & | ||
| & .or. (cbc_loc == 1 .and. bc${XYZ}$e == BC_CHAR_NR_SUB_INFLOW)) then | ||
| else if ((cbc_loc == -1 .and. bc_${XYZ}$%beg == BC_CHAR_NR_SUB_INFLOW) & | ||
| & .or. (cbc_loc == 1 .and. bc_${XYZ}$%end == BC_CHAR_NR_SUB_INFLOW)) then | ||
| call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, dpres_ds, dvel_ds) |
Contributor
There was a problem hiding this comment.
2. Bc beg/end missing on gpu 🐞 Bug ≡ Correctness
OpenACC device code in m_cbc now branches on bc_x/bc_y/bc_z%beg/%end, but those members are not included in the OpenACC GPU_DECLARE list and are never GPU_UPDATE’d, so kernels can read undefined values and apply incorrect characteristic boundary conditions.
Agent Prompt
### Issue description
CBC OpenACC kernels now read `bc_x%beg/%end` (and y/z) on-device, but these fields are not part of the OpenACC `GPU_DECLARE` set and are not copied to the device during GPU initialization. This can make boundary-condition selection nondeterministic/wrong in GPU runs.
### Issue Context
- `m_cbc` uses `bc_${XYZ}$%beg/%end` inside device loops to choose the CBC formulation.
- Under OpenACC, `m_global_parameters` only declares `bc_*%vb*`/`%ve*` on the device.
- `s_initialize_gpu_vars` updates `bc_*%vb*`/`%ve*` and `bc_*%grcbc_*` but not `bc_*%beg/%end`.
### Fix Focus Areas
- src/simulation/m_global_parameters.fpp[211-223]
- src/simulation/m_start_up.fpp[1035-1061]
- src/simulation/m_cbc.fpp[725-786]
### Expected fix
- Add `bc_x%beg, bc_x%end, bc_y%beg, bc_y%end, bc_z%beg, bc_z%end` to the OpenACC `GPU_DECLARE(create=...)` list (or declare the full `bc_x/bc_y/bc_z` structs for OpenACC).
- Add a matching `GPU_UPDATE(device='[...]')` for those beg/end members in `s_initialize_gpu_vars` after they are finalized on the host.
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
Contributor
Code Review by Qodo
|
Contributor
📝 WalkthroughWalkthroughThis pull request performs a large-scale refactoring of equation index management across the codebase. The changes introduce two new derived types— 🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
Contributor
There was a problem hiding this comment.
Actionable comments posted: 11
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (4)
src/simulation/m_time_steppers.fpp (1)
916-959:⚠️ Potential issue | 🟡 MinorMissing deallocations for
eqn_idx%nandeqn_idx%callocations.The allocation block allocates
q_prim_vf(eqn_idx%n)%sfwhenadv_nis true (Line 239) andq_prim_vf(eqn_idx%c)%sfwhensurface_tensionis true (Line 290), but the finalization subroutine does not deallocate these fields.🔧 Proposed fix to add missing deallocations
if (bubbles_euler) then do i = eqn_idx%bub%beg, eqn_idx%bub%end @:DEALLOCATE(q_prim_vf(i)%sf) end do + if (adv_n) then + @:DEALLOCATE(q_prim_vf(eqn_idx%n)%sf) + end if end if if (model_eqns == 3) then do i = eqn_idx%int_en%beg, eqn_idx%int_en%end @:DEALLOCATE(q_prim_vf(i)%sf) end do end if + + if (surface_tension) then + @:DEALLOCATE(q_prim_vf(eqn_idx%c)%sf) + end if end ifAs per coding guidelines: "Every
@:ALLOCATE(...)Fypp macro MUST have a matching@:DEALLOCATE(...)."src/simulation/m_rhs.fpp (1)
1727-1747:⚠️ Potential issue | 🟠 Major
s_finalize_rhs_moduleno longer mirrors theq_cons_qp/q_prim_qpallocation paths.
q_cons_qp%vf(l)%sfis allocated forl = 1:sys_size, andq_prim_qpis also allocated pasteqn_idx%Eins_initialize_rhs_module, but this finalizer only freesmom:Eand does alias cleanup forcont/adv. That leaves the conservative continuity/advection targets and every primitive field aboveeqn_idx%Elive at shutdown.Proposed direction
- if (.not. igr) then - do j = eqn_idx%cont%beg, eqn_idx%cont%end - if (relativity) then - @:DEALLOCATE(q_cons_qp%vf(j)%sf) - @:DEALLOCATE(q_prim_qp%vf(j)%sf) - else - nullify (q_prim_qp%vf(j)%sf) - end if - end do - - do j = eqn_idx%adv%beg, eqn_idx%adv%end - nullify (q_prim_qp%vf(j)%sf) - end do - - do j = eqn_idx%mom%beg, eqn_idx%E - @:DEALLOCATE(q_cons_qp%vf(j)%sf) - @:DEALLOCATE(q_prim_qp%vf(j)%sf) - end do - end if + if (.not. igr) then + do j = 1, sys_size + if (associated(q_prim_qp%vf(j)%sf, q_cons_qp%vf(j)%sf)) then + nullify (q_prim_qp%vf(j)%sf) + end if + end do + + do j = 1, sys_size + if (associated(q_cons_qp%vf(j)%sf)) then + @:DEALLOCATE(q_cons_qp%vf(j)%sf) + end if + if (associated(q_prim_qp%vf(j)%sf)) then + @:DEALLOCATE(q_prim_qp%vf(j)%sf) + end if + end do + end ifAs per coding guidelines "Every
@:ALLOCATE(...)Fypp macro MUST have a matching@:DEALLOCATE(...)."src/pre_process/m_data_output.fpp (1)
619-675:⚠️ Potential issue | 🟡 MinorUse unlimited-width integer formatting in
indices.dat.The new file still formats indices with
I3/I2. Once a case pushes any equation index past those widths, Fortran will print***/**, which corrupts the mapping you just made easier to consume.Suggested fix
- write (iu, '(I3,A20,A20)') i, "\alpha_{" // trim(temp) // "} \rho_{" // trim(temp) // "}", & + write (iu, '(I0,1X,A20,1X,A20)') i, "\alpha_{" // trim(temp) // "} \rho_{" // trim(temp) // "}", & & "\alpha_{" // trim(temp) // "} \rho" ... - write (iu, '(I3,A20,A20)') i, "\rho u_" // coord(i - momxb + 1), "u_" // coord(i - momxb + 1) + write (iu, '(I0,1X,A20,1X,A20)') i, "\rho u_" // coord(i - momxb + 1), "u_" // coord(i - momxb + 1) ... - if (eqn_idx%E /= 0) write (iu, '(I3,A20,A20)') eqn_idx%E, "\rho U", "p" + if (eqn_idx%E /= 0) write (iu, '(I0,1X,A20,1X,A20)') eqn_idx%E, "\rho U", "p" ... - write (iu, '(I3,A20,A20)') i, "\alpha_{" // trim(temp) // "}", "\alpha_{" // trim(temp) // "}" + write (iu, '(I0,1X,A20,1X,A20)') i, "\alpha_{" // trim(temp) // "}", "\alpha_{" // trim(temp) // "}" ... - write (iu, '(I3,A20,A20)') chemxb + i - 1, "Y_{" // trim(species_names(i)) // "} \rho", & + write (iu, '(I0,1X,A20,1X,A20)') chemxb + i - 1, "Y_{" // trim(species_names(i)) // "} \rho", & & "Y_{" // trim(species_names(i)) // "}" ... - if (beg /= 0) write (iu, '("[",I2,",",I2,"]",A)') beg, end, label + if (beg /= 0) write (iu, '("[",I0,",",I0,"]",A)') beg, end, labelBased on learnings "Flag modifications to public subroutine signatures, parameter defaults, or output formats as they affect downstream code and user workflows."
src/simulation/m_cbc.fpp (1)
727-785:⚠️ Potential issue | 🔴 CriticalThe BC struct members
begandendare not GPU-declared for OpenACC builds.The GPU branches in lines 727–785 read
bc_${XYZ}$%begandbc_${XYZ}$%enddirectly within GPU kernels. However,src/simulation/m_global_parameters.fpplines 218–222 show that for OpenACC, only the vector fields (bc_x%vb1,bc_x%vb2,bc_x%vb3,bc_x%ve1,bc_x%ve2,bc_x%ve3) are GPU_DECLARE'd; the scalarbegandendmembers are left host-only. OpenMP correctly declares the full struct. This mismatch means OpenACC GPU kernels may silently read stale or uninitialized host memory when branching on BC metadata, causing divergence from CPU behavior.Add GPU declarations for
bc_x%beg,bc_x%end,bc_y%beg,bc_y%end,bc_z%beg,bc_z%endinsrc/simulation/m_global_parameters.fppfor the OpenACC path, or refactor to use the already-declared vector fields.
🧹 Nitpick comments (6)
claude-code-review.yml (1)
145-164: Context fetching handles errors silently.The
2>/dev/null || trueat line 151 and similar patterns silently swallow API errors. While acceptable for "optional" context, consider adding minimal logging (e.g.,echo "Skipped: $path" >&2) when files fail to fetch for easier debugging during workflow development.src/pre_process/m_perturbation.fpp (1)
46-52: Useeqn_idx%advhere instead ofeqn_idx%E + offset.These lines work only because the advection block currently sits immediately after energy. That reintroduces the same implicit layout coupling this refactor is trying to remove.
♻️ Proposed refactor
- perturb_alpha = q_prim_vf(eqn_idx%E + perturb_sph_fluid)%sf(i, j, k) + perturb_alpha = q_prim_vf(eqn_idx%adv%beg + perturb_sph_fluid - 1)%sf(i, j, k) ... - q_prim_vf(l)%sf(i, j, k) = q_prim_vf(eqn_idx%E + l)%sf(i, j, k)*fluid_rho(l) + q_prim_vf(eqn_idx%cont%beg + l - 1)%sf(i, j, k) = & + q_prim_vf(eqn_idx%adv%beg + l - 1)%sf(i, j, k)*fluid_rho(l)src/common/m_boundary_common.fpp (1)
1054-1058: Run the full shared-target test matrix for this API change.Line 1054 changes a
src/commonboundary routine signature on the capillary path. Build coverage is good, but I’d still want./mfc.sh test -j 8across pre_process, simulation, and post_process before merging.As per coding guidelines, "
src/common/**/*.{fpp,f90}: Test changes tosrc/common/comprehensively by running ALL three targets (./mfc.sh test -j 8) since this shared code affects all three executables".src/simulation/m_start_up.fpp (1)
1035-1049: Duplicate GPU_UPDATE with overlapping variables.Lines 1035-1038 and line 1049 both update
adv_n,adap_dt,adap_dt_tol,adap_dt_max_iters,eqn_idx%n,pi_fac, andlow_Machto the device. This appears to be unintentional duplication that could be removed for clarity and minor performance improvement.♻️ Suggested fix to remove duplicate GPU_UPDATE
$:GPU_UPDATE(device='[R0ref, p0ref, rho0ref, ss, pv, vd, mu_l, mu_v, mu_g, gam_v, gam_g, M_v, M_g, R_v, R_g, Tw, cp_v, & & cp_g, k_vl, k_gl, gam, gam_m, Eu, Ca, Web, Re_inv, Pe_c, phi_vg, phi_gv, omegaN, bubbles_euler, & & polytropic, polydisperse, qbmm, ptil, bubble_model, thermal, poly_sigma, adv_n, adap_dt, adap_dt_tol, & & adap_dt_max_iters, eqn_idx%n, pi_fac, low_Mach]') if (bubbles_euler) then $:GPU_UPDATE(device='[weight, R0]') if (.not. polytropic) then $:GPU_UPDATE(device='[pb0, Pe_T, k_g, k_v, mass_g0, mass_v0, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c]') else if (qbmm) then $:GPU_UPDATE(device='[pb0]') end if end if - $:GPU_UPDATE(device='[adv_n, adap_dt, adap_dt_tol, adap_dt_max_iters, eqn_idx%n, pi_fac, low_Mach]') - $:GPU_UPDATE(device='[acoustic_source, num_source]')src/simulation/m_riemann_solvers.fpp (2)
253-258: Prefereqn_idx%alfovereqn_idx%E + ifor alpha slots.These accesses still assume the alpha block starts immediately after pressure. Since the new index holder already exposes
eqn_idx%alf, using that explicitly would finish the decoupling this refactor is aiming for and avoid baking the old layout back into the solver.Representative change
- alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, eqn_idx%E + i) - alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, eqn_idx%E + i) + alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, eqn_idx%alf - 1 + i) + alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, eqn_idx%alf - 1 + i)Also applies to: 935-940, 1821-1859, 2238-2261, 2443-2467, 2840-2887, 3383-3388
2236-2245: Drop the duplicated alpha extraction loop.The second loop reloads the same
alpha_L/alpha_Rvalues and does no extra work. In this GPU hot path it only adds instruction count and register pressure.Proposed fix
$:GPU_LOOP(parallelism='[seq]') do i = 1, num_fluids alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, eqn_idx%E + i) alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, eqn_idx%E + i) end do - $:GPU_LOOP(parallelism='[seq]') - do i = 1, num_fluids - alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, eqn_idx%E + i) - alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, eqn_idx%E + i) - end doBased on learnings: Avoid unnecessary host/device transfers in hot loops, redundant allocations, and algorithmic inefficiency.
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: CHILL
Plan: Pro
Run ID: ac3fd33f-76bc-45a2-90ba-13bfdcedc3ed
📒 Files selected for processing (41)
claude-code-review.ymlsrc/common/include/1dHardcodedIC.fppsrc/common/include/2dHardcodedIC.fppsrc/common/include/3dHardcodedIC.fppsrc/common/m_boundary_common.fppsrc/common/m_chemistry.fppsrc/common/m_derived_types.fppsrc/common/m_phase_change.fppsrc/common/m_variables_conversion.fppsrc/post_process/m_data_output.fppsrc/post_process/m_derived_variables.fppsrc/post_process/m_global_parameters.fppsrc/post_process/m_start_up.fppsrc/pre_process/m_assign_variables.fppsrc/pre_process/m_data_output.fppsrc/pre_process/m_global_parameters.fppsrc/pre_process/m_icpp_patches.fppsrc/pre_process/m_initial_condition.fppsrc/pre_process/m_perturbation.fppsrc/pre_process/m_start_up.fppsrc/simulation/m_acoustic_src.fppsrc/simulation/m_body_forces.fppsrc/simulation/m_bubbles_EE.fppsrc/simulation/m_bubbles_EL.fppsrc/simulation/m_cbc.fppsrc/simulation/m_compute_cbc.fppsrc/simulation/m_data_output.fppsrc/simulation/m_global_parameters.fppsrc/simulation/m_hyperelastic.fppsrc/simulation/m_hypoelastic.fppsrc/simulation/m_ibm.fppsrc/simulation/m_igr.fppsrc/simulation/m_pressure_relaxation.fppsrc/simulation/m_qbmm.fppsrc/simulation/m_rhs.fppsrc/simulation/m_riemann_solvers.fppsrc/simulation/m_sim_helpers.fppsrc/simulation/m_start_up.fppsrc/simulation/m_surface_tension.fppsrc/simulation/m_time_steppers.fppsrc/simulation/m_viscous.fpp
Comment on lines
+672
to
+673
| call s_compute_pressure(qK_cons_vf(eqn_idx%E)%sf(j, k, l), qK_cons_vf(eqn_idx%alf)%sf(j, k, l), dyn_pres_K, & | ||
| & pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres, T, pres_mag=pres_mag) |
Contributor
There was a problem hiding this comment.
Don't index qK_cons_vf with eqn_idx%alf unconditionally.
alf is unused for most EOS branches, but this call still evaluates qK_cons_vf(eqn_idx%alf) first. eqn_idx%alf is not initialized for every layout, so model-1/pre-process paths can read an undefined slot before s_compute_pressure even runs.
Suggested fix
- call s_compute_pressure(qK_cons_vf(eqn_idx%E)%sf(j, k, l), qK_cons_vf(eqn_idx%alf)%sf(j, k, l), dyn_pres_K, &
- & pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres, T, pres_mag=pres_mag)
+ if (model_eqns == 4 .or. bubbles_euler) then
+ vftmp = qK_cons_vf(eqn_idx%alf)%sf(j, k, l)
+ else
+ vftmp = 0._wp
+ end if
+
+ call s_compute_pressure(qK_cons_vf(eqn_idx%E)%sf(j, k, l), vftmp, dyn_pres_K, &
+ & pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres, T, pres_mag=pres_mag)
Comment on lines
+563
to
568
| allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb)) | ||
| allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb)) | ||
|
|
||
| if (qbmm) then | ||
| allocate (bub_idx%moms(nb, nmom)) | ||
| allocate (qbmm_idx%moms(nb, nmom)) | ||
| do i = 1, nb |
Contributor
There was a problem hiding this comment.
Finalize the new qbmm_idx buffers as well.
These arrays are now allocated during initialization, but the module finalizer still never deallocates them. That leaves the QBMM index storage leaking across repeated post-process runs in the same process.
Suggested fix
impure subroutine s_finalize_global_parameters_module
@@
+ if (allocated(qbmm_idx%moms)) deallocate(qbmm_idx%moms)
+ if (allocated(qbmm_idx%ms)) deallocate(qbmm_idx%ms)
+ if (allocated(qbmm_idx%ps)) deallocate(qbmm_idx%ps)
+ if (allocated(qbmm_idx%vs)) deallocate(qbmm_idx%vs)
+ if (allocated(qbmm_idx%rs)) deallocate(qbmm_idx%rs)Also applies to: 643-644
Comment on lines
+614
to
+619
| allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb)) | ||
| allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb)) | ||
|
|
||
| if (qbmm) then | ||
| allocate (bub_idx%moms(nb, nmom)) | ||
| allocate (bub_idx%fullmom(nb,0:nmom,0:nmom)) | ||
| allocate (qbmm_idx%moms(nb, nmom)) | ||
| allocate (qbmm_idx%fullmom(nb,0:nmom,0:nmom)) |
Contributor
There was a problem hiding this comment.
Add a matching teardown for the new qbmm_idx allocations.
These new module allocations never get released in s_finalize_global_parameters_module, so repeated initialize/finalize paths will leak the QBMM index tables.
Suggested fix
impure subroutine s_finalize_global_parameters_module
@@
+ if (allocated(qbmm_idx%fullmom)) deallocate(qbmm_idx%fullmom)
+ if (allocated(qbmm_idx%moms)) deallocate(qbmm_idx%moms)
+ if (allocated(qbmm_idx%ms)) deallocate(qbmm_idx%ms)
+ if (allocated(qbmm_idx%ps)) deallocate(qbmm_idx%ps)
+ if (allocated(qbmm_idx%vs)) deallocate(qbmm_idx%vs)
+ if (allocated(qbmm_idx%rs)) deallocate(qbmm_idx%rs)Also applies to: 696-697
Comment on lines
+382
to
383
| q_prim_vf(eqn_idx%alf)%sf(i, j, & | ||
| & 0) = patch_icpp(patch_id)%alpha(1)*exp(-0.5_wp*((myr - radius)**2._wp)/(thickness/3._wp)**2._wp) |
Contributor
There was a problem hiding this comment.
Guard eqn_idx%alf before writing the annulus profile.
Line 382 and Line 444 dereference q_prim_vf(eqn_idx%alf) unconditionally. eqn_idx%alf is model-specific and inactive outside model_eqns == 4, so these paths can fall through to q_prim_vf(0) on other models.
Suggested fix
subroutine s_icpp_varcircle(patch_id, patch_id_fp, q_prim_vf)
+ if (eqn_idx%alf <= 0) then
+ call s_mpi_abort('s_icpp_varcircle requires model_eqns == 4')
+ end if
...
subroutine s_icpp_3dvarcircle(patch_id, patch_id_fp, q_prim_vf)
+ if (eqn_idx%alf <= 0) then
+ call s_mpi_abort('s_icpp_3dvarcircle requires model_eqns == 4')
+ end ifAlso applies to: 444-445
Comment on lines
+75
to
79
| 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) | ||
| if (bubbles_euler) then | ||
| q_prim_vf(alf_idx)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(alf_idx)%sf(i, j, k) | ||
| q_prim_vf(eqn_idx%alf)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(eqn_idx%alf)%sf(i, j, k) | ||
| end if |
Contributor
There was a problem hiding this comment.
Guard this momentum perturbation for 1D.
On Line 75, eqn_idx%mom%end aliases eqn_idx%mom%beg when num_vels == 1. That means the first assignment overwrites the only velocity component before Line 76 scales it, so the 1D path gets a different perturbation than intended.
🐛 Proposed fix
- q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
+ if (num_vels > 1) then
+ q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k)
+ end if
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)📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
Suggested change
| 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) | |
| if (bubbles_euler) then | |
| q_prim_vf(alf_idx)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(alf_idx)%sf(i, j, k) | |
| q_prim_vf(eqn_idx%alf)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(eqn_idx%alf)%sf(i, j, k) | |
| end if | |
| if (num_vels > 1) then | |
| q_prim_vf(eqn_idx%mom%end)%sf(i, j, k) = rand_real*q_prim_vf(eqn_idx%mom%beg)%sf(i, j, k) | |
| end if | |
| 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 (bubbles_euler) then | |
| q_prim_vf(eqn_idx%alf)%sf(i, j, k) = (1._wp + rand_real)*q_prim_vf(eqn_idx%alf)%sf(i, j, k) | |
| end if |
Comment on lines
+922
to
929
| @:ALLOCATE(qbmm_idx%rs(nb), qbmm_idx%vs(nb)) | ||
| @:ALLOCATE(qbmm_idx%ps(nb), qbmm_idx%ms(nb)) | ||
|
|
||
| gam = bub_pp%gam_g | ||
|
|
||
| if (qbmm) then | ||
| @:ALLOCATE(bub_idx%moms(nb, nmom)) | ||
| @:ALLOCATE(qbmm_idx%moms(nb, nmom)) | ||
| do i = 1, nb |
Contributor
There was a problem hiding this comment.
Pair the new qbmm_idx allocations with @:DEALLOCATE in teardown.
qbmm_idx%rs/vs/ps/ms and qbmm_idx%moms are now allocated here, but s_finalize_global_parameters_module never releases them. That leaves the new QBMM mapping storage leaking across reinitialization and violates the module’s existing allocate/finalize contract. As per coding guidelines, every @:ALLOCATE(...) Fypp macro MUST have a matching @:DEALLOCATE(...).
Suggested fix
impure subroutine s_finalize_global_parameters_module
@@
+ if (allocated(qbmm_idx%moms)) @:DEALLOCATE(qbmm_idx%moms)
+ if (allocated(qbmm_idx%ms)) @:DEALLOCATE(qbmm_idx%ms)
+ if (allocated(qbmm_idx%ps)) @:DEALLOCATE(qbmm_idx%ps)
+ if (allocated(qbmm_idx%vs)) @:DEALLOCATE(qbmm_idx%vs)
+ if (allocated(qbmm_idx%rs)) @:DEALLOCATE(qbmm_idx%rs)Also applies to: 995-996
Comment on lines
+939
to
+950
| local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i & | ||
| & + 1, j, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, & | ||
| & k))/(2._wp*dx) ! force is the negative pressure gradient | ||
| local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(E_idx + fluid_idx)%sf(i, & | ||
| & j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy) | ||
| local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, & | ||
| & j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy) | ||
| cell_volume = abs(dx*dy) | ||
| ! add the 3D component of the pressure gradient, if we are working in 3 dimensions | ||
| if (num_dims == 3) then | ||
| dz = z_cc(k + 1) - z_cc(k) | ||
| local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(E_idx + fluid_idx)%sf(i, & | ||
| & j, k + 1) - q_prim_vf(E_idx + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz) | ||
| local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E + fluid_idx) & | ||
| & %sf(i, j, k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, & | ||
| & k - 1))/(2._wp*dz) |
Contributor
There was a problem hiding this comment.
Pressure-force loop is reading advection fields as pressure.
Only eqn_idx%E is the pressure slot. For fluid_idx > 0, eqn_idx%E + fluid_idx steps into the advection block, so this loop starts subtracting grad(alpha) into the IBM force in multi-fluid cases. Compute the pressure gradient once from q_prim_vf(eqn_idx%E) here, or switch to the actual per-fluid pressure index if that was the intent.
| do i = 1, num_fluids | ||
| alpha_rho(i) = q_cons_vf(i)%sf(j, k, l) | ||
| alpha(i) = q_cons_vf(E_idx + i)%sf(j, k, l) | ||
| alpha(i) = q_cons_vf(eqn_idx%E + i)%sf(j, k, l) |
Contributor
There was a problem hiding this comment.
AMD fallback still overruns alpha_rho/alpha at 3 fluids.
This loop runs to num_fluids, but under the USING_AMD branch above both arrays are still length 2. Since the AMD startup path allows num_fluids == 3, the third iteration writes past both buffers here.
💡 Minimal fix
- #:if not MFC_CASE_OPTIMIZATION and USING_AMD
- real(wp), dimension(2) :: alpha_rho, alpha
+ #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+ real(wp), dimension(3) :: alpha_rho, alphaBased on learnings: when compiling with USING_AMD and not MFC_CASE_OPTIMIZATION, similar simulation arrays are fixed at length 3, and s_check_amd allows num_fluids <= 3.
📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
Suggested change
| alpha(i) = q_cons_vf(eqn_idx%E + i)%sf(j, k, l) | |
| #:if not MFC_CASE_OPTIMIZATION and USING_AMD | |
| real(wp), dimension(3) :: alpha_rho, alpha |
Comment on lines
1821
to
1826
| if (cyl_coord) then | ||
| do i = 1, num_dims | ||
| @:DEALLOCATE(tau_re_vf(cont_idx%end + i)%sf) | ||
| @:DEALLOCATE(tau_re_vf(eqn_idx%cont%end + i)%sf) | ||
| end do | ||
| @:DEALLOCATE(tau_re_vf(e_idx)%sf) | ||
| @:DEALLOCATE(tau_re_vf(eqn_idx%E)%sf) | ||
| @:DEALLOCATE(tau_re_vf) |
Contributor
There was a problem hiding this comment.
Deallocate tau_Re_vf for every viscous run, not only cylindrical ones.
tau_Re_vf is allocated under if (viscous) in s_initialize_rhs_module, but this finalizer frees it only when cyl_coord is true. Cartesian viscous runs will leak these buffers.
Minimal fix
- if (cyl_coord) then
+ if (viscous) then
do i = 1, num_dims
@:DEALLOCATE(tau_re_vf(eqn_idx%cont%end + i)%sf)
end do
@:DEALLOCATE(tau_re_vf(eqn_idx%E)%sf)
@:DEALLOCATE(tau_re_vf)
end ifAs per coding guidelines "If an array is allocated inside an if block, its deallocation must follow the same condition to avoid memory leaks."
Comment on lines
+449
to
+450
| G_L = G_L*max((1._wp - qL_prim_rs${XYZ}$_vf(j, k, l, eqn_idx%damage)), 0._wp) | ||
| G_R = G_R*max((1._wp - qR_prim_rs${XYZ}$_vf(j, k, l, eqn_idx%damage)), 0._wp) |
Contributor
There was a problem hiding this comment.
Use j + 1 for the right-state damage sample.
G_R is computed from the right Riemann state everywhere else, but these reads take eqn_idx%damage from j instead of j + 1. That applies the left-cell damage to the right shear modulus and perturbs the elastic energy/wave-speed calculation across damaged interfaces.
Proposed fix
- G_R = G_R*max((1._wp - qR_prim_rs${XYZ}$_vf(j, k, l, eqn_idx%damage)), 0._wp)
+ G_R = G_R*max((1._wp - qR_prim_rs${XYZ}$_vf(j + 1, k, l, eqn_idx%damage)), 0._wp)Based on learnings: Riemann solver indexing: left states at index j, right states at j+1; off-by-one errors corrupt fluxes.
Also applies to: 1129-1130
Comment on lines
+922
to
+928
| @:ALLOCATE(qbmm_idx%rs(nb), qbmm_idx%vs(nb)) | ||
| @:ALLOCATE(qbmm_idx%ps(nb), qbmm_idx%ms(nb)) | ||
|
|
||
| gam = bub_pp%gam_g | ||
|
|
||
| if (qbmm) then | ||
| @:ALLOCATE(bub_idx%moms(nb, nmom)) | ||
| @:ALLOCATE(qbmm_idx%moms(nb, nmom)) |
Contributor
There was a problem hiding this comment.
1. qbmm_idx allocations not deallocated 📘 Rule violation ☼ Reliability
qbmm_idx allocatables are allocated via @:ALLOCATE(...) but there is no corresponding @:DEALLOCATE(...) for these members, which breaks the project’s required allocation/deallocation pairing and risks leaks across repeated init/finalize cycles. Add matching deallocations in the module finalization path (and any early-return/error paths if applicable).
Agent Prompt
## Issue description
`qbmm_idx` allocatable members are allocated with `@:ALLOCATE(...)` but are never freed with `@:DEALLOCATE(...)`, violating the required pairing and risking leaks when initialization/finalization occurs more than once.
## Issue Context
Allocations occur in `s_initialize_global_parameters_module` (or equivalent init logic) for `qbmm_idx%rs`, `qbmm_idx%vs`, `qbmm_idx%ps`, `qbmm_idx%ms`, and (when `qbmm`) `qbmm_idx%moms`. The module already has a `s_finalize_global_parameters_module` where other `@:DEALLOCATE(...)` calls are performed, but it does not currently deallocate `qbmm_idx` members.
## Fix Focus Areas
- src/simulation/m_global_parameters.fpp[922-952]
- src/simulation/m_global_parameters.fpp[995-1012]
- src/simulation/m_global_parameters.fpp[1299-1340]
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
Comment on lines
+1187
to
+1202
| momxb = eqn_idx%mom%beg | ||
| momxe = eqn_idx%mom%end | ||
| advxb = eqn_idx%adv%beg | ||
| advxe = eqn_idx%adv%end | ||
| contxb = eqn_idx%cont%beg | ||
| contxe = eqn_idx%cont%end | ||
| bubxb = eqn_idx%bub%beg | ||
| bubxe = eqn_idx%bub%end | ||
| strxb = eqn_idx%stress%beg | ||
| strxe = eqn_idx%stress%end | ||
| intxb = eqn_idx%int_en%beg | ||
| intxe = eqn_idx%int_en%end | ||
| xibeg = eqn_idx%xi%beg | ||
| xiend = eqn_idx%xi%end | ||
| chemxb = eqn_idx%species%beg | ||
| chemxe = eqn_idx%species%end |
Contributor
There was a problem hiding this comment.
3. Uninitialized eqn_idx fields read 🐞 Bug ☼ Reliability
eqn_idx sub-fields like stress/xi/B/species are only assigned under feature flags, but multiple targets read them unconditionally (e.g., copying to strxb/xibeg or writing indices.dat ranges). When a feature is disabled, this is a read of undefined values and can produce incorrect indices.dat output or propagate garbage indices into downstream logic.
Agent Prompt
### Issue description
Several `eqn_idx` members (e.g., `stress`, `xi`, `B`, `int_en`, `species`, and scalars like `c`, `damage`, `psi`) are conditionally assigned, but are read unconditionally later (copied into `strxb/strxe/...` and written to `indices.dat`). This reads undefined values when the related feature is disabled.
### Issue Context
- In `m_global_parameters` (simulation/pre_process/post_process), optional ranges are only set in feature-gated blocks.
- Later code unconditionally does `strxb = eqn_idx%stress%beg` etc.
- `m_data_output` (pre_process) unconditionally calls `write_range(eqn_idx%stress%beg, ...)`, `write_range(eqn_idx%xi%beg, ...)`, etc.
### Fix Focus Areas
- src/simulation/m_global_parameters.fpp[1050-1105]
- src/simulation/m_global_parameters.fpp[1187-1206]
- src/pre_process/m_global_parameters.fpp[730-806]
- src/pre_process/m_data_output.fpp[641-652]
### What to change
1. At the start of the equation-index construction routine in each target’s `m_global_parameters`, explicitly initialize all `eqn_idx` members to a known inactive state:
- Set every `int_bounds_info` range’s `%beg`/`%end` to `0`.
- Set scalar indices like `E`, `n`, `gamma`, `pi_inf`, `c`, `damage`, `psi` to `0` (or the established sentinel if the codebase uses one consistently).
2. Only then compute/overwrite the active indices under feature flags.
3. Optionally encapsulate this in a small helper like `call s_reset_eqn_idx(eqn_idx)` to avoid duplicated initialization across targets.
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
- Add GPU_UPDATE for bc_x/y/z%beg/end so CBC device kernels see correct BC codes (was broken: old code extracted to scalars + updated, new code accessed struct members that were never synced to device) - Add matching @:DEALLOCATE/@deallocate for qbmm_idx%rs/vs/ps/ms/moms in s_finalize_global_parameters_module for all three targets - Remove GPU_DECLARE(create='[qbmm_idx]'): struct has allocatable members and is never accessed directly on device; members are always copied to local arrays before GPU use - Remove redundant eqn_idx%n from secondary GPU_UPDATE (already covered by full eqn_idx struct update in s_initialize_global_parameters_module) - Remove stale claude-code-review.yml from repo root (live copy is in .github/workflows/)
This was referenced Apr 10, 2026
Open
Single-line 'if (qbmm) @:DEALLOCATE(...)' is invalid because the macro expands to multiple statements. Use explicit then/end if block instead.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1365 +/- ##
=======================================
Coverage 64.90% 64.90%
=======================================
Files 70 70
Lines 18254 18259 +5
Branches 1508 1505 -3
=======================================
+ Hits 11847 11851 +4
- Misses 5444 5448 +4
+ Partials 963 960 -3 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
E_idx, stress_idx%beg, c_idx, B_idx%beg/end were renamed as part of the eqn_idx refactor but QPVF_IDX_VARS in case.py (used to generate Fortran index expressions for case-optimization builds) was not updated, causing pre_process build failures in case-opt CI jobs.
…ferences in QPVF_IDX_VARS
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #828.
Summary
BC scalars (
bc_xyz_info)bc_xyz_infoderived type groupingbc_x/y/zasbc%x/y/zs_populate_capillary_buffersnow takes an explicitbcargument instead of implicit globalsbcxb/bcxe/bcyb/bcye/bczb/bczefromm_cbc.fppbc${XYZ}$b/bc${XYZ}$ewithbc_${XYZ}$%beg/bc_${XYZ}$%endGPU_UPDATEforbc_x/y/z%beg/endso CBC device kernels see correct valuesUnified equation index type (
eqn_idx)eqn_idx_infoderived type inm_derived_types.fppconsolidating all equation indices:int_bounds_info):cont,mom,adv,bub,stress,xi,B,int_en,speciesinteger):E,n,alf,gamma,pi_inf,c,damage,psicont_idx,mom_idx,adv_idx,bub_idx,stress_idx,xi_idx,B_idx,internalEnergies_idx,species_idx) with a singleeqn_idxin all three targetseqn_idxinstead of 6+ separate scalars and rangeseqn_idxcontains no allocatable members — safe forGPU_DECLAREas a single struct in all casesQBMM index split (
qbmm_idx)rs,vs,ps,ms,moms,fullmom) frombub_bounds_info, leaving it with justbeg/end(nowint_bounds_info)qbmm_idx_infotype andqbmm_idxmodule variable for QBMM moment index mappings with its own@:ALLOCATE/@:DEALLOCATElifecycle in all three targetsGPU_DECLAREremoved forqbmm_idx— members are always accessed via host-side local copies before GPU use, so no device struct neededindices.datimprovements (pre_process)1withnewunit=iudo i = E, Eloop with a guarded writeeqn_idxfields directly via awrite_rangeinternal helperxi,B, and color functionTest plan
./mfc.sh precheck -j 8— all 6 checks pass./mfc.sh build -j 8 --no-gpu— all 3 targets build./mfc.sh build -j 8 --gpu acc— all 3 targets build with nvfortran OpenACC