Add rotating mibms #1014

danieljvickers · 2025-10-07T16:34:04Z

User description

Description

This PR improves on the previous MIBM PR in that it adds higher-order time stepping of the IB centroids/angles (supporting RK1-3, as opposed to only RK1) and it adds rotation of IBs.

Type of change

Please delete options that are not relevant.

New feature (non-breaking change which adds functionality)

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?

Please describe the tests that you ran to verify your changes.
Provide instructions so we can reproduce.
Please also list any relevant details for your test configuration

I generated a tumbling rectangle which rotates, moves in a stationary fluid, starts with a rotated initial angle, and performs RK3 integration. The case file is as given:

import json
import math

Mu = 1.84e-05
gam_a = 1.4

# Configuring case dictionary
print(
    json.dumps(
        {
            # Logistics
            "run_time_info": "T",
            # Computational Domain Parameters
            # For these computations, the cylinder is placed at the (0,0,0)
            # domain origin.
            # axial direction
            "x_domain%beg": 0.0e00,
            "x_domain%end": 6.0e-03,
            # r direction
            "y_domain%beg": 0.0e00,
            "y_domain%end": 6.0e-03,
            "cyl_coord": "F",
            "m": 75,
            "n": 75,
            "p": 0,
            "dt": 6.0e-6,
            "t_step_start": 0,
            "t_step_stop": 10000,
            "t_step_save": 100,
            # Simulation Algorithm Parameters
            # Only one patches are necessary, the air tube
            "num_patches": 1,
            # Use the 5 equation model
            "model_eqns": 2,
            "alt_soundspeed": "F",
            # One fluids: air
            "num_fluids": 1,
            # time step
            "mpp_lim": "F",
            # Correct errors when computing speed of sound
            "mixture_err": "T",
            # Use TVD RK3 for time marching
            "time_stepper": 3,
            # Use WENO5
            "weno_order": 5,
            "weno_eps": 1.0e-16,
            "weno_Re_flux": "T",
            "weno_avg": "T",
            "avg_state": 2,
            "mapped_weno": "T",
            "null_weights": "F",
            "mp_weno": "T",
            "riemann_solver": 2,
            "wave_speeds": 1,
            # We use ghost-cell
            "bc_x%beg": -3,
            "bc_x%end": -3,
            "bc_y%beg": -3,
            "bc_y%end": -3,
            # Set IB to True and add 1 patch
            "ib": "T",
            "num_ibs": 1,
            "viscous": "T",
            # Formatted Database Files Structure Parameters
            "format": 1,
            "precision": 2,
            "prim_vars_wrt": "T",
            "E_wrt": "T",
            "parallel_io": "T",
            # Patch: Constant Tube filled with air
            # Specify the cylindrical air tube grid geometry
            "patch_icpp(1)%geometry": 3,
            "patch_icpp(1)%x_centroid": 3.0e-03,
            # Uniform medium density, centroid is at the center of the domain
            "patch_icpp(1)%y_centroid": 3.0e-03,
            "patch_icpp(1)%length_x": 6.0e-03,
            "patch_icpp(1)%length_y": 6.0e-03,
            # Specify the patch primitive variables
            "patch_icpp(1)%vel(1)": 0.00e00,
            "patch_icpp(1)%vel(2)": 0.0e00,
            "patch_icpp(1)%pres": 1.0e00,
            "patch_icpp(1)%alpha_rho(1)": 1.0e00,
            "patch_icpp(1)%alpha(1)": 10.0e00,
            # Patch: Cylinder Immersed Boundary
            "patch_ib(1)%geometry": 3,
            "patch_ib(1)%x_centroid": 4.5e-03,
            "patch_ib(1)%y_centroid": 3.0e-03,
            "patch_ib(1)%length_x": 0.6e-03,
            "patch_ib(1)%length_y": 0.6e-03,
            "patch_ib(1)%slip": "F",
            "patch_ib(1)%moving_ibm": 1,
            "patch_ib(1)%vel(1)": -0.05,
            "patch_ib(1)%angles(1)": 0.0, # x-axis rotation in radians
            "patch_ib(1)%angles(2)": 0.0, # y-axis rotation
            "patch_ib(1)%angles(3)": 0.78539816339, # z-axis rotation
            "patch_ib(1)%angular_vel(1)": 0.0, # x-axis rotational velcoity in radians per second
            "patch_ib(1)%angular_vel(2)": 0.0, # y-axis rotation
            "patch_ib(1)%angular_vel(3)": 100.0, # z-axis rotation
            # Fluids Physical Parameters
            "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50(Not 1.40)
            "fluid_pp(1)%pi_inf": 0,
            "fluid_pp(1)%Re(1)": 2500000,
        }
    )
)

Test Configuration:

What computers and compilers did you use to test this:

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

Description

Add rotation support for immersed boundaries (IBs)
Implement higher-order time stepping (RK1-3) for IB motion
Add angular velocity and rotation matrix calculations
Update ghost point velocity calculations with rotational effects

