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 1b41bfc62c..fc632fa691 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -456,9 +456,6 @@ contains 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) @@ -468,21 +465,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 @@ -556,19 +554,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) @@ -580,7 +573,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) @@ -588,17 +582,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 @@ -617,21 +614,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) @@ -642,7 +634,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) @@ -650,20 +642,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 @@ -683,40 +678,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 @@ -895,6 +887,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 67e98193bd..fd1a006145 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -76,14 +76,17 @@ contains call nvtxStartRange("SETUP-IBM-MODULE") + $: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() ! allocate some arrays for MPI communication, if required by this simulation #ifdef MFC_MPI @@ -868,13 +871,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 +936,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 @@ -1122,6 +1125,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 +1163,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 +1184,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