periodic BCs mostly working

Author: wilfonba

src/common/m_boundary_common.fpp

@@ -1942,6 +1942,21 @@ contains
         offset_x%beg = buff_size; offset_x%end = buff_size
         offset_y%beg = buff_size; offset_y%end = buff_size
         offset_z%beg = buff_size; offset_z%end = buff_size
+
+#ifdef MFC_MPI
+        ! Populate global domain boundaries with stretched grids
+        call s_mpi_allreduce_min(x_cb(-1), glb_bounds(1)%beg)
+        call s_mpi_allreduce_max(x_cb(m), glb_bounds(1)%end)
+
+        if (n > 0) then
+            call s_mpi_allreduce_min(y_cb(-1), glb_bounds(2)%beg)
+            call s_mpi_allreduce_max(y_cb(n), glb_bounds(2)%end)
+            if (p > 0) then
+                call s_mpi_allreduce_min(z_cb(-1), glb_bounds(3)%beg)
+                call s_mpi_allreduce_max(z_cb(p), glb_bounds(3)%end)
+            end if
+        end if
+#endif
 #endif
 
 #ifndef MFC_PRE_PROCESS

src/simulation/m_bubbles_EL.fpp

@@ -79,8 +79,11 @@ module m_bubbles_EL
     $:GPU_DECLARE(create='[nBubs,Rmax_glb,Rmin_glb,q_beta,q_beta_idx]')
 
     integer, allocatable, dimension(:) :: keep_bubble, prefix_sum
+    integer, allocatable, dimension(:) :: wrap_bubble_loc
+    integer, allocatable, dimension(:,:) :: wrap_bubble_dir
     integer :: active_bubs
     $:GPU_DECLARE(create='[keep_bubble, prefix_sum, active_bubs]')
+    $:GPU_DECLARE(create='[wrap_bubble_loc, wrap_bubble_dir]')
 
 contains
 
@@ -110,6 +113,17 @@ contains
             call s_mpi_abort('Please check the lag_params%solver_approach input')
         end if
 
+        pcomm_coords(1)%beg = x_cb(buff_size - fd_number - 1)
+        pcomm_coords(1)%end = x_cb(m - buff_size + fd_number)
+        if (n > 0) then
+            pcomm_coords(2)%beg = y_cb(buff_size - fd_number - 1)
+            pcomm_coords(2)%end = y_cb(n - buff_size + fd_number)
+            if (p > 0) then
+                pcomm_coords(3)%beg = z_cb(buff_size - fd_number - 1)
+                pcomm_coords(3)%end = z_cb(p - buff_size + fd_number)
+            end if
+        end if
+
         $:GPU_UPDATE(device='[lag_num_ts, q_beta_idx]')
 
         @:ALLOCATE(q_beta%vf(1:q_beta_idx))
@@ -149,6 +163,7 @@ contains
         @:ALLOCATE(mtn_dveldt(1:nBubs_glb, 1:3, 1:lag_num_ts))
 
         @:ALLOCATE(keep_bubble(1:nBubs_glb), prefix_sum(1:nBubs_glb))
+        @:ALLOCATE(wrap_bubble_loc(1:nBubs_glb), wrap_bubble_dir(1:nBubs_glb, 1:num_dims))
 
         if (adap_dt .and. f_is_default(adap_dt_tol)) adap_dt_tol = dflt_adap_dt_tol
 
@@ -274,7 +289,8 @@ contains
                                           gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, &
                                           intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, &
                                           mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, &
-                                          gas_dmvdt, mtn_dposdt, mtn_dveldt, lag_num_ts, nBubs)
+                                          gas_dmvdt, mtn_dposdt, mtn_dveldt, lag_num_ts, nBubs, &
+                                          dest=1)
         end if
 
         $:GPU_UPDATE(device='[bubbles_lagrange, lag_params]')
@@ -576,8 +592,6 @@ contains
             call MPI_TYPE_FREE(view, ierr)
         end if
 