Diagram Walkthrough

flowchart LR
  A["IB Parameters"] --> B["Rotation Matrix"]
  B --> C["Ghost Point Velocity"]
  C --> D["Time Stepping"]
  D --> E["Updated IB Position/Angle"]
  E --> A

File Walkthrough

Relevant files

Enhancement

6 files

m_derived_types.fpp `Add rotation and time stepping fields`	+10/-0
m_ib_patches.fpp `Add rotation matrix calculation and coordinate transformation`	+60/-9
m_global_parameters.fpp `Initialize rotation matrices and angular velocity`	+10/-3
m_initial_condition.fpp `Add rotation matrix initialization call`	+5/-0
m_ibm.fpp `Update ghost point velocity with rotation`	+26/-28
m_time_steppers.fpp `Implement RK1-3 time stepping for IB motion`	+134/-0

Miscellaneous

2 files

m_start_up.fpp `Remove redundant MIB update call`	+0/-4
p_main.fpp `Add debug print for IB parameters`	+1/-0

Configuration changes

1 files

case_dicts.py `Add angular velocity parameters to case dictionaries`	+8/-4

…ifying rotated markers and velocity calculations

qodo-merge-pro · 2025-10-07T16:34:53Z

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪
🧪 No relevant tests
🔒 No security concerns identified
⚡ Recommended focus areas for review Duplicate Logic The MIBM update/integration logic for velocities, angles, centroids is copy-pasted across multiple RK stages with slight variations, increasing maintenance risk and potential inconsistencies. Consider refactoring into shared routines for RK1/2/3 stages and centralizing the Euler/step storage updates. if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(1)%vf) if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) if (ib) then ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms if (moving_immersed_boundary_flag) then do i = 1, num_ibs ! start by using euler's method naiively, but eventually incorporate more sophistocation if (patch_ib(i)%moving_ibm == 1) then do j = 1, 3 patch_ib(i)%vel(j) = patch_ib(i)%vel(j) + 0.0dt patch_ib(i)%angular_vel(j) = patch_ib(i)%angular_vel(j) + 0.0dt ! Update the angle of the IB patch_ib(i)%angles(j) = patch_ib(i)%angles(j) + patch_ib(i)%angular_vel(j)dt end do ! Update the position of the IB patch_ib(i)%x_centroid = patch_ib(i)%x_centroid + patch_ib(i)%vel(1)dt patch_ib(i)%y_centroid = patch_ib(i)%y_centroid + patch_ib(i)%vel(2)dt patch_ib(i)%z_centroid = patch_ib(i)%z_centroid + patch_ib(i)%vel(3)dt end if end do call s_update_mib(num_ibs, levelset, levelset_norm) end if if (qbmm .and. .not. polytropic) then call s_ibm_correct_state(q_cons_ts(1)%vf, q_prim_vf, pb_ts(1)%sf, mv_ts(1)%sf) else call s_ibm_correct_state(q_cons_ts(1)%vf, q_prim_vf) end if end if call nvtxEndRange call cpu_time(finish) wall_time = abs(finish - start) if (t_step >= 2) then wall_time_avg = (wall_time + (t_step - 2)wall_time_avg)/(t_step - 1) else wall_time_avg = 0._wp end if end subroutine s_1st_order_tvd_rk !> 2nd order TVD RK time-stepping algorithm !! @param t_step Current time-step impure subroutine s_2nd_order_tvd_rk(t_step, time_avg) #ifdef _CRAYFTN !DIR$ OPTIMIZE (-haggress) #endif integer, intent(in) :: t_step real(wp), intent(inout) :: time_avg integer :: i, j, k, l, q!< Generic loop iterator real(wp) :: start, finish integer :: dest ! Stage 1 of 2 call cpu_time(start) call nvtxStartRange("TIMESTEP") call s_compute_rhs(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, t_step, time_avg, 1) if (run_time_info) then if (igr) then call s_write_run_time_information(q_cons_ts(1)%vf, t_step) else call s_write_run_time_information(q_prim_vf, t_step) end if end if if (probe_wrt) then call s_time_step_cycling(t_step) end if if (cfl_dt) then if (mytime >= t_stop) return else if (t_step == t_step_stop) return end if if (bubbles_lagrange .and. .not. adap_dt) call s_update_lagrange_tdv_rk(stage=1) #if defined(__NVCOMPILER_GPU_UNIFIED_MEM) \|\| defined(FRONTIER_UNIFIED) $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(2)%vf(i)%sf(j, k, l) = & q_cons_ts(1)%vf(i)%sf(j, k, l) q_cons_ts(1)%vf(i)%sf(j, k, l) = & q_cons_ts(1)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l) end do end do end do end do dest = 1 ! Result in q_cons_ts(1)%vf #else $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(2)%vf(i)%sf(j, k, l) = & q_cons_ts(1)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l) end do end do end do end do dest = 2 ! Result in q_cons_ts(2)%vf #endif !Evolve pb and mv for non-polytropic qbmm if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode pb_ts(2)%sf(j, k, l, q, i) = & pb_ts(1)%sf(j, k, l, q, i) & + dtrhs_pb(j, k, l, q, i) end do end do end do end do end do end if if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode mv_ts(2)%sf(j, k, l, q, i) = & mv_ts(1)%sf(j, k, l, q, i) & + dtrhs_mv(j, k, l, q, i) end do end do end do end do end do end if if (bodyForces) call s_apply_bodyforces(q_cons_ts(dest)%vf, q_prim_vf, rhs_vf, dt) if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(dest)%vf) if (model_eqns == 3 .and. (.not. relax)) then call s_pressure_relaxation_procedure(q_cons_ts(dest)%vf) end if if (adv_n) call s_comp_alpha_from_n(q_cons_ts(dest)%vf) if (ib) then ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms if (moving_immersed_boundary_flag) then do i = 1, num_ibs ! start by using euler's method naiively, but eventually incorporate more sophistocation if (patch_ib(i)%moving_ibm == 1) then patch_ib(i)%step_vel = patch_ib(i)%vel patch_ib(i)%step_angular_vel = patch_ib(i)%angular_vel patch_ib(i)%step_angles = patch_ib(i)%angles do j = 1, 3 patch_ib(i)%vel(j) = patch_ib(i)%step_vel(j) + 0.0dt patch_ib(i)%angular_vel(j) = patch_ib(i)%step_angular_vel(j) + 0.0dt ! Update the angle of the IB patch_ib(i)%angles(j) = patch_ib(i)%step_angles(j) + patch_ib(i)%angular_vel(j)dt end do patch_ib(i)%step_x_centroid = patch_ib(i)%x_centroid patch_ib(i)%step_y_centroid = patch_ib(i)%y_centroid patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid patch_ib(i)%x_centroid = patch_ib(i)%step_x_centroid + patch_ib(i)%vel(1)dt patch_ib(i)%y_centroid = patch_ib(i)%step_y_centroid + patch_ib(i)%vel(2)dt patch_ib(i)%z_centroid = patch_ib(i)%step_z_centroid + patch_ib(i)%vel(3)dt end if end do call s_update_mib(num_ibs, levelset, levelset_norm) end if if (qbmm .and. .not. polytropic) then call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf, pb_ts(2)%sf, mv_ts(2)%sf) else call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf) end if end if ! Stage 2 of 2 call s_compute_rhs(q_cons_ts(dest)%vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_ts(2)%sf, rhs_pb, mv_ts(2)%sf, rhs_mv, t_step, time_avg, 2) if (bubbles_lagrange .and. .not. adap_dt) call s_update_lagrange_tdv_rk(stage=2) #if defined(__NVCOMPILER_GPU_UNIFIED_MEM) \|\| defined(FRONTIER_UNIFIED) $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(1)%vf(i)%sf(j, k, l) = & (q_cons_ts(2)%vf(i)%sf(j, k, l) & + q_cons_ts(1)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l))/2._wp end do end do end do end do dest = 1 ! Result in q_cons_ts(1)%vf #else $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(1)%vf(i)%sf(j, k, l) = & (q_cons_ts(1)%vf(i)%sf(j, k, l) & + q_cons_ts(2)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l))/2._wp end do end do end do end do dest = 1 ! Result in q_cons_ts(1)%vf #endif if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode pb_ts(1)%sf(j, k, l, q, i) = & (pb_ts(1)%sf(j, k, l, q, i) & + pb_ts(2)%sf(j, k, l, q, i) & + dtrhs_pb(j, k, l, q, i))/2._wp end do end do end do end do end do end if if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode mv_ts(1)%sf(j, k, l, q, i) = & (mv_ts(1)%sf(j, k, l, q, i) & + mv_ts(2)%sf(j, k, l, q, i) & + dtrhs_mv(j, k, l, q, i))/2._wp end do end do end do end do end do end if if (bodyForces) call s_apply_bodyforces(q_cons_ts(dest)%vf, q_prim_vf, rhs_vf, 2._wpdt/3._wp) if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(dest)%vf) if (model_eqns == 3 .and. (.not. relax)) then call s_pressure_relaxation_procedure(q_cons_ts(dest)%vf) end if if (adv_n) call s_comp_alpha_from_n(q_cons_ts(dest)%vf) if (ib) then if (moving_immersed_boundary_flag) then do i = 1, num_ibs ! start by using euler's method naiively, but eventually incorporate more sophistocation if (patch_ib(i)%moving_ibm == 1) then do j = 1, 3 patch_ib(i)%vel(j) = (patch_ib(i)%vel(j) + patch_ib(i)%step_vel(j) + 0.0dt)/2._wp patch_ib(i)%angular_vel(j) = (patch_ib(i)%angular_vel(j) + patch_ib(i)%step_angular_vel(j) + 0.0dt)/2._wp patch_ib(i)%angles(j) = (patch_ib(i)%angles(j) + patch_ib(i)%step_angles(j) + patch_ib(i)%step_angular_vel(j)dt)/2._wp end do patch_ib(i)%x_centroid = (patch_ib(i)%x_centroid + patch_ib(i)%step_x_centroid + patch_ib(i)%step_vel(1)dt)/2._wp patch_ib(i)%y_centroid = (patch_ib(i)%y_centroid + patch_ib(i)%step_y_centroid + patch_ib(i)%step_vel(2)dt)/2._wp patch_ib(i)%z_centroid = (patch_ib(i)%z_centroid + patch_ib(i)%step_z_centroid + patch_ib(i)%step_vel(3)dt)/2._wp end if end do call s_update_mib(num_ibs, levelset, levelset_norm) end if if (qbmm .and. .not. polytropic) then call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf, pb_ts(1)%sf, mv_ts(1)%sf) else call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf) end if end if call nvtxEndRange call cpu_time(finish) wall_time = abs(finish - start) if (t_step >= 2) then wall_time_avg = (wall_time + (t_step - 2)wall_time_avg)/(t_step - 1) else wall_time_avg = 0._wp end if end subroutine s_2nd_order_tvd_rk !> 3rd order TVD RK time-stepping algorithm !! @param t_step Current time-step impure subroutine s_3rd_order_tvd_rk(t_step, time_avg) #ifdef _CRAYFTN !DIR$ OPTIMIZE (-haggress) #endif integer, intent(IN) :: t_step real(wp), intent(INOUT) :: time_avg integer :: i, j, k, l, q !< Generic loop iterator real(wp) :: start, finish integer :: dest ! Stage 1 of 3 if (.not. adap_dt) then call cpu_time(start) call nvtxStartRange("TIMESTEP") end if call s_compute_rhs(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_ts(1)%sf, rhs_pb, mv_ts(1)%sf, rhs_mv, t_step, time_avg, 1) if (run_time_info) then if (igr) then call s_write_run_time_information(q_cons_ts(1)%vf, t_step) else call s_write_run_time_information(q_prim_vf, t_step) end if end if if (probe_wrt) then call s_time_step_cycling(t_step) end if if (cfl_dt) then if (mytime >= t_stop) return else if (t_step == t_step_stop) return end if if (bubbles_lagrange .and. .not. adap_dt) call s_update_lagrange_tdv_rk(stage=1) #if defined(__NVCOMPILER_GPU_UNIFIED_MEM) \|\| defined(FRONTIER_UNIFIED) $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(2)%vf(i)%sf(j, k, l) = & q_cons_ts(1)%vf(i)%sf(j, k, l) q_cons_ts(1)%vf(i)%sf(j, k, l) = & q_cons_ts(1)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l) end do end do end do end do dest = 1 ! result in q_cons_ts(1)%vf #else $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(2)%vf(i)%sf(j, k, l) = & q_cons_ts(1)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l) end do end do end do end do dest = 2 ! result in q_cons_ts(2)%vf #endif !Evolve pb and mv for non-polytropic qbmm if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode pb_ts(2)%sf(j, k, l, q, i) = & pb_ts(1)%sf(j, k, l, q, i) & + dtrhs_pb(j, k, l, q, i) end do end do end do end do end do end if if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode mv_ts(2)%sf(j, k, l, q, i) = & mv_ts(1)%sf(j, k, l, q, i) & + dtrhs_mv(j, k, l, q, i) end do end do end do end do end do end if if (bodyForces) call s_apply_bodyforces(q_cons_ts(dest)%vf, q_prim_vf, rhs_vf, dt) if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(dest)%vf) if (model_eqns == 3 .and. (.not. relax)) then call s_pressure_relaxation_procedure(q_cons_ts(dest)%vf) end if if (adv_n) call s_comp_alpha_from_n(q_cons_ts(dest)%vf) if (ib) then if (moving_immersed_boundary_flag) then do i = 1, num_ibs if (patch_ib(i)%moving_ibm == 1) then patch_ib(i)%step_vel = patch_ib(i)%vel patch_ib(i)%step_angular_vel = patch_ib(i)%angular_vel patch_ib(i)%step_angles = patch_ib(i)%angles do j = 1, 3 patch_ib(i)%vel(j) = patch_ib(i)%step_vel(j) + 0.0dt patch_ib(i)%angular_vel(j) = patch_ib(i)%step_angular_vel(j) + 0.0dt ! Update the angle of the IB patch_ib(i)%angles(j) = patch_ib(i)%step_angles(j) + patch_ib(i)%angular_vel(j)dt end do patch_ib(i)%step_x_centroid = patch_ib(i)%x_centroid patch_ib(i)%step_y_centroid = patch_ib(i)%y_centroid patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid patch_ib(i)%x_centroid = patch_ib(i)%step_x_centroid + patch_ib(i)%vel(1)dt patch_ib(i)%y_centroid = patch_ib(i)%step_y_centroid + patch_ib(i)%vel(2)dt patch_ib(i)%z_centroid = patch_ib(i)%step_z_centroid + patch_ib(i)%vel(3)dt end if end do call s_update_mib(num_ibs, levelset, levelset_norm) end if if (qbmm .and. .not. polytropic) then call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf, pb_ts(2)%sf, mv_ts(2)%sf) else call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf) end if end if ! Stage 2 of 3 call s_compute_rhs(q_cons_ts(dest)%vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_ts(2)%sf, rhs_pb, mv_ts(2)%sf, rhs_mv, t_step, time_avg, 2) if (bubbles_lagrange .and. .not. adap_dt) call s_update_lagrange_tdv_rk(stage=2) #if defined(__NVCOMPILER_GPU_UNIFIED_MEM) \|\| defined(FRONTIER_UNIFIED) $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(1)%vf(i)%sf(j, k, l) = & (3._wpq_cons_ts(2)%vf(i)%sf(j, k, l) & + q_cons_ts(1)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l))/4._wp end do end do end do end do dest = 1 ! Result in q_cons_ts(1)%vf #else $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(2)%vf(i)%sf(j, k, l) = & (3._wpq_cons_ts(1)%vf(i)%sf(j, k, l) & + q_cons_ts(2)%vf(i)%sf(j, k, l) & + dtrhs_vf(i)%sf(j, k, l))/4._wp end do end do end do end do dest = 2 ! Result in q_cons_ts(2)%vf #endif if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode pb_ts(2)%sf(j, k, l, q, i) = & (3._wppb_ts(1)%sf(j, k, l, q, i) & + pb_ts(2)%sf(j, k, l, q, i) & + dtrhs_pb(j, k, l, q, i))/4._wp end do end do end do end do end do end if if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode mv_ts(2)%sf(j, k, l, q, i) = & (3._wpmv_ts(1)%sf(j, k, l, q, i) & + mv_ts(2)%sf(j, k, l, q, i) & + dtrhs_mv(j, k, l, q, i))/4._wp end do end do end do end do end do end if if (bodyForces) call s_apply_bodyforces(q_cons_ts(dest)%vf, q_prim_vf, rhs_vf, dt/4._wp) if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(dest)%vf) if (model_eqns == 3 .and. (.not. relax)) then call s_pressure_relaxation_procedure(q_cons_ts(dest)%vf) end if if (adv_n) call s_comp_alpha_from_n(q_cons_ts(dest)%vf) if (ib) then if (moving_immersed_boundary_flag) then do i = 1, num_ibs if (patch_ib(i)%moving_ibm == 1) then do j = 1, 3 patch_ib(i)%vel(j) = (patch_ib(i)%vel(j) + 3._wppatch_ib(i)%step_vel(j) + 0.0dt)/4._wp patch_ib(i)%angular_vel(j) = (patch_ib(i)%angular_vel(j) + 3._wppatch_ib(i)%step_angular_vel(j) + 0.0dt)/4._wp patch_ib(i)%angles(j) = (patch_ib(i)%angles(j) + 3._wppatch_ib(i)%step_angles(j) + patch_ib(i)%angular_vel(j)dt)/4._wp end do patch_ib(i)%x_centroid = (patch_ib(i)%x_centroid + 3._wppatch_ib(i)%step_x_centroid + patch_ib(i)%vel(1)dt)/4._wp patch_ib(i)%y_centroid = (patch_ib(i)%y_centroid + 3._wppatch_ib(i)%step_y_centroid + patch_ib(i)%vel(2)dt)/4._wp patch_ib(i)%z_centroid = (patch_ib(i)%z_centroid + 3._wppatch_ib(i)%step_z_centroid + patch_ib(i)%vel(3)dt)/4._wp end if end do call s_update_mib(num_ibs, levelset, levelset_norm) end if if (qbmm .and. .not. polytropic) then call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf, pb_ts(2)%sf, mv_ts(2)%sf) else call s_ibm_correct_state(q_cons_ts(dest)%vf, q_prim_vf) end if end if ! Stage 3 of 3 call s_compute_rhs(q_cons_ts(dest)%vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_ts(2)%sf, rhs_pb, mv_ts(2)%sf, rhs_mv, t_step, time_avg, 3) if (bubbles_lagrange .and. .not. adap_dt) call s_update_lagrange_tdv_rk(stage=3) #if defined(__NVCOMPILER_GPU_UNIFIED_MEM) \|\| defined(FRONTIER_UNIFIED) $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(1)%vf(i)%sf(j, k, l) = & (q_cons_ts(2)%vf(i)%sf(j, k, l) & + 2._wpq_cons_ts(1)%vf(i)%sf(j, k, l) & + 2._wpdtrhs_vf(i)%sf(j, k, l))/3._wp end do end do end do end do dest = 1 ! Result in q_cons_ts(1)%vf #else $:GPU_PARALLEL_LOOP(collapse=4) do i = 1, sys_size do l = 0, p do k = 0, n do j = 0, m q_cons_ts(1)%vf(i)%sf(j, k, l) = & (q_cons_ts(1)%vf(i)%sf(j, k, l) & + 2._wpq_cons_ts(2)%vf(i)%sf(j, k, l) & + 2._wpdtrhs_vf(i)%sf(j, k, l))/3._wp end do end do end do end do dest = 1 ! Result in q_cons_ts(2)%vf #endif if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode pb_ts(1)%sf(j, k, l, q, i) = & (pb_ts(1)%sf(j, k, l, q, i) & + 2._wppb_ts(2)%sf(j, k, l, q, i) & + 2._wpdtrhs_pb(j, k, l, q, i))/3._wp end do end do end do end do end do end if if (qbmm .and. (.not. polytropic)) then $:GPU_PARALLEL_LOOP(collapse=5) do i = 1, nb do l = 0, p do k = 0, n do j = 0, m do q = 1, nnode mv_ts(1)%sf(j, k, l, q, i) = & (mv_ts(1)%sf(j, k, l, q, i) & + 2._wpmv_ts(2)%sf(j, k, l, q, i) & + 2._wpdtrhs_mv(j, k, l, q, i))/3._wp end do end do end do end do end do end if if (bodyForces) call s_apply_bodyforces(q_cons_ts(dest)%vf, q_prim_vf, rhs_vf, 2._wpdt/3._wp) if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(dest)%vf) if (model_eqns == 3 .and. (.not. relax)) then call s_pressure_relaxation_procedure(q_cons_ts(dest)%vf) end if call nvtxStartRange("RHS-ELASTIC") if (hyperelasticity) call s_hyperelastic_rmt_stress_update(q_cons_ts(dest)%vf, q_prim_vf) call nvtxEndRange if (adv_n) call s_comp_alpha_from_n(q_cons_ts(dest)%vf) if (ib) then if (moving_immersed_boundary_flag) then do i = 1, num_ibs if (patch_ib(i)%moving_ibm == 1) then do j = 1, 3 patch_ib(i)%vel(j) = (2._wppatch_ib(i)%vel(j) + patch_ib(i)%step_vel(j) + 0.0dt)/3._wp patch_ib(i)%angular_vel(j) = (2._wppatch_ib(i)%angular_vel(j) + patch_ib(i)%step_angular_vel(j) + 0.0dt)/3._wp patch_ib(i)%angles(j) = (2._wppatch_ib(i)%angles(j) + patch_ib(i)%step_angles(j) + 2._wppatch_ib(i)%angular_vel(j)dt)/3._wp end do patch_ib(i)%x_centroid = (2._wppatch_ib(i)%x_centroid + patch_ib(i)%step_x_centroid + 2._wppatch_ib(i)%vel(1)dt)/3._wp patch_ib(i)%y_centroid = (2._wppatch_ib(i)%y_centroid + patch_ib(i)%step_y_centroid + 2._wppatch_ib(i)%vel(2)dt)/3._wp patch_ib(i)%z_centroid = (2._wppatch_ib(i)%z_centroid + patch_ib(i)%step_z_centroid + 2._wppatch_ib(i)%vel(3)dt)/3._wp end if end do call s_update_mib(num_ibs, levelset, levelset_norm) end if if (qbmm .and. .not. polytropic) then Possible Issue The new cross_product uses hard-coded real(8) while the rest of the module uses real(wp). This type mismatch can lead to precision inconsistencies or implicit conversions. Align the routine and its interface with real(wp). impure subroutine s_finalize_ibm_module() @:DEALLOCATE(ib_markers%sf) @:DEALLOCATE(levelset%sf) @:DEALLOCATE(levelset_norm%sf) end subroutine s_finalize_ibm_module function cross_product(a, b) result(c) implicit none real(8), intent(in) :: a(3), b(3) real(8) :: c(3) c(1) = a(2)b(3) - a(3)b(2) c(2) = a(3)b(1) - a(1)b(3) c(3) = a(1)b(2) - a(2)b(1) end function cross_product end module m_ibm Rotation Order/2D-3D Consistency The rotation matrix composition assumes a fixed X→Y→Z order; inverse is built via transposed order. Verify this matches the intended convention and dimensionality. In 2D, only Z rotation is used; ensure rectangles and angular_vel usage are consistent with this choice and document expected angle units and order.

