Add rotating mibms

examples/2D_ibm_airfoil/case.py

@@ -94,12 +94,12 @@
             "patch_ib(1)%t": 0.15,
             "patch_ib(1)%p": 0.4,
             "patch_ib(1)%m": 0.02,
-            "patch_ib(1)%theta": 30,
+            "patch_ib(1)%angles(3)": -0.5235987756,  # 30 degrees clockwise rotation, in radians
             # Fluids Physical Parameters
             # Use the same stiffness as the air bubble
-            "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50(Not 1.40)
+            "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50 (Not 1.40)
             "fluid_pp(1)%pi_inf": 0,
-            "fluid_pp(2)%gamma": 1.0e00 / (gam_b - 1.0e00),  # 2.50(Not 1.40)
+            "fluid_pp(2)%gamma": 1.0e00 / (gam_b - 1.0e00),  # 2.50 (Not 1.40)
             "fluid_pp(2)%pi_inf": 0,
         }
     )

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/pre_process/m_data_output.fpp

@@ -550,8 +550,6 @@ contains
                         write (2, FMT) airfoil_grid_l(j)%x, airfoil_grid_l(j)%y
                     end do
                     close (2)
-
-                    print *, "Np", Np
                 end if
             end do
         end if

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/pre_process/m_mpi_proxy.fpp

@@ -80,6 +80,10 @@ contains
                 call MPI_BCAST(patch_bc(i)%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 
+            #:for VAR in ['vel', 'angular_vel', 'angles']
+                call MPI_BCAST(patch_ib(i)%${VAR}$, 3, mpi_p, 0, MPI_COMM_WORLD, ierr)
+            #:endfor
+
             call MPI_BCAST(patch_bc(i)%radius, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
 
             #:for VAR in ['centroid', 'length']

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
@@ -191,8 +194,8 @@ contains
         $:GPU_PARALLEL_LOOP(private='[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, &
-            & j,k,l,q]')
+            & gamma,pi_inf,Re_K,G_K,Gs,gp,innerp,norm,buf, radial_vector, &
+            & rotation_velocity, j,k,l,q]')
         do i = 1, num_gps
 
             gp = ghost_points(i)
@@ -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
@@ -381,7 +390,7 @@ contains
         type(ghost_point) :: gp
 
         integer :: q, dim !< Iterator variables
-        integer :: i, j, k !< Location indexes
+        integer :: i, j, k, l !< Location indexes
         integer :: patch_id !< IB Patch ID
         integer :: dir
         integer :: index
@@ -436,7 +445,6 @@ contains
                                .or. temp_loc > s_cc(index + 1)))
                         index = index + dir
                         if (index < -buff_size .or. index > bound) then
-                            print *, q, index, bound, buff_size
                             print *, "temp_loc=", temp_loc, " s_cc(index)=", s_cc(index), " s_cc(index+1)=", s_cc(index + 1)
                             print *, "Increase buff_size further in m_helper_basic (currently set to a minimum of 10)"
                             error stop "Increase buff_size"
@@ -896,26 +904,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)
@@ -924,14 +912,15 @@ contains
         type(levelset_field), intent(inout) :: levelset
         type(levelset_norm_field), intent(inout) :: levelset_norm
 
-        integer :: i
+        integer :: i, ierr
 
         ! 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
+                call s_update_ib_rotation_matrix(i)
             end if
         end do
 
@@ -963,4 +952,14 @@ contains
 
     end subroutine s_finalize_ibm_module
 
+    function cross_product(a, b) result(c)
+        implicit none
+        real(wp), intent(in) :: a(3), b(3)
+        real(wp) :: 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_mpi_proxy.fpp

@@ -203,11 +203,14 @@ contains
 
         do i = 1, num_ibs
             #:for VAR in [ 'radius', 'length_x', 'length_y', &
-                & 'x_centroid', 'y_centroid', 'c', 'm', 'p', 't', 'theta', 'slip', &
-                'moving_ibm', 'vel',]
+                & 'x_centroid', 'y_centroid', 'c', 'm', 'p', 't', 'theta', 'slip']
                 call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             #:endfor
