Skip to content

Levelset refactor#1123

Open
danieljvickers wants to merge 49 commits intoMFlowCode:masterfrom
danieljvickers:levelset-refactor
Open

Levelset refactor#1123
danieljvickers wants to merge 49 commits intoMFlowCode:masterfrom
danieljvickers:levelset-refactor

Conversation

@danieljvickers
Copy link
Member

@danieljvickers danieljvickers commented Feb 4, 2026

User description

User description

Description

This PR is a significant refactor of the levelset code. It contains

  1. The removal of IB markers and levelset values from pre-processing. Now all values are computed in simulation
  2. The deletion of the levelset and levelset norm global arrays. The are now tied to ghost points. The levelset code is now parallelized over ghost points instead of grid cells. This reduces the total amount of memory that needs to be allocated for levelsets by an order of magnitude in simple cases, and by several orders of magnitude when using multiple immersed boundaries.
  3. GPU compute of STL IB markers and levelst values. This is done by loading STL models into global memory, and they are computed separately for IB markers and levelset values. Then each separate loop was parallelized.

Fixes #1011

Type of change

  • New feature (non-breaking change which adds functionality)
  • Refactor

Scope

  • This PR comprises a set of related changes with a common goal

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?

This currently passes all tests on MacOS and Ubuntu operating systms using the GNU compiler.

Checklist

  • I have added comments for the new code
  • I added Doxygen docstrings to the new code
  • I have made corresponding changes to the documentation (docs/)
  • I have added regression tests to the test suite so that people can verify in the future that the feature is behaving as expected
  • I have added example cases in examples/ that demonstrate my new feature performing as expected.
    They run to completion and demonstrate "interesting physics"
  • I ran ./mfc.sh format before committing my code
  • New and existing tests pass locally with my changes, including with GPU capability enabled (both NVIDIA hardware with NVHPC compilers and AMD hardware with CRAY compilers) and disabled
  • This PR does not introduce any repeated code (it follows the DRY principle)
  • I cannot think of a way to condense this code and reduce any introduced additional line count

If 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:

  • Checked that the code compiles using NVHPC compilers
  • Checked that the code compiles using CRAY compilers
  • Ran the code on either V100, A100, or H100 GPUs and ensured the new feature performed as expected (the GPU results match the CPU results)
  • Ran the code on MI200+ GPUs and ensure the new features performed as expected (the GPU results match the CPU results)
  • Enclosed the new feature via nvtx ranges so that they can be identified in profiles
  • Ran a Nsight Systems profile using ./mfc.sh run XXXX --gpu -t simulation --nsys, and have attached the output file (.nsys-rep) and plain text results to this PR
  • Ran a Rocprof Systems profile using ./mfc.sh run XXXX --gpu -t simulation --rsys --hip-trace, and have attached the output file and plain text results to this PR.
  • Ran my code using various numbers of different GPUs (1, 2, and 8, for example) in parallel and made sure that the results scale similarly to what happens if you run without the new code/feature

PR Type

Enhancement, Refactor


Description

  • Removed levelset computation from preprocessing: IB markers and levelset values are no longer computed during pre-processing; all computations now occur during simulation runtime

  • Eliminated global levelset arrays: Replaced global levelset and levelset_norm arrays with local storage in the ghost_point derived type, reducing memory allocation by orders of magnitude

  • Parallelized levelset computation over ghost points: Refactored levelset code to parallelize over ghost points instead of grid cells, improving memory efficiency and GPU performance

  • GPU-accelerated STL model processing: Implemented GPU compute for STL IB markers and levelset values by loading STL models into global memory with separate parallelized loops

  • New m_compute_levelset module: Created dedicated module for all levelset computations with geometry-specific functions (circles, rectangles, ellipses, spheres, cylinders, airfoils, STL models)

  • Added STL model configuration parameters: Extended IB patch parameters to support model_filepath, model_spc, model_threshold, and transformation parameters (model_translate, model_rotate)

  • Simplified data I/O: Removed levelset and levelset norm from preprocessing output; added dedicated IB data output subroutines for simulation runtime

  • Updated test infrastructure: Adjusted grid resolution for circle test cases and updated golden test files to reflect removal of pre-computed IB markers


Diagram Walkthrough

flowchart LR
  A["Preprocessing<br/>Input Data"] -->|"No levelset<br/>computation"| B["Simulation<br/>Initialization"]
  B -->|"s_instantiate_STL_models"| C["Load STL Models<br/>to Global Memory"]
  C -->|"GPU compute"| D["Ghost Point<br/>Levelset Storage"]
  D -->|"s_apply_levelset"| E["IBM Updates<br/>During Simulation"]
  E -->|"s_write_ib_data_file"| F["Output IB Data<br/>at Runtime"]
Loading

File Walkthrough

Relevant files
Refactor
8 files
m_ib_patches.fpp
Refactor levelset computation from preprocessing to simulation runtime

src/simulation/m_ib_patches.fpp

  • Removed m_compute_levelset module dependency and levelset computation
    calls from s_apply_ib_patches
  • Added new s_instantiate_STL_models() subroutine to load and preprocess
    STL models during initialization
  • Refactored s_ib_model() to use pre-instantiated models from global
    models array instead of computing levelsets inline
  • Removed local eta variable declarations and replaced with local scope
    usage in geometry-specific functions
  • Added offset variable to airfoil subroutines for cleaner centroid
    offset handling
+165/-234
m_data_output.fpp
Remove levelset data output from preprocessing                     

src/pre_process/m_data_output.fpp

  • Removed ib_markers, levelset, and levelset_norm parameters from data
    output subroutine signatures
  • Removed file I/O operations for writing IB markers, levelset, and
    levelset norm data
  • Simplified abstract interface and implementation for
    s_write_abstract_data_files()
+6/-252 
m_ibm.fpp
Refactor IBM module to use ghost-point-based levelsets     

src/simulation/m_ibm.fpp

  • Removed global levelset and levelset_norm arrays and replaced with
    ghost-point-local storage
  • Added call to s_instantiate_STL_models() during IBM setup
  • Replaced s_apply_ib_patches() calls to remove levelset parameters
  • Changed s_compute_image_points() to use levelset values from
    ghost_point structure instead of global arrays
  • Updated s_update_mib() to call new s_apply_levelset() instead of
    computing levelsets inline
  • Removed patch_id_fp array allocation
+38/-44 
m_start_up.fpp
Remove levelset restart file I/O from simulation startup 

src/simulation/m_start_up.fpp

  • Removed code that reads IB markers, levelset, and levelset norm data
    from restart files
  • Removed airfoil grid data reading from preprocessing
  • Added call to s_write_ib_data_file() after IBM setup during
    initialization
  • Simplified s_read_ic_data_files() call signature by removing IB marker
    parameter
+6/-212 
m_start_up.fpp
Remove levelset handling from preprocessing startup           

src/pre_process/m_start_up.fpp

  • Removed m_ib_patches module import
  • Removed ib_markers, levelset, and levelset_norm variable declarations
    and allocations
  • Removed IB marker reading from preprocessing startup
  • Simplified s_read_ic_data_files() interface by removing IB marker
    parameter
+6/-70   
m_mpi_common.fpp
Simplify MPI data initialization for IBM                                 

src/common/m_mpi_common.fpp

  • Removed levelset and levelset_norm parameters from
    s_initialize_mpi_data() subroutine
  • Removed MPI type creation for levelset and levelset norm data
    structures
  • Removed airfoil grid MPI I/O setup code
  • Simplified MPI data initialization for IB markers only
+2/-66   
m_initial_condition.fpp
Remove IBM-related variables from preprocessing                   

src/pre_process/m_initial_condition.fpp

  • Removed m_ib_patches module import
  • Removed ib_markers, levelset, and levelset_norm variable declarations
  • Removed IB patch allocation and initialization code
  • Removed s_apply_ib_patches() call from preprocessing
+1/-31   
m_time_steppers.fpp
Remove levelset parameters from immersed boundary update 

src/simulation/m_time_steppers.fpp

  • Simplified the s_update_mib subroutine call by removing levelset and
    levelset_norm parameters
  • Reflects the refactoring where levelset values are now computed in
    simulation rather than passed as arguments
+1/-1     
Enhancement
5 files
m_compute_levelset.fpp
New module for ghost-point-based levelset computation       

src/simulation/m_compute_levelset.fpp

  • New module created to handle all levelset computations for immersed
    boundaries
  • Implements s_apply_levelset() subroutine that computes levelsets for
    ghost points instead of grid cells
  • Contains geometry-specific levelset functions for circles, rectangles,
    ellipses, spheres, cylinders, airfoils, and STL models
  • Levelset and normal vector computations now tied to ghost_point
    derived type instead of global arrays
  • Parallelized over ghost points using GPU macros for improved memory
    efficiency
+723/-0 
m_data_output.fpp
Add dedicated IB data output subroutines                                 

src/simulation/m_data_output.fpp

  • Added new subroutines s_write_serial_ib_data() and
    s_write_parallel_ib_data() for IB data output
  • Created s_write_ib_data_file() wrapper to handle both serial and
    parallel I/O
  • Removed levelset and levelset norm output from main data writing
    routines
  • Simplified MPI data initialization calls by removing levelset
    parameters
+96/-21 
m_global_parameters.fpp
Add default initialization for IB patch parameters             

src/simulation/m_global_parameters.fpp

  • Added initialization loop for patch_ib() array with default values for
    all IB patch parameters
  • Sets proper defaults for STL model transformation parameters (scale,
    translate, rotate)
  • Initializes rotation matrices and moving boundary parameters
+41/-0   
m_derived_types.fpp
Add derived types for model storage and ghost point levelsets

src/common/m_derived_types.fpp

  • Added new t_model_array derived type to store STL model data with
    boundary vertices and interpolation information
  • Extended ghost_point derived type with levelset and levelset_norm
    fields for local storage
  • Removed global levelset array dependencies by moving data to ghost
    point structure
+13/-0   
case_dicts.py
Add STL model and transformation parameters to IB configuration

toolchain/mfc/run/case_dicts.py

  • Added three new immersed boundary patch parameters: model_filepath,
    model_spc, and model_threshold
  • Added three new model transformation parameters: model_translate and
    model_rotate for each spatial direction
  • Commented out model_scale parameter for potential future use
+6/-1     
Tests
7 files
cases.py
Adjust grid resolution for circle test cases                         

toolchain/mfc/test/cases.py

  • Updated Circle test case to use n=49 grid resolution for improved
    precision
  • Added comment explaining machine-level precision sensitivity for
    circular geometries
  • Minor formatting adjustment to boundary condition test setup
+6/-2     
golden.txt
Remove IB markers from golden test output                               

tests/7FA04E95/golden.txt

  • Removed the ib_markers output line from golden test data
  • Reflects the refactoring where IB markers are no longer pre-computed
    and stored globally
+1/-2     
golden.txt
Remove IB markers output from test golden file                     

tests/5600D63B/golden.txt

  • Removed the last line containing D/ib_markers.00.dat with all zero
    values
  • Kept the D/cons.5.00.000050.dat line with updated values
  • Test golden file updated to reflect removal of IB markers from output
+1/-2     
golden-metadata.txt
Update test metadata with new build environment details   

tests/7F70E665/golden-metadata.txt

  • Updated timestamp from 2025-01-20 to 2026-02-03
  • Changed Git commit hash and branch information
  • Updated CMake version and compiler information (AppleClang to GNU)
  • Changed system from macOS (Apple M1 Pro) to Linux (Intel i7-12700K)
  • Reordered build configuration sections (post_process moved to top)
  • Added OpenMP : OFF configuration line
  • Expanded CPU information with detailed lscpu output
+84/-41 
golden-metadata.txt
Update test metadata with new build environment details   

tests/F60D6594/golden-metadata.txt

  • Updated timestamp from 2025-01-20 to 2026-02-03
  • Changed Git commit hash and branch information
  • Updated CMake version and compiler information (AppleClang to GNU)
  • Changed system from macOS (Apple M1 Pro) to Linux (Intel i7-12700K)
  • Reordered build configuration sections
  • Added OpenMP : OFF configuration line
  • Expanded CPU information with detailed lscpu output
+84/-41 
golden-metadata.txt
Update test metadata with new build environment details   

tests/4F5A5E32/golden-metadata.txt

  • Updated timestamp from 2025-01-20 to 2026-02-03
  • Changed Git commit hash and branch information
  • Updated CMake version and compiler information (AppleClang to GNU)
  • Changed system from macOS (Apple M1 Pro) to Linux (Intel i7-12700K)
  • Reordered build configuration sections
  • Added OpenMP : OFF configuration line
  • Expanded CPU information with detailed lscpu output
+81/-38 
golden-metadata.txt
Update test metadata with new build environment details   