qodo-merge-pro · 2025-10-07T16:36:13Z

src/simulation/m_time_steppers.fpp

High-level Suggestion

The logic for updating Immersed Boundary (IB) motion is duplicated across several time-stepping subroutines in m_time_steppers.fpp. This should be refactored into a single, parameterized subroutine to improve code maintainability. [High-level, importance: 7]

Solution Walkthrough:

Before:

subroutine s_tvd_rk2(...) ... ! Step 1 if (moving_immersed_boundary_flag) then do i = 1, num_ibs if (patch_ib(i)%moving_ibm == 1) then ! ... save state and perform Euler step ... end if end do call s_update_mib(...) end if ... ! Step 2 if (moving_immersed_boundary_flag) then do i = 1, num_ibs if (patch_ib(i)%moving_ibm == 1) then ! ... perform RK2 update step ... end if end do call s_update_mib(...) end if end subroutine

After:

subroutine s_update_ib_position(rk_stage, rk_order) ! Logic to update IB position based on RK stage and order if (moving_immersed_boundary_flag) then do i = 1, num_ibs if (patch_ib(i)%moving_ibm == 1) then ! ... unified update logic using rk_stage/order ... end if end do call s_update_mib(...) end if end subroutine subroutine s_tvd_rk2(...) ... call s_update_ib_position(rk_stage=1, rk_order=2) ... call s_update_ib_position(rk_stage=2, rk_order=2) end subroutine