-        print *, proc_rank, lag_id
-
         if (proc_rank == 0) then
             write (*, '(A,I0,A,I0)') 'Read ', file_tot_part, ' particles from restart file at t_step = ', save_count
             write (*, '(A,E15.7,A,E15.7)') 'Restart time = ', mytime, ', dt = ', dt
@@ -1358,12 +1372,16 @@ contains
 
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
         integer, intent(in) :: dest
-        integer :: k, i, patch_id, offset
+        integer :: k, i, q
+        integer ::  patch_id, newBubs, new_idx
+        real(wp) :: offset
         integer, dimension(3) :: cell
 
         $:GPU_PARALLEL_LOOP(private='[cell]')
         do k = 1, nBubs
             keep_bubble(k) = 1
+            wrap_bubble_loc(k) = 0
+            wrap_bubble_dir(k,:) = 0
 
             ! Relocate bubbles at solid boundaries and delete bubbles that leave
             ! buffer regions
@@ -1373,6 +1391,16 @@ contains
             elseif (any(bc_x%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) &
                     .and. mtn_pos(k, 1, dest) > x_cb(m) - intfc_rad(k, dest)) then
                 mtn_pos(k, 1, dest) = x_cb(m) - intfc_rad(k, dest)
+            elseif (bc_x%beg == BC_PERIODIC .and. mtn_pos(k, 1, dest) < pcomm_coords(1)%beg .and. &
+                    mtn_posPrev(k, 1, dest) > pcomm_coords(1)%beg) then
+                !print*, "1, -1", mtn_pos(k, 1, dest), mtn_posPrev(k, 1, dest), pcomm_coords(1)%beg
+                wrap_bubble_dir(k,1) = 1
+                wrap_bubble_loc(k) = -1
+            elseif (bc_x%end == BC_PERIODIC .and. mtn_pos(k, 1, dest) > pcomm_coords(1)%end .and. &
+                    mtn_posPrev(k, 1, dest) < pcomm_coords(1)%end) then
+                !print*, "1, 1", mtn_pos(k, 1, dest), mtn_posPrev(k, 1, dest), pcomm_coords(1)%end
+                wrap_bubble_dir(k,1) = 1
+                wrap_bubble_loc(k) = 1
             elseif (mtn_pos(k, 1, dest) >= x_cb(m + buff_size - fd_number)) then
                 keep_bubble(k) = 0
             elseif (mtn_pos(k, 1, dest) < x_cb(fd_number - buff_size - 1)) then
@@ -1385,6 +1413,16 @@ contains
             else if (any(bc_y%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) &
                      .and. mtn_pos(k, 2, dest) > y_cb(n) - intfc_rad(k, dest)) then
                 mtn_pos(k, 2, dest) = y_cb(n) - intfc_rad(k, dest)
+            elseif (bc_y%beg == BC_PERIODIC .and. mtn_pos(k, 2, dest) < pcomm_coords(2)%beg .and. &
+                    mtn_posPrev(k, 2, dest) > pcomm_coords(2)%beg) then
+                !print*, "2, -1", mtn_pos(k, 2, dest), mtn_posPrev(k, 2, dest), pcomm_coords(2)%beg
+                wrap_bubble_dir(k,2) = 1
+                wrap_bubble_loc(k) = -1
+            elseif (bc_y%end == BC_PERIODIC .and. mtn_pos(k, 2, dest) > pcomm_coords(2)%end .and. &
+                    mtn_posPrev(k, 2, dest) < pcomm_coords(2)%end) then
+                !print*, "2, 1", mtn_pos(k, 2, dest), mtn_posPrev(k, 2, dest), pcomm_coords(2)%end
+                wrap_bubble_dir(k,2) = 1
+                wrap_bubble_loc(k) = 1
             elseif (mtn_pos(k, 2, dest) >= y_cb(n + buff_size - fd_number)) then
                 keep_bubble(k) = 0
             elseif (mtn_pos(k, 2, dest) < y_cb(fd_number - buff_size - 1)) then