tests/8D8F6424/golden-metadata.txt

  • Updated timestamp from 2025-01-20 to 2026-02-03
  • Changed Git commit hash and branch information
  • Updated CMake version and compiler information (AppleClang to GNU)
  • Changed system from macOS (Apple M1 Pro) to Linux (Intel i7-12700K)
  • Reordered build configuration sections
  • Added OpenMP : OFF configuration line
  • Expanded CPU information with detailed lscpu output
+81/-38 
Miscellaneous
3 files
golden-metadata.txt
Update golden test metadata and build environment               

tests/B0CE19C5/golden-metadata.txt

  • Updated test metadata with new generation timestamp and git commit
    hash
  • Changed build system configuration details (CMake version, compiler
    paths, system information)
  • Updated CPU information from macOS to Linux system specifications
+93/-102
golden-metadata.txt
Update golden test metadata and build environment               

tests/7DCE34B4/golden-metadata.txt

  • Updated test metadata with new generation timestamp and git commit
    hash
  • Changed build system configuration details (CMake version, compiler
    paths, system information)
  • Updated CPU information from macOS to Linux system specifications
+90/-99 
golden-metadata.txt
Update golden test metadata and build environment               

tests/6171E9D4/golden-metadata.txt

  • Updated test metadata with new generation timestamp and git commit
    hash
  • Changed build system configuration details (CMake version, compiler
    paths, system information)
  • Updated CPU information from macOS M1 Pro to Linux x86_64 system
    specifications
+84/-41 
Additional files
63 files
m_compute_levelset.fpp +0/-625 
m_model.fpp +0/-1     
m_global_parameters.fpp +0/-10   
m_icpp_patches.fpp +0/-6     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-3     
golden.txt +0/-1     
golden.txt +0/-1     
golden-metadata.txt +78/-35 
golden.txt +10/-11 
golden.txt +0/-1     
golden.txt +1/-2     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +1/-2     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +10/-11 
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +10/-11 
golden-metadata.txt +81/-38 
golden.txt +14/-15 
golden-metadata.txt +81/-38 
golden.txt +14/-15 
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +18/-19 
golden.txt +12/-13 
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +18/-19 
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +1/-2     
golden.txt +0/-1     
golden.txt +12/-13 
golden-metadata.txt +81/-38 
golden.txt +10/-11 
golden.txt +1/-2     
golden.txt +1/-2     
golden.txt +0/-1     
golden-metadata.txt +81/-38 
golden.txt +14/-15 
golden.txt +1/-2     
golden.txt +1/-2     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +0/-1     
golden.txt +14/-15 
golden.txt +1/-2     

Summary by CodeRabbit

  • New Features

    • Improved STL/model-based immersed-boundary support and a new dedicated IB data writer; per-ghost-point level-set data now driven by models.
  • Refactor

    • Simplified startup and I/O: removed legacy per-patch level-set/marker file plumbing and centralized model-driven IB initialization with sane defaults.
  • Bug Fixes

    • Streamlined MPI/I/O and reduced redundant branches for IB workflows.
  • Tests

    • Golden outputs and metadata updated to reflect the new IB/model data formats.

CodeAnt-AI Description

Compute levelset and normals per ghost point and preload STL models; remove global levelset arrays

What Changed

  • Levelset distances and normals are computed for each ghost point on demand and stored on the ghost point structure; global levelset and levelset_norm arrays and their allocations/IO parameters were removed
  • STL models are loaded and transformed once during setup and reused for IB marking and levelset queries instead of re-reading per-cell; IB marker assignment now uses the preloaded models
  • IB marker data is written at startup and included in restart output (serial and MPI paths); many read/write interfaces no longer require separate levelset arguments
  • Image-point computation now reads distance and normal from the ghost point; when an image point falls outside the computational box the code prints a detailed diagnostic and aborts

Impact

✅ Lower memory during simulation
✅ IB markers written at startup and in restarts
✅ Shorter STL setup time on startup

💡 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:

@codeant-ai ask: Your question here

This lets you have a chat with CodeAnt AI about your pull request, making it easier to understand and improve your code.

Example

@codeant-ai ask: Can you suggest a safer alternative to storing this secret?

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:

@codeant-ai: Your feedback here

This helps CodeAnt AI learn and adapt to your team's coding style and standards.

Example

@codeant-ai: Do not flag unused imports.

Retrigger review

Ask CodeAnt AI to review the PR again, by typing:

@codeant-ai: review

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.

danieljvickers and others added 30 commits December 12, 2025 13:48
Copy link
Contributor

Choose a reason for hiding this comment

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

High-level Suggestion

To improve modularity, move the STL model handling logic (reading, transforming, interpolating) from m_ib_patches to the m_model module. This will decouple the immersed boundary implementation from the specifics of STL models. [High-level, importance: 7]

Solution Walkthrough:

Before:

! In src/simulation/m_ib_patches.fpp
module m_ib_patches
    ...
    type(t_model_array), allocatable, target :: models(:)
    ...
contains
    subroutine s_instantiate_STL_models()
        ! Reads STL file
        model = f_model_read(...)
        ...
        ! Transforms model
        call s_transform_model(model, ...)
        ...
        ! Interpolates model
        call f_interpolate_3D(model, ...)
        ...
        ! Stores prepared model data in the global `models` array
        models(patch_id)%model = model
        ...
    end subroutine s_instantiate_STL_models
    ...
end module m_ib_patches

After:

! In src/simulation/m_model.fpp (conceptual)
module m_model
    ...
    ! The global models array is now owned by the model module
    type(t_model_array), allocatable, target :: models(:)
    ...
contains
    subroutine s_instantiate_STL_models()
        ! All logic for reading, transforming, and preparing
        ! STL models is now encapsulated within the m_model module.
        ...
    end subroutine s_instantiate_STL_models
    ...
end module m_model

! In src/simulation/m_ib_patches.fpp
module m_ib_patches
    use m_model, only: models ! Use the global models array
    ...
end module m_ib_patches

Comment on lines 642 to 721
subroutine s_model_levelset(gp)
$:GPU_ROUTINE(parallelism='[seq]')

type(ghost_point), intent(inout) :: gp

integer :: i, j, k, patch_id, boundary_edge_count, total_vertices
type(t_model) :: model
real(wp), pointer, dimension(:, :, :) :: boundary_v
real(wp), pointer, dimension(:, :) :: interpolated_boundary_v
logical :: interpolate
real(wp), dimension(1:3) :: point
real(wp) :: normals(1:3) !< Boundary normal buffer
real(wp) :: distance

patch_id = gp%ib_patch_id
i = gp%loc(1)
j = gp%loc(2)
k = gp%loc(3)

! load in model values
model = models(patch_id)%model
interpolate = models(patch_id)%interpolate
boundary_edge_count = models(patch_id)%boundary_edge_count
total_vertices = models(patch_id)%total_vertices
boundary_v => models(patch_id)%boundary_v
interpolated_boundary_v => models(patch_id)%interpolated_boundary_v

! determine where we are located in space
point = (/x_cc(i), y_cc(j), 0._wp/)
if (p > 0) then
point(3) = z_cc(k)
end if

if (grid_geometry == 3) then
point = f_convert_cyl_to_cart(point)
end if

! 3D models
if (p > 0) then

! Get the boundary normals and shortest distance between the cell center and the model boundary
call f_distance_normals_3D(model, point, normals, distance)

! Get the shortest distance between the cell center and the interpolated model boundary
if (interpolate) then
gp%levelset = f_interpolated_distance(interpolated_boundary_v, total_vertices, point)
else
gp%levelset = distance
end if

! Correct the sign of the levelset
gp%levelset = -abs(gp%levelset)

! Assign the levelset_norm
gp%levelset_norm = normals(1:3)
else
! 2D models
if (interpolate) then
! Get the shortest distance between the cell center and the model boundary
gp%levelset = f_interpolated_distance(interpolated_boundary_v, total_vertices, point)
else
! Get the shortest distance between the cell center and the interpolated model boundary
gp%levelset = f_distance(boundary_v, boundary_edge_count, point)
end if

! Correct the sign of the levelset
gp%levelset = -abs(gp%levelset)

! Get the boundary normals
call f_normals(boundary_v, &
boundary_edge_count, &
point, &
normals)

! Assign the levelset_norm
gp%levelset_norm = normals(1:3)

end if

end subroutine s_model_levelset
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: In s_model_levelset, nullify the local pointers boundary_v and interpolated_boundary_v before the subroutine exits to prevent potential dangling pointers. [general, importance: 4]

Suggested change
subroutine s_model_levelset(gp)
$:GPU_ROUTINE(parallelism='[seq]')
type(ghost_point), intent(inout) :: gp
integer :: i, j, k, patch_id, boundary_edge_count, total_vertices
type(t_model) :: model
real(wp), pointer, dimension(:, :, :) :: boundary_v
real(wp), pointer, dimension(:, :) :: interpolated_boundary_v
logical :: interpolate
real(wp), dimension(1:3) :: point
real(wp) :: normals(1:3) !< Boundary normal buffer
real(wp) :: distance
patch_id = gp%ib_patch_id
i = gp%loc(1)
j = gp%loc(2)
k = gp%loc(3)
! load in model values
model = models(patch_id)%model
interpolate = models(patch_id)%interpolate
boundary_edge_count = models(patch_id)%boundary_edge_count
total_vertices = models(patch_id)%total_vertices
boundary_v => models(patch_id)%boundary_v
interpolated_boundary_v => models(patch_id)%interpolated_boundary_v
! determine where we are located in space
point = (/x_cc(i), y_cc(j), 0._wp/)
if (p > 0) then
point(3) = z_cc(k)
end if
if (grid_geometry == 3) then
point = f_convert_cyl_to_cart(point)
end if
! 3D models
if (p > 0) then
! Get the boundary normals and shortest distance between the cell center and the model boundary
call f_distance_normals_3D(model, point, normals, distance)
! Get the shortest distance between the cell center and the interpolated model boundary
if (interpolate) then
gp%levelset = f_interpolated_distance(interpolated_boundary_v, total_vertices, point)
else
gp%levelset = distance
end if
! Correct the sign of the levelset
gp%levelset = -abs(gp%levelset)
! Assign the levelset_norm
gp%levelset_norm = normals(1:3)
else
! 2D models
if (interpolate) then
! Get the shortest distance between the cell center and the model boundary
gp%levelset = f_interpolated_distance(interpolated_boundary_v, total_vertices, point)
else
! Get the shortest distance between the cell center and the interpolated model boundary
gp%levelset = f_distance(boundary_v, boundary_edge_count, point)
end if
! Correct the sign of the levelset
gp%levelset = -abs(gp%levelset)
! Get the boundary normals
call f_normals(boundary_v, &
boundary_edge_count, &
point, &
normals)
! Assign the levelset_norm
gp%levelset_norm = normals(1:3)
end if
end subroutine s_model_levelset
subroutine s_model_levelset(gp)
$:GPU_ROUTINE(parallelism='[seq]')
type(ghost_point), intent(inout) :: gp
...
! load in model values
model = models(patch_id)%model
interpolate = models(patch_id)%interpolate
boundary_edge_count = models(patch_id)%boundary_edge_count
total_vertices = models(patch_id)%total_vertices
boundary_v => models(patch_id)%boundary_v
interpolated_boundary_v => models(patch_id)%interpolated_boundary_v
! determine where we are located in space
point = (/x_cc(i), y_cc(j), 0._wp/)
...
! Assign the levelset_norm
gp%levelset_norm = normals(1:3)
end if
! De-associate pointers
nullify(boundary_v)
nullify(interpolated_boundary_v)
end subroutine s_model_levelset

Comment on lines 917 to 918
$:GPU_PARALLEL_LOOP(private='[i,j,k,point,cell_num]', copy='[ib_markers_sf]',&
& copyin='[patch_id,x_cc,y_cc,ncells,model]', collapse=3)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Add z_cc, dx, dy, and dz to the copyin clause of the $:GPU_PARALLEL_LOOP in s_ib_model to make them available on the GPU. [possible issue, importance: 9]

Suggested change
$:GPU_PARALLEL_LOOP(private='[i,j,k,point,cell_num]', copy='[ib_markers_sf]',&
& copyin='[patch_id,x_cc,y_cc,ncells,model]', collapse=3)
$:GPU_PARALLEL_LOOP(private='[i,j,k,point,cell_num]', copy='[ib_markers_sf]',&
& copyin='[patch_id,x_cc,y_cc,z_cc,dx,dy,dz,ncells,model]', collapse=3)

Comment on lines 1157 to 1165
impure subroutine s_finalize_ibm_module()

@:DEALLOCATE(ib_markers%sf)
@:DEALLOCATE(levelset%sf)
@:DEALLOCATE(levelset_norm%sf)
if (allocated(airfoil_grid_u)) then
@:DEALLOCATE(airfoil_grid_u)
@:DEALLOCATE(airfoil_grid_l)
end if

