diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index baf9d0591c..dad7035583 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -34,32 +34,35 @@ contains integer, intent(IN) :: ib_patch_id real(wp) :: radius, dist - real(wp) :: x_centroid, y_centroid + real(wp), dimension(2) :: center real(wp), dimension(3) :: dist_vec integer :: i, j !< Loop index variables radius = patch_ib(ib_patch_id)%radius - x_centroid = patch_ib(ib_patch_id)%x_centroid - y_centroid = patch_ib(ib_patch_id)%y_centroid - - do i = 0, m - do j = 0, n - - dist_vec(1) = x_cc(i) - x_centroid - dist_vec(2) = y_cc(j) - y_centroid - dist_vec(3) = 0 - dist = sqrt(sum(dist_vec**2)) - levelset%sf(i, j, 0, ib_patch_id) = dist - radius - if (f_approx_equal(dist, 0._wp)) then - levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0 - else - levelset_norm%sf(i, j, 0, ib_patch_id, :) = & - dist_vec(:)/dist - end if + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid + #:call GPU_PARALLEL_LOOP(private='[i,j,dist_vec,dist]', & + & copyin='[ib_patch_id,center,radius]', collapse=2) + do i = 0, m + do j = 0, n + + dist_vec(1) = x_cc(i) - center(1) + dist_vec(2) = y_cc(j) - center(2) + dist_vec(3) = 0._wp + dist = sqrt(sum(dist_vec**2)) + levelset%sf(i, j, 0, ib_patch_id) = dist - radius + if (f_approx_equal(dist, 0._wp)) then + levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0 + else + levelset_norm%sf(i, j, 0, ib_patch_id, :) = & + dist_vec(:)/dist + end if + + end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_circle_levelset @@ -71,78 +74,81 @@ contains real(wp) :: dist, global_dist integer :: global_id - real(wp) :: x_centroid, y_centroid real(wp), dimension(3) :: dist_vec real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame + real(wp), dimension(1:2) :: center real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation integer :: i, j, k !< Loop index variables - x_centroid = patch_ib(ib_patch_id)%x_centroid - y_centroid = patch_ib(ib_patch_id)%y_centroid + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) - do i = 0, m - do j = 0, n - xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp] ! get coordinate frame centered on IB - xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate - - if (xy_local(2) >= 0._wp) then - ! finds the location on the airfoil grid with the minimum distance (closest) - do k = 1, Np - dist_vec(1) = xy_local(1) - airfoil_grid_u(k)%x - dist_vec(2) = xy_local(2) - airfoil_grid_u(k)%y - dist_vec(3) = 0._wp - dist = sqrt(sum(dist_vec**2)) - if (k == 1) then - global_dist = dist - global_id = k - else - if (dist < global_dist) then + #:call GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,dist_vec,dist,global_dist,global_id]', & + & copyin='[ib_patch_id,center,rotation,inverse_rotation,airfoil_grid_u,airfoil_grid_l]', collapse=2) + do i = 0, m + do j = 0, n + 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 coordinate + + if (xy_local(2) >= 0._wp) then + ! finds the location on the airfoil grid with the minimum distance (closest) + do k = 1, Np + dist_vec(1) = xy_local(1) - airfoil_grid_u(k)%x + dist_vec(2) = xy_local(2) - airfoil_grid_u(k)%y + dist_vec(3) = 0._wp + dist = sqrt(sum(dist_vec**2)) + if (k == 1) then global_dist = dist global_id = k + else + if (dist < global_dist) then + global_dist = dist + global_id = k + end if end if - end if - end do - dist_vec(1) = xy_local(1) - airfoil_grid_u(global_id)%x - dist_vec(2) = xy_local(2) - airfoil_grid_u(global_id)%y - dist_vec(3) = 0 - dist = global_dist - else - ! TODO :: This looks identical to the above code but using the lower array. Refactor this. - do k = 1, Np - dist_vec(1) = xy_local(1) - airfoil_grid_l(k)%x - dist_vec(2) = xy_local(2) - airfoil_grid_l(k)%y + end do + dist_vec(1) = xy_local(1) - airfoil_grid_u(global_id)%x + dist_vec(2) = xy_local(2) - airfoil_grid_u(global_id)%y dist_vec(3) = 0 - dist = sqrt(sum(dist_vec**2)) - if (k == 1) then - global_dist = dist - global_id = k - else - if (dist < global_dist) then + dist = global_dist + else + ! TODO :: This looks identical to the above code but using the lower array. Refactor this. + do k = 1, Np + dist_vec(1) = xy_local(1) - airfoil_grid_l(k)%x + dist_vec(2) = xy_local(2) - airfoil_grid_l(k)%y + dist_vec(3) = 0 + dist = sqrt(sum(dist_vec**2)) + if (k == 1) then global_dist = dist global_id = k + else + if (dist < global_dist) then + global_dist = dist + global_id = k + end if end if - end if - end do - dist_vec(1) = xy_local(1) - airfoil_grid_l(global_id)%x - dist_vec(2) = xy_local(2) - airfoil_grid_l(global_id)%y - dist_vec(3) = 0 - dist = global_dist - end if - - levelset%sf(i, j, 0, ib_patch_id) = dist - if (f_approx_equal(dist, 0._wp)) then - levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp - else - levelset_norm%sf(i, j, 0, ib_patch_id, :) = & - matmul(rotation, dist_vec(:))/dist ! convert the normal vector back to global grid coordinates - end if + end do + dist_vec(1) = xy_local(1) - airfoil_grid_l(global_id)%x + dist_vec(2) = xy_local(2) - airfoil_grid_l(global_id)%y + dist_vec(3) = 0 + dist = global_dist + end if + + levelset%sf(i, j, 0, ib_patch_id) = dist + if (f_approx_equal(dist, 0._wp)) then + levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp + else + levelset_norm%sf(i, j, 0, ib_patch_id, :) = & + matmul(rotation, dist_vec(:))/dist ! convert the normal vector back to global grid coordinates + end if + end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_airfoil_levelset @@ -154,99 +160,102 @@ contains real(wp) :: dist, dist_surf, dist_side, global_dist integer :: global_id - real(wp) :: x_centroid, y_centroid, z_centroid, lz, z_max, z_min + real(wp) :: lz, z_max, z_min real(wp), dimension(3) :: dist_vec - real(wp), dimension(1:3) :: xyz_local !< x, y, z coordinates in local IB frame + real(wp), dimension(1:3) :: xyz_local, center !< x, y, z coordinates in local IB frame real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation real(wp) :: length_z integer :: i, j, k, l !< Loop index variables - x_centroid = patch_ib(ib_patch_id)%x_centroid - y_centroid = patch_ib(ib_patch_id)%y_centroid - z_centroid = patch_ib(ib_patch_id)%z_centroid + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid + center(3) = patch_ib(ib_patch_id)%z_centroid lz = patch_ib(ib_patch_id)%length_z inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) - z_max = z_centroid + lz/2 - z_min = z_centroid - lz/2 - - do l = 0, p - do j = 0, n - do i = 0, m - - xyz_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_cc(l) - z_centroid] ! get coordinate frame centered on IB - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates - - if (xyz_local(2) >= y_centroid) then - do k = 1, Np - dist_vec(1) = xyz_local(1) - airfoil_grid_u(k)%x - dist_vec(2) = xyz_local(2) - airfoil_grid_u(k)%y - dist_vec(3) = 0 - dist_surf = sqrt(sum(dist_vec**2)) - if (k == 1) then - global_dist = dist_surf - global_id = k - else - if (dist_surf < global_dist) then + z_max = center(3) + lz/2 + z_min = center(3) - lz/2 + + #:call GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,dist_vec,dist,global_dist,global_id,dist_side,dist_surf]', & + & copyin='[ib_patch_id,center,rotation,inverse_rotation,airfoil_grid_u,airfoil_grid_l,z_min,z_max]', collapse=3) + do l = 0, p + do j = 0, n + do i = 0, m + + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB + xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates + + if (xyz_local(2) >= center(2)) then + do k = 1, Np + dist_vec(1) = xyz_local(1) - airfoil_grid_u(k)%x + dist_vec(2) = xyz_local(2) - airfoil_grid_u(k)%y + dist_vec(3) = 0 + dist_surf = sqrt(sum(dist_vec**2)) + if (k == 1) then global_dist = dist_surf global_id = k + else + if (dist_surf < global_dist) then + global_dist = dist_surf + global_id = k + end if end if - end if - end do - dist_vec(1) = xyz_local(1) - airfoil_grid_u(global_id)%x - dist_vec(2) = xyz_local(2) - airfoil_grid_u(global_id)%y - dist_vec(3) = 0 - dist_surf = global_dist - else - do k = 1, Np - dist_vec(1) = xyz_local(1) - airfoil_grid_l(k)%x - dist_vec(2) = xyz_local(2) - airfoil_grid_l(k)%y + end do + dist_vec(1) = xyz_local(1) - airfoil_grid_u(global_id)%x + dist_vec(2) = xyz_local(2) - airfoil_grid_u(global_id)%y dist_vec(3) = 0 - dist_surf = sqrt(sum(dist_vec**2)) - if (k == 1) then - global_dist = dist_surf - global_id = k - else - if (dist_surf < global_dist) then + dist_surf = global_dist + else + do k = 1, Np + dist_vec(1) = xyz_local(1) - airfoil_grid_l(k)%x + dist_vec(2) = xyz_local(2) - airfoil_grid_l(k)%y + dist_vec(3) = 0 + dist_surf = sqrt(sum(dist_vec**2)) + if (k == 1) then global_dist = dist_surf global_id = k + else + if (dist_surf < global_dist) then + global_dist = dist_surf + global_id = k + end if end if - end if - end do - dist_vec(1) = xyz_local(1) - airfoil_grid_l(global_id)%x - dist_vec(2) = xyz_local(2) - airfoil_grid_l(global_id)%y - dist_vec(3) = 0 - dist_surf = global_dist - end if + end do + dist_vec(1) = xyz_local(1) - airfoil_grid_l(global_id)%x + dist_vec(2) = xyz_local(2) - airfoil_grid_l(global_id)%y + dist_vec(3) = 0 + dist_surf = global_dist + end if - dist_side = min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l))) + dist_side = min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l))) - if (dist_side < dist_surf) then - levelset%sf(i, j, l, ib_patch_id) = dist_side - if (f_approx_equal(dist_side, abs(z_cc(l) - z_min))) then - levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, -1/) - else - levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, 1/) - end if - levelset_norm%sf(i, j, l, ib_patch_id, :) = & - matmul(rotation, levelset_norm%sf(i, j, l, ib_patch_id, :)/dist_surf) - else - levelset%sf(i, j, l, ib_patch_id) = dist_surf - if (f_approx_equal(dist_surf, 0._wp)) then - levelset_norm%sf(i, j, l, ib_patch_id, :) = 0 - else + if (dist_side < dist_surf) then + levelset%sf(i, j, l, ib_patch_id) = dist_side + if (f_approx_equal(dist_side, abs(z_cc(l) - z_min))) then + levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, -1/) + else + levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, 1/) + end if levelset_norm%sf(i, j, l, ib_patch_id, :) = & - matmul(rotation, dist_vec(:)/dist_surf) + matmul(rotation, levelset_norm%sf(i, j, l, ib_patch_id, :)/dist_surf) + else + levelset%sf(i, j, l, ib_patch_id) = dist_surf + if (f_approx_equal(dist_surf, 0._wp)) then + levelset_norm%sf(i, j, l, ib_patch_id, :) = 0 + else + levelset_norm%sf(i, j, l, ib_patch_id, :) = & + matmul(rotation, dist_vec(:)/dist_surf) + end if end if - end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_3D_airfoil_levelset @@ -261,9 +270,9 @@ contains real(wp) :: min_dist real(wp) :: side_dists(4) - real(wp) :: x_centroid, y_centroid real(wp) :: length_x, length_y - real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame + real(wp), dimension(1:3) :: xy_local, dist_vec !< x and y coordinates in local IB frame + real(wp), dimension(2) :: center !< x and y coordinates in local IB frame real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation integer :: i, j, k !< Loop index variables @@ -271,8 +280,8 @@ contains length_x = patch_ib(ib_patch_id)%length_x length_y = patch_ib(ib_patch_id)%length_y - x_centroid = patch_ib(ib_patch_id)%x_centroid - y_centroid = patch_ib(ib_patch_id)%y_centroid + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) @@ -281,49 +290,49 @@ contains bottom_left(1) = -length_x/2 bottom_left(2) = -length_y/2 - do i = 0, m - do j = 0, n - xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp] - xy_local = matmul(inverse_rotation, xy_local) - - if ((xy_local(1) > bottom_left(1) .and. xy_local(1) < top_right(1)) .or. & - (xy_local(2) > bottom_left(2) .and. xy_local(2) < top_right(2))) then - - side_dists(1) = bottom_left(1) - xy_local(1) - side_dists(2) = top_right(1) - xy_local(1) - side_dists(3) = bottom_left(2) - xy_local(2) - side_dists(4) = top_right(2) - xy_local(2) - min_dist = initial_distance_buffer - idx = 1 - - do k = 1, 4 - if (abs(side_dists(k)) < abs(min_dist)) then - idx = k - min_dist = side_dists(idx) - end if - end do + #:call GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', & + & copyin='[ib_patch_id,center,bottom_left,top_right,inverse_rotation,rotation]', collapse=2) + do i = 0, m + do j = 0, n + xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] + xy_local = matmul(inverse_rotation, xy_local) + + if ((xy_local(1) > bottom_left(1) .and. xy_local(1) < top_right(1)) .or. & + (xy_local(2) > bottom_left(2) .and. xy_local(2) < top_right(2))) then + + side_dists(1) = bottom_left(1) - xy_local(1) + side_dists(2) = top_right(1) - xy_local(1) + side_dists(3) = bottom_left(2) - xy_local(2) + side_dists(4) = top_right(2) - xy_local(2) + min_dist = side_dists(1) + idx = 1 + + do k = 2, 4 + if (abs(side_dists(k)) < abs(min_dist)) then + idx = k + min_dist = side_dists(idx) + end if + end do - levelset%sf(i, j, 0, ib_patch_id) = side_dists(idx) - levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp - if (.not. f_approx_equal(side_dists(idx), 0._wp)) then - if (idx == 1 .or. idx == 2) then - ! vector points along the x axis - levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(idx)/ & - abs(side_dists(idx)) + levelset%sf(i, j, 0, ib_patch_id) = side_dists(idx) + dist_vec = 0._wp + if (.not. f_approx_equal(side_dists(idx), 0._wp)) then + if (idx == 1 .or. idx == 2) then + ! vector points along the x axis + dist_vec(1) = side_dists(idx)/abs(side_dists(idx)) + else + ! vector points along the y axis + dist_vec(2) = side_dists(idx)/abs(side_dists(idx)) + end if + ! convert the normal vector back into the global coordinate system + levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, dist_vec) else - ! vector points along the y axis - levelset_norm%sf(i, j, 0, ib_patch_id, 2) = side_dists(idx)/ & - abs(side_dists(idx)) + levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp end if - ! convert the normal vector back into the global coordinate system - levelset_norm%sf(i, j, 0, ib_patch_id, :) = & - matmul(rotation, levelset_norm%sf(i, j, 0, ib_patch_id, :)) - else - levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp end if - end if + end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_rectangle_levelset @@ -334,111 +343,105 @@ contains integer, intent(IN) :: ib_patch_id real(wp) :: Right, Left, Bottom, Top, Front, Back - real(wp) :: x, y, z, min_dist + real(wp) :: min_dist real(wp) :: side_dists(6) - real(wp) :: x_centroid, y_centroid, z_centroid + real(wp), dimension(3) :: center real(wp) :: length_x, length_y, length_z - real(wp), dimension(1:3) :: xyz_local !< x and y coordinates in local IB frame + real(wp), dimension(1:3) :: xyz_local, dist_vec !< x and y coordinates in local IB frame real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation - integer :: i, j, k, l, idx !< Loop index variables + integer :: i, j, k !< Loop index variables length_x = patch_ib(ib_patch_id)%length_x length_y = patch_ib(ib_patch_id)%length_y length_z = patch_ib(ib_patch_id)%length_z - x_centroid = patch_ib(ib_patch_id)%x_centroid - y_centroid = patch_ib(ib_patch_id)%y_centroid - z_centroid = patch_ib(ib_patch_id)%z_centroid + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid + center(3) = patch_ib(ib_patch_id)%z_centroid inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) Right = length_x/2 Left = -length_x/2 - Top = length_y/2 Bottom = -length_y/2 - Front = length_z/2 Back = -length_z/2 - do i = 0, m - do j = 0, n - do k = 0, p - - xyz_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_cc(k) - z_centroid] ! get coordinate frame centered on IB - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinate - - if ((xyz_local(1) > Left .and. xyz_local(1) < Right) .or. & - (xyz_local(2) > Bottom .and. xyz_local(2) < Top) .or. & - (xyz_local(3) > Back .and. xyz_local(3) < Front)) then - - side_dists(1) = Left - xyz_local(1) - side_dists(2) = xyz_local(1) - Right - side_dists(3) = Bottom - xyz_local(2) - side_dists(4) = xyz_local(2) - Top - side_dists(5) = Back - xyz_local(3) - side_dists(6) = xyz_local(3) - Front - min_dist = minval(abs(side_dists)) - - ! TODO :: The way that this is written, it looks like we will - ! trigger at the first size that is close to the minimum distance, - ! meaning corners where side_dists are the same will - ! trigger on what may not actually be the minimum, - ! leading to undesired behavior. This should be resolved - ! and this code should be cleaned up. I verified this behavior - ! with tests. - levelset_norm%sf(i, j, k, ib_patch_id, :) = 0._wp - if (f_approx_equal(min_dist, abs(side_dists(1)))) then - levelset%sf(i, j, k, ib_patch_id) = side_dists(1) - if (.not. f_approx_equal(side_dists(1), 0._wp)) then - levelset_norm%sf(i, j, k, ib_patch_id, 1) = side_dists(1)/ & - abs(side_dists(1)) - end if + #:call GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local,dist_vec]', & + & copyin='[ib_patch_id,center,inverse_rotation,rotation,Right,Left,Top,Bottom,Front,Back]', collapse=3) + do i = 0, m + do j = 0, n + do k = 0, p + + 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 coordinate + + if ((xyz_local(1) > Left .and. xyz_local(1) < Right) .or. & + (xyz_local(2) > Bottom .and. xyz_local(2) < Top) .or. & + (xyz_local(3) > Back .and. xyz_local(3) < Front)) then + + side_dists(1) = Left - xyz_local(1) + side_dists(2) = xyz_local(1) - Right + side_dists(3) = Bottom - xyz_local(2) + side_dists(4) = xyz_local(2) - Top + side_dists(5) = Back - xyz_local(3) + side_dists(6) = xyz_local(3) - Front + min_dist = minval(abs(side_dists)) + + ! TODO :: The way that this is written, it looks like we will + ! trigger at the first size that is close to the minimum distance, + ! meaning corners where side_dists are the same will + ! trigger on what may not actually be the minimum, + ! leading to undesired behavior. This should be resolved + ! and this code should be cleaned up. It also means that + ! rotating the box 90 degrees will cause tests to fail. + dist_vec = 0._wp + if (f_approx_equal(min_dist, abs(side_dists(1)))) then + levelset%sf(i, j, k, ib_patch_id) = side_dists(1) + if (.not. f_approx_equal(side_dists(1), 0._wp)) then + dist_vec(1) = side_dists(1)/abs(side_dists(1)) + end if - else if (f_approx_equal(min_dist, abs(side_dists(2)))) then - levelset%sf(i, j, k, ib_patch_id) = side_dists(2) - if (.not. f_approx_equal(side_dists(2), 0._wp)) then - levelset_norm%sf(i, j, k, ib_patch_id, 1) = -side_dists(2)/ & - abs(side_dists(2)) - end if + else if (f_approx_equal(min_dist, abs(side_dists(2)))) then + levelset%sf(i, j, k, ib_patch_id) = side_dists(2) + if (.not. f_approx_equal(side_dists(2), 0._wp)) then + dist_vec(1) = -side_dists(2)/abs(side_dists(2)) + end if - else if (f_approx_equal(min_dist, abs(side_dists(3)))) then - levelset%sf(i, j, k, ib_patch_id) = side_dists(3) - if (.not. f_approx_equal(side_dists(3), 0._wp)) then - levelset_norm%sf(i, j, k, ib_patch_id, 2) = side_dists(3)/ & - abs(side_dists(3)) - end if + else if (f_approx_equal(min_dist, abs(side_dists(3)))) then + levelset%sf(i, j, k, ib_patch_id) = side_dists(3) + if (.not. f_approx_equal(side_dists(3), 0._wp)) then + dist_vec(2) = side_dists(3)/abs(side_dists(3)) + end if - else if (f_approx_equal(min_dist, abs(side_dists(4)))) then - levelset%sf(i, j, k, ib_patch_id) = side_dists(4) - if (.not. f_approx_equal(side_dists(4), 0._wp)) then - levelset_norm%sf(i, j, k, ib_patch_id, 2) = -side_dists(4)/ & - abs(side_dists(4)) - end if + else if (f_approx_equal(min_dist, abs(side_dists(4)))) then + levelset%sf(i, j, k, ib_patch_id) = side_dists(4) + if (.not. f_approx_equal(side_dists(4), 0._wp)) then + dist_vec(2) = -side_dists(4)/abs(side_dists(4)) + end if - else if (f_approx_equal(min_dist, abs(side_dists(5)))) then - levelset%sf(i, j, k, ib_patch_id) = side_dists(5) - if (.not. f_approx_equal(side_dists(5), 0._wp)) then - levelset_norm%sf(i, j, k, ib_patch_id, 3) = side_dists(5)/ & - abs(side_dists(5)) - end if + else if (f_approx_equal(min_dist, abs(side_dists(5)))) then + levelset%sf(i, j, k, ib_patch_id) = side_dists(5) + if (.not. f_approx_equal(side_dists(5), 0._wp)) then + dist_vec(3) = side_dists(5)/abs(side_dists(5)) + end if - else if (f_approx_equal(min_dist, abs(side_dists(6)))) then - levelset%sf(i, j, k, ib_patch_id) = side_dists(6) - if (.not. f_approx_equal(side_dists(6), 0._wp)) then - levelset_norm%sf(i, j, k, ib_patch_id, 3) = -side_dists(6)/ & - abs(side_dists(6)) + else if (f_approx_equal(min_dist, abs(side_dists(6)))) then + levelset%sf(i, j, k, ib_patch_id) = side_dists(6) + if (.not. f_approx_equal(side_dists(6), 0._wp)) then + dist_vec(3) = -side_dists(6)/abs(side_dists(6)) + end if end if + levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, dist_vec) end if - levelset_norm%sf(i, j, 0, ib_patch_id, :) = & - matmul(rotation, levelset_norm%sf(i, j, 0, ib_patch_id, :)) - end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_cuboid_levelset @@ -449,33 +452,34 @@ contains integer, intent(IN) :: ib_patch_id real(wp) :: radius, dist - real(wp) :: x_centroid, y_centroid, z_centroid - real(wp), dimension(3) :: dist_vec + real(wp), dimension(3) :: dist_vec, center integer :: i, j, k !< Loop index variables radius = patch_ib(ib_patch_id)%radius - x_centroid = patch_ib(ib_patch_id)%x_centroid - y_centroid = patch_ib(ib_patch_id)%y_centroid - z_centroid = patch_ib(ib_patch_id)%z_centroid - - do i = 0, m - do j = 0, n - do k = 0, p - dist_vec(1) = x_cc(i) - x_centroid - dist_vec(2) = y_cc(j) - y_centroid - dist_vec(3) = z_cc(k) - z_centroid - dist = sqrt(sum(dist_vec**2)) - levelset%sf(i, j, k, ib_patch_id) = dist - radius - if (f_approx_equal(dist, 0._wp)) then - levelset_norm%sf(i, j, k, ib_patch_id, :) = (/1, 0, 0/) - else - levelset_norm%sf(i, j, k, ib_patch_id, :) = & - dist_vec(:)/dist - end if + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid + center(3) = patch_ib(ib_patch_id)%z_centroid + + #:call GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,dist]', & + & copyin='[ib_patch_id,center,radius]', collapse=3) + do i = 0, m + do j = 0, n + do k = 0, p + dist_vec(1) = x_cc(i) - center(1) + dist_vec(2) = y_cc(j) - center(2) + dist_vec(3) = z_cc(k) - center(3) + dist = sqrt(sum(dist_vec**2)) + levelset%sf(i, j, k, ib_patch_id) = dist - radius + if (f_approx_equal(dist, 0._wp)) then + levelset_norm%sf(i, j, k, ib_patch_id, :) = (/1, 0, 0/) + else + levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_vec(:)/dist + end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_sphere_levelset @@ -486,80 +490,77 @@ contains integer, intent(IN) :: ib_patch_id real(wp) :: radius - real(wp) :: x_centroid, y_centroid, z_centroid - real(wp) :: length_x, length_y, length_z - real(wp), dimension(3) :: centroid_vec, dist_sides_vec, dist_surface_vec + real(wp), dimension(3) :: dist_sides_vec, dist_surface_vec, length + real(wp), dimension(2) :: boundary real(wp) :: dist_side, dist_surface, side_pos - type(bounds_info) :: boundary integer :: i, j, k !< Loop index variables - real(wp), dimension(1:3) :: xyz_local !< x and y coordinates in local IB frame + real(wp), dimension(1:3) :: xyz_local, center !< x and y coordinates in local IB frame real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation radius = patch_ib(ib_patch_id)%radius - x_centroid = patch_ib(ib_patch_id)%x_centroid - y_centroid = patch_ib(ib_patch_id)%y_centroid - z_centroid = patch_ib(ib_patch_id)%z_centroid - length_x = patch_ib(ib_patch_id)%length_x - length_y = patch_ib(ib_patch_id)%length_y - length_z = patch_ib(ib_patch_id)%length_z + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid + center(3) = patch_ib(ib_patch_id)%z_centroid + length(1) = patch_ib(ib_patch_id)%length_x + length(2) = patch_ib(ib_patch_id)%length_y + length(3) = patch_ib(ib_patch_id)%length_z inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) - if (.not. f_approx_equal(length_x, 0._wp)) then - boundary%beg = -0.5_wp*length_x - boundary%end = 0.5_wp*length_x + if (.not. f_approx_equal(length(1), 0._wp)) then + boundary(1) = -0.5_wp*length(1) + boundary(2) = 0.5_wp*length(1) dist_sides_vec = (/1, 0, 0/) dist_surface_vec = (/0, 1, 1/) - else if (.not. f_approx_equal(length_y, 0._wp)) then - boundary%beg = -0.5_wp*length_y - boundary%end = 0.5_wp*length_y + else if (.not. f_approx_equal(length(2), 0._wp)) then + boundary(1) = -0.5_wp*length(2) + boundary(2) = 0.5_wp*length(2) dist_sides_vec = (/0, 1, 0/) dist_surface_vec = (/1, 0, 1/) - else if (.not. f_approx_equal(length_z, 0._wp)) then - boundary%beg = -0.5_wp*length_z - boundary%end = 0.5_wp*length_z + else if (.not. f_approx_equal(length(3), 0._wp)) then + boundary(1) = -0.5_wp*length(3) + boundary(2) = 0.5_wp*length(3) dist_sides_vec = (/0, 0, 1/) dist_surface_vec = (/1, 1, 0/) end if - do i = 0, m - do j = 0, n - do k = 0, p - xyz_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_cc(k) - z_centroid] ! get coordinate frame centered on IB - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates - - ! get distance to flat edge of cylinder - side_pos = dot_product(xyz_local, dist_sides_vec) - dist_side = min(abs(side_pos - boundary%beg), & - abs(boundary%end - side_pos)) - ! get distance to curved side of cylinder - dist_surface = norm2(xyz_local*dist_surface_vec) & - - radius - - if (dist_side < abs(dist_surface)) then - ! if the closest edge is flat - levelset%sf(i, j, k, ib_patch_id) = -dist_side - if (f_approx_equal(dist_side, abs(side_pos - boundary%beg))) then - levelset_norm%sf(i, j, k, ib_patch_id, :) = -dist_sides_vec + #:call GPU_PARALLEL_LOOP(private='[i,j,k,side_pos,dist_side,dist_surface,xyz_local]', & + & copyin='[ib_patch_id,center,radius,inverse_rotation,rotation,dist_sides_vec,dist_surface_vec]', collapse=3) + do i = 0, m + do j = 0, n + do k = 0, p + 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 + + ! get distance to flat edge of cylinder + side_pos = dot_product(xyz_local, dist_sides_vec) + dist_side = min(abs(side_pos - boundary(1)), & + abs(boundary(2) - side_pos)) + ! get distance to curved side of cylinder + dist_surface = norm2(xyz_local*dist_surface_vec) & + - radius + + if (dist_side < abs(dist_surface)) then + ! if the closest edge is flat + levelset%sf(i, j, k, ib_patch_id) = -dist_side + if (f_approx_equal(dist_side, abs(side_pos - boundary(1)))) then + levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, -dist_sides_vec) + else + levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, dist_sides_vec) + end if else - levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_sides_vec - end if - else - levelset%sf(i, j, k, ib_patch_id) = dist_surface + levelset%sf(i, j, k, ib_patch_id) = dist_surface - levelset_norm%sf(i, j, k, ib_patch_id, :) = & - xyz_local*dist_surface_vec - levelset_norm%sf(i, j, k, ib_patch_id, :) = & - levelset_norm%sf(i, j, k, ib_patch_id, :)/ & - norm2(levelset_norm%sf(i, j, k, ib_patch_id, :)) - end if - levelset_norm%sf(i, j, k, ib_patch_id, :) = & - matmul(rotation, levelset_norm%sf(i, j, k, ib_patch_id, :)) + xyz_local = xyz_local*dist_surface_vec + xyz_local = xyz_local/norm2(xyz_local) + levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, xyz_local) + end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_cylinder_levelset diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index de1b1f1908..c9f422b973 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -140,6 +140,7 @@ contains integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + real(wp), dimension(1:2) :: center real(wp) :: radius integer :: i, j, k !< Generic loop iterators @@ -147,8 +148,8 @@ contains ! Transferring the circular patch's radius, centroid, smearing patch ! identity and smearing coefficient information - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid radius = patch_ib(patch_id)%radius ! Initializing the pseudo volume fraction value to 1. The value will @@ -161,15 +162,18 @@ contains ! that cell. If both queries check out, the primitive variables of ! the current patch are assigned to this cell. - do j = 0, n - do i = 0, m - if ((x_cc(i) - x_centroid)**2 & - + (y_cc(j) - y_centroid)**2 <= radius**2) & - then - ib_markers_sf(i, j, 0) = patch_id - end if + #:call GPU_PARALLEL_LOOP(private='[i,j]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,radius]', collapse=2) + do j = 0, n + do i = 0, m + if ((x_cc(i) - center(1))**2 & + + (y_cc(j) - center(2))**2 <= radius**2) & + then + ib_markers_sf(i, j, 0) = patch_id + end if + end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_ib_circle @@ -186,10 +190,11 @@ contains integer :: Np1, Np2 real(wp), dimension(1:3) :: xy_local !< 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 - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid ca_in = patch_ib(patch_id)%c pa = patch_ib(patch_id)%p ma = patch_ib(patch_id)%m @@ -209,114 +214,118 @@ contains if (.not. allocated(airfoil_grid_u)) then allocate (airfoil_grid_u(1:Np)) allocate (airfoil_grid_l(1:Np)) - end if - - ! TODO :: The below instantiations are already handles by the loop below - airfoil_grid_u(1)%x = 0._wp - airfoil_grid_u(1)%y = 0._wp - airfoil_grid_l(1)%x = 0._wp - airfoil_grid_l(1)%y = 0._wp - - eta = 1._wp - - do i = 1, Np1 + Np2 - 1 - ! TODO :: This allocated the upper and lower airfoil arrays, and does not need to be performed each time the IB markers are updated. Place this as a separate subroutine. - if (i <= Np1) then - xc = i*(pa*ca_in/Np1) - xa = xc/ca_in - yc = (ma/pa**2)*(2*pa*xa - xa**2) - dycdxc = (2*ma/pa**2)*(pa - xa) - else - xc = pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) - xa = xc/ca_in - yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) - dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) - end if + ! TODO :: The below instantiations are already handles by the loop below + airfoil_grid_u(1)%x = 0._wp + airfoil_grid_u(1)%y = 0._wp + + airfoil_grid_l(1)%x = 0._wp + airfoil_grid_l(1)%y = 0._wp + + eta = 1._wp + + do i = 1, Np1 + Np2 - 1 + ! TODO :: This allocated the upper and lower airfoil arrays, and does not need to be performed each time the IB markers are updated. Place this as a separate subroutine. + if (i <= Np1) then + xc = i*(pa*ca_in/Np1) + xa = xc/ca_in + yc = (ma/pa**2)*(2*pa*xa - xa**2) + dycdxc = (2*ma/pa**2)*(pa - xa) + else + xc = pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) + xa = xc/ca_in + yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) + dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) + end if - yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) - sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp - cos_c = 1/(1 + dycdxc**2)**0.5_wp + yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) + sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp + cos_c = 1/(1 + dycdxc**2)**0.5_wp - xu = xa - yt*sin_c - yu = yc + yt*cos_c + xu = xa - yt*sin_c + yu = yc + yt*cos_c - xl = xa + yt*sin_c - yl = yc - yt*cos_c + xl = xa + yt*sin_c + yl = yc - yt*cos_c - xu = xu*ca_in - yu = yu*ca_in + xu = xu*ca_in + yu = yu*ca_in - xl = xl*ca_in - yl = yl*ca_in + xl = xl*ca_in + yl = yl*ca_in - airfoil_grid_u(i + 1)%x = xu - airfoil_grid_u(i + 1)%y = yu + airfoil_grid_u(i + 1)%x = xu + airfoil_grid_u(i + 1)%y = yu - airfoil_grid_l(i + 1)%x = xl - airfoil_grid_l(i + 1)%y = yl + airfoil_grid_l(i + 1)%x = xl + airfoil_grid_l(i + 1)%y = yl - end do + end do - airfoil_grid_u(Np)%x = ca_in - airfoil_grid_u(Np)%y = 0._wp + airfoil_grid_u(Np)%x = ca_in + airfoil_grid_u(Np)%y = 0._wp - airfoil_grid_l(Np)%x = ca_in - airfoil_grid_l(Np)%y = 0._wp + airfoil_grid_l(Np)%x = ca_in + airfoil_grid_l(Np)%y = 0._wp - do j = 0, n - do i = 0, m - xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp] ! get coordinate frame centered on IB - xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates + end if - 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 (airfoil_grid_u(k)%x < xy_local(1) .and. k <= Np) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_u(k)%x, xy_local(1))) then - if (xy_local(2) <= airfoil_grid_u(k)%y) then - !!IB - ib_markers_sf(i, j, 0) = patch_id - end if + #:call GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,f]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,inverse_rotation,ma,ca_in,airfoil_grid_u,airfoil_grid_l]', collapse=2) + do j = 0, n + do i = 0, m + 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 + + 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 - f = (airfoil_grid_u(k)%x - xy_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) - if (xy_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then - !!IB - ib_markers_sf(i, j, 0) = patch_id - end if + yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) + dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) end if - else - k = 1 - do while (airfoil_grid_l(k)%x < xy_local(1)) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_l(k)%x, xy_local(1))) then - if (xy_local(2) >= airfoil_grid_l(k)%y) then + if (xy_local(2) >= 0._wp) then + k = 1 + do while (airfoil_grid_u(k)%x < xy_local(1) .and. k <= Np) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_u(k)%x, xy_local(1))) then + if (xy_local(2) <= airfoil_grid_u(k)%y) then + !!IB + ib_markers_sf(i, j, 0) = patch_id + end if + else + f = (airfoil_grid_u(k)%x - xy_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) + if (xy_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then !!IB - ib_markers_sf(i, j, 0) = patch_id + ib_markers_sf(i, j, 0) = patch_id + end if end if else - f = (airfoil_grid_l(k)%x - xy_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) + k = 1 + do while (airfoil_grid_l(k)%x < xy_local(1)) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_l(k)%x, xy_local(1))) then + if (xy_local(2) >= airfoil_grid_l(k)%y) then + !!IB + ib_markers_sf(i, j, 0) = patch_id + end if + else + f = (airfoil_grid_l(k)%x - xy_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) - if (xy_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then + if (xy_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then !!IB - ib_markers_sf(i, j, 0) = patch_id + ib_markers_sf(i, j, 0) = patch_id + end if end if end if end if - end if + end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_ib_airfoil @@ -332,12 +341,12 @@ contains integer :: i, j, k, l integer :: Np1, Np2 - real(wp), dimension(1:3) :: xyz_local !< x, y, z coordinates in local IB frame + real(wp), dimension(1:3) :: xyz_local, center !< x, y, z coordinates in local IB frame real(wp), dimension(1:3, 1:3) :: inverse_rotation - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - z_centroid = patch_ib(patch_id)%z_centroid + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid + center(3) = patch_ib(patch_id)%z_centroid lz = patch_ib(patch_id)%length_z ca_in = patch_ib(patch_id)%c pa = patch_ib(patch_id)%p @@ -355,121 +364,118 @@ contains #endif Np = Np1 + Np2 + 1 - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) - - airfoil_grid_u(1)%x = 0._wp - airfoil_grid_u(1)%y = 0._wp - - airfoil_grid_l(1)%x = 0._wp - airfoil_grid_l(1)%y = 0._wp + if (.not. allocated(airfoil_grid_u)) then + allocate (airfoil_grid_u(1:Np)) + allocate (airfoil_grid_l(1:Np)) - z_max = lz/2 - z_min = -lz/2 + airfoil_grid_u(1)%x = 0._wp + airfoil_grid_u(1)%y = 0._wp - eta = 1._wp + airfoil_grid_l(1)%x = 0._wp + airfoil_grid_l(1)%y = 0._wp - do i = 1, Np1 + Np2 - 1 - if (i <= Np1) then - xc = i*(pa*ca_in/Np1) - xa = xc/ca_in - yc = (ma/pa**2)*(2*pa*xa - xa**2) - dycdxc = (2*ma/pa**2)*(pa - xa) - else - xc = pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) - xa = xc/ca_in - yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) - dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) - end if + z_max = lz/2 + z_min = -lz/2 - yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) - sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp - cos_c = 1/(1 + dycdxc**2)**0.5_wp + eta = 1._wp - xu = xa - yt*sin_c - yu = yc + yt*cos_c + do i = 1, Np1 + Np2 - 1 + if (i <= Np1) then + xc = i*(pa*ca_in/Np1) + xa = xc/ca_in + yc = (ma/pa**2)*(2*pa*xa - xa**2) + dycdxc = (2*ma/pa**2)*(pa - xa) + else + xc = pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) + xa = xc/ca_in + yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) + dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) + end if - xl = xa + yt*sin_c - yl = yc - yt*cos_c + yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4) + sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp + cos_c = 1/(1 + dycdxc**2)**0.5_wp - xu = xu*ca_in - yu = yu*ca_in + xu = xa - yt*sin_c + yu = yc + yt*cos_c - xl = xl*ca_in - yl = yl*ca_in + xl = xa + yt*sin_c + yl = yc - yt*cos_c - airfoil_grid_u(i + 1)%x = xu - airfoil_grid_u(i + 1)%y = yu + xu = xu*ca_in + yu = yu*ca_in - airfoil_grid_l(i + 1)%x = xl - airfoil_grid_l(i + 1)%y = yl + xl = xl*ca_in + yl = yl*ca_in - end do + airfoil_grid_u(i + 1)%x = xu + airfoil_grid_u(i + 1)%y = yu - airfoil_grid_u(Np)%x = ca_in - airfoil_grid_u(Np)%y = 0._wp + airfoil_grid_l(i + 1)%x = xl + airfoil_grid_l(i + 1)%y = yl - airfoil_grid_l(Np)%x = ca_in - airfoil_grid_l(Np)%y = 0._wp + end do - do l = 0, p - do j = 0, n - do i = 0, m - xyz_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_cc(l) - z_centroid] ! get coordinate frame centered on IB - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates + airfoil_grid_u(Np)%x = ca_in + airfoil_grid_u(Np)%y = 0._wp - if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then + airfoil_grid_l(Np)%x = ca_in + airfoil_grid_l(Np)%y = 0._wp + end if - if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then - xa = xyz_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 (xyz_local(2) >= 0._wp) then - k = 1 - do while (airfoil_grid_u(k)%x < xyz_local(1)) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_u(k)%x, xyz_local(1))) then - if (xyz_local(2) <= airfoil_grid_u(k)%y) then - !!IB - ib_markers_sf(i, j, l) = patch_id - end if - else - f = (airfoil_grid_u(k)%x - xyz_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) - if (xyz_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then + #:call GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,f]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,inverse_rotation,ma,ca_in,airfoil_grid_u,airfoil_grid_l]', collapse=3) + do l = 0, p + do j = 0, n + do i = 0, m + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB + xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates + + if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then + + if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then + if (xyz_local(2) >= 0._wp) then + k = 1 + do while (airfoil_grid_u(k)%x < xyz_local(1)) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_u(k)%x, xyz_local(1))) then + if (xyz_local(2) <= airfoil_grid_u(k)%y) then !!IB - ib_markers_sf(i, j, l) = patch_id - end if - end if - else - k = 1 - do while (airfoil_grid_l(k)%x < xyz_local(1)) - k = k + 1 - end do - if (f_approx_equal(airfoil_grid_l(k)%x, xyz_local(1))) then - if (xyz_local(2) >= airfoil_grid_l(k)%y) then + ib_markers_sf(i, j, l) = patch_id + end if + else + f = (airfoil_grid_u(k)%x - xyz_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x) + if (xyz_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then !!IB - ib_markers_sf(i, j, l) = patch_id + ib_markers_sf(i, j, l) = patch_id + end if end if else - f = (airfoil_grid_l(k)%x - xyz_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) + k = 1 + do while (airfoil_grid_l(k)%x < xyz_local(1)) + k = k + 1 + end do + if (f_approx_equal(airfoil_grid_l(k)%x, xyz_local(1))) then + if (xyz_local(2) >= airfoil_grid_l(k)%y) then + !!IB + ib_markers_sf(i, j, l) = patch_id + end if + else + f = (airfoil_grid_l(k)%x - xyz_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x) - if (xyz_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then + if (xyz_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then !!IB - ib_markers_sf(i, j, l) = patch_id + ib_markers_sf(i, j, l) = patch_id + end if end if end if end if end if - end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_ib_3D_airfoil @@ -492,6 +498,7 @@ contains integer :: i, j, k !< generic loop iterators real(wp) :: pi_inf, gamma, lit_gamma !< 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 pi_inf = fluid_pp(1)%pi_inf @@ -499,19 +506,12 @@ contains lit_gamma = (1._wp + gamma)/gamma ! Transferring the rectangle's centroid and length information - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - length_x = patch_ib(patch_id)%length_x - length_y = patch_ib(patch_id)%length_y + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid + 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(:, :) - ! Computing the beginning and the end x- and y-coordinates of the - ! rectangle based on its centroid and lengths - x_boundary%beg = -0.5_wp*length_x - x_boundary%end = 0.5_wp*length_x - y_boundary%beg = -0.5_wp*length_y - y_boundary%end = 0.5_wp*length_y - ! Since the rectangular patch does not allow for its boundaries to ! be smoothed out, the pseudo volume fraction is set to 1 to ensure ! that only the current patch contributes to the fluid state in the @@ -522,22 +522,26 @@ contains ! 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. - do j = 0, n - do i = 0, m - ! get the x and y coordinates in the local IB frame - xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp] - xy_local = matmul(inverse_rotation, xy_local) - if (x_boundary%beg <= xy_local(1) .and. & - x_boundary%end >= xy_local(1) .and. & - y_boundary%beg <= xy_local(2) .and. & - y_boundary%end >= xy_local(2)) then - - ! Updating the patch identities bookkeeping variable - ib_markers_sf(i, j, 0) = patch_id + #:call GPU_PARALLEL_LOOP(private='[i,j, xy_local]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,length,inverse_rotation,x_cc,y_cc]', collapse=2) + do j = 0, n + do i = 0, m + ! 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) - end if + 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 + + ! Updating the patch identities bookkeeping variable + ib_markers_sf(i, j, 0) = patch_id + + end if + end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_ib_rectangle @@ -557,15 +561,16 @@ contains ! Generic loop iterators integer :: i, j, k real(wp) :: radius + real(wp), dimension(1:3) :: center !! Variables to initialize the pressure field that corresponds to the !! bubble-collapse test case found in Tiwari et al. (2013) ! Transferring spherical patch's radius, centroid, smoothing patch ! identity and smoothing coefficient information - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - z_centroid = patch_ib(patch_id)%z_centroid + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid + center(3) = patch_ib(patch_id)%z_centroid radius = patch_ib(patch_id)%radius ! Initializing the pseudo volume fraction value to 1. The value will @@ -577,24 +582,27 @@ contains ! and verifying whether the current patch has permission to write to ! that cell. If both queries check out, the primitive variables of ! the current patch are assigned to this cell. - do k = 0, p - do j = 0, n - do i = 0, m - if (grid_geometry == 3) then - call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) - else - cart_y = y_cc(j) - cart_z = z_cc(k) - end if - ! Updating the patch identities bookkeeping variable - if (((x_cc(i) - x_centroid)**2 & - + (cart_y - y_centroid)**2 & - + (cart_z - z_centroid)**2 <= radius**2)) then - ib_markers_sf(i, j, k) = patch_id - end if + #:call GPU_PARALLEL_LOOP(private='[i,j,k,cart_y,cart_z]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,radius]', collapse=3) + do k = 0, p + do j = 0, n + do i = 0, m + if (grid_geometry == 3) then + call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) + else + cart_y = y_cc(j) + cart_z = z_cc(k) + end if + ! Updating the patch identities bookkeeping variable + if (((x_cc(i) - center(1))**2 & + + (cart_y - center(2))**2 & + + (cart_z - center(3))**2 <= radius**2)) then + ib_markers_sf(i, j, k) = patch_id + end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_ib_sphere @@ -614,27 +622,18 @@ contains integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf integer :: i, j, k !< Generic loop iterators - real(wp), dimension(1:3) :: xyz_local !< x and y coordinates in local IB frame + 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 ! Transferring the cuboid's centroid and length information - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - z_centroid = patch_ib(patch_id)%z_centroid - length_x = patch_ib(patch_id)%length_x - length_y = patch_ib(patch_id)%length_y - length_z = patch_ib(patch_id)%length_z + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid + center(3) = patch_ib(patch_id)%z_centroid + 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(:, :) - ! Computing the beginning and the end x-, y- and z-coordinates of - ! the cuboid based on its centroid and lengths - x_boundary%beg = -0.5_wp*length_x - x_boundary%end = 0.5_wp*length_x - y_boundary%beg = -0.5_wp*length_y - y_boundary%end = 0.5_wp*length_y - z_boundary%beg = -0.5_wp*length_z - z_boundary%end = 0.5_wp*length_z - ! Since the cuboidal patch does not allow for its boundaries to get ! smoothed out, the pseudo volume fraction is set to 1 to make sure ! that only the current patch contributes to the fluid state in the @@ -645,33 +644,36 @@ contains ! 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. - do k = 0, p - do j = 0, n - do i = 0, m - - if (grid_geometry == 3) then - ! TODO :: This does not work and is not covered by any tests. This should be fixed - call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) - else - cart_y = y_cc(j) - cart_z = z_cc(k) - end if - xyz_local = [x_cc(i) - x_centroid, cart_y - y_centroid, cart_z - z_centroid] ! get coordinate frame centered on IB - xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates - - if (x_boundary%beg <= xyz_local(1) .and. & - x_boundary%end >= xyz_local(1) .and. & - y_boundary%beg <= xyz_local(2) .and. & - y_boundary%end >= xyz_local(2) .and. & - z_boundary%beg <= xyz_local(3) .and. & - z_boundary%end >= xyz_local(3)) then - - ! Updating the patch identities bookkeeping variable - ib_markers_sf(i, j, k) = patch_id - end if + #:call GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local,cart_y,cart_z]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,length,inverse_rotation]', collapse=3) + do k = 0, p + do j = 0, n + do i = 0, m + + if (grid_geometry == 3) then + ! TODO :: This does not work and is not covered by any tests. This should be fixed + call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) + else + cart_y = y_cc(j) + cart_z = z_cc(k) + end if + xyz_local = [x_cc(i), cart_y, cart_z] - 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 + + ! Updating the patch identities bookkeeping variable + ib_markers_sf(i, j, k) = patch_id + end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_ib_cuboid @@ -693,30 +695,21 @@ contains integer :: i, j, k !< Generic loop iterators real(wp) :: radius - real(wp), dimension(1:3) :: xyz_local !< x and y coordinates in local IB frame + 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 ! Transferring the cylindrical patch's centroid, length, radius, ! smoothing patch identity and smoothing coefficient information - x_centroid = patch_ib(patch_id)%x_centroid - y_centroid = patch_ib(patch_id)%y_centroid - z_centroid = patch_ib(patch_id)%z_centroid - length_x = patch_ib(patch_id)%length_x - length_y = patch_ib(patch_id)%length_y - length_z = patch_ib(patch_id)%length_z + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid + center(3) = patch_ib(patch_id)%z_centroid + 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(:, :) - ! Computing the beginning and the end x-, y- and z-coordinates of - ! the cylinder based on its centroid and lengths - x_boundary%beg = -0.5_wp*length_x - x_boundary%end = 0.5_wp*length_x - y_boundary%beg = -0.5_wp*length_y - y_boundary%end = 0.5_wp*length_y - z_boundary%beg = -0.5_wp*length_z - z_boundary%end = 0.5_wp*length_z - ! Initializing the pseudo volume fraction value to 1. The value will ! be modified as the patch is laid out on the grid, but only in the ! case that smearing of the cylindrical patch's boundary is enabled. @@ -726,43 +719,46 @@ contains ! 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. - do k = 0, p - do j = 0, n - do i = 0, m - - if (grid_geometry == 3) then - call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) - else - cart_y = y_cc(j) - cart_z = z_cc(k) - end if - xyz_local = [x_cc(i) - x_centroid, cart_y - y_centroid, cart_z - z_centroid] ! 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_x) .and. & - xyz_local(2)**2 & - + xyz_local(3)**2 <= radius**2 .and. & - x_boundary%beg <= xyz_local(1) .and. & - x_boundary%end >= xyz_local(1)) & - .or. & - (.not. f_is_default(length_y) .and. & - xyz_local(1)**2 & - + xyz_local(3)**2 <= radius**2 .and. & - y_boundary%beg <= xyz_local(2) .and. & - y_boundary%end >= xyz_local(2)) & - .or. & - (.not. f_is_default(length_z) .and. & - xyz_local(1)**2 & - + xyz_local(2)**2 <= radius**2 .and. & - z_boundary%beg <= xyz_local(3) .and. & - z_boundary%end >= xyz_local(3)))) then - - ! Updating the patch identities bookkeeping variable - ib_markers_sf(i, j, k) = patch_id - end if + #:call GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local,cart_y,cart_z]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,length,radius,inverse_rotation]', collapse=3) + do k = 0, p + do j = 0, n + do i = 0, m + + if (grid_geometry == 3) then + call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k)) + else + cart_y = y_cc(j) + cart_z = z_cc(k) + end if + xyz_local = [x_cc(i), cart_y, cart_z] - 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 + + ! Updating the patch identities bookkeeping variable + ib_markers_sf(i, j, k) = patch_id + end if + end do end do end do - end do + #:endcall GPU_PARALLEL_LOOP end subroutine s_ib_cylinder diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index de1517e079..f82cca7466 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -145,17 +145,20 @@ end subroutine s_setup_mpi_io_params !> Helper subroutine to read IB data files !! @param file_loc_base Base file location for IB data - impure subroutine s_read_ib_data_files(file_loc_base) + impure subroutine s_read_ib_data_files(file_loc_base, t_step) character(len=*), intent(in) :: file_loc_base + integer, intent(in), optional :: t_step character(LEN=len_trim(file_loc_base) + 20) :: file_loc logical :: file_exist - integer :: ifile, ierr, data_size + integer :: ifile, ierr, data_size, var_MOK #ifdef MFC_MPI integer, dimension(MPI_STATUS_SIZE) :: status integer(KIND=MPI_OFFSET_KIND) :: disp + integer :: m_MOK, n_MOK, p_MOK, MOK, WP_MOK, save_index + #endif if (.not. ib) return @@ -171,8 +174,20 @@ impure subroutine s_read_ib_data_files(file_loc_base) #ifdef MFC_MPI call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) + n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) + p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) + MOK = int(1._wp, MPI_OFFSET_KIND) + WP_MOK = int(8._wp, MPI_OFFSET_KIND) + save_index = t_step/t_step_save ! get the number of saves done to this point + data_size = (m + 1)*(n + 1)*(p + 1) - disp = 0 + var_MOK = int(sys_size + 1, MPI_OFFSET_KIND) + if (t_step == 0) then + disp = 0 + else + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1 + int(save_index, MPI_OFFSET_KIND)) + end if call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & 'native', mpi_info_int, ierr) @@ -253,9 +268,6 @@ impure subroutine s_read_serial_data_files(t_step) write (t_step_dir, '(A,I0,A,I0)') '/p_all/p', proc_rank, '/', t_step t_step_dir = trim(case_dir)//trim(t_step_dir) - write (t_step_ib_dir, '(A,I0,A,I0)') '/p_all/p', proc_rank, '/', 0 - t_step_ib_dir = trim(case_dir)//trim(t_step_ib_dir) - ! Inquiring as to the existence of the time-step directory file_loc = trim(t_step_dir)//'/.' call my_inquire(file_loc, dir_check) @@ -266,15 +278,6 @@ impure subroutine s_read_serial_data_files(t_step) ' is missing. Exiting.') end if - file_loc = trim(t_step_ib_dir)//'/.' - call my_inquire(file_loc, dir_check) - - ! If the time-step directory is missing, the post-process exits. - if (dir_check .neqv. .true.) then - call s_mpi_abort('Time-step folder '//trim(t_step_ib_dir)// & - ' is missing. Exiting.') - end if - if (bc_io) then call s_read_serial_boundary_condition_files(t_step_dir, bc_type) else @@ -317,7 +320,7 @@ impure subroutine s_read_serial_data_files(t_step) end do ! Reading IB data using helper subroutine - call s_read_ib_data_files(t_step_ib_dir) + call s_read_ib_data_files(t_step_dir) end subroutine s_read_serial_data_files @@ -542,7 +545,7 @@ impure subroutine s_read_parallel_conservative_data(t_step, m_MOK, n_MOK, p_MOK, end do end if - call s_read_ib_data_files(trim(case_dir)//'/restart_data'//trim(mpiiofs)) + call s_read_ib_data_files(trim(case_dir)//'/restart_data'//trim(mpiiofs), t_step) else call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting.') end if @@ -573,7 +576,7 @@ impure subroutine s_read_parallel_conservative_data(t_step, m_MOK, n_MOK, p_MOK, call s_mpi_barrier() call MPI_FILE_CLOSE(ifile, ierr) - call s_read_ib_data_files(trim(case_dir)//'/restart_data'//trim(mpiiofs)) + call s_read_ib_data_files(trim(case_dir)//'/restart_data'//trim(mpiiofs), t_step) else call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting.') end if diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index e5e6dc9b4f..7ae10f4408 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -503,7 +503,7 @@ contains FORM='unformatted', & STATUS='new') - write (2) ib_markers%sf; close (2) + write (2) ib_markers%sf(0:m, 0:n, 0:p); close (2) end if gamma = fluid_pp(1)%gamma @@ -1013,6 +1013,24 @@ contains end if call MPI_FILE_CLOSE(ifile, ierr) + + !Write ib data + if (ib) then + write (file_loc, '(A)') 'ib.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & + mpi_info_int, ifile, ierr) + + var_MOK = int(sys_size + 1, MPI_OFFSET_KIND) + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1 + int(t_step/t_step_save)) + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & + 'native', mpi_info_int, ierr) + call MPI_FILE_WRITE_ALL(ifile, MPI_IO_IB_DATA%var%sf, data_size, & + MPI_INTEGER, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + end if + end if #endif diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 760bdda8d8..1a159c74d6 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -121,9 +121,10 @@ contains call s_mpi_allreduce_integer_sum(num_gps, max_num_gps) call s_mpi_allreduce_integer_sum(num_inner_gps, max_num_inner_gps) + ! set the size of the ghost point arrays to be the amount of points total, plus a factor of 2 buffer $:GPU_UPDATE(device='[num_gps, num_inner_gps]') - @:ALLOCATE(ghost_points(1:int(max_num_gps * 2.0))) - @:ALLOCATE(inner_points(1:int(max_num_inner_gps * 2.0))) + @:ALLOCATE(ghost_points(1:int((max_num_gps + max_num_inner_gps) * 2.0))) + @:ALLOCATE(inner_points(1:int((max_num_gps + max_num_inner_gps) * 2.0))) $:GPU_ENTER_DATA(copyin='[ghost_points,inner_points]') @@ -262,6 +263,15 @@ contains norm = norm/buf vel_norm_IP = sum(vel_IP*norm)*norm vel_g = vel_IP - vel_norm_IP + if (patch_ib(patch_id)%moving_ibm /= 0) then + ! compute the linear velocity of the ghost point due to rotation + radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, & + patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid] + rotation_velocity = cross_product(matmul(patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector) + + ! add only the component of the IB's motion that is normal to the surface + vel_g = vel_g + sum((patch_ib(patch_id)%vel + rotation_velocity)*norm)*norm + end if else if (patch_ib(patch_id)%moving_ibm == 0) then ! we know the object is not moving if moving_ibm is 0 (false) @@ -922,7 +932,9 @@ contains integer :: i, ierr ! Clears the existing immersed boundary indices - ib_markers%sf = 0 + ib_markers%sf = 0._wp + levelset%sf = 0._wp + levelset_norm%sf = 0._wp ! recalulcate the rotation matrix based upon the new angles do i = 1, num_ibs @@ -936,7 +948,8 @@ contains ! recompute the new ib_patch locations and broadcast them. call s_apply_ib_patches(ib_markers%sf(0:m, 0:n, 0:p), levelset, levelset_norm) call s_populate_ib_buffers() ! transmits the new IB markers via MPI - $:GPU_UPDATE(device='[ib_markers%sf, levelset%sf, levelset_norm%sf]') + $:GPU_UPDATE(device='[ib_markers%sf]') + $:GPU_UPDATE(host='[levelset%sf, levelset_norm%sf]') ! recalculate the ghost point locations and coefficients call s_find_num_ghost_points(num_gps, num_inner_gps) diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 5ecc2c707a..d578eb7ff7 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -202,8 +202,8 @@ contains end do do i = 1, num_ibs - #:for VAR in [ 'radius', 'length_x', 'length_y', & - & 'x_centroid', 'y_centroid', 'c', 'm', 'p', 't', 'theta', 'slip'] + #:for VAR in [ 'radius', 'length_x', 'length_y', 'length_z', & + & 'x_centroid', 'y_centroid', 'z_centroid', 'c', 'm', 'p', 't', 'theta', 'slip'] call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor #:for VAR in ['vel', 'angular_vel', 'angles']