Skip to content

Commit 6446782

Browse files
author
Daniel Vickers
committed
Successfully ran 128k IBs
1 parent 0f79876 commit 6446782

3 files changed

Lines changed: 105 additions & 87 deletions

File tree

‎src/common/m_constants.fpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ module m_constants
2424
integer, parameter :: num_fluids_max = 10 !< Maximum number of fluids in the simulation
2525
integer, parameter :: num_probes_max = 10 !< Maximum number of flow probes in the simulation
2626
integer, parameter :: num_patches_max = 10 !< Maximum number of IC patches
27-
integer, parameter :: num_ib_patches_max = 50000 !< Maximum number of immersed boundary patches (patch_ib)
27+
integer, parameter :: num_ib_patches_max = 300000 !< Maximum number of immersed boundary patches (patch_ib)
2828
integer, parameter :: num_local_ibs_max = 2000 !< Maximum number of immersed boundary patches (patch_ib)
2929
integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches
3030
integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13)

‎src/simulation/m_collisions.fpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,8 @@ contains
6262

6363
! get is distance used in the force calculation with each IB and each wall
6464
call s_detect_wall_collisions()
65-
call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions)
65+
! call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions)
66+
call s_detect_ib_collisions_n2(num_considered_collisions)
6667

6768
select case (collision_model)
6869
case (1) ! soft sphere model
@@ -98,6 +99,8 @@ contains
9899
encoded_pid2 = collision_lookup(i, 4)
99100
call s_decode_patch_periodicity(encoded_pid1, pid1, xp1, yp1, zp1)
100101
call s_decode_patch_periodicity(encoded_pid2, pid2, xp2, yp2, zp2)
102+
! call s_get_neighborhood_idx(pid1, pid1) ! global patch ID -> local index call s_get_neighborhood_idx(pid2, pid2)
103+
if (pid1 <= 0 .or. pid2 <= 0) cycle
101104

102105
centroid_1(1) = patch_ib(pid1)%x_centroid + real(xp1, wp)*(x_domain%end - x_domain%beg)
103106
centroid_1(2) = patch_ib(pid1)%y_centroid + real(yp1, wp)*(y_domain%end - y_domain%beg)

‎src/simulation/m_ibm.fpp‎

Lines changed: 100 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -171,17 +171,21 @@ contains
171171
do j = 0, m
172172
patch_id = ib_markers%sf(j, k, l)
173173
if (patch_id /= 0) then
174-
q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp
175-
rho = 0._wp
176-
do i = 1, num_fluids
177-
rho = rho + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l)
178-
end do
174+
call s_decode_patch_periodicity(patch_id, patch_id)
175+
call s_get_neighborhood_idx(patch_id, patch_id)
176+
if (patch_id > 0) then
177+
q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp
178+
rho = 0._wp
179+
do i = 1, num_fluids
180+
rho = rho + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l)
181+
end do
179182

180-
! Sets the momentum
181-
do i = 1, num_dims
182-
q_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho
183-
q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)
184-
end do
183+
! Sets the momentum
184+
do i = 1, num_dims
185+
q_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho
186+
q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)
187+
end do
188+
end if ! patch_id > 0
185189
end if
186190
end do
187191
end do
@@ -925,84 +929,92 @@ contains
925929
encoded_ib_idx = ib_markers%sf(i, j, k)
926930
if (encoded_ib_idx /= 0) then
927931
call s_decode_patch_periodicity(encoded_ib_idx, ib_idx)
928-
929-
! get the vector pointing to the grid cell from the IB centroid
930-
if (num_dims == 3) then
931-
radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, &
932-
& patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid]
933-
else
934-
radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, &
935-
& patch_ib(ib_idx)%y_centroid, 0._wp]
936-
end if
937-
dx = x_cc(i + 1) - x_cc(i)
938-
dy = y_cc(j + 1) - y_cc(j)
939-
940-
local_force_contribution(:) = 0._wp
941-
do fluid_idx = 0, num_fluids - 1
942-
! Get the pressure contribution to force via a finite difference to compute the 2D components of the
943-
! gradient of the pressure and cell volume
944-
local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i &
945-
& + 1, j, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, &
946-
& k))/(2._wp*dx) ! force is the negative pressure gradient
947-
local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, &
948-
& j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
949-
cell_volume = abs(dx*dy)
950-
! add the 3D component of the pressure gradient, if we are working in 3 dimensions
932+
call s_get_neighborhood_idx(ib_idx, ib_idx) ! global patch ID -> local index
933+
if (ib_idx > 0) then
934+
! get the vector pointing to the grid cell from the IB centroid
951935
if (num_dims == 3) then
952-
dz = z_cc(k + 1) - z_cc(k)
953-
local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E + fluid_idx) &
954-
& %sf(i, j, k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, &
955-
& k - 1))/(2._wp*dz)
956-
cell_volume = abs(cell_volume*dz)
936+
radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, &
937+
& patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid]
938+
else
939+
radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, &
940+
& patch_ib(ib_idx)%y_centroid, 0._wp]
957941
end if
958-
end do
959-
960-
! get the viscous stress and add its contribution if that is considered
961-
if (viscous) then
962-
! compute the volume-weighted local dynamic viscosity
963-
dynamic_viscosity = 0._wp
964-
do fluid_idx = 1, num_fluids
965-
! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
966-
dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + eqn_idx%adv%beg - 1)%sf(i, j, &
967-
& k)*dynamic_viscosities(fluid_idx))
942+
dx = x_cc(i + 1) - x_cc(i)
943+
dy = y_cc(j + 1) - y_cc(j)
944+
945+
local_force_contribution(:) = 0._wp
946+
do fluid_idx = 0, num_fluids - 1
947+
! Get the pressure contribution to force via a finite difference to compute the 2D components of the
948+
! gradient of the pressure and cell volume
949+
local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E + fluid_idx) &
950+
& %sf(i + 1, j, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, &
951+
& k))/(2._wp*dx) ! force is the negative pressure gradient
952+
local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E + fluid_idx) &
953+
& %sf(i, j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j - 1, &
954+
& k))/(2._wp*dy)
955+
cell_volume = abs(dx*dy)
956+
! add the 3D component of the pressure gradient, if we are working in 3 dimensions
957+
if (num_dims == 3) then
958+
dz = z_cc(k + 1) - z_cc(k)
959+
local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E + fluid_idx) &
960+
& %sf(i, j, k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, &
961+
& k - 1))/(2._wp*dz)
962+
cell_volume = abs(cell_volume*dz)
963+
end if
968964
end do
969965