end subroutine s_finalize_ibm_module
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: In s_finalize_ibm_module, deallocate the globally allocated models array to prevent a memory leak. [general, importance: 7]

Suggested change
impure subroutine s_finalize_ibm_module()
@:DEALLOCATE(ib_markers%sf)
@:DEALLOCATE(levelset%sf)
@:DEALLOCATE(levelset_norm%sf)
if (allocated(airfoil_grid_u)) then
@:DEALLOCATE(airfoil_grid_u)
@:DEALLOCATE(airfoil_grid_l)
end if
end subroutine s_finalize_ibm_module
impure subroutine s_finalize_ibm_module()
@:DEALLOCATE(ib_markers%sf)
if (allocated(airfoil_grid_u)) then
@:DEALLOCATE(airfoil_grid_u)
@:DEALLOCATE(airfoil_grid_l)
end if
if (allocated(models)) then
! Free each STL model if needed, e.g.: call f_model_free(models(i)%model)
deallocate(models)
end if
end subroutine s_finalize_ibm_module

SIMULATION[f"patch_ib({ib_id})%{real_attr}"] = ty

for dir_id in range(1, 4):
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic()
# SIMULATION[f"patch_ib({ib_id})%model_scale({dir_id})"] = ParamType.REAL
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Uncomment the model_scale parameter in the case dictionary to allow users to configure this value from the Python interface, aligning it with the new Fortran defaults. [general, importance: 7]

Suggested change
# SIMULATION[f"patch_ib({ib_id})%model_scale({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%model_scale({dir_id})"] = ParamType.REAL.analytic()

!! @param levelset Closest distance from each grid cell to IB
!! @param levelset_norm Vector pointing in the direction of the closest distance
impure subroutine s_compute_image_points(ghost_points_in, levelset_in, levelset_norm_in)
impure subroutine s_compute_image_points(ghost_points_in)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Refactor the error handling in s_compute_image_points to move print and error stop calls outside the $:GPU_PARALLEL_LOOP to ensure safety and stability on GPUs. [possible issue, importance: 8]

end type t_model

type :: t_model_array
type(t_model), allocatable :: model
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: The component model in t_model_array is declared as a scalar allocatable rather than an allocatable array, which prevents this type from correctly representing multiple models; any calling code that expects to store or iterate over several models in this container will either silently only store one model or malfunction when trying to index it. [logic error]

Severity Level: Major ⚠️
- ❌ STL model loading routines rely on arrays.
- ⚠️ Preprocessing multiple-model workflows misbehave.
- ⚠️ ic_patch_parameters model handling impacted during setup.
Suggested change
type(t_model), allocatable :: model
type(t_model), allocatable, dimension(:) :: model
Steps of Reproduction ✅
1. Open the type definition at src/common/m_derived_types.fpp lines 186-193 (type
t_model_array definition).

2. Create a small test or locate a caller that intends to hold multiple models. Example
test (new file src/tests/allocate_model_test.f90):

   - use m_derived_types

   - type(t_model_array) :: t

   - integer :: n

   - n = 2

   - allocate(t%model(n)) ! attempt to allocate an array of models

3. Compile the test with the project's Fortran compiler (e.g., gfortran or nvfortran). The
compiler will reject the allocation or indexing:

   - Expected failure: "Error: Cannot allocate an object with bounds; 'model' is a scalar
   allocatable" or "subscripted object 'model' is not an array".

4. Alternatively, try to reference t%model(1) in code that iterates models; this will
produce a compile-time error "subscripted object is not an array" (observable in build
logs). This reproduces the malfunction concretely from the type definition at the cited
lines.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/m_derived_types.fpp
**Line:** 187:187
**Comment:**
	*Logic Error: The component `model` in `t_model_array` is declared as a scalar allocatable rather than an allocatable array, which prevents this type from correctly representing multiple models; any calling code that expects to store or iterate over several models in this container will either silently only store one model or malfunction when trying to index it.

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.

Comment on lines +451 to +452
real(wp) :: levelset
real(wp), dimension(1:3) :: levelset_norm
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: The newly added levelset and levelset_norm fields in ghost_point are not initialized, so any code that reads these values before they are explicitly set (e.g., for classification or interpolation) will operate on garbage data and may misclassify ghost points. [possible bug]

Severity Level: Critical 🚨
- ❌ Ghost-point classification may be incorrect.
- ❌ Levelset-based IB forces/markers miscomputed.
- ⚠️ Simulation results nondeterministic across runs.
Suggested change
real(wp) :: levelset
real(wp), dimension(1:3) :: levelset_norm
real(wp) :: levelset = 0.0_wp
real(wp), dimension(1:3) :: levelset_norm = 0.0_wp
Steps of Reproduction ✅
1. Inspect the ghost_point type in src/common/m_derived_types.fpp lines 448-455; lines
451-452 add levelset and levelset_norm fields.

2. Create a minimal reproducer (or find IB code that allocates ghost points) that declares
a ghost_point variable and reads the fields before any assignment:

   - use m_derived_types

   - type(ghost_point) :: g

   - print *, g%levelset, g%levelset_norm

3. Compile and run the reproducer. The printed values will be indeterminate
(uninitialized), demonstrating that reading these fields without prior initialization
yields garbage or non-deterministic values.

4. In a real simulation run, any IB/ghost-point routine that reads g%levelset or
g%levelset_norm for classification or interpolation before explicitly writing them (for
example, ghost-point setup or levelset-based boundary condition code paths) will make
wrong decisions; verify by locating IB routines that use ghost_point fields after
allocation and before explicit assignment.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/m_derived_types.fpp
**Line:** 451:452
**Comment:**
	*Possible Bug: The newly added `levelset` and `levelset_norm` fields in `ghost_point` are not initialized, so any code that reads these values before they are explicitly set (e.g., for classification or interpolation) will operate on garbage data and may misclassify ghost points.

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.

@@ -428,13 +535,13 @@ contains
end if

$:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,f]', copy='[ib_markers_sf]',&
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: In s_ib_3D_airfoil, the variables z_min and z_max are only assigned inside the if (.not. allocated(airfoil_grid_u)) block, so on subsequent calls when the grids are already allocated they remain uninitialized but are still used in the GPU loop condition, causing undefined behavior and incorrect IB marker placement; they should be set on every call. [logic error]

Severity Level: Major ⚠️
- ❌ Incorrect 3D airfoil IB markers placement.
- ⚠️ Time-stepped updates may produce nondeterministic results.
- ⚠️ Affects simulations using 3D airfoil patches.
- ⚠️ Hard-to-debug geometry errors during runs.
Suggested change
$:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,f]', copy='[ib_markers_sf]',&
z_max = lz/2
z_min = -lz/2
if (.not. allocated(airfoil_grid_u)) then
allocate (airfoil_grid_u(1:Np))
allocate (airfoil_grid_l(1:Np))
airfoil_grid_u(1)%x = 0._wp
airfoil_grid_u(1)%y = 0._wp
airfoil_grid_l(1)%x = 0._wp
airfoil_grid_l(1)%y = 0._wp
Steps of Reproduction ✅
1. Open src/simulation/m_ib_patches.fpp and find s_ib_3D_airfoil (subroutine header
earlier in the file). Local variables include "real(wp) :: lz, z_max, z_min" declared at
the subroutine start.

2. First call to s_ib_3D_airfoil (when airfoil grids unallocated) executes the block
starting at line 535: the code allocates airfoil_grid_u/l and sets z_max = lz/2 and z_min
= -lz/2 inside that allocation block.

3. On subsequent calls (common: s_apply_ib_patches calls per timestep), (.not.
allocated(airfoil_grid_u)) is false, so the allocation block is skipped and z_min/z_max
are not re-assigned in this call. Because z_min and z_max are local (stack) variables,
they are effectively uninitialized on these later calls.

4. The GPU parallel loop (starting around line 537) then uses z_min/z_max in the axial
extent check (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max). This leads to undefined
behavior or incorrect IB marker placement for 3D airfoils on repeated updates.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_ib_patches.fpp
**Line:** 537:537
**Comment:**
	*Logic Error: In `s_ib_3D_airfoil`, the variables `z_min` and `z_max` are only assigned inside the `if (.not. allocated(airfoil_grid_u))` block, so on subsequent calls when the grids are already allocated they remain uninitialized but are still used in the GPU loop condition, causing undefined behavior and incorrect IB marker placement; they should be set on every call.

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
Copy link
Contributor

codeant-ai bot commented Feb 4, 2026

CodeAnt AI finished reviewing your PR.

Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot left a comment

Choose a reason for hiding this comment

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

5 issues found across 86 files

Note: This PR contains a large number of files. cubic only reviews up to 75 files per PR, so some files may not have been reviewed.

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_compute_levelset.fpp">

<violation number="1" location="src/simulation/m_compute_levelset.fpp:239">
P1: Incorrect comparison in local coordinate system. After transforming to local coordinates (centered at origin), `xyz_local(2)` should be compared to `0._wp`, not `center(2)`. The 2D version `s_airfoil_levelset` correctly uses `xy_local(2) >= 0._wp` for this check.</violation>
</file>

<file name="toolchain/mfc/run/case_dicts.py">

<violation number="1" location="toolchain/mfc/run/case_dicts.py:378">
P2: SIMULATION schema omits model_scale for IB models. Simulation uses patch_ib%model_scale when reading STL models, so validation will reject valid configs that set model_scale. Add model_scale to the SIMULATION schema alongside translate/rotate.</violation>
</file>

<file name="src/simulation/m_ibm.fpp">

<violation number="1" location="src/simulation/m_ibm.fpp:126">
P0: This `GPU_UPDATE` overwrites the `ghost_points` data on the device (calculated in `s_apply_levelset`) with stale data from the host. 

As with `s_ibm_setup`, `s_apply_levelset` computes values on the GPU that are needed by `s_compute_image_points` on the host. The direction of the update should be from device to host.</violation>

<violation number="2" location="src/simulation/m_ibm.fpp:594">
P2: Remove this debug print statement. It is inside a loop over the entire domain and will generate excessive output, severely impacting performance and cluttering the logs.</violation>
</file>

<file name="src/simulation/m_data_output.fpp">

<violation number="1" location="src/simulation/m_data_output.fpp:1049">
P0: This `write` statement uses a format string `(A,I0,A,I0)` that requires 4 arguments, but only 1 argument is provided (`trim(case_dir)//'/p_all'`). This will likely cause a runtime error. Additionally, the line appears to be redundant as `t_step_dir` is overwritten immediately after.</violation>
</file>

Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.

character(LEN=path_len + 2*name_len) :: t_step_dir

! Creating or overwriting the time-step root directory
write (t_step_dir, '(A,I0,A,I0)') trim(case_dir)//'/p_all'
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Feb 4, 2026

Choose a reason for hiding this comment

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