src/simulation/m_ibm.fpp

codecov · 2025-10-07T19:50:37Z

Codecov Report

❌ Patch coverage is 50.57034% with 130 lines in your changes missing coverage. Please review.
✅ Project coverage is 42.11%. Comparing base (5c9d069) to head (21d0594).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/common/m_ib_patches.fpp	44.18%	63 Missing and 9 partials ⚠️
src/common/m_compute_levelset.fpp	46.87%	41 Missing and 10 partials ⚠️
src/simulation/m_time_steppers.fpp	66.66%	4 Missing and 2 partials ⚠️
src/simulation/m_ibm.fpp	83.33%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1014      +/-   ##
==========================================
+ Coverage   41.82%   42.11%   +0.28%     
==========================================
  Files          70       70              
  Lines       19921    19936      +15     
  Branches     2490     2485       -5     
==========================================
+ Hits         8332     8396      +64     
+ Misses      10043     9996      -47     
+ Partials     1546     1544       -2

☔ 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.

danieljvickers · 2025-10-08T01:26:37Z

It should be noted that even if this passes the final test (which I think it will), it is not ready to be merged. I realized that there are two changes for the actual MIBM code corrections that must be implemented which won't affect current test suite. I will make another comment when this is ready to merge.

… shift and rotate the grid cell into their local frame and then compute validity from the local frame. This is true for all IB patches except for circles and spheres, as rotational symetry solves the problem for me in those cases.