@@ -1398,6 +1436,14 @@ contains
                 else if (any(bc_z%end == (/BC_REFLECTIVE, BC_CHAR_SLIP_WALL, BC_SLIP_WALL, BC_NO_SLIP_WALL/)) &
                          .and. mtn_pos(k, 3, dest) > z_cb(p) - intfc_rad(k, dest)) then
                     mtn_pos(k, 3, dest) = z_cb(p) - intfc_rad(k, dest)
+                elseif (bc_z%beg == BC_PERIODIC .and. mtn_pos(k, 3, dest) < pcomm_coords(3)%beg .and. &
+                        mtn_posPrev(k, 3, dest) > pcomm_coords(3)%beg) then
+                    wrap_bubble_dir(k,3) = 1
+                    wrap_bubble_loc(k) = -1
+                elseif (bc_z%end == BC_PERIODIC .and. mtn_pos(k, 3, dest) > pcomm_coords(3)%end .and. &
+                        mtn_posPrev(k, 3, dest) < pcomm_coords(3)%end) then
+                    wrap_bubble_dir(k,3) = 1
+                    wrap_bubble_loc(k) = 1
                 elseif (mtn_pos(k, 3, dest) >= z_cb(p + buff_size - fd_number)) then
                     keep_bubble(k) = 0
                 elseif (mtn_pos(k, 3, dest) < z_cb(fd_number - buff_size - 1)) then
@@ -1436,6 +1482,7 @@ contains
             end if
         end do
 
+        call nvtxStartRange("LAG-BC")
         call nvtxStartRange("LAG-BC-DEV2HOST")
         $:GPU_UPDATE(host='[bub_R0, Rmax_stats, Rmin_stats, gas_mg, gas_betaT, &
             & gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, intfc_rad, intfc_vel, &
@@ -1444,6 +1491,7 @@ contains
             & nBubs]')
         call nvtxEndRange
 
+        ! Handle deletion of bubbles leaving local domain
         do k = 1, nBubs
             if (k == 1) then
                 prefix_sum(k) = keep_bubble(k)
@@ -1456,33 +1504,41 @@ contains
         do k = 1, nBubs
             if (keep_bubble(k) == 1) then
                 if (prefix_sum(k) /= k) then
-                    bub_R0(prefix_sum(k)) = bub_R0(k)
-                    Rmax_stats(prefix_sum(k)) = Rmax_stats(k)
-                    Rmin_stats(prefix_sum(k)) = Rmin_stats(k)
-                    gas_mg(prefix_sum(k)) = gas_mg(k)
-                    gas_betaT(prefix_sum(k)) = gas_betaT(k)
-                    gas_betaC(prefix_sum(k)) = gas_betaC(k)
-                    bub_dphidt(prefix_sum(k)) = bub_dphidt(k)
-                    lag_id(prefix_sum(k), 1) = lag_id(k, 1)
-                    gas_p(prefix_sum(k), 1:2) = gas_p(k, 1:2)
-                    gas_mv(prefix_sum(k), 1:2) = gas_mv(k, 1:2)
-                    intfc_rad(prefix_sum(k), 1:2) = intfc_rad(k, 1:2)
-                    intfc_vel(prefix_sum(k), 1:2) = intfc_vel(k, 1:2)
-                    mtn_pos(prefix_sum(k), 1:3, 1:2) = mtn_pos(k, 1:3, 1:2)
-                    mtn_posPrev(prefix_sum(k), 1:3, 1:2) = mtn_posPrev(k, 1:3, 1:2)
-                    mtn_vel(prefix_sum(k), 1:3, 1:2) = mtn_vel(k, 1:3, 1:2)
-                    mtn_s(prefix_sum(k), 1:3, 1:2) = mtn_s(k, 1:3, 1:2)
-                    intfc_draddt(prefix_sum(k), 1:lag_num_ts) = intfc_draddt(k, 1:lag_num_ts)
-                    intfc_dveldt(prefix_sum(k), 1:lag_num_ts) = intfc_dveldt(k, 1:lag_num_ts)
-                    gas_dpdt(prefix_sum(k), 1:lag_num_ts) = gas_dpdt(k, 1:lag_num_ts)
-                    gas_dmvdt(prefix_sum(k), 1:lag_num_ts) = gas_dmvdt(k, 1:lag_num_ts)
-                    mtn_dposdt(prefix_sum(k), 1:3, 1:lag_num_ts) = mtn_dposdt(k, 1:3, 1:lag_num_ts)
-                    mtn_dveldt(prefix_sum(k), 1:3, 1:lag_num_ts) = mtn_dveldt(k, 1:3, 1:lag_num_ts)
+                    call s_copy_lag_bubble(k, prefix_sum(k))
                 end if
             end if
             if (k == nBubs) nBubs = active_bubs
         end do
 