P0: This write statement uses a format string (A,I0,A,I0) that requires 4 arguments, but only 1 argument is provided (trim(case_dir)//'/p_all'). This will likely cause a runtime error. Additionally, the line appears to be redundant as t_step_dir is overwritten immediately after.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_data_output.fpp, line 1049:

<comment>This `write` statement uses a format string `(A,I0,A,I0)` that requires 4 arguments, but only 1 argument is provided (`trim(case_dir)//'/p_all'`). This will likely cause a runtime error. Additionally, the line appears to be redundant as `t_step_dir` is overwritten immediately after.</comment>

<file context>
@@ -1016,26 +1018,99 @@ contains
+        character(LEN=path_len + 2*name_len) :: t_step_dir
+
+        ! Creating or overwriting the time-step root directory
+        write (t_step_dir, '(A,I0,A,I0)') trim(case_dir)//'/p_all'
+        write (t_step_dir, '(a,i0,a,i0)') trim(case_dir)//'/p_all/p', &
+            proc_rank, '/', time_step
</file context>
Fix with Cubic

patch_id = ib_markers%sf(i, j, 0)
ghost_points_in(count)%ib_patch_id = &
patch_id
print *, "patch_id: ", patch_id
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Feb 4, 2026

Choose a reason for hiding this comment

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

P2: Remove this debug print statement. It is inside a loop over the entire domain and will generate excessive output, severely impacting performance and cluttering the logs.

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 594:

<comment>Remove this debug print statement. It is inside a loop over the entire domain and will generate excessive output, severely impacting performance and cluttering the logs.</comment>

<file context>
@@ -598,6 +591,7 @@ contains
                             patch_id = ib_markers%sf(i, j, 0)
                             ghost_points_in(count)%ib_patch_id = &
                                 patch_id
+                            print *, "patch_id: ", patch_id
 
                             ghost_points_in(count)%slip = patch_ib(patch_id)%slip
</file context>
Fix with Cubic

Copy link
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 16

Caution

Some comments are outside the diff and can’t be posted inline due to platform limitations.

⚠️ Outside diff range comments (2)
tests/F60D6594/golden-metadata.txt (1)

9-143: ⚠️ Potential issue | 🟠 Major

Remove personal developer path from golden metadata.

The Fypp path contains a personal home directory: /home/dan/Documents/repos/MFC/build/venv/bin/fypp (lines 29, 63, 97, 131). This makes the golden metadata machine-specific and non-portable. Other test directories use system paths (e.g., /opt/homebrew/bin/gfortran, /usr/bin/cc), or the path should be resolved relative to the repository at test runtime.

Update the golden metadata to use the actual system path or a normalized reference that works across CI/CD environments and different developer machines.

src/simulation/m_ib_patches.fpp (1)

908-932: ⚠️ Potential issue | 🟠 Major

Add eta to private list and ensure f_model_is_inside has GPU-callable annotation.

eta is written and read within each GPU loop iteration but is not declared private, causing potential data races. Additionally, f_model_is_inside is called from GPU device code but has no $:GPU_ROUTINE annotation, risking device compilation issues.

🔧 Proposed fix
-        $:GPU_PARALLEL_LOOP(private='[i,j,k,point]', copy='[ib_markers_sf]',&
+        $:GPU_PARALLEL_LOOP(private='[i,j,k,point,eta]', copy='[ib_markers_sf]',&
                   & copyin='[patch_id,x_cc,y_cc,ncells,model]', collapse=3)

Additionally, verify that f_model_is_inside in src/common/m_model.fpp (line 488) is marked with $:GPU_ROUTINE(function_name='f_model_is_inside', parallelism='seq') immediately after its declaration, or mark it now if missing.

🤖 Fix all issues with AI agents
In `@src/common/m_derived_types.fpp`:
- Around line 186-193: The allocatable array components in type t_model_array
(boundary_v and interpolated_boundary_v) must be declared with the target
attribute so they can be safely associated by pointers in m_compute_levelset;
update the type t_model_array to add the target attribute to the allocatable
declarations for boundary_v and interpolated_boundary_v (keeping their existing
kinds/dimensions) so pointer association from models (which is declared target)
works correctly and is portable to device compilers.

In `@src/simulation/m_compute_levelset.fpp`:
- Around line 430-527: The cuboid levelset has ambiguity at corners and under
rotation; replace the cascade of side_dists checks in subroutine
s_cuboid_levelset with a robust oriented-box signed-distance calculation:
compute half_lengths = (length_x/2,length_y/2,length_z/2), q = abs(xyz_local) -
half_lengths, outside = norm(max(q,0)), inside = min(max(q(1),max(q(2),q(3))),
0), and set gp%levelset = outside + inside; compute the local normal as (if
outside > 0) normalize(max(q,0)) * sign(xyz_local) else unit vector along the
component with largest q (use sign(xyz_local) for direction), then rotate that
normal back with matmul(rotation, normal) into gp%levelset_norm; replace
references to side_dists/min_dist and ensure zero-division is guarded when
normal magnitude is ~0.0.
- Around line 638-721: s_model_levelset copies the whole derived-model into
local variable "model" which triggers deep copies of allocatable triangle data
(e.g., model%trs) per ghost point; change the local scalar declaration to a
pointer/associate to the target in "models" and use that pointer instead (e.g.,
replace "type(t_model) :: model" and assignments like "model =
models(patch_id)%model" with a pointer association such as "type(t_model),
pointer :: model_ptr" and "model_ptr => models(patch_id)%model") so the routine
reads the model by reference and avoids device-side deep copies; update all
subsequent uses of "model" to use the pointer name (model_ptr) or an equivalent
associate.
- Around line 29-88: The GPU parallel loops in s_apply_levelset declare patch_id
and patch_geometry as loop-local but they are not listed in the
GPU_PARALLEL_LOOP private clause, causing potential data races; update both
GPU_PARALLEL_LOOP directives (the one before the 3D loop and the one before the
2D loop) to include patch_id and patch_geometry in the private list so each
thread gets its own copy when executing the body that reads/writes
gps(i)%ib_patch_id and uses patch_ib(patch_id)%geometry.
- Around line 121-300: The GPU kernel routines s_airfoil_levelset and
s_3D_airfoil_levelset read the arrays airfoil_grid_u, airfoil_grid_l and the
scalar Np but the GPU_PARALLEL_LOOP in s_apply_levelset only copyin patch_ib
(with default(present)), causing missing-device-data faults; fix by adding
airfoil_grid_u, airfoil_grid_l and Np to the GPU data movement (e.g., include
them in the copyin clause of the GPU_PARALLEL_LOOP that calls s_apply_levelset)
or alternatively ensure they are moved to the device at startup via a persistent
$:GPU_ENTER_DATA(...) region (for example in m_start_up.fpp) so the kernels see
these symbols on the device.

In `@src/simulation/m_data_output.fpp`:
- Around line 1026-1033: The code writes ib.dat into a per-rank timestep
directory without ensuring the directory exists; update s_write_serial_ib_data
(and the parallel writer s_write_parallel_ib_data / caller s_write_ib_data_file)
to create the target timestep directory before opening the file: compute the
output path used (e.g. /p_all/p<rank>/<time_step> or the parallel path built in
s_write_parallel_ib_data), check for existence and create missing directories
including parents (mkdir -p semantics) and handle/report errors before
attempting MPI_FILE_OPEN or open; ensure the same directory-creation logic is
applied in both the serial and parallel paths so s_write_ib_data_file(0) cannot
fail due to a missing timestep directory.

In `@src/simulation/m_ibm.fpp`:
- Line 602: Remove the leftover debug output that prints patch_id for every 2D
ghost point: delete the print *, "patch_id: ", patch_id statement in the routine
(the line referencing patch_id) so normal simulation runs no longer flood the
console; if retaining optional diagnostics is desired instead, wrap that print
behind an existing debug/verbose flag check rather than leaving it
unconditional.
- Around line 501-509: The error stop call must be replaced with a call to
s_mpi_abort(...) passing a descriptive message; locate the block that prints
diagnostics (references: ghost_points_in, x_cc, y_cc, z_cc, buff_size,
m_helper_basic, dist, norm) and remove the final `error stop "Ghost Point and
Image Point on Different Processors"` and instead call s_mpi_abort with a clear
message (e.g., include "Ghost Point and Image Point on Different Processors"
plus any concise context you want), ensuring the s_mpi_abort call is used for
termination per guidelines.

In `@tests/1B218CF1/golden.txt`:
- Line 14: The golden baseline changed because the D/cons.7.00.000050.dat block
was regenerated from the new runtime IB/levelset path; update the PR and test
metadata to state that D/cons.7.00.000050.dat was intentionally regenerated from
the new IB/levelset runtime (include the regeneration command or script and the
runtime/version used), add/update focused tests around the cons.7 block
(targeting the code paths that produce D/cons.7.00.000050.dat, e.g., the
levelset/IB processing functions), and include a short justification entry in
the test/golden manifest so future diffs reference this change and prevent
unintended regressions.

In `@tests/34626101/golden.txt`:
- Line 10: The golden baseline change is a deliberate regeneration of the
D/cons.5.00.000050.dat block; update the PR and test metadata to record that
this file was regenerated from the updated runtime IB pipeline, add a short
justification note in tests/34626101/golden.txt (or the test's metadata section)
naming D/cons.5.00.000050.dat and the pipeline version, and add focused tests
that assert the specific behavior/values introduced by the new cons.5 block (so
regression coverage exists) — ensure the PR description explicitly states the
regeneration and links the new focused tests to the golden update.

In `@tests/6171E9D4/golden-metadata.txt`:
- Around line 145-193: The PR regenerated all golden metadata files on different
hardware (Intel x86_64, branch generate-new-golden-temp) without explaining why;
update the PR to state whether the regeneration was intentional, why it was
necessary for "Added 3D models to pregeneration step", and whether you
standardised environment metadata or accidentally committed local snapshots
(refer to golden-metadata.txt and the generate-new-golden-temp branch); amend
the commit message/PR description to document intent, add a note to the
test-generation script (e.g., generate_golden_metadata or whatever script
produced these files) describing when and how to regenerate metadata, and if
needed add CI/CONTRIBUTING guidance about hardware-dependent metadata or make
metadata generation machine-independent to avoid future churn.

In `@tests/7DCE34B4/golden-metadata.txt`:
- Around line 1-7: The golden-metadata.txt records a dirty working tree and a
temporary branch (git commit d987bdb5a38c912183f371ddff56c6b1642d05d7 on
generate-new-golden-temp), which reduces reproducibility; either (A) confirm
this rebaseline is intentional by updating the test metadata to remove or
annotate the “dirty” and “temp branch” warnings and add a short justification
note, or (B) regenerate the golden from a clean, stable ref: checkout a clean
commit/branch, run the same mfc.sh invocation (test --generate --only 7DCE34B4),
and commit the regenerated tests so golden-metadata.txt reflects a clean,
reproducible state.

In `@tests/C7BE8FEA/golden.txt`:
- Line 10: Regenerate the golden named by the entry "D/cons.5.00.000050.dat"
from the updated test case using the new runtime IB path and attach a brief note
in the PR describing the regeneration command and inputs; add focused
unit/integration tests that exercise the runtime IB/levelset changes (targeting
the code paths that produce "cons.5.00.000050" outputs) so the golden reflects
real behavior, and ensure the PR includes a short justification that the golden
change is expected (or revert if unintended) to prevent silent regressions.

In `@tests/C87A466A/golden.txt`:
- Line 10: The golden snapshot line for D/cons.5.00.000050.dat was changed
without documenting how it was regenerated; update the PR/test to include the
exact regeneration command and rationale (e.g., script or test harness command,
data source, random seed or versioned dependency) and state why the new numeric
values for D/cons.5.00.000050.dat are expected (what code change produced them
and why this is not a regression); reference the golden file name
D/cons.5.00.000050.dat and include the provenance (command, input fixtures,
seed, and expected behavior) in the commit message or test metadata so reviewers
can reproduce the golden file and verify there is no silent drift.

In `@toolchain/mfc/run/case_dicts.py`:
- Around line 371-373: PRE_PROCESS and SIMULATION schemas are inconsistent
because model_scale is commented out in the SIMULATION tuple; update the
SIMULATION parameter list (the tuple that currently includes ("moving_ibm",
ParamType.INT), ("mass", ParamType.REAL), ("model_filepath", ParamType.STR),
("model_spc", ParamType.INT), ("model_threshold", ParamType.REAL)) to include
("model_scale", ParamType.REAL) so the SIMULATION handling matches PRE_PROCESS
(or alternatively remove model_scale from PRE_PROCESS if you prefer that
direction); locate the tuple in case_dicts.py around the PRE_PROCESS/SIMULATION
definitions and add the ("model_scale", ParamType.REAL) entry to SIMULATION to
keep validation consistent.

In `@toolchain/mfc/test/cases.py`:
- Around line 415-418: The formatting for the dictionary entry in the
define_case_d call for Circle{suffix} is inconsistent: change the key/value pair
`'n' :49` to follow the file style (e.g., `'n': 49`) in the call to
define_case_d so spacing around the colon matches other cases and maintains
consistent formatting.
🧹 Nitpick comments (3)
tests/F60D6594/golden-metadata.txt (1)

145-193: Consider portability implications of expanded CPU metadata.

The CPU section has been significantly expanded to include detailed hardware information (architecture, cache topology, CPU flags, vulnerability mitigations). While this provides comprehensive system documentation, it introduces machine-specific details that may cause golden file validation issues when tests run on different hardware configurations.

Consider whether:

  1. This level of CPU detail is necessary for test F60D6594's validation logic
  2. Some fields (especially vulnerability mitigations, specific CPU flags) should be excluded or normalized to improve portability
  3. The test framework has tolerance for hardware-specific fields when comparing against golden metadata

If hardware-specific metadata is required for debugging but not validation, consider separating it into a different file or marking it as non-validated.

tests/8D8F6424/golden-metadata.txt (1)

147-193: Reconsider including exhaustive machine-specific CPU details in golden metadata.

The CPU section was completely replaced with 46 lines of detailed x86_64-specific information (cache topology, CPU flags, vulnerability mitigations, etc.). While capturing build environment is valuable, this level of machine-specific detail causes golden file churn whenever:

  • Tests regenerate on different developer machines
  • CI infrastructure changes
  • OS/kernel updates affect vulnerability reporting

Consider limiting golden metadata to test-relevant configuration:

  • Architecture (x86_64 vs ARM)
  • CPU vendor/model (if tests are performance-sensitive)
  • Core/thread count (if parallelism affects outputs)

Exclude volatile details like:

  • Specific vulnerability mitigations (change with kernel updates)
  • BogoMIPS, CPU flags, cache topologies
  • lscpu's full output

This reduces maintenance burden while preserving test validity. Based on learnings, golden files should capture meaningful test conditions without unnecessary machine-specific noise.

src/simulation/m_ib_patches.fpp (1)

144-176: Consider avoiding per-rank STL file reads to reduce filesystem contention.

Reading large STL files on every MPI rank can overwhelm shared storage. A rank‑0 read + broadcast (or MPI‑IO) would scale better.

Comment on lines +186 to +193
type :: t_model_array
type(t_model), allocatable :: model
real(wp), allocatable, dimension(:, :, :) :: boundary_v
real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v
integer :: boundary_edge_count
integer :: total_vertices
logical :: interpolate
end type t_model_array
Copy link
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

🧩 Analysis chain

🏁 Script executed:

# First, let's examine the exact code at lines 186-193 in m_derived_types.fpp
sed -n '186,193p' src/common/m_derived_types.fpp

Repository: MFlowCode/MFC

Length of output: 395


🏁 Script executed:

# Find and examine m_compute_levelset to see how it uses these components
fd -e fpp -e f90 | xargs grep -l "m_compute_levelset" | head -5

Repository: MFlowCode/MFC

Length of output: 119


🏁 Script executed:

# Search for pointer association patterns with boundary_v or interpolated_boundary_v
rg "boundary_v\s*=>" src/

Repository: MFlowCode/MFC

Length of output: 258


🏁 Script executed:

# Check if TARGET is already on allocatable components in the actual file
rg -A 10 "type :: t_model_array" src/common/m_derived_types.fpp

Repository: MFlowCode/MFC

Length of output: 553


🏁 Script executed:

# Examine the context around pointer associations in m_compute_levelset
rg -B 5 -A 5 "boundary_v => models" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 577


🏁 Script executed:

# Look for pointer declarations in m_compute_levelset
rg "pointer.*::.*boundary_v" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 186


🏁 Script executed:

# Get broader context - find the subroutine containing these assignments
rg -B 20 "boundary_v => models" src/simulation/m_compute_levelset.fpp | head -60

Repository: MFlowCode/MFC

Length of output: 926


🏁 Script executed:

# Check if other allocatable components in derived types use TARGET
rg "allocatable.*target|target.*allocatable" src/common/m_derived_types.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Search for other derived type definitions with allocatable components
rg -A 8 "type ::" src/common/m_derived_types.fpp | grep -A 8 "allocatable" | head -40

Repository: MFlowCode/MFC

Length of output: 479


🏁 Script executed:

# Check the coding guidelines file to see if TARGET is required for allocatable in derived types
cat .cursor/rules/mfc-agent-rules.mdc | head -100

Repository: MFlowCode/MFC

Length of output: 3968


🏁 Script executed:

# Check if models (the array of t_model_array) is declared with TARGET
rg -B 5 -A 5 "models.*t_model_array|allocatable.*models" src/simulation/m_compute_levelset.fpp | head -40

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Search for how models is declared - find the variable declaration
rg "type.*t_model_array.*::|:.*models" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Look for the full context of where models is declared/allocated
rg -B 10 "models(patch_id)" src/simulation/m_compute_levelset.fpp | head -60

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Find where models is passed as a parameter or declared in m_compute_levelset
rg "subroutine.*models|function.*models" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Look for all subroutine/function signatures that use models
ast-grep --pattern 'subroutine $FUNC($$$models$$$)
  $$$
end subroutine'

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Find the subroutine containing the pointer associations
rg -B 50 "boundary_v => models" src/simulation/m_compute_levelset.fpp | grep "subroutine\|function" | tail -1

Repository: MFlowCode/MFC

Length of output: 92


🏁 Script executed:

# Get the full subroutine signature to see if models is a parameter
sed -n '/subroutine.*compute_levelset/,/end subroutine/p' src/simulation/m_compute_levelset.fpp | head -80

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Find the full subroutine s_model_levelset with its signature and initial context
rg -A 60 "subroutine s_model_levelset" src/simulation/m_compute_levelset.fpp | head -80

Repository: MFlowCode/MFC

Length of output: 2385


🏁 Script executed:

# Check what type gp is and how models is accessed (likely via module-level variable)
rg "type.*::|intent.*models" src/simulation/m_compute_levelset.fpp | head -20

Repository: MFlowCode/MFC

Length of output: 572


🏁 Script executed:

# Search for module-level declarations in m_compute_levelset
rg "^[[:space:]]*(type|real|integer|logical|character).*models|^[[:space:]]*allocatable.*models" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Look at the top of m_compute_levelset module for variable declarations
head -100 src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 3212


🏁 Script executed:

# Search for where models is declared or imported - it's used but not a parameter
rg "models|allocate.*models|type.*models" src/simulation/m_compute_levelset.fpp | grep -v "patch_models\|s_model_\|f_model_\|t_model"

Repository: MFlowCode/MFC

Length of output: 440


🏁 Script executed:

# Search in m_ib_patches which is used in m_compute_levelset
rg "models" src/common/m_ib_patches.fpp | head -20

Repository: MFlowCode/MFC

Length of output: 124


🏁 Script executed:

# Look at m_global_parameters to see if models is declared there
rg "models" src/common/m_global_parameters.fpp | head -20

Repository: MFlowCode/MFC

Length of output: 131


🏁 Script executed:

# Search broadly for where models array is declared as module variable
rg "^[[:space:]]*(allocatable|target).*models|type.*t_model_array.*models" src/

Repository: MFlowCode/MFC

Length of output: 146


🏁 Script executed:

# Verify the coding guideline reference about pointer encapsulation
grep -A 5 "pointer.*encapsulation\|encapsulation.*pointer" .cursor/rules/mfc-agent-rules.mdc

Repository: MFlowCode/MFC

Length of output: 39


Add target attribute to allocatable components in t_model_array.

The boundary_v and interpolated_boundary_v allocatable components are pointer targets in m_compute_levelset. Since the parent models variable itself is declared with target, these components should also be explicitly marked as target for consistency, portability, and safety—especially for device code where some compilers enforce strict pointer-association rules.

🛠️ Suggested fix
 type :: t_model_array
-    type(t_model), allocatable :: model
-    real(wp), allocatable, dimension(:, :, :) :: boundary_v
-    real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v
+    type(t_model), allocatable, target :: model
+    real(wp), allocatable, target, dimension(:, :, :) :: boundary_v
+    real(wp), allocatable, target, dimension(:, :) :: interpolated_boundary_v
     integer :: boundary_edge_count
     integer :: total_vertices
     logical :: interpolate
 end type t_model_array
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
type :: t_model_array
type(t_model), allocatable :: model
real(wp), allocatable, dimension(:, :, :) :: boundary_v
real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v
integer :: boundary_edge_count
integer :: total_vertices
logical :: interpolate
end type t_model_array
type :: t_model_array
type(t_model), allocatable, target :: model
real(wp), allocatable, target, dimension(:, :, :) :: boundary_v
real(wp), allocatable, target, dimension(:, :) :: interpolated_boundary_v
integer :: boundary_edge_count
integer :: total_vertices
logical :: interpolate
end type t_model_array
🤖 Prompt for AI Agents
In `@src/common/m_derived_types.fpp` around lines 186 - 193, The allocatable array
components in type t_model_array (boundary_v and interpolated_boundary_v) must
be declared with the target attribute so they can be safely associated by
pointers in m_compute_levelset; update the type t_model_array to add the target
attribute to the allocatable declarations for boundary_v and
interpolated_boundary_v (keeping their existing kinds/dimensions) so pointer
association from models (which is declared target) works correctly and is
portable to device compilers.

Comment on lines 29 to 88
impure subroutine s_apply_levelset(gps, num_gps)

type(ghost_point), dimension(:), intent(inout) :: gps
integer, intent(in) :: num_gps

integer :: i, patch_id, patch_geometry

! 3D Patch Geometries
if (p > 0) then

$:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib]')
do i = 1, num_gps

patch_id = gps(i)%ib_patch_id
patch_geometry = patch_ib(patch_id)%geometry

if (patch_geometry == 8) then
call s_sphere_levelset(gps(i))
elseif (patch_geometry == 9) then
call s_cuboid_levelset(gps(i))
elseif (patch_geometry == 10) then
call s_cylinder_levelset(gps(i))
elseif (patch_geometry == 11) then
call s_3D_airfoil_levelset(gps(i))
! STL+IBM patch
elseif (patch_geometry == 12) then
call s_model_levelset(gps(i))
end if
end do
$:END_GPU_PARALLEL_LOOP()
!> @}

! 2D Patch Geometries
elseif (n > 0) then

$:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib]')
do i = 1, num_gps

patch_id = gps(i)%ib_patch_id
patch_geometry = patch_ib(patch_id)%geometry

if (patch_geometry == 2) then
call s_circle_levelset(gps(i))
elseif (patch_geometry == 3) then
call s_rectangle_levelset(gps(i))
elseif (patch_geometry == 4) then
call s_airfoil_levelset(gps(i))
! STL+IBM patch
elseif (patch_geometry == 5) then
call s_model_levelset(gps(i))
elseif (patch_geometry == 6) then
call s_ellipse_levelset(gps(i))
end if
end do
$:END_GPU_PARALLEL_LOOP()
!> @}