-            call MPI_BCAST(patch_ib(i)%geometry, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            #:for VAR in ['vel', 'angular_vel', 'angles']
+                call MPI_BCAST(patch_ib(i)%${VAR}$, 3, mpi_p, 0, MPI_COMM_WORLD, ierr)
+            #:endfor
+            call MPI_BCAST(patch_ib(i)%geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(patch_ib(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
         end do
 
         do j = 1, num_probes_max

src/simulation/m_start_up.fpp

@@ -581,7 +581,6 @@ contains
                 if (patch_ib(i)%c > 0) then
                     Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1
                     allocate (MPI_IO_airfoil_IB_DATA%var(1:2*Np))
-                    print *, "HERE Np", Np
                 end if
             end do
         end if
@@ -943,8 +942,6 @@ contains
             do j = 1, num_ibs
                 if (patch_ib(j)%c > 0) then
 
-                    print *, "HERE Np", Np
-
                     allocate (airfoil_grid_u(1:Np))
                     allocate (airfoil_grid_l(1:Np))
 
@@ -1135,10 +1132,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)
 
 

src/simulation/m_time_steppers.fpp

@@ -566,6 +566,35 @@ contains
             if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf)
 
             if (ib) then
+                ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms
+                if (moving_immersed_boundary_flag) then
+                    do i = 1, num_ibs
+                        if (s == 1) then
+                            patch_ib(i)%step_vel = patch_ib(i)%vel
+                            patch_ib(i)%step_angular_vel = patch_ib(i)%angular_vel
+                            patch_ib(i)%step_angles = patch_ib(i)%angles
+                            patch_ib(i)%step_x_centroid = patch_ib(i)%x_centroid
+                            patch_ib(i)%step_y_centroid = patch_ib(i)%y_centroid
+                            patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid
+                        end if
+
+                        if (patch_ib(i)%moving_ibm == 1) then
+                            do j = 1, 3
+                                patch_ib(i)%vel(j) = (rk_coef(s, 1)*patch_ib(i)%step_vel(j) + rk_coef(s, 2)*patch_ib(i)%vel(j) + rk_coef(s, 3)*0._wp*dt)/rk_coef(s, 4) ! 0.0 is a placeholder for accelerations
+                                patch_ib(i)%angular_vel(j) = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel(j) + rk_coef(s, 2)*patch_ib(i)%angular_vel(j) + rk_coef(s, 3)*0._wp*dt)/rk_coef(s, 4)
+
+                                ! Update the angle of the IB
+                                patch_ib(i)%angles(j) = (rk_coef(s, 1)*patch_ib(i)%step_angles(j) + rk_coef(s, 2)*patch_ib(i)%angles(j) + rk_coef(s, 3)*patch_ib(i)%angular_vel(j)*dt)/rk_coef(s, 4)
+                            end do
+
+                            ! Update the position of the IB
+                            patch_ib(i)%x_centroid = (rk_coef(s, 1)*patch_ib(i)%step_x_centroid + rk_coef(s, 2)*patch_ib(i)%x_centroid + rk_coef(s, 3)*patch_ib(i)%vel(1)*dt)/rk_coef(s, 4)
+                            patch_ib(i)%y_centroid = (rk_coef(s, 1)*patch_ib(i)%step_y_centroid + rk_coef(s, 2)*patch_ib(i)%y_centroid + rk_coef(s, 3)*patch_ib(i)%vel(2)*dt)/rk_coef(s, 4)
+                            patch_ib(i)%z_centroid = (rk_coef(s, 1)*patch_ib(i)%step_z_centroid + rk_coef(s, 2)*patch_ib(i)%z_centroid + rk_coef(s, 3)*patch_ib(i)%vel(3)*dt)/rk_coef(s, 4)
+                        end if
+                    end do
+                    call s_update_mib(num_ibs, levelset, levelset_norm)
+                end if
                 if (qbmm .and. .not. polytropic) then
                     call s_ibm_correct_state(q_cons_ts(1)%vf, q_prim_vf, pb_ts(1)%sf, mv_ts(1)%sf)
                 else