Fix periodic ib issues

docs/documentation/case.md

@@ -197,7 +197,7 @@ MPI topology is automatically optimized to maximize the parallel efficiency for
 
 The Table lists the patch parameters.
 The parameters define the geometries and physical parameters of fluid components (patch) in the domain at initial condition.
-Note that the domain must be fully filled with patche(s).
+Note that the domain must be fully filled with patches.
 The code outputs error messages when an empty region is left in the domain.
 
 - `tau_e(i)` is the `i`-th component of the elastic stress tensor, ordered as `tau_xx`, `tau_xy`, `tau_yy`, `tau_xz`, `tau_yz`, and `tau_zz`. 1D simulation requires `tau(1)`, 2D `tau(1:3)`, and 3D `tau(1:6)`.

examples/3D_mibm_periodic_collision/case.py

@@ -0,0 +1,164 @@
+"""
+Case file to demonstrate and test the periodic collision of immersed boundaries.
+This file was generated by Conrad Delgado as a minimum reproducer to prevent regression
+of periodic immersed boundaries.
+"""
+
+import json
+import math
+
+import numpy as np
+
+gam_a = 1.4
+
+# particle diameter
+D = 0.1
+R = D / 2.0
+
+# particle params
+rho_s = 10.0
+vol_s = 4.0 / 3.0 * np.pi * R**3
+mass_s = rho_s * vol_s
+N_s = 2
+
+# fluid params
+M = 2.0
+Re = 500.0
+P = 101325
+rho = 1.225
+v1 = M * np.sqrt(gam_a * P / rho)
+mu = rho * v1 * D / Re
+
+# timestep
+dt = 1.0e-05
+Nt = 50
+t_save = 5
+
+# grid
+Nx = 31
+Ny = 63
+Nz = 31
+
+# immersed boundary dictionary
+ib_dict = {}
+ib_dict.update(
+    {
+        f"patch_ib({1})%geometry": 8,
+        f"patch_ib({1})%x_centroid": 0.0,
+        f"patch_ib({1})%y_centroid": -2.1 * D,
+        f"patch_ib({1})%z_centroid": 0.0,
+        f"patch_ib({1})%vel(2)": -100.0,
+        f"patch_ib({1})%radius": D / 2,
+        f"patch_ib({1})%slip": "F",
+        f"patch_ib({1})%moving_ibm": 2,
+        f"patch_ib({1})%mass": mass_s,
+        f"patch_ib({2})%geometry": 8,
+        f"patch_ib({2})%x_centroid": 0.0,
+        f"patch_ib({2})%y_centroid": 1.8 * D,
+        f"patch_ib({2})%z_centroid": 0.0,
+        f"patch_ib({2})%vel(2)": +100.0,
+        f"patch_ib({2})%radius": D / 2,
+        f"patch_ib({2})%slip": "F",
+        f"patch_ib({2})%moving_ibm": 2,
+        f"patch_ib({2})%mass": mass_s,
+    }
+)
+
+# Configuring case dictionary
+case_dict = {
+    # Logistics
+    "run_time_info": "T",
+    # Computational Domain Parameters
+    # x direction
+    "x_domain%beg": -1.25 * D,
+    "x_domain%end": 1.25 * D,
+    # y direction
+    "y_domain%beg": -2.5 * D,
+    "y_domain%end": 2.5 * D,
+    # z direction
+    "z_domain%beg": -1.25 * D,
+    "z_domain%end": 1.25 * D,
+    "cyl_coord": "F",
+    "m": Nx,
+    "n": Ny,
+    "p": Nz,
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": Nt,
+    "t_step_save": t_save,
+    # Simulation Algorithm Parameters
+    # Only one patches are necessary, the air tube
+    "num_patches": 1,
+    # Use the 5 equation model
+    "model_eqns": 2,
+    # 6 equations model does not need the K \div(u) term
+    "alt_soundspeed": "F",
+    # One fluids: air
+    "num_fluids": 1,
+    # time step
+    "mpp_lim": "F",
+    # Correct errors when computing speed of sound
+    "mixture_err": "T",
+    # Use TVD RK3 for time marching
+    "time_stepper": 3,
+    # Reconstruct the primitive variables to minimize spurious
+    # Use WENO5
+    "weno_order": 5,
+    "weno_eps": 1.0e-16,
+    "weno_Re_flux": "T",
+    "weno_avg": "T",
+    "avg_state": 2,
+    "null_weights": "F",
+    "mp_weno": "T",
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    # periodic bc
+    "bc_x%beg": -3,
+    "bc_x%end": -3,
+    "bc_y%beg": -1,
+    "bc_y%end": -1,
+    "bc_z%beg": -3,
+    "bc_z%end": -3,
+    # Set IB to True and add 1 patch
+    "ib": "T",
+    "num_ibs": N_s,
+    "viscous": "T",
+    # Formatted Database Files Structure Parameters
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "E_wrt": "T",
+    "parallel_io": "T",
+    "ib_state_wrt": "T",
+    "fd_order": 4,
+    # Patch: Constant Tube filled with air
+    # Specify the cylindrical air tube grid geometry
+    "patch_icpp(1)%geometry": 9,
+    "patch_icpp(1)%x_centroid": 0.0,
+    # Uniform medium density, centroid is at the center of the domain
+    "patch_icpp(1)%y_centroid": 0.0,
+    "patch_icpp(1)%z_centroid": 0.0,
+    "patch_icpp(1)%length_x": 2.5 * D,
+    "patch_icpp(1)%length_y": 5.0 * D,
+    "patch_icpp(1)%length_z": 2.5 * D,
+    # Specify the patch primitive variables
+    "patch_icpp(1)%vel(1)": 0.0e00,
+    "patch_icpp(1)%vel(2)": 0.0e00,
+    "patch_icpp(1)%vel(3)": 0.0e00,
+    "patch_icpp(1)%pres": P,
+    "patch_icpp(1)%alpha_rho(1)": rho,
+    "patch_icpp(1)%alpha(1)": 1.0e00,
+    # Patch: Sphere Immersed Boundary
+    # Fluids Physical Parameters
+    "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50(Not 1.40)
+    "fluid_pp(1)%pi_inf": 0,
+    "fluid_pp(1)%Re(1)": 1.0 / mu,
+    "collision_model": 1,  # soft-sphere collision model
+    "ib_coefficient_of_friction": 0.1,
+    "collision_time": 20.0 * dt,
+    "coefficient_of_restitution": 0.9,  # almost perfectly elastic
+}
+
+case_dict.update(ib_dict)
+
+print(json.dumps(case_dict))