end if

end subroutine s_apply_levelset
Copy link
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

🧩 Analysis chain

🏁 Script executed:

cat -n src/simulation/m_compute_levelset.fpp | head -100

Repository: MFlowCode/MFC

Length of output: 3912


🏁 Script executed:

find . -name "parallel_macros.fpp" -o -name "macros.fpp" | head -5

Repository: MFlowCode/MFC

Length of output: 129


🏁 Script executed:

rg "GPU_PARALLEL_LOOP" --type fpp -A 3 | head -50

Repository: MFlowCode/MFC

Length of output: 84


🏁 Script executed:

rg "GPU_PARALLEL_LOOP" -A 5 | head -80

Repository: MFlowCode/MFC

Length of output: 5525


🏁 Script executed:

cat src/common/include/parallel_macros.fpp | head -150

Repository: MFlowCode/MFC

Length of output: 6753


🏁 Script executed:

cat src/common/include/acc_macros.fpp | grep -A 20 "def ACC_PARALLEL_LOOP"

Repository: MFlowCode/MFC

Length of output: 1484


🏁 Script executed:

cat src/common/include/acc_macros.fpp | grep -A 10 "def GEN_PRIVATE_STR"

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

cat src/common/include/shared_parallel_macros.fpp | grep -A 15 "def GEN_PRIVATE_STR"

Repository: MFlowCode/MFC

Length of output: 716


🏁 Script executed:

cat src/common/include/shared_parallel_macros.fpp | grep -A 15 "def GEN_PARENTHESES_CLAUSE"

Repository: MFlowCode/MFC

Length of output: 621


🏁 Script executed:

cat src/common/include/shared_parallel_macros.fpp | grep -A 20 "def GEN_CLAUSE"

Repository: MFlowCode/MFC

Length of output: 1002


🏁 Script executed:

cat -n src/simulation/m_compute_levelset.fpp | sed -n '29,88p'

Repository: MFlowCode/MFC

Length of output: 2577


Privatize patch_id and patch_geometry in both GPU loops.

These variables are assigned per-iteration but aren't in the private list, creating data races across GPU threads. Add them to the private parameter in both GPU_PARALLEL_LOOP directives (lines 39 and 64).

🛠️ Suggested fix
-            $:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib]')
+            $:GPU_PARALLEL_LOOP(private='[i,patch_id,patch_geometry]', copy='[gps]', copyin='[patch_ib]')

Apply the same change at line 64 for the 2D loop.

🤖 Prompt for AI Agents
In `@src/simulation/m_compute_levelset.fpp` around lines 29 - 88, The GPU parallel
loops in s_apply_levelset declare patch_id and patch_geometry as loop-local but
they are not listed in the GPU_PARALLEL_LOOP private clause, causing potential
data races; update both GPU_PARALLEL_LOOP directives (the one before the 3D loop
and the one before the 2D loop) to include patch_id and patch_geometry in the
private list so each thread gets its own copy when executing the body that
reads/writes gps(i)%ib_patch_id and uses patch_ib(patch_id)%geometry.

Comment on lines 121 to 300
subroutine s_airfoil_levelset(gp)
$:GPU_ROUTINE(parallelism='[seq]')

type(ghost_point), intent(inout) :: gp

real(wp) :: dist, global_dist
integer :: global_id
real(wp), dimension(3) :: dist_vec

real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
real(wp), dimension(1:2) :: center
real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation

integer :: i, j, k, ib_patch_id !< Loop index variables

ib_patch_id = gp%ib_patch_id
i = gp%loc(1)
j = gp%loc(2)

center(1) = patch_ib(ib_patch_id)%x_centroid
center(2) = patch_ib(ib_patch_id)%y_centroid
inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)

xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset

if (xy_local(2) >= 0._wp) then
! finds the location on the airfoil grid with the minimum distance (closest)
do k = 1, Np
dist_vec(1) = xy_local(1) - airfoil_grid_u(k)%x
dist_vec(2) = xy_local(2) - airfoil_grid_u(k)%y
dist_vec(3) = 0._wp
dist = sqrt(sum(dist_vec**2))
if (k == 1) then
global_dist = dist
global_id = k
else
if (dist < global_dist) then
global_dist = dist
global_id = k
end if
end if
end do
dist_vec(1) = xy_local(1) - airfoil_grid_u(global_id)%x
dist_vec(2) = xy_local(2) - airfoil_grid_u(global_id)%y
dist_vec(3) = 0
dist = global_dist
else
do k = 1, Np
dist_vec(1) = xy_local(1) - airfoil_grid_l(k)%x
dist_vec(2) = xy_local(2) - airfoil_grid_l(k)%y
dist_vec(3) = 0
dist = sqrt(sum(dist_vec**2))
if (k == 1) then
global_dist = dist
global_id = k
else
if (dist < global_dist) then
global_dist = dist
global_id = k
end if
end if
end do
dist_vec(1) = xy_local(1) - airfoil_grid_l(global_id)%x
dist_vec(2) = xy_local(2) - airfoil_grid_l(global_id)%y
dist_vec(3) = 0
dist = global_dist
end if

gp%levelset = dist
if (f_approx_equal(dist, 0._wp)) then
gp%levelset_norm = 0._wp
else
gp%levelset_norm = matmul(rotation, dist_vec(:))/dist ! convert the normal vector back to global grid coordinates
end if

end subroutine s_airfoil_levelset

subroutine s_3D_airfoil_levelset(gp)
$:GPU_ROUTINE(parallelism='[seq]')

type(ghost_point), intent(inout) :: gp

real(wp) :: dist, dist_surf, dist_side, global_dist
integer :: global_id
real(wp) :: lz, z_max, z_min
real(wp), dimension(3) :: dist_vec

real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame
real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation

real(wp) :: length_z

integer :: i, j, k, l, ib_patch_id !< Loop index variables

ib_patch_id = gp%ib_patch_id
i = gp%loc(1)
j = gp%loc(2)
l = gp%loc(3)

center(1) = patch_ib(ib_patch_id)%x_centroid
center(2) = patch_ib(ib_patch_id)%y_centroid
center(3) = patch_ib(ib_patch_id)%z_centroid
lz = patch_ib(ib_patch_id)%length_z
inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)

z_max = center(3) + lz/2
z_min = center(3) - lz/2

xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset

if (xyz_local(2) >= center(2)) then
do k = 1, Np
dist_vec(1) = xyz_local(1) - airfoil_grid_u(k)%x
dist_vec(2) = xyz_local(2) - airfoil_grid_u(k)%y
dist_vec(3) = 0
dist_surf = sqrt(sum(dist_vec**2))
if (k == 1) then
global_dist = dist_surf
global_id = k
else
if (dist_surf < global_dist) then
global_dist = dist_surf
global_id = k
end if
end if
end do
dist_vec(1) = xyz_local(1) - airfoil_grid_u(global_id)%x
dist_vec(2) = xyz_local(2) - airfoil_grid_u(global_id)%y
dist_vec(3) = 0
dist_surf = global_dist
else
do k = 1, Np
dist_vec(1) = xyz_local(1) - airfoil_grid_l(k)%x
dist_vec(2) = xyz_local(2) - airfoil_grid_l(k)%y
dist_vec(3) = 0
dist_surf = sqrt(sum(dist_vec**2))
if (k == 1) then
global_dist = dist_surf
global_id = k
else
if (dist_surf < global_dist) then
global_dist = dist_surf
global_id = k
end if
end if
end do
dist_vec(1) = xyz_local(1) - airfoil_grid_l(global_id)%x
dist_vec(2) = xyz_local(2) - airfoil_grid_l(global_id)%y
dist_vec(3) = 0
dist_surf = global_dist
end if

dist_side = min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l)))

if (dist_side < dist_surf) then
gp%levelset = dist_side
if (f_approx_equal(dist_side, abs(z_cc(l) - z_min))) then
gp%levelset_norm = (/0._wp, 0._wp, -1._wp/)
else
gp%levelset_norm = (/0._wp, 0._wp, 1._wp/)
end if
gp%levelset_norm = matmul(rotation, gp%levelset_norm)
else
gp%levelset = dist_surf
if (f_approx_equal(dist_surf, 0._wp)) then
gp%levelset_norm = 0._wp
else
gp%levelset_norm = matmul(rotation, dist_vec(:)/dist_surf)
end if
end if

end subroutine s_3D_airfoil_levelset
Copy link
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

🧩 Analysis chain

🏁 Script executed:

# First, let's find the file and check its overall structure
find . -name "m_compute_levelset.fpp" -type f

Repository: MFlowCode/MFC

Length of output: 96


🏁 Script executed:

# Let's examine the full m_compute_levelset.fpp file to find declarations and allocations
cat -n src/simulation/m_compute_levelset.fpp | head -150

Repository: MFlowCode/MFC

Length of output: 6070


🏁 Script executed:

# Search for where airfoil_grid_u, airfoil_grid_l, and Np are declared
rg "airfoil_grid_u|airfoil_grid_l|Np" src/simulation/m_compute_levelset.fpp -B 2 -A 2

Repository: MFlowCode/MFC

Length of output: 2174


🏁 Script executed:

# Look for allocate statements for these variables
rg "allocate.*airfoil" src/simulation/m_compute_levelset.fpp -B 2 -A 5

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Check if there are any enter data regions in the file
rg "enter data" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Look at the parallel macros to understand GPU_ROUTINE expansion
cat src/common/include/parallel_macros.fpp | head -200

Repository: MFlowCode/MFC

Length of output: 8148


🏁 Script executed:

# Get the full size of the file first
wc -l src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 98


🏁 Script executed:

# Search for where these variables might be declared in this file
cat src/simulation/m_compute_levelset.fpp | tail -n +150 | head -n 100

Repository: MFlowCode/MFC

Length of output: 3973


🏁 Script executed:

# Search for where airfoil_grid is declared - could be in another module
rg "airfoil_grid_u|airfoil_grid_l" src/ -l

Repository: MFlowCode/MFC

Length of output: 260


🏁 Script executed:

# Check the m_ib_patches module
find . -name "*ib_patches*" -type f

Repository: MFlowCode/MFC

Length of output: 131


🏁 Script executed:

# Look for Np declaration
rg "^.*Np" src/simulation/ -A 2 -B 2 | head -50

Repository: MFlowCode/MFC

Length of output: 3737


🏁 Script executed:

# Check ACC_ROUTINE macro expansion
rg "ACC_ROUTINE" src/common/include/ -A 5 | head -80

Repository: MFlowCode/MFC

Length of output: 1181


🏁 Script executed:

# Check the full declaration in m_global_parameters
rg "airfoil_grid_u|airfoil_grid_l|^.*integer.*Np" src/simulation/m_global_parameters.fpp -B 3 -A 3

Repository: MFlowCode/MFC

Length of output: 482


🏁 Script executed:

# Check if there's a managed attribute anywhere
rg "managed" src/simulation/ -i

Repository: MFlowCode/MFC

Length of output: 343


🏁 Script executed:

# Look for enter data regions in m_compute_levelset.fpp
rg "GPU_ENTER_DATA|GPU_EXIT_DATA|enter data" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Check the allocate calls to see if they use managed or other attributes
rg "allocate.*airfoil" src/simulation/m_start_up.fpp -B 2 -A 2

Repository: MFlowCode/MFC

Length of output: 558


🏁 Script executed:

# Look at acc_macros to understand ACC_ROUTINE expansion
cat src/common/include/acc_macros.fpp | head -150

Repository: MFlowCode/MFC

Length of output: 5881


🏁 Script executed:

# Verify the allocation doesn't use managed
cat src/simulation/m_start_up.fpp | grep -A 5 -B 5 "allocate (airfoil_grid"

Repository: MFlowCode/MFC

Length of output: 521


🏁 Script executed:

# Check if there's any device-aware allocate macro being used
rg "@:ALLOCATE" src/simulation/m_start_up.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Let's look at the actual default clause in ACC_PARALLEL_LOOP macro to confirm it uses default(present)
grep -A 30 "def ACC_PARALLEL_LOOP" src/common/include/acc_macros.fpp | grep -i "default"

Repository: MFlowCode/MFC

Length of output: 277


🏁 Script executed:

# Verify which variables are in copyin for the calling loop
grep -B 2 -A 10 "GPU_PARALLEL_LOOP.*copy=" src/simulation/m_compute_levelset.fpp | head -40

Repository: MFlowCode/MFC

Length of output: 1055


🏁 Script executed:

# Check if there are ANY enter data regions in the entire simulation directory
rg "GPU_ENTER_DATA|acc enter data" src/simulation/ | head -20

Repository: MFlowCode/MFC

Length of output: 1863


🏁 Script executed:

# Let's see the complete context around s_apply_levelset to understand the full data region
sed -n '29,88p' src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 2157


🏁 Script executed:

# Check if airfoil_grid variables are used anywhere else with GPU directives
rg "airfoil_grid" src/simulation/ -B 2 -A 2 | grep -i "gpu\|acc\|omp" | head -20

Repository: MFlowCode/MFC

Length of output: 1793


🏁 Script executed:

# Final verification: check if there's any GPU_DATA wrapper around the entire s_apply_levelset call
rg "GPU_DATA" src/simulation/m_compute_levelset.fpp

Repository: MFlowCode/MFC

Length of output: 39


🏁 Script executed:

# Check where s_apply_levelset is called from to see if there's a data region there
rg "s_apply_levelset" src/simulation/ -B 3 -A 3 | head -40

Repository: MFlowCode/MFC

Length of output: 2003


Add airfoil_grid_u, airfoil_grid_l, and Np to GPU data movement clauses.

These arrays and scalar are read inside s_airfoil_levelset and s_3D_airfoil_levelset, but the GPU_PARALLEL_LOOP in s_apply_levelset only includes patch_ib in its copyin clause. With default(present), this causes a runtime fault on GPU. Add them to the copyin clause or move them into a persistent $:GPU_ENTER_DATA(...) region at startup in m_start_up.fpp.

🤖 Prompt for AI Agents
In `@src/simulation/m_compute_levelset.fpp` around lines 121 - 300, The GPU kernel
routines s_airfoil_levelset and s_3D_airfoil_levelset read the arrays
airfoil_grid_u, airfoil_grid_l and the scalar Np but the GPU_PARALLEL_LOOP in
s_apply_levelset only copyin patch_ib (with default(present)), causing
missing-device-data faults; fix by adding airfoil_grid_u, airfoil_grid_l and Np
to the GPU data movement (e.g., include them in the copyin clause of the
GPU_PARALLEL_LOOP that calls s_apply_levelset) or alternatively ensure they are
moved to the device at startup via a persistent $:GPU_ENTER_DATA(...) region
(for example in m_start_up.fpp) so the kernels see these symbols on the device.

Comment on lines 430 to 527
subroutine s_cuboid_levelset(gp)
$:GPU_ROUTINE(parallelism='[seq]')

type(ghost_point), intent(inout) :: gp

real(wp) :: Right, Left, Bottom, Top, Front, Back
real(wp) :: min_dist
real(wp) :: side_dists(6)

real(wp), dimension(3) :: center
real(wp) :: length_x, length_y, length_z
real(wp), dimension(1:3) :: xyz_local, dist_vec !< x and y coordinates in local IB frame
real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation

integer :: i, j, k !< Loop index variables
integer :: ib_patch_id !< patch ID

ib_patch_id = gp%ib_patch_id
i = gp%loc(1)
j = gp%loc(2)
k = gp%loc(3)

length_x = patch_ib(ib_patch_id)%length_x
length_y = patch_ib(ib_patch_id)%length_y
length_z = patch_ib(ib_patch_id)%length_z

center(1) = patch_ib(ib_patch_id)%x_centroid
center(2) = patch_ib(ib_patch_id)%y_centroid
center(3) = patch_ib(ib_patch_id)%z_centroid

inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)

Right = length_x/2
Left = -length_x/2
Top = length_y/2
Bottom = -length_y/2
Front = length_z/2
Back = -length_z/2

xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinate

side_dists(1) = Left - xyz_local(1)
side_dists(2) = xyz_local(1) - Right
side_dists(3) = Bottom - xyz_local(2)
side_dists(4) = xyz_local(2) - Top
side_dists(5) = Back - xyz_local(3)
side_dists(6) = xyz_local(3) - Front
min_dist = minval(abs(side_dists))

