Ran test case and everything appears to be working

Author: danieljvickers

src/common/m_derived_types.fpp

@@ -296,6 +296,7 @@ module m_derived_types
         real(wp), dimension(1:3) :: angles
         real(wp), dimension(1:3) :: step_angles
         real(wp), dimension(1:3, 1:3) :: rotation_matrix !< matrix that converts from IB reference frame to fluid reference frame
+        real(wp), dimension(1:3, 1:3) :: rotation_matrix_inverse !< matrix that converts from fluid reference frame to IB reference frame
 
         real(wp) :: c, p, t, m
 

src/common/m_ib_patches.fpp

@@ -27,7 +27,7 @@ module m_ib_patches
 
     implicit none
 
-    private; public :: s_apply_ib_patches
+    private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix
 
     real(wp) :: x_centroid, y_centroid, z_centroid
     real(wp) :: length_x, length_y, length_z
@@ -546,7 +546,7 @@ contains
         do j = 0, n
             do i = 0, m
                 ! get the x and y coodinates in the local IB frame
-                xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid]
+                xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp]
                 xy_local = matmul(inverse_rotation, xy_local)
                 if (x_boundary%beg <= xy_local(1) .and. &
                     x_boundary%end >= xy_local(1) .and. &
@@ -993,6 +993,51 @@ contains
 
     end subroutine s_ib_model
 
+    !> Subroutine that computes a rotation matrix for converting to the rotating frame of the boundary
+    subroutine s_update_ib_rotation_matrix(patch_id)
+
+        integer, intent(in) :: patch_id
+        integer :: i
+
+        real(wp), dimension(3, 3, 3) :: rotation
+        real(wp) :: angle
+
+        ! construct the x, y, and z rotation matrices
+        if (num_dims == 3) then
+          ! also compute the x and y axes in 3D
+          angle = patch_ib(patch_id)%angles(1)
+          rotation(1, 1, :) = [1._wp, 0._wp     , 0._wp      ]
+          rotation(1, 2, :) = [0._wp, cos(angle), -sin(angle)]
+          rotation(1, 3, :) = [0._wp, sin(angle), cos(angle) ]
+
+          angle = patch_ib(patch_id)%angles(2)
+          rotation(2, 1, :) = [cos(angle) , 0._wp, sin(angle)]
+          rotation(2, 2, :) = [0._wp      , 1._wp, 0._wp     ]
+          rotation(2, 3, :) = [-sin(angle), 0._wp, cos(angle)]
+
+          ! apply the y rotation to the x rotation
+          patch_ib(patch_id)%rotation_matrix(:, :) = matmul(rotation(1, :, :), rotation(2, :, :))
+          patch_ib(patch_id)%rotation_matrix_inverse(:, :) = matmul(transpose(rotation(2, :, :)), transpose(rotation(1, :, :)))
+        end if
+
+        ! z component first, since it applies in 2D and 3D
+        angle = patch_ib(patch_id)%angles(3)
+        rotation(3, 1, :) = [cos(angle), -sin(angle), 0._wp]
+        rotation(3, 2, :) = [sin(angle), cos(angle) , 0._wp]
+        rotation(3, 3, :) = [0._wp     , 0._wp      , 1._wp]
+
+        if (num_dims == 3) then
+          ! apply the z rotation to the xy rotation in 3D
+          patch_ib(patch_id)%rotation_matrix(:, :) = matmul(patch_ib(patch_id)%rotation_matrix(:, :), rotation(3, :, :))
+          patch_ib(patch_id)%rotation_matrix_inverse(:, :) = matmul(transpose(rotation(3, :, :)), patch_ib(patch_id)%rotation_matrix_inverse(:, :))
+        else
+          ! write out only the z rotation in 2D
+          patch_ib(patch_id)%rotation_matrix(:, :) = rotation(3, :, :)
+          patch_ib(patch_id)%rotation_matrix_inverse(:, :) = transpose(rotation(3, :, :))
+        end if
+
+    end subroutine s_update_ib_rotation_matrix
+
     subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z)
         $:GPU_ROUTINE(parallelism='[seq]')
 

src/pre_process/m_global_parameters.fpp

@@ -550,6 +550,7 @@ contains
             patch_ib(i)%rotation_matrix(1, 1) = 1._wp
             patch_ib(i)%rotation_matrix(2, 2) = 1._wp
             patch_ib(i)%rotation_matrix(3, 3) = 1._wp
+            patch_ib(i)%rotation_matrix_inverse = patch_ib(i)%rotation_matrix
         end do
 
         ! Fluids physical parameters

src/pre_process/m_initial_condition.fpp

@@ -180,6 +180,8 @@ contains
         !!              primitive variables are converted to conservative ones.
     impure subroutine s_generate_initial_condition
 
+        integer :: i
+
         ! Converting the conservative variables to the primitive ones given
         ! preexisting initial condition data files were read in on start-up
         if (old_ic) then
@@ -190,6 +192,9 @@ contains
         end if
 
         if (ib) then
+            do i = 1, num_ibs
+              ! call s_update_ib_rotation_matrix(i)
+            end do
             call s_apply_ib_patches(ib_markers%sf, levelset, levelset_norm)
         end if
         call s_apply_icpp_patches(patch_id_fp, q_prim_vf)