src/common/m_global_parameters_common.fpp

@@ -109,8 +109,11 @@ module m_global_parameters_common
 
     !> @name Processor coordinates and parallel-IO addressing (identical declaration across all three targets)
     !> @{
-    integer, allocatable, dimension(:) :: proc_coords  !< Processor coordinates in MPI_CART_COMM
-    integer, allocatable, dimension(:) :: start_idx    !< Starting cell-center index of local processor in global grid
+    integer, allocatable, dimension(:) :: proc_coords      !< Processor coordinates in MPI_CART_COMM
+    integer, allocatable, dimension(:) :: start_idx        !< Starting cell-center index of local processor in global grid
+    integer                            :: num_procs_x = 1  !< Number of MPI ranks in x-direction
+    integer                            :: num_procs_y = 1  !< Number of MPI ranks in y-direction
+    integer                            :: num_procs_z = 1  !< Number of MPI ranks in z-direction
     !> @}
 
     !> @name MPI info for parallel IO with Lustre file systems (identical across all three targets)

src/common/m_mpi_common.fpp

@@ -1021,7 +1021,6 @@ contains
     subroutine s_mpi_decompose_computational_domain
 
 #ifdef MFC_MPI
-        integer :: num_procs_x, num_procs_y, num_procs_z  !< Optimal number of processors in the x-, y- and z-directions
         !> Non-optimal number of processors in the x-, y- and z-directions
         real(wp) :: tmp_num_procs_x, tmp_num_procs_y, tmp_num_procs_z
         real(wp) :: fct_min        !< Processor factorization (fct) minimization parameter

src/post_process/m_start_up.fpp

@@ -41,7 +41,6 @@ module m_start_up
     complex(c_double_complex), allocatable  :: data_cmplx(:,:,:), data_cmplx_y(:,:,:), data_cmplx_z(:,:,:)
     real(wp), allocatable, dimension(:,:,:) :: En_real
     real(wp), allocatable, dimension(:)     :: En