! TODO :: The way that this is written, it looks like we will
! trigger at the first size that is close to the minimum distance,
! meaning corners where side_dists are the same will
! trigger on what may not actually be the minimum,
! leading to undesired behavior. This should be resolved
! and this code should be cleaned up. It also means that
! rotating the box 90 degrees will cause tests to fail.
dist_vec = 0._wp
if (f_approx_equal(min_dist, abs(side_dists(1)))) then
gp%levelset = side_dists(1)
if (.not. f_approx_equal(side_dists(1), 0._wp)) then
dist_vec(1) = side_dists(1)/abs(side_dists(1))
end if

else if (f_approx_equal(min_dist, abs(side_dists(2)))) then
gp%levelset = side_dists(2)
if (.not. f_approx_equal(side_dists(2), 0._wp)) then
dist_vec(1) = -side_dists(2)/abs(side_dists(2))
end if

else if (f_approx_equal(min_dist, abs(side_dists(3)))) then
gp%levelset = side_dists(3)
if (.not. f_approx_equal(side_dists(3), 0._wp)) then
dist_vec(2) = side_dists(3)/abs(side_dists(3))
end if

else if (f_approx_equal(min_dist, abs(side_dists(4)))) then
gp%levelset = side_dists(4)
if (.not. f_approx_equal(side_dists(4), 0._wp)) then
dist_vec(2) = -side_dists(4)/abs(side_dists(4))
end if

else if (f_approx_equal(min_dist, abs(side_dists(5)))) then
gp%levelset = side_dists(5)
if (.not. f_approx_equal(side_dists(5), 0._wp)) then
dist_vec(3) = side_dists(5)/abs(side_dists(5))
end if

else if (f_approx_equal(min_dist, abs(side_dists(6)))) then
gp%levelset = side_dists(6)
if (.not. f_approx_equal(side_dists(6), 0._wp)) then
dist_vec(3) = -side_dists(6)/abs(side_dists(6))
end if
end if
gp%levelset_norm = matmul(rotation, dist_vec)

end subroutine s_cuboid_levelset
Copy link
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

The TODO indicates a known correctness risk in cuboid levelset.

The comment notes corner ambiguity and rotation-dependent failures. Please resolve this before release, or gate it to avoid incorrect IB geometry results.

If you want, I can propose an oriented-box signed-distance implementation to replace the TODO block.

🤖 Prompt for AI Agents
In `@src/simulation/m_compute_levelset.fpp` around lines 430 - 527, The cuboid
levelset has ambiguity at corners and under rotation; replace the cascade of
side_dists checks in subroutine s_cuboid_levelset with a robust oriented-box
signed-distance calculation: compute half_lengths =
(length_x/2,length_y/2,length_z/2), q = abs(xyz_local) - half_lengths, outside =
norm(max(q,0)), inside = min(max(q(1),max(q(2),q(3))), 0), and set gp%levelset =
outside + inside; compute the local normal as (if outside > 0)
normalize(max(q,0)) * sign(xyz_local) else unit vector along the component with
largest q (use sign(xyz_local) for direction), then rotate that normal back with
matmul(rotation, normal) into gp%levelset_norm; replace references to
side_dists/min_dist and ensure zero-division is guarded when normal magnitude is
~0.0.

Copy link
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 0

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_ib_patches.fpp (1)

899-943: ⚠️ Potential issue | 🔴 Critical

f_model_is_inside cannot be called from GPU device code—it is impure and lacks the required GPU routine marker.

This function (defined in src/common/m_model.fpp:488) contains call random_number(...) and is declared impure, making it incompatible with GPU device execution. Either:

  1. Refactor the function to be pure and mark it with $:GPU_ROUTINE(function_name='f_model_is_inside', parallelism='[seq]'), or
  2. Move the model-inside check outside the GPU loop and pre-compute the results before the parallel region, or
  3. Replace it with a GPU-compatible version without random sampling.

Per coding guidelines, GPU-callable helpers must be explicitly marked with the GPU routine macro, and impure operations (random number generation) cannot execute in device code.

@codeant-ai
Copy link
Contributor

codeant-ai bot commented Feb 5, 2026

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 ·
Reddit ·
LinkedIn

@codeant-ai codeant-ai bot added size:XXL This PR changes 1000+ lines, ignoring generated files and removed size:XXL This PR changes 1000+ lines, ignoring generated files labels Feb 5, 2026
Copy link
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 2

🤖 Fix all issues with AI agents
In `@src/simulation/m_compute_levelset.fpp`:
- Around line 64-83: The GPU loop over gps must privatize loop-local variables
and include missing read-only data: add patch_id and patch_geometry to the
private clause of the $:GPU_PARALLEL_LOOP for the 2D loop and add
airfoil_grid_u, airfoil_grid_l, and Np to the copyin clause (similar to the 3D
loop) so s_airfoil_levelset(gps(i)) has the needed data; ensure gps and patch_ib
remain in the appropriate copy/copyin clauses and verify the private/copyin
changes reference the same symbols (patch_id, patch_geometry, airfoil_grid_u,
airfoil_grid_l, Np, gps, patch_ib).

In `@src/simulation/m_ib_patches.fpp`:
- Line 913: Fix the typo in the comment for the f_model_is_inside note: change
"requries" to "requires" in the comment that references the non-GPU-compatible
random_number() subroutine and the TODO about commented macro calls; ensure the
comment reads "...requires the non-GPU-compatible random_number() subroutine..."
so spell-checkers no longer flag it while preserving the intended note about
f_model_is_inside and random_number().
🧹 Nitpick comments (3)
src/simulation/m_compute_levelset.fpp (2)

561-562: Use real literals for consistency with real(wp) arrays.

gp%levelset_norm is a real(wp) array, but the assignment uses integer literals. While Fortran will implicitly convert, using explicit real literals (1._wp, 0._wp) is more consistent with the codebase style.

♻️ Suggested fix
         if (f_approx_equal(dist, 0._wp)) then
-            gp%levelset_norm = (/1, 0, 0/)
+            gp%levelset_norm = (/1._wp, 0._wp, 0._wp/)
         else

169-172: Use real literals for dist_vec(3) assignments.

Multiple assignments like dist_vec(3) = 0 use integer literals. For consistency with other code in this file that uses 0._wp, consider using real literals throughout.

Similar occurrences at lines 171, 177, 191, 246, 260, 266, 280.

src/simulation/m_ib_patches.fpp (1)

117-243: Consider adding bounds check and idempotency guard for s_instantiate_STL_models.

  1. If called multiple times, the subroutine will reallocate models(patch_id)%model without deallocating previous allocations, causing memory leaks.
  2. Local arrays boundary_v and interpolated_boundary_v are allocated conditionally but the deallocation path is not shown within this subroutine.
🛡️ Suggested guard for re-entry
     subroutine s_instantiate_STL_models()
+        ! Guard against multiple calls
+        if (allocated(models)) then
+            return
+        end if
+
+        @:ALLOCATE(models(1:num_ibs))
+
         ! Variables for IBM+STL

Comment on lines +64 to +83
$:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib]')
do i = 1, num_gps

patch_id = gps(i)%ib_patch_id
patch_geometry = patch_ib(patch_id)%geometry

if (patch_geometry == 2) then
call s_circle_levelset(gps(i))
elseif (patch_geometry == 3) then
call s_rectangle_levelset(gps(i))
elseif (patch_geometry == 4) then
call s_airfoil_levelset(gps(i))
! STL+IBM patch
elseif (patch_geometry == 5) then
call s_model_levelset(gps(i))
elseif (patch_geometry == 6) then
call s_ellipse_levelset(gps(i))
end if
end do
$:END_GPU_PARALLEL_LOOP()
Copy link
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Privatize loop-local variables and add missing data for 2D GPU loop.

  1. patch_id and patch_geometry need to be in the private clause (same issue as 3D loop).
  2. The 2D loop calls s_airfoil_levelset (line 75) which uses airfoil_grid_u, airfoil_grid_l, and Np, but these are not in the copyin clause unlike the 3D loop at line 39.
🛠️ Suggested fix
-            $:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib]')
+            $:GPU_PARALLEL_LOOP(private='[i,patch_id,patch_geometry]', copy='[gps]', copyin='[patch_ib,airfoil_grid_u,airfoil_grid_l,Np]')
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
$:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib]')
do i = 1, num_gps
patch_id = gps(i)%ib_patch_id
patch_geometry = patch_ib(patch_id)%geometry
if (patch_geometry == 2) then
call s_circle_levelset(gps(i))
elseif (patch_geometry == 3) then
call s_rectangle_levelset(gps(i))
elseif (patch_geometry == 4) then
call s_airfoil_levelset(gps(i))
! STL+IBM patch
elseif (patch_geometry == 5) then
call s_model_levelset(gps(i))
elseif (patch_geometry == 6) then
call s_ellipse_levelset(gps(i))
end if
end do
$:END_GPU_PARALLEL_LOOP()
$:GPU_PARALLEL_LOOP(private='[i,patch_id,patch_geometry]', copy='[gps]', copyin='[patch_ib,airfoil_grid_u,airfoil_grid_l,Np]')
do i = 1, num_gps
patch_id = gps(i)%ib_patch_id
patch_geometry = patch_ib(patch_id)%geometry
if (patch_geometry == 2) then
call s_circle_levelset(gps(i))
elseif (patch_geometry == 3) then
call s_rectangle_levelset(gps(i))
elseif (patch_geometry == 4) then
call s_airfoil_levelset(gps(i))
! STL+IBM patch
elseif (patch_geometry == 5) then
call s_model_levelset(gps(i))
elseif (patch_geometry == 6) then
call s_ellipse_levelset(gps(i))
end if
end do
$:END_GPU_PARALLEL_LOOP()
🤖 Prompt for AI Agents
In `@src/simulation/m_compute_levelset.fpp` around lines 64 - 83, The GPU loop
over gps must privatize loop-local variables and include missing read-only data:
add patch_id and patch_geometry to the private clause of the $:GPU_PARALLEL_LOOP
for the 2D loop and add airfoil_grid_u, airfoil_grid_l, and Np to the copyin
clause (similar to the 3D loop) so s_airfoil_levelset(gps(i)) has the needed
data; ensure gps and patch_ib remain in the appropriate copy/copyin clauses and
verify the private/copyin changes reference the same symbols (patch_id,
patch_geometry, airfoil_grid_u, airfoil_grid_l, Np, gps, patch_ib).

Comment on lines +183 to +202
call f_check_interpolation_3D(model, (/dx, dy, dz/), interpolate)
else
call f_check_interpolation_2D(boundary_v, boundary_edge_count, (/minval(dx), minval(dy), 0._wp/), interpolate)
end if

! Show the number of edges and boundary edges in 2D STL models
if (proc_rank == 0 .and. p == 0) then
print *, ' * Number of 2D model boundary edges:', boundary_edge_count
end if

! Interpolate the STL model along the edges (2D) and on triangle facets (3D)
if (interpolate) then
if (proc_rank == 0) then
print *, ' * Interpolating STL vertices.'
end if

