diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index a4099937b6..1bfd59af1c 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -889,6 +889,8 @@ contains integer, intent(in) :: time_step + $:GPU_UPDATE(host='[patch_ib(1:num_ibs)]') + if (parallel_io) then call s_write_parallel_ib_data(time_step) else diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index fcfce42fc1..41f8524b87 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -207,23 +207,18 @@ contains integer, intent(in) :: patch_id type(integer_field), intent(inout) :: ib_markers integer, intent(in) :: xp, yp !< integers containing the periodicity projection information - real(wp) :: f, ca_in, pa, ma - real(wp) :: xa, yc, dycdxc + real(wp) :: f, ca_in integer :: i, j, k, il, ir, jl, jr integer :: Np_local, airfoil_id integer :: encoded_patch_id real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame - real(wp), dimension(1:3,1:3) :: inverse_rotation airfoil_id = patch_ib(patch_id)%airfoil_id center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) ca_in = ib_airfoil(airfoil_id)%c - pa = ib_airfoil(airfoil_id)%p - ma = ib_airfoil(airfoil_id)%m Np_local = ib_airfoil_grids(airfoil_id)%Np - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) offset(:) = patch_ib(patch_id)%centroid_offset(:) ! encode the periodicity information into the patch_id @@ -238,24 +233,15 @@ contains call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir) call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr) - $:GPU_PARALLEL_LOOP(private='[i, j, xy_local, k, f, xa, yc, dycdxc]', copyin='[encoded_patch_id, center, & - & inverse_rotation, offset, ma, pa, ca_in, airfoil_id, Np_local, ib_airfoil_grids(airfoil_id)%upper, & - & ib_airfoil_grids(airfoil_id)%lower]', collapse=2) + $:GPU_PARALLEL_LOOP(private='[i, j, xy_local, k, f]', copyin='[encoded_patch_id, center, offset, ca_in, airfoil_id, & + & Np_local, ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]', collapse=2) do j = jl, jr do i = il, ir xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB - xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates + xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xy_local) ! rotate the frame into the IB's coordinates xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then - xa = xy_local(1)/ca_in - if (xa <= pa) then - yc = (ma/pa**2)*(2*pa*xa - xa**2) - dycdxc = (2*ma/pa**2)*(pa - xa) - else - yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) - dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) - end if if (xy_local(2) >= 0._wp) then k = 1 do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < xy_local(1) .and. k <= Np_local) @@ -304,12 +290,11 @@ contains integer, intent(in) :: patch_id type(integer_field), intent(inout) :: ib_markers integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information - real(wp) :: lz, z_max, z_min, f, ca_in, pa, ma + real(wp) :: lz, z_max, z_min, f, ca_in integer :: i, j, k, l, il, ir, jl, jr, ll, lr - integer :: Np_local, airfoil_id + integer :: airfoil_id integer :: encoded_patch_id real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame - real(wp), dimension(1:3,1:3) :: inverse_rotation airfoil_id = patch_ib(patch_id)%airfoil_id center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) @@ -317,10 +302,6 @@ contains center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg) lz = patch_ib(patch_id)%length_z ca_in = ib_airfoil(airfoil_id)%c - pa = ib_airfoil(airfoil_id)%p - ma = ib_airfoil(airfoil_id)%m - Np_local = ib_airfoil_grids(airfoil_id)%Np - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) offset(:) = patch_ib(patch_id)%centroid_offset(:) z_max = lz/2 @@ -341,15 +322,15 @@ contains call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr) call get_bounding_indices(center(3) - ca_in, center(3) + ca_in, z_cc, ll, lr) - $:GPU_PARALLEL_LOOP(private='[i, j, l, xyz_local, k, f]', copyin='[encoded_patch_id, center, inverse_rotation, offset, & - & ma, pa, ca_in, airfoil_id, Np_local, ib_airfoil_grids(airfoil_id)%upper, & - & ib_airfoil_grids(airfoil_id)%lower, z_min, z_max]', collapse=3) + $:GPU_PARALLEL_LOOP(private='[i, j, l, xyz_local, k, f]', copyin='[encoded_patch_id, center, offset, ca_in, airfoil_id, & + & ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower, z_min, z_max]', collapse=3) do l = ll, lr do j = jl, jr do i = il, ir ! get coordinate frame centered on IB xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates + ! rotate the frame into the IB's coordinates + xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then @@ -410,15 +391,11 @@ contains real(wp) :: corner_distance !< Equation of state parameters real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame - real(wp), dimension(1:3,1:3) :: inverse_rotation ! Transferring the rectangle's centroid and length information center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) - length(1) = patch_ib(patch_id)%length_x - length(2) = patch_ib(patch_id)%length_y - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) ! encode the periodicity information into the patch_id call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) @@ -428,21 +405,22 @@ contains jl = -gp_layers - 1 ir = m + gp_layers + 1 jr = n + gp_layers + 1 - corner_distance = sqrt(dot_product(length, length))/2._wp ! maximum distance any marker can be from the center + ! maximum distance any marker can be from the center + corner_distance = 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2) call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir) call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr) ! Assign primitive variables if rectangle covers cell and patch has write permission - $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center, length, inverse_rotation, x_cc, & - & y_cc]', collapse=2) + $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center]', collapse=2) do j = jl, jr do i = il, ir ! get the x and y coordinates in the local IB frame xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] - xy_local = matmul(inverse_rotation, xy_local) + xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xy_local) - if (-0.5_wp*length(1) <= xy_local(1) .and. 0.5_wp*length(1) >= xy_local(1) .and. -0.5_wp*length(2) <= xy_local(2) & - & .and. 0.5_wp*length(2) >= xy_local(2)) then + if (-0.5_wp*patch_ib(patch_id)%length_x <= xy_local(1) .and. 0.5_wp*patch_ib(patch_id)%length_x >= xy_local(1) & + & .and. -0.5_wp*patch_ib(patch_id)%length_y <= xy_local(2) & + & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xy_local(2)) then ! Updating the patch identities bookkeeping variable ib_markers%sf(i, j, 0) = encoded_patch_id end if @@ -516,19 +494,14 @@ contains integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information integer :: i, j, k, ir, il, jr, jl, kr, kl !< Generic loop iterators integer :: encoded_patch_id - real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame - real(wp), dimension(1:3,1:3) :: inverse_rotation + real(wp), dimension(1:3) :: xyz_local, center !< x and y coordinates in local IB frame real(wp) :: corner_distance - ! Transferring the cuboid's centroid and length information + ! Transferring the cuboid's centroid center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg) - length(1) = patch_ib(patch_id)%length_x - length(2) = patch_ib(patch_id)%length_y - length(3) = patch_ib(patch_id)%length_z - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) ! encode the periodicity information into the patch_id call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) @@ -540,7 +513,8 @@ contains ir = m + gp_layers + 1 jr = n + gp_layers + 1 kr = p + gp_layers + 1 - corner_distance = sqrt(dot_product(length, length))/2._wp ! maximum distance any marker can be from the center + corner_distance = 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2 & + & + patch_ib(patch_id)%length_z**2) ! maximum distance any marker can be from the center call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir) call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr) call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr) @@ -548,17 +522,20 @@ contains ! Checking whether the cuboid covers a particular cell in the domain and verifying whether the current patch has permission ! to write to to that cell. If both queries check out, the primitive variables of the current patch are assigned to this ! cell. - $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center, length, inverse_rotation]', & - & collapse=3) + $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center]', collapse=3) do k = kl, kr do j = jl, jr do i = il, ir xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates - - if (-0.5*length(1) <= xyz_local(1) .and. 0.5*length(1) >= xyz_local(1) .and. -0.5*length(2) <= xyz_local(2) & - & .and. 0.5*length(2) >= xyz_local(2) .and. -0.5*length(3) <= xyz_local(3) .and. 0.5*length(3) & - & >= xyz_local(3)) then + ! rotate the frame into the IB's coordinates + xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) + + if (-0.5_wp*patch_ib(patch_id)%length_x <= xyz_local(1) & + & .and. 0.5_wp*patch_ib(patch_id)%length_x >= xyz_local(1) .and. & + & -0.5_wp*patch_ib(patch_id)%length_y <= xyz_local(2) & + & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xyz_local(2) .and. & + & -0.5_wp*patch_ib(patch_id)%length_z <= xyz_local(3) & + & .and. 0.5_wp*patch_ib(patch_id)%length_z >= xyz_local(3)) then ! Updating the patch identities bookkeeping variable ib_markers%sf(i, j, k) = encoded_patch_id end if @@ -577,21 +554,16 @@ contains integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators integer :: encoded_patch_id - real(wp) :: radius real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame real(wp), dimension(1:3,1:3) :: inverse_rotation real(wp) :: corner_distance - ! Transferring the cylindrical patch's centroid, length, radius, + ! Transferring the cylindrical patch's centroid center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg) - length(1) = patch_ib(patch_id)%length_x - length(2) = patch_ib(patch_id)%length_y - length(3) = patch_ib(patch_id)%length_z - radius = patch_ib(patch_id)%radius - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) + length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, patch_ib(patch_id)%length_z] ! encode the periodicity information into the patch_id call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id) @@ -602,7 +574,7 @@ contains ir = m + gp_layers + 1 jr = n + gp_layers + 1 kr = p + gp_layers + 1 - corner_distance = sqrt(radius**2 + maxval(length)**2) ! distance to rim of cylinder + corner_distance = sqrt(patch_ib(patch_id)%radius**2 + maxval(length)**2) ! distance to rim of cylinder call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir) call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr) call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr) @@ -610,20 +582,23 @@ contains ! Checking whether the cylinder covers a particular cell in the domain and verifying whether the current patch has the ! permission to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to ! this cell. - $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center, length, radius, & - & inverse_rotation]', collapse=3) + $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center]', collapse=3) do k = kl, kr do j = jl, jr do i = il, ir xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates - - if (((.not. f_is_default(length(1)) .and. xyz_local(2)**2 + xyz_local(3)**2 <= radius**2 .and. & - & -0.5_wp*length(1) <= xyz_local(1) .and. 0.5_wp*length(1) >= xyz_local(1)) & - & .or. (.not. f_is_default(length(2)) .and. xyz_local(1)**2 + xyz_local(3)**2 <= radius**2 .and. & - & -0.5_wp*length(2) <= xyz_local(2) .and. 0.5_wp*length(2) >= xyz_local(2)) & - & .or. (.not. f_is_default(length(3)) .and. xyz_local(1)**2 + xyz_local(2)**2 <= radius**2 .and. & - & -0.5_wp*length(3) <= xyz_local(3) .and. 0.5_wp*length(3) >= xyz_local(3)))) then + ! rotate the frame into the IB's coordinates + xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) + + if (((.not. f_is_default(patch_ib(patch_id)%length_x) .and. xyz_local(2)**2 + xyz_local(3) & + & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_x <= xyz_local(1) & + & .and. 0.5_wp*patch_ib(patch_id)%length_x >= xyz_local(1)) & + & .or. (.not. f_is_default(patch_ib(patch_id)%length_y) .and. xyz_local(1)**2 + xyz_local(3) & + & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_y <= xyz_local(2) & + & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xyz_local(2)) & + & .or. (.not. f_is_default(patch_ib(patch_id)%length_z) .and. xyz_local(1)**2 + xyz_local(2) & + & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_z <= xyz_local(3) & + & .and. 0.5_wp*patch_ib(patch_id)%length_z >= xyz_local(3)))) then ! Updating the patch identities bookkeeping variable ib_markers%sf(i, j, k) = encoded_patch_id end if @@ -643,40 +618,37 @@ contains integer :: i, j, il, ir, jl, jr !< Generic loop iterators integer :: encoded_patch_id real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame - real(wp), dimension(1:2) :: ellipse_coeffs !< a and b in the ellipse coefficients real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame - real(wp), dimension(1:3,1:3) :: inverse_rotation + real(wp) :: bounding_radius - ! Transferring the ellipse's centroid and length information + ! Transferring the ellipse's centroid center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) - ellipse_coeffs(1) = 0.5_wp*patch_ib(patch_id)%length_x - ellipse_coeffs(2) = 0.5_wp*patch_ib(patch_id)%length_y - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) ! encode the periodicity information into the patch_id call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id) ! find the indices to the left and right of the IB in i, j, k + bounding_radius = 0.5_wp*max(patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y) il = -gp_layers - 1 jl = -gp_layers - 1 ir = m + gp_layers + 1 jr = n + gp_layers + 1 - call get_bounding_indices(center(1) - maxval(ellipse_coeffs)*2._wp, center(1) + maxval(ellipse_coeffs)*2._wp, x_cc, il, ir) - call get_bounding_indices(center(2) - maxval(ellipse_coeffs)*2._wp, center(2) + maxval(ellipse_coeffs)*2._wp, y_cc, jl, jr) + call get_bounding_indices(center(1) - bounding_radius*2._wp, center(1) + bounding_radius*2._wp, x_cc, il, ir) + call get_bounding_indices(center(2) - bounding_radius*2._wp, center(2) + bounding_radius*2._wp, y_cc, jl, jr) ! Checking whether the ellipse covers a particular cell in the domain - $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center, ellipse_coeffs, inverse_rotation, & - & x_cc, y_cc]', collapse=2) + $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center]', collapse=2) do j = jl, jr do i = il, ir ! get the x and y coordinates in the local IB frame xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] - xy_local = matmul(inverse_rotation, xy_local) + xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse(:,:), xy_local) ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1 - if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then + if ((xy_local(1)/(0.5_wp*patch_ib(patch_id)%length_x))**2 + (xy_local(2)/(0.5_wp*patch_ib(patch_id)%length_y)) & + & **2 <= 1._wp) then ! Updating the patch identities bookkeeping variable ib_markers%sf(i, j, 0) = encoded_patch_id end if @@ -853,6 +825,8 @@ contains !> Compute a rotation matrix for converting to the rotating frame of the boundary subroutine s_update_ib_rotation_matrix(patch_id) + $:GPU_ROUTINE(parallelism='[seq]') + integer, intent(in) :: patch_id real(wp), dimension(3, 3, 3) :: rotation real(wp) :: angle diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index b710583965..253123184d 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -76,14 +76,23 @@ contains call nvtxStartRange("SETUP-IBM-MODULE") + ! GPU routines require updated cell centers + $:GPU_UPDATE(device='[num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, x_domain, y_domain, ib_bc_x%beg, ib_bc_y%beg]') + if (p /= 0) then + $:GPU_UPDATE(device='[z_cc, dz, z_domain, ib_bc_z%beg]') + end if + $:GPU_UPDATE(device='[patch_ib(1:num_ibs)]') + ! do all set up for moving immersed boundaries + $:GPU_PARALLEL_LOOP(private='[i]') do i = 1, num_ibs if (patch_ib(i)%moving_ibm /= 0) then call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) end if call s_update_ib_rotation_matrix(i) end do - $:GPU_UPDATE(device='[patch_ib(1:num_ibs)]') + $:END_GPU_PARALLEL_LOOP() + $:GPU_UPDATE(host='[patch_ib(1:num_ibs)]') ! allocate some arrays for MPI communication, if required by this simulation #ifdef MFC_MPI @@ -94,11 +103,6 @@ contains end if #endif - ! GPU routines require updated cell centers - $:GPU_UPDATE(device='[num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, x_domain, y_domain, ib_bc_x%beg, ib_bc_y%beg]') - if (p /= 0) then - $:GPU_UPDATE(device='[z_cc, dz, z_domain, ib_bc_z%beg]') - end if call s_update_ib_lookup() ! recompute the new ib_patch locations @@ -868,13 +872,13 @@ contains $:END_GPU_PARALLEL_LOOP() ! recalulcate the rotation matrix based upon the new angles + $:GPU_PARALLEL_LOOP(private='[i]') do i = 1, num_ibs if (patch_ib(i)%moving_ibm /= 0) then call s_update_ib_rotation_matrix(i) end if end do - - $:GPU_UPDATE(device='[patch_ib]') + $:END_GPU_PARALLEL_LOOP() ! recompute the new ib_patch locations and broadcast them. call nvtxStartRange("APPLY-IB-PATCHES") @@ -933,7 +937,7 @@ contains $:GPU_PARALLEL_LOOP(private='[ib_idx, ib_idx_temp, encoded_ib_idx, fluid_idx, radial_vector, local_force_contribution, & & cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress_div, & & viscous_stress_div_1, viscous_stress_div_2, dx, dy, dz]', copy='[forces, torques]', & - & copyin='[patch_ib, dynamic_viscosities]', collapse=3) + & copyin='[dynamic_viscosities]', collapse=3) do i = 0, m do j = 0, n do k = 0, p @@ -1050,10 +1054,12 @@ contains end do ! apply the summed forces + $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') do i = 1, num_ibs patch_ib(i)%force(:) = forces(i,:) patch_ib(i)%torque(:) = torques(i,:) end do + $:END_GPU_PARALLEL_LOOP() call nvtxEndRange @@ -1122,6 +1128,8 @@ contains !> Computes the moment of inertia for an immersed boundary subroutine s_compute_moment_of_inertia(ib_idx, axis) + $:GPU_ROUTINE(parallelism='[seq]') + real(wp), dimension(3), intent(in) :: axis !< the axis about which we compute the moment. Only required in 3D. integer, intent(in) :: ib_idx real(wp) :: moment, distance_to_axis, cell_volume @@ -1158,13 +1166,10 @@ contains ib_marker = patch_ib(ib_idx)%gbl_patch_id - $:GPU_PARALLEL_LOOP(private='[position, closest_point_along_axis, vector_to_axis, distance_to_axis]', copy='[moment, & - & count]', copyin='[ib_marker, cell_volume, normal_axis]', collapse=3) do i = 0, m do j = 0, n do k = 0, p if (ib_markers%sf(i, j, k) == ib_marker) then - $:GPU_ATOMIC(atomic='update') count = count + 1 ! increment the count of total cells in the boundary ! get the position in local coordinates so that the axis passes through 0, 0, 0 @@ -1182,13 +1187,11 @@ contains distance_to_axis = dot_product(vector_to_axis, vector_to_axis) ! saves the distance to the axis squared ! compute the position component of the moment - $:GPU_ATOMIC(atomic='update') moment = moment + distance_to_axis end if end do end do end do - $:END_GPU_PARALLEL_LOOP() ! write the final moment assuming the points are all uniform density patch_ib(ib_idx)%moment = moment*patch_ib(ib_idx)%mass/(count*cell_volume) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 0f92271a0f..252301892a 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -714,6 +714,7 @@ contains if (moving_immersed_boundary_flag) call s_compute_ib_forces(q_prim_vf, fluid_pp) + $:GPU_PARALLEL_LOOP(private='[i, gbl_id]', copyin='[s]') do i = 1, num_ibs if (s == 1) then patch_ib(i)%step_vel = patch_ib(i)%vel @@ -759,8 +760,8 @@ contains & 2)*patch_ib(i)%z_centroid + rk_coef(s, 3)*patch_ib(i)%vel(3)*dt)/rk_coef(s, 4) end if end do + $:END_GPU_PARALLEL_LOOP() - $:GPU_UPDATE(device='[patch_ib]') call s_update_mib(num_ibs) call nvtxEndRange