-    integer                                 :: num_procs_x, num_procs_y, num_procs_z
     integer                                 :: Nx, Ny, Nz, Nxloc, Nyloc, Nyloc2, Nzloc, Nf
     integer                                 :: ierr
     integer                                 :: MPI_COMM_CART, MPI_COMM_CART12, MPI_COMM_CART13

src/simulation/m_collisions.fpp

@@ -20,14 +20,17 @@ module m_collisions
     implicit none
 
     private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module, &
-        & f_local_rank_owns_location, f_neighborhood_ranks_own_location
+        & f_local_rank_owns_location, f_neighborhood_ranks_own_location, ib_gbl_idx_lookup
     ! overlap distances for computing collisions
     integer, allocatable, dimension(:,:)  :: collision_lookup
     real(wp), allocatable, dimension(:,:) :: wall_overlap_distances
     real(wp)                              :: spring_stiffness, damping_parameter
     $:GPU_DECLARE(create='[spring_stiffness, damping_parameter]')
     $:GPU_DECLARE(create='[collision_lookup, wall_overlap_distances]')
 
+    integer, dimension(:), allocatable :: ib_gbl_idx_lookup
+    $:GPU_DECLARE(create='[ib_gbl_idx_lookup]')
+
 contains
 
     subroutine s_initialize_collisions_module()
@@ -99,6 +102,9 @@ contains
             encoded_pid2 = collision_lookup(i, 4)
             call s_decode_patch_periodicity(encoded_pid1, pid1, xp1, yp1, zp1)
             call s_decode_patch_periodicity(encoded_pid2, pid2, xp2, yp2, zp2)
+            pid1 = collision_lookup(i, 1)
+            pid2 = collision_lookup(i, 2)
+
             ! call s_get_neighborhood_idx(pid1, pid1) ! global patch ID -> local index call s_get_neighborhood_idx(pid2, pid2)
             if (pid1 <= 0 .or. pid2 <= 0) cycle
 
@@ -291,6 +297,8 @@ contains
             ! get the decoded pairs for checking if they exist, using ii,jj,kk as dummy indices
             call s_decode_patch_periodicity(raw_pairs(pair_idx, 1), decoded_pairs(1), ii, jj, kk)
             call s_decode_patch_periodicity(raw_pairs(pair_idx, 2), decoded_pairs(2), ii, jj, kk)
+            decoded_pairs(1) = ib_gbl_idx_lookup(decoded_pairs(1))
+            decoded_pairs(2) = ib_gbl_idx_lookup(decoded_pairs(2))
 
             ! skip self-collisions (an IB cannot collide with its own periodic image)
             if (decoded_pairs(1) == decoded_pairs(2)) cycle
@@ -329,34 +337,46 @@ contains
     subroutine s_detect_ib_collisions_n2(num_considered_collisions)
 
         integer, intent(out)   :: num_considered_collisions
-        integer                :: pid1, pid2, encoded_pid1, encoded_pid2, current_collisions
+        integer                :: pid1, pid2, encoded_pid2, current_collisions
+        integer                :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper, xp, yp, zp
         real(wp), dimension(3) :: centroid_1, centroid_2, distance_vec
 
         num_considered_collisions = 0
 
-        $:GPU_PARALLEL_LOOP(private='[pid1, pid2, centroid_1, centroid_2, distance_vec, current_collisions]', &
-                            & copy='[num_considered_collisions]')
+        call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
+
+        $:GPU_PARALLEL_LOOP(private='[pid1, pid2, encoded_pid2, centroid_1, centroid_2, xp, yp, zp, distance_vec, &
+                            & current_collisions]', copyin='[xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper]', copy='[num_considered_collisions]')
         do pid1 = 1, num_ibs - 1
+            centroid_1 = [patch_ib(pid1)%x_centroid, patch_ib(pid1)%y_centroid, 0._wp]
+            if (num_dims == 3) centroid_1(3) = patch_ib(pid1)%z_centroid
             do pid2 = pid1 + 1, num_ibs
