Add rotating mibms

src/common/m_derived_types.fpp

@@ -290,6 +290,13 @@ module m_derived_types
         real(wp) :: x_centroid, y_centroid, z_centroid !<
         !! Location of the geometric center, i.e. the centroid, of the patch. It
         !! is specified through its x-, y- and z-coordinates, respectively.
+        real(wp) :: step_x_centroid, step_y_centroid, step_z_centroid !<
+        !! Centroid locations of intermediate steps in the time_stepper module
+
+        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
 
@@ -322,6 +329,9 @@ module m_derived_types
         !! Patch conditions for moving imersed boundaries
         integer :: moving_ibm ! 0 for no moving, 1 for moving, 2 for moving on forced path
         real(wp), dimension(1:3) :: vel
+        real(wp), dimension(1:3) :: step_vel ! velocity array used to store intermediate steps in the time_stepper module
+        real(wp), dimension(1:3) :: angular_vel
+        real(wp), dimension(1:3) :: step_angular_vel ! velocity array used to store intermediate steps in the time_stepper module
 
     end type ib_patch_parameters
 

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
@@ -512,6 +512,8 @@ contains
 
         integer :: i, j, k !< generic loop iterators
         real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters
+        real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
+        real(wp), dimension(1:3, 1:3) :: inverse_rotation
 
         pi_inf = fluid_pp(1)%pi_inf
         gamma = fluid_pp(1)%gamma
@@ -522,13 +524,14 @@ contains
         y_centroid = patch_ib(patch_id)%y_centroid
         length_x = patch_ib(patch_id)%length_x
         length_y = patch_ib(patch_id)%length_y
+        inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
 
         ! Computing the beginning and the end x- and y-coordinates of the
         ! rectangle based on its centroid and lengths
-        x_boundary%beg = x_centroid - 0.5_wp*length_x
-        x_boundary%end = x_centroid + 0.5_wp*length_x
-        y_boundary%beg = y_centroid - 0.5_wp*length_y
-        y_boundary%end = y_centroid + 0.5_wp*length_y
+        x_boundary%beg = -0.5_wp*length_x
+        x_boundary%end = 0.5_wp*length_x
+        y_boundary%beg = -0.5_wp*length_y
+        y_boundary%end = 0.5_wp*length_y
 
         ! Since the rectangular patch does not allow for its boundaries to
         ! be smoothed out, the pseudo volume fraction is set to 1 to ensure
@@ -542,10 +545,13 @@ contains
         ! variables of the current patch are assigned to this cell.
         do j = 0, n
             do i = 0, m
-                if (x_boundary%beg <= x_cc(i) .and. &
-                    x_boundary%end >= x_cc(i) .and. &
-                    y_boundary%beg <= y_cc(j) .and. &
-                    y_boundary%end >= y_cc(j)) then
+                ! get the x and y coordinates in the local IB frame
+                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. &
+                    y_boundary%beg <= xy_local(2) .and. &
+                    y_boundary%end >= xy_local(2)) then
 
                     ! Updating the patch identities bookkeeping variable
                     ib_markers_sf(i, j, 0) = patch_id
@@ -987,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

@@ -541,9 +541,16 @@ contains
 
             ! Variables to handle moving imersed boundaries, defaulting to no movement
             patch_ib(i)%moving_ibm = 0
-            patch_ib(i)%vel(1) = 0._wp
-            patch_ib(i)%vel(2) = 0._wp
-            patch_ib(i)%vel(3) = 0._wp
+            patch_ib(i)%vel(:) = 0._wp
+            patch_ib(i)%angles(:) = 0._wp
+            patch_ib(i)%angular_vel(:) = 0._wp
+
+            ! sets values of a rotation matrix which can be used when calculating rotations
+            patch_ib(i)%rotation_matrix = 0._wp
+            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

@@ -94,8 +94,9 @@ contains
         do i = 1, num_ibs
             if (patch_ib(i)%moving_ibm /= 0) then
                 moving_immersed_boundary_flag = .true.