src/simulation/m_ibm.fpp

@@ -905,50 +905,6 @@ contains
 
     end subroutine s_interpolate_image_point
 
-    !> Subroutine that computes a rotation matrix for converting to the rotating frame of the boundary
-    impure subroutine s_update_ib_rotation_matrix(patch_id)
-
-        integer, intent(in) :: patch_id
-        integer :: i
-
-        real(wp), dimension(3, 3, 3) :: rotation
-        real(wp) :: angle
-
-        ! construct the x, y, and z rotation matrices
-        if (num_dims == 3) then
-          ! also compute the x and y axes in 3D
-          angle = patch_ib(patch_id)%angles(1)
-          rotation(1, 1, :) = [1._wp, 0._wp     , 0._wp      ]
-          rotation(1, 2, :) = [0._wp, cos(angle), -sin(angle)]
-          rotation(1, 3, :) = [0._wp, sin(angle), cos(angle) ]
-
-          angle = patch_ib(patch_id)%angles(2)
-          rotation(2, 1, :) = [cos(angle) , 0._wp, sin(angle)]
-          rotation(2, 2, :) = [0._wp      , 1._wp, 0._wp     ]
-          rotation(2, 3, :) = [-sin(angle), 0._wp, cos(angle)]
-
-          ! apply the y rotation to the x rotation
-          patch_ib(patch_id)%rotation_matrix(:, :) = matmul(rotation(1, :, :), rotation(2, :, :))
-          patch_ib(patch_id)%rotation_matrix_inverse(:, :) = matmul(tranpose(rotation(2, :, :)), transpose(rotation(1, :, :)))
-        end if
-        ! z component first, since it applies in 2D and 3D
-        angle = patch_ib(patch_id)%angles(3)
-        rotation(3, 1, :) = [cos(angle), -sin(angle), 0._wp]
-        rotation(3, 2, :) = [sin(angle), cos(angle) , 0._wp]
-        rotation(3, 3, :) = [0._wp     , 0._wp      , 1._wp]
-
-        if (num_dims == 3) then
-          ! apply the z rotation to the xy rotation in 3D
-          patch_ib(patch_id)%rotation_matrix(:, :) = matmul(patch_ib(patch_id)%rotation_matrix(:, :), rotation(3, :, :))
-          patch_ib(patch_id)%rotation_matrix_inverse(:, :) = matmul(transpose(rotation(3, :, :)), patch_ib(patch_id)%rotation_matrix_inverse(:, :))
-        else
-          ! write out only the z rotation in 2D
-          patch_ib(patch_id)%rotation_matrix(:, :) = rotation(3, :, :)
-          patch_ib(patch_id)%rotation_matrix_inverse(:, :) = transpose(rotation(3, :, :))
-        end if
-
-    end subroutine s_update_ib_rotation_matrix
-
     !> Resets the current indexes of immersed boundaries and replaces them after updating
     !> the position of each moving immersed boundary
     impure subroutine s_update_mib(num_ibs, levelset, levelset_norm)
@@ -964,8 +920,8 @@ contains
 
         ! recalulcate the rotation matrix based upon the new angles
         do i = 1, num_ibs
-          if patch_ib(i)%moving_ibm call s_update_ib_rotation_matrix(i)
-        end if
+          if (patch_ib(i)%moving_ibm .ne. 0) call s_update_ib_rotation_matrix(i)
+        end do
 
         ! recompute the new ib_patch locations and broadcast them.
         call s_apply_ib_patches(ib_markers%sf(0:m, 0:n, 0:p), levelset, levelset_norm)

src/simulation/m_time_steppers.fpp

@@ -836,7 +836,7 @@ contains
                 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.0*dt)/2._wp
-                        patch_ib(i)%angular_vel(j) = (patch_ib(i)%angular_vel(j) + patch_ib(i)%angular_step_vel(j) + 0.0*dt)/2._wp
+                        patch_ib(i)%angular_vel(j) = (patch_ib(i)%angular_vel(j) + patch_ib(i)%step_angular_vel(j) + 0.0*dt)/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
@@ -1235,7 +1235,7 @@ contains
                         patch_ib(i)%vel(j) = (2._wp*patch_ib(i)%vel(j) + patch_ib(i)%step_vel(j) + 0.0*dt)/3._wp
                         patch_ib(i)%angular_vel(j) = (2._wp*patch_ib(i)%angular_vel(j) + patch_ib(i)%step_angular_vel(j) + 0.0*dt)/3._wp
 
-                        patch_ib(i)%angles(j) = (2._wp*patch_ib(i)%angles(j) + patch_ib(i)%step_angles(j) + patch_ib(i)%angular_vel(j)*dt)/3._wp
+                        patch_ib(i)%angles(j) = (2._wp*patch_ib(i)%angles(j) + patch_ib(i)%step_angles(j) + 2._wp*patch_ib(i)%angular_vel(j)*dt)/3._wp
                     end do
 
                     patch_ib(i)%x_centroid = (2._wp*patch_ib(i)%x_centroid + patch_ib(i)%step_x_centroid + 2._wp*patch_ib(i)%vel(1)*dt)/3._wp