…ment has been left.

… with multiple MPI ranks

examples/2D_ibm_airfoil/case.py

…upport

…PU build failing

danieljvickers · 2025-10-12T21:22:48Z

I checked twice, and the this most-recent commit does build on frontier for CPU when I sign into a login node and run it manually. The failure message is rather non-descript. Does anyone have thoughtss on why the CPU build would suddenly fail like this?

…to add-rotating-mibms

danieljvickers · 2025-10-15T07:52:00Z

CPU frontier test suite timed out after 23 hours without printing logs. Not really sure what to make of that one. I will try to run something locally tomorrow if I can get an allocation. I wasn't able to get an interactive session on frontier all day today.

wilfonba · 2025-10-15T12:35:50Z

CPU frontier test suite timed out after 23 hours without printing logs. Not really sure what to make of that one. I will try to run something locally tomorrow if I can get an allocation. I wasn't able to get an interactive session on frontier all day today.

That just means it wasn't able to get a node to run the tests. I restarted it.

sbryngelson · 2025-10-15T19:24:24Z

We can have 2D and 3D IBMs, so I suppose we need test cases for both?

We also need examples for moving IBMs in 2D and 3D, also an update to the readme.

danieljvickers · 2025-10-15T19:27:41Z

I was planning on adding that (examples and documentation) with a follow up PR getting 2D and 3D cases from literature. Would you prefer that I do that into this PR?