970-
! get the linear force components first
971-
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, j, k)
972-
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, j, k)
973-
viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1, &
974-
& 1:3))/(2._wp*dx) ! get x derivative of the first-row of viscous stress tensor
975-
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, &
976-
& 1:3) ! add the x components of the divergence to the force
977-
978-
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j - 1, k)
979-
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j + 1, k)
980-
viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2, &
981-
& 1:3))/(2._wp*dy) ! get y derivative of the second-row of viscous stress tensor
982-
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, &
983-
& 1:3) ! add the y components of the divergence to the force
966+
! get the viscous stress and add its contribution if that is considered
967+
if (viscous) then
968+
! compute the volume-weighted local dynamic viscosity
969+
dynamic_viscosity = 0._wp
970+
do fluid_idx = 1, num_fluids
971+
! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
972+
dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + eqn_idx%adv%beg - 1)%sf(i, j, &
973+
& k)*dynamic_viscosities(fluid_idx))
974+
end do
984975

985-
if (num_dims == 3) then
986-
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j, &
987-
& k - 1)
988-
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j, &
989-
& k + 1)
990-
viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3, &
991-
& 1:3))/(2._wp*dz) ! get z derivative of the third-row of viscous stress tensor
992-
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, &
993-
& 1:3) ! add the z components of the divergence to the force
976+
! get the linear force components first
977+
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, &
978+
& j, k)
979+
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, &
980+
& j, k)
981+
viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1, &
982+
& 1:3))/(2._wp*dx) ! get x derivative of the first-row of viscous stress tensor
983+
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, &
984+
& 1:3) ! add the x components of the divergence to the force
985+
986+
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
987+
& j - 1, k)
988+
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
989+
& j + 1, k)
990+
viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2, &
991+
& 1:3))/(2._wp*dy) ! get y derivative of the second-row of viscous stress tensor
992+
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, &
993+
& 1:3) ! add the y components of the divergence to the force
994+
995+
if (num_dims == 3) then
996+
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
997+
& j, k - 1)
998+
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
999+
& j, k + 1)
1000+
viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3, &
1001+
& 1:3))/(2._wp*dz) &
1002+
& ! get z derivative of the third-row of viscous stress tensor
1003+
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, &
1004+
& 1:3) ! add the z components of the divergence to the force
1005+
end if
9941006
end if
995-
end if
9961007

997-
call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
1008+
call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
9981009

999-
! Update the force and torque values atomically to prevent race conditions
1000-
do l = 1, 3
1001-
$:GPU_ATOMIC(atomic='update')
1002-
forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
1003-
$:GPU_ATOMIC(atomic='update')
1004-
torques(ib_idx, l) = torques(ib_idx, l) + local_torque_contribution(l)*cell_volume
1005-
end do
1010+
! Update the force and torque values atomically to prevent race conditions
1011+
do l = 1, 3
1012+
$:GPU_ATOMIC(atomic='update')
1013+
forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
1014+
$:GPU_ATOMIC(atomic='update')
1015+
torques(ib_idx, l) = torques(ib_idx, l) + local_torque_contribution(l)*cell_volume
1016+
end do
1017+
end if ! ib_idx > 0
10061018
end if
10071019
end do
10081020
end do
@@ -1055,7 +1067,7 @@ contains
10551067
subroutine s_compute_centroid_offset(ib_marker)
10561068

10571069
integer, intent(in) :: ib_marker
1058-
integer :: i, j, k, num_cells, num_cells_local
1070+
integer :: i, j, k, num_cells, num_cells_local, decoded_gbl_id
10591071
real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local
10601072

10611073
! Offset only needs to be computes for specific geometries
@@ -1069,10 +1081,13 @@ contains
10691081
do i = 0, m
10701082
do j = 0, n
10711083
do k = 0, p
1072-
if (ib_markers%sf(i, j, k) == ib_marker) then
1073-
num_cells_local = num_cells_local + 1
1074-
center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp]
1075-
if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k)
1084+
if (ib_markers%sf(i, j, k) /= 0) then
1085+
call s_decode_patch_periodicity(ib_markers%sf(i, j, k), decoded_gbl_id)
1086+
if (decoded_gbl_id == patch_ib(ib_marker)%gbl_patch_id) then
1087+
num_cells_local = num_cells_local + 1
1088+
center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp]
1089+
if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k)
1090+
end if
10761091
end if
10771092
end do
10781093
end do

0 commit comments

Comments
 (0)