+        ! Handle periodic wrapping of bubbles on same processor
+        newBubs = 0
+        do k = 1, nBubs
+            if (any(wrap_bubble_dir(k, :) == 1)) then
+                newBubs = newBubs + 1
+                new_idx = nBubs + newBubs
+                call s_copy_lag_bubble(k, new_idx)
+                do i = 1, num_dims
+                    if (wrap_bubble_dir(k, i) == 1) then
+                        offset = glb_bounds(i)%end - glb_bounds(i)%beg
+                        if (wrap_bubble_loc(k) == 1) then
+                            do q = 1, 2
+                                mtn_pos(new_idx, i, q) = mtn_pos(new_idx, i, q) - offset
+                                mtn_posPrev(new_idx, i, q) = mtn_posPrev(new_idx, i, q) - offset
+                            end do
+                        else if (wrap_bubble_loc(k) == -1) then
+                            do q = 1, 2
+                                mtn_pos(new_idx, i, q) = mtn_pos(new_idx, i, q) + offset
+                                mtn_posPrev(new_idx, i, q) = mtn_posPrev(new_idx, i, q) + offset
+                            end do
+                        end if
+                    end if
+                end do
+            end if
+        end do
+        print*, nBubs, newBubs
+        nBubs = nBubs + newBubs
+
+        ! Handle MPI transfer of bubbles going to another processor's local domain
         if (num_procs > 1) then
             call nvtxStartRange("LAG-BC-TRANSFER-LIST")
             call s_add_particles_to_transfer_list(nBubs, mtn_pos(:, :, dest), mtn_posPrev(:, :, dest))
@@ -1493,7 +1549,8 @@ contains
                                           gas_betaC, bub_dphidt, lag_id, gas_p, gas_mv, &
                                           intfc_rad, intfc_vel, mtn_pos, mtn_posPrev, mtn_vel, &
                                           mtn_s, intfc_draddt, intfc_dveldt, gas_dpdt, &
-                                          gas_dmvdt, mtn_dposdt, mtn_dveldt, lag_num_ts, nBubs)
+                                          gas_dmvdt, mtn_dposdt, mtn_dveldt, lag_num_ts, nBubs, &
+                                          dest)
             call nvtxEndRange
         end if
 
@@ -1503,6 +1560,7 @@ contains
             & mtn_pos, mtn_posPrev, mtn_vel, mtn_s, intfc_draddt, intfc_dveldt, &
             & gas_dpdt, gas_dmvdt, mtn_dposdt, mtn_dveldt, nBubs]')
         call nvtxEndRange
+        call nvtxEndRange
 
     end subroutine s_enforce_EL_bubbles_boundary_conditions
 
@@ -2034,42 +2092,34 @@ contains
 
     !> The purpose of this subroutine is to remove one specific particle if dt is too small.
           !! @param bub_id Particle id