…

On Wed, Oct 15, 2025, 3:24 PM Spencer Bryngelson ***@***.***> wrote: *sbryngelson* left a comment (MFlowCode/MFC#1014) <#1014 (comment)> We can have 2D and 3D IBMs, so I suppose we need test cases for both? We also need examples for moving IBMs in 2D and 3D, also an update to the readme. — Reply to this email directly, view it on GitHub <#1014 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BMWYI3XSNSTCJ3IU2FFRMP33X2NQLAVCNFSM6AAAAACIRBVD4GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTIMBXHE2DKNRTG4> . You are receiving this because you authored the thread.Message ID: ***@***.***>

sbryngelson · 2025-10-15T19:31:06Z

A simple example in this PR is needed (and tests). You can modify the example to be 'better' later if you want, or just add another example. I just don't want dangling features.

danieljvickers · 2025-10-15T19:37:34Z

That sounds fair enough. I have a couple that I have been playing with anyways. I can add a couple.

danieljvickers added 11 commits October 7, 2025 11:59

Pushing from local

d953d5b

Succesfully added RK time stepping

9f6af97

Finalizing time-step module

8d3941d

Addded handlers for new case file parameters

734b649

added hooks for a rotation matrix

e1c3b3f

Added the rotation matrix calculation, which will be needed for ident…

effd324

…ifying rotated markers and velocity calculations

Finished adding velocity update, but now I need to fix compiler errors

ac0013d

Added some inverse rotation calculations to the square case

50799ab

Ready to test with a rotated square

c04cfae

Ran test case and everything appears to be working

e80952b

Spelling and formatting run

6ddc795

danieljvickers requested review from a team as code owners October 7, 2025 16:34

qodo-merge-pro bot added the Review effort 4/5 label Oct 7, 2025

qodo-merge-pro bot reviewed Oct 7, 2025

View reviewed changes

src/simulation/m_ibm.fpp Show resolved Hide resolved

Fixed single-precision building

ce5f071

danieljvickers added 10 commits October 8, 2025 11:13

Updated patch_IB rotation values

007949d

Updated 2D airfoil

ddfcd0f

Fixed compiler issue with airfoil changes. All tests pass locally

372ddf0

Updated rectangle levelset

603b42f

Swapped out Id values

ab29c8a

Modified the cuboid levelset to use local coords

71437fd

Final levelset calculation updated

fcccef7

Fixed compiling errors'

d13de48

Temporary fix for issue involving cuboid normal vector selection. Com…

91e0f50

…ment has been left.

Fixed issue with pre-processing not applying correct initial rotation…

78fd281

… with multiple MPI ranks

wilfonba reviewed Oct 10, 2025

View reviewed changes

examples/2D_ibm_airfoil/case.py Show resolved Hide resolved

danieljvickers added 2 commits October 10, 2025 11:14

Formatting

7cf947e

Adding OpenMP/OpenACC tags for vecturs used in paralel loop for GPU s…

63c87c5

…upport

danieljvickers requested a review from wilfonba October 11, 2025 19:18

danieljvickers added 2 commits October 11, 2025 21:21

found an error in 3d Airfoil case that is unrelated to the frontier C…

1af3670

…PU build failing

Merge branch 'master' into add-rotating-mibms

5a60a7f

wilfonba and others added 7 commits October 13, 2025 08:20

Build with one task for testing

272ca02

.github

8fc4902

Found my error. Thanks Ben for being a rubber ducky :)