if (p > 0) then
call f_interpolate_3D(model, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
else
call f_interpolate_2D(boundary_v, boundary_edge_count, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: The STL interpolation and interpolation-check routines are called with array constructors built from the full dx, dy, and dz arrays instead of scalar characteristic spacings, which likely violates their expected interface (a 3-element spacing vector) and can lead to incorrect interpolation behavior or subtle shape mismatches. [possible bug]

Severity Level: Major ⚠️
- ❌ STL vertex interpolation may be incorrect.
- ⚠️ Geometry approximation density may be wrong.
- ⚠️ Affects 3D STL-imported immersed boundaries.
Suggested change
call f_check_interpolation_3D(model, (/dx, dy, dz/), interpolate)
else
call f_check_interpolation_2D(boundary_v, boundary_edge_count, (/minval(dx), minval(dy), 0._wp/), interpolate)
end if
! Show the number of edges and boundary edges in 2D STL models
if (proc_rank == 0 .and. p == 0) then
print *, ' * Number of 2D model boundary edges:', boundary_edge_count
end if
! Interpolate the STL model along the edges (2D) and on triangle facets (3D)
if (interpolate) then
if (proc_rank == 0) then
print *, ' * Interpolating STL vertices.'
end if
if (p > 0) then
call f_interpolate_3D(model, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
else
call f_interpolate_2D(boundary_v, boundary_edge_count, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
! Pass characteristic grid spacings in each direction
call f_check_interpolation_3D(model, (/minval(dx), minval(dy), minval(dz)/), interpolate)
else
call f_check_interpolation_2D(boundary_v, boundary_edge_count, (/minval(dx), minval(dy), 0._wp/), interpolate)
end if
! Show the number of edges and boundary edges in 2D STL models
if (proc_rank == 0 .and. p == 0) then
print *, ' * Number of 2D model boundary edges:', boundary_edge_count
end if
! Interpolate the STL model along the edges (2D) and on triangle facets (3D)
if (interpolate) then
if (proc_rank == 0) then
print *, ' * Interpolating STL vertices.'
end if
if (p > 0) then
call f_interpolate_3D(model, (/minval(dx), minval(dy), minval(dz)/), interpolated_boundary_v, total_vertices)
else
call f_interpolate_2D(boundary_v, boundary_edge_count, (/minval(dx), minval(dy), 0._wp/), interpolated_boundary_v, total_vertices)
Steps of Reproduction ✅
1. Open src/simulation/m_ib_patches.fpp and locate subroutine s_instantiate_STL_models
(definition begins at line ~117). The interpolation decision is performed at the check
beginning on line 181.

2. Run a 3D simulation (p > 0) that includes an STL-based immersed boundary
(patch_ib(...)%geometry == 5 or 12). The code path in s_instantiate_STL_models will reach
the interpolation decision at lines 181-186 and then call f_interpolate_3D at lines
~199-202.

3. Observe that the current code passes array constructors (/dx, dy, dz/) where dx,dy,dz
are full per-cell arrays instead of a 3-element spacing vector (e.g. characteristic
spacing per direction). If f_check_interpolation_3D or f_interpolate_3D expect a 3-element
vector (spacing in x,y,z) this mismatch can produce incorrect interpolation decisions or
incorrect spacing used when generating interpolated vertices.

4. Reproduce by running an STL case where interpolation is required (interpolate .true.);
compare interpolated vertex counts or geometry against a run where min-spacings are
explicitly supplied — differences in interpolation density or mis-sized geometry indicate
the issue.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_ib_patches.fpp
**Line:** 183:202
**Comment:**
	*Possible Bug: The STL interpolation and interpolation-check routines are called with array constructors built from the full `dx`, `dy`, and `dz` arrays instead of scalar characteristic spacings, which likely violates their expected interface (a 3-element spacing vector) and can lead to incorrect interpolation behavior or subtle shape mismatches.

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.

xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset

if (xyz_local(2) >= center(2)) then
Copy link
Contributor

@codeant-ai codeant-ai bot Feb 5, 2026

Choose a reason for hiding this comment

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

Suggestion: In the 3D airfoil levelset routine, the choice between the upper and lower airfoil surfaces incorrectly compares a local y-coordinate (xyz_local(2)) with a global centroid value (center(2)), so for rotated patches or non-zero centroids, many points will be classified against the wrong surface, yielding wrong distances and normals. [logic error]

Severity Level: Major ⚠️
- ❌ Levelset for 3D airfoil patches (s_3d_airfoil_levelset).
- ⚠️ Normals for rotated/offset airfoils incorrect.
- ⚠️ Ghost-point classification wrong for many points.
Suggested change
if (xyz_local(2) >= center(2)) then
if (xyz_local(2) >= 0._wp) then
Steps of Reproduction ✅
1. Ensure a 3D airfoil IBM patch is present so s_apply_levelset dispatches to
s_3d_airfoil_levelset. See s_apply_levelset dispatch in
src/simulation/m_compute_levelset.fpp: lines 40-52 where patch_geometry == 11 calls
s_3d_airfoil_levelset (call at line 51 in the PR diff).

2. Run the simulation path that invokes s_apply_levelset with ghost points for that
patch_id (s_apply_levelset entry at line 29 in the same file). The loop at lines 40-57
iterates ghost points and calls s_3d_airfoil_levelset(gps(i)).

3. Inside s_3d_airfoil_levelset (subroutine header at line 204), the code computes a local
coordinate xyz_local (line 238) and then decides which airfoil surface (upper/lower) to
search with the test at line 242: if (xyz_local(2) >= center(2)) then.

4. For any patch with a non-zero centroid or a rotation matrix (patch_ib(...)%x_centroid,
%y_centroid, %rotation_matrix set), comparing xyz_local(2) to center(2) is incorrect:
center(2) is the global centroid (set earlier at lines 227-229), not the local coordinate
origin. This causes points that should be classified by sign of the local y (xyz_local(2)
relative to 0) to be routed to the wrong surface. Reproduce by creating an airfoil patch
with non-zero y_centroid or non-identity rotation and observing incorrect gp%levelset /
gp%levelset_norm values emitted by s_3d_airfoil_levelset at lines 284-301 and 295-301.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_compute_levelset.fpp
**Line:** 242:242
**Comment:**
	*Logic Error: In the 3D airfoil levelset routine, the choice between the upper and lower airfoil surfaces incorrectly compares a local y-coordinate (`xyz_local(2)`) with a global centroid value (`center(2)`), so for rotated patches or non-zero centroids, many points will be classified against the wrong surface, yielding wrong distances and normals.

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.

Comment on lines 284 to 288
dist_side = min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l)))

if (dist_side < dist_surf) then
gp%levelset = dist_side
if (f_approx_equal(dist_side, abs(z_cc(l) - z_min))) then
Copy link
Contributor

@codeant-ai codeant-ai bot Feb 5, 2026

Choose a reason for hiding this comment

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

Suggestion: The distance to the end caps of the 3D airfoil is computed using global z_cc and z_min/z_max, which ignores the patch rotation and local coordinate system, so for rotated or offset airfoils the extrusion distance and normals along the spanwise direction are wrong. [logic error]

Severity Level: Major ⚠️
- ❌ End-cap distances for 3D airfoils (s_3d_airfoil_levelset).
- ⚠️ Spanwise normals for rotated patches incorrect.
- ⚠️ Levelset sign/magnitude errors near end faces.
Suggested change
dist_side = min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l)))
if (dist_side < dist_surf) then
gp%levelset = dist_side
if (f_approx_equal(dist_side, abs(z_cc(l) - z_min))) then
dist_side = min(abs(xyz_local(3) + 0.5_wp*lz), &
abs(0.5_wp*lz - xyz_local(3)))
if (dist_side < dist_surf) then
gp%levelset = dist_side
if (f_approx_equal(dist_side, abs(xyz_local(3) + 0.5_wp*lz))) then
Steps of Reproduction ✅
1. Trigger s_apply_levelset with a 3D airfoil patch (patch_geometry == 11).
s_apply_levelset invokes s_3d_airfoil_levelset in src/simulation/m_compute_levelset.fpp
(call at line 51 in the PR diff).

2. In s_3d_airfoil_levelset (subroutine at line 204), the code computes xyz_local by
translating and rotating the global coordinates into the IB-local frame (lines 238-240).

3. Later the code computes dist_side using global z coordinates at line 284: dist_side =
min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l))). This ignores the local rotated frame and
centroid offset computed earlier (lines 229-236 and 238-240).

4. For any patch with non-zero rotation or centroid offset in z
(patch_ib(...)%rotation_matrix, %z_centroid, %length_z non-zero), dist_side will be wrong.
Reproduce by creating a rotated/offset 3D airfoil patch and inspecting gp%levelset and
gp%levelset_norm after s_3d_airfoil_levelset returns — the end-cap distances and normals
(assigned on lines 286-294) will be inconsistent with the rotated local geometry.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_compute_levelset.fpp
**Line:** 284:288
**Comment:**
	*Logic Error: The distance to the end caps of the 3D airfoil is computed using global `z_cc` and `z_min/z_max`, which ignores the patch rotation and local coordinate system, so for rotated or offset airfoils the extrusion distance and normals along the spanwise direction are wrong.

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.

xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
xy_local = matmul(inverse_rotation, xy_local)

if ((xy_local(1) > bottom_left(1) .and. xy_local(1) < top_right(1)) .or. &
Copy link
Contributor

@codeant-ai codeant-ai bot Feb 5, 2026

Choose a reason for hiding this comment

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

Suggestion: In the rectangle levelset routine, the inside-region test uses a logical OR between the x and y ranges, so points that are outside in one direction but inside in the other are incorrectly treated as inside the rectangle, leading to wrong levelset values and normals in those regions. [logic error]

Severity Level: Major ⚠️
- ❌ Levelset for rectangle patches (s_rectangle_levelset).
- ⚠️ Normals near rectangle edges incorrect.
- ⚠️ Misclassification of outside points as inside.
Suggested change
if ((xy_local(1) > bottom_left(1) .and. xy_local(1) < top_right(1)) .or. &
if ((xy_local(1) > bottom_left(1) .and. xy_local(1) < top_right(1)) .and. &
Steps of Reproduction ✅
1. Ensure a 2D rectangle IBM patch is present so s_apply_levelset dispatches to
s_rectangle_levelset. In src/simulation/m_compute_levelset.fpp, s_apply_levelset calls
s_rectangle_levelset when patch_geometry == 3 (call location around lines 70-76 in the PR
diff).

2. Execute s_rectangle_levelset (subroutine header at line 306). The code computes
xy_local (lines 342-343) and then enters the bounds test at lines 345-346 using a logical
OR.

3. Pick a ghost point whose local x is inside the rectangle but local y is outside (or
vice-versa). Because the condition uses OR, the code treats that point as "inside" and
proceeds to compute side distances and normals (lines 348-376), producing an incorrect
gp%levelset and gp%levelset_norm for that point.

4. Reproduce by constructing a rectangle patch and sampling points that are
inside-x/outside-y (or inside-y/outside-x) and observing that gp%levelset is set (should
be outside) — evidence that the OR is incorrect and the intended test is an AND.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_compute_levelset.fpp
**Line:** 345:345
**Comment:**
	*Logic Error: In the rectangle levelset routine, the inside-region test uses a logical OR between the x and y ranges, so points that are outside in one direction but inside in the other are incorrectly treated as inside the rectangle, leading to wrong levelset values and normals in those regions.

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
Copy link
Contributor

codeant-ai bot commented Feb 5, 2026

CodeAnt AI Incremental review completed.

Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot left a comment

Choose a reason for hiding this comment

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

1 issue found across 3 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_compute_levelset.fpp">

<violation number="1" location="src/simulation/m_compute_levelset.fpp:39">
P2: The 3D GPU loop now always copyin’s airfoil_grid_u/l, but those arrays are only allocated when an airfoil patch exists. In 3D runs without airfoils this can attempt a GPU copy of unallocated arrays and fail at runtime. Guard the copyin with allocation checks or ensure these arrays are allocated before entering the loop.</violation>
</file>

Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.

! 3D Patch Geometries
if (p > 0) then

$:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib,airfoil_grid_u,airfoil_grid_l]')
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Feb 5, 2026

Choose a reason for hiding this comment

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

P2: The 3D GPU loop now always copyin’s airfoil_grid_u/l, but those arrays are only allocated when an airfoil patch exists. In 3D runs without airfoils this can attempt a GPU copy of unallocated arrays and fail at runtime. Guard the copyin with allocation checks or ensure these arrays are allocated before entering the loop.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_compute_levelset.fpp, line 39:

<comment>The 3D GPU loop now always copyin’s airfoil_grid_u/l, but those arrays are only allocated when an airfoil patch exists. In 3D runs without airfoils this can attempt a GPU copy of unallocated arrays and fail at runtime. Guard the copyin with allocation checks or ensure these arrays are allocated before entering the loop.</comment>

<file context>
@@ -36,7 +36,7 @@ contains
         if (p > 0) then
 
-            $:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib]')
+            $:GPU_PARALLEL_LOOP(private='[i]', copy='[gps]', copyin='[patch_ib,airfoil_grid_u,airfoil_grid_l]')
             do i = 1, num_gps
 
</file context>
Fix with Cubic

Daniel J Vickers and others added 3 commits February 4, 2026 23:23
@codecov
Copy link

codecov bot commented Feb 5, 2026

Codecov Report

❌ Patch coverage is 58.50622% with 200 lines in your changes missing coverage. Please review.
✅ Project coverage is 44.09%. Comparing base (1d3286a) to head (05c069e).

Files with missing lines Patch % Lines
src/simulation/m_compute_levelset.fpp 53.77% 136 Missing and 5 partials ⚠️
src/simulation/m_ib_patches.fpp 53.94% 31 Missing and 4 partials ⚠️
src/simulation/m_ibm.fpp 43.47% 13 Missing ⚠️
src/common/m_model.fpp 0.00% 4 Missing ⚠️
src/pre_process/m_start_up.fpp 25.00% 3 Missing ⚠️
src/simulation/m_data_output.fpp 90.32% 3 Missing ⚠️
src/simulation/m_start_up.fpp 50.00% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1123      +/-   ##
==========================================
+ Coverage   44.03%   44.09%   +0.06%     
==========================================
  Files          70       70              
  Lines       20649    20447     -202     
  Branches     2053     1973      -80     
==========================================
- Hits         9093     9017      -76     
+ Misses      10368    10295      -73     
+ Partials     1188     1135      -53     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Review effort 4/5 size:XXL This PR changes 1000+ lines, ignoring generated files

Development

Successfully merging this pull request may close these issues.

Remove Levelset and Levelset Norm Arrays

1 participant