Skip to content

Commit ac0013d

Browse files
Finished adding velocity update, but now I need to fix compiler errors
1 parent effd324 commit ac0013d

File tree

1 file changed

+17
-1
lines changed

1 file changed

+17
-1
lines changed

src/simulation/m_ibm.fpp

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -182,6 +182,8 @@ contains
182182
real(wp), dimension(3) :: norm !< Normal vector from GP to IP
183183
real(wp), dimension(3) :: physical_loc !< Physical loc of GP
184184
real(wp), dimension(3) :: vel_g !< Velocity of GP
185+
real(wp), dimension(3) :: radial_vector !< vector from centroid to ghost point
186+
real(wp), dimension(3) :: rotation_velocity !< speed of the ghost point due to rotation
185187

186188
real(wp) :: nbub
187189
real(wp) :: buf
@@ -265,9 +267,13 @@ contains
265267
! we know the object is not moving if moving_ibm is 0 (false)
266268
vel_g = 0._wp
267269
else
270+
radial_vector = physical_loc - (patch_ib(patch_id)%x_centroid, &
271+
patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid)
272+
rotational_velocity = cross_product(matmul(patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector)
268273
do q = 1, 3
269274
! if mibm is 1 or 2, then the boundary may be moving
270-
vel_g(q) = patch_ib(patch_id)%vel(q)
275+
vel_g(q) = patch_ib(patch_id)%vel(q) ! add the linear velocity
276+
vel_g(q) = vel_g(q) + rotational_velocity(q) ! add the rotational velocity
271277
end do
272278
end if
273279
end if
@@ -978,4 +984,14 @@ contains
978984
979985
end subroutine s_finalize_ibm_module
980986
987+
function cross_product(a, b) result(c)
988+
implicit none
989+
real(8), intent(in) :: a(3), b(3)
990+
real(8) :: c(3)
991+
992+
c(1) = a(2)*b(3) - a(3)*b(2)
993+
c(2) = a(3)*b(1) - a(1)*b(3)
994+
c(3) = a(1)*b(2) - a(2)*b(1)
995+
end function cross_product
996+
981997
end module m_ibm

0 commit comments

Comments
 (0)