-    impure subroutine s_remove_lag_bubble(bub_id)
-
-        integer, intent(in) :: bub_id
-
-        integer :: i
-
-        $:GPU_LOOP(parallelism='[seq]')
-        do i = bub_id, nBubs - 1
-            bub_R0(i) = bub_R0(i + 1)
-            Rmax_stats(i) = Rmax_stats(i + 1)
-            Rmin_stats(i) = Rmin_stats(i + 1)
-            gas_mg(i) = gas_mg(i + 1)
-            gas_betaT(i) = gas_betaT(i + 1)
-            gas_betaC(i) = gas_betaC(i + 1)
-            bub_dphidt(i) = bub_dphidt(i + 1)
-            lag_id(i, 1) = lag_id(i + 1, 1)
-            gas_p(i, 1:2) = gas_p(i + 1, 1:2)
-            gas_mv(i, 1:2) = gas_mv(i + 1, 1:2)
-            intfc_rad(i, 1:2) = intfc_rad(i + 1, 1:2)
-            intfc_vel(i, 1:2) = intfc_vel(i + 1, 1:2)
-            mtn_pos(i, 1:3, 1:2) = mtn_pos(i + 1, 1:3, 1:2)
-            mtn_posPrev(i, 1:3, 1:2) = mtn_posPrev(i + 1, 1:3, 1:2)
-            mtn_vel(i, 1:3, 1:2) = mtn_vel(i + 1, 1:3, 1:2)
-            mtn_s(i, 1:3, 1:2) = mtn_s(i + 1, 1:3, 1:2)
-            intfc_draddt(i, 1:lag_num_ts) = intfc_draddt(i + 1, 1:lag_num_ts)
-            intfc_dveldt(i, 1:lag_num_ts) = intfc_dveldt(i + 1, 1:lag_num_ts)
-            gas_dpdt(i, 1:lag_num_ts) = gas_dpdt(i + 1, 1:lag_num_ts)
-            gas_dmvdt(i, 1:lag_num_ts) = gas_dmvdt(i + 1, 1:lag_num_ts)
-            mtn_dposdt(i, 1:3, 1:lag_num_ts) = mtn_dposdt(i + 1, 1:3, 1:lag_num_ts)
-            mtn_dveldt(i, 1:3, 1:lag_num_ts) = mtn_dveldt(i + 1, 1:3, 1:lag_num_ts)
-        end do
-
-        nBubs = nBubs - 1
-        $:GPU_UPDATE(device='[nBubs]')
-
-    end subroutine s_remove_lag_bubble
+    impure subroutine s_copy_lag_bubble(src, dest)
+
+        integer, intent(in) :: src, dest
+
+        bub_R0(dest) = bub_R0(src)
+        Rmax_stats(dest) = Rmax_stats(src)
+        Rmin_stats(dest) = Rmin_stats(src)
+        gas_mg(dest) = gas_mg(src)
+        gas_betaT(dest) = gas_betaT(src)
+        gas_betaC(dest) = gas_betaC(src)
+        bub_dphidt(dest) = bub_dphidt(src)
+        lag_id(dest, 1) = lag_id(src, 1)
+        gas_p(dest, 1:2) = gas_p(src, 1:2)
+        gas_mv(dest, 1:2) = gas_mv(src, 1:2)
+        intfc_rad(dest, 1:2) = intfc_rad(src, 1:2)
+        intfc_vel(dest, 1:2) = intfc_vel(src, 1:2)
+        mtn_vel(dest, 1:3, 1:2) = mtn_vel(src, 1:3, 1:2)
+        mtn_s(dest, 1:3, 1:2) = mtn_s(src, 1:3, 1:2)
+        intfc_draddt(dest, 1:lag_num_ts) = intfc_draddt(src, 1:lag_num_ts)
+        intfc_dveldt(dest, 1:lag_num_ts) = intfc_dveldt(src, 1:lag_num_ts)
+        gas_dpdt(dest, 1:lag_num_ts) = gas_dpdt(src, 1:lag_num_ts)
+        gas_dmvdt(dest, 1:lag_num_ts) = gas_dmvdt(src, 1:lag_num_ts)
+        mtn_dposdt(dest, 1:3, 1:lag_num_ts) = mtn_dposdt(src, 1:3, 1:lag_num_ts)
+        mtn_dveldt(dest, 1:3, 1:lag_num_ts) = mtn_dveldt(src, 1:3, 1:lag_num_ts)
+        mtn_pos(dest, 1:3, 1:2) = mtn_pos(src, 1:3, 1:2)
+        mtn_posPrev(dest, 1:3, 1:2) = mtn_posPrev(src, 1:3, 1:2)
+
+    end subroutine s_copy_lag_bubble
 
     !> The purpose of this subroutine is to deallocate variables
     impure subroutine s_finalize_lagrangian_solver()

