Volume filtered Euler-Lagrange framework for particle-laden flow explicit unclosed term computation #1043
Conversation
2hrs for Frontier CI (MFlowCode#832)
PR Reviewer Guide 🔍Here are some key observations to aid the review process:
|
![qodo-code-review[bot]](https://avatars.githubusercontent.com/in/484649?s=60&v=4){kind=small}
![cubic-dev-ai[bot]](https://avatars.githubusercontent.com/in/1082092?s=60&v=4){kind=small}
There was a problem hiding this comment.
12 issues found across 36 files
Prompt for AI agents (all issues)
Check if these issues are valid — if so, understand the root cause of each and fix them.
<file name="voronoi/gen_voronoi_2D.py">
<violation number="1" location="voronoi/gen_voronoi_2D.py:35">
P2: The Lloyd iteration recomputes Voronoi cells using initial_points, so every iteration uses the original tessellation and ignores updated point positions. Use the current points array so the relaxation actually updates.</violation>
</file>
<file name="src/simulation/m_compute_statistics.fpp">
<violation number="1" location="src/simulation/m_compute_statistics.fpp:92">
P2: The OpenACC setup is applied to the whole `Msn_q_cons_filtered` array inside a per-element loop. This likely skips proper setup for each element. Use the indexed element in the setup call.</violation>
<violation number="2" location="src/simulation/m_compute_statistics.fpp:150">
P2: The OpenACC setup is applied to the whole `stat_q_cons_filtered` array inside a per-element loop. Use the indexed element to ensure each field is initialized for GPU use.</violation>
</file>
<file name="src/post_process/m_data_input.f90">
<violation number="1" location="src/post_process/m_data_input.f90:679">
P2: alt_sys is hardcoded for a 3D/sys_size=7 layout, but the allocated filtered arrays are sized using num_dims and sys_size. In 1D/2D or other sys_size values this miscounts variables and can break MPI IO (wrong NVARS_MOK / loop bounds). Compute alt_sys from num_dims and sys_size instead of fixed constants.</violation>
</file>
<file name="src/pre_process/m_data_output.fpp">
<violation number="1" location="src/pre_process/m_data_output.fpp:840">
P2: Levelset MPI-IO should use mpi_io_p and scale the element count by mpi_io_type; using mpi_p/data_size*num_ibs breaks mixed-precision (half) output sizes.</violation>
<violation number="2" location="src/pre_process/m_data_output.fpp:860">
P2: Levelset_norm MPI-IO should use mpi_io_p and multiply by mpi_io_type; using mpi_p with raw element counts breaks mixed-precision output sizing.</violation>
</file>
<file name="voronoi/gen_voronoi_3D.py">
<violation number="1" location="voronoi/gen_voronoi_3D.py:21">
P2: `lloyd_relaxation_3d` mutates `initial_points` because `points` is just a reference. This makes the later “initial points” plot incorrect. Copy the array before modifying it.</violation>
<violation number="2" location="voronoi/gen_voronoi_3D.py:62">
P2: `os.mkdir` will fail if `../runs` doesn’t already exist. Use `os.makedirs(..., exist_ok=True)` to create parent directories safely.</violation>
</file>
<file name="src/simulation/m_ibm.fpp">
<violation number="1" location="src/simulation/m_ibm.fpp:620">
P1: The bound limit now allows `index` to reach the last valid coordinate index, but the loop still accesses `s_cc(index + 1)`. This can read past the end of `x_cc/y_cc/z_cc` when `index == bound`. Keep the bound at `m/n/p + buff_size - 1` so `index+1` stays in range.</violation>
</file>
<file name="runs/phi01/case.py">
<violation number="1" location="runs/phi01/case.py:35">
P2: `np.loadtxt` returns a 1D array for a single sphere, but the loop indexes `sphere_loc[i, 0]`, which raises an index error. Normalize to 2D so single-sphere files work.</violation>
<violation number="2" location="runs/phi01/case.py:158">
P2: `periodic_forcing` with `ib` requires `fluid_volume_fraction`, but the case dict sets `particle_vf` instead, so validation will fail and the forcing code won’t have the required volume fraction. Rename or replace this entry with `fluid_volume_fraction`.</violation>
</file>
<file name="src/common/m_mpi_common.fpp">
<violation number="1" location="src/common/m_mpi_common.fpp:2276">
P2: Deallocate the scalarfield buffers under the same condition used for allocation. Otherwise runs with periodic_forcing enabled (but volume_filtering_momentum_eqn disabled) will leak these allocations.</violation>
</file>
Since this is your first cubic review, here's how it works:
- cubic automatically reviews your code and comments on bugs and improvements
- Teach cubic by replying to its comments. cubic learns from your replies and gets better over time
- Ask questions if you need clarification on any suggestion
Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 18
Note
Due to the large number of review comments, Critical, Major severity comments were prioritized as inline comments.
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/simulation/m_data_output.fpp (1)
85-123: Add guard for filtered output in serial mode and refactor argument listThe serial I/O path ignores filtered fields even when
q_filtered_wrtis enabled, leading to silent omission from output. Additionally,s_write_data_fileshas 12 parameters, exceeding the ≤6 argument guideline. Bundle the filtered statistics (filtered_fluid_indicator_function, stat_q_cons_filtered, stat_filtered_pressure, stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch) into a derived-type struct, then add a guard to abort if filtered output is requested without parallel I/O support.Guard to prevent silent filtered data loss
if (.not. parallel_io) then + if (q_filtered_wrt) call s_mpi_abort('q_filtered_wrt requires parallel_io or serial filtered I/O support') call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) else call s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, & filtered_fluid_indicator_function, & stat_q_cons_filtered, stat_filtered_pressure, & stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch) end if
🤖 Fix all issues with AI agents
In `@src/post_process/m_data_input.f90`:
- Around line 679-709: Replace the hard-coded 3D-only expression that sets
alt_sys with a computation based on num_dims and sys_size: compute the extra
variables (filtered indicator, stats per conserved variable, pressure stats,
Reynolds stress components, effective viscosity, interfacial momentum exchange)
by deriving their counts from num_dims (e.g., number of velocity components =
num_dims, number of independent Reynolds stress components =
num_dims*(num_dims+1)/2, and any per-component multiplicative factors currently
hard-coded as 4) and then set alt_sys = sys_size + derived_extra_count; update
the assignment that sets alt_sys (and any comments) so all reads use the new
computed value instead of the 3D-specific literal.
In `@src/simulation/m_additional_forcing.fpp`:
- Around line 192-201: The finalizer s_finalize_additional_forcing_module is
missing deallocation of the window arrays allocated during init; add safe
deallocation calls for rho_window, u_window, and eps_window inside
s_finalize_additional_forcing_module (similar to the existing
DEALLOCATE(q_periodic_force) pattern) ensuring you check for ALLOCATED (or
ASSOCIATED for pointers) before DEALLOCATE and also release any device/copy
resources if these arrays were offloaded (match whatever device deallocation
pattern your codebase uses).
- Around line 84-88: The subroutine s_compute_periodic_forcing currently
declares rhs_vf as intent(in) but mutates it when adding forcing terms; change
the parameter intent to intent(inout) in the subroutine signature for rhs_vf so
the declaration matches its usage, and update any matching interface or module
declaration/prototype for s_compute_periodic_forcing to the same intent(inout)
to keep signatures consistent across the codebase.
- Around line 105-113: The kinetic-energy term is hardcoded to q_cons_vf(2:5)
which breaks non-3D or varying-sys_size layouts; replace the explicit indices
with layout-agnostic helpers by summing momentum components using the mom_idx
int_bounds_info (loop from mom_idx%lower to mom_idx%upper, accumulating
q_cons_vf(contxe + current_mom_idx)%sf(i,j,k)**2) and use E_idx
(q_cons_vf(E_idx)%sf(i,j,k)) for the energy term; update the spatial_eps
accumulation to compute (q_cons_vf(E_idx)%sf(...) - 0.5_wp * sum_mom_sq /
q_cons_vf(1)%sf(...)) * dVol so the code works for 1D/2D/multi-fluid layouts.
In `@src/simulation/m_compute_statistics.fpp`:
- Around line 84-93: ACC_SETUP_VFs is being invoked on the entire arrays
Msn_q_cons_filtered and stat_q_cons_filtered instead of the current element,
which can break macro expansion; update the calls to target the element by
replacing @:ACC_SETUP_VFs(Msn_q_cons_filtered) with
@:ACC_SETUP_VFs(Msn_q_cons_filtered(i)) (and similarly use
stat_q_cons_filtered(i)) inside the loop where you allocate vf and sf, and apply
the same change for the other occurrence around lines 142-151 so the macro runs
on each element instance.
- Around line 239-247: The formulas in the GPU loop of m_compute_statistics.fpp
can divide by zero when ns <= 3 or when the variance Msn(2)%sf(i,j,k) is zero;
modify the loop (the block computing q_stat(1..4)%sf(i,j,k) from Msn) to first
guard using the local ns and the local variance value: always assign q_stat(1) =
Msn(1), then if ns <= 1 set q_stat(2..4) to a safe default (e.g., 0 or a
sentinel), if ns <= 2 set q_stat(3) default, if ns <= 3 set q_stat(4) default,
and additionally check if Msn(2)%sf(i,j,k) is <= eps (small positive threshold)
and in that case avoid dividing by it by returning defaults for q_stat(2..4) or
using a stabilized expression; ensure these checks are done per element inside
the same collapsed GPU loop so no invalid divides occur for any (i,j,k).
In `@src/simulation/m_data_output.fpp`:
- Around line 848-852: The code uses a hard-coded constant (9*4 + 9*4 + 3*4 +
6*4) when setting alt_sys, which assumes num_dims==3 and a fixed sys_size
layout; change this to compute the extra variable count from num_dims and
sys_size and use that computed value in the alt_sys assignment. Specifically,
replace the magic constant with a computed extra_count (e.g., sum the actual
per-variable-array sizes derived from num_dims and sys_size for the arrays
represented by those terms) and set alt_sys = sys_size + 1 + extra_count in the
same conditional branch that checks beta, q_filtered_wrt, t_step, and
t_step_stop so the offset is correct for 1D/2D/3D and different equation sets.
- Around line 957-963: The code only calls s_initialize_mpi_data_filtered inside
the ib branch, leaving filtered buffers uninitialized when q_filtered_wrt is
true but ib is false or when file_per_process=true; update the logic to either
(A) add a validation that if q_filtered_wrt is true and parallel I/O is used
(file_per_process true) then ib must be true and raise an error/abort early, or
(B) move the call to s_initialize_mpi_data_filtered out of the ib-only block so
it runs whenever q_filtered_wrt is true (guarded by the same
file_per_process/alt_sys checks currently used to allocate alt_sys), ensuring
s_initialize_mpi_data_filtered is invoked for all code paths that set up
filtered MPI buffers; reference q_filtered_wrt, s_initialize_mpi_data_filtered,
ib, file_per_process and alt_sys when making the change.
In `@src/simulation/m_ibm.fpp`:
- Around line 1313-1317: The deallocation block conditionally skips freeing
levelset%sf and levelset_norm%sf when store_levelset is false, causing a leak
because these arrays are allocated unconditionally in s_initialize_ibm_module;
update the deallocation sequence to always DEALLOCATE levelset%sf and
levelset_norm%sf (remove them from the store_levelset conditional) or match the
allocation logic in s_initialize_ibm_module by only allocating them when
store_levelset is true—modify the code around the DEALLOCATE calls for
ib_markers%sf, levelset%sf, and levelset_norm%sf to ensure deallocation mirrors
allocation (refer to s_initialize_ibm_module, levelset%sf, levelset_norm%sf, and
store_levelset).
- Around line 644-648: Replace the Fortran intrinsic error stop with the
MPI-safe abort by calling s_mpi_abort(...) instead of "error stop";
specifically, in the block that currently prints the diagnostic (using
proc_rank, temp_loc, index, patch_id, dim, dir, i, j, k) and then calls error
stop "Increase buff_size", remove or keep the prints as desired and replace the
error stop with a call to s_mpi_abort that passes a concise diagnostic string
(e.g. include proc_rank, index, buff_size, bound and patch_id) so all ranks are
terminated cleanly; ensure you call the existing s_mpi_abort routine (same
spelling) and pass the message rather than using error stop.
- Around line 554-610: The code path that computes levelset on-the-fly assumes
spherical IBs when store_levelset is false; add a guard before the "else !
compute levelset..." block to detect non-spherical bodies and abort with
s_mpi_abort (or alternatively set store_levelset = .true. and fall back to
stored levelsets). Concretely, check a unique identifier for the IB shape (e.g.,
patch_ib(patch_id)%shape, or if only spheres have patch_ib(patch_id)%radius
available, test that) right before the on-the-fly computation and call
s_mpi_abort("store_levelset=F requires spherical IBs for patch_id=X") if the IB
is not a sphere so we do not silently compute incorrect image points.
In `@src/simulation/m_start_up.fpp`:
- Around line 1669-1671: The statistics module initialized when
volume_filtering_momentum_eqn is enabled (m_compute_statistics) is not being
finalized, causing host/device allocation leaks; add a finalize call alongside
the existing finalize for the FFTW filter by invoking the statistics finalizer
(call s_finalize_statistics_module()) inside the conditional for
volume_filtering_momentum_eqn so that when volume_filtering_momentum_eqn is true
both s_finalize_fftw_explicit_filter_module and s_finalize_statistics_module are
called.
In `@src/simulation/m_volume_filtering.fpp`:
- Around line 776-779: The division by q_cons_filtered(1)%sf(i,j,k) in the
reynolds_stress update can hit zero or be tiny inside particles; fix by guarding
the denominator—use a small floor (e.g. eps_rho) or max(...) against a filtered
fluid indicator before dividing so you never divide by values below eps_rho;
update the line in the loop that updates reynolds_stress(l,q)%sf(i,j,k) to
divide by a safe_denominator (max(q_cons_filtered(1)%sf(i,j,k), eps_rho) or
conditioned on filtered_fluid_indicator_function) to ensure numerical stability.
- Around line 1545-1552: Several places in the subroutine compute
q_cons_filtered(l)%sf(i,j,k) (and similar assignments later) by dividing by
filtered_fluid_indicator_function%sf(i,j,k) without checking for zero; replicate
the protection used in the Reynolds stress computation: compute a safe
denominator = max(filtered_fluid_indicator_function%sf(i,j,k), small_eps) (or
test if > small_eps) and use that safe denominator (or conditionally set the
result to zero) when assigning q_cons_filtered, and apply the same pattern for
the other divisions that reference filtered_fluid_indicator_function%sf to avoid
NaNs; locate the divisions via the symbols q_cons_filtered, data_real_3D_slabz
and filtered_fluid_indicator_function%sf and replace them with the guarded
version.
- Around line 571-590: The grad_fluid_indicator arrays are allocated but never
zeroed before the accumulation loop, so initialize each
grad_fluid_indicator(i)%sf to zero after allocation (or allocate with a zero
source) to avoid reading/adding uninitialized memory; specifically, after the
@:ALLOCATE and @:ACC_SETUP_SFs calls for grad_fluid_indicator in the block that
allocates 1:3 and grad_fluid_indicator(i)%sf, set grad_fluid_indicator(i)%sf =
0.0 (or an equivalent memset/fill) for i = 1,num_dims so the GPU_PARALLEL_LOOP
accumulation over fd_coeffs and fluid_indicator_function uses a defined starting
value.
In `@src/simulation/p_main.fpp`:
- Around line 59-64: Add corresponding finalizer routines for each initializer
called under volume_filtering_momentum_eqn: implement
s_finalize_filtering_kernel(), s_finalize_filtered_fluid_indicator_function(),
and s_finalize_fluid_indicator_gradient() that free/clean resources allocated by
s_initialize_filtering_kernel(),
s_initialize_filtered_fluid_indicator_function(), and
s_initialize_fluid_indicator_gradient(). Then update s_finalize_modules() to
call these new finalizers when volume_filtering_momentum_eqn is enabled; also
move the existing call to s_finalize_fftw_explicit_filter_module() out of the
special-case in m_start_up.fpp and invoke it from s_finalize_modules() under the
same volume_filtering_momentum_eqn guard so all four modules are finalized via
the standard cleanup path.
In `@voronoi/gen_voronoi_2D.py`:
- Around line 33-36: The loop recreates the Voronoi tessellation with the
immutable initial_points instead of the evolving points array, preventing Lloyd
relaxation from progressing; change the call in the iterations loop so
freud.locality.Voronoi().compute(...) uses the current points variable (not
initial_points) — update the voro.compute((box, initial_points)) invocation to
voro.compute((box, points)) so subsequent centroid moves operate on the newly
tessellated cells in functions/variables around voro, voro_data, vertices and
the centroid-update block.
- Around line 21-26: compute_simplex_centroid and compute_simplex_area currently
slice vertices as simplex_vertices[:, :, 0] which pulls coordinate components
instead of per-vertex vectors; change the indexing to extract each vertex vector
as simplex_vertices[:, 0, :], simplex_vertices[:, 1, :], simplex_vertices[:, 2,
:] in both compute_simplex_centroid and compute_simplex_area so centroids are
computed as the mean of the three vertex vectors and area is computed from
cross(v2 - v1, v3 - v1) with axis=1 for the norm.
🟡 Minor comments (7)
voronoi/gen_voronoi_3D.py-81-87 (1)
81-87: Fix spelling to unblock spell-check.The pipeline flags “overlaping” as a typo in both the comment and log message.
✅ Proposed fix
- # check no spheres are overlaping + # check no spheres are overlapping @@ - print(f'spheres overlaping, dist={dist}, spheres #: {i}, {j}') + print(f'spheres overlapping, dist={dist}, spheres #: {i}, {j}')src/common/m_mpi_common.fpp-2275-2284 (1)
2275-2284: Deallocate scalarfield buffers when periodic_forcing is enabled.
Buffers are allocated for periodic_forcing OR volume_filtering_momentum_eqn, but finalized only for volume_filtering_momentum_eqn today.🐛 Proposed fix
- if (volume_filtering_momentum_eqn) then + if (periodic_forcing .or. volume_filtering_momentum_eqn) then @:DEALLOCATE(buff_send_scalarfield) @:DEALLOCATE(buff_recv_scalarfield) end ifsrc/common/m_boundary_common.fpp-289-396 (1)
289-396: Fix deallocation condition to match allocation.The buffers
buff_send_scalarfieldandbuff_recv_scalarfieldare allocated when(periodic_forcing .or. volume_filtering_momentum_eqn)(m_mpi_common.fpp:106), but deallocated only whenvolume_filtering_momentum_eqn(m_mpi_common.fpp:2276). This creates a memory leak ifperiodic_forcingis enabled without volume filtering. The deallocation condition at line 2276 should be updated to match the allocation condition.src/simulation/m_time_steppers.fpp-49-52 (1)
49-52: Remove unused module importm_compute_statistics.The import of
m_volume_filteringis correct and is used at line 610 (s_apply_fourier_filter). However,m_compute_statisticsis imported but none of its exported routines are called in this file. Remove this unused dependency.toolchain/mfc/case_validator.py-360-366 (1)
360-366: Resolve lint warnings (superfluous parens, trailing whitespace).CI reports superfluous parentheses in the
if (periodic_ibs):line and trailing whitespace in the new checks. Please clean these to keep the lint job green.✏️ Suggested fix
- if (periodic_ibs): + if periodic_ibs:Also applies to: 1104-1169
toolchain/mfc/case_validator.py-356-357 (1)
356-357: Update num_ibs error message to match the new limit.The constraint now allows up to 1000, but the message still mentions 10.
✏️ Suggested fix
- self.prohibit(ib and (num_ibs <= 0 or num_ibs > 1000), - "num_ibs must be between 1 and num_patches_max (10)") + self.prohibit(ib and (num_ibs <= 0 or num_ibs > 1000), + "num_ibs must be between 1 and num_patches_max (1000)")src/post_process/m_data_input.f90-243-299 (1)
243-299: Wrap long allocate lines to satisfy the linterStatic analysis flags multiple lines >100 chars in the filtered-array allocations. Please wrap these with continuation lines to avoid CI failures.
✂️ Example wrapping
- allocate (stat_reynolds_stress(i, j)%vf(k)%sf(local_start_idx:end_x, local_start_idx:end_y, local_start_idx:end_z)) + allocate (stat_reynolds_stress(i, j)%vf(k)%sf( & + local_start_idx:end_x, local_start_idx:end_y, local_start_idx:end_z))
🧹 Nitpick comments (18)
voronoi/gen_voronoi_3D.py (1)
61-63: Preferos.makedirs(..., exist_ok=True)for output dir creation.This avoids the
== Falsecomparison (Ruff E712) and is more robust.♻️ Proposed refactor
- if os.path.exists(output_dir) == False: - os.mkdir(output_dir) + os.makedirs(output_dir, exist_ok=True)Please confirm this aligns with your linting/style expectations.
runs/phi01/case.py (2)
1-3: Unused import:mathThe
mathmodule is imported but never used in this script. Consider removing it.import json -import math import numpy as np
34-36: Relative path dependency for sphere locations.The script loads sphere positions from a relative path
'sphere_array_locations.txt'. This assumes the script is executed from its containing directory (runs/phi01/). Consider using__file__to make the path resolution more robust.+import os + +# load initial sphere locations +script_dir = os.path.dirname(os.path.abspath(__file__)) +sphere_loc = np.loadtxt(os.path.join(script_dir, 'sphere_array_locations.txt')) -sphere_loc = np.loadtxt('sphere_array_locations.txt')runs/3d_drag_test/case.py (1)
1-3: Unused import:mathThe
mathmodule is imported but never used.import json -import math import numpy as npsrc/common/m_compute_levelset.fpp (2)
518-519: Unused variable:ierrThe variable
ierris declared but never used in this subroutine. Remove it to avoid compiler warnings about unused variables.- integer :: i, j, k, ix, iy, iz !< Loop index variables - integer :: ierr + integer :: i, j, k, ix, iy, iz !< Loop index variables
526-548: Periodic image loop always iterates over 8 combinations.When the sphere center is not near any domain boundary,
center(d, 2) = center(d, 1)for all dimensions. The inner 2×2×2 loop (lines 562-575) then redundantly computes the same distance 8 times. While correctness is preserved, this is inefficient.Consider either:
- Skipping the inner loop when
center(:,1) == center(:,2)for all dimensions, or- Tracking which dimensions actually have distinct periodic images and only iterating over those.
This is a minor performance concern for a per-timestep levelset computation but could matter with many IBs.
Also applies to: 562-575
runs/3d_periodic_ibs_test/centered/case.py (1)
1-3: Unused import:mathThe
mathmodule is imported but never used.import json -import math import numpy as npruns/3d_periodic_ibs_test/off-center/case.py (1)
1-3: Unused import:mathThe
mathmodule is imported but never used.import json -import math import numpy as npruns/3d_1sphere_filtering/case.py (1)
1-3: Unused import:mathThe
mathmodule is imported but never used.import json -import math import numpy as npsrc/simulation/m_global_parameters.fpp (1)
1208-1210: Avoid hard-coded “109” sizing for filtered stats.
The extra MPI_IO_DATA slots are fixed at 109 in both allocation and teardown. Consider deriving this count from num_dims and sys_size (or centralizing it as a named constant) to keep it synchronized with s_initialize_mpi_data_filtered and to prevent config drift.Also applies to: 1233-1236, 1433-1436
src/common/m_mpi_common.fpp (1)
340-372: Confirm filtered-stat indexing assumptions.
The offsets (37/73/85/105) and alt_sys = sys_size + 109 imply a fixed layout (likely 3D single‑fluid). If other configurations are possible, derive offsets from num_dims/sys_size or assert supported cases to avoid out‑of‑bounds indexing.src/common/m_ib_patches.fpp (1)
617-640: Add ix/iy/iz to the GPU private list.These inner-loop indices are loop-local; include them in
privateto align with GPU loop guidelines.♻️ Suggested fix
- $:GPU_PARALLEL_LOOP(collapse=3, private='[cart_y, cart_z]', copy='[ib_markers_sf]', copyin='[patch_id, center, r2]') + $:GPU_PARALLEL_LOOP(collapse=3, private='[cart_y, cart_z, ix, iy, iz]', copy='[ib_markers_sf]', copyin='[patch_id, center, r2]')As per coding guidelines, declare loop-local variables with
private='[...]'in GPU parallel loop macros.src/simulation/m_ibm.fpp (1)
154-228: Declare loop indices as private in GPU parallel loops.These loops introduce j/k/l but don’t specify
private, which the GPU macros expect.♻️ Suggested fix (apply pattern to all periodic buffer loops)
- $:GPU_PARALLEL_LOOP(collapse=3) + $:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l]')As per coding guidelines, declare loop-local variables with
private='[...]'in GPU parallel loop macros.src/simulation/m_volume_filtering.fpp (5)
333-337: Consider usingnewunitfor file I/O.Hard-coded unit number
100can conflict with other I/O operations. Modern Fortran providesnewunitfor automatic unit assignment.♻️ Suggested improvement
integer :: particle_force_unit ! Add as module variable ! Then in the open statement: open(newunit=particle_force_unit, file='particle_force.bin', status='replace', & form='unformatted', access='stream', action='write')
1029-1038: Subroutine exceeds size and argument limits.
s_filter_batchis approximately 728 lines (limit: 500) and has 8 arguments (limit: 6). Per coding guidelines, consider:
- Grouping related arguments into a derived type
- Extracting repetitive FFT forward/backward blocks into helper subroutines
The X-Y FFT forward sequence (lines 1041-1088, 1091-1136, 1139-1188, etc.) and the Z FFT + convolution + inverse sequence are repeated for each field type. Extracting these into helpers like
s_fft_xy_forwardands_fft_z_convolvewould significantly reduce code size and improve maintainability.
713-721: Subroutine has 8 arguments (limit: 6).Consider grouping the tensor arguments (
reynolds_stress,visc_stress,pres_visc_stress,div_pres_visc_stress) into a derived type to reduce the argument count and improve readability.
863-863: Remove commented-out debug code.This commented-out print statement appears to be a debug artifact and should be removed before merging.
♻️ Proposed fix
write (100) force_glb flush (100) - !print *, force_glb(1, 1) / (0.5_wp * rho_inf_ref * u_inf_ref**2 * pi * patch_ib(1)%radius**2) end if
1-2007: File exceeds 1000-line limit.At ~2007 lines, this module significantly exceeds the 1000-line guideline. Consider splitting into smaller, focused modules:
m_volume_filtering_types— type definitions and public interfacesm_volume_filtering_fft— FFT initialization, forward/backward transforms, transpose routinesm_volume_filtering_compute— stress tensor, Reynolds stress, effective viscosity computationsm_volume_filtering_particles— particle force calculations
| do q = 1, num_dims | ||
| reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) & | ||
| - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k)/q_cons_filtered(1)%sf(i, j, k)) | ||
| end do |
There was a problem hiding this comment.
Potential division by zero on filtered density.
The Reynolds stress computation divides by q_cons_filtered(1)%sf(i, j, k) (filtered density). In regions where the fluid indicator function approaches zero (inside particles), the filtered density could be very small, causing numerical instability or division by zero.
Consider adding a small epsilon floor or checking against filtered_fluid_indicator_function:
real(wp), parameter :: eps_rho = 1.0e-12_wp
...
reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) &
- (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k) &
/ max(q_cons_filtered(1)%sf(i, j, k), eps_rho))🤖 Prompt for AI Agents
In `@src/simulation/m_volume_filtering.fpp` around lines 776 - 779, The division
by q_cons_filtered(1)%sf(i,j,k) in the reynolds_stress update can hit zero or be
tiny inside particles; fix by guarding the denominator—use a small floor (e.g.
eps_rho) or max(...) against a filtered fluid indicator before dividing so you
never divide by values below eps_rho; update the line in the loop that updates
reynolds_stress(l,q)%sf(i,j,k) to divide by a safe_denominator
(max(q_cons_filtered(1)%sf(i,j,k), eps_rho) or conditioned on
filtered_fluid_indicator_function) to ensure numerical stability.
| do i = 0, m | ||
| do j = 0, n | ||
| do k = 0, p | ||
| q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)/(real(Nx*Ny*Nz, wp)*filtered_fluid_indicator_function%sf(i, j, k)) | ||
| end do | ||
| end do | ||
| end do | ||
| $:END_GPU_PARALLEL_LOOP() |
There was a problem hiding this comment.
Division by filtered_fluid_indicator_function without zero-check.
Multiple locations in this subroutine divide by filtered_fluid_indicator_function%sf(i, j, k) (lines 1548, 1597, 1647, 1699). Same concern as with filtered density—this could cause NaN when the fluid indicator approaches zero.
Apply the same protection pattern as suggested for the Reynolds stress computation.
🤖 Prompt for AI Agents
In `@src/simulation/m_volume_filtering.fpp` around lines 1545 - 1552, Several
places in the subroutine compute q_cons_filtered(l)%sf(i,j,k) (and similar
assignments later) by dividing by filtered_fluid_indicator_function%sf(i,j,k)
without checking for zero; replicate the protection used in the Reynolds stress
computation: compute a safe denominator =
max(filtered_fluid_indicator_function%sf(i,j,k), small_eps) (or test if >
small_eps) and use that safe denominator (or conditionally set the result to
zero) when assigning q_cons_filtered, and apply the same pattern for the other
divisions that reference filtered_fluid_indicator_function%sf to avoid NaNs;
locate the divisions via the symbols q_cons_filtered, data_real_3D_slabz and
filtered_fluid_indicator_function%sf and replace them with the guarded version.
|
CodeAnt AI is reviewing your PR. Thanks for using CodeAnt! 🎉We're free for open-source projects. if you're enjoying it, help us grow by sharing. Share on X · |
codeant-ai
bot
added
the
size:XXL
![codeant-ai[bot]](https://avatars.githubusercontent.com/in/646884?s=60&v=4){kind=small}
src/simulation/m_ibm.fpp
Outdated
| subsection_x = ib_markers%sf(i - gp_layers:i + gp_layers, j, k) | ||
| subsection_y = ib_markers%sf(i, j - gp_layers:j + gp_layers, k) | ||
| subsection_z = ib_markers%sf(i, j, k - gp_layers:k + gp_layers) | ||
| if (any(subsection_x == 0) .or. & | ||
| any(subsection_y == 0) .or. & |
There was a problem hiding this comment.
✅ Suggestion: Logic error: the new code checks three 1D cross-sections (subsection_x, subsection_y, subsection_z) instead of the full 3D neighborhood; this can miss zero entries located off those cross-sections and misclassify ghost vs inner points. Restore the original 3D neighborhood test using subsection_3D = ib_markers%sf(i - gp_layers:i + gp_layers, j - gp_layers:j + gp_layers, k - gp_layers:k + gp_layers) and any(subsection_3D == 0). [logic error]
Severity Level: Critical 🚨
- ❌ Ghost-point counting may be incorrect.
- ⚠️ s_ibm_setup ghost allocations affected.
- ❌ Downstream interpolation uses wrong point sets.| subsection_x = ib_markers%sf(i - gp_layers:i + gp_layers, j, k) | |
| subsection_y = ib_markers%sf(i, j - gp_layers:j + gp_layers, k) | |
| subsection_z = ib_markers%sf(i, j, k - gp_layers:k + gp_layers) | |
| if (any(subsection_x == 0) .or. & | |
| any(subsection_y == 0) .or. & | |
| subsection_3D = ib_markers%sf( & | |
| i - gp_layers:i + gp_layers, & | |
| j - gp_layers:j + gp_layers, & | |
| k - gp_layers:k + gp_layers) | |
| if (any(subsection_3D == 0)) then |
Steps of Reproduction ✅
1. Trigger s_ibm_setup() which calls s_find_num_ghost_points(num_gps, num_inner_gps).
s_ibm_setup is the IBM initialization path used during module setup (call present in
s_ibm_setup in src/simulation/m_ibm.fpp).
2. Within s_find_num_ghost_points (src/simulation/m_ibm.fpp), iterate a 3D immersed
boundary region where gp_layers > 0 and the IB geometry is not axis-aligned.
3. Execution reaches the block starting at line 707 where the code samples three 1D slices
(subsection_x/y/z) instead of the full 3D block.
4. If a zero-valued marker exists in the 3D neighborhood but not on any of the three
central cross-sections, the current condition (lines 707–715) will not detect it and will
misclassify a true ghost point as an inner point. Reproduce by constructing an IB with
zeros off the central planes and observing num_gps_out undercount.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/simulation/m_ibm.fpp
**Line:** 707:711
**Comment:**
*Logic Error: Logic error: the new code checks three 1D cross-sections (`subsection_x`, `subsection_y`, `subsection_z`) instead of the full 3D neighborhood; this can miss zero entries located off those cross-sections and misclassify ghost vs inner points. Restore the original 3D neighborhood test using `subsection_3D = ib_markers%sf(i - gp_layers:i + gp_layers, j - gp_layers:j + gp_layers, k - gp_layers:k + gp_layers)` and `any(subsection_3D == 0)`.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.There was a problem hiding this comment.
There are issues if an immersed boundary is periodically wrapped around the domain and near a domain/processor corner. Since the MPI buffer routines only populate buffers in the 3 cardinal directions, buffer regions located off a diagonal are not populated, which can result in cells being marked as outer ghost points when they are not.
There was a problem hiding this comment.
Question:
There are issues if an immersed boundary is periodically wrapped around the domain and near a domain/processor corner. Since the MPI buffer routines only populate buffers in the 3 cardinal directions, buffer regions located off a diagonal are not populated, which can result in cells being marked as outer ghost points when they are not.
Answer:
This looks like a context/assumption mismatch rather than something you can “fix” purely with Qodo Merge settings: the bot’s suggestion focuses on the local 3D neighborhood test (1D slices vs full 3D), while your comment points to a distributed-memory halo/buffer-population limitation (missing diagonal buffers near periodic + corner), which can also cause ghost/inner misclassification even if the neighborhood test is correct.
What you can do via configuration is steer /improve to explicitly consider these HPC/MPI corner cases and avoid over-confident “restore old logic” recommendations unless it can rule out the buffer/halo issue.
Add targeted extra_instructions (and optionally bump suggestion depth) like this:
[pr_code_suggestions]
extra_instructions = """\
When reviewing ghost/inner point classification in IBM/MPI code:
- Consider periodic wrapping and processor-corner cases.
- Explicitly check whether halo/buffer exchange populates diagonal/corner neighbors; if only 3 cardinal directions are exchanged, call out that the neighborhood data may be incomplete.
- Do NOT recommend switching to a full 3D neighborhood `any(subsection_3D==0)` unless you can confirm the 3D block is actually populated on-rank (including corner halos).
- Prefer suggestions that add/adjust MPI exchange (or document the limitation) and/or add tests that reproduce periodic+corner scenarios.
- If the correctness depends on runtime decomposition/halo policy, state the needed assumptions instead of asserting a single “correct” fix.
"""
suggestions_depth = "exhaustive" # optional, for more edge-case coverageIf you maintain repo best-practices, you can also encode this as a project rule so it’s consistently checked:
- Add a
best_practices.mdentry describing “MPI halo exchange must populate corner/diagonal regions (or code must not assume it)” and “periodic wraparound + corner decomposition must be tested”.
If you share which platform you use (GitHub/GitLab/Bitbucket) and whether you run /improve automatically or manually, I can tailor where to put this (wiki .pr_agent.toml vs local .pr_agent.toml) and how to scope it to only the IBM/MPI files.
Relevant Sources:
- https://docs.qodo.ai/qodo-documentation/qodo-merge/tools/improve#extra-instructions-and-best-practices
- https://docs.qodo.ai/qodo-documentation/qodo-merge/usage-guide/configuration_options#local-configuration-file
- https://docs.qodo.ai/qodo-documentation/qodo-merge/usage-guide/configuration_options#wiki-configuration-file
- https://docs.qodo.ai/qodo-documentation/qodo-merge/tools/improve#controlling-suggestions-depth
src/simulation/m_ibm.fpp
Outdated
| subsection_x = ib_markers%sf(i - gp_layers:i + gp_layers, j, k) | ||
| subsection_y = ib_markers%sf(i, j - gp_layers:j + gp_layers, k) | ||
| subsection_z = ib_markers%sf(i, j, k - gp_layers:k + gp_layers) | ||
| if (any(subsection_x == 0) .or. & | ||
| any(subsection_y == 0) .or. & | ||
| any(subsection_z == 0)) then |
There was a problem hiding this comment.
✅ Suggestion: Logic error (same pattern as above) in the ghost-point discovery routine: using three 1D slices to test for neighboring zeros (subsection_x, subsection_y, subsection_z) is not equivalent to checking the full 3D block and will miss zeros off the central cross-sections; restore the full 3D block test with subsection_3D and any(subsection_3D == 0). [logic error]
Severity Level: Critical 🚨
- ❌ Ghost-point detection fails in 3D cases.
- ⚠️ Interpolation coefficients become invalid.
- ❌ Image-point calculations use wrong points.| subsection_x = ib_markers%sf(i - gp_layers:i + gp_layers, j, k) | |
| subsection_y = ib_markers%sf(i, j - gp_layers:j + gp_layers, k) | |
| subsection_z = ib_markers%sf(i, j, k - gp_layers:k + gp_layers) | |
| if (any(subsection_x == 0) .or. & | |
| any(subsection_y == 0) .or. & | |
| any(subsection_z == 0)) then | |
| subsection_3D = ib_markers%sf( & | |
| i - gp_layers:i + gp_layers, & | |
| j - gp_layers:j + gp_layers, & | |
| k - gp_layers:k + gp_layers) | |
| if (any(subsection_3D == 0)) then |
Steps of Reproduction ✅
1. Call s_find_ghost_points(ghost_points, inner_points) during IBM setup/update (called
from s_ibm_setup and s_update_mib in src/simulation/m_ibm.fpp).
2. During the 3D loop (do k = 0, p) the code reaches the section starting at line 799 that
constructs subsection_x/y/z instead of the full 3D block.
3. If an IB neighborhood contains ib_markers%sf == 0 at an off-center location (not on the
three central cross-sections), the if condition at lines 799–804 will evaluate false and a
true ghost point will be missed.
4. Reproduce by preparing a 3D IB configuration with a zero in the 3D neighborhood away
from the central planes, run s_find_ghost_points, and observe that ghost_points is missing
that location (misclassification).Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/simulation/m_ibm.fpp
**Line:** 799:804
**Comment:**
*Logic Error: Logic error (same pattern as above) in the ghost-point discovery routine: using three 1D slices to test for neighboring zeros (`subsection_x`, `subsection_y`, `subsection_z`) is not equivalent to checking the full 3D block and will miss zeros off the central cross-sections; restore the full 3D block test with `subsection_3D` and `any(subsection_3D == 0)`.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.
Nitpicks 🔍
|
|
CodeAnt AI finished reviewing your PR. |
There was a problem hiding this comment.
Actionable comments posted: 1
🤖 Fix all issues with AI agents
In `@voronoi/gen_voronoi_3D.py`:
- Around line 28-37: The current loop over points builds tetrahedra by taking
consecutive vertex triples (in the block that constructs tetrahedra then
computes centroids via compute_tetrahedron_centroid and volumes via
compute_tetrahedron_volume), which can produce invalid tetrahedra; replace that
tetrahedralization with a proper Delaunay decomposition of the vertex set (use
scipy.spatial.Delaunay on vertices[i] to get tetrahedra indices and then compute
centroids/volumes from those tetrahedra), or simpler and preferred, use freud's
precomputed values by replacing the per-cell centroid/volume computation with
voro_data.centroid[i] and voro_data.volumes[i] to get correct weighted centroids
and volumes. Ensure you remove the consecutive-triple logic and adapt downstream
code to consume the Delaunay-based tetrahedra or voro_data fields.
♻️ Duplicate comments (2)
voronoi/gen_voronoi_3D.py (2)
20-22: Avoid mutating the caller’s input array.
points = initial_pointsaliases the input, so the relaxation mutatesinitial_pointsand the “initial points” plot becomes incorrect. Copy the array before iterating.🛠️ Suggested fix
- points = initial_points + points = np.array(initial_points, copy=True)
61-63: Create the output directory recursively.
os.mkdirfails if../runsdoesn’t exist. Useos.makedirs(..., exist_ok=True)instead.🛠️ Suggested fix
- if os.path.exists(output_dir) == False: - os.mkdir(output_dir) + os.makedirs(output_dir, exist_ok=True)
🧹 Nitpick comments (2)
examples/3D_random_sphere_array/case.py (2)
1-3: Unused import:mathThe
mathmodule is imported but never used. All mathematical operations usenumpyinstead.Suggested fix
import json -import math import numpy as np
36-38: Consider adding basic validation for the sphere location file.If
sphere_array_locations.txtis missing or has fewer than 3 columns, the script will fail with an unclear error when accessingsphere_loc[i, 0:3]in the loop below. A brief check could improve usability:Suggested validation
# load initial sphere locations -sphere_loc = np.loadtxt("sphere_array_locations.txt") +sphere_loc = np.loadtxt("sphere_array_locations.txt", ndmin=2) +if sphere_loc.ndim != 2 or sphere_loc.shape[1] < 3: + raise ValueError("sphere_array_locations.txt must contain at least 3 columns (x, y, z)") N_sphere = len(sphere_loc)
|
CodeAnt AI is running Incremental review Thanks for using CodeAnt! 🎉We're free for open-source projects. if you're enjoying it, help us grow by sharing. Share on X · |
codeant-ai
bot
added
size:XXL
|
CodeAnt AI Incremental review completed. |
There was a problem hiding this comment.
3 issues found across 19 files (changes from recent commits).
Prompt for AI agents (all issues)
Check if these issues are valid — if so, understand the root cause of each and fix them.
<file name="src/post_process/m_global_parameters.fpp">
<violation number="1" location="src/post_process/m_global_parameters.fpp:866">
P1: volume_filter_size is only computed for model_eqns == 2, but it is used for MPI allocations whenever q_filtered_wrt is true. If q_filtered_wrt is enabled with a different model, this uses an uninitialized size and can trigger invalid allocations.</violation>
</file>
<file name="src/simulation/m_global_parameters.fpp">
<violation number="1" location="src/simulation/m_global_parameters.fpp:1214">
P1: volume_filter_size is only computed for model_eqns == 2, but it is used for allocations whenever q_filtered_wrt is enabled. For other model_eqns values this leaves volume_filter_size uninitialized and can cause invalid allocation sizes at runtime.</violation>
</file>
<file name="src/simulation/m_additional_forcing.fpp">
<violation number="1" location="src/simulation/m_additional_forcing.fpp:168">
P2: Use total density for the energy forcing denominator. Dividing by only the first continuity field breaks multi-fluid cases and can cause divide-by-zero when that phase is absent.</violation>
</file>
Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.
There was a problem hiding this comment.
Actionable comments posted: 4
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/simulation/m_start_up.fpp (1)
1279-1352: Guard filtered-stat I/O to avoid uninitialized buffers.
q_filtered_wrtcan be enabled withoutvolume_filter_momentum_eqn, yet the save path updates/writesstat_*arrays regardless. That can dereference uninitialized/unallocated buffers or write garbage whent_step < t_step_start_stats. Gate the filtered I/O onvolume_filter_momentum_eqn(and the stats start step).🐛 Proposed fix
- if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop)) then + if (q_filtered_wrt .and. volume_filter_momentum_eqn .and. & + t_step >= t_step_start_stats .and. (t_step == 0 .or. t_step == t_step_stop)) then $:GPU_UPDATE(host='[filtered_fluid_indicator_function%sf]') do i = 1, num_dims do j = 1, num_dims do k = 1, 4 $:GPU_UPDATE(host='[stat_reynolds_stress(i, j)%vf(k)%sf]') $:GPU_UPDATE(host='[stat_eff_visc(i, j)%vf(k)%sf]') end do end do end do do i = 1, num_dims do j = 1, 4 $:GPU_UPDATE(host='[stat_int_mom_exch(i)%vf(j)%sf]') end do end do do i = 1, E_idx do j = 1, 4 $:GPU_UPDATE(host='[stat_q_cons_filtered(i)%vf(j)%sf]') end do end do do i = 1, 4 $:GPU_UPDATE(host='[stat_filtered_pressure(i)%sf]') end do end if @@ - else if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop)) then + else if (q_filtered_wrt .and. volume_filter_momentum_eqn .and. & + t_step >= t_step_start_stats .and. (t_step == 0 .or. t_step == t_step_stop)) then call s_write_data_files(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, save_count, bc_type, & filtered_fluid_indicator_function=filtered_fluid_indicator_function, & stat_q_cons_filtered=stat_q_cons_filtered, stat_filtered_pressure=stat_filtered_pressure, & stat_reynolds_stress=stat_reynolds_stress, stat_eff_visc=stat_eff_visc, stat_int_mom_exch=stat_int_mom_exch)
🤖 Fix all issues with AI agents
In `@src/common/m_mpi_common.fpp`:
- Around line 346-389: The code uses hardcoded index offsets when assigning
MPI_IO_DATA%var (e.g., sys_size + 1, +37, +73, +85, +105) which is fragile and
error-prone; replace those magic numbers by computing named offsets from
sys_size (or define constants) based on the sizes of the preceding blocks
(referencing filtered_fluid_indicator_function%sf, stat_reynolds_stress,
stat_eff_visc, stat_int_mom_exch, stat_q_cons_filtered, stat_filtered_pressure)
so the indexing is explicit and maintainable, and fix the view-creation loop to
start at sys_size + 1 (change do i = sys_size, alt_sys to do i = sys_size + 1,
alt_sys) to avoid overwriting the last base variable's view for
MPI_IO_DATA%view.
In `@src/post_process/m_data_input.f90`:
- Around line 243-301: The long Allocate statements in subroutine
s_allocate_filtered_arrays (e.g., allocations for
filtered_fluid_indicator_function%sf, stat_q_cons_filtered%vf(j)%sf,
stat_reynolds_stress%vf(k)%sf, stat_eff_visc%vf(k)%sf,
stat_int_mom_exch%vf(j)%sf) exceed line-length limits; break each long
allocate(...) line into multiple lines using Fortran continuation (&) after
commas or before the continued token, indent continuations for readability, and
ensure parentheses and array bounds (local_start_idx:end_x,
local_start_idx:end_y, local_start_idx:end_z) remain intact on the continued
lines so the allocate calls compile cleanly and satisfy Fortitude.
In `@src/simulation/m_volume_filtering.fpp`:
- Around line 214-216: The abort message incorrectly references "p" while the
check uses mod(Ny, num_procs); update the s_mpi_abort call in this conditional
(the call to s_mpi_abort inside the block that checks mod(Ny, num_procs) /= 0)
so the message refers to Ny (or "y-dimension" / "number of y points") being
divisible by num_procs rather than "p", e.g. change the error string to mention
Ny/rows or y-dimension divisibility by num_procs.
In `@toolchain/mfc/case_validator.py`:
- Around line 1104-1192: Strip all trailing whitespace characters from the
affected block (inside the methods check_periodic_forcing and
check_volume_filtering) so no lines end with extra spaces; edit the lines in
those functions' bodies and docstrings to remove trailing spaces and re-run lint
to confirm the Pylint trailing-whitespace warnings are resolved.
♻️ Duplicate comments (9)
src/simulation/m_data_output.fpp (3)
1018-1054: Avoid writing the beta slot twice in the filtered MPI loop.
Whenq_filtered_wrtis true andbetais present, the loop writessys_size+1and then the “Correction” block writes it again, risking overlap/corruption.🐛 Suggested fix
else if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop)) then do i = 1, alt_sys + if (present(beta) .and. i == sys_size + 1) cycle var_MOK = int(i, MPI_OFFSET_KIND) disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, MPI_IO_DATA%view(i), & 'native', mpi_info_int, ierr) call MPI_FILE_WRITE_ALL(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, & mpi_io_p, status, ierr) end do
117-123: Serial path drops filtered outputs.
Whenparallel_iois false, the filtered/statistics arguments are ignored, soq_filtered_wrtcannot produce outputs in serial runs. Please extend the serial path or explicitly reject this configuration.💡 Suggested fix (also update s_write_serial_data_files signature)
- if (.not. parallel_io) then - call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) - else + if (.not. parallel_io) then + if (q_filtered_wrt .and. present(filtered_fluid_indicator_function)) then + call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta, & + filtered_fluid_indicator_function, & + stat_q_cons_filtered, stat_filtered_pressure, & + stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch) + else + call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) + end if + else call s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, & filtered_fluid_indicator_function, & stat_q_cons_filtered, stat_filtered_pressure, & stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch) end if
957-963: Filtered MPI buffers can remain uninitialized.
s_initialize_mpi_data_filteredis only called underib, soq_filtered_wrtwithib = .false.(or file-per-process) leaves filtered views/buffers uninitialized. Move the call outside theibbranch or assert unsupported combinations.src/simulation/m_volume_filtering.fpp (3)
585-594: Finite-difference stencil can read out of bounds.
The loop runsi=0..m,j=0..n,k=0..pwhile indexingi+lwithl=-4..4, which can access beyond bounds near edges. Restrict loop bounds or use buffered indices.🐛 Suggested fix (bounds-safe loop)
- $:GPU_PARALLEL_LOOP(collapse=3, copyin='[fd_coeffs]') - do i = 0, m - do j = 0, n - do k = 0, p + $:GPU_PARALLEL_LOOP(collapse=3, copyin='[fd_coeffs]') + do i = 4, m-4 + do j = 4, n-4 + do k = 4, p-4 $:GPU_LOOP(parallelism='[seq]') do l = -4, 4 grad_fluid_indicator(1)%sf(i, j, k) = grad_fluid_indicator(1)%sf(i, j, k) + fd_coeffs(l)*fluid_indicator_function%sf(i + l, j, k)/dx(i) grad_fluid_indicator(2)%sf(i, j, k) = grad_fluid_indicator(2)%sf(i, j, k) + fd_coeffs(l)*fluid_indicator_function%sf(i, j + l, k)/dy(j) grad_fluid_indicator(3)%sf(i, j, k) = grad_fluid_indicator(3)%sf(i, j, k) + fd_coeffs(l)*fluid_indicator_function%sf(i, j, k + l)/dz(k) end do end do end do end do
785-786: Guard against zero filtered density.
q_cons_filtered(1)%sfcan be zero inside particles, causing division-by-zero/NaNs in Reynolds stress.🐛 Suggested fix
+ real(wp), parameter :: eps_rho = 1.0e-12_wp ... - reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) & - - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k)/q_cons_filtered(1)%sf(i, j, k)) + reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) & + - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k) / & + max(q_cons_filtered(1)%sf(i, j, k), eps_rho))
1553-1560: Filter back-transform divides by near-zero fluid indicator.
Several assignments divide byfiltered_fluid_indicator_function%sf, which can be tiny near/inside particles. Add an epsilon floor (or conditional zeroing) for stability.🐛 Suggested fix (apply to all similar divisions)
+ real(wp), parameter :: eps_phi = 1.0e-12_wp + real(wp) :: denom ... - q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/filtered_fluid_indicator_function%sf(i, j, k) + denom = max(filtered_fluid_indicator_function%sf(i, j, k), eps_phi) + q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/denomAlso applies to: 1602-1607, 1652-1656, 1704-1708
src/simulation/p_main.fpp (1)
59-64: Ensure volume-filtering resources are finalized.
Initialization is added here, but the matching finalization calls are not visible. Please confirm these are invoked ins_finalize_modules()whenvolume_filter_momentum_eqnis enabled.#!/bin/bash # Verify volume-filtering finalizers are called in shutdown paths. rg -n "s_finalize_.*filter|s_finalize_fftw_explicit_filter_module" src/simulation rg -n "volume_filter_momentum_eqn" src/simulation/m_start_up.fpp src/simulation/p_main.fppsrc/common/m_ib_patches.fpp (1)
617-633: Avoid shared module state inside the GPU loop (thread-safety).
s_convert_cylindrical_to_cartesian_coordwrites module-levelcart_y/cart_zinside a parallel region. This is the same race previously flagged; please make the conversion thread‑local (e.g., usingf_convert_cyl_to_cartand local temps) instead of module variables.src/common/m_mpi_common.fpp (1)
1604-1613: Slab decomposition validation relies on generic error message.The validation for slab decomposition divisibility is present via the
ierrflag check at lines 1731-1734, but the error message is generic ("Unsupported combination of values..."). Consider adding a more specific error message for the slab decomposition case to help users diagnose the issue. This is a duplicate of a past review suggestion.
🧹 Nitpick comments (6)
src/simulation/m_data_output.fpp (2)
85-115: Reduce argument count ins_write_data_files.
This signature now exceeds the ≤6-argument guideline; consider a derived-type params struct for the filtered/statistics payload. As per coding guidelines, ...
801-814: Reduce argument count ins_write_parallel_data_files.
This signature exceeds the ≤6-argument guideline; bundle the optional filtered/statistics fields into a derived-type params struct. As per coding guidelines, ...src/simulation/m_volume_filtering.fpp (1)
721-729: Trim parameter lists for filtering routines.
s_setup_terms_filteringands_filter_batchexceed the ≤6-argument guideline; wrap these inputs into a derived-type params object for maintainability. As per coding guidelines, ...Also applies to: 1037-1046
src/post_process/m_data_input.f90 (1)
664-669: Reduce argument count via a params struct.
s_read_parallel_filtered_datatakes 8 arguments, exceeding the project guideline (≤6). Consider a small derived type (e.g.,type(mpi_io_params)) to pass the MPI offsets as a bundle.As per coding guidelines, keep routine arguments ≤ 6 by grouping related parameters.
toolchain/mfc/run/case_dicts.py (1)
118-143: Large parameter dictionary expansion for IB patches.Expanding
ib_idrange from 10 to 1000 generates thousands of dictionary entries. While this aligns with the PR objective to support up to 1000 IB patches, consider the impact on:
- Schema validation performance (fastjsonschema compilation)
- Memory footprint of the parameter dictionaries
This may be acceptable for the use case, but worth noting for future optimization if startup time becomes an issue.
src/simulation/m_global_parameters.fpp (1)
1048-1050: Consider documenting the volume_filter_size formula.The formula
sys_size + 1 + 4*(2*num_dims**2 + num_dims + E_idx + 1)accounts for multiple filtered quantities with 4 statistical moments each, but the calculation is not self-documenting. Consider adding a comment explaining the breakdown:
sys_size: base conservative variables+1: filtered fluid indicator function4 * (...): 4 statistical moments for each filtered quantity
There was a problem hiding this comment.
Actionable comments posted: 7
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/simulation/m_data_output.fpp (1)
85-116: Reduce argument count via a params struct.This subroutine now exceeds the ≤6 argument guideline. Consider bundling the filtered/stat fields into a derived-type params struct (or module-level context) to keep the API stable and easier to evolve. As per coding guidelines, routine arguments should be ≤6.
🤖 Fix all issues with AI agents
In `@src/post_process/m_data_input.f90`:
- Around line 251-299: Several allocate(...) statements in
s_allocate_filtered_arrays (e.g. filtered_fluid_indicator_function%sf(...),
stat_q_cons_filtered(i)%vf(j)%sf(...), stat_reynolds_stress(i,j)%vf(k)%sf(...),
stat_eff_visc(i,j)%vf(k)%sf(...), stat_int_mom_exch(i)%vf(j)%sf(...)) exceed 100
chars and must be wrapped using Fortran continuation; split the long extent list
across lines with an ampersand at the end of the previous line and a leading
ampersand on the continued line so each allocate call stays under 100 chars,
preserving the full allocate(...) expression and indices exactly and keeping the
existing loop structure and indices intact.
In `@src/simulation/m_compute_statistics.fpp`:
- Around line 161-169: The subroutine
s_compute_statistics_momentum_unclosed_terms currently exceeds the ≤6-argument
guideline; create a derived type (e.g., type :: momentum_unclosed_params) that
bundles the unclosed-term fields (reynolds_stress, eff_visc, int_mom_exch,
q_cons_filtered, filtered_pressure) and replace those five arguments in
s_compute_statistics_momentum_unclosed_terms with a single intent(in) parameter
of that new type, keeping only the remaining scalar arguments (t_step,
t_step_start_stats) as separate arguments; update all call sites to construct
and pass the new params struct and adjust any internal references to use the
struct members.
In `@src/simulation/m_data_output.fpp`:
- Around line 801-814: The subroutine s_write_parallel_data_files currently
exceeds the ≤6-argument guideline; refactor by introducing a single derived type
(e.g. a params/holder type) to carry the optional filtered/stat fields
(stat_q_cons_filtered, stat_filtered_pressure, stat_reynolds_stress,
stat_eff_visc, stat_int_mom_exch and optionally
filtered_fluid_indicator_function and beta) and replace those separate optional
arguments with one optional instance of that derived type in the
s_write_parallel_data_files signature; update the type definition, adjust the
intent and optional attributes accordingly, and update all call sites to
construct/pass the new params object so the routine arguments are reduced to ≤6
while preserving intent(inout) semantics for the contained fields.
In `@src/simulation/m_global_parameters.fpp`:
- Around line 1048-1050: The computation of volume_filter_size is done too early
(when q_filtered_wrt is true) before additional contributions are added to
sys_size (surface tension, damage, chemistry), which can under-size MPI buffers;
move or recompute the volume_filter_size assignment so it runs after sys_size is
fully finalized (i.e., right before the MPI_IO_DATA allocation), keeping the
same conditional on q_filtered_wrt and using the final sys_size value; update
any duplicate early computations (also present around the other occurrence) to
the same post-finalization location to ensure buffers are sized correctly.
In `@src/simulation/m_start_up.fpp`:
- Around line 1152-1157: Before passing 1._wp/fluid_pp(1)%Re(1) into
s_volume_filter_momentum_eqn, validate that fluid_pp(1)%Re(1) is set and
strictly positive; if it is missing or <= 0, call s_mpi_abort with a clear error
message and stop execution. Update the block guarded by
volume_filter_momentum_eqn (the call to s_volume_filter_momentum_eqn) to perform
this check on fluid_pp(1)%Re(1) and only compute the reciprocal
1._wp/fluid_pp(1)%Re(1) when the value is valid; reference the symbols
fluid_pp(1)%Re(1), s_volume_filter_momentum_eqn and s_mpi_abort when making the
change.
In `@src/simulation/m_volume_filtering.fpp`:
- Around line 214-216: The abort message is misleading: the check compares Ny
against num_procs but the s_mpi_abort text refers to "p"; update the error text
to reference Ny (or "Ny" and "number of ranks") so it matches the condition.
Locate the conditional that uses Ny and num_procs in the routine (the if
statement and the call to s_mpi_abort) and change the message passed to
s_mpi_abort to something like "Volume filtering requires Ny to be divisible by
the number of ranks" (or similar phrasing referencing Ny) so logs accurately
reflect the checked dimension.
In `@toolchain/mfc/case_validator.py`:
- Around line 1104-1192: Trim all trailing whitespace characters from the newly
added validator methods check_periodic_forcing and check_volume_filtering (and
any adjacent lines in that diff); locate occurrences inside those functions
(including lines with calls to self.get(...), self.prohibit(...), and the
docstrings) and remove trailing spaces at line ends so the linter warning is
resolved.
♻️ Duplicate comments (18)
src/simulation/p_main.fpp (1)
25-64: Add finalization for volume-filtering resources.The new initializers are invoked under
volume_filter_momentum_eqn, but there’s no corresponding cleanup in the finalize path. Please add finalizers (or confirms_finalize_modulesalready handles them) to avoid leaking large allocations.src/common/m_mpi_common.fpp (1)
341-389: Fix alt_sys so filtered MPI views cover the full variable range.
alt_sysis set tovolume_filter_size, so the view-creation loop can miss appended variables. It should besys_size + volume_filter_sizeand the loop should start atsys_size + 1.🐛 Proposed fix
- alt_sys = volume_filter_size + alt_sys = sys_size + volume_filter_size ... - do i = sys_size, alt_sys + do i = sys_size + 1, alt_syssrc/post_process/m_data_input.f90 (1)
679-709: Fix alt_sys so MPI filtered reads include base + filtered variables.
alt_sysis set tovolume_filter_size, soNVARS_MOKand the read loop undercount when base variables are present. It should besys_size + volume_filter_size.🐛 Proposed fix
- alt_sys = volume_filter_size ! filtered indicator, stats of: R_u, R_mu, F_imet, q_cons_filtered, pressure + alt_sys = sys_size + volume_filter_size ! filtered indicator + stats appended to base varssrc/simulation/m_ibm.fpp (3)
620-630: Fix off-by-one in image-point search bounds.
index+1is accessed, so allowingindex==m+buff_sizecan read past the coordinate arrays.🐛 Proposed fix
- bound = m + buff_size + bound = m + buff_size - 1 ... - bound = n + buff_size + bound = n + buff_size - 1 ... - bound = p + buff_size + bound = p + buff_size - 1
648-651: Uses_mpi_abortinstead oferror stop.
error stopcan hang MPI ranks; switch to the MPI-aware abort. As per coding guidelines, uses_mpi_abortfor error termination.🐛 Suggested fix
- print *, "Increase buff_size further in m_helper_basic (currently set to a minimum of 10)" - error stop "Increase buff_size" + call s_mpi_abort("Increase buff_size further in m_helper_basic (currently set to a minimum of 10)")
1316-1320: Deallocate levelset fields even whenstore_levelset=F.
They are allocated unconditionally in initialization, so conditional deallocation leaks memory. As per coding guidelines, use the deallocation macro.🐛 Proposed fix
@:DEALLOCATE(ib_markers%sf) - if (store_levelset) then - @:DEALLOCATE(levelset%sf) - @:DEALLOCATE(levelset_norm%sf) - end if + @:DEALLOCATE(levelset%sf) + @:DEALLOCATE(levelset_norm%sf)src/post_process/m_global_parameters.fpp (2)
123-124: Ensurevolume_filter_sizeis initialized for allq_filtered_wrtcases.
It’s computed only undermodel_eqns == 2, but used for allocation wheneverq_filtered_wrtis true; ifq_filtered_wrtis enabled for another model, the size is undefined. Consider computing it aftersys_sizeis finalized for all models, or validate thatq_filtered_wrtis only used with model_eqns=2.♻️ Suggested placement (after sys_size is finalized)
- if (q_filtered_wrt) then - volume_filter_size = sys_size + 1 + 4*(2*num_dims**2 + num_dims + E_idx + 1) - end if ... + if (q_filtered_wrt) then + volume_filter_size = sys_size + 1 + 4*(2*num_dims**2 + num_dims + E_idx + 1) + end ifAlso applies to: 698-700, 865-868
1058-1062: Deallocate extra MPI_IO_DATA entries for filtered outputs.
Nulling pointers alone leaks the allocatedsfbuffers. As per coding guidelines, deallocate with the macro when associated.🐛 Proposed fix
if (q_filtered_wrt) then do i = sys_size + 1, volume_filter_size - MPI_IO_DATA%var(i)%sf => null() + if (associated(MPI_IO_DATA%var(i)%sf)) then + @:DEALLOCATE(MPI_IO_DATA%var(i)%sf) + end if end do end ifsrc/simulation/m_data_output.fpp (3)
117-123: Serial path drops filtered outputs whenparallel_iois false.
q_filtered_wrtdata never reachess_write_serial_data_files, so filtered outputs are silently skipped in serial runs. This matches a prior finding.🔧 Suggested fix
- if (.not. parallel_io) then - call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) + if (.not. parallel_io) then + if (q_filtered_wrt .and. present(filtered_fluid_indicator_function)) then + call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta, & + filtered_fluid_indicator_function, stat_q_cons_filtered, stat_filtered_pressure, & + stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch) + else + call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) + end if else call s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, & filtered_fluid_indicator_function, & stat_q_cons_filtered, stat_filtered_pressure, & stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch) end if
957-963: Filtered MPI buffers are still only initialized inside theibbranch.When
q_filtered_wrtis true butibis false (or whenfile_per_processis true), the filtered MPI buffers/views remain uninitialized. This mirrors a prior finding.
1018-1054: Potential duplicate write ofbetaslot whenq_filtered_wrtis enabled.The
q_filtered_wrtloop writes indices1..alt_sys, which includessys_size+1whenpresent(beta)setsalt_sys = sys_size + 1, and thenbetais written again in the correction block. This is a previously reported overlap risk.src/simulation/m_additional_forcing.fpp (2)
111-115: Momentum-energy term still assumes 3 components.
momxb,momxb+1,momxb+2hard-code a 3D layout and break 1D/2D or nonstandard layouts. This matches a prior review note.
166-168: Energy forcing denominator is single-fluid and unguarded.Division by
q_cons_vf(contxb)can be zero and ignores multi-fluid mixture density, consistent with earlier review findings.src/simulation/m_global_parameters.fpp (1)
1048-1050:volume_filter_sizeonly set in model_eqns==2.If
q_filtered_wrtis enabled under other equation sets,volume_filter_sizecan be uninitialized and lead to invalid allocations. This repeats a prior review finding.src/simulation/m_compute_statistics.fpp (1)
239-246: Moments divide by zero whenns <= 3or variance is zero.The formulas can divide by
ns-1/2/3and byMsn(2)when the variance is zero. This matches a prior review finding.src/simulation/m_volume_filtering.fpp (3)
585-594: Stencil access can go out of bounds at domain edges.
i+l,j+l,k+lwill underflow/overflow near boundaries unless halos are ≥4. Restrict loop ranges or guard indices.
785-786: Guard division by filtered density.
q_cons_filtered(1)%sf(i,j,k)can be ~0 in solid regions, causing NaNs. Add a small floor.🐛 Suggested pattern
- reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) & - - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k)/q_cons_filtered(1)%sf(i, j, k)) + reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) & + - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k) & + / max(q_cons_filtered(1)%sf(i, j, k), eps_rho))
1553-1557: Protect divisions byfiltered_fluid_indicator_function.
The denominator can be 0 in solid cells, leading to NaNs forq_cons_filtered,filtered_pressure,reynolds_stress, andeff_visc. Add an epsilon or conditional guard.🐛 Suggested pattern (apply to all four sites)
- q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/filtered_fluid_indicator_function%sf(i, j, k) + q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm / & + max(filtered_fluid_indicator_function%sf(i, j, k), eps_phi)Also applies to: 1602-1606, 1652-1656, 1704-1708
🧹 Nitpick comments (4)
toolchain/mfc/case_validator.py (1)
1104-1155: Reduce local variable count incheck_periodic_forcingto satisfy pylint.
A small table-driven validation can drop several locals while keeping the same error messages.♻️ Possible refactor sketch
- u_inf_ref = self.get('u_inf_ref') - rho_inf_ref = self.get('rho_inf_ref') - P_inf_ref = self.get('P_inf_ref') - forcing_window = self.get('forcing_window') - forcing_dt = self.get('forcing_dt') + required = { + 'u_inf_ref': "u_inf_ref (desired bulk velocity) must be specified for periodic_forcing", + 'rho_inf_ref': "rho_inf_ref (desired bulk density) must be specified for periodic_forcing", + 'P_inf_ref': "P_inf_ref (desired pressure based on bulk internal energy) must be specified for periodic_forcing", + 'forcing_window': "forcing_window must be specified for periodic_forcing", + 'forcing_dt': "forcing_dt must be specified for periodic_forcing", + } + for key, msg in required.items(): + self.prohibit(self.get(key) is None, msg)src/simulation/m_ibm.fpp (1)
154-228: Add explicit private lists to periodic-buffer GPU loops.
If the macro doesn’t implicitly privatizej/k/l, these loops can race. Please verify and addprivate='[j,k,l]'(apply to all six loops). As per coding guidelines, declare loop-local variables in GPU parallel macros.♻️ Example (apply to all six loops)
- $:GPU_PARALLEL_LOOP(collapse=3) + $:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l]')src/simulation/m_volume_filtering.fpp (2)
721-729: Consider grouping arguments to keep subroutines ≤6 params.
Boths_setup_terms_filteringands_filter_batchexceed 6 arguments. A derived-type params struct would keep signatures manageable. As per coding guidelines, consider consolidating related fields.Also applies to: 1037-1045
92-95: Verify mixed‑precision slab buffers are intentional.
Slab buffers and MPI traffic usecomplex(c_float_complex)/MPI_COMPLEXwhile FFTs operate on double‑complex. This introduces implicit casts and potential precision loss. If double precision is required end‑to‑end, align buffer kinds and MPI datatypes.🔧 Possible alignment (if double precision is desired)
- complex(c_float_complex), allocatable :: data_cmplx_slabz(:, :, :), data_cmplx_slaby(:, :, :) + complex(c_double_complex), allocatable :: data_cmplx_slabz(:, :, :), data_cmplx_slaby(:, :, :) - call MPI_Alltoall(sendbuf_sf, NxC*Nyloc*Nzloc, MPI_COMPLEX, & - recvbuf_sf, NxC*Nyloc*Nzloc, MPI_COMPLEX, MPI_COMM_WORLD, ierr) + call MPI_Alltoall(sendbuf_sf, NxC*Nyloc*Nzloc, MPI_DOUBLE_COMPLEX, & + recvbuf_sf, NxC*Nyloc*Nzloc, MPI_DOUBLE_COMPLEX, MPI_COMM_WORLD, ierr)Also applies to: 895-897, 933-935, 973-975, 1015-1017
| subroutine s_compute_statistics_momentum_unclosed_terms(t_step, t_step_start_stats, reynolds_stress, eff_visc, int_mom_exch, q_cons_filtered, filtered_pressure) | ||
| type(scalar_field), dimension(num_dims, num_dims), intent(in) :: reynolds_stress | ||
| type(scalar_field), dimension(num_dims, num_dims), intent(in) :: eff_visc | ||
| type(scalar_field), dimension(num_dims), intent(in) :: int_mom_exch | ||
| type(scalar_field), dimension(sys_size), intent(in) :: q_cons_filtered | ||
| type(scalar_field), intent(in) :: filtered_pressure | ||
| integer, intent(in) :: t_step | ||
| integer, intent(in) :: t_step_start_stats | ||
|
|
There was a problem hiding this comment.
🛠️ Refactor suggestion | 🟠 Major
Reduce argument count via a params struct.
This routine now exceeds the ≤6 argument guideline; consider packing the unclosed-term fields into a derived-type params struct. As per coding guidelines, routine arguments should be ≤6.
🤖 Prompt for AI Agents
In `@src/simulation/m_compute_statistics.fpp` around lines 161 - 169, The
subroutine s_compute_statistics_momentum_unclosed_terms currently exceeds the
≤6-argument guideline; create a derived type (e.g., type ::
momentum_unclosed_params) that bundles the unclosed-term fields
(reynolds_stress, eff_visc, int_mom_exch, q_cons_filtered, filtered_pressure)
and replace those five arguments in s_compute_statistics_momentum_unclosed_terms
with a single intent(in) parameter of that new type, keeping only the remaining
scalar arguments (t_step, t_step_start_stats) as separate arguments; update all
call sites to construct and pass the new params struct and adjust any internal
references to use the struct members.
| impure subroutine s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, & | ||
| filtered_fluid_indicator_function, & | ||
| stat_q_cons_filtered, stat_filtered_pressure, & | ||
| stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch) | ||
|
|
||
| type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf | ||
| integer, intent(in) :: t_step | ||
| type(scalar_field), intent(inout), optional :: beta | ||
| type(scalar_field), intent(inout), optional :: filtered_fluid_indicator_function | ||
| type(vector_field), dimension(E_idx), intent(inout), optional :: stat_q_cons_filtered | ||
| type(scalar_field), dimension(4), intent(inout), optional :: stat_filtered_pressure | ||
| type(vector_field), dimension(num_dims, num_dims), intent(inout), optional :: stat_reynolds_stress | ||
| type(vector_field), dimension(num_dims, num_dims), intent(inout), optional :: stat_eff_visc | ||
| type(vector_field), dimension(num_dims), intent(inout), optional :: stat_int_mom_exch |
There was a problem hiding this comment.
🛠️ Refactor suggestion | 🟠 Major
Reduce argument count via a params struct.
This signature now exceeds the ≤6 argument guideline. Consider a single derived-type carrying filtered/stat fields to keep signatures concise and consistent. As per coding guidelines, routine arguments should be ≤6.
🤖 Prompt for AI Agents
In `@src/simulation/m_data_output.fpp` around lines 801 - 814, The subroutine
s_write_parallel_data_files currently exceeds the ≤6-argument guideline;
refactor by introducing a single derived type (e.g. a params/holder type) to
carry the optional filtered/stat fields (stat_q_cons_filtered,
stat_filtered_pressure, stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch
and optionally filtered_fluid_indicator_function and beta) and replace those
separate optional arguments with one optional instance of that derived type in
the s_write_parallel_data_files signature; update the type definition, adjust
the intent and optional attributes accordingly, and update all call sites to
construct/pass the new params object so the routine arguments are reduced to ≤6
while preserving intent(inout) semantics for the contained fields.
|
CodeAnt AI is running Incremental review Thanks for using CodeAnt! 🎉We're free for open-source projects. if you're enjoying it, help us grow by sharing. Share on X · |
codeant-ai
bot
added
size:XXL
| do i = 1, num_dims | ||
| do j = 1, 4 | ||
| MPI_IO_DATA%var(sys_size + volume_filter_dt%stat_mom_exch_idx + (i - 1)*4 + (j - 1))%sf => stat_int_mom_exch(i)%vf(j)%sf(0:m, 0:n, 0:p) | ||
| end do |
There was a problem hiding this comment.
Suggestion: Off-by-one error when creating MPI data type views: the loop currently starts at sys_size, which will recreate/override the view for index sys_size (already created earlier) and may leave the intended first new variable without a dedicated view; change the loop to start at sys_size + 1. [off-by-one]
Severity Level: Critical 🚨
- ❌ MPI parallel I/O writes mismatched fields.
- ⚠️ Volume-filtered statistics output corrupted.
- ⚠️ Post-processing reads incorrect variables.| end do | |
| ! Define the view for each newly added variable (avoid recreating existing sys_size view) | |
| do i = sys_size + 1, alt_sys |
Steps of Reproduction ✅
1. Build/run the code with MPI enabled and volume-filter features enabled so filtered
initialization runs.
2. During startup the program executes s_initialize_mpi_data (src/common/m_mpi_common.fpp)
which earlier defines MPI views for 1..alt_sys (see the loop creating MPI_IO_DATA%view(i)
for i = 1, alt_sys in the same file).
3. Later s_initialize_mpi_data_filtered (src/common/m_mpi_common.fpp) executes the shown
loop (lines starting at 366) which begins at i = sys_size and therefore
recreates/overwrites the already-committed view at index sys_size instead of creating the
next new view.
4. Result: the view-to-variable mapping becomes misaligned for subsequent MPI I/O calls
(MPI_IO_DATA%view(index) no longer corresponds to the intended variable). This is
observable as incorrect fields being written/read during parallel IO at
initialization/write phases that use MPI_IO_DATA%view.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/common/m_mpi_common.fpp
**Line:** 366:366
**Comment:**
*Off By One: Off-by-one error when creating MPI data type views: the loop currently starts at `sys_size`, which will recreate/override the view for index `sys_size` (already created earlier) and may leave the intended first new variable without a dedicated view; change the loop to start at `sys_size + 1`.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.| file_loc = trim(t_step_dir)//'/levelset.dat' | ||
|
|
||
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | ||
| write (1) levelset%sf | ||
| close (1) | ||
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | ||
| write (1) levelset%sf | ||
| close (1) |
There was a problem hiding this comment.
Suggestion: Missing presence check for optional levelset: the subroutine accepts levelset as OPTIONAL but the new code writes levelset%sf unconditionally; if the caller did not pass levelset this will cause a runtime error. Wrap the write in if (present(levelset)) then ... end if. [null pointer]
Severity Level: Critical 🚨
- ❌ Serial data output can crash during levelset file writing.
- ⚠️ Pre-processing step writing initial state affected.| file_loc = trim(t_step_dir)//'/levelset.dat' | |
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | |
| write (1) levelset%sf | |
| close (1) | |
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | |
| write (1) levelset%sf | |
| close (1) | |
| if (present(levelset)) then | |
| file_loc = trim(t_step_dir)//'/levelset.dat' | |
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | |
| write (1) levelset%sf | |
| close (1) | |
| end if |
Steps of Reproduction ✅
1. Inspect subroutine s_write_serial_data_files in src/pre_process/m_data_output.fpp (the
subroutine signature includes optional arguments `levelset, levelset_norm`). Execution
enters the block guarded by ib and store_levelset at the lines shown in the PR hunk (if
(store_levelset) then ... — the levelset write begins at line 226).
2. Ensure the code executes the serial data-output path: run with parallel_io disabled so
s_write_data_files => s_write_serial_data_files is used (assignment occurs in
s_initialize_data_output_module). This leads to the serial writer being invoked.
3. Trigger a run where ib and store_levelset are true (so the block at lines 224–237 is
entered), but the caller of s_write_serial_data_files did not pass the optional argument
`levelset` (i.e., present(levelset) is .false.). When execution reaches line 229 (write
(1) levelset%sf) the code dereferences the absent optional argument.
4. Observe a Fortran runtime error (attempt to reference absent OPTIONAL actual) or crash
at run-time when write (1) levelset%sf executes (failure point:
src/pre_process/m_data_output.fpp lines ~226-230). If callers always supply `levelset`
this path won't fail — but the current code makes that assumption unguarded.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/pre_process/m_data_output.fpp
**Line:** 226:230
**Comment:**
*Null Pointer: Missing presence check for optional `levelset`: the subroutine accepts `levelset` as OPTIONAL but the new code writes `levelset%sf` unconditionally; if the caller did not pass `levelset` this will cause a runtime error. Wrap the write in `if (present(levelset)) then ... end if`.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.| file_loc = trim(t_step_dir)//'/levelset_norm.dat' | ||
|
|
||
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | ||
| write (1) levelset_norm%sf | ||
| close (1) | ||
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | ||
| write (1) levelset_norm%sf | ||
| close (1) |
There was a problem hiding this comment.
Suggestion: Missing presence check for optional levelset_norm: the subroutine accepts levelset_norm as OPTIONAL but the new code writes levelset_norm%sf unconditionally; this will fail at runtime if levelset_norm was not provided. Guard the write with if (present(levelset_norm)) then ... end if. [null pointer]
Severity Level: Critical 🚨
- ❌ Serial data output can crash when writing levelset normals.
- ⚠️ Initial condition levelset normalization output unreliable.| file_loc = trim(t_step_dir)//'/levelset_norm.dat' | |
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | |
| write (1) levelset_norm%sf | |
| close (1) | |
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | |
| write (1) levelset_norm%sf | |
| close (1) | |
| if (present(levelset_norm)) then | |
| file_loc = trim(t_step_dir)//'/levelset_norm.dat' | |
| open (1, FILE=trim(file_loc), FORM='unformatted', STATUS=status) | |
| write (1) levelset_norm%sf | |
| close (1) | |
| end if |
Steps of Reproduction ✅
1. In src/pre_process/m_data_output.fpp, locate the same s_write_serial_data_files
subroutine and the store_levelset block added in the PR (the levelset_norm write starts at
line 232 in the new hunk).
2. Execute the serial data-output path (s_write_data_files => s_write_serial_data_files)
with ib and store_levelset enabled so the code enters the block encompassing lines
224–237.
3. If the caller did not supply the optional argument `levelset_norm`
(present(levelset_norm) is .false.), execution reaches line 235 (write (1)
levelset_norm%sf) and dereferences the absent optional argument.
4. Observe a runtime failure or crash at the write of levelset_norm
(src/pre_process/m_data_output.fpp lines ~232-236). If all callers are already passing
levelset_norm this code will not fail, but the subroutine signature allows omission so the
unguarded access is unsafe.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/pre_process/m_data_output.fpp
**Line:** 232:236
**Comment:**
*Null Pointer: Missing presence check for optional `levelset_norm`: the subroutine accepts `levelset_norm` as OPTIONAL but the new code writes `levelset_norm%sf` unconditionally; this will fail at runtime if `levelset_norm` was not provided. Guard the write with `if (present(levelset_norm)) then ... end if`.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.| if (store_levelset) then | ||
| allocate (levelset%sf(0:m, 0:n, 0:p, 1:num_ibs)) | ||
| allocate (levelset_norm%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)) | ||
| end if |
There was a problem hiding this comment.
Suggestion: The new code allocates levelset%sf and levelset_norm%sf only when store_levelset is true but leaves them uninitialized; if other code expects these arrays to exist (or reads them) this will cause unallocated/undefined-data runtime errors — allocate them unconditionally (or ensure caller is guarded) and initialize them to zero right after allocation to avoid uninitialized reads. [possible bug]
Severity Level: Critical 🚨
- ❌ Initialization crash when IB enabled without stored levelsets.
- ⚠️ Immersed-boundary setup fails preventing IB features.
- ❌ Job aborts during start-up for IB-enabled cases.| if (store_levelset) then | |
| allocate (levelset%sf(0:m, 0:n, 0:p, 1:num_ibs)) | |
| allocate (levelset_norm%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)) | |
| end if | |
| allocate (levelset%sf(0:m, 0:n, 0:p, 1:num_ibs)) | |
| allocate (levelset_norm%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)) | |
| levelset%sf = 0 | |
| levelset_norm%sf = 0 |
Steps of Reproduction ✅
1. Build and run the code that calls module m_initial_condition (file
src/pre_process/m_initial_condition.fpp).
2. During startup, s_initialize_initial_condition_module allocates IB bookkeeping: see
allocation block at lines 104-107 in this file (if (store_levelset) then ... end if). If
store_levelset is false, levelset%sf and levelset_norm%sf are not allocated.
3. Later the initialization sequence invokes s_generate_initial_condition which (when ib
is true) calls s_apply_ib_patches(ib_markers%sf, levelset, levelset_norm). The call site
is in s_generate_initial_condition in the same file (function name
s_generate_initial_condition; call present in the subroutine body).
4. With runtime configuration where ib = .true. and store_levelset = .false.,
s_apply_ib_patches receives unallocated components levelset%sf and/or levelset_norm%sf and
will trigger a Fortran runtime error (unallocated array or invalid memory access) when
those components are accessed.
Note: This is not a hypothetical path — the allocation decision is at lines 104-107 and
the call into s_apply_ib_patches is unconditional for ib in s_generate_initial_condition;
therefore leaving levelset arrays unallocated is a concrete reproducible failure scenario.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/pre_process/m_initial_condition.fpp
**Line:** 104:107
**Comment:**
*Possible Bug: The new code allocates `levelset%sf` and `levelset_norm%sf` only when `store_levelset` is true but leaves them uninitialized; if other code expects these arrays to exist (or reads them) this will cause unallocated/undefined-data runtime errors — allocate them unconditionally (or ensure caller is guarded) and initialize them to zero right after allocation to avoid uninitialized reads.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.|
|
||
| # load initial sphere locations | ||
| sphere_loc = np.loadtxt("sphere_array_locations.txt") | ||
| N_sphere = len(sphere_loc) |
There was a problem hiding this comment.
Suggestion: np.loadtxt returns a 1-D array for a single-row file which makes subsequent indexing like sphere_loc[i, 0] crash; convert the result to a 2-D array and compute the number of spheres from shape[0] so single-row and multi-row inputs both work. [possible bug]
Severity Level: Critical 🚨
- ❌ Case generation crashes for single-row sphere file.
- ❌ Immersed-boundary dictionary construction fails.
- ⚠️ Single-sphere simulation setup blocked.| N_sphere = len(sphere_loc) | |
| sphere_loc = np.atleast_2d(sphere_loc) | |
| N_sphere = sphere_loc.shape[0] |
Steps of Reproduction ✅
1. Create a single-line locations file:
`examples/3D_random_sphere_array/sphere_array_locations.txt` containing "0.0 0.0 0.0".
2. Run the script: `python examples/3D_random_sphere_array/case.py` (top-level script,
executes from file start; see print at line 164).
3. Script executes load at line 37 (`sphere_loc = np.loadtxt(...)`) producing a 1-D array
of shape (3,) instead of (1,3).
4. Loop at line 42 (`for i in range(N_sphere):`) computes N_sphere = len(sphere_loc) at
line 38 as 3, then at line 46 the code uses `sphere_loc[i, 0]`, causing an IndexError: too
many indices for array when i==0 (access at lines 46-48).
5. Observe script crash before JSON case output at line 164; immersed-boundary dictionary
`ib_dict` construction fails.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** examples/3D_random_sphere_array/case.py
**Line:** 38:38
**Comment:**
*Possible Bug: np.loadtxt returns a 1-D array for a single-row file which makes subsequent indexing like `sphere_loc[i, 0]` crash; convert the result to a 2-D array and compute the number of spheres from `shape[0]` so single-row and multi-row inputs both work.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.| "patch_icpp(1)%alpha_rho(1)": rho, | ||
| "patch_icpp(1)%alpha(1)": 1.0e00, | ||
| # Fluids Physical Parameters | ||
| "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00), # 2.50(Not 1.40) |
There was a problem hiding this comment.
Suggestion: The value assigned to "fluid_pp(1)%gamma" is 1.0 / (gam_a - 1.0) which equals 1/(gamma-1) (e.g. 2.5 for gamma=1.4); if the field expects the physical ratio of specific heats it should be gam_a instead — the current assignment likely mixes up gamma with (gamma-1) and is a logic error. [logic error]
Severity Level: Critical 🚨
- ❌ Case JSON contains incorrect gamma value.
- ⚠️ Thermodynamic properties computed incorrectly.
- ⚠️ Simulation equation-of-state behavior altered.| "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00), # 2.50(Not 1.40) | |
| "fluid_pp(1)%gamma": gam_a, |
Steps of Reproduction ✅
1. Open `examples/3D_random_sphere_array/case.py` and note gamma variable at line 5
(`gam_a = 1.4`).
2. Run `python examples/3D_random_sphere_array/case.py` to print the case JSON (print at
line 164).
3. Inspect printed JSON field `fluid_pp(1)%gamma` populated from line 140; it will be `1.0
/ (gam_a - 1.0)` → 2.5 for gam_a=1.4.
4. Confirm that the produced value differs from the physical specific-heat ratio
(gam_a=1.4) by inspecting lines 5 and 140. This shows the script sets 2.5 into
`fluid_pp(1)%gamma` rather than 1.4.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** examples/3D_random_sphere_array/case.py
**Line:** 140:140
**Comment:**
*Logic Error: The value assigned to `"fluid_pp(1)%gamma"` is `1.0 / (gam_a - 1.0)` which equals 1/(gamma-1) (e.g. 2.5 for gamma=1.4); if the field expects the physical ratio of specific heats it should be `gam_a` instead — the current assignment likely mixes up gamma with (gamma-1) and is a logic error.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.| "P_inf_ref": P, | ||
| "mom_f_idx": 1, | ||
| "forcing_window": 1, | ||
| "forcing_dt": 1.0 / (0.5 * dt), |
There was a problem hiding this comment.
Suggestion: forcing_dt is computed as the reciprocal of 0.5 * dt (i.e. 2/dt) which yields an extremely large value for small dt; it's almost certainly intended to be half the timestep (0.5 * dt) rather than its inverse. [logic error]
Severity Level: Critical 🚨
- ❌ Forcing timestep misconfigured in case JSON.
- ⚠️ Periodic forcing frequency will be incorrect.
- ⚠️ Simulation bulk forcing behavior altered.| "forcing_dt": 1.0 / (0.5 * dt), | |
| "forcing_dt": 0.5 * dt, |
Steps of Reproduction ✅
1. Inspect dt definition at line 27 (`dt = 1.0e-06`) in
`examples/3D_random_sphere_array/case.py`.
2. Run script: `python examples/3D_random_sphere_array/case.py` to print case JSON (print
at line 164).
3. In produced JSON, observe `case_dict["forcing_dt"]` printed from line 150 with value
`1.0 / (0.5 * dt)` → for dt=1e-6 this equals 2_000_000.0, not the likely intended 5e-07.
4. Confirm mismatch between intended "half timestep" and the reciprocal value by
inspecting line 150 and dt at line 27.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** examples/3D_random_sphere_array/case.py
**Line:** 150:150
**Comment:**
*Logic Error: `forcing_dt` is computed as the reciprocal of `0.5 * dt` (i.e. 2/dt) which yields an extremely large value for small `dt`; it's almost certainly intended to be half the timestep (`0.5 * dt`) rather than its inverse.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.|
|
||
| if (present(beta)) then | ||
| alt_sys = sys_size + 1 | ||
| else if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop)) then |
There was a problem hiding this comment.
Suggestion: Logic mismatch: the code increases alt_sys to include volume_filter_dt%stat_size when q_filtered_wrt is true, but the MPI I/O data structures for those extra variables are only initialized via s_initialize_mpi_data_filtered inside the if (ib) branch earlier; if ib is false this will cause out-of-bounds access when the code later loops to alt_sys. Restrict adding the extra variables to the same condition that initializes them (i.e., require ib), so the number of variables matches the MPI I/O initialization. [logic error]
Severity Level: Critical 🚨
- ❌ MPI restart/restart_data writes can crash.
- ❌ Parallel MPI I/O may write corrupted restart files.
- ⚠️ Restarts using volume-filter stats become unreliable.
- ⚠️ Affects s_write_parallel_data_files restart output.| else if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop)) then | |
| else if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop) .and. ib) then |
Steps of Reproduction ✅
1. Build and run the MPI-enabled simulation where data output uses parallel I/O (call
path: s_write_data_files -> s_write_parallel_data_files in
src/simulation/m_data_output.fpp). The alt_sys computation is at lines 848-854 in
s_write_parallel_data_files.
2. Ensure file_per_process = .false. so the non-file-per-process branch is executed (this
is the branch that uses alt_sys later to size MPI view/writes).
3. Configure the run with q_filtered_wrt = .true. and trigger a write at t_step == 0 (or
t_step == t_step_stop). Also ensure ib = .false. so IBM-specific initialization is
skipped.
4. Execution reaches the block at src/simulation/m_data_output.fpp:848-854 and sets
alt_sys = sys_size + volume_filter_dt%stat_size, but because ib = .false. the code path
that registers extra MPI variables (call s_initialize_mpi_data_filtered) is not executed
(the s_initialize_mpi_data_filtered call is only invoked inside the ib branch earlier in
s_write_parallel_data_files). Later the code loops up to alt_sys and uses
MPI_IO_DATA%view(i)/var(i) entries that were never initialized, causing out-of-bounds
reads/writes or MPI I/O errors (crash or corrupted restart file).Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/simulation/m_data_output.fpp
**Line:** 850:850
**Comment:**
*Logic Error: Logic mismatch: the code increases `alt_sys` to include `volume_filter_dt%stat_size` when `q_filtered_wrt` is true, but the MPI I/O data structures for those extra variables are only initialized via `s_initialize_mpi_data_filtered` inside the `if (ib)` branch earlier; if `ib` is false this will cause out-of-bounds access when the code later loops to `alt_sys`. Restrict adding the extra variables to the same condition that initializes them (i.e., require `ib`), so the number of variables matches the MPI I/O initialization.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.| if (store_levelset) then | ||
| @:ALLOCATE(levelset%sf(-buff_size:m+buff_size, & | ||
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs)) | ||
| @:ALLOCATE(levelset_norm%sf(-buff_size:m+buff_size, & | ||
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs, 1:3)) | ||
| end if |
There was a problem hiding this comment.
Suggestion: Inconsistent conditional allocation: the new conditional allocation of levelset%sf and levelset_norm%sf is only performed when p > 0 and store_levelset is true, but other code paths (e.g. unguarded GPU updates and the unconditional deallocation/zeroing in other routines) expect those arrays to exist. This can lead to unallocated-array accesses or resource leaks when store_levelset is false or when p > 0. Make the allocation behavior consistent with usage — either allocate unconditionally (as the rest of the code expects) or ensure all uses and deallocations are guarded by the same store_levelset condition. [resource leak]
Severity Level: Critical 🚨
- ❌ s_ibm_setup may access unallocated arrays causing crashes.
- ❌ s_update_mib writes levelset arrays when unallocated.
- ⚠️ GPU data macros assume arrays exist unconditionally.
- ⚠️ Simulations with store_levelset=.false. unstable.| if (store_levelset) then | |
| @:ALLOCATE(levelset%sf(-buff_size:m+buff_size, & | |
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs)) | |
| @:ALLOCATE(levelset_norm%sf(-buff_size:m+buff_size, & | |
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs, 1:3)) | |
| end if | |
| @:ALLOCATE(levelset%sf(-buff_size:m+buff_size, & | |
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs)) | |
| @:ALLOCATE(levelset_norm%sf(-buff_size:m+buff_size, & | |
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs, 1:3)) |
Steps of Reproduction ✅
1. Inspect initialization in src/simulation/m_ibm.fpp: the allocation guard exists at
lines 70-75 inside s_initialize_ibm_module (`if (store_levelset) then ... end if`). This
prevents levelset%sf and levelset_norm%sf allocation when store_levelset is .false.
2. Observe s_ibm_setup (called during module initialization) performs device updates which
expect these arrays. The ACC setup lines for SFs are conditionally added at lines 86-88,
but other code paths reference levelset arrays unconditionally (e.g.,
$:GPU_UPDATE(device='[levelset%sf]') appears in the IBM setup sequence after buffers are
populated).
3. Run a configuration with p>0 and store_levelset = .false. Start the simulation so
s_initialize_ibm_module() runs (allocation skipped). When s_ibm_setup proceeds it will
execute GPU update / host routines referencing levelset%sf or levelset_norm%sf that are
not allocated.
4. Expected outcome: Fortran runtime reports "unallocated variable" or the GPU pragmas
produce undefined behavior or crash. If instead store_levelset is true, behavior is
unchanged. This demonstrates a real mismatch between conditional allocation and
unconditional use; the pattern is not merely stylistic.Prompt for AI Agent 🤖
This is a comment left during a code review.
**Path:** src/simulation/m_ibm.fpp
**Line:** 70:75
**Comment:**
*Resource Leak: Inconsistent conditional allocation: the new conditional allocation of `levelset%sf` and `levelset_norm%sf` is only performed when `p > 0` and `store_levelset` is true, but other code paths (e.g. unguarded GPU updates and the unconditional deallocation/zeroing in other routines) expect those arrays to exist. This can lead to unallocated-array accesses or resource leaks when `store_levelset` is false or when `p > 0`. Make the allocation behavior consistent with usage — either allocate unconditionally (as the rest of the code expects) or ensure all uses and deallocations are guarded by the same `store_levelset` condition.
Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.
|
CodeAnt AI Incremental review completed. |
There was a problem hiding this comment.
1 issue found across 16 files (changes from recent commits).
Prompt for AI agents (all issues)
Check if these issues are valid — if so, understand the root cause of each and fix them.
<file name="src/simulation/m_ibm.fpp">
<violation number="1" location="src/simulation/m_ibm.fpp:70">
P1: Conditional allocation of levelset/levelset_norm leaves later unconditional uses on unallocated arrays when store_levelset is false, which will crash or corrupt memory in 3D runs.</violation>
</file>
Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs)) | ||
| @:ALLOCATE(levelset_norm%sf(-buff_size:m+buff_size, & | ||
| -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs, 1:3)) | ||
| if (store_levelset) then |
There was a problem hiding this comment.
P1: Conditional allocation of levelset/levelset_norm leaves later unconditional uses on unallocated arrays when store_levelset is false, which will crash or corrupt memory in 3D runs.
Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_ibm.fpp, line 70:
<comment>Conditional allocation of levelset/levelset_norm leaves later unconditional uses on unallocated arrays when store_levelset is false, which will crash or corrupt memory in 3D runs.</comment>
<file context>
@@ -67,10 +67,12 @@ contains
- -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs))
- @:ALLOCATE(levelset_norm%sf(-buff_size:m+buff_size, &
- -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs, 1:3))
+ if (store_levelset) then
+ @:ALLOCATE(levelset%sf(-buff_size:m+buff_size, &
+ -buff_size:n+buff_size, -buff_size:p+buff_size, 1:num_ibs))
</file context>
There was a problem hiding this comment.
Actionable comments posted: 7
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/simulation/m_ibm.fpp (1)
76-83: 2D case allocates levelset unconditionally, inconsistent with 3D and deallocation.When
p == 0(2D), lines 79-82 allocatelevelset%sfandlevelset_norm%sfunconditionally, but lines 1321-1324 only deallocate them whenstore_levelsetis true. This creates a memory leak in 2D whenstore_levelset = .false..🐛 Suggested fix
else @:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, & -buff_size:n+buff_size, 0:0)) - @:ALLOCATE(levelset%sf(-buff_size:m+buff_size, & - -buff_size:n+buff_size, 0:0, 1:num_ibs)) - @:ALLOCATE(levelset_norm%sf(-buff_size:m+buff_size, & - -buff_size:n+buff_size, 0:0, 1:num_ibs, 1:3)) + if (store_levelset) then + @:ALLOCATE(levelset%sf(-buff_size:m+buff_size, & + -buff_size:n+buff_size, 0:0, 1:num_ibs)) + @:ALLOCATE(levelset_norm%sf(-buff_size:m+buff_size, & + -buff_size:n+buff_size, 0:0, 1:num_ibs, 1:3)) + end if end if
🤖 Fix all issues with AI agents
In `@examples/3D_random_sphere_array/case.py`:
- Around line 36-38: The code assumes sphere_loc is 2-D but np.loadtxt can
return a 1-D array for single-row files; change the loading to ensure a 2-D
array and compute count from rows by wrapping the result with np.atleast_2d
(e.g., sphere_loc = np.atleast_2d(np.loadtxt(...))) and replace len(sphere_loc)
with sphere_loc.shape[0] so references like sphere_loc[i, 0] work safely.
- Around line 164-166: case_dict may contain NumPy scalar values from
computations (e.g., np.sqrt or array indexing) which json.dumps cannot
serialize; update the serialization to convert NumPy types to native Python
types by passing a default handler or by casting values before dumping. For
example, after case_dict.update(ib_dict) call json.dumps(case_dict,
default=lambda o: o.item() if hasattr(o, "item") else o) or pre-walk
case_dict/ib_dict to replace NumPy scalars with float/int via .item(); reference
the variables case_dict, ib_dict and the json.dumps call to locate where to
apply this change.
In `@src/common/m_ib_patches.fpp`:
- Around line 87-89: The initialization currently calls levelset routines
unconditionally for several geometry types and must be protected by the same
store_levelset guard used for spheres: wrap each call that computes levelset and
levelset_norm (e.g., calls like s_cuboid_levelset, s_cylinder_levelset,
s_airfoil3d_levelset and the 2D variants such as s_circle_levelset,
s_rectangle_levelset, s_airfoil2d_levelset, s_ellipse_levelset) in if
(store_levelset) then ... end if, so levelset/levelset_norm are only referenced
when they have been allocated; apply this change to every geometry branch in
m_ib_patches.fpp where levelset routines are invoked.
In `@src/post_process/m_data_input.f90`:
- Around line 535-537: Several statements in src/post_process/m_data_input.f90
(notably the conditional invoking s_read_parallel_filtered_data with
q_filtered_wrt, t_step and t_step_stop) exceed 100 characters; wrap those long
Fortran source lines using the continuation character (&) so they satisfy
Fortitude S001. Locate the long conditional and the long call to
s_read_parallel_filtered_data (and the other long lines around the ranges
mentioned: the block from ~671–703 and the line near ~762) and split argument
lists and logical expressions onto subsequent lines with trailing & on the
broken line and leading &/aligned indentation on the continuation lines,
preserving the original order of tokens and parentheses. Ensure no token is
split across lines and keep code readable by grouping related arguments (e.g.,
t_step checks on one line, the call arguments on subsequent lines) while keeping
each resulting line under 100 chars.
In `@src/pre_process/m_initial_condition.fpp`:
- Around line 102-107: The code allocates levelset and levelset_norm only when
store_levelset is true but later calls s_apply_ib_patches and other geometry
routines that unconditionally pass levelset/levelset_norm, causing potential
dereference of unallocated arrays; update callers (e.g., the s_apply_ib_patches
invocation that currently passes ib_markers%sf, levelset, levelset_norm) to
conditionally pass those optional arguments or use the alternate call when
store_levelset is false, and also wrap any computations that write/read
levelset/levelset_norm in an if (store_levelset) guard so non-spherical IBs
never touch unallocated levelset or levelset_norm storage.
In `@src/simulation/m_additional_forcing.fpp`:
- Around line 111-115: The kinetic energy term currently hardcodes three
momentum components (momxb, momxb+1, momxb+2); change it to compute ke by
looping over offsets 0..(num_dims-1) and summing
q_cons_vf(momxb+offset)%sf(i,j,k)**2, then divide by rho and use 0.5*ke in the
spatial_eps update (keep E_idx, dVol and rho usage). Declare and initialize a
local ke variable inside the GPU loop, add it to the loop's private clause, and
replace the explicit three-component expression with the loop-based ke
calculation so the code works for 1D/2D/3D.
- Around line 48-52: The calculation for domain_vol unconditionally uses
domain_glb(3,*) which causes out-of-bounds in 2D; update the block around
domain_vol to compute volume based on num_dims (use 3D product when num_dims==3,
2D product when num_dims==2, and 1D when num_dims==1) before computing
avg_coeff, ensuring domain_glb indices are only accessed for existing dimensions
and preserving the subsequent avg_coeff = 1._wp/(domain_vol*(1._wp -
volfrac_phi)).
♻️ Duplicate comments (15)
src/pre_process/m_data_output.fpp (1)
842-865: Fix MPI datatype/count for levelset and levelset_norm writes.
Usingmpi_pand raw element counts bypasses mixed-precision sizing; usempi_io_pand scale bympi_io_type.🐛 Proposed fix
- call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, MPI_IO_levelset_DATA%view, & + call MPI_FILE_SET_VIEW(ifile, disp, mpi_io_p, MPI_IO_levelset_DATA%view, & 'native', mpi_info_int, ierr) - call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelset_DATA%var%sf, data_size*num_ibs, & - mpi_p, status, ierr) + call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelset_DATA%var%sf, data_size*num_ibs*mpi_io_type, & + mpi_io_p, status, ierr) … - call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, MPI_IO_levelsetnorm_DATA%view, & + call MPI_FILE_SET_VIEW(ifile, disp, mpi_io_p, MPI_IO_levelsetnorm_DATA%view, & 'native', mpi_info_int, ierr) - call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelsetnorm_DATA%var%sf, data_size*num_ibs*3, & - mpi_p, status, ierr) + call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelsetnorm_DATA%var%sf, data_size*num_ibs*3*mpi_io_type, & + mpi_io_p, status, ierr)src/post_process/m_global_parameters.fpp (1)
1066-1070: Deallocate q_filtered_wrt MPI_IO buffers before nulling.
The extra stats slots are allocated but only nullified here, which leaks memory. Based on learnings, please use the device-aware deallocation macro.🐛 Proposed fix
if (q_filtered_wrt) then do i = sys_size + 1, sys_size + volume_filter_dt%stat_size - MPI_IO_DATA%var(i)%sf => null() + if (associated(MPI_IO_DATA%var(i)%sf)) then + @:DEALLOCATE(MPI_IO_DATA%var(i)%sf) + end if + MPI_IO_DATA%var(i)%sf => null() end do end iftoolchain/mfc/case_validator.py (2)
1111-1124: Remove trailing whitespace in the new validators (lint failure).
Pylint is failing on trailing whitespace in these added blocks; trimming resolves the pipeline error.Also applies to: 1162-1200
1194-1195: Fix the even-dimension parity check.
The current(dim + 1) % 2test rejects valid even sizes; check the dimension directly.🐛 Proposed fix
- self.prohibit((m_glb + 1) % 2 != 0 or (n_glb + 1) % 2 != 0 or (p_glb + 1) % 2 != 0, + self.prohibit(m_glb % 2 != 0 or n_glb % 2 != 0 or p_glb % 2 != 0, "volume_filter_momentum_eqn requires global dimensions divisible by 2 in every direction")src/common/m_mpi_common.fpp (1)
386-391: Avoid overwriting base MPI views when initializing filtered views.
The loop starts atsys_size, which can reinitialize the last base variable’s view. Start atsys_size + 1instead.🐛 Proposed fix
- do i = sys_size, alt_sys + do i = sys_size + 1, alt_sys call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, & MPI_ORDER_FORTRAN, mpi_p, MPI_IO_DATA%view(i), ierr) call MPI_TYPE_COMMIT(MPI_IO_DATA%view(i), ierr) end dosrc/common/m_ib_patches.fpp (1)
619-665: Avoid shared-state coordinate conversion inside the GPU loop.
s_convert_cylindrical_to_cartesian_coordwrites module-levelcart_y/cart_z, which is unsafe in parallel loops. Use the pure function with thread-local temporaries instead.🐛 Suggested fix (thread‑local conversion)
- if (grid_geometry == 3) then - call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) - else - cart_y = y_cc(j) - cart_z = z_cc(k) - end if + if (grid_geometry == 3) then + tmp_cyl = (/0._wp, y_cc(j), z_cc(k)/) + tmp_cart = f_convert_cyl_to_cart(tmp_cyl) + cart_y = tmp_cart(2) + cart_z = tmp_cart(3) + else + cart_y = y_cc(j) + cart_z = z_cc(k) + end ifsrc/simulation/m_data_output.fpp (2)
85-123: Serial I/O path does not forward filtered data arguments.When
parallel_iois false (Line 117-118), the call tos_write_serial_data_filesdoes not pass the new filtered arguments (filtered_fluid_indicator_function,stat_q_cons_filtered, etc.). Ifq_filtered_wrtis enabled with serial I/O, the filtered data will not be written.Either extend
s_write_serial_data_filesto accept and handle these arguments, or add a runtime check/abort whenq_filtered_wrtis true andparallel_iois false.
801-814: Argument count exceeds the ≤6 guideline.
s_write_parallel_data_filesnow accepts 10 arguments. As per coding guidelines, consider introducing a derived-type params struct to bundle the filtered/stat fields and reduce the signature complexity.src/simulation/m_ibm.fpp (2)
625-634: Off-by-one risk: loop accessess_cc(index + 1)but bound allowsindex == bound.The bound calculation at lines 627, 630, 633 sets
bound = m/n/p + buff_size, but the while-loop at line 648-656 accessess_cc(index + 1). Whenindex == bound,s_cc(index + 1)reads past the allocated array bounds.🐛 Suggested fix
if (dim == 1) then s_cc => x_cc - bound = m + buff_size + bound = m + buff_size - 1 elseif (dim == 2) then s_cc => y_cc - bound = n + buff_size + bound = n + buff_size - 1 else s_cc => z_cc - bound = p + buff_size + bound = p + buff_size - 1 end if
651-655: Replaceerror stopwiths_mpi_abortfor MPI safety.Using
error stopcan leave other MPI ranks hanging. As per coding guidelines, uses_mpi_abortfor error termination.🐛 Suggested fix
if (index < -buff_size .or. index > bound) then - print *, "proc_rank=", proc_rank, "temp_loc=", temp_loc, " index=", index, "ib=", patch_id, "dim", dim, "dir", dir, "i, j, k", i, j, k - print *, "Increase buff_size further in m_helper_basic (currently set to a minimum of 10)" - error stop "Increase buff_size" + call s_mpi_abort("Image point search exceeded buffer bounds. Increase buff_size (minimum 10).") end ifsrc/simulation/m_additional_forcing.fpp (1)
170-172: Division by zero risk when local density is zero.Line 172 divides by
rho, which is computed by summing fluid densities. At phase interfaces or in regions where all fluids are absent,rhocould be zero or near-zero, causing division by zero or numerical instability.🐛 Suggested fix
! f_E - q_periodic_force(3)%sf(i, j, k) = (P_inf_ref*gammas(1) - phase_eps)*forcing_dt & - + q_cons_vf(contxe + mom_f_idx)%sf(i, j, k)*q_periodic_force(2)%sf(i, j, k)/rho + if (rho > 1.e-12_wp) then + q_periodic_force(3)%sf(i, j, k) = (P_inf_ref*gammas(1) - phase_eps)*forcing_dt & + + q_cons_vf(contxe + mom_f_idx)%sf(i, j, k)*q_periodic_force(2)%sf(i, j, k)/rho + else + q_periodic_force(3)%sf(i, j, k) = (P_inf_ref*gammas(1) - phase_eps)*forcing_dt + end ifsrc/simulation/m_volume_filtering.fpp (4)
214-216: Misleading error message: references 'p' but checks 'Ny'.The check validates
mod(Ny, num_procs)but the error message incorrectly states the requirement is forp.✏️ Suggested fix
if (mod(Ny, num_procs) /= 0) then - call s_mpi_abort('Volume filtering requires p to be divisible by the number of ranks') + call s_mpi_abort('Volume filtering requires Ny to be divisible by the number of ranks') end if
585-598: Critical: Out-of-bounds array access in finite-difference stencil.The loop iterates
i = 0, m,j = 0, n,k = 0, p, but the stencil accessesi + lwherel ∈ [-4, 4]. At boundary cells (e.g.,i = 0), this accessesfluid_indicator_function%sf(-4, j, k), causing out-of-bounds memory access unlessidwbuffprovides sufficient ghost cells.Verify that
idwbuffbounds provide at least 4 ghost cells in each direction, or restrict the loop ranges:🐛 Option 1: Restrict loop ranges (safer)
$:GPU_PARALLEL_LOOP(collapse=3, copyin='[fd_coeffs]') - do i = 0, m - do j = 0, n - do k = 0, p + do i = 4, m - 4 + do j = 4, n - 4 + do k = 4, p - 4 $:GPU_LOOP(parallelism='[seq]') do l = -4, 4🔍 Option 2: Verify buffer allocation
#!/bin/bash # Check idwbuff allocation bounds to verify sufficient ghost cells rg -n 'idwbuff' --type-add 'fpp:*.fpp' --type=fpp -C3 | head -60
785-787: Potential division by zero on filtered density.When
q_cons_filtered(1)%sf(i, j, k)(filtered density) approaches zero inside solid regions, this division produces NaN/Inf values.🐛 Proposed fix with epsilon floor
+ real(wp), parameter :: eps_rho = 1.0e-12_wp reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) & - - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k)/q_cons_filtered(1)%sf(i, j, k)) + - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k) & + / max(q_cons_filtered(1)%sf(i, j, k), eps_rho))
1552-1560: Multiple divisions byfiltered_fluid_indicator_functionrisk NaN in solid regions.Lines 1556, 1605, 1655, and 1707 all divide by
filtered_fluid_indicator_function%sf(i, j, k), which is zero inside particles. This systematically produces NaN/Inf values throughout the filtered fields.Apply a consistent epsilon guard across all four locations:
🐛 Proposed pattern for all divisions
Define a module-level constant and apply uniformly:
real(wp), parameter :: eps_indicator = 1.0e-12_wpThen replace each division:
-q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/filtered_fluid_indicator_function%sf(i, j, k) +q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm & + / max(filtered_fluid_indicator_function%sf(i, j, k), eps_indicator)Apply the same pattern at lines 1605, 1655, and 1707.
🧹 Nitpick comments (8)
src/pre_process/m_global_parameters.fpp (1)
949-954: Use @:ALLOCATE for the new levelset MPI buffers.
These allocations are in a.fppmodule and should use the device-aware macro. As per coding guidelines, please prefer@:ALLOCATEhere.♻️ Suggested change
- allocate (MPI_IO_levelset_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs)) - allocate (MPI_IO_levelsetnorm_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)) + @:ALLOCATE(MPI_IO_levelset_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs), & + MPI_IO_levelsetnorm_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3))src/post_process/m_global_parameters.fpp (1)
872-878: Use @:ALLOCATE for q_filtered_wrt MPI_IO_DATA buffers.
These new allocations are in a.fppmodule; using the device-aware macro keeps GPU builds consistent. As per coding guidelines, please prefer@:ALLOCATEhere.♻️ Suggested change
- allocate (MPI_IO_DATA%var(i)%sf(0:m, 0:n, 0:p)) + @:ALLOCATE(MPI_IO_DATA%var(i)%sf(0:m, 0:n, 0:p))src/post_process/m_data_input.f90 (1)
671-680: Consider bundling MPI I/O params into a derived type.
This signature has 8 parameters; grouping the MPI offsets into a small derived type would align with the ≤6-argument guideline and simplify call sites. As per coding guidelines, ...src/simulation/m_ibm.fpp (2)
713-723: Remove or enable commented-out code.Lines 714-723 contain commented-out alternative logic for 1D subsection checks. If this code is no longer needed, remove it to reduce clutter. If it's intended for future use, consider adding a TODO comment explaining its purpose.
805-811: Remove or enable commented-out code.Lines 806-811 contain the same commented-out 1D subsection logic pattern. This duplicates the commented code in
s_find_num_ghost_points. Either remove both or document the intent.src/simulation/m_volume_filtering.fpp (3)
339-343: Consider usingnewunit=for automatic unit assignment.Hardcoded unit number
100may conflict with other file operations in the codebase. Usingnewunit=provides safer automatic assignment.♻️ Suggested improvement
+ integer, save :: particle_force_unit ! file for particle forces if (compute_particle_drag) then if (proc_rank == 0) then - open (unit=100, file='particle_force.bin', status='replace', form='unformatted', access='stream', action='write') + open (newunit=particle_force_unit, file='particle_force.bin', status='replace', form='unformatted', access='stream', action='write') end if end ifThen update the
writeandclosestatements to useparticle_force_unitinstead of100.
1037-1765: Subroutines_filter_batchexceeds size guidelines (~730 lines).Per coding guidelines, subroutines should be ≤ 500 lines. This subroutine contains highly repetitive FFT transform patterns for different field types (conservative vars, pressure, Reynolds stress, effective viscosity, interphase momentum exchange).
Consider extracting reusable helper routines:
s_fft_forward_xy: Handles X FFT → reorder → Y FFT → store to batch slabs_fft_z_convolve: Handles Z FFT → convolution → inverse Z FFTs_fft_backward_yx: Handles Y inverse FFT → reorder → X inverse FFT → store resultThis would reduce
s_filter_batchto ~100-150 lines of orchestration logic and improve maintainability.
3-2015: File exceeds 1000-line guideline (~2015 lines).The coding guidelines recommend files ≤ 1000 lines. Consider splitting into focused submodules:
m_volume_filtering: Public API, orchestration, core filtering workflowm_volume_filtering_fft: FFT plans, forward/backward transforms, slab transposesm_volume_filtering_terms: Stress tensor, Reynolds stress, effective viscosity computations
| # load initial sphere locations | ||
| sphere_loc = np.loadtxt("sphere_array_locations.txt") | ||
| N_sphere = len(sphere_loc) |
There was a problem hiding this comment.
Handle single-row sphere files safely.
np.loadtxt returns a 1‑D array when the file has one row, so len(sphere_loc) becomes 3 and sphere_loc[i, 0] will fail. Use np.atleast_2d and shape[0] to keep consistent 2‑D shape.
✅ Proposed fix
-# load initial sphere locations
-sphere_loc = np.loadtxt("sphere_array_locations.txt")
-N_sphere = len(sphere_loc)
+# load initial sphere locations
+sphere_loc = np.atleast_2d(np.loadtxt("sphere_array_locations.txt"))
+N_sphere = sphere_loc.shape[0]📝 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.
| # load initial sphere locations | |
| sphere_loc = np.loadtxt("sphere_array_locations.txt") | |
| N_sphere = len(sphere_loc) | |
| # load initial sphere locations | |
| sphere_loc = np.atleast_2d(np.loadtxt("sphere_array_locations.txt")) | |
| N_sphere = sphere_loc.shape[0] |
🤖 Prompt for AI Agents
In `@examples/3D_random_sphere_array/case.py` around lines 36 - 38, The code
assumes sphere_loc is 2-D but np.loadtxt can return a 1-D array for single-row
files; change the loading to ensure a 2-D array and compute count from rows by
wrapping the result with np.atleast_2d (e.g., sphere_loc =
np.atleast_2d(np.loadtxt(...))) and replace len(sphere_loc) with
sphere_loc.shape[0] so references like sphere_loc[i, 0] work safely.
| case_dict.update(ib_dict) | ||
|
|
||
| print(json.dumps(case_dict)) |
There was a problem hiding this comment.
Prevent JSON serialization errors from NumPy scalars.
json.dumps can’t serialize NumPy scalar types (e.g., np.float64) produced by np.sqrt and array indexing. Add a default handler (or cast to float) to avoid a runtime TypeError.
✅ Proposed fix
-case_dict.update(ib_dict)
-
-print(json.dumps(case_dict))
+case_dict.update(ib_dict)
+
+def _json_default(obj):
+ if isinstance(obj, np.generic):
+ return obj.item()
+ raise TypeError(f"Object of type {type(obj).__name__} is not JSON serializable")
+
+print(json.dumps(case_dict, default=_json_default))📝 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.
| case_dict.update(ib_dict) | |
| print(json.dumps(case_dict)) | |
| case_dict.update(ib_dict) | |
| def _json_default(obj): | |
| if isinstance(obj, np.generic): | |
| return obj.item() | |
| raise TypeError(f"Object of type {type(obj).__name__} is not JSON serializable") | |
| print(json.dumps(case_dict, default=_json_default)) |
🤖 Prompt for AI Agents
In `@examples/3D_random_sphere_array/case.py` around lines 164 - 166, case_dict
may contain NumPy scalar values from computations (e.g., np.sqrt or array
indexing) which json.dumps cannot serialize; update the serialization to convert
NumPy types to native Python types by passing a default handler or by casting
values before dumping. For example, after case_dict.update(ib_dict) call
json.dumps(case_dict, default=lambda o: o.item() if hasattr(o, "item") else o)
or pre-walk case_dict/ib_dict to replace NumPy scalars with float/int via
.item(); reference the variables case_dict, ib_dict and the json.dumps call to
locate where to apply this change.
| if (store_levelset) then | ||
| call s_sphere_levelset(i, levelset, levelset_norm) | ||
| end if |
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
# First, get the full context of the file to understand structure
wc -l src/common/m_ib_patches.fppRepository: MFlowCode/MFC
Length of output: 89
🏁 Script executed:
# Look at the relevant section around lines 87-89 to understand the geometry handling
sed -n '70,110p' src/common/m_ib_patches.fppRepository: MFlowCode/MFC
Length of output: 1686
🏁 Script executed:
# Find all store_levelset usages in the codebase
rg -n "store_levelset" src -g '!**/test/**' -A 2 -B 2Repository: MFlowCode/MFC
Length of output: 12569
🏁 Script executed:
# Search for other geometry types (cylinder, cuboid, etc.) and their levelset handling
rg -n "s_.*_levelset|levelset.*cylinder|levelset.*cuboid" src -g '!**/test/**'Repository: MFlowCode/MFC
Length of output: 3102
🏁 Script executed:
# Look for levelset allocation and initialization
rg -n "allocate.*levelset|levelset.*allocat" src -g '!**/test/**' -A 1 -B 1Repository: MFlowCode/MFC
Length of output: 1084
🏁 Script executed:
# Check the 2D section to see if similar guards are missing
sed -n '100,130p' src/common/m_ib_patches.fppRepository: MFlowCode/MFC
Length of output: 1356
🏁 Script executed:
# Verify the validation is actually in the code path that matters
sed -n '556,575p' src/simulation/m_ibm.fppRepository: MFlowCode/MFC
Length of output: 1245
🏁 Script executed:
# Check when m_ib_patches is called in relation to m_ibm initialization
rg -n "s_ib_patches|s_initialize_levelset_field" src -g '!**/test/**' -B 2 -A 2Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Search for where the validation check in m_ibm is called (is it in initialization or runtime?)
rg -n "image_point_locations|gp%ib_patch_id" src -g '!**/test/**' | head -20Repository: MFlowCode/MFC
Length of output: 330
🏁 Script executed:
# Check if levelset is passed by optional argument and if it can be unallocated
sed -n '30,150p' src/common/m_ib_patches.fppRepository: MFlowCode/MFC
Length of output: 5497
🏁 Script executed:
# Check if the levelset subroutines handle unallocated optional arguments gracefully
sed -n '31,80p' src/common/m_compute_levelset.fppRepository: MFlowCode/MFC
Length of output: 1893
🏁 Script executed:
# Check if the levelset subroutines check for present() before accessing optional arguments
rg -n "if.*present.*levelset" src/common/m_compute_levelset.fppRepository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
# Get the full signature and start of cuboid_levelset to see if it checks present
sed -n '396,420p' src/common/m_compute_levelset.fppRepository: MFlowCode/MFC
Length of output: 1055
🏁 Script executed:
# Find where s_apply_ib_patches is called to see if there's an outer guard
rg -n "s_apply_ib_patches" src -g '!**/test/**' -B 5 -A 5Repository: MFlowCode/MFC
Length of output: 3945
🏁 Script executed:
# Check the signature of levelset subroutines to see if they safely handle unallocated optional
sed -n '505,530p' src/common/m_compute_levelset.fppRepository: MFlowCode/MFC
Length of output: 1237
🏁 Script executed:
# Check when s_apply_ib_patches is called - does it only happen when store_levelset=true?
sed -n '200,210p' src/pre_process/m_initial_condition.fppRepository: MFlowCode/MFC
Length of output: 412
🏁 Script executed:
# Check if there's any validation BEFORE calling s_apply_ib_patches
sed -n '1150,1170p' src/simulation/m_ibm.fppRepository: MFlowCode/MFC
Length of output: 882
🏁 Script executed:
# Final verification: Check if any levelset subroutine checks present() before accessing
sed -n '505,570p' src/common/m_compute_levelset.fpp | head -60Repository: MFlowCode/MFC
Length of output: 2887
🏁 Script executed:
# Check cylinder_levelset for the same pattern
sed -n '611,650p' src/common/m_compute_levelset.fppRepository: MFlowCode/MFC
Length of output: 1896
🏁 Script executed:
# Verify: are levelset subroutines called during initialization (s_apply_ib_patches)?
# And does initialization pass allocated or unallocated optional arguments?
sed -n '100,110p' src/pre_process/m_initial_condition.fppRepository: MFlowCode/MFC
Length of output: 418
🏁 Script executed:
# Critical: Check if levelset is allocated before s_apply_ib_patches call during init
sed -n '100,210p' src/pre_process/m_initial_condition.fppRepository: MFlowCode/MFC
Length of output: 4813
Add store_levelset guards to all geometry types, not just spheres.
Currently, only the spherical geometry (lines 87–89) is guarded by if (store_levelset). The cuboid, cylinder, 3D airfoil, and all 2D geometries (circle, rectangle, airfoil, ellipse) call their respective levelset subroutines unconditionally. Since levelset and levelset_norm are only allocated when store_levelset=.true. (see m_initial_condition.fpp:104–107), these unguarded calls will attempt to access unallocated optional arrays during initialization, causing a crash if store_levelset=false with non-spherical geometries.
While a runtime check in m_ibm.fpp:563 prevents non-spherical geometries when store_levelset=false, this validation occurs too late—after initialization has already failed. Apply the guard to all geometry branches: cuboid, cylinder, 3D airfoil, and all 2D types.
🤖 Prompt for AI Agents
In `@src/common/m_ib_patches.fpp` around lines 87 - 89, The initialization
currently calls levelset routines unconditionally for several geometry types and
must be protected by the same store_levelset guard used for spheres: wrap each
call that computes levelset and levelset_norm (e.g., calls like
s_cuboid_levelset, s_cylinder_levelset, s_airfoil3d_levelset and the 2D variants
such as s_circle_levelset, s_rectangle_levelset, s_airfoil2d_levelset,
s_ellipse_levelset) in if (store_levelset) then ... end if, so
levelset/levelset_norm are only referenced when they have been allocated; apply
this change to every geometry branch in m_ib_patches.fpp where levelset routines
are invoked.
| if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop)) then | ||
| call s_read_parallel_filtered_data(t_step, m_MOK, n_MOK, p_MOK, WP_MOK, MOK, str_MOK, NVARS_MOK) | ||
| end if |
There was a problem hiding this comment.
Wrap long lines to satisfy Fortitude S001.
Several lines exceed the 100‑char limit. Please wrap with continuation &.
✂️ Example wrap
- call s_read_parallel_filtered_data(t_step, m_MOK, n_MOK, p_MOK, WP_MOK, MOK, str_MOK, NVARS_MOK)
+ call s_read_parallel_filtered_data(t_step, m_MOK, n_MOK, p_MOK, &
+ WP_MOK, MOK, str_MOK, NVARS_MOK)Also applies to: 671-703, 762-762
🧰 Tools
🪛 Fortitude (0.7.5)
[error] 536-536: line length of 108, exceeds maximum 100
(S001)
🤖 Prompt for AI Agents
In `@src/post_process/m_data_input.f90` around lines 535 - 537, Several statements
in src/post_process/m_data_input.f90 (notably the conditional invoking
s_read_parallel_filtered_data with q_filtered_wrt, t_step and t_step_stop)
exceed 100 characters; wrap those long Fortran source lines using the
continuation character (&) so they satisfy Fortitude S001. Locate the long
conditional and the long call to s_read_parallel_filtered_data (and the other
long lines around the ranges mentioned: the block from ~671–703 and the line
near ~762) and split argument lists and logical expressions onto subsequent
lines with trailing & on the broken line and leading &/aligned indentation on
the continuation lines, preserving the original order of tokens and parentheses.
Ensure no token is split across lines and keep code readable by grouping related
arguments (e.g., t_step checks on one line, the call arguments on subsequent
lines) while keeping each resulting line under 100 chars.
| if (ib) then | ||
| allocate (ib_markers%sf(0:m, 0:n, 0:p)) | ||
| allocate (levelset%sf(0:m, 0:n, 0:p, 1:num_ibs)) | ||
| allocate (levelset_norm%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)) | ||
| if (store_levelset) then | ||
| allocate (levelset%sf(0:m, 0:n, 0:p, 1:num_ibs)) | ||
| allocate (levelset_norm%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)) | ||
| end if |
There was a problem hiding this comment.
Guard levelset usage when store_levelset is false.
levelset/levelset_norm are now conditionally allocated, but IB patch application still passes them unconditionally and some geometries compute levelsets without a guard. If store_levelset is false with non-spherical IBs, this will dereference unallocated storage. Consider gating those calls or omitting the optional args when storage is disabled.
if (store_levelset) then
call s_apply_ib_patches(ib_markers%sf, levelset, levelset_norm)
else
call s_apply_ib_patches(ib_markers%sf)
end if🤖 Prompt for AI Agents
In `@src/pre_process/m_initial_condition.fpp` around lines 102 - 107, The code
allocates levelset and levelset_norm only when store_levelset is true but later
calls s_apply_ib_patches and other geometry routines that unconditionally pass
levelset/levelset_norm, causing potential dereference of unallocated arrays;
update callers (e.g., the s_apply_ib_patches invocation that currently passes
ib_markers%sf, levelset, levelset_norm) to conditionally pass those optional
arguments or use the alternate call when store_levelset is false, and also wrap
any computations that write/read levelset/levelset_norm in an if
(store_levelset) guard so non-spherical IBs never touch unallocated levelset or
levelset_norm storage.
| ! total cartesian domain volume | ||
| domain_vol = (domain_glb(1, 2) - domain_glb(1, 1))*(domain_glb(2, 2) - domain_glb(2, 1))*(domain_glb(3, 2) - domain_glb(3, 1)) | ||
|
|
||
| ! coefficient used for phase averages | ||
| avg_coeff = 1._wp/(domain_vol*(1._wp - volfrac_phi)) |
There was a problem hiding this comment.
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
# First, check if the file exists and read the full context around the flagged lines
cat -n src/simulation/m_additional_forcing.fpp | head -100Repository: MFlowCode/MFC
Length of output: 4067
🏁 Script executed:
#!/bin/bash
# Search for domain_glb initialization and how it's set up
rg "domain_glb\s*=" src/ -B 2 -A 2Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
#!/bin/bash
# Check for num_dims or ndims references to understand dimensionality handling
rg "num_dims|ndims" src/simulation/m_additional_forcing.fpp -B 2 -A 2Repository: MFlowCode/MFC
Length of output: 39
🏁 Script executed:
#!/bin/bash
# Search for how domain_vol is used to understand if 2D vs 3D is handled elsewhere
rg "domain_vol" src/ -B 2 -A 2 | head -100Repository: MFlowCode/MFC
Length of output: 1056
🏁 Script executed:
#!/bin/bash
# Search for domain_glb definition - likely in m_global_parameters or similar
rg "domain_glb" src/ -B 1 -A 1 | head -80Repository: MFlowCode/MFC
Length of output: 6093
🏁 Script executed:
#!/bin/bash
# Look for dimension-related parameters
rg "num_dims|n_dims|ndims" src/ --include="*.fpp" --include="*.f90" | head -30Repository: MFlowCode/MFC
Length of output: 498
🏁 Script executed:
#!/bin/bash
# Check if m_global_parameters exists and search for domain_glb there
fd "m_global_parameters" src/ -type fRepository: MFlowCode/MFC
Length of output: 227
🏁 Script executed:
#!/bin/bash
# Find where domain_glb is declared
rg "domain_glb" src/ | grep -E "(dimension|allocate|real)" | head -20Repository: MFlowCode/MFC
Length of output: 149
🏁 Script executed:
#!/bin/bash
# Look for global_parameters module file
find src -name "*global*" -type fRepository: MFlowCode/MFC
Length of output: 176
🏁 Script executed:
#!/bin/bash
# Check for dimension-related code without the problematic flag
rg "dimension.*3" src/ | head -30Repository: MFlowCode/MFC
Length of output: 3235
🏁 Script executed:
#!/bin/bash
# Check how domain_glb is allocated in m_mpi_common
cat -n src/common/m_mpi_common.fpp | grep -A 5 -B 5 "domain_glb"Repository: MFlowCode/MFC
Length of output: 2165
🏁 Script executed:
#!/bin/bash
# Look at m_global_parameters for the simulation to see dimension handling
cat -n src/simulation/m_global_parameters.fpp | head -200Repository: MFlowCode/MFC
Length of output: 11355
🏁 Script executed:
#!/bin/bash
# Search for where domain_glb is set/initialized
rg "domain_glb\(" src/ -B 3 -A 3 | grep -E "(allocate|=|dimension)" | head -40Repository: MFlowCode/MFC
Length of output: 3891
Fix domain volume calculation to handle 2D simulations.
Line 49 unconditionally accesses domain_glb(3,1) and domain_glb(3,2), but for 2D simulations domain_glb is allocated with size (num_dims, 2), causing out-of-bounds access. Guard the calculation based on num_dims:
Suggested fix
! total cartesian domain volume
if (num_dims == 3) then
domain_vol = (domain_glb(1, 2) - domain_glb(1, 1))*(domain_glb(2, 2) - domain_glb(2, 1))*(domain_glb(3, 2) - domain_glb(3, 1))
else if (num_dims == 2) then
domain_vol = (domain_glb(1, 2) - domain_glb(1, 1))*(domain_glb(2, 2) - domain_glb(2, 1))
else
domain_vol = (domain_glb(1, 2) - domain_glb(1, 1))
end if🤖 Prompt for AI Agents
In `@src/simulation/m_additional_forcing.fpp` around lines 48 - 52, The
calculation for domain_vol unconditionally uses domain_glb(3,*) which causes
out-of-bounds in 2D; update the block around domain_vol to compute volume based
on num_dims (use 3D product when num_dims==3, 2D product when num_dims==2, and
1D when num_dims==1) before computing avg_coeff, ensuring domain_glb indices are
only accessed for existing dimensions and preserving the subsequent avg_coeff =
1._wp/(domain_vol*(1._wp - volfrac_phi)).
| spatial_eps = spatial_eps + ((q_cons_vf(E_idx)%sf(i, j, k) - 0.5_wp*( & | ||
| q_cons_vf(momxb)%sf(i, j, k)**2 + & | ||
| q_cons_vf(momxb + 1)%sf(i, j, k)**2 + & | ||
| q_cons_vf(momxb + 2)%sf(i, j, k)**2)/ & | ||
| rho)*dVol) ! rho*e |
There was a problem hiding this comment.
Kinetic energy calculation assumes exactly 3 momentum components.
Lines 111-115 explicitly reference momxb, momxb + 1, and momxb + 2, which assumes a 3D simulation. For 1D or 2D cases, this will access invalid array indices.
🐛 Suggested fix using num_dims
- spatial_eps = spatial_eps + ((q_cons_vf(E_idx)%sf(i, j, k) - 0.5_wp*( &
- q_cons_vf(momxb)%sf(i, j, k)**2 + &
- q_cons_vf(momxb + 1)%sf(i, j, k)**2 + &
- q_cons_vf(momxb + 2)%sf(i, j, k)**2)/ &
- rho)*dVol) ! rho*e
+ ke = 0._wp
+ do l = 1, num_dims
+ ke = ke + q_cons_vf(momxb + l - 1)%sf(i, j, k)**2
+ end do
+ spatial_eps = spatial_eps + ((q_cons_vf(E_idx)%sf(i, j, k) - 0.5_wp*ke/rho)*dVol) ! rho*eNote: ke would need to be declared and added to the private clause of the GPU loop.
🤖 Prompt for AI Agents
In `@src/simulation/m_additional_forcing.fpp` around lines 111 - 115, The kinetic
energy term currently hardcodes three momentum components (momxb, momxb+1,
momxb+2); change it to compute ke by looping over offsets 0..(num_dims-1) and
summing q_cons_vf(momxb+offset)%sf(i,j,k)**2, then divide by rho and use 0.5*ke
in the spatial_eps update (keep E_idx, dVol and rho usage). Declare and
initialize a local ke variable inside the GPU loop, add it to the loop's private
clause, and replace the explicit three-component expression with the loop-based
ke calculation so the code works for 1D/2D/3D.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1043 +/- ##
==========================================
- Coverage 44.07% 42.92% -1.15%
==========================================
Files 71 74 +3
Lines 20369 21876 +1507
Branches 1986 2064 +78
==========================================
+ Hits 8978 9391 +413
- Misses 10246 11300 +1054
- Partials 1145 1185 +40 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
3 participants
User description
User description
Description
Please include a summary of the changes and the related issue(s) if they exist.
Please also include relevant motivation and context.
Fixes #(issue) [optional]
-Overarching goal of investigating unclosed terms present in the volume filtered Euler-Lagrange framework
-Implements explicit computation of the unclosed terms present in the volume filtered momentum equation in the Euler-Lagrange framework for particle laden flows
-Explicit immersed boundary drag force computation
-General field statistics computation (1st-4th order statistical moments)
-Periodic spherical immersed boundaries
-Bulk flow periodic forcing
Type of change
Please delete options that are not relevant.
Scope
If you cannot check the above box, please split your PR into multiple PRs that each have a common goal.
How Has This Been Tested?
Please describe the tests that you ran to verify your changes.
Provide instructions so we can reproduce.
Please also list any relevant details for your test configuration
Test Configuration:
-MacOs
-Linux (Ubuntu)
-NCSA Delta Nvidia A100 GPUs
Checklist
docs/)examples/that demonstrate my new feature performing as expected.They run to completion and demonstrate "interesting physics"
./mfc.sh formatbefore committing my codeIf your code changes any code source files (anything in
src/simulation)To make sure the code is performing as expected on GPU devices, I have:
nvtxranges so that they can be identified in profiles./mfc.sh run XXXX --gpu -t simulation --nsys, and have attached the output file (.nsys-rep) and plain text results to this PR./mfc.sh run XXXX --gpu -t simulation --rsys --hip-trace, and have attached the output file and plain text results to this PR.PR Type
Enhancement, Tests
Description
Implements comprehensive volume filtering framework for particle-laden flows with explicit computation of unclosed terms in the Euler-Lagrange framework
Adds FFT-based filtering infrastructure supporting both CPU (FFTW) and GPU (cuFFT) backends with MPI parallelization
Computes key quantities including Reynolds stress tensors, effective viscosity, and interphase momentum exchange terms
Implements periodic spherical immersed boundaries with on-the-fly levelset computation when storage is disabled
Adds periodic bulk flow forcing based on space-time averaged quantities with fluid indicator function weighting
Extends MPI communication infrastructure with scalar field send/receive buffers and global domain bounds computation
Implements statistical moment computation (1st-4th order) using Welford's online algorithm with GPU parallelization
Adds support for writing volume-filtered conservative variables, pressure, Reynolds stress, effective viscosity, and interphase momentum exchange
Includes 5 new test cases demonstrating volume filtering, periodic forcing, and particle drag computation
Increases maximum immersed boundary patches from 10 to 1000 to support larger particle-laden systems
Adds comprehensive case validation for volume filtering and periodic forcing requirements
Includes utility scripts for generating 3D Voronoi-based sphere distributions for uniform particle packing
Diagram Walkthrough
File Walkthrough
26 files
m_volume_filtering.fpp
Volume filtering module for Euler-Lagrange particle-laden flowssrc/simulation/m_volume_filtering.fpp
flows with explicit computation of unclosed terms in the
Euler-Lagrange framework
and GPU (cuFFT) backends with MPI parallelization
viscosity, and interphase momentum exchange terms
management for immersed boundary methods
m_data_output.fpp
Add volume-filtered data output capabilities to I/O modulesrc/simulation/m_data_output.fpp
s_write_data_filesands_write_parallel_data_filessubroutineswith optional parameters for filtered quantities
pressure, Reynolds stress, effective viscosity, and interphase
momentum exchange
and final time steps when
q_filtered_wrtflag is enableds_initialize_mpi_data_filteredfor MPI-based initializationof filtered field data structures
m_time_steppers.fpp
Add volume filtering module dependencies to time stepperssrc/simulation/m_time_steppers.fpp
m_volume_filteringandm_compute_statisticscapabilities into time-stepping routines
m_mpi_common.fpp
MPI scalar field communication and global domain boundssrc/common/m_mpi_common.fpp
memory management
domain_glbarray to store global domain bounds across allprocessors
s_mpi_sendrecv_variables_buffers_scalarfieldsubroutinefor scalar field MPI communication with pack/unpack operations
s_mpi_global_domain_boundssubroutine to compute and broadcastglobal domain boundaries
distribution
m_start_up.fpp
Volume filtering and statistics computation integrationsrc/simulation/m_start_up.fpp
additional forcing
store_levelsetflag isenabled
statistics calculation in main time-stepping loop
filtered data
statistics modules
m_compute_statistics.fpp
Statistical moments computation for volume filtered variablessrc/simulation/m_compute_statistics.fpp
for unclosed terms
interphase momentum exchange, and filtered conservative variables
moments with GPU parallelization
module
m_ibm.fpp
Periodic spherical immersed boundaries and on-the-fly levelsetsrc/simulation/m_ibm.fpp
all three directions
spherical boundaries when
store_levelsetis disabledboundaries
of full 3D subsection
m_data_input.f90
Post-processing support for filtered variables and statisticssrc/post_process/m_data_input.f90
statistical moment arrays
s_allocate_filtered_arrayssubroutine to allocate filteredvariable arrays
s_read_parallel_filtered_datasubroutine to read volume filtereddata from restart files
routine
m_additional_forcing.fpp
Periodic bulk flow forcing for particle-laden flowssrc/simulation/m_additional_forcing.fpp
averaged quantities
indicator function weighting
m_boundary_common.fpp
Scalar field boundary condition populationsrc/common/m_boundary_common.fpp
s_populate_scalarfield_bufferssubroutine for populating scalarfield boundary buffers
s_periodic_scalarfieldpure subroutine for periodicboundary conditions on scalar fields
scalar fields
m_ib_patches.fpp
Periodic sphere initialization with domain wrappingsrc/common/m_ib_patches.fpp
of sphere centroids
mark grid cells accordingly
projections
m_compute_levelset.fpp
Periodic levelset computation for spherical boundariessrc/common/m_compute_levelset.fpp
for spheres
domain boundaries
nearest surface
computation
m_mpi_proxy.fpp
MPI parameter broadcasting for new featuressrc/simulation/m_mpi_proxy.fpp
t_step_start_stats,mom_f_idx,forcing_windowperiodic_ibs,compute_particle_drag,periodic_forcing,forcing_wrt,volume_filtering_momentum_eqn,store_levelset,slab_domain_decomposition,q_filtered_wrtu_inf_ref,rho_inf_ref,P_inf_ref,filter_width,forcing_dt,fluid_volume_fractionm_start_up.fpp
Pre-processing domain bounds initializationsrc/pre_process/m_start_up.fpp
s_mpi_global_domain_boundsto compute global domainbounds during preprocessing
and levelset storage flags
m_global_parameters.fpp
Global parameters for volume filtering and periodic forcingsrc/simulation/m_global_parameters.fpp
and immersed boundary features
s_assign_default_values_to_user_inputsfor volume filtering
m_start_up.fpp
Post-processing output for filtered quantities and statisticssrc/post_process/m_start_up.fpp
q_filtered_wrtis enabled(Reynolds stress, effective viscosity, interfacial momentum exchange,
filtered conserved variables, filtered pressure)
m_data_output.fpp
Conditional levelset storage in restart filessrc/pre_process/m_data_output.fpp
store_levelsetflagwriting
store_levelsetis falsem_global_parameters.fpp
Post-process global parameters for filtering and IB featuressrc/post_process/m_global_parameters.fpp
levelset storage, domain decomposition, and filtered output
q_filtered_wrtis enabledm_mpi_proxy.fpp
MPI broadcast for new immersed boundary parameterssrc/pre_process/m_mpi_proxy.fpp
periodic_ibs,store_levelset,slab_domain_decompositionm_mpi_proxy.fpp
MPI broadcast for post-process filtering parameterssrc/post_process/m_mpi_proxy.fpp
periodic_ibs,store_levelset,slab_domain_decomposition,q_filtered_wrtm_rhs.fpp
Periodic forcing integration into RHS computationsrc/simulation/m_rhs.fpp
m_additional_forcingprofiling ranges
s_compute_periodic_forcingwhenperiodic_forcingis enabledm_global_parameters.fpp
Pre-process global parameters for periodic IB supportsrc/pre_process/m_global_parameters.fpp
domain decomposition
s_assign_default_values_to_user_inputsp_main.fpp
Volume filtering initialization in main simulationsrc/simulation/p_main.fpp
m_volume_filteringforcing
when volume filtering is enabled
m_constants.fpp
Increased maximum patches limit for IB supportsrc/common/m_constants.fpp
num_patches_maxfrom 10 to 1000 to support larger numbers ofimmersed boundaries
case_validator.py
Case validation for volume filtering and periodic forcingtoolchain/mfc/case_validator.py
num_ibsvalidation limit from 10 to 1000periodic BCs
check_periodic_forcingmethod validating required parameters forperiodic forcing
check_volume_filteringmethod validating filter width, domaindecomposition, and grid requirements
validate_commonmethodcase_dicts.py
Case dictionary parameter definitions for new featurestoolchain/mfc/run/case_dicts.py
support
SIMULATION
and periodic forcing
q_filtered_wrtto POST_PROCESS dictionary1 files
delta.mako
GPU memory allocation increase for Delta systemtoolchain/templates/delta.mako
NCSA Delta
6 files
case.py
Test case for volume filtered particle-laden flowruns/phi01/case.py
computation
case.py
Single sphere volume filtering test caseruns/3d_1sphere_filtering/case.py
writing
case.py
Periodic immersed boundary test caseruns/3d_1sphere_periodic/case.py
case.py
Off-center periodic IB test configurationruns/3d_periodic_ibs_test/off-center/case.py
case.py
Centered periodic IB test configurationruns/3d_periodic_ibs_test/centered/case.py
case.py
Drag force computation test caseruns/3d_drag_test/case.py
3 files
gen_voronoi_3D.py
3D Voronoi sphere distribution generatorvoronoi/gen_voronoi_3D.py
distributions
gen_voronoi_2D.py
2D Voronoi circle distribution generatorvoronoi/gen_voronoi_2D.py
distributions
sphere_array_locations.txt
Sphere location data for phi=0.1 test caseruns/phi01/sphere_array_locations.txt
flow test case
distribution
Summary by CodeRabbit
New Features
Improvements
✏️ Tip: You can customize this high-level summary in your review settings.
CodeAnt-AI Description
Add volume-filtered momentum term computation, FFT-based filtering, periodic IBs and particle force output
What Changed
Impact
✅ Volume-filtered momentum terms available in outputs✅ Clearer particle force reports (particle_force.bin)✅ Support for periodic immersed boundaries and periodic forcing💡 Usage Guide
Checking Your Pull Request
Every time you make a pull request, our system automatically looks through it. We check for security issues, mistakes in how you're setting up your infrastructure, and common code problems. We do this to make sure your changes are solid and won't cause any trouble later.
Talking to CodeAnt AI
Got a question or need a hand with something in your pull request? You can easily get in touch with CodeAnt AI right here. Just type the following in a comment on your pull request, and replace "Your question here" with whatever you want to ask:
This lets you have a chat with CodeAnt AI about your pull request, making it easier to understand and improve your code.
Example
Preserve Org Learnings with CodeAnt
You can record team preferences so CodeAnt AI applies them in future reviews. Reply directly to the specific CodeAnt AI suggestion (in the same thread) and replace "Your feedback here" with your input:
This helps CodeAnt AI learn and adapt to your team's coding style and standards.
Example
Retrigger review
Ask CodeAnt AI to review the PR again, by typing:
Check Your Repository Health
To analyze the health of your code repository, visit our dashboard at https://app.codeant.ai. This tool helps you identify potential issues and areas for improvement in your codebase, ensuring your repository maintains high standards of code health.