-                exit
+
             end if
+            call s_update_ib_rotation_matrix(i)
         end do
 
         ! Allocating the patch identities bookkeeping variable
@@ -116,8 +117,8 @@ contains
         call s_mpi_allreduce_integer_sum(num_inner_gps, max_num_inner_gps)
 
         $:GPU_UPDATE(device='[num_gps, num_inner_gps]')
-        @:ALLOCATE(ghost_points(1:int(max_num_gps * 1.2)))
-        @:ALLOCATE(inner_points(1:int(max_num_inner_gps * 1.2)))
+        @:ALLOCATE(ghost_points(1:int(max_num_gps * 1.5)))
+        @:ALLOCATE(inner_points(1:int(max_num_inner_gps * 1.5)))
 
         $:GPU_ENTER_DATA(copyin='[ghost_points,inner_points]')
 
@@ -149,7 +150,7 @@ contains
         !!  @param q_prim_vf Primitive variables
         !!  @param pb Internal bubble pressure
         !!  @param mv Mass of vapor in bubble
-    pure subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in)
+    subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in)
 
         type(scalar_field), &
             dimension(sys_size), &
@@ -182,6 +183,8 @@ contains
         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) :: nbub
         real(wp) :: buf
@@ -265,9 +268,15 @@ contains
                     ! we know the object is not moving if moving_ibm is 0 (false)
                     vel_g = 0._wp
                 else
+                    ! get the vector that points from the centroid to the ghost
+                    radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, &
+                                                    patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid]
+                    ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear velocity
+                    rotation_velocity = cross_product(matmul(patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector)
                     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)
+                        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
@@ -896,26 +905,6 @@ contains
 
     end subroutine s_interpolate_image_point
 
-    !> Subroutine the updates the moving imersed boundary positions via Euler's method
-    impure subroutine s_propagate_mib(patch_id)
-
-        integer, intent(in) :: patch_id
-        integer :: i
-
-        ! start by using euler's method naiively, but eventually incorporate more sophistocation
-        if (patch_ib(patch_id)%moving_ibm == 1) then
-            ! this continues with euler's method, which is obviously not that great and we need to add acceleration
-            do i = 1, 3
-                patch_ib(patch_id)%vel(i) = patch_ib(patch_id)%vel(i) + 0.0*dt ! TODO :: ADD EXTERNAL FORCES HERE
-            end do
-
-            patch_ib(patch_id)%x_centroid = patch_ib(patch_id)%x_centroid + patch_ib(patch_id)%vel(1)*dt
-            patch_ib(patch_id)%y_centroid = patch_ib(patch_id)%y_centroid + patch_ib(patch_id)%vel(2)*dt
-            patch_ib(patch_id)%z_centroid = patch_ib(patch_id)%z_centroid + patch_ib(patch_id)%vel(3)*dt
-        end if
-
-    end subroutine s_propagate_mib
-
     !> 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)
@@ -929,10 +918,9 @@ contains
         ! Clears the existing immersed boundary indices
         ib_markers%sf = 0
 
+        ! recalulcate the rotation matrix based upon the new angles
         do i = 1, num_ibs
-            if (patch_ib(i)%moving_ibm /= 0) then
-                call s_propagate_mib(i)  ! TODO :: THIS IS DONE TERRIBLY WITH EULER METHOD
-            end if
+            if (patch_ib(i)%moving_ibm /= 0) call s_update_ib_rotation_matrix(i)
         end do
 
         ! recompute the new ib_patch locations and broadcast them.
@@ -963,4 +951,14 @@ contains
 
     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

src/simulation/m_start_up.fpp

@@ -1135,10 +1135,6 @@ contains
             end do
         end if
 
-        if (moving_immersed_boundary_flag) then
-            call s_update_mib(num_ibs, levelset, levelset_norm)
-        end if
-
         call s_compute_derived_variables(t_step)