Fix periodic ib issues

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,35 +337,47 @@ 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
-            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
-            end do
+            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
+            collision_partner_loop: do pid2 = pid1 + 1, num_ibs
+                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
+                                cycle collision_partner_loop
+                            end if
+                        end do
+                    end do
+                end do
+            end do collision_partner_loop
         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