src/simulation/m_global_parameters.fpp

@@ -63,14 +63,14 @@ module m_global_parameters
     !> @{
 
     real(wp), target, allocatable, dimension(:) :: x_cb, y_cb, z_cb
+    type(bounds_info), dimension(3) :: glb_bounds !<
     !> @}
 
     !> @name Cell-center (CC) locations in the x-, y- and z-directions, respectively
     !> @{
 
     real(wp), target, allocatable, dimension(:) :: x_cc, y_cc, z_cc
     !> @}
-    !type(bounds_info) :: x_domain, y_domain, z_domain !<
     !! Locations of the domain bounds in the x-, y- and z-coordinate directions
     !> @name Cell-width distributions in the x-, y- and z-directions, respectively
     !> @{
@@ -231,6 +231,7 @@ module m_global_parameters
 
     integer :: num_bc_patches
     logical :: bc_io
+    logical, dimension(3) :: periodic_bc
     !> @name Boundary conditions (BC) in the x-, y- and z-directions, respectively
     !> @{
     type(int_bounds_info) :: bc_x, bc_y, bc_z
@@ -239,10 +240,6 @@ module m_global_parameters
     $:GPU_DECLARE(create='[bc_y%vb1, bc_y%vb2, bc_y%vb3, bc_y%ve1, bc_y%ve2, bc_y%ve3]')
     $:GPU_DECLARE(create='[bc_z%vb1, bc_z%vb2, bc_z%vb3, bc_z%ve1, bc_z%ve2, bc_z%ve3]')
 
-    type(bounds_info) :: x_domain, y_domain, z_domain
-    real(wp) :: x_a, y_a, z_a
-    real(wp) :: x_b, y_b, z_b
-
     logical :: parallel_io !< Format of the data files
     logical :: file_per_process !< shared file or not when using parallel io
     integer :: precision !< Precision of output files
@@ -252,6 +249,9 @@ module m_global_parameters
     integer, allocatable, dimension(:) :: proc_coords !<
     !! Processor coordinates in MPI_CART_COMM
 
+    type(bounds_info), allocatable, dimension(:) :: pcomm_coords
+    !! Coordinates for EL particle transfer
+
     type(int_bounds_info), dimension(3) :: nidx !< Indices for neighboring processors
 
     integer, allocatable, dimension(:, :, :) :: neighbor_ranks
@@ -652,6 +652,7 @@ contains
 
         num_bc_patches = 0
         bc_io = .false.
+        periodic_bc = .false.
 
         bc_x%beg = dflt_int; bc_x%end = dflt_int
         bc_y%beg = dflt_int; bc_y%end = dflt_int
@@ -664,9 +665,9 @@ contains
             #:endfor
         #:endfor
 
-        x_domain%beg = dflt_int; x_domain%end = dflt_int
-        y_domain%beg = dflt_int; y_domain%end = dflt_int
-        z_domain%beg = dflt_int; z_domain%end = dflt_int
+        glb_bounds(1)%beg = dflt_int; glb_bounds(1)%end = dflt_int
+        glb_bounds(2)%beg = dflt_int; glb_bounds(2)%end = dflt_int
+        glb_bounds(3)%beg = dflt_int; glb_bounds(3)%end = dflt_int
 
         ! Fluids physical parameters
         do i = 1, num_fluids_max
@@ -1388,6 +1389,7 @@ contains
         #:endif
 
         allocate (proc_coords(1:num_dims))
+        allocate (pcomm_coords(1:num_dims))
 
         if (parallel_io .neqv. .true.) return
 
@@ -1424,6 +1426,7 @@ contains
         end if
 
         deallocate (proc_coords)
+        deallocate (pcomm_coords)
         if (parallel_io) then
             deallocate (start_idx)