-                centroid_1 = [patch_ib(pid1)%x_centroid, patch_ib(pid1)%y_centroid, 0._wp]
-                centroid_2 = [patch_ib(pid2)%x_centroid, patch_ib(pid2)%y_centroid, 0._wp]
-                if (num_dims == 3) then
-                    centroid_1(3) = patch_ib(pid1)%z_centroid
-                    centroid_2(3) = patch_ib(pid2)%z_centroid
-                end if
-                distance_vec = centroid_2 - centroid_1
-
-                if (norm2(distance_vec) < patch_ib(pid1)%radius + patch_ib(pid2)%radius) then
-                    $:GPU_ATOMIC(atomic='capture')
-                    num_considered_collisions = num_considered_collisions + 1
-                    current_collisions = num_considered_collisions
-                    $:END_GPU_ATOMIC_CAPTURE()
-
-                    collision_lookup(current_collisions, 1) = pid1
-                    collision_lookup(current_collisions, 2) = pid2
-                    collision_lookup(current_collisions, 3) = pid1
-                    collision_lookup(current_collisions, 4) = pid2
-                end if
+                periodic_search: do xp = xp_lower, xp_upper
+                    do yp = yp_lower, yp_upper
+                        do zp = zp_lower, zp_upper
+                            centroid_2(1) = patch_ib(pid2)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
+                            centroid_2(2) = patch_ib(pid2)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
+                            if (num_dims == 3) centroid_2(3) = patch_ib(pid2)%z_centroid + real(zp, &
+                                & wp)*(z_domain%end - z_domain%beg)
+                            distance_vec = centroid_2 - centroid_1
+
+                            if (norm2(distance_vec) < patch_ib(pid1)%radius + patch_ib(pid2)%radius) then
+                                $:GPU_ATOMIC(atomic='capture')
+                                num_considered_collisions = num_considered_collisions + 1
+                                current_collisions = num_considered_collisions
+                                $:END_GPU_ATOMIC_CAPTURE()
+
+                                call s_encode_patch_periodicity(patch_ib(pid2)%gbl_patch_id, xp, yp, zp, encoded_pid2)
+
+                                collision_lookup(current_collisions, 1) = pid1
+                                collision_lookup(current_collisions, 2) = pid2
+                                collision_lookup(current_collisions, 3) = patch_ib(pid1)%gbl_patch_id
+                                collision_lookup(current_collisions, 4) = encoded_pid2
+                                exit periodic_search
+                            end if
+                        end do
+                    end do
+                end do periodic_search
             end do
         end do
         $:END_GPU_PARALLEL_LOOP()

src/simulation/m_derived_variables.fpp

@@ -21,8 +21,6 @@ module m_derived_variables
     private; public :: s_initialize_derived_variables_module, s_initialize_derived_variables, s_compute_derived_variables, &
         & s_finalize_derived_variables_module
 
-    ! fd_coeff_x, fd_coeff_y, fd_coeff_z: declared in m_global_parameters so m_viscous can see them too
-
     ! @name Variables for computing acceleration
     !> @{
     real(wp), public, allocatable, dimension(:,:,:) :: accel_mag

src/simulation/m_ib_patches.fpp

@@ -23,7 +23,7 @@ module m_ib_patches
     implicit none
 
     private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_STL_models, s_decode_patch_periodicity, &
-        & s_initialize_ib_airfoils
+        & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities
 
 contains
 
@@ -48,8 +48,9 @@ contains
         real(wp), dimension(3) :: center, xyz_local, length
         real(wp) :: radius, eta
 
-        !  3D Patch Geometries
+        $:GPU_UPDATE(host='[patch_ib(1:num_ibs)]')
 
+        !  3D Patch Geometries
         if (num_dims == 3) then
             call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
             do xp = xp_lower, xp_upper
@@ -678,14 +679,17 @@ contains
         #:endfor
 
         ! z only if 3D
-        if (present(zp_lower) .and. p /= 0) then
+        if (present(zp_lower) .and. num_dims == 3) then
             if (ib_bc_z%beg == BC_PERIODIC) then
                 zp_lower = -1
                 zp_upper = 1
             else
                 zp_lower = 0
                 zp_upper = 0
             end if
+        else if (present(zp_lower)) then
+            zp_lower = 0
+            zp_upper = 0
         end if
 
     end subroutine s_get_periodicities