From 81ed41197f42bacf692c8e9c7216a17702520628 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Thu, 10 Sep 2026 08:05:18 -0400 Subject: [PATCH 1/3] Initial commit --- src/simulation/m_ibm.fpp | 189 +++++++++++++++++++++++---------------- 1 file changed, 112 insertions(+), 77 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 99a212260..f34415836 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -139,6 +139,99 @@ contains end subroutine s_ibm_setup + subroutine s_compute_ghost_point_pressure(gp, gp_patch_id, alpha_rho_IP, pres_IP, pres_GP) + + $:GPU_ROUTINE(parallelism='[seq]') + + type(ghost_point), intent(in) :: gp + integer, intent(in) :: gp_patch_id + #:if not MFC_CASE_OPTIMIZATION and USING_AMD + real(wp), dimension(3), intent(in) :: alpha_rho_IP + #:else + real(wp), dimension(num_fluids), intent(in) :: alpha_rho_IP + #:endif + real(wp), intent(in) :: pres_IP + real(wp), intent(out) :: pres_GP + integer :: q !< Iterator variable + + pres_GP = 0._wp + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_fluids + select case (eoss(q)) + case (eos_ideal_gas, eos_stiffened_gas) + ! Pressure correction for moving IB: accounts for acceleration of IB surface + pres_GP = pres_GP + pres_IP/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & + & *dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, & + & gp%levelset_norm)) + case default + ! TODO: moving-IB pressure correction is not yet derived for state-dependent EOS + ! (Mie-Gruneisen/JWL/Vinet); wire up the correct formula in a follow-up PR. +#ifndef MFC_GPU + call s_mpi_abort('s_compute_ghost_point_pressure: moving IB pressure correction is only ' & + & // 'implemented for eos_ideal_gas/eos_stiffened_gas') +#endif + end select + end do + + end subroutine s_compute_ghost_point_pressure + + subroutine s_compute_ghost_point_velocity(gp, gp_patch_id, radial_vector, vel_IP, vel_GP) + + $:GPU_ROUTINE(parallelism='[seq]') + + type(ghost_point), intent(in) :: gp + integer, intent(in) :: gp_patch_id + real(wp), dimension(3), intent(in) :: radial_vector, vel_IP + real(wp), dimension(3), intent(out) :: vel_GP + real(wp), dimension(3) :: norm, vel_norm_IP, rotation_velocity + real(wp) :: buf, v_blow_eff + + ! Calculate velocity of ghost cell + if (gp%slip) then + norm(1:3) = gp%levelset_norm + buf = sqrt(sum(norm**2)) + norm = norm/buf + vel_norm_IP = sum(vel_IP*norm)*norm + vel_GP = vel_IP - vel_norm_IP + if (patch_ib(gp_patch_id)%moving_ibm /= 0) then + ! compute the linear velocity of the ghost point due to rotation + call s_cross_product(patch_ib(gp_patch_id)%angular_vel, radial_vector, rotation_velocity) + + ! add only the component of the IB's motion that is normal to the surface + vel_GP = vel_GP + sum((patch_ib(gp_patch_id)%vel + rotation_velocity)*norm)*norm + end if + else + if (patch_ib(gp_patch_id)%moving_ibm == 0) then + ! we know the object is not moving if moving_ibm is 0 (false) + vel_GP = 0._wp + else + ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear + ! velocity + call s_cross_product(patch_ib(gp_patch_id)%angular_vel, radial_vector, rotation_velocity) + do q = 1, 3 + ! if mibm is 1 or 2, then the boundary may be moving + vel_g(q) = patch_ib(gp_patch_id)%vel(q) ! add the linear velocity + vel_GP(q) = vel_GP(q) + rotation_velocity(q) ! add the rotational velocity + end do + end if + end if + + ! Burning/injecting surface: superimpose wall-normal (outward) blowing on the + ! ghost velocity so the immersed surface transpires/injects gas into the flow. + if (patch_ib(patch_id)%v_blow > 0._wp) then + v_blow_eff = patch_ib(patch_id)%v_blow + ! Pressure-coupled burn rate (Vieille's law r_dot ~ p^n) + if (patch_ib(patch_id)%burn_rate_pref > 0._wp) then + ! max(pres_IP, 0) guards the fractional power against a transient negative + v_blow_eff = v_blow_eff*(max(pres_IP, 0._wp)/patch_ib(patch_id)%burn_rate_pref)**patch_ib(patch_id)%burn_rate_exp + end if + norm(1:3) = gp%levelset_norm + buf = sqrt(sum(norm**2)) + if (buf > 0._wp) vel_g = vel_g + v_blow_eff*norm/buf + end if + + end subroutine s_compute_ghost_point_velocity + !> Update the conservative variables at the ghost points subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in) @@ -152,7 +245,7 @@ contains real(wp), dimension(2) :: Re_K real(wp) :: G_K real(wp) :: qv_K - real(wp) :: pres_IP + real(wp) :: pres_IP, pres_GP real(wp), dimension(3) :: vel_IP, vel_norm_IP real(wp) :: c_IP @@ -221,10 +314,10 @@ contains $:END_GPU_PARALLEL_LOOP() if (num_gps > 0) then - $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, & - & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & - & innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id, & - & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, vel_sum_g, E_ghost, alpha_q, alpha_rho_q, e_q]') + $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, pres_GP, vel_IP, vel_g, & + & vel_norm_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, & + & G_K, Gs, gp, innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, & + & patch_id, Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, vel_sum_g, E_ghost, alpha_q, alpha_rho_q, e_q]') do i = 1, num_gps gp = ghost_points(i) j = gp%loc(1) @@ -233,11 +326,8 @@ contains patch_id = ghost_points(i)%ib_patch_id ! Calculate physical location of GP - if (p > 0) then - physical_loc = [x_cc(j), y_cc(k), z_cc(l)] - else - physical_loc = [x_cc(j), y_cc(k), 0._wp] - end if + physical_loc = [x_cc(j), y_cc(k), 0._wp] + if (num_dims == 3) physical_loc(3) = z_cc(l) ! Interpolate primitive variables at image point associated w/ GP if (bubbles_euler .and. .not. qbmm) then @@ -285,14 +375,8 @@ contains if (patch_ib(patch_id)%moving_ibm <= 1) then q_prim_vf(eqn_idx%E)%sf(j, k, l) = pres_IP else - q_prim_vf(eqn_idx%E)%sf(j, k, l) = 0._wp - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_fluids - ! Pressure correction for moving IB: accounts for acceleration of IB surface - q_prim_vf(eqn_idx%E)%sf(j, k, l) = q_prim_vf(eqn_idx%E)%sf(j, k, & - & l) + pres_IP/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & - & *dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, gp%levelset_norm)) - end do + call s_compute_ghost_point_pressure(gp, patch_id, alpha_rho_IP, pres_IP, pres_GP) + q_prim_vf(eqn_idx%E)%sf(j, k, l) = pres_GP end if ! If in simulation, use acc mixture subroutines @@ -303,64 +387,17 @@ contains call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_K, alpha_IP, alpha_rho_IP, Re_K) end if - if (patch_ib(patch_id)%moving_ibm /= 0) then - ! get the vector that points from the centroid to the ghost - radial_vector(1) = physical_loc(1) - (patch_ib(patch_id)%x_centroid + real(ghost_points(i)%x_periodicity, & - & wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)) - radial_vector(2) = physical_loc(2) - (patch_ib(patch_id)%y_centroid + real(ghost_points(i)%y_periodicity, & - & wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)) - radial_vector(3) = 0._wp - if (num_dims == 3) radial_vector(3) = physical_loc(3) - (patch_ib(patch_id)%z_centroid & - & + real(ghost_points(i)%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)) - end if + ! get the vector that points from the centroid to the ghost + radial_vector(1) = physical_loc(1) - (patch_ib(patch_id)%x_centroid + real(ghost_points(i)%x_periodicity, & + & wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)) + radial_vector(2) = physical_loc(2) - (patch_ib(patch_id)%y_centroid + real(ghost_points(i)%y_periodicity, & + & wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)) + radial_vector(3) = 0._wp + if (num_dims == 3) radial_vector(3) = physical_loc(3) - (patch_ib(patch_id)%z_centroid & + & + real(ghost_points(i)%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)) ! Calculate velocity of ghost cell - if (gp%slip) then - norm(1:3) = gp%levelset_norm - buf = sqrt(sum(norm**2)) - norm = norm/buf - vel_norm_IP = sum(vel_IP*norm)*norm - vel_g = vel_IP - vel_norm_IP - if (patch_ib(patch_id)%moving_ibm /= 0) then - ! compute the linear velocity of the ghost point due to rotation - call s_cross_product(patch_ib(patch_id)%angular_vel, radial_vector, rotation_velocity) - - ! add only the component of the IB's motion that is normal to the surface - vel_g = vel_g + sum((patch_ib(patch_id)%vel + rotation_velocity)*norm)*norm - end if - else - if (patch_ib(patch_id)%moving_ibm == 0) then - ! we know the object is not moving if moving_ibm is 0 (false) - vel_g = 0._wp - else - ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear - ! velocity - call s_cross_product(patch_ib(patch_id)%angular_vel, radial_vector, rotation_velocity) - do q = 1, 3 - ! if mibm is 1 or 2, then the boundary may be moving - vel_g(q) = patch_ib(patch_id)%vel(q) ! add the linear velocity - vel_g(q) = vel_g(q) + rotation_velocity(q) ! add the rotational velocity - end do - end if - end if - - ! Burning/injecting surface: superimpose wall-normal (outward) blowing on the - ! ghost velocity so the immersed surface transpires/injects gas into the flow. - if (patch_ib(patch_id)%v_blow > 0._wp) then - v_blow_eff = patch_ib(patch_id)%v_blow - ! Pressure-coupled burn rate (Vieille's law r_dot ~ p^n): the local surface - ! pressure scales the blowing speed, giving chamber-pressure feedback (internal - ! ballistics) in a closed chamber. Off (constant) when burn_rate_pref <= 0. - if (patch_ib(patch_id)%burn_rate_pref > 0._wp) then - ! max(pres_IP, 0) guards the fractional power against a transient negative - ! interpolated pressure, which would otherwise return NaN and poison the field. - v_blow_eff = v_blow_eff*(max(pres_IP, & - & 0._wp)/patch_ib(patch_id)%burn_rate_pref)**patch_ib(patch_id)%burn_rate_exp - end if - norm(1:3) = gp%levelset_norm - buf = sqrt(sum(norm**2)) - if (buf > 0._wp) vel_g = vel_g + v_blow_eff*norm/buf - end if + call s_compute_ghost_point_velocity(gp, patch_id, radial_vector, vel_IP, vel_g) ! Set momentum vel_sum_g = 0._wp @@ -379,9 +416,7 @@ contains end do ! Set color function - if (surface_tension) then - q_cons_vf(eqn_idx%c)%sf(j, k, l) = c_IP - end if + if (surface_tension) q_cons_vf(eqn_idx%c)%sf(j, k, l) = c_IP ! Set Energy if (chemistry) then From 9ae77c5538ff91ba68bfe4a8e4eeae08d4355ac7 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Thu, 10 Sep 2026 08:10:26 -0400 Subject: [PATCH 2/3] Claude initial review of the changes for syntactical bugs --- src/simulation/m_ibm.fpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index f34415836..1895b7ade 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -175,16 +175,18 @@ contains end subroutine s_compute_ghost_point_pressure - subroutine s_compute_ghost_point_velocity(gp, gp_patch_id, radial_vector, vel_IP, vel_GP) + subroutine s_compute_ghost_point_velocity(gp, gp_patch_id, radial_vector, vel_IP, pres_IP, vel_GP) $:GPU_ROUTINE(parallelism='[seq]') type(ghost_point), intent(in) :: gp integer, intent(in) :: gp_patch_id real(wp), dimension(3), intent(in) :: radial_vector, vel_IP + real(wp), intent(in) :: pres_IP real(wp), dimension(3), intent(out) :: vel_GP real(wp), dimension(3) :: norm, vel_norm_IP, rotation_velocity real(wp) :: buf, v_blow_eff + integer :: q !< Iterator variable ! Calculate velocity of ghost cell if (gp%slip) then @@ -210,7 +212,7 @@ contains call s_cross_product(patch_ib(gp_patch_id)%angular_vel, radial_vector, rotation_velocity) do q = 1, 3 ! if mibm is 1 or 2, then the boundary may be moving - vel_g(q) = patch_ib(gp_patch_id)%vel(q) ! add the linear velocity + vel_GP(q) = patch_ib(gp_patch_id)%vel(q) ! add the linear velocity vel_GP(q) = vel_GP(q) + rotation_velocity(q) ! add the rotational velocity end do end if @@ -218,16 +220,17 @@ contains ! Burning/injecting surface: superimpose wall-normal (outward) blowing on the ! ghost velocity so the immersed surface transpires/injects gas into the flow. - if (patch_ib(patch_id)%v_blow > 0._wp) then - v_blow_eff = patch_ib(patch_id)%v_blow + if (patch_ib(gp_patch_id)%v_blow > 0._wp) then + v_blow_eff = patch_ib(gp_patch_id)%v_blow ! Pressure-coupled burn rate (Vieille's law r_dot ~ p^n) - if (patch_ib(patch_id)%burn_rate_pref > 0._wp) then + if (patch_ib(gp_patch_id)%burn_rate_pref > 0._wp) then ! max(pres_IP, 0) guards the fractional power against a transient negative - v_blow_eff = v_blow_eff*(max(pres_IP, 0._wp)/patch_ib(patch_id)%burn_rate_pref)**patch_ib(patch_id)%burn_rate_exp + v_blow_eff = v_blow_eff*(max(pres_IP, & + & 0._wp)/patch_ib(gp_patch_id)%burn_rate_pref)**patch_ib(gp_patch_id)%burn_rate_exp end if norm(1:3) = gp%levelset_norm buf = sqrt(sum(norm**2)) - if (buf > 0._wp) vel_g = vel_g + v_blow_eff*norm/buf + if (buf > 0._wp) vel_GP = vel_GP + v_blow_eff*norm/buf end if end subroutine s_compute_ghost_point_velocity @@ -397,7 +400,7 @@ contains & + real(ghost_points(i)%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)) ! Calculate velocity of ghost cell - call s_compute_ghost_point_velocity(gp, patch_id, radial_vector, vel_IP, vel_g) + call s_compute_ghost_point_velocity(gp, patch_id, radial_vector, vel_IP, pres_IP, vel_g) ! Set momentum vel_sum_g = 0._wp From 3699ed84699ccf3479106f3055daf567b62bb578 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Sun, 13 Sep 2026 00:35:49 -0400 Subject: [PATCH 3/3] Addressed review comments --- src/simulation/m_ibm.fpp | 38 ++++++++++----------------------- toolchain/mfc/case_validator.py | 2 +- 2 files changed, 12 insertions(+), 28 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 1895b7ade..dcbf16247 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -26,7 +26,7 @@ module m_ibm implicit none private :: s_compute_image_points, s_compute_interpolation_coeffs, s_interpolate_image_point, s_find_ghost_points, & - & s_find_num_ghost_points + & s_find_num_ghost_points, s_compute_ghost_point_pressure, s_compute_ghost_point_velocity ; public :: s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module type(integer_field), public :: ib_markers @@ -157,20 +157,9 @@ contains pres_GP = 0._wp $:GPU_LOOP(parallelism='[seq]') do q = 1, num_fluids - select case (eoss(q)) - case (eos_ideal_gas, eos_stiffened_gas) - ! Pressure correction for moving IB: accounts for acceleration of IB surface - pres_GP = pres_GP + pres_IP/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & - & *dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, & - & gp%levelset_norm)) - case default - ! TODO: moving-IB pressure correction is not yet derived for state-dependent EOS - ! (Mie-Gruneisen/JWL/Vinet); wire up the correct formula in a follow-up PR. -#ifndef MFC_GPU - call s_mpi_abort('s_compute_ghost_point_pressure: moving IB pressure correction is only ' & - & // 'implemented for eos_ideal_gas/eos_stiffened_gas') -#endif - end select + ! Pressure correction for moving IB: accounts for acceleration of IB surface + pres_GP = pres_GP + pres_IP/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & + & *dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, gp%levelset_norm)) end do end subroutine s_compute_ghost_point_pressure @@ -249,7 +238,7 @@ contains real(wp) :: G_K real(wp) :: qv_K real(wp) :: pres_IP, pres_GP - real(wp), dimension(3) :: vel_IP, vel_norm_IP + real(wp), dimension(3) :: vel_IP real(wp) :: c_IP #:if not MFC_CASE_OPTIMIZATION and USING_AMD @@ -269,18 +258,13 @@ contains #:endif real(wp) :: alpha_q, alpha_rho_q, e_q real(wp) :: T_IP, mw_IP, e_IP !< Image-point temperature, mixture MW, and mass-specific internal energy (chemistry) - real(wp) :: v_blow_eff !< Effective surface blowing speed (after any pressure-coupled burn-rate scaling) ! Primitive variables at the image point associated with a ghost point, interpolated from surrounding fluid cells. - real(wp), dimension(3) :: norm !< Normal vector from GP to IP - real(wp), dimension(3) :: physical_loc !< Physical loc of GP - real(wp), dimension(3) :: vel_g !< Velocity of GP - real(wp), dimension(3) :: radial_vector !< vector from centroid to ghost point - real(wp), dimension(3) :: rotation_velocity !< speed of the ghost point due to rotation + real(wp), dimension(3) :: physical_loc !< Physical loc of GP + real(wp), dimension(3) :: vel_g !< Velocity of GP + real(wp), dimension(3) :: radial_vector !< vector from centroid to ghost point real(wp) :: nbub - real(wp) :: buf type(ghost_point) :: gp - type(ghost_point) :: innerp ! set the Moving IBM interior conservative variables $:GPU_PARALLEL_LOOP(private='[i, j, k, patch_id, rho]', collapse=3) @@ -318,9 +302,9 @@ contains if (num_gps > 0) then $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, pres_GP, vel_IP, vel_g, & - & vel_norm_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, & - & G_K, Gs, gp, innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, & - & patch_id, Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, vel_sum_g, E_ghost, alpha_q, alpha_rho_q, e_q]') + & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & + & radial_vector, j, k, l, q, qv_K, c_IP, nbub, patch_id, Ys_IP, T_IP, mw_IP, e_IP, vel_sum_g, & + & E_ghost, alpha_q, alpha_rho_q, e_q]') do i = 1, num_gps gp = ghost_points(i) j = gp%loc(1) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index b15cbffb9..956a10d16 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -1069,7 +1069,7 @@ def check_eos_selector(self): self.prohibit(self.get("wave_speeds") == 2, "a state-dependent eos (mie_gruneisen, jwl, vinet) requires wave_speeds = 1 (the PVRS estimate is stiffened-gas only)") for j in range(1, (self.get("num_patches") or 0) + 1): self.prohibit(self.get(f"patch_icpp({j})%hcid") in (202, 203), f"patch_icpp({j})%hcid = 202/203 read fluid_pp(1)%gamma, which a state-dependent eos does not set") - for flag in ("bubbles_euler", "bubbles_lagrange", "igr", "relativity", "mhd", "chemistry", "relax"): + for flag in ("bubbles_euler", "bubbles_lagrange", "igr", "relativity", "mhd", "chemistry", "relax", "ib"): self.prohibit(self.get(flag, "F") == "T", f"a state-dependent eos (mie_gruneisen, jwl, vinet) is not supported with {flag} = T") def check_stiffened_eos(self):