Skip to content

Commit 9ae77c5

Browse files
Claude initial review of the changes for syntactical bugs
1 parent 81ed411 commit 9ae77c5

1 file changed

Lines changed: 11 additions & 8 deletions

File tree

src/simulation/m_ibm.fpp

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -175,16 +175,18 @@ contains
175175

176176
end subroutine s_compute_ghost_point_pressure
177177

178-
subroutine s_compute_ghost_point_velocity(gp, gp_patch_id, radial_vector, vel_IP, vel_GP)
178+
subroutine s_compute_ghost_point_velocity(gp, gp_patch_id, radial_vector, vel_IP, pres_IP, vel_GP)
179179

180180
$:GPU_ROUTINE(parallelism='[seq]')
181181

182182
type(ghost_point), intent(in) :: gp
183183
integer, intent(in) :: gp_patch_id
184184
real(wp), dimension(3), intent(in) :: radial_vector, vel_IP
185+
real(wp), intent(in) :: pres_IP
185186
real(wp), dimension(3), intent(out) :: vel_GP
186187
real(wp), dimension(3) :: norm, vel_norm_IP, rotation_velocity
187188
real(wp) :: buf, v_blow_eff
189+
integer :: q !< Iterator variable
188190

189191
! Calculate velocity of ghost cell
190192
if (gp%slip) then
@@ -210,24 +212,25 @@ contains
210212
call s_cross_product(patch_ib(gp_patch_id)%angular_vel, radial_vector, rotation_velocity)
211213
do q = 1, 3
212214
! if mibm is 1 or 2, then the boundary may be moving
213-
vel_g(q) = patch_ib(gp_patch_id)%vel(q) ! add the linear velocity
215+
vel_GP(q) = patch_ib(gp_patch_id)%vel(q) ! add the linear velocity
214216
vel_GP(q) = vel_GP(q) + rotation_velocity(q) ! add the rotational velocity
215217
end do
216218
end if
217219
end if
218220
219221
! Burning/injecting surface: superimpose wall-normal (outward) blowing on the
220222
! ghost velocity so the immersed surface transpires/injects gas into the flow.
221-
if (patch_ib(patch_id)%v_blow > 0._wp) then
222-
v_blow_eff = patch_ib(patch_id)%v_blow
223+
if (patch_ib(gp_patch_id)%v_blow > 0._wp) then
224+
v_blow_eff = patch_ib(gp_patch_id)%v_blow
223225
! Pressure-coupled burn rate (Vieille's law r_dot ~ p^n)
224-
if (patch_ib(patch_id)%burn_rate_pref > 0._wp) then
226+
if (patch_ib(gp_patch_id)%burn_rate_pref > 0._wp) then
225227
! max(pres_IP, 0) guards the fractional power against a transient negative
226-
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
228+
v_blow_eff = v_blow_eff*(max(pres_IP, &
229+
& 0._wp)/patch_ib(gp_patch_id)%burn_rate_pref)**patch_ib(gp_patch_id)%burn_rate_exp
227230
end if
228231
norm(1:3) = gp%levelset_norm
229232
buf = sqrt(sum(norm**2))
230-
if (buf > 0._wp) vel_g = vel_g + v_blow_eff*norm/buf
233+
if (buf > 0._wp) vel_GP = vel_GP + v_blow_eff*norm/buf
231234
end if
232235

233236
end subroutine s_compute_ghost_point_velocity
@@ -397,7 +400,7 @@ contains
397400
& + real(ghost_points(i)%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg))
398401

399402
! Calculate velocity of ghost cell
400-
call s_compute_ghost_point_velocity(gp, patch_id, radial_vector, vel_IP, vel_g)
403+
call s_compute_ghost_point_velocity(gp, patch_id, radial_vector, vel_IP, pres_IP, vel_g)
401404

402405
! Set momentum
403406
vel_sum_g = 0._wp

0 commit comments

Comments
 (0)