64fa336

Integrated with new time stepper code

29db19a

Resolved error. Ready to merge

c29ec88

Merge branch 'add-rotating-mibms' of github.com:danieljvickers/MFC in…

604e1db

…to add-rotating-mibms

Formatting

cbfb2da

sbryngelson approved these changes Oct 15, 2025

View reviewed changes

danieljvickers added 6 commits October 15, 2025 16:39

Added a tumbling rectangle case

7560a02

Fixed GPU errors on NVHPC

8906927

Generated golden files for totating sphere

d4a1280

Updated case file to be longer

876f16a

Formatting

f1b6fb7

Spelling

21d0594

sbryngelson merged commit 1609707 into MFlowCode:master Oct 17, 2025
29 checks passed

danieljvickers deleted the add-rotating-mibms branch October 17, 2025 14:37

Add rotating mibms #1014

Add rotating mibms #1014

Uh oh!

Conversation

danieljvickers commented Oct 7, 2025 • edited by qodo-merge-pro bot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

User description

Description

Type of change

Scope

How Has This Been Tested?

Checklist

If your code changes any code source files (anything in src/simulation)

PR Type

Description

Diagram Walkthrough

File Walkthrough

Uh oh!

qodo-merge-pro bot commented Oct 7, 2025

PR Reviewer Guide 🔍

Uh oh!

qodo-merge-pro bot Oct 7, 2025

Choose a reason for hiding this comment

High-level Suggestion

Solution Walkthrough:

Before:

After:

Uh oh!

Uh oh!

codecov bot commented Oct 7, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

danieljvickers commented Oct 8, 2025

Uh oh!

Uh oh!

danieljvickers commented Oct 12, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

danieljvickers commented Oct 15, 2025

Uh oh!

wilfonba commented Oct 15, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

sbryngelson commented Oct 15, 2025

Uh oh!

danieljvickers commented Oct 15, 2025 via email

Uh oh!

sbryngelson commented Oct 15, 2025

Uh oh!

danieljvickers commented Oct 15, 2025

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Milestone

Development

Uh oh!

3 participants

danieljvickers commented Oct 7, 2025 •

edited by qodo-merge-pro bot

Loading

If your code changes any code source files (anything in `src/simulation`)

codecov bot commented Oct 7, 2025 •

edited

Loading

danieljvickers commented Oct 12, 2025 •

edited

Loading

wilfonba commented Oct 15, 2025 •

edited

Loading