From 450ff9fa995ba3cba9d94fbf8fdb432cc9847f76 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Fri, 17 Oct 2025 11:46:25 -0400 Subject: [PATCH 01/20] Replaced ib step dir to always read the correct IB files for the time step --- src/post_process/m_data_input.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index de1517e079..eed7a3ed34 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -253,7 +253,7 @@ 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 + write (t_step_ib_dir, '(A,I0,A,I0)') '/p_all/p', proc_rank, '/', t_step t_step_ib_dir = trim(case_dir)//trim(t_step_ib_dir) ! Inquiring as to the existence of the time-step directory From 1257732f44617f0eef30a04f72459e84a6b711a3 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Fri, 17 Oct 2025 12:55:59 -0400 Subject: [PATCH 02/20] Resolved serial IO issue: --- src/post_process/m_data_input.f90 | 14 +------------- src/simulation/m_data_output.fpp | 2 +- src/simulation/m_ibm.fpp | 5 +++-- 3 files changed, 5 insertions(+), 16 deletions(-) diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index eed7a3ed34..c149774045 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -253,9 +253,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, '/', t_step - 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 +263,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 +305,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 diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 5f5b79c481..682574c1ab 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -509,7 +509,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 diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 6fe4932574..6fd7bd689e 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -118,9 +118,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]') From d6e197611fa92b0583541f44af4934cef493936d Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Mon, 20 Oct 2025 13:36:05 -0400 Subject: [PATCH 03/20] Successfully writing out markers with parallel io except for initial step --- src/post_process/m_data_input.f90 | 21 ++++++++++++++++----- src/simulation/m_data_output.fpp | 19 +++++++++++++++++++ src/simulation/m_ibm.fpp | 9 +++++++++ 3 files changed, 44 insertions(+), 5 deletions(-) diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index c149774045..4588387254 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,16 @@ 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) + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1 + int(save_index, MPI_OFFSET_KIND)) call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & 'native', mpi_info_int, ierr) @@ -530,7 +541,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 @@ -561,7 +572,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 682574c1ab..442a401240 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1019,6 +1019,25 @@ contains end if call MPI_FILE_CLOSE(ifile, ierr) + + !Writeib data + if (ib) then + write (file_loc, '(A)') 'ib.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + print *, "Writing file ", 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 6fd7bd689e..33042798f9 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -266,6 +266,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) From a24ff394055f5ee1d861c3f782d80b9806678a1a Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Mon, 20 Oct 2025 13:44:55 -0400 Subject: [PATCH 04/20] IO is working --- src/simulation/m_data_output.fpp | 1 - src/simulation/p_main.fpp | 3 +++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 442a401240..b578acfbba 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1024,7 +1024,6 @@ contains if (ib) then write (file_loc, '(A)') 'ib.dat' file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - print *, "Writing file ", file_loc call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & mpi_info_int, ifile, ierr) diff --git a/src/simulation/p_main.fpp b/src/simulation/p_main.fpp index 29cb3b8281..3eefefd88f 100644 --- a/src/simulation/p_main.fpp +++ b/src/simulation/p_main.fpp @@ -70,6 +70,9 @@ program p_main call nvtxEndRange ! INIT + ! save the data at t_step=0 + call s_save_data(t_step, start, finish, io_time_avg, nt) + call nvtxStartRange("SIMULATION-TIME-MARCH") ! Time-stepping Loop do From 64a57bb777c27f6579f2f6b67389a3579b58d317 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Tue, 21 Oct 2025 17:42:28 -0400 Subject: [PATCH 05/20] mostly working, but original levelset not changing --- src/common/m_ib_patches.fpp | 19 ++++++++++++++----- src/simulation/m_ibm.fpp | 2 +- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index de1b1f1908..d3c4a78aff 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -492,6 +492,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 @@ -512,6 +513,11 @@ contains y_boundary%beg = -0.5_wp*length_y y_boundary%end = 0.5_wp*length_y + length(1) = length_x + length(2) = length_y + center(1) = x_centroid + center(2) = y_centroid + ! 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,15 +528,18 @@ 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. + $:GPU_PARALLEL_LOOP(private='[i,j, xy_local]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,length,inverse_rotation]', 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) - x_centroid, y_cc(j) - y_centroid, 0._wp] + xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 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 + + 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 diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 33042798f9..1a65607fd6 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -941,7 +941,7 @@ 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='[levelset%sf, levelset_norm%sf]') ! recalculate the ghost point locations and coefficients call s_find_num_ghost_points(num_gps, num_inner_gps) From 71dc5567b173701c4c71a7e4d1a5309a3880a8a8 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Tue, 21 Oct 2025 17:45:37 -0400 Subject: [PATCH 06/20] Back to working rectangle case with moving IBs --- src/simulation/m_ibm.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 1a65607fd6..33042798f9 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -941,7 +941,7 @@ 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='[levelset%sf, levelset_norm%sf]') + $:GPU_UPDATE(device='[ib_markers%sf, levelset%sf, levelset_norm%sf]') ! recalculate the ghost point locations and coefficients call s_find_num_ghost_points(num_gps, num_inner_gps) From f02636d456dd029d220bf101d451037215f299e0 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Wed, 22 Oct 2025 11:20:18 -0400 Subject: [PATCH 07/20] Finished porting all IB marker compute to the GPU --- src/common/m_ib_patches.fpp | 296 ++++++++++++++++++------------------ 1 file changed, 148 insertions(+), 148 deletions(-) diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index d3c4a78aff..1245d09221 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -139,7 +139,8 @@ 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,10 +162,12 @@ contains ! 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]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,radius]', collapse=2) do j = 0, n do i = 0, m - if ((x_cc(i) - x_centroid)**2 & - + (y_cc(j) - y_centroid)**2 <= radius**2) & + 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 @@ -186,10 +189,11 @@ contains integer :: Np1, Np2 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame + real(wp), dimension(1:3) :: 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,64 +213,67 @@ 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 + ! 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 - 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 + end if + $: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) - x_centroid, y_cc(j) - y_centroid, 0._wp] ! get coordinate frame centered on IB + 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 @@ -332,12 +339,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,67 +362,71 @@ contains #endif Np = Np1 + Np2 + 1 - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) + if (.not. allocated(airfoil_grid_u)) then + 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_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 + airfoil_grid_l(1)%x = 0._wp + airfoil_grid_l(1)%y = 0._wp - z_max = lz/2 - z_min = -lz/2 + z_max = lz/2 + z_min = -lz/2 - eta = 1._wp + eta = 1._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 + 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 - 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 + end if + $: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) - x_centroid, y_cc(j) - y_centroid, z_cc(l) - z_centroid] ! get coordinate frame centered on IB + 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 @@ -562,6 +573,7 @@ contains integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + integer, dimension(1:3):: center ! Generic loop iterators integer :: i, j, k @@ -572,9 +584,9 @@ contains ! 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 @@ -586,6 +598,8 @@ 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. + $:GPU_PARALLEL_LOOP(private='[i,j,k]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,radius,inverse_rotation]', collapse=3) do k = 0, p do j = 0, n do i = 0, m @@ -596,9 +610,9 @@ contains 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 + 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 @@ -623,27 +637,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 @@ -654,6 +659,8 @@ 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. + $:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local]', 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 @@ -665,15 +672,15 @@ contains 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 = [x_cc(i) - center(1), cart_y - center(2), cart_z - center(3)] ! 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 + if (-0.5length(1) <= xyz_local(1) .and. & + 0.5length(1) >= xyz_local(1) .and. & + -0.5length(2) <= xyz_local(2) .and. & + 0.5length(2) >= xyz_local(2) .and. & + -0.5length(3) <= xyz_local(3) .and. & + 0.5length(3) >= xyz_local(3)) then ! Updating the patch identities bookkeeping variable ib_markers_sf(i, j, k) = patch_id @@ -702,30 +709,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. @@ -735,6 +733,8 @@ 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. + $:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local]', 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 @@ -745,26 +745,26 @@ contains 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 = [x_cc(i) - center(1), cart_y - center(2), cart_z - center(3)] ! 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)) & + -0.5*length(1) <= xyz_local(1) .and. & + 0.5*length(1) >= 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)) & + -0.5*length(2) <= xyz_local(2) .and. & + 0.5*length(2) >= 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 + -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 From 8965ae288b4003419655cc2157b8565f0d5b9b00 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Wed, 22 Oct 2025 11:30:54 -0400 Subject: [PATCH 08/20] Full GPU support working for airfoils --- src/common/m_compute_levelset.fpp | 10 ++++++---- src/common/m_ib_patches.fpp | 16 ++++++++-------- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index cc473ea544..df651574ba 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -71,22 +71,24 @@ 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(:, :) + $:GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,dist_vec,dist,global_dist,global_id]', copy='[levelset,levelset_norm]',& + & 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) - x_centroid, y_cc(j) - y_centroid, 0._wp] ! get coordinate frame centered on IB + 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 diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index 1245d09221..df6464e676 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -189,7 +189,7 @@ contains integer :: Np1, Np2 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame - real(wp), dimension(1:3) :: center !< 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 center(1) = patch_ib(patch_id)%x_centroid @@ -599,7 +599,7 @@ contains ! 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]', copy='[ib_markers_sf]',& - & copyin='[patch_id,center,radius,inverse_rotation]', collapse=3) + & copyin='[patch_id,center,radius]', collapse=3) do k = 0, p do j = 0, n do i = 0, m @@ -675,12 +675,12 @@ contains xyz_local = [x_cc(i) - center(1), cart_y - center(2), cart_z - center(3)] ! get coordinate frame centered on IB xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates - if (-0.5length(1) <= xyz_local(1) .and. & - 0.5length(1) >= xyz_local(1) .and. & - -0.5length(2) <= xyz_local(2) .and. & - 0.5length(2) >= xyz_local(2) .and. & - -0.5length(3) <= xyz_local(3) .and. & - 0.5length(3) >= xyz_local(3)) then + 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 From e87fbff90d3de6b73eef4956d9dc4086493abd31 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Wed, 22 Oct 2025 11:37:13 -0400 Subject: [PATCH 09/20] Removed depricated code --- src/common/m_ib_patches.fpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index df6464e676..690a4b018d 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -432,14 +432,6 @@ contains 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 - 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)) From 240c58867971e1fe0a6b6c1c1b8efa7020aa3196 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Fri, 24 Oct 2025 08:18:39 -0400 Subject: [PATCH 10/20] Removed some unused variables and parallelized more levelset calculations --- src/common/m_compute_levelset.fpp | 64 +++++++++++++++++-------------- 1 file changed, 35 insertions(+), 29 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index df651574ba..daa0599a3b 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -34,21 +34,23 @@ 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 + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid + $:GPU_PARALLEL_LOOP(private='[i,j,dist_vec,dist]', copy='[levelset,levelset_norm]',& + & copyin='[ib_patch_id,center,radius]', collapse=2) 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_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 @@ -156,31 +158,33 @@ 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 + z_max = center(3) + lz/2 + z_min = center(3) - lz/2 + $:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,dist_vec,dist,global_dist,global_id,dist_side,dist_surf]', copy='[levelset,levelset_norm]',& + & 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) - x_centroid, y_cc(j) - y_centroid, z_cc(l) - z_centroid] ! get coordinate frame centered on IB + 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) >= y_centroid) then @@ -263,9 +267,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(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 @@ -273,8 +277,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(:, :) @@ -283,9 +287,11 @@ contains bottom_left(1) = -length_x/2 bottom_left(2) = -length_y/2 + $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local]', copy='[levelset,levelset_norm]',& + & copyin='[ib_patch_id,center,bottom_left,top_right,initial_distance_buffer,inverse_rotation,rotation]', collapse=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 = [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. & @@ -336,41 +342,41 @@ 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, 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 + $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local]', copy='[levelset,levelset_norm]',& + & 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) - x_centroid, y_cc(j) - y_centroid, z_cc(k) - z_centroid] ! get coordinate frame centered on IB + 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. & @@ -390,8 +396,8 @@ contains ! 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. + ! and this code should be cleaned up. It also means that + ! rotating the box 90 degrees will cause tests to fail. 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) From 77c14321ac7e7ef01f77e2fee3a0b49db0789624 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Fri, 24 Oct 2025 08:36:32 -0400 Subject: [PATCH 11/20] All levelsets updated for GPU compute --- src/common/m_compute_levelset.fpp | 59 ++++++++++++++++--------------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index daa0599a3b..9aebf92455 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -457,22 +457,23 @@ 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 + 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 + $:GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,dist]', copy='[levelset,levelset_norm]',& + & 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) - x_centroid - dist_vec(2) = y_cc(j) - y_centroid - dist_vec(3) = z_cc(k) - z_centroid + 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 @@ -494,54 +495,54 @@ 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) :: centroid_vec, 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 + 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 + 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 + 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 + $:GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset,levelset_norm]',& + & 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) - x_centroid, y_cc(j) - y_centroid, z_cc(k) - z_centroid] ! get coordinate frame centered on IB + 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%beg), & - abs(boundary%end - side_pos)) + 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 @@ -549,7 +550,7 @@ contains 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 + if (f_approx_equal(dist_side, abs(side_pos - boundary(1)))) then levelset_norm%sf(i, j, k, ib_patch_id, :) = -dist_sides_vec else levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_sides_vec From bd6096b87e5b0967338eeeb24919d10ffba2e430 Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Fri, 24 Oct 2025 13:27:07 -0400 Subject: [PATCH 12/20] Passed test suite of IBM cases --- src/common/m_compute_levelset.fpp | 71 ++++++++++++------------------- src/common/m_ib_patches.fpp | 30 ++++++------- 2 files changed, 43 insertions(+), 58 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index 9aebf92455..fa85808713 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -187,7 +187,7 @@ contains 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) >= y_centroid) then + 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 @@ -268,7 +268,7 @@ contains real(wp) :: side_dists(4) 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 @@ -287,7 +287,7 @@ contains bottom_left(1) = -length_x/2 bottom_left(2) = -length_y/2 - $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', copy='[levelset,levelset_norm]',& & copyin='[ib_patch_id,center,bottom_left,top_right,initial_distance_buffer,inverse_rotation,rotation]', collapse=2) do i = 0, m do j = 0, n @@ -312,20 +312,17 @@ contains end do levelset%sf(i, j, 0, ib_patch_id) = side_dists(idx) - levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp + 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 - levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(idx)/ & - abs(side_dists(idx)) + dist_vec(1) = side_dists(idx)/abs(side_dists(idx)) else ! vector points along the y axis - levelset_norm%sf(i, j, 0, ib_patch_id, 2) = side_dists(idx)/ & - abs(side_dists(idx)) + 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, levelset_norm%sf(i, j, 0, ib_patch_id, :)) + levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, dist_vec) else levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0._wp end if @@ -347,7 +344,7 @@ contains 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 !< Loop index variables @@ -370,7 +367,7 @@ contains Front = length_z/2 Back = -length_z/2 - $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local,dist_vec]', copy='[levelset,levelset_norm]',& & copyin='[ib_patch_id,center,inverse_rotation,rotation,Right,Left,Top,Bottom,Front,Back]', collapse=3) do i = 0, m do j = 0, n @@ -398,51 +395,44 @@ contains ! 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. - levelset_norm%sf(i, j, k, ib_patch_id, :) = 0._wp + 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 - levelset_norm%sf(i, j, k, ib_patch_id, 1) = side_dists(1)/ & - abs(side_dists(1)) + 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)) + 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)) + 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)) + 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)) + 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)) + dist_vec(3) = -side_dists(6)/abs(side_dists(6)) end if end if - levelset_norm%sf(i, j, 0, ib_patch_id, :) = & - matmul(rotation, levelset_norm%sf(i, j, 0, ib_patch_id, :)) + levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, dist_vec) end if end do end do @@ -479,8 +469,7 @@ contains 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 + levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_vec(:)/dist end if end do end do @@ -495,7 +484,7 @@ contains integer, intent(IN) :: ib_patch_id real(wp) :: radius - real(wp), dimension(3) :: centroid_vec, dist_sides_vec, dist_surface_vec, length + real(wp), dimension(3) :: dist_sides_vec, dist_surface_vec, length real(wp), dimension(2) :: boundary real(wp) :: dist_side, dist_surface, side_pos integer :: i, j, k !< Loop index variables @@ -514,24 +503,24 @@ contains 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 + 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 + 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 + 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 - $:GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,k,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset,levelset_norm]',& & 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 @@ -551,21 +540,17 @@ contains ! 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, :) = -dist_sides_vec + levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, -dist_sides_vec) else - levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_sides_vec + levelset_norm%sf(i, j, k, ib_patch_id, :) = matmul(rotation, dist_sides_vec) end if else 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, :)) + 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 - levelset_norm%sf(i, j, k, ib_patch_id, :) = & - matmul(rotation, levelset_norm%sf(i, j, k, ib_patch_id, :)) end do end do end do diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index 690a4b018d..9e863137ca 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -565,11 +565,11 @@ contains integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf - integer, dimension(1:3):: center ! 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) @@ -590,7 +590,7 @@ 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. - $:GPU_PARALLEL_LOOP(private='[i,j,k]', copy='[ib_markers_sf]',& + $: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 @@ -651,7 +651,7 @@ 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. - $:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local]', copy='[ib_markers_sf]',& + $: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 @@ -664,7 +664,7 @@ contains cart_y = y_cc(j) cart_z = z_cc(k) end if - xyz_local = [x_cc(i) - center(1), cart_y - center(2), cart_z - center(3)] ! get coordinate frame centered on IB + 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. & @@ -725,7 +725,7 @@ 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. - $:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local]', copy='[ib_markers_sf]',& + $: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 @@ -737,26 +737,26 @@ contains cart_y = y_cc(j) cart_z = z_cc(k) end if - xyz_local = [x_cc(i) - center(1), cart_y - center(2), cart_z - center(3)] ! get coordinate frame centered on IB + 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_x) .and. & + if (((.not. f_is_default(length(1)) .and. & xyz_local(2)**2 & + xyz_local(3)**2 <= radius**2 .and. & - -0.5*length(1) <= xyz_local(1) .and. & - 0.5*length(1) >= xyz_local(1)) & + -0.5_wp*length(1) <= xyz_local(1) .and. & + 0.5_wp*length(1) >= xyz_local(1)) & .or. & - (.not. f_is_default(length_y) .and. & + (.not. f_is_default(length(2)) .and. & xyz_local(1)**2 & + xyz_local(3)**2 <= radius**2 .and. & - -0.5*length(2) <= xyz_local(2) .and. & - 0.5*length(2) >= xyz_local(2)) & + -0.5_wp*length(2) <= xyz_local(2) .and. & + 0.5_wp*length(2) >= xyz_local(2)) & .or. & - (.not. f_is_default(length_z) .and. & + (.not. f_is_default(length(3)) .and. & xyz_local(1)**2 & + xyz_local(2)**2 <= radius**2 .and. & - -0.5*length(3) <= xyz_local(3) .and. & - 0.5*length(3) >= xyz_local(3)))) then + -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 From 91b256f2127d9e9e0b01a39592ac7f259b58725f Mon Sep 17 00:00:00 2001 From: Daniel J Vickers Date: Fri, 24 Oct 2025 13:39:08 -0400 Subject: [PATCH 13/20] Ran formatting --- src/common/m_ib_patches.fpp | 8 ++++---- src/post_process/m_data_input.f90 | 4 ++-- src/simulation/m_data_output.fpp | 14 +++++++------- src/simulation/m_ibm.fpp | 4 ++-- 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index 9e863137ca..52d36a17d4 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -139,7 +139,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 @@ -540,9 +540,9 @@ contains xy_local = matmul(inverse_rotation, 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(1) >= xy_local(1) .and. & -0.5_wp*length(2) <= xy_local(2) .and. & - 0.5_wp*length(2) >= xy_local(2)) then + 0.5_wp*length(2) >= xy_local(2)) then ! Updating the patch identities bookkeeping variable ib_markers_sf(i, j, 0) = patch_id @@ -569,7 +569,7 @@ contains ! Generic loop iterators integer :: i, j, k real(wp) :: radius - real(wp), dimension(1:3):: center + 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) diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index 4588387254..4429be273a 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -158,7 +158,7 @@ impure subroutine s_read_ib_data_files(file_loc_base, t_step) 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 @@ -179,7 +179,7 @@ impure subroutine s_read_ib_data_files(file_loc_base, t_step) 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 + save_index = t_step/t_step_save ! get the number of saves done to this point data_size = (m + 1)*(n + 1)*(p + 1) var_MOK = int(sys_size + 1, MPI_OFFSET_KIND) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index b578acfbba..532a5ebff5 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1022,16 +1022,16 @@ contains !Writeib 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) - + 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)) + 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) + '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) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 33042798f9..ffcb50fbed 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -271,9 +271,9 @@ contains 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 + 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 From 0389b67b308943ce3cece5061b1991bf01fcfb14 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Fri, 24 Oct 2025 16:27:00 -0400 Subject: [PATCH 14/20] The premature data save fixes MIBM frame 0 data, but ruins test suite. Need another solution --- src/simulation/p_main.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/p_main.fpp b/src/simulation/p_main.fpp index 3eefefd88f..6bc1ab5790 100644 --- a/src/simulation/p_main.fpp +++ b/src/simulation/p_main.fpp @@ -71,7 +71,7 @@ program p_main call nvtxEndRange ! INIT ! save the data at t_step=0 - call s_save_data(t_step, start, finish, io_time_avg, nt) + ! call s_save_data(t_step, start, finish, io_time_avg, nt) call nvtxStartRange("SIMULATION-TIME-MARCH") ! Time-stepping Loop From eae197453474767220f6d76d934d0708f0a6168a Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Mon, 27 Oct 2025 13:59:57 -0400 Subject: [PATCH 15/20] Added something that may fix the cray issue --- src/common/m_compute_levelset.fpp | 7 ++++--- src/simulation/m_data_output.fpp | 2 +- src/simulation/p_main.fpp | 3 --- 3 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index fa85808713..9aab9c2f86 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -264,7 +264,7 @@ contains integer, intent(in) :: ib_patch_id real(wp) :: top_right(2), bottom_left(2) - real(wp) :: min_dist + real(wp) :: min_dist, initial_min_dist real(wp) :: side_dists(4) real(wp) :: length_x, length_y @@ -286,9 +286,10 @@ contains top_right(2) = length_y/2 bottom_left(1) = -length_x/2 bottom_left(2) = -length_y/2 + initial_min_dist = initial_distance_buffer $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', copy='[levelset,levelset_norm]',& - & copyin='[ib_patch_id,center,bottom_left,top_right,initial_distance_buffer,inverse_rotation,rotation]', collapse=2) + & copyin='[ib_patch_id,center,bottom_left,top_right,initial_min_dist,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] @@ -301,7 +302,7 @@ contains 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 + min_dist = initial_min_dist idx = 1 do k = 1, 4 diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 532a5ebff5..a3a422abbb 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1020,7 +1020,7 @@ contains call MPI_FILE_CLOSE(ifile, ierr) - !Writeib data + !Write ib data if (ib) then write (file_loc, '(A)') 'ib.dat' file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) diff --git a/src/simulation/p_main.fpp b/src/simulation/p_main.fpp index 6bc1ab5790..29cb3b8281 100644 --- a/src/simulation/p_main.fpp +++ b/src/simulation/p_main.fpp @@ -70,9 +70,6 @@ program p_main call nvtxEndRange ! INIT - ! save the data at t_step=0 - ! call s_save_data(t_step, start, finish, io_time_avg, nt) - call nvtxStartRange("SIMULATION-TIME-MARCH") ! Time-stepping Loop do From 3dc61e0a49f0be8fea3db4cfe86c4a88b4943213 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Wed, 29 Oct 2025 10:42:15 -0400 Subject: [PATCH 16/20] Intermittend commit --- src/common/m_compute_levelset.fpp | 25 ++++++++++++------------- src/common/m_ib_patches.fpp | 22 +++++----------------- 2 files changed, 17 insertions(+), 30 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index 9aab9c2f86..ef58c01faf 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -43,7 +43,7 @@ contains center(1) = patch_ib(ib_patch_id)%x_centroid center(2) = patch_ib(ib_patch_id)%y_centroid - $:GPU_PARALLEL_LOOP(private='[i,j,dist_vec,dist]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,dist_vec,dist]', copy='[levelset%sf,levelset_norm%sf]',& & copyin='[ib_patch_id,center,radius]', collapse=2) do i = 0, m do j = 0, n @@ -86,7 +86,7 @@ contains inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) - $:GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,dist_vec,dist,global_dist,global_id]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,dist_vec,dist,global_dist,global_id]', copy='[levelset%sf,levelset_norm%sf]',& & copyin='[ib_patch_id,center,rotation,inverse_rotation,airfoil_grid_u,airfoil_grid_l]', collapse=2) do i = 0, m do j = 0, n @@ -178,7 +178,7 @@ contains z_max = center(3) + lz/2 z_min = center(3) - lz/2 - $:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,dist_vec,dist,global_dist,global_id,dist_side,dist_surf]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,dist_vec,dist,global_dist,global_id,dist_side,dist_surf]', copy='[levelset%sf,levelset_norm%sf]',& & 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 @@ -257,14 +257,14 @@ contains end subroutine s_3D_airfoil_levelset !> Initialize IBM module - pure subroutine s_rectangle_levelset(ib_patch_id, levelset, levelset_norm) + subroutine s_rectangle_levelset(ib_patch_id, levelset, levelset_norm) type(levelset_field), intent(INOUT), optional :: levelset type(levelset_norm_field), intent(INOUT), optional :: levelset_norm integer, intent(in) :: ib_patch_id real(wp) :: top_right(2), bottom_left(2) - real(wp) :: min_dist, initial_min_dist + real(wp) :: min_dist real(wp) :: side_dists(4) real(wp) :: length_x, length_y @@ -286,10 +286,9 @@ contains top_right(2) = length_y/2 bottom_left(1) = -length_x/2 bottom_left(2) = -length_y/2 - initial_min_dist = initial_distance_buffer - $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', copy='[levelset,levelset_norm]',& - & copyin='[ib_patch_id,center,bottom_left,top_right,initial_min_dist,inverse_rotation,rotation]', collapse=2) + $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', copy='[levelset%sf,levelset_norm%sf]',& + & copyin='[ib_patch_id,center,bottom_left,top_right,inverse_rotation,rotation,x_cc,y_cc]', collapse=2) do i = 0, m do j = 0, n xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] @@ -302,10 +301,10 @@ contains 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_min_dist + min_dist = side_dists(1) idx = 1 - do k = 1, 4 + do k = 2, 4 if (abs(side_dists(k)) < abs(min_dist)) then idx = k min_dist = side_dists(idx) @@ -368,7 +367,7 @@ contains Front = length_z/2 Back = -length_z/2 - $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local,dist_vec]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local,dist_vec]', copy='[levelset%sf,levelset_norm%sf]',& & copyin='[ib_patch_id,center,inverse_rotation,rotation,Right,Left,Top,Bottom,Front,Back]', collapse=3) do i = 0, m do j = 0, n @@ -457,7 +456,7 @@ contains center(2) = patch_ib(ib_patch_id)%y_centroid center(3) = patch_ib(ib_patch_id)%z_centroid - $:GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,dist]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,dist]', copy='[levelset%sf,levelset_norm%sf]',& & copyin='[ib_patch_id,center,radius]', collapse=3) do i = 0, m do j = 0, n @@ -521,7 +520,7 @@ contains dist_surface_vec = (/1, 1, 0/) end if - $:GPU_PARALLEL_LOOP(private='[i,j,k,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset,levelset_norm]',& + $:GPU_PARALLEL_LOOP(private='[i,j,k,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset%sf,levelset_norm%sf]',& & 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 diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index 52d36a17d4..9bdd16e828 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -503,24 +503,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 - - length(1) = length_x - length(2) = length_y - center(1) = x_centroid - center(2) = y_centroid - ! 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 @@ -532,7 +520,7 @@ contains ! 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, xy_local]', copy='[ib_markers_sf]',& - & copyin='[patch_id,center,length,inverse_rotation]', collapse=2) + & 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 From e1591e33b9d50c0a500bbef8958da056ab9ebab7 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Wed, 29 Oct 2025 11:35:01 -0400 Subject: [PATCH 17/20] Passes cray compiler tests --- src/common/m_compute_levelset.fpp | 7 +++++-- src/simulation/m_ibm.fpp | 7 +++++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index ef58c01faf..8dd13c3230 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -287,8 +287,10 @@ contains bottom_left(1) = -length_x/2 bottom_left(2) = -length_y/2 - $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', copy='[levelset%sf,levelset_norm%sf]',& - & copyin='[ib_patch_id,center,bottom_left,top_right,inverse_rotation,rotation,x_cc,y_cc]', collapse=2) + ! $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,idx,side_dists,xy_local,dist_vec]', copy='[levelset%sf,levelset_norm%sf, dbg]',& + ! & copyin='[ib_patch_id,center,bottom_left,top_right,inverse_rotation,rotation,x_cc,y_cc]', collapse=2) + $: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] @@ -330,6 +332,7 @@ contains end do end do + end subroutine s_rectangle_levelset pure subroutine s_cuboid_levelset(ib_patch_id, levelset, levelset_norm) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index ffcb50fbed..6a4aa4f58d 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -927,7 +927,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 @@ -941,7 +943,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) From a301fbc738db96c278e2ba3b04ac36a4e0d453d0 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Wed, 29 Oct 2025 11:44:23 -0400 Subject: [PATCH 18/20] Modified all GPU parallel statements and ran formatting --- src/common/m_compute_levelset.fpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index 8dd13c3230..a2c224dad5 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -43,7 +43,7 @@ contains center(1) = patch_ib(ib_patch_id)%x_centroid center(2) = patch_ib(ib_patch_id)%y_centroid - $:GPU_PARALLEL_LOOP(private='[i,j,dist_vec,dist]', copy='[levelset%sf,levelset_norm%sf]',& + $: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 @@ -86,7 +86,7 @@ contains inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) - $:GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,dist_vec,dist,global_dist,global_id]', copy='[levelset%sf,levelset_norm%sf]',& + $: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 @@ -178,7 +178,7 @@ contains z_max = center(3) + lz/2 z_min = center(3) - lz/2 - $:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,dist_vec,dist,global_dist,global_id,dist_side,dist_surf]', copy='[levelset%sf,levelset_norm%sf]',& + $: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 @@ -332,7 +332,6 @@ contains end do end do - end subroutine s_rectangle_levelset pure subroutine s_cuboid_levelset(ib_patch_id, levelset, levelset_norm) @@ -370,7 +369,7 @@ contains Front = length_z/2 Back = -length_z/2 - $:GPU_PARALLEL_LOOP(private='[i,j,k,min_dist,side_dists,xyz_local,dist_vec]', copy='[levelset%sf,levelset_norm%sf]',& + $: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 @@ -459,7 +458,7 @@ contains center(2) = patch_ib(ib_patch_id)%y_centroid center(3) = patch_ib(ib_patch_id)%z_centroid - $:GPU_PARALLEL_LOOP(private='[i,j,k,dist_vec,dist]', copy='[levelset%sf,levelset_norm%sf]',& + $: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 @@ -523,7 +522,7 @@ contains dist_surface_vec = (/1, 1, 0/) end if - $:GPU_PARALLEL_LOOP(private='[i,j,k,side_pos,dist_side,dist_surface,xyz_local]', copy='[levelset%sf,levelset_norm%sf]',& + $: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 From 7a0553987574f4aef6b60e67ae39f76a1057b08a Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 29 Oct 2025 14:15:09 -0400 Subject: [PATCH 19/20] Fixed IB marker write issue with something more sustainable --- src/post_process/m_data_input.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index 4429be273a..f82cca7466 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -183,7 +183,11 @@ impure subroutine s_read_ib_data_files(file_loc_base, t_step) data_size = (m + 1)*(n + 1)*(p + 1) 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(save_index, 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) From 76ec5eded1b6a8f6063f1039bff2e5c6af68e9e2 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 29 Oct 2025 20:06:48 -0400 Subject: [PATCH 20/20] Fixed bug with centroid and length not being comunicated over MPI --- src/simulation/m_mpi_proxy.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 81460b5c70..440a2eaad4 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']