From f85f9831b5b1c380cd8dec0fea05475b2910d548 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Tue, 9 Jun 2026 16:32:16 -0400 Subject: [PATCH 01/14] Initial modifications to allow us to use a unified boudns check on patch geometries --- examples/2D_mibm_particle_cloud/case.py | 1 + src/common/m_patch_geometries.fpp | 175 ++++++++++++++++++++++++ src/simulation/m_checker.fpp | 1 + src/simulation/m_global_parameters.fpp | 1 + src/simulation/m_ibm.fpp | 1 + toolchain/mfc/params/definitions.py | 2 + 6 files changed, 181 insertions(+) create mode 100644 src/common/m_patch_geometries.fpp diff --git a/examples/2D_mibm_particle_cloud/case.py b/examples/2D_mibm_particle_cloud/case.py index 68abf5badc..77cc908b1b 100644 --- a/examples/2D_mibm_particle_cloud/case.py +++ b/examples/2D_mibm_particle_cloud/case.py @@ -79,6 +79,7 @@ "ib": "T", "num_ibs": 0, "viscous": "T", + "many_ib_patch_parallelism": "T", # Collision model (soft-sphere, from 3D_mibm_sphere_head_on_collision) "collision_model": 1, "coefficient_of_restitution": 0.9, diff --git a/src/common/m_patch_geometries.fpp b/src/common/m_patch_geometries.fpp new file mode 100644 index 0000000000..658116a8ca --- /dev/null +++ b/src/common/m_patch_geometries.fpp @@ -0,0 +1,175 @@ +!> +!! @file +!! @brief Contains module m_ibm + +#:include 'macros.fpp' + +!> @brief Contains helper functions specific to various patch gemoetries for determining if a grid cell lies inside of or outside of a patch geometry +module m_patch_geometries + + use m_derived_types + use m_global_parameters + use m_variables_conversion + use m_helper + use m_helper_basic + use m_constants + use m_model + + implicit none + + public :: f_is_inside_circle, f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_rectangle, f_is_inside_cuboid, & + & f_is_inside_airfoil, f_is_inside_airfoil_3D, f_is_inside_ellipse + +contains + + !> Check if the x, and y coordinates would be located inside a circle with the patch_id's radius + function f_is_inside_circle(patch_id, x, y) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y + logical :: is_inside + + is_inside = x**2 + y**2 <= patch_ib(patch_id)%radius**2 + + end function f_is_inside_circle + + !> Check if the x, y, and z coordinates would be located inside a sphere with the patch_id's radius + function f_is_inside_sphere(patch_id, x, y, z) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside + + is_inside = x**2 + y**2 + z**2 <= patch_ib(patch_id)%radius**2 + + end function f_is_inside_sphere + + !> Check which length of the cylinder is not default. Use that direction as the height and the other two coordinate + ! values as the radius check + function f_is_inside_cylinder(patch_id, x, y, z) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside + + ! check if the cylinder is extended in the x direction + is_inside = (.not. f_is_default(patch_ib(patch_id)%length_x)) .and. y**2 + z & + & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_x <= x & + & .and. 0.5_wp*patch_ib(patch_id)%length_x >= x + + ! check if the cylinder is extended in the z direction + is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_y)) .and. x**2 + z & + & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_y <= y & + & .and. 0.5_wp*patch_ib(patch_id)%length_y >= y) + + ! check if the cylinder is extended in the z direction + is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_z)) .and. x**2 + y & + & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_z <= z & + & .and. 0.5_wp*patch_ib(patch_id)%length_z >= z) + + ! TODO :: This can easily reuse the f_is_inside_circle function with an additional length requirement + + end function f_is_inside_cylinder + + !> Check if the x, y, and possibly z coordinates would be located inside a cuboid with the patch_id's lengths + function f_is_inside_cuboid(patch_id, x, y, z) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside + + ! check if x and y are inside the rectangle plane at z=0 + is_inside = -0.5_wp*patch_ib(patch_id)%length_x <= x .and. 0.5_wp*patch_ib(patch_id)%length_x >= x & + is_inside = is_inside .and. -0.5_wp*patch_ib(patch_id)%length_y <= y.and. 0.5_wp*patch_ib(patch_id)%length_y >= y + + ! if we are in 3D, this is a cuboid and so we must also check the z axis + if (num_dims == 3) is_inside = is_inside .and. -0.5_wp*patch_ib(patch_id)%length_z <= z .and. 0.5_wp*patch_ib(patch_id)%length_z >= z + + end function f_is_inside_cuboid + + !> Check if the x, y, are bounded by a NACA airfoil. Check if the z coordinate is inside the left and right edges of the airfoil, if set. + function f_is_inside_airfoil(patch_id, x, y, z) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside + + integer :: k, airfoil_id + real(wp) :: f + + airfoil_id = patch_ib(patch_id)%airfoil_id + + is_inside = .false. + + ! check the initial x bounds of the grid cell + if (.not. (x >= 0._wp .and. x <= ib_airfoil(airfoil_id)%c)) return + + ! if we are in 3D, we must also check the z axis + if (num_dims == 3 .and. (.not. (-0.5_wp*patch_ib(patch_id)%length_z <= z .and. 0.5_wp*patch_ib(patch_id)%length_z >= z))) return + + ! our check branches for the upper and lower half of the airfoil + if (y >= 0._wp) then + + ! incriment the iterator so we know where in the airfoil arrays to look + k = 1 + do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < x) + k = k + 1 + end do + + ! If the values are approximately equivilent, skip the next check + if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, x)) then + if (y <= ib_airfoil_grids(airfoil_id)%upper(k)%y) is_inside = .true. + else + + ! check if the y value is below the upper edge of the airfoil + f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - x) & + & /(ib_airfoil_grids(airfoil_id)%upper(k)%x - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x) + if (y <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y & + & + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) is_inside = .true. + + end if + else + + ! incriment the iterator so we know where in the airfoil arrays to look + k = 1 + do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < x) + k = k + 1 + end do + + ! If the values are approximately equivilent, skip the next check + if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, x)) then + if (y >= ib_airfoil_grids(airfoil_id)%lower(k)%y) is_inside = .true. + else + + ! check if the y value is above the lower edge of the airfoil + f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - x) & + & /(ib_airfoil_grids(airfoil_id)%lower(k)%x - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x) + if (y >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y & + & + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) is_inside = .true. + + end if + end if + + end function f_is_inside_airfoil + + function f_is_inside_ellipse(patch_id, x, y) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y + logical :: is_inside + + end function f_is_inside_ellipse + + end module m_patch_geometries diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index 2ed8aca3c3..d1586ec337 100644 --- a/src/simulation/m_checker.fpp +++ b/src/simulation/m_checker.fpp @@ -37,6 +37,7 @@ contains call s_check_inputs_time_stepping @:PROHIBIT(ib_state_wrt .and. .not. ib, "ib_state_wrt requires ib to be enabled") + @:PROHIBIT(many_ib_patch_parallelism .and. .not. ib, "many_ib_patch_parallelism requires ib to be enabled") if (num_particle_clouds > 0) then call s_check_inputs_particle_clouds diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 1bd6ea26f8..3f4eaea555 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -534,6 +534,7 @@ contains collision_time = dflt_real ib_coefficient_of_friction = dflt_real ib_state_wrt = .false. + many_ib_patch_parallelism = .false. ! Bubble modeling bubbles_euler = .false. diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 253123184d..f95a332c3d 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -19,6 +19,7 @@ module m_ibm use m_ib_patches use m_viscous use m_model + use m_patch_geometries use m_collisions implicit none diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index e22787089b..10eb666d15 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -617,6 +617,7 @@ def _load(): _r("coefficient_of_restitution", REAL, {"ib"}) _r("collision_time", REAL, {"ib"}) _r("ib_coefficient_of_friction", REAL, {"ib"}) + _r("many_ib_patch_parallelism", LOG, {"ib"}) # Probes for n in ["num_probes", "num_integrals"]: @@ -1254,6 +1255,7 @@ def _nv(targets: set, *names: str) -> None: "ib_coefficient_of_friction", "num_particle_clouds", "ib_neighborhood_radius", + "many_ib_patch_parallelism", "particle_cloud", "tau_star", "cont_damage_s", From 4c944543bc610ec4e3d2d5c19de34aeae7016759 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Tue, 9 Jun 2026 21:44:17 -0400 Subject: [PATCH 02/14] Removed all of the old outdated subroutines and merged them into a monolithic subroutine for grid-cell parallelism --- src/common/m_patch_geometries.fpp | 87 +++- src/simulation/m_ib_patches.fpp | 703 +++++------------------------- 2 files changed, 191 insertions(+), 599 deletions(-) diff --git a/src/common/m_patch_geometries.fpp b/src/common/m_patch_geometries.fpp index 658116a8ca..7cbf56bc0c 100644 --- a/src/common/m_patch_geometries.fpp +++ b/src/common/m_patch_geometries.fpp @@ -17,8 +17,8 @@ module m_patch_geometries implicit none - public :: f_is_inside_circle, f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_rectangle, f_is_inside_cuboid, & - & f_is_inside_airfoil, f_is_inside_airfoil_3D, f_is_inside_ellipse + public :: f_is_inside_circle, f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_cuboid, & + & f_is_inside_airfoil, f_is_inside_ellipse contains @@ -172,4 +172,87 @@ contains end function f_is_inside_ellipse + subroutine s_get_bounding_box_corner_distance(patch_id, bounding_distance) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(out) :: bounding_distance + + ! TODO :: FILL ME IN + + end subroutine s_get_bounding_box_corner_distance + + !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid) + !! and before s_apply_ib_patches or s_apply_levelset. + subroutine s_initialize_ib_airfoils() + + integer :: i, j, airfoil_id + integer :: Np, Np1, Np2 + real(wp) :: ca_in, pa, ma, ta + real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c + + do i = 1, num_ibs + if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle + + airfoil_id = patch_ib(i)%airfoil_id + ca_in = ib_airfoil(airfoil_id)%c + pa = ib_airfoil(airfoil_id)%p + ma = ib_airfoil(airfoil_id)%m + ta = ib_airfoil(airfoil_id)%t + + Np1 = int((pa*ca_in/dx(0))*20) + Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) + Np = Np1 + Np2 + 1 + ib_airfoil_grids(airfoil_id)%Np = Np + $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%Np]') + + if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then + @:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np)) + @:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np)) + + ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp + ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp + ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp + ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp + + do j = 1, Np1 + Np2 - 1 + if (j <= Np1) then + xc = j*(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 + (j - 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 + + xu = (xa - yt*sin_c)*ca_in + yu = (yc + yt*cos_c)*ca_in + xl = (xa + yt*sin_c)*ca_in + yl = (yc - yt*cos_c)*ca_in + + ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu + ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu + ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl + ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl + end do + + ib_airfoil_grids(airfoil_id)%upper(Np)%x = ca_in + ib_airfoil_grids(airfoil_id)%upper(Np)%y = 0._wp + ib_airfoil_grids(airfoil_id)%lower(Np)%x = ca_in + ib_airfoil_grids(airfoil_id)%lower(Np)%y = 0._wp + + $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]') + end if + end do + + end subroutine s_initialize_ib_airfoils + end module m_patch_geometries diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 41f8524b87..18ce5d188f 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -26,106 +26,81 @@ module m_ib_patches contains - !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid) - !! and before s_apply_ib_patches or s_apply_levelset. - subroutine s_initialize_ib_airfoils() - - integer :: i, j, airfoil_id - integer :: Np, Np1, Np2 - real(wp) :: ca_in, pa, ma, ta - real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c - - do i = 1, num_ibs - if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle - - airfoil_id = patch_ib(i)%airfoil_id - ca_in = ib_airfoil(airfoil_id)%c - pa = ib_airfoil(airfoil_id)%p - ma = ib_airfoil(airfoil_id)%m - ta = ib_airfoil(airfoil_id)%t - - Np1 = int((pa*ca_in/dx(0))*20) - Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) - Np = Np1 + Np2 + 1 - ib_airfoil_grids(airfoil_id)%Np = Np - $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%Np]') - - if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then - @:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np)) - @:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np)) - - ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp - ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp - ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp - ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp - - do j = 1, Np1 + Np2 - 1 - if (j <= Np1) then - xc = j*(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 + (j - 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 - - xu = (xa - yt*sin_c)*ca_in - yu = (yc + yt*cos_c)*ca_in - xl = (xa + yt*sin_c)*ca_in - yl = (yc - yt*cos_c)*ca_in - - ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu - ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu - ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl - ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl - end do - - ib_airfoil_grids(airfoil_id)%upper(Np)%x = ca_in - ib_airfoil_grids(airfoil_id)%upper(Np)%y = 0._wp - ib_airfoil_grids(airfoil_id)%lower(Np)%x = ca_in - ib_airfoil_grids(airfoil_id)%lower(Np)%y = 0._wp - - $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]') - end if - end do - - end subroutine s_initialize_ib_airfoils - !> Apply all immersed boundary patch geometries to mark interior cells in the IB marker array impure subroutine s_apply_ib_patches(ib_markers) type(integer_field), intent(inout) :: ib_markers - integer :: i, xp, yp, zp !< iterators + + if (many_ib_patch_parallelism) then + call s_apply_ib_patches_ib_parallelism(ib_markers) + else + call s_apply_ib_patches_grid_cell_parallelism(ib_markers) + end if + + end subroutine s_apply_ib_patches + + subroutine s_apply_ib_patches_grid_cell_parallelism(ib_markers) + + type(integer_field), intent(inout) :: ib_markers + integer :: patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds + real(wp), dimension(3) :: center, xyz_local + real(wp) :: bounding_box_corner_distance ! 3D Patch Geometries - - if (p > 0) then + if (num_dims == 3) then !> IB Patches !> @{ call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper) do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper do zp = zp_lower, zp_upper - do i = 1, num_ibs - if (patch_ib(i)%geometry == 8) then - call s_ib_sphere(i, ib_markers, xp, yp, zp) - else if (patch_ib(i)%geometry == 9) then - call s_ib_cuboid(i, ib_markers, xp, yp, zp) - else if (patch_ib(i)%geometry == 10) then - call s_ib_cylinder(i, ib_markers, xp, yp, zp) - else if (patch_ib(i)%geometry == 11) then - call s_ib_3D_airfoil(i, ib_markers, xp, yp, zp) - else if (patch_ib(i)%geometry == 12) then - call s_ib_3d_model(i, ib_markers, xp, yp, zp) - end if + do patch_id = 1, num_ibs + 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) + call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) + + ! 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) + + ! find the indices to the left and right of the IB in i, j, k + il = -gp_layers - 1 + jl = -gp_layers - 1 + kl = -gp_layers - 1 + ir = m + gp_layers + 1 + jr = n + gp_layers + 1 + kr = p + gp_layers + 1 + call get_bounding_indices(center(1) - bounding_box_corner_distance, center(1) + bounding_box_corner_distance, x_cc, il, ir) + call get_bounding_indices(center(2) - bounding_box_corner_distance, center(2) + bounding_box_corner_distance, y_cc, jl, jr) + call get_bounding_indices(center(3) - bounding_box_corner_distance, center(3) + bounding_box_corner_distance, z_cc, kl, kr) + + $:GPU_PARALLEL_LOOP(private='[i, j, k]', copyin='[patch_id, encoded_patch_id, center]', collapse=3) + do k = kl, kr + do j = jl, jr + do i = il, ir + ! get coordinate frame centered on IB + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] + ! rotate the frame into the IB's coordinates + xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) + xyz_local = xyz_local - patch_ib(patch_id)%offset ! airfoils are a patch that require a centroid offset + + ! perform the interior check for the patch geometry of this IB + if (patch_ib(patch_id)%geometry == 8) then + if (f_is_inside_sphere(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 9) then + if (f_is_inside_cuboid(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(ipatch_id)%geometry == 10) then + if (f_is_inside_cylinder(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 11) then + if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 12) then + if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + end if + end do + end do + end do + $:END_GPU_PARALLEL_LOOP() end do end do end do @@ -139,524 +114,58 @@ contains call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper) do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper - do i = 1, num_ibs - if (patch_ib(i)%geometry == 2) then - call s_ib_circle(i, ib_markers, xp, yp) - else if (patch_ib(i)%geometry == 3) then - call s_ib_rectangle(i, ib_markers, xp, yp) - else if (patch_ib(i)%geometry == 4) then - call s_ib_airfoil(i, ib_markers, xp, yp) - else if (patch_ib(i)%geometry == 5) then - call s_ib_model(i, ib_markers, xp, yp) - else if (patch_ib(i)%geometry == 6) then - call s_ib_ellipse(i, ib_markers, xp, yp) - end if + do patch_id = 1, num_ibs + 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) = 0._wp + call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) + + ! 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 + il = -gp_layers - 1 + jl = -gp_layers - 1 + ir = m + gp_layers + 1 + jr = n + gp_layers + 1 + call get_bounding_indices(center(1) - bounding_box_corner_distance, center(1) + bounding_box_corner_distance, x_cc, il, ir) + call get_bounding_indices(center(2) - bounding_box_corner_distance, center(2) + bounding_box_corner_distance, y_cc, jl, jr) + + $:GPU_PARALLEL_LOOP(private='[i, j, k]', copyin='[patch_id, encoded_patch_id, center]', collapse=2) + do j = jl, jr + do i = il, ir + ! get coordinate frame centered on IB + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] + ! rotate the frame into the IB's coordinates + xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) + xyz_local = xyz_local - patch_ib(patch_id)%offset ! airfoils are a patch that require a centroid offset + + ! perform the interior check for the patch geometry of this IB + if (patch_ib(patch_id)%geometry == 2) then + if (f_is_inside_circle(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 3) then + if (f_is_inside_rectangle(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(ipatch_id)%geometry == 4) then + if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 5) then + if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 6) then + if (f_is_inside_ellipse(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + end if + end do + end do + $:END_GPU_PARALLEL_LOOP() end do end do end do !> @} end if - end subroutine s_apply_ib_patches - - !> Mark cells inside a circular immersed boundary - subroutine s_ib_circle(patch_id, ib_markers, xp, yp) - - integer, intent(in) :: patch_id - integer, intent(in) :: xp, yp !< integers containing the periodicity projection information - type(integer_field), intent(inout) :: ib_markers - real(wp), dimension(1:2) :: center - real(wp) :: radius - integer :: i, j, il, ir, jl, jr !< Generic loop iterators - integer :: encoded_patch_id - - ! Transferring the circular patch's radius, centroid, smearing patch identity and smearing coefficient information + end subroutine s_apply_ib_patches_grid_cell_parallelism - 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) - radius = patch_ib(patch_id)%radius - - ! 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 - il = -gp_layers - 1 - jl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir) - call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr) - - ! Assign primitive variables if circle covers cell and patch has write permission - - $:GPU_PARALLEL_LOOP(private='[i, j]', copyin='[encoded_patch_id, center, radius]', collapse=2) - do j = jl, jr - do i = il, ir - if ((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 <= radius**2) then - ib_markers%sf(i, j, 0) = encoded_patch_id - end if - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_circle - - !> Mark cells inside a 2D NACA 4-digit airfoil immersed boundary - subroutine s_ib_airfoil(patch_id, ib_markers, xp, yp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - integer, intent(in) :: xp, yp !< integers containing the periodicity projection information - real(wp) :: f, ca_in - integer :: i, j, k, il, ir, jl, jr - integer :: Np_local, airfoil_id - integer :: encoded_patch_id - real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame - real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame - - airfoil_id = patch_ib(patch_id)%airfoil_id - center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) - center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) - ca_in = ib_airfoil(airfoil_id)%c - Np_local = ib_airfoil_grids(airfoil_id)%Np - offset(:) = patch_ib(patch_id)%centroid_offset(:) - - ! 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 - il = -gp_layers - 1 - jl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - ! maximum distance any marker can be from the center is the length of the airfoil - call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir) - call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr) - - $:GPU_PARALLEL_LOOP(private='[i, j, xy_local, k, f]', copyin='[encoded_patch_id, center, offset, ca_in, airfoil_id, & - & Np_local, ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]', collapse=2) - do j = jl, jr - do i = il, ir - xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB - xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xy_local) ! rotate the frame into the IB's coordinates - xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset - - if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then - if (xy_local(2) >= 0._wp) then - k = 1 - do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < xy_local(1) .and. k <= Np_local) - k = k + 1 - end do - if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, xy_local(1))) then - if (xy_local(2) <= ib_airfoil_grids(airfoil_id)%upper(k)%y) then - ib_markers%sf(i, j, 0) = encoded_patch_id - end if - else - f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - xy_local(1))/(ib_airfoil_grids(airfoil_id)%upper(k)%x & - & - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x) - if (xy_local(2) <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y & - & + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) then - ib_markers%sf(i, j, 0) = encoded_patch_id - end if - end if - else - k = 1 - do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < xy_local(1)) - k = k + 1 - end do - if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, xy_local(1))) then - if (xy_local(2) >= ib_airfoil_grids(airfoil_id)%lower(k)%y) then - ib_markers%sf(i, j, 0) = encoded_patch_id - end if - else - f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - xy_local(1))/(ib_airfoil_grids(airfoil_id)%lower(k)%x & - & - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x) - if (xy_local(2) >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y & - & + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) then - ib_markers%sf(i, j, 0) = encoded_patch_id - end if - end if - end if - end if - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_airfoil - - !> Mark cells inside a 3D extruded NACA 4-digit airfoil immersed boundary with finite span - subroutine s_ib_3D_airfoil(patch_id, ib_markers, xp, yp, zp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information - real(wp) :: lz, z_max, z_min, f, ca_in - integer :: i, j, k, l, il, ir, jl, jr, ll, lr - integer :: airfoil_id - integer :: encoded_patch_id - real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame - - airfoil_id = patch_ib(patch_id)%airfoil_id - center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) - center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) - center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg) - lz = patch_ib(patch_id)%length_z - ca_in = ib_airfoil(airfoil_id)%c - offset(:) = patch_ib(patch_id)%centroid_offset(:) - - z_max = lz/2 - z_min = -lz/2 - - ! 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) - - ! find the indices to the left and right of the IB in i, j, k - il = -gp_layers - 1 - jl = -gp_layers - 1 - ll = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - lr = p + gp_layers + 1 - ! maximum distance any marker can be from the center is the length of the airfoil - call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir) - call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr) - call get_bounding_indices(center(3) - ca_in, center(3) + ca_in, z_cc, ll, lr) - - $:GPU_PARALLEL_LOOP(private='[i, j, l, xyz_local, k, f]', copyin='[encoded_patch_id, center, offset, ca_in, airfoil_id, & - & ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower, z_min, z_max]', collapse=3) - do l = ll, lr - do j = jl, jr - do i = il, ir - ! get coordinate frame centered on IB - xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] - ! rotate the frame into the IB's coordinates - xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset - - if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then - if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then - if (xyz_local(2) >= 0._wp) then - k = 1 - do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < xyz_local(1)) - k = k + 1 - end do - if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, xyz_local(1))) then - if (xyz_local(2) <= ib_airfoil_grids(airfoil_id)%upper(k)%y) then - ! IB - ib_markers%sf(i, j, l) = encoded_patch_id - end if - else - f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - xyz_local(1)) & - & /(ib_airfoil_grids(airfoil_id)%upper(k)%x - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x) - if (xyz_local(2) <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y & - & + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) then - ib_markers%sf(i, j, l) = encoded_patch_id - end if - end if - else - k = 1 - do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < xyz_local(1)) - k = k + 1 - end do - if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, xyz_local(1))) then - if (xyz_local(2) >= ib_airfoil_grids(airfoil_id)%lower(k)%y) then - ib_markers%sf(i, j, l) = encoded_patch_id - end if - else - f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - xyz_local(1)) & - & /(ib_airfoil_grids(airfoil_id)%lower(k)%x - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x) - if (xyz_local(2) >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y & - & + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) then - ib_markers%sf(i, j, l) = encoded_patch_id - end if - end if - end if - end if - end if - end do - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_3D_airfoil - - !> Mark cells inside a rectangular immersed boundary - subroutine s_ib_rectangle(patch_id, ib_markers, xp, yp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - integer, intent(in) :: xp, yp !< integers containing the periodicity projection information - integer :: i, j, il, ir, jl, jr !< generic loop iterators - integer :: encoded_patch_id - real(wp) :: corner_distance !< Equation of state parameters - real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame - real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame - - ! Transferring the rectangle's centroid and length information - - center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg) - center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg) - - ! 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 - il = -gp_layers - 1 - jl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - ! 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]', 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(patch_ib(patch_id)%rotation_matrix_inverse, xy_local) - - 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 - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_rectangle - - !> Mark cells inside a spherical immersed boundary - subroutine s_ib_sphere(patch_id, ib_markers, xp, yp, zp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information - ! Generic loop iterators - integer :: i, j, k - integer :: il, ir, jl, jr, kl, kr - integer :: encoded_patch_id - real(wp) :: radius - real(wp), dimension(1:3) :: center - - 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) - radius = patch_ib(patch_id)%radius - - ! completely skip particles not in the domain - if (center(1) - radius > x_cc(m + gp_layers + 1) .or. center(1) + radius < x_cc(-gp_layers - 1) .or. center(2) & - & - radius > y_cc(n + gp_layers + 1) .or. center(2) + radius < y_cc(-gp_layers - 1) .or. center(3) - radius > z_cc(p & - & + gp_layers + 1) .or. center(3) + radius < z_cc(-gp_layers - 1)) then - return - end if - - ! 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) - - ! find the indices to the left and right of the IB in i, j, k - il = -gp_layers - 1 - jl = -gp_layers - 1 - kl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - kr = p + gp_layers + 1 - call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir) - call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr) - call get_bounding_indices(center(3) - radius, center(3) + radius, z_cc, kl, kr) - - ! Checking whether the sphere covers a particular cell in the domain 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]', copyin='[encoded_patch_id, center, radius]', collapse=3) - do k = kl, kr - do j = jl, jr - do i = il, ir - ! Updating the patch identities bookkeeping variable - if (((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 + (z_cc(k) - center(3))**2 <= radius**2)) then - ib_markers%sf(i, j, k) = encoded_patch_id - end if - end do - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_sphere - - !> Mark cells inside a cuboidal immersed boundary - subroutine s_ib_cuboid(patch_id, ib_markers, xp, yp, zp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - 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 !< x and y coordinates in local IB frame - real(wp) :: corner_distance - - ! 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) - - ! 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) - - ! find the indices to the left and right of the IB in i, j, k - il = -gp_layers - 1 - jl = -gp_layers - 1 - kl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - kr = p + gp_layers + 1 - 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) - - ! 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]', 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 - ! 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 - end do - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_cuboid - - !> Mark cells inside a cylindrical immersed boundary - subroutine s_ib_cylinder(patch_id, ib_markers, xp, yp, zp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - 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), 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 - - 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 = [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) - - il = -gp_layers - 1 - jl = -gp_layers - 1 - kl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - kr = p + gp_layers + 1 - 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) - - ! 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]', 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 - ! 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 - end do - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_cylinder - - !> Mark cells inside a 2D elliptical immersed boundary - subroutine s_ib_ellipse(patch_id, ib_markers, xp, yp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - integer, intent(in) :: xp, yp !< integers containing the periodicity projection information - 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) :: center !< x and y coordinates in local IB frame - real(wp) :: bounding_radius - - ! 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) - - ! 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) - 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]', 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(patch_ib(patch_id)%rotation_matrix_inverse(:,:), xy_local) - - ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1 - 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 - end do - end do - $:END_GPU_PARALLEL_LOOP() + subroutine s_apply_ib_patches_ib_parallelism() - end subroutine s_ib_ellipse + end subroutine s_apply_ib_patches_ib_parallelism !> The STL patch is a 2D geometry that is imported from an STL file. subroutine s_ib_model(patch_id, ib_markers, xp, yp) From 6f763746a0b8ff2d0b452db449680eb2fe0b5595 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Tue, 9 Jun 2026 21:51:04 -0400 Subject: [PATCH 03/14] part of the way through adding a secondary loop --- src/common/m_patch_geometries.fpp | 354 +++++++++++++++--------------- src/simulation/m_ib_patches.fpp | 138 +++++++++--- 2 files changed, 286 insertions(+), 206 deletions(-) diff --git a/src/common/m_patch_geometries.fpp b/src/common/m_patch_geometries.fpp index 7cbf56bc0c..3734eb8b55 100644 --- a/src/common/m_patch_geometries.fpp +++ b/src/common/m_patch_geometries.fpp @@ -4,7 +4,8 @@ #:include 'macros.fpp' -!> @brief Contains helper functions specific to various patch gemoetries for determining if a grid cell lies inside of or outside of a patch geometry +!> @brief Contains helper functions specific to various patch gemoetries for determining if a grid cell lies inside of or outside of +!! a patch geometry module m_patch_geometries use m_derived_types @@ -17,34 +18,34 @@ module m_patch_geometries implicit none - public :: f_is_inside_circle, f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_cuboid, & - & f_is_inside_airfoil, f_is_inside_ellipse + public :: f_is_inside_circle, f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_cuboid, f_is_inside_airfoil, & + & f_is_inside_ellipse contains !> Check if the x, and y coordinates would be located inside a circle with the patch_id's radius function f_is_inside_circle(patch_id, x, y) result(is_inside) - $:GPU_ROUTINE(parallelism='[seq]') + $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y - logical :: is_inside + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y + logical :: is_inside - is_inside = x**2 + y**2 <= patch_ib(patch_id)%radius**2 + is_inside = x**2 + y**2 <= patch_ib(patch_id)%radius**2 end function f_is_inside_circle !> Check if the x, y, and z coordinates would be located inside a sphere with the patch_id's radius function f_is_inside_sphere(patch_id, x, y, z) result(is_inside) - $:GPU_ROUTINE(parallelism='[seq]') + $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z - logical :: is_inside + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside - is_inside = x**2 + y**2 + z**2 <= patch_ib(patch_id)%radius**2 + is_inside = x**2 + y**2 + z**2 <= patch_ib(patch_id)%radius**2 end function f_is_inside_sphere @@ -52,207 +53,202 @@ contains ! values as the radius check function f_is_inside_cylinder(patch_id, x, y, z) result(is_inside) - $:GPU_ROUTINE(parallelism='[seq]') + $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z - logical :: is_inside + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside - ! check if the cylinder is extended in the x direction - is_inside = (.not. f_is_default(patch_ib(patch_id)%length_x)) .and. y**2 + z & - & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_x <= x & - & .and. 0.5_wp*patch_ib(patch_id)%length_x >= x + ! check if the cylinder is extended in the x direction + is_inside = (.not. f_is_default(patch_ib(patch_id)%length_x)) .and. y**2 + z**2 <= patch_ib(patch_id)%radius**2 .and. & + & -0.5_wp*patch_ib(patch_id)%length_x <= x .and. 0.5_wp*patch_ib(patch_id)%length_x >= x - ! check if the cylinder is extended in the z direction - is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_y)) .and. x**2 + z & - & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_y <= y & - & .and. 0.5_wp*patch_ib(patch_id)%length_y >= y) + ! check if the cylinder is extended in the z direction + is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_y)) .and. x**2 & + & + z**2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_y <= y & + & .and. 0.5_wp*patch_ib(patch_id)%length_y >= y) - ! check if the cylinder is extended in the z direction - is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_z)) .and. x**2 + y & - & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_z <= z & - & .and. 0.5_wp*patch_ib(patch_id)%length_z >= z) + ! check if the cylinder is extended in the z direction + is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_z)) .and. x**2 & + & + y**2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_z <= z & + & .and. 0.5_wp*patch_ib(patch_id)%length_z >= z) - ! TODO :: This can easily reuse the f_is_inside_circle function with an additional length requirement + ! TODO :: This can easily reuse the f_is_inside_circle function with an additional length requirement end function f_is_inside_cylinder !> Check if the x, y, and possibly z coordinates would be located inside a cuboid with the patch_id's lengths function f_is_inside_cuboid(patch_id, x, y, z) result(is_inside) - $:GPU_ROUTINE(parallelism='[seq]') + $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z - logical :: is_inside + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside - ! check if x and y are inside the rectangle plane at z=0 - is_inside = -0.5_wp*patch_ib(patch_id)%length_x <= x .and. 0.5_wp*patch_ib(patch_id)%length_x >= x & - is_inside = is_inside .and. -0.5_wp*patch_ib(patch_id)%length_y <= y.and. 0.5_wp*patch_ib(patch_id)%length_y >= y + ! check if x and y are inside the rectangle plane at z=0 + is_inside = -0.5_wp*patch_ib(patch_id)%length_x <= x .and. 0.5_wp*patch_ib(patch_id)%length_x >= x is_inside = is_inside & + & .and. -0.5_wp*patch_ib(patch_id)%length_y <= y .and. 0.5_wp*patch_ib(patch_id)%length_y >= y - ! if we are in 3D, this is a cuboid and so we must also check the z axis - if (num_dims == 3) is_inside = is_inside .and. -0.5_wp*patch_ib(patch_id)%length_z <= z .and. 0.5_wp*patch_ib(patch_id)%length_z >= z + ! if we are in 3D, this is a cuboid and so we must also check the z axis + if (num_dims == 3) is_inside = is_inside .and. -0.5_wp*patch_ib(patch_id)%length_z <= z & + & .and. 0.5_wp*patch_ib(patch_id)%length_z >= z end function f_is_inside_cuboid - !> Check if the x, y, are bounded by a NACA airfoil. Check if the z coordinate is inside the left and right edges of the airfoil, if set. + !> Check if the x, y, are bounded by a NACA airfoil. Check if the z coordinate is inside the left and right edges of the + !! airfoil, if set. function f_is_inside_airfoil(patch_id, x, y, z) result(is_inside) - $:GPU_ROUTINE(parallelism='[seq]') - - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z - logical :: is_inside - - integer :: k, airfoil_id - real(wp) :: f - - airfoil_id = patch_ib(patch_id)%airfoil_id - - is_inside = .false. - - ! check the initial x bounds of the grid cell - if (.not. (x >= 0._wp .and. x <= ib_airfoil(airfoil_id)%c)) return - - ! if we are in 3D, we must also check the z axis - if (num_dims == 3 .and. (.not. (-0.5_wp*patch_ib(patch_id)%length_z <= z .and. 0.5_wp*patch_ib(patch_id)%length_z >= z))) return - - ! our check branches for the upper and lower half of the airfoil - if (y >= 0._wp) then - - ! incriment the iterator so we know where in the airfoil arrays to look - k = 1 - do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < x) - k = k + 1 - end do - - ! If the values are approximately equivilent, skip the next check - if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, x)) then - if (y <= ib_airfoil_grids(airfoil_id)%upper(k)%y) is_inside = .true. - else - - ! check if the y value is below the upper edge of the airfoil - f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - x) & - & /(ib_airfoil_grids(airfoil_id)%upper(k)%x - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x) - if (y <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y & - & + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) is_inside = .true. - - end if - else - - ! incriment the iterator so we know where in the airfoil arrays to look - k = 1 - do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < x) - k = k + 1 - end do - - ! If the values are approximately equivilent, skip the next check - if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, x)) then - if (y >= ib_airfoil_grids(airfoil_id)%lower(k)%y) is_inside = .true. - else - - ! check if the y value is above the lower edge of the airfoil - f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - x) & - & /(ib_airfoil_grids(airfoil_id)%lower(k)%x - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x) - if (y >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y & - & + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) is_inside = .true. - - end if - end if + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y, z + logical :: is_inside + integer :: k, airfoil_id + real(wp) :: f + + airfoil_id = patch_ib(patch_id)%airfoil_id + + is_inside = .false. + + ! check the initial x bounds of the grid cell + if (.not. (x >= 0._wp .and. x <= ib_airfoil(airfoil_id)%c)) return + + ! if we are in 3D, we must also check the z axis + if (num_dims == 3 .and. (.not. (-0.5_wp*patch_ib(patch_id)%length_z <= z .and. 0.5_wp*patch_ib(patch_id)%length_z >= z))) & + & return + + ! our check branches for the upper and lower half of the airfoil + if (y >= 0._wp) then + ! incriment the iterator so we know where in the airfoil arrays to look + k = 1 + do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < x) + k = k + 1 + end do + + ! If the values are approximately equivilent, skip the next check + if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, x)) then + if (y <= ib_airfoil_grids(airfoil_id)%upper(k)%y) is_inside = .true. + else + ! check if the y value is below the upper edge of the airfoil + f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - x)/(ib_airfoil_grids(airfoil_id)%upper(k)%x & + & - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x) + if (y <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) & + & is_inside = .true. + end if + else + ! incriment the iterator so we know where in the airfoil arrays to look + k = 1 + do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < x) + k = k + 1 + end do + + ! If the values are approximately equivilent, skip the next check + if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, x)) then + if (y >= ib_airfoil_grids(airfoil_id)%lower(k)%y) is_inside = .true. + else + ! check if the y value is above the lower edge of the airfoil + f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - x)/(ib_airfoil_grids(airfoil_id)%lower(k)%x & + & - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x) + if (y >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) & + & is_inside = .true. + end if + end if end function f_is_inside_airfoil function f_is_inside_ellipse(patch_id, x, y) result(is_inside) - $:GPU_ROUTINE(parallelism='[seq]') + $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y - logical :: is_inside + integer, intent(in) :: patch_id + real(wp), intent(in) :: x, y + logical :: is_inside end function f_is_inside_ellipse subroutine s_get_bounding_box_corner_distance(patch_id, bounding_distance) - $:GPU_ROUTINE(parallelism='[seq]') + $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(out) :: bounding_distance + integer, intent(in) :: patch_id + real(wp), intent(out) :: bounding_distance - ! TODO :: FILL ME IN + ! TODO :: FILL ME IN end subroutine s_get_bounding_box_corner_distance - !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid) + !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid) !! and before s_apply_ib_patches or s_apply_levelset. subroutine s_initialize_ib_airfoils() - integer :: i, j, airfoil_id - integer :: Np, Np1, Np2 - real(wp) :: ca_in, pa, ma, ta - real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c - - do i = 1, num_ibs - if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle - - airfoil_id = patch_ib(i)%airfoil_id - ca_in = ib_airfoil(airfoil_id)%c - pa = ib_airfoil(airfoil_id)%p - ma = ib_airfoil(airfoil_id)%m - ta = ib_airfoil(airfoil_id)%t - - Np1 = int((pa*ca_in/dx(0))*20) - Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) - Np = Np1 + Np2 + 1 - ib_airfoil_grids(airfoil_id)%Np = Np - $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%Np]') - - if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then - @:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np)) - @:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np)) - - ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp - ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp - ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp - ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp - - do j = 1, Np1 + Np2 - 1 - if (j <= Np1) then - xc = j*(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 + (j - 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 - - xu = (xa - yt*sin_c)*ca_in - yu = (yc + yt*cos_c)*ca_in - xl = (xa + yt*sin_c)*ca_in - yl = (yc - yt*cos_c)*ca_in - - ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu - ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu - ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl - ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl - end do - - ib_airfoil_grids(airfoil_id)%upper(Np)%x = ca_in - ib_airfoil_grids(airfoil_id)%upper(Np)%y = 0._wp - ib_airfoil_grids(airfoil_id)%lower(Np)%x = ca_in - ib_airfoil_grids(airfoil_id)%lower(Np)%y = 0._wp - - $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]') - end if - end do - - end subroutine s_initialize_ib_airfoils - - end module m_patch_geometries + integer :: i, j, airfoil_id + integer :: Np, Np1, Np2 + real(wp) :: ca_in, pa, ma, ta + real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c + + do i = 1, num_ibs + if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle + + airfoil_id = patch_ib(i)%airfoil_id + ca_in = ib_airfoil(airfoil_id)%c + pa = ib_airfoil(airfoil_id)%p + ma = ib_airfoil(airfoil_id)%m + ta = ib_airfoil(airfoil_id)%t + + Np1 = int((pa*ca_in/dx(0))*20) + Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) + Np = Np1 + Np2 + 1 + ib_airfoil_grids(airfoil_id)%Np = Np + $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%Np]') + + if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then + @:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np)) + @:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np)) + + ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp + ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp + ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp + ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp + + do j = 1, Np1 + Np2 - 1 + if (j <= Np1) then + xc = j*(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 + (j - 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 + + xu = (xa - yt*sin_c)*ca_in + yu = (yc + yt*cos_c)*ca_in + xl = (xa + yt*sin_c)*ca_in + yl = (yc - yt*cos_c)*ca_in + + ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu + ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu + ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl + ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl + end do + + ib_airfoil_grids(airfoil_id)%upper(Np)%x = ca_in + ib_airfoil_grids(airfoil_id)%upper(Np)%y = 0._wp + ib_airfoil_grids(airfoil_id)%lower(Np)%x = ca_in + ib_airfoil_grids(airfoil_id)%lower(Np)%y = 0._wp + + $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]') + end if + end do + + end subroutine s_initialize_ib_airfoils + +end module m_patch_geometries diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 18ce5d188f..af074a3724 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -30,24 +30,25 @@ contains impure subroutine s_apply_ib_patches(ib_markers) type(integer_field), intent(inout) :: ib_markers - + if (many_ib_patch_parallelism) then - call s_apply_ib_patches_ib_parallelism(ib_markers) + call s_apply_ib_patches_ib_parallelism(ib_markers) else - call s_apply_ib_patches_grid_cell_parallelism(ib_markers) + call s_apply_ib_patches_grid_cell_parallelism(ib_markers) end if end subroutine s_apply_ib_patches subroutine s_apply_ib_patches_grid_cell_parallelism(ib_markers) - type(integer_field), intent(inout) :: ib_markers - integer :: patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators + type(integer_field), intent(inout) :: ib_markers + integer :: patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds - real(wp), dimension(3) :: center, xyz_local - real(wp) :: bounding_box_corner_distance + real(wp), dimension(3) :: center, xyz_local + real(wp) :: bounding_box_corner_distance ! 3D Patch Geometries + if (num_dims == 3) then !> IB Patches !> @{ @@ -71,11 +72,15 @@ contains ir = m + gp_layers + 1 jr = n + gp_layers + 1 kr = p + gp_layers + 1 - call get_bounding_indices(center(1) - bounding_box_corner_distance, center(1) + bounding_box_corner_distance, x_cc, il, ir) - call get_bounding_indices(center(2) - bounding_box_corner_distance, center(2) + bounding_box_corner_distance, y_cc, jl, jr) - call get_bounding_indices(center(3) - bounding_box_corner_distance, center(3) + bounding_box_corner_distance, z_cc, kl, kr) - - $:GPU_PARALLEL_LOOP(private='[i, j, k]', copyin='[patch_id, encoded_patch_id, center]', collapse=3) + call get_bounding_indices(center(1) - bounding_box_corner_distance, & + & center(1) + bounding_box_corner_distance, x_cc, il, ir) + call get_bounding_indices(center(2) - bounding_box_corner_distance, & + & center(2) + bounding_box_corner_distance, y_cc, jl, jr) + call get_bounding_indices(center(3) - bounding_box_corner_distance, & + & center(3) + bounding_box_corner_distance, z_cc, kl, kr) + + $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[patch_id, encoded_patch_id, center]', & + & collapse=3) do k = kl, kr do j = jl, jr do i = il, ir @@ -83,19 +88,25 @@ contains xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - xyz_local = xyz_local - patch_ib(patch_id)%offset ! airfoils are a patch that require a centroid offset + ! airfoils are a patch that require a centroid offset + xyz_local = xyz_local - patch_ib(patch_id)%offset ! perform the interior check for the patch geometry of this IB if (patch_ib(patch_id)%geometry == 8) then - if (f_is_inside_sphere(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_sphere(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 9) then - if (f_is_inside_cuboid(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_cuboid(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(ipatch_id)%geometry == 10) then - if (f_is_inside_cylinder(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_cylinder(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 11) then - if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 12) then - if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id end if end do end do @@ -108,7 +119,7 @@ contains !> @} ! 2D Patch Geometries - else if (n > 0) then + else if (num_dims == 2) then !> IB Patches !> @{ call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper) @@ -128,29 +139,35 @@ contains jl = -gp_layers - 1 ir = m + gp_layers + 1 jr = n + gp_layers + 1 - call get_bounding_indices(center(1) - bounding_box_corner_distance, center(1) + bounding_box_corner_distance, x_cc, il, ir) - call get_bounding_indices(center(2) - bounding_box_corner_distance, center(2) + bounding_box_corner_distance, y_cc, jl, jr) + call get_bounding_indices(center(1) - bounding_box_corner_distance, & + & center(1) + bounding_box_corner_distance, x_cc, il, ir) + call get_bounding_indices(center(2) - bounding_box_corner_distance, & + & center(2) + bounding_box_corner_distance, y_cc, jl, jr) - $:GPU_PARALLEL_LOOP(private='[i, j, k]', copyin='[patch_id, encoded_patch_id, center]', collapse=2) + $:GPU_PARALLEL_LOOP(private='[i, j, xyz_local]', copyin='[patch_id, encoded_patch_id, center]', collapse=2) do j = jl, jr do i = il, ir ! get coordinate frame centered on IB xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - xyz_local = xyz_local - patch_ib(patch_id)%offset ! airfoils are a patch that require a centroid offset + ! airfoils are a patch that require a centroid offset + xyz_local = xyz_local - patch_ib(patch_id)%offset ! perform the interior check for the patch geometry of this IB if (patch_ib(patch_id)%geometry == 2) then if (f_is_inside_circle(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 3) then - if (f_is_inside_rectangle(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_rectangle(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(ipatch_id)%geometry == 4) then - if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then - if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 6) then - if (f_is_inside_ellipse(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_ellipse(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id end if end do end do @@ -165,6 +182,73 @@ contains subroutine s_apply_ib_patches_ib_parallelism() + type(integer_field), intent(inout) :: ib_markers + integer :: patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators + integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds + real(wp), dimension(3) :: center, xyz_local + real(wp) :: bounding_box_corner_distance + + if (num_dims == 3) then + ! get the periodicities + call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper) + else if (num_dims == 2) then + ! get the periodicities considered + call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper) + + $:GPU_PARALLEL_LOOP(private='[xp, yp, patch_id, center]') + do xp = xp_lower, xp_upper + do yp = yp_lower, yp_upper + do patch_id = 1, num_ibs + 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) = 0._wp + call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) + + ! 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 + il = -gp_layers - 1 + jl = -gp_layers - 1 + ir = m + gp_layers + 1 + jr = n + gp_layers + 1 + call get_bounding_indices(center(1) - bounding_box_corner_distance, & + & center(1) + bounding_box_corner_distance, x_cc, il, ir) + call get_bounding_indices(center(2) - bounding_box_corner_distance, & + & center(2) + bounding_box_corner_distance, y_cc, jl, jr) + + do j = jl, jr + do i = il, ir + ! get coordinate frame centered on IB + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] + ! rotate the frame into the IB's coordinates + xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) + ! airfoils are a patch that require a centroid offset + xyz_local = xyz_local - patch_ib(patch_id)%offset + + ! perform the interior check for the patch geometry of this IB + if (patch_ib(patch_id)%geometry == 2) then + if (f_is_inside_circle(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 3) then + if (f_is_inside_rectangle(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id + else if (patch_ib(ipatch_id)%geometry == 4) then + if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + & k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 5) then + if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + & k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 6) then + if (f_is_inside_ellipse(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + end if + end do + end do + end do + end do + end do + $:END_GPU_PARALLEL_LOOP + end if + end subroutine s_apply_ib_patches_ib_parallelism !> The STL patch is a 2D geometry that is imported from an STL file. From ec363dfd10079df5c0dca015b6680c0b37f27850 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 09:01:18 -0400 Subject: [PATCH 04/14] Unified some of the patch geometry functions. Finished out the first draft of the grid parallelism subroutine. Deprecated cylinder value to only favor length_x being used --- src/common/m_patch_geometries.fpp | 81 +++++---------- src/simulation/m_ib_patches.fpp | 166 ++++++++++++++++++++++++------ 2 files changed, 158 insertions(+), 89 deletions(-) diff --git a/src/common/m_patch_geometries.fpp b/src/common/m_patch_geometries.fpp index 3734eb8b55..f1eaeea974 100644 --- a/src/common/m_patch_geometries.fpp +++ b/src/common/m_patch_geometries.fpp @@ -18,106 +18,75 @@ module m_patch_geometries implicit none - public :: f_is_inside_circle, f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_cuboid, f_is_inside_airfoil, & - & f_is_inside_ellipse + public :: f_is_inside_sphere, f_is_inside_cylinder, f_is_inside_cuboid, f_is_inside_airfoil, f_is_inside_ellipse contains - !> Check if the x, and y coordinates would be located inside a circle with the patch_id's radius - function f_is_inside_circle(patch_id, x, y) result(is_inside) - - $:GPU_ROUTINE(parallelism='[seq]') - - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y - logical :: is_inside - - is_inside = x**2 + y**2 <= patch_ib(patch_id)%radius**2 - - end function f_is_inside_circle - !> Check if the x, y, and z coordinates would be located inside a sphere with the patch_id's radius - function f_is_inside_sphere(patch_id, x, y, z) result(is_inside) + function f_is_inside_sphere(x, y, z, radius) result(is_inside) $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z + real(wp), intent(in) :: radius, x, y, z logical :: is_inside - is_inside = x**2 + y**2 + z**2 <= patch_ib(patch_id)%radius**2 + is_inside = x**2 + y**2 + z**2 <= radius**2 end function f_is_inside_sphere !> Check which length of the cylinder is not default. Use that direction as the height and the other two coordinate ! values as the radius check - function f_is_inside_cylinder(patch_id, x, y, z) result(is_inside) + function f_is_inside_cylinder(polar_x, polar_y, height, radius, length) result(is_inside) $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z + real(wp), intent(in) :: x, y, z, radius, length logical :: is_inside - ! check if the cylinder is extended in the x direction - is_inside = (.not. f_is_default(patch_ib(patch_id)%length_x)) .and. y**2 + z**2 <= patch_ib(patch_id)%radius**2 .and. & - & -0.5_wp*patch_ib(patch_id)%length_x <= x .and. 0.5_wp*patch_ib(patch_id)%length_x >= x - - ! check if the cylinder is extended in the z direction - is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_y)) .and. x**2 & - & + z**2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_y <= y & - & .and. 0.5_wp*patch_ib(patch_id)%length_y >= y) - - ! check if the cylinder is extended in the z direction - is_inside = is_inside .or. ((.not. f_is_default(patch_ib(patch_id)%length_z)) .and. x**2 & - & + y**2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_z <= z & - & .and. 0.5_wp*patch_ib(patch_id)%length_z >= z) + ! check if the circular component of the cylinder is correct + is_inside = polar_x**2 + polar_y**2 <= patch_ib(patch_id)%radius**2 - ! TODO :: This can easily reuse the f_is_inside_circle function with an additional length requirement + ! in 3D, also check the length of the cylinder + if (num_dims == 3) is_inside = is_inside .and. -0.5_wp*length <= height .and. 0.5_wp*length >= height end function f_is_inside_cylinder !> Check if the x, y, and possibly z coordinates would be located inside a cuboid with the patch_id's lengths - function f_is_inside_cuboid(patch_id, x, y, z) result(is_inside) + function f_is_inside_cuboid(x, y, z, length) result(is_inside) $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z - logical :: is_inside + real(wp), intent(in) :: x, y, z + real(wp), dimension(3), intent(in) :: length + logical :: is_inside ! check if x and y are inside the rectangle plane at z=0 - is_inside = -0.5_wp*patch_ib(patch_id)%length_x <= x .and. 0.5_wp*patch_ib(patch_id)%length_x >= x is_inside = is_inside & - & .and. -0.5_wp*patch_ib(patch_id)%length_y <= y .and. 0.5_wp*patch_ib(patch_id)%length_y >= y + is_inside = -0.5_wp*length(1) <= x .and. 0.5_wp*length(1) >= x .and. -0.5_wp*length(2) <= y .and. 0.5_wp*length(2) >= y ! if we are in 3D, this is a cuboid and so we must also check the z axis - if (num_dims == 3) is_inside = is_inside .and. -0.5_wp*patch_ib(patch_id)%length_z <= z & - & .and. 0.5_wp*patch_ib(patch_id)%length_z >= z + if (num_dims == 3) is_inside = is_inside .and. -0.5_wp*length(3) <= z .and. 0.5_wp*length(3) >= z end function f_is_inside_cuboid !> Check if the x, y, are bounded by a NACA airfoil. Check if the z coordinate is inside the left and right edges of the !! airfoil, if set. - function f_is_inside_airfoil(patch_id, x, y, z) result(is_inside) + function f_is_inside_airfoil(x, y, z, length, airfoil_id) result(is_inside) $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y, z + real(wp), intent(in) :: x, y, z, length + integer, intent(in) :: airfoil_id logical :: is_inside - integer :: k, airfoil_id + integer :: k real(wp) :: f - airfoil_id = patch_ib(patch_id)%airfoil_id - is_inside = .false. ! check the initial x bounds of the grid cell if (.not. (x >= 0._wp .and. x <= ib_airfoil(airfoil_id)%c)) return ! if we are in 3D, we must also check the z axis - if (num_dims == 3 .and. (.not. (-0.5_wp*patch_ib(patch_id)%length_z <= z .and. 0.5_wp*patch_ib(patch_id)%length_z >= z))) & - & return + if (num_dims == 3 .and. (.not. (-0.5_wp*length(3) <= z .and. 0.5_wp*length(3) >= z))) return ! our check branches for the upper and lower half of the airfoil if (y >= 0._wp) then @@ -158,13 +127,13 @@ contains end function f_is_inside_airfoil - function f_is_inside_ellipse(patch_id, x, y) result(is_inside) + function f_is_inside_ellipse(x, y, length) result(is_inside) $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(in) :: x, y - logical :: is_inside + real(wp), intent(in) :: x, y + real(wp), dimension(3), intent(in) :: length + logical :: is_inside end function f_is_inside_ellipse diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index af074a3724..3871927b1f 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -44,14 +44,12 @@ contains type(integer_field), intent(inout) :: ib_markers integer :: patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds - real(wp), dimension(3) :: center, xyz_local - real(wp) :: bounding_box_corner_distance + real(wp), dimension(3) :: center, xyz_local, length + real(wp) :: bounding_box_corner_distance, radius ! 3D Patch Geometries if (num_dims == 3) then - !> IB Patches - !> @{ call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper) do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper @@ -79,8 +77,8 @@ contains call get_bounding_indices(center(3) - bounding_box_corner_distance, & & center(3) + bounding_box_corner_distance, z_cc, kl, kr) - $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[patch_id, encoded_patch_id, center]', & - & collapse=3) + $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, length, radius]', copyin='[patch_id, & + & encoded_patch_id, center]', collapse=3) do k = kl, kr do j = jl, jr do i = il, ir @@ -88,23 +86,33 @@ contains xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - ! airfoils are a patch that require a centroid offset - xyz_local = xyz_local - patch_ib(patch_id)%offset ! perform the interior check for the patch geometry of this IB if (patch_ib(patch_id)%geometry == 8) then - if (f_is_inside_sphere(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + ! sphere geometry + radius = patch_ib(patch_id)%radius + if (f_is_inside_sphere(x_cc(i), y_cc(j), z_cc(k), radius)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 9) then - if (f_is_inside_cuboid(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + ! cuboid geometry + length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, & + & patch_ib(patch_id)%length_z] + if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(ipatch_id)%geometry == 10) then - if (f_is_inside_cylinder(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & - & k) = encoded_patch_id + ! cylinder geometry + radius = radius = patch_ib(patch_id)%radius + if (f_is_inside_cylinder(x_cc(i), y_cc(j), z_cc(k), radius, & + & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 11) then - if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & - & k) = encoded_patch_id + ! 3D airfoil geometry + airfoil_id = patch_ib(patch_id)%airfoil_id + xyz_local = xyz_local - patch_ib(patch_id)%offset + if (f_is_inside_airfoil(x_cc(i), y_cc(j), z_cc(k), patch_ib(patch_id)%length_z, & + & airoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 12) then + ! STL model geometry + xyz_local = xyz_local - patch_ib(patch_id)%offset if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & & k) = encoded_patch_id end if @@ -116,12 +124,9 @@ contains end do end do end do - !> @} ! 2D Patch Geometries else if (num_dims == 2) then - !> IB Patches - !> @{ call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper) do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper @@ -151,23 +156,33 @@ contains xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - ! airfoils are a patch that require a centroid offset - xyz_local = xyz_local - patch_ib(patch_id)%offset ! perform the interior check for the patch geometry of this IB if (patch_ib(patch_id)%geometry == 2) then - if (f_is_inside_circle(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + ! circular geometries + radius = patch_ib(patch_id)%radius + if (f_is_inside_cylinder(0._wp, x_cc(i), y_cc(j), radius, 0._wp)) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 3) then - if (f_is_inside_rectangle(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + ! rectangular geometries + length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] + if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(ipatch_id)%geometry == 4) then - if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + ! 2D airfoil geometry + airfoil_id = patch_ib(patch_id)%airfoil_id + xyz_local = xyz_local - patch_ib(patch_id)%offset + if (f_is_inside_airfoil(x_cc(i), y_cc(j), 0._wp, 0._wp, airfoil_id)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then + ! STL model geometry + xyz_local = xyz_local - patch_ib(patch_id)%offset if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 6) then - if (f_is_inside_ellipse(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + ! ellipse geometry + length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] + if (f_is_inside_ellipse(x_cc(i), y_cc(j), length)) ib_markers%sf(i, j, k) = encoded_patch_id end if end do end do @@ -175,7 +190,6 @@ contains end do end do end do - !> @} end if end subroutine s_apply_ib_patches_grid_cell_parallelism @@ -191,13 +205,88 @@ contains if (num_dims == 3) then ! get the periodicities call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper) + + do xp = xp_lower, xp_upper + do yp = yp_lower, yp_upper + do zp = zp_lower, zp_upper + $:GPU_PARALLEL_LOOP(private='[xp, yp, zp, i, il, ir, j, jl, jr, k, kl, kr, xyz_local, length, radius, & + & bounding_box_corner_distance, patch_id, encoded_patch_id, center]') + do patch_id = 1, num_ibs + 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) + call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) + + ! 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) + + ! find the indices to the left and right of the IB in i, j, k + il = -gp_layers - 1 + jl = -gp_layers - 1 + kl = -gp_layers - 1 + ir = m + gp_layers + 1 + jr = n + gp_layers + 1 + kr = p + gp_layers + 1 + call get_bounding_indices(center(1) - bounding_box_corner_distance, & + & center(1) + bounding_box_corner_distance, x_cc, il, ir) + call get_bounding_indices(center(2) - bounding_box_corner_distance, & + & center(2) + bounding_box_corner_distance, y_cc, jl, jr) + call get_bounding_indices(center(3) - bounding_box_corner_distance, & + & center(3) + bounding_box_corner_distance, z_cc, kl, kr) + do k = kl, kr + do j = jl, jr + do i = il, ir + ! get coordinate frame centered on IB + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] + ! rotate the frame into the IB's coordinates + xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) + + ! perform the interior check for the patch geometry of this IB + if (patch_ib(patch_id)%geometry == 8) then + ! sphere geometry + radius = patch_ib(patch_id)%radius + if (f_is_inside_sphere(x_cc(i), y_cc(j), z_cc(k), radius)) ib_markers%sf(i, j, & + & k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 9) then + ! cuboid geometry + length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, & + & patch_ib(patch_id)%length_z] + if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & + & k) = encoded_patch_id + else if (patch_ib(ipatch_id)%geometry == 10) then + ! cylinder geometry + radius = radius = patch_ib(patch_id)%radius + if (f_is_inside_cylinder(x_cc(i), y_cc(j), z_cc(k), radius, & + & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 11) then + ! 3D airfoil geometry + airfoil_id = patch_ib(patch_id)%airfoil_id + xyz_local = xyz_local - patch_ib(patch_id)%offset + if (f_is_inside_airfoil(x_cc(i), y_cc(j), z_cc(k), patch_ib(patch_id)%length_z, & + & airoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 12) then + ! STL model geometry + xyz_local = xyz_local - patch_ib(patch_id)%offset + if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + & k) = encoded_patch_id + end if + end do + end do + end do + end do + $:END_GPU_PARALLEL_LOOP() + end do + end do + end do else if (num_dims == 2) then - ! get the periodicities considered + ! get the periodicities call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper) - $:GPU_PARALLEL_LOOP(private='[xp, yp, patch_id, center]') + $:GPU_PARALLEL_LOOP(private='[xp, yp, patch_id, center]', collapse=3) do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper + $:GPU_PARALLEL_LOOP(private='[xp, yp, i, il, ir, j, jl, jr, xyz_local, length, radius, & + & bounding_box_corner_distance, patch_id, encoded_patch_id, center]') do patch_id = 1, num_ibs 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) @@ -223,30 +312,41 @@ contains xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - ! airfoils are a patch that require a centroid offset - xyz_local = xyz_local - patch_ib(patch_id)%offset ! perform the interior check for the patch geometry of this IB if (patch_ib(patch_id)%geometry == 2) then - if (f_is_inside_circle(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + ! circular geometries + radius = patch_ib(patch_id)%radius + if (f_is_inside_cylinder(x_cc(i), y_cc(j), 0._wp, radius, 0._wp)) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 3) then - if (f_is_inside_rectangle(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + ! rectangular geometries + length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] + if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(ipatch_id)%geometry == 4) then - if (f_is_inside_airfoil(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + ! 2D airfoil geometry + airfoil_id = patch_ib(patch_id)%airfoil_id + xyz_local = xyz_local - patch_ib(patch_id)%offset + if (f_is_inside_airfoil(x_cc(i), y_cc(j), 0._wp, 0._wp, airfoil_id)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then + ! STL model geometry + model_id = patch_ib(patch_id)%model_id + xyz_local = xyz_local - patch_ib(patch_id)%offset if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 6) then - if (f_is_inside_ellipse(patch_id, x_cc(i), y_cc(j))) ib_markers%sf(i, j, k) = encoded_patch_id + ! ellipse geometry + length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] + if (f_is_inside_ellipse(x_cc(i), y_cc(j), length)) ib_markers%sf(i, j, k) = encoded_patch_id end if end do end do end do + $:END_GPU_PARALLEL_LOOP() end do end do - $:END_GPU_PARALLEL_LOOP end if end subroutine s_apply_ib_patches_ib_parallelism From 72529c9724d9566878834cdb899713f941c161bb Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 10:29:27 -0400 Subject: [PATCH 05/14] Passes circle tests --- src/common/m_patch_geometries.fpp | 92 +----------- src/post_process/m_global_parameters.fpp | 3 + src/pre_process/m_global_parameters.fpp | 5 +- src/simulation/m_ib_patches.fpp | 183 ++++++++++++++++++----- 4 files changed, 154 insertions(+), 129 deletions(-) diff --git a/src/common/m_patch_geometries.fpp b/src/common/m_patch_geometries.fpp index f1eaeea974..2709e89a6b 100644 --- a/src/common/m_patch_geometries.fpp +++ b/src/common/m_patch_geometries.fpp @@ -40,11 +40,11 @@ contains $:GPU_ROUTINE(parallelism='[seq]') - real(wp), intent(in) :: x, y, z, radius, length + real(wp), intent(in) :: polar_x, polar_y, height, radius, length logical :: is_inside ! check if the circular component of the cylinder is correct - is_inside = polar_x**2 + polar_y**2 <= patch_ib(patch_id)%radius**2 + is_inside = polar_x**2 + polar_y**2 <= radius**2 ! in 3D, also check the length of the cylinder if (num_dims == 3) is_inside = is_inside .and. -0.5_wp*length <= height .and. 0.5_wp*length >= height @@ -86,7 +86,7 @@ contains if (.not. (x >= 0._wp .and. x <= ib_airfoil(airfoil_id)%c)) return ! if we are in 3D, we must also check the z axis - if (num_dims == 3 .and. (.not. (-0.5_wp*length(3) <= z .and. 0.5_wp*length(3) >= z))) return + if (num_dims == 3 .and. (.not. (-0.5_wp*length <= z .and. 0.5_wp*length >= z))) return ! our check branches for the upper and lower half of the airfoil if (y >= 0._wp) then @@ -135,89 +135,9 @@ contains real(wp), dimension(3), intent(in) :: length logical :: is_inside - end function f_is_inside_ellipse - - subroutine s_get_bounding_box_corner_distance(patch_id, bounding_distance) - - $:GPU_ROUTINE(parallelism='[seq]') + ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1 + is_inside = (x/(0.5_wp*length(1)))**2 + (y/(0.5_wp*length(2)))**2 <= 1._wp - integer, intent(in) :: patch_id - real(wp), intent(out) :: bounding_distance - - ! TODO :: FILL ME IN - - end subroutine s_get_bounding_box_corner_distance - - !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid) - !! and before s_apply_ib_patches or s_apply_levelset. - subroutine s_initialize_ib_airfoils() - - integer :: i, j, airfoil_id - integer :: Np, Np1, Np2 - real(wp) :: ca_in, pa, ma, ta - real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c - - do i = 1, num_ibs - if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle - - airfoil_id = patch_ib(i)%airfoil_id - ca_in = ib_airfoil(airfoil_id)%c - pa = ib_airfoil(airfoil_id)%p - ma = ib_airfoil(airfoil_id)%m - ta = ib_airfoil(airfoil_id)%t - - Np1 = int((pa*ca_in/dx(0))*20) - Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) - Np = Np1 + Np2 + 1 - ib_airfoil_grids(airfoil_id)%Np = Np - $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%Np]') - - if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then - @:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np)) - @:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np)) - - ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp - ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp - ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp - ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp - - do j = 1, Np1 + Np2 - 1 - if (j <= Np1) then - xc = j*(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 + (j - 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 - - xu = (xa - yt*sin_c)*ca_in - yu = (yc + yt*cos_c)*ca_in - xl = (xa + yt*sin_c)*ca_in - yl = (yc - yt*cos_c)*ca_in - - ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu - ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu - ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl - ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl - end do - - ib_airfoil_grids(airfoil_id)%upper(Np)%x = ca_in - ib_airfoil_grids(airfoil_id)%upper(Np)%y = 0._wp - ib_airfoil_grids(airfoil_id)%lower(Np)%x = ca_in - ib_airfoil_grids(airfoil_id)%lower(Np)%y = 0._wp - - $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]') - end if - end do - - end subroutine s_initialize_ib_airfoils + end function f_is_inside_ellipse end module m_patch_geometries diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index b9673f354c..574b01cd0b 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -106,6 +106,9 @@ module m_global_parameters integer, dimension(3, 2) :: shear_BC_flip_indices !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num) integer, allocatable, dimension(:) :: proc_coords !< Processor coordinates in MPI_CART_COMM integer, allocatable, dimension(:) :: start_idx !< Starting cell-center index of local processor in global grid + type(ib_airfoil_parameters), allocatable, dimension(:) :: ib_airfoil !< Per-airfoil NACA parameters (unused in post_process) + type(ib_airfoil_grid), allocatable, dimension(:) :: ib_airfoil_grids !< Per-airfoil computed surface grids (unused in post_process) + #ifdef MFC_MPI type(mpi_io_var), public :: MPI_IO_DATA type(mpi_io_ib_var), public :: MPI_IO_IB_DATA diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index bd28e73326..2c4bd9d4b3 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -94,8 +94,9 @@ module m_global_parameters !> @name Immersed Boundaries !> @{ type(ib_patch_parameters), dimension(num_ib_patches_max_namelist) :: patch_ib !< Immersed boundary patch parameters - type(ib_airfoil_parameters), dimension(num_ib_airfoils_max) :: ib_airfoil !< Per-airfoil NACA parameters - type(ib_stl_parameters), dimension(num_stl_models_max) :: stl_models !< Per-STL model parameters (namelist) + type(ib_airfoil_parameters), dimension(num_ib_airfoils_max) :: ib_airfoil !< Per-airfoil NACA parameters + type(ib_airfoil_grid), allocatable, dimension(:) :: ib_airfoil_grids !< Per-airfoil computed surface grids (unused in pre_process) + type(ib_stl_parameters), dimension(num_stl_models_max) :: stl_models !< Per-STL model parameters (namelist) !> @} !> @name Non-polytropic bubble gas compression diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 3871927b1f..081267c7d4 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -12,6 +12,7 @@ !> @brief Immersed boundary patch geometry constructors for 2D and 3D shapes module m_ib_patches + use m_patch_geometries use m_model ! Subroutine(s) related to STL files use m_derived_types ! Definitions of the derived types use m_global_parameters @@ -42,7 +43,7 @@ contains subroutine s_apply_ib_patches_grid_cell_parallelism(ib_markers) type(integer_field), intent(inout) :: ib_markers - integer :: patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators + integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds real(wp), dimension(3) :: center, xyz_local, length real(wp) :: bounding_box_corner_distance, radius @@ -77,13 +78,13 @@ contains call get_bounding_indices(center(3) - bounding_box_corner_distance, & & center(3) + bounding_box_corner_distance, z_cc, kl, kr) - $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, length, radius]', copyin='[patch_id, & + $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, length, radius, airfoil_id]', copyin='[patch_id, & & encoded_patch_id, center]', collapse=3) do k = kl, kr do j = jl, jr do i = il, ir ! get coordinate frame centered on IB - xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)] ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) @@ -99,22 +100,22 @@ contains & patch_ib(patch_id)%length_z] if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & & k) = encoded_patch_id - else if (patch_ib(ipatch_id)%geometry == 10) then + else if (patch_ib(patch_id)%geometry == 10) then ! cylinder geometry - radius = radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(x_cc(i), y_cc(j), z_cc(k), radius, & + radius = patch_ib(patch_id)%radius + if (f_is_inside_cylinder(y_cc(j), z_cc(k), x_cc(i), radius, & & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 11) then ! 3D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id - xyz_local = xyz_local - patch_ib(patch_id)%offset + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset if (f_is_inside_airfoil(x_cc(i), y_cc(j), z_cc(k), patch_ib(patch_id)%length_z, & - & airoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id + & airfoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry - xyz_local = xyz_local - patch_ib(patch_id)%offset - if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & - & k) = encoded_patch_id + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset + ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + ! & k) = encoded_patch_id end if end do end do @@ -148,8 +149,9 @@ contains & center(1) + bounding_box_corner_distance, x_cc, il, ir) call get_bounding_indices(center(2) - bounding_box_corner_distance, & & center(2) + bounding_box_corner_distance, y_cc, jl, jr) + print *, "Bounding indices: ", il, ir, jl, jr - $:GPU_PARALLEL_LOOP(private='[i, j, xyz_local]', copyin='[patch_id, encoded_patch_id, center]', collapse=2) + $:GPU_PARALLEL_LOOP(private='[i, j, xyz_local, airfoil_id]', copyin='[patch_id, encoded_patch_id, center]', collapse=2) do j = jl, jr do i = il, ir ! get coordinate frame centered on IB @@ -161,28 +163,29 @@ contains if (patch_ib(patch_id)%geometry == 2) then ! circular geometries radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(0._wp, x_cc(i), y_cc(j), radius, 0._wp)) ib_markers%sf(i, j, & + print *, i, j + if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 3) then ! rectangular geometries length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & - & k) = encoded_patch_id - else if (patch_ib(ipatch_id)%geometry == 4) then + & 0) = encoded_patch_id + else if (patch_ib(patch_id)%geometry == 4) then ! 2D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id - xyz_local = xyz_local - patch_ib(patch_id)%offset + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset if (f_is_inside_airfoil(x_cc(i), y_cc(j), 0._wp, 0._wp, airfoil_id)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + & 0) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry - xyz_local = xyz_local - patch_ib(patch_id)%offset - if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset + ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + ! & 0) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_ellipse(x_cc(i), y_cc(j), length)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_ellipse(x_cc(i), y_cc(j), length)) ib_markers%sf(i, j, 0) = encoded_patch_id end if end do end do @@ -194,13 +197,13 @@ contains end subroutine s_apply_ib_patches_grid_cell_parallelism - subroutine s_apply_ib_patches_ib_parallelism() + subroutine s_apply_ib_patches_ib_parallelism(ib_markers) type(integer_field), intent(inout) :: ib_markers - integer :: patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators + integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds - real(wp), dimension(3) :: center, xyz_local - real(wp) :: bounding_box_corner_distance + real(wp), dimension(3) :: center, xyz_local, length + real(wp) :: bounding_box_corner_distance, radius if (num_dims == 3) then ! get the periodicities @@ -210,7 +213,7 @@ contains do yp = yp_lower, yp_upper do zp = zp_lower, zp_upper $:GPU_PARALLEL_LOOP(private='[xp, yp, zp, i, il, ir, j, jl, jr, k, kl, kr, xyz_local, length, radius, & - & bounding_box_corner_distance, patch_id, encoded_patch_id, center]') + & bounding_box_corner_distance, patch_id, airfoil_id, model_id, encoded_patch_id, center]') do patch_id = 1, num_ibs 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) @@ -237,7 +240,7 @@ contains do j = jl, jr do i = il, ir ! get coordinate frame centered on IB - xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] + xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)] ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) @@ -253,22 +256,22 @@ contains & patch_ib(patch_id)%length_z] if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & & k) = encoded_patch_id - else if (patch_ib(ipatch_id)%geometry == 10) then + else if (patch_ib(patch_id)%geometry == 10) then ! cylinder geometry - radius = radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(x_cc(i), y_cc(j), z_cc(k), radius, & + radius = patch_ib(patch_id)%radius + if (f_is_inside_cylinder(y_cc(j), z_cc(k), x_cc(i), radius, & & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 11) then ! 3D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id - xyz_local = xyz_local - patch_ib(patch_id)%offset + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset if (f_is_inside_airfoil(x_cc(i), y_cc(j), z_cc(k), patch_ib(patch_id)%length_z, & - & airoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id + & airfoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry - xyz_local = xyz_local - patch_ib(patch_id)%offset - if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & - & k) = encoded_patch_id + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset + ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & + ! & k) = encoded_patch_id end if end do end do @@ -286,7 +289,7 @@ contains do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper $:GPU_PARALLEL_LOOP(private='[xp, yp, i, il, ir, j, jl, jr, xyz_local, length, radius, & - & bounding_box_corner_distance, patch_id, encoded_patch_id, center]') + & bounding_box_corner_distance, patch_id, airfoil_id, model_id, encoded_patch_id, center]') do patch_id = 1, num_ibs 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) @@ -324,18 +327,18 @@ contains length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & & k) = encoded_patch_id - else if (patch_ib(ipatch_id)%geometry == 4) then + else if (patch_ib(patch_id)%geometry == 4) then ! 2D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id - xyz_local = xyz_local - patch_ib(patch_id)%offset + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset if (f_is_inside_airfoil(x_cc(i), y_cc(j), 0._wp, 0._wp, airfoil_id)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry model_id = patch_ib(patch_id)%model_id - xyz_local = xyz_local - patch_ib(patch_id)%offset - if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset + ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & + ! & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] @@ -515,6 +518,78 @@ contains end subroutine s_ib_3d_model + !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid) + !! and before s_apply_ib_patches or s_apply_levelset. + subroutine s_initialize_ib_airfoils() + + integer :: i, j, airfoil_id + integer :: Np, Np1, Np2 + real(wp) :: ca_in, pa, ma, ta + real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c + + do i = 1, num_ibs + if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle + + airfoil_id = patch_ib(i)%airfoil_id + ca_in = ib_airfoil(airfoil_id)%c + pa = ib_airfoil(airfoil_id)%p + ma = ib_airfoil(airfoil_id)%m + ta = ib_airfoil(airfoil_id)%t + + Np1 = int((pa*ca_in/dx(0))*20) + Np2 = int(((ca_in - pa*ca_in)/dx(0))*20) + Np = Np1 + Np2 + 1 + ib_airfoil_grids(airfoil_id)%Np = Np + $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%Np]') + + if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then + @:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np)) + @:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np)) + + ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp + ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp + ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp + ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp + + do j = 1, Np1 + Np2 - 1 + if (j <= Np1) then + xc = j*(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 + (j - 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 + + xu = (xa - yt*sin_c)*ca_in + yu = (yc + yt*cos_c)*ca_in + xl = (xa + yt*sin_c)*ca_in + yl = (yc - yt*cos_c)*ca_in + + ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu + ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu + ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl + ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl + end do + + ib_airfoil_grids(airfoil_id)%upper(Np)%x = ca_in + ib_airfoil_grids(airfoil_id)%upper(Np)%y = 0._wp + ib_airfoil_grids(airfoil_id)%lower(Np)%x = ca_in + ib_airfoil_grids(airfoil_id)%lower(Np)%y = 0._wp + + $:GPU_UPDATE(device='[ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower]') + end if + end do + + end subroutine s_initialize_ib_airfoils + !> Compute a rotation matrix for converting to the rotating frame of the boundary subroutine s_update_ib_rotation_matrix(patch_id) @@ -562,6 +637,32 @@ contains end subroutine s_update_ib_rotation_matrix + subroutine s_get_bounding_box_corner_distance(patch_id, bounding_distance) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: patch_id + real(wp), intent(out) :: bounding_distance + + if (patch_ib(patch_id)%geometry == 2) then + ! circular geometries + bounding_distance = patch_ib(patch_id)%radius + else if (patch_ib(patch_id)%geometry == 3) then + ! rectangular geometries + bounding_distance = 0._wp + else if (patch_ib(patch_id)%geometry == 4) then + ! 2D airfoil geometry + bounding_distance = 0._wp + else if (patch_ib(patch_id)%geometry == 5) then + ! STL model geometry + bounding_distance = 0._wp + else if (patch_ib(patch_id)%geometry == 6) then + ! ellipse geometry + bounding_distance = 0._wp + end if + + end subroutine s_get_bounding_box_corner_distance + subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index) real(wp), intent(in) :: left_bound, right_bound From ce800e40c6918abc2356ac2f67242ad83dab0258 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 11:19:54 -0400 Subject: [PATCH 06/14] All tests pass, except for STL mdoels --- src/post_process/m_global_parameters.fpp | 13 +- src/pre_process/m_global_parameters.fpp | 7 +- src/simulation/m_ib_patches.fpp | 286 +++++++++++------------ 3 files changed, 145 insertions(+), 161 deletions(-) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 574b01cd0b..818f792a10 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -100,14 +100,15 @@ module m_global_parameters type(int_bounds_info) :: bc_x, bc_y, bc_z !> @} - integer :: shear_num !< Number of shear stress components - integer, dimension(3) :: shear_indices !< Indices of the stress components that represent shear stress - integer :: shear_BC_flip_num !< Number of shear stress components to reflect for boundary conditions - integer, dimension(3, 2) :: shear_BC_flip_indices !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num) + integer :: shear_num !< Number of shear stress components + integer, dimension(3) :: shear_indices !< Indices of the stress components that represent shear stress + integer :: shear_BC_flip_num !< Number of shear stress components to reflect for boundary conditions + integer, dimension(3, 2) :: shear_BC_flip_indices !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num) integer, allocatable, dimension(:) :: proc_coords !< Processor coordinates in MPI_CART_COMM integer, allocatable, dimension(:) :: start_idx !< Starting cell-center index of local processor in global grid - type(ib_airfoil_parameters), allocatable, dimension(:) :: ib_airfoil !< Per-airfoil NACA parameters (unused in post_process) - type(ib_airfoil_grid), allocatable, dimension(:) :: ib_airfoil_grids !< Per-airfoil computed surface grids (unused in post_process) + type(ib_airfoil_parameters), allocatable, dimension(:) :: ib_airfoil !< Per-airfoil NACA parameters (unused in post_process) + !> Per-airfoil computed surface grids (unused in post_process) + type(ib_airfoil_grid), allocatable, dimension(:) :: ib_airfoil_grids #ifdef MFC_MPI type(mpi_io_var), public :: MPI_IO_DATA diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 2c4bd9d4b3..804578a84c 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -94,9 +94,10 @@ module m_global_parameters !> @name Immersed Boundaries !> @{ type(ib_patch_parameters), dimension(num_ib_patches_max_namelist) :: patch_ib !< Immersed boundary patch parameters - type(ib_airfoil_parameters), dimension(num_ib_airfoils_max) :: ib_airfoil !< Per-airfoil NACA parameters - type(ib_airfoil_grid), allocatable, dimension(:) :: ib_airfoil_grids !< Per-airfoil computed surface grids (unused in pre_process) - type(ib_stl_parameters), dimension(num_stl_models_max) :: stl_models !< Per-STL model parameters (namelist) + type(ib_airfoil_parameters), dimension(num_ib_airfoils_max) :: ib_airfoil !< Per-airfoil NACA parameters + !> Per-airfoil computed surface grids (unused in pre_process) + type(ib_airfoil_grid), allocatable, dimension(:) :: ib_airfoil_grids + type(ib_stl_parameters), dimension(num_stl_models_max) :: stl_models !< Per-STL model parameters (namelist) !> @} !> @name Non-polytropic bubble gas compression diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 081267c7d4..30a16c12ff 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -43,10 +43,10 @@ contains subroutine s_apply_ib_patches_grid_cell_parallelism(ib_markers) type(integer_field), intent(inout) :: ib_markers - integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators - integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds - real(wp), dimension(3) :: center, xyz_local, length - real(wp) :: bounding_box_corner_distance, radius + integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators + integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds + real(wp), dimension(3) :: center, xyz_local, length + real(wp) :: radius ! 3D Patch Geometries @@ -59,24 +59,12 @@ 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) center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg) - call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) ! 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) ! find the indices to the left and right of the IB in i, j, k - il = -gp_layers - 1 - jl = -gp_layers - 1 - kl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - kr = p + gp_layers + 1 - call get_bounding_indices(center(1) - bounding_box_corner_distance, & - & center(1) + bounding_box_corner_distance, x_cc, il, ir) - call get_bounding_indices(center(2) - bounding_box_corner_distance, & - & center(2) + bounding_box_corner_distance, y_cc, jl, jr) - call get_bounding_indices(center(3) - bounding_box_corner_distance, & - & center(3) + bounding_box_corner_distance, z_cc, kl, kr) + call get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, length, radius, airfoil_id]', copyin='[patch_id, & & encoded_patch_id, center]', collapse=3) @@ -92,25 +80,26 @@ contains if (patch_ib(patch_id)%geometry == 8) then ! sphere geometry radius = patch_ib(patch_id)%radius - if (f_is_inside_sphere(x_cc(i), y_cc(j), z_cc(k), radius)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), & + & radius)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 9) then ! cuboid geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, & & patch_ib(patch_id)%length_z] - if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), & + & length)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 10) then ! cylinder geometry radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(y_cc(j), z_cc(k), x_cc(i), radius, & + if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, & & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 11) then ! 3D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - if (f_is_inside_airfoil(x_cc(i), y_cc(j), z_cc(k), patch_ib(patch_id)%length_z, & - & airfoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), & + & patch_ib(patch_id)%length_z, airfoil_id)) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset @@ -135,23 +124,15 @@ 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) center(3) = 0._wp - call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) ! 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 - il = -gp_layers - 1 - jl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - call get_bounding_indices(center(1) - bounding_box_corner_distance, & - & center(1) + bounding_box_corner_distance, x_cc, il, ir) - call get_bounding_indices(center(2) - bounding_box_corner_distance, & - & center(2) + bounding_box_corner_distance, y_cc, jl, jr) - print *, "Bounding indices: ", il, ir, jl, jr - - $:GPU_PARALLEL_LOOP(private='[i, j, xyz_local, airfoil_id]', copyin='[patch_id, encoded_patch_id, center]', collapse=2) + call get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) + + $:GPU_PARALLEL_LOOP(private='[i, j, xyz_local, airfoil_id]', copyin='[patch_id, encoded_patch_id, & + & center]', collapse=2) do j = jl, jr do i = il, ir ! get coordinate frame centered on IB @@ -163,20 +144,19 @@ contains if (patch_ib(patch_id)%geometry == 2) then ! circular geometries radius = patch_ib(patch_id)%radius - print *, i, j - if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, & + & j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 3) then ! rectangular geometries length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) ib_markers%sf(i, j, & & 0) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 4) then ! 2D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - if (f_is_inside_airfoil(x_cc(i), y_cc(j), 0._wp, 0._wp, airfoil_id)) ib_markers%sf(i, j, & - & 0) = encoded_patch_id + if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, & + & airfoil_id)) ib_markers%sf(i, j, 0) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset @@ -185,7 +165,8 @@ contains else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_ellipse(x_cc(i), y_cc(j), length)) ib_markers%sf(i, j, 0) = encoded_patch_id + if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) ib_markers%sf(i, j, & + & 0) = encoded_patch_id end if end do end do @@ -200,10 +181,10 @@ contains subroutine s_apply_ib_patches_ib_parallelism(ib_markers) type(integer_field), intent(inout) :: ib_markers - integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators - integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds - real(wp), dimension(3) :: center, xyz_local, length - real(wp) :: bounding_box_corner_distance, radius + integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators + integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds + real(wp), dimension(3) :: center, xyz_local, length + real(wp) :: radius if (num_dims == 3) then ! get the periodicities @@ -213,29 +194,18 @@ contains do yp = yp_lower, yp_upper do zp = zp_lower, zp_upper $:GPU_PARALLEL_LOOP(private='[xp, yp, zp, i, il, ir, j, jl, jr, k, kl, kr, xyz_local, length, radius, & - & bounding_box_corner_distance, patch_id, airfoil_id, model_id, encoded_patch_id, center]') + & patch_id, airfoil_id, model_id, encoded_patch_id, center]') do patch_id = 1, num_ibs 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) - call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) ! 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) ! find the indices to the left and right of the IB in i, j, k - il = -gp_layers - 1 - jl = -gp_layers - 1 - kl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - kr = p + gp_layers + 1 - call get_bounding_indices(center(1) - bounding_box_corner_distance, & - & center(1) + bounding_box_corner_distance, x_cc, il, ir) - call get_bounding_indices(center(2) - bounding_box_corner_distance, & - & center(2) + bounding_box_corner_distance, y_cc, jl, jr) - call get_bounding_indices(center(3) - bounding_box_corner_distance, & - & center(3) + bounding_box_corner_distance, z_cc, kl, kr) + call get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) + do k = kl, kr do j = jl, jr do i = il, ir @@ -248,25 +218,26 @@ contains if (patch_ib(patch_id)%geometry == 8) then ! sphere geometry radius = patch_ib(patch_id)%radius - if (f_is_inside_sphere(x_cc(i), y_cc(j), z_cc(k), radius)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), & + & radius)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 9) then ! cuboid geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, & & patch_ib(patch_id)%length_z] - if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), & + & length)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 10) then ! cylinder geometry radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(y_cc(j), z_cc(k), x_cc(i), radius, & + if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, & & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 11) then ! 3D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - if (f_is_inside_airfoil(x_cc(i), y_cc(j), z_cc(k), patch_ib(patch_id)%length_z, & - & airfoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), & + & patch_ib(patch_id)%length_z, airfoil_id)) ib_markers%sf(i, j, & + & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset @@ -294,20 +265,13 @@ 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) center(3) = 0._wp - call s_get_bounding_box_corner_distance(patch_id, bounding_box_corner_distance) ! 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 - il = -gp_layers - 1 - jl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - call get_bounding_indices(center(1) - bounding_box_corner_distance, & - & center(1) + bounding_box_corner_distance, x_cc, il, ir) - call get_bounding_indices(center(2) - bounding_box_corner_distance, & - & center(2) + bounding_box_corner_distance, y_cc, jl, jr) + ! find the indices to the left and right of the IB in i, j, k + call get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) do j = jl, jr do i = il, ir @@ -320,19 +284,19 @@ contains if (patch_ib(patch_id)%geometry == 2) then ! circular geometries radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(x_cc(i), y_cc(j), 0._wp, radius, 0._wp)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, & + & j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 3) then ! rectangular geometries length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_cuboid(x_cc(i), y_cc(j), z_cc(k), length)) ib_markers%sf(i, j, & + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) ib_markers%sf(i, j, & & k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 4) then ! 2D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - if (f_is_inside_airfoil(x_cc(i), y_cc(j), 0._wp, 0._wp, airfoil_id)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, & + & airfoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry model_id = patch_ib(patch_id)%model_id @@ -342,7 +306,8 @@ contains else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_ellipse(x_cc(i), y_cc(j), length)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) ib_markers%sf(i, j, & + & k) = encoded_patch_id end if end do end do @@ -387,30 +352,6 @@ contains jl = -gp_layers - 1 ir = m + gp_layers + 1 jr = n + gp_layers + 1 - - ! Local-space bounding box extents (min=1, max=2 in the third index) - lx(1) = stl_bounding_boxes(model_id, 1, 1) + offset(1) - lx(2) = stl_bounding_boxes(model_id, 1, 3) + offset(1) - ly(1) = stl_bounding_boxes(model_id, 2, 1) + offset(2) - ly(2) = stl_bounding_boxes(model_id, 2, 3) + offset(2) - - bbox_min = 1e12 - bbox_max = -1e12 - ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB - do cx = 1, 2 - do cy = 1, 2 - local_corner = [lx(cx), ly(cy), 0._wp] - world_corner = matmul(rotation, local_corner) + center - bbox_min(1) = min(bbox_min(1), world_corner(1)) - bbox_min(2) = min(bbox_min(2), world_corner(2)) - bbox_max(1) = max(bbox_max(1), world_corner(1)) - bbox_max(2) = max(bbox_max(2), world_corner(2)) - end do - end do - - call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir) - call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr) - $:GPU_PARALLEL_LOOP(private='[i, j, xy_local, eta]', copyin='[patch_id, model_id, encoded_patch_id, center, & & inverse_rotation, offset, threshold]', collapse=2) do i = il, ir @@ -467,36 +408,6 @@ contains jr = n + gp_layers + 1 kr = p + gp_layers + 1 - ! Local-space bounding box extents (min=1, max=2 in the third index) - lx(1) = stl_bounding_boxes(model_id, 1, 1) + offset(1) - lx(2) = stl_bounding_boxes(model_id, 1, 3) + offset(1) - ly(1) = stl_bounding_boxes(model_id, 2, 1) + offset(2) - ly(2) = stl_bounding_boxes(model_id, 2, 3) + offset(2) - lz(1) = stl_bounding_boxes(model_id, 3, 1) + offset(3) - lz(2) = stl_bounding_boxes(model_id, 3, 3) + offset(3) - - bbox_min = 1e12 - bbox_max = -1e12 - ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB - do cx = 1, 2 - do cy = 1, 2 - do cz = 1, 2 - local_corner = [lx(cx), ly(cy), lz(cz)] - world_corner = matmul(rotation, local_corner) + center - bbox_min(1) = min(bbox_min(1), world_corner(1)) - bbox_min(2) = min(bbox_min(2), world_corner(2)) - bbox_min(3) = min(bbox_min(3), world_corner(3)) - bbox_max(1) = max(bbox_max(1), world_corner(1)) - bbox_max(2) = max(bbox_max(2), world_corner(2)) - bbox_max(3) = max(bbox_max(3), world_corner(3)) - end do - end do - end do - - call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir) - call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr) - call get_bounding_indices(bbox_min(3), bbox_max(3), z_cc, kl, kr) - $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, eta]', copyin='[patch_id, model_id, encoded_patch_id, center, & & inverse_rotation, offset, threshold]', collapse=3) do i = il, ir @@ -637,33 +548,104 @@ contains end subroutine s_update_ib_rotation_matrix - subroutine s_get_bounding_box_corner_distance(patch_id, bounding_distance) + subroutine get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) $:GPU_ROUTINE(parallelism='[seq]') - integer, intent(in) :: patch_id - real(wp), intent(out) :: bounding_distance + real(wp), dimension(3), intent(in) :: center + integer, intent(in) :: patch_id + integer, intent(out) :: il, ir, jl, jr, kl, kr + real(wp), dimension(3) :: bbox_min, bbox_max, local_corner, world_corner + real(wp), dimension(2) :: lx, ly, lz + integer :: cx, cy, cz - if (patch_ib(patch_id)%geometry == 2) then - ! circular geometries - bounding_distance = patch_ib(patch_id)%radius + if (patch_ib(patch_id)%geometry == 2 .or. patch_ib(patch_id)%geometry == 8) then + ! circle and sphere geometries + bbox_min = center - patch_ib(patch_id)%radius + bbox_max = center + patch_ib(patch_id)%radius else if (patch_ib(patch_id)%geometry == 3) then ! rectangular geometries - bounding_distance = 0._wp - else if (patch_ib(patch_id)%geometry == 4) then - ! 2D airfoil geometry - bounding_distance = 0._wp + bbox_min = center - 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2) + bbox_max = center + 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2) + else if (patch_ib(patch_id)%geometry == 4 .or. patch_ib(patch_id)%geometry == 11) then + ! airfoil geometries TODO :: This can be better optimized, since airfoils are typically very long in one dimension + bbox_min = center - ib_airfoil(patch_ib(patch_id)%airfoil_id)%c + bbox_max = center + ib_airfoil(patch_ib(patch_id)%airfoil_id)%c else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry - bounding_distance = 0._wp + lx(1) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 1, 1) + patch_ib(patch_id)%centroid_offset(1) + lx(2) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 1, 3) + patch_ib(patch_id)%centroid_offset(1) + ly(1) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 2, 1) + patch_ib(patch_id)%centroid_offset(2) + ly(2) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 2, 3) + patch_ib(patch_id)%centroid_offset(2) + + bbox_min = 1e12 + bbox_max = -1e12 + ! Enumerate all 4 corners of the local bounding box, rotate to world space, track world-space AABB + do cx = 1, 2 + do cy = 1, 2 + local_corner = [lx(cx), ly(cy), 0._wp] + world_corner = matmul(patch_ib(patch_id)%rotation_matrix, local_corner) + center + bbox_min(1) = min(bbox_min(1), world_corner(1)) + bbox_min(2) = min(bbox_min(2), world_corner(2)) + bbox_max(1) = max(bbox_max(1), world_corner(1)) + bbox_max(2) = max(bbox_max(2), world_corner(2)) + end do + end do else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry - bounding_distance = 0._wp + bbox_min = center - 0.5_wp*max(patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y) + bbox_max = center + 0.5_wp*max(patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y) + else if (patch_ib(patch_id)%geometry == 9) then + ! cuboid geometries + bbox_min = center - 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) + bbox_max = center + 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) + else if (patch_ib(patch_id)%geometry == 10) then + ! cylinder geometry + bbox_min = center - sqrt(patch_ib(patch_id)%radius**2 + patch_ib(patch_id)%length_x**2) + bbox_max = center + sqrt(patch_ib(patch_id)%radius**2 + patch_ib(patch_id)%length_x**2) + else if (patch_ib(patch_id)%geometry == 12) then + ! Local-space bounding box extents (min=1, max=2 in the third index) + lx(1) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 1, 1) + patch_ib(patch_id)%centroid_offset(1) + lx(2) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 1, 3) + patch_ib(patch_id)%centroid_offset(1) + ly(1) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 2, 1) + patch_ib(patch_id)%centroid_offset(2) + ly(2) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 2, 3) + patch_ib(patch_id)%centroid_offset(2) + lz(1) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 3, 1) + patch_ib(patch_id)%centroid_offset(3) + lz(2) = stl_bounding_boxes(patch_ib(patch_id)%model_id, 3, 3) + patch_ib(patch_id)%centroid_offset(3) + + bbox_min = 1e12 + bbox_max = -1e12 + ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB + do cx = 1, 2 + do cy = 1, 2 + do cz = 1, 2 + local_corner = [lx(cx), ly(cy), lz(cz)] + world_corner = matmul(patch_ib(patch_id)%rotation_matrix, local_corner) + center + bbox_min(1) = min(bbox_min(1), world_corner(1)) + bbox_min(2) = min(bbox_min(2), world_corner(2)) + bbox_min(3) = min(bbox_min(3), world_corner(3)) + bbox_max(1) = max(bbox_max(1), world_corner(1)) + bbox_max(2) = max(bbox_max(2), world_corner(2)) + bbox_max(3) = max(bbox_max(3), world_corner(3)) + end do + end do + end do end if - end subroutine s_get_bounding_box_corner_distance + il = -gp_layers - 1 + jl = -gp_layers - 1 + kl = -gp_layers - 1 + ir = m + gp_layers + 1 + jr = n + gp_layers + 1 + kr = p + gp_layers + 1 + call get_indices_from_bounds(bbox_min(1), bbox_max(1), x_cc, il, ir) + call get_indices_from_bounds(bbox_min(2), bbox_max(2), y_cc, jl, jr) + if (num_dims == 3) call get_indices_from_bounds(bbox_min(3), bbox_max(3), z_cc, kl, kr) + + end subroutine get_bounding_indices - subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index) + subroutine get_indices_from_bounds(left_bound, right_bound, cell_centers, left_index, right_index) real(wp), intent(in) :: left_bound, right_bound integer, intent(inout) :: left_index, right_index @@ -700,7 +682,7 @@ contains end do right_index = itr_right - end subroutine get_bounding_indices + end subroutine get_indices_from_bounds !> Encode the patch ID with a unique offset containing periodicity information subroutine s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id) From f5dc0d7dd9449dd130fd2017cad87dd2dc73c4c7 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 11:25:53 -0400 Subject: [PATCH 07/14] All tests pass on GNU compilers --- src/simulation/m_ib_patches.fpp | 153 ++++++-------------------------- 1 file changed, 27 insertions(+), 126 deletions(-) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 30a16c12ff..dbc2fd92ee 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -46,7 +46,7 @@ contains integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds real(wp), dimension(3) :: center, xyz_local, length - real(wp) :: radius + real(wp) :: radius, eta ! 3D Patch Geometries @@ -66,8 +66,8 @@ contains ! find the indices to the left and right of the IB in i, j, k call get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) - $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, length, radius, airfoil_id]', copyin='[patch_id, & - & encoded_patch_id, center]', collapse=3) + $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, length, radius, airfoil_id, eta]', & + & copyin='[patch_id, encoded_patch_id, center]', collapse=3) do k = kl, kr do j = jl, jr do i = il, ir @@ -103,8 +103,11 @@ contains else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & - ! & k) = encoded_patch_id + model_id = patch_ib(patch_id)%model_id + eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xyz_local) + if (eta > stl_models(model_id)%model_threshold) then + ib_markers%sf(i, j, k) = encoded_patch_id + end if end if end do end do @@ -131,7 +134,7 @@ contains ! find the indices to the left and right of the IB in i, j, k call get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) - $:GPU_PARALLEL_LOOP(private='[i, j, xyz_local, airfoil_id]', copyin='[patch_id, encoded_patch_id, & + $:GPU_PARALLEL_LOOP(private='[i, j, xyz_local, airfoil_id, eta]', copyin='[patch_id, encoded_patch_id, & & center]', collapse=2) do j = jl, jr do i = il, ir @@ -160,8 +163,11 @@ contains else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & - ! & 0) = encoded_patch_id + model_id = patch_ib(patch_id)%model_id + eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xyz_local) + if (eta > stl_models(model_id)%model_threshold) then + ib_markers%sf(i, j, k) = encoded_patch_id + end if else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] @@ -184,7 +190,7 @@ contains integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds real(wp), dimension(3) :: center, xyz_local, length - real(wp) :: radius + real(wp) :: radius, eta if (num_dims == 3) then ! get the periodicities @@ -194,7 +200,7 @@ contains do yp = yp_lower, yp_upper do zp = zp_lower, zp_upper $:GPU_PARALLEL_LOOP(private='[xp, yp, zp, i, il, ir, j, jl, jr, k, kl, kr, xyz_local, length, radius, & - & patch_id, airfoil_id, model_id, encoded_patch_id, center]') + & patch_id, airfoil_id, model_id, encoded_patch_id, center, eta]') do patch_id = 1, num_ibs 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) @@ -241,8 +247,11 @@ contains else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), z_cc(k))) ib_markers%sf(i, j, & - ! & k) = encoded_patch_id + model_id = patch_ib(patch_id)%model_id + eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xyz_local) + if (eta > stl_models(model_id)%model_threshold) then + ib_markers%sf(i, j, k) = encoded_patch_id + end if end if end do end do @@ -260,7 +269,7 @@ contains do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper $:GPU_PARALLEL_LOOP(private='[xp, yp, i, il, ir, j, jl, jr, xyz_local, length, radius, & - & bounding_box_corner_distance, patch_id, airfoil_id, model_id, encoded_patch_id, center]') + & bounding_box_corner_distance, patch_id, airfoil_id, model_id, encoded_patch_id, center, eta]') do patch_id = 1, num_ibs 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) @@ -299,10 +308,12 @@ contains & airfoil_id)) ib_markers%sf(i, j, k) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry - model_id = patch_ib(patch_id)%model_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - ! if (f_is_inside_model(patch_id, x_cc(i), y_cc(j), 0._wp)) ib_markers%sf(i, j, & - ! & k) = encoded_patch_id + model_id = patch_ib(patch_id)%model_id + eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xyz_local) + if (eta > stl_models(model_id)%model_threshold) then + ib_markers%sf(i, j, k) = encoded_patch_id + end if else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] @@ -319,116 +330,6 @@ contains end subroutine s_apply_ib_patches_ib_parallelism - !> The STL patch is a 2D geometry that is imported from an STL file. - subroutine s_ib_model(patch_id, ib_markers, xp, yp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - integer, intent(in) :: xp, yp !< integers containing the periodicity projection information - integer :: i, j, il, ir, jl, jr !< Generic loop iterators - integer :: model_id, encoded_patch_id - integer :: cx, cy - real(wp) :: lx(2), ly(2) - real(wp), dimension(1:2) :: bbox_min, bbox_max - real(wp), dimension(1:3) :: local_corner, world_corner - real(wp) :: eta, threshold - real(wp), dimension(1:3) :: offset - real(wp), dimension(1:3) :: center, xy_local - real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation - - center = 0._wp - 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) - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) - rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:) - offset(:) = patch_ib(patch_id)%centroid_offset(:) - model_id = patch_ib(patch_id)%model_id - threshold = stl_models(model_id)%model_threshold - - ! 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) - - il = -gp_layers - 1 - jl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - $:GPU_PARALLEL_LOOP(private='[i, j, xy_local, eta]', copyin='[patch_id, model_id, encoded_patch_id, center, & - & inverse_rotation, offset, threshold]', collapse=2) - do i = il, ir - do j = jl, jr - xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] - xy_local = matmul(inverse_rotation, xy_local) - xy_local = xy_local - offset - - eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xy_local) - - ! Reading STL boundary vertices and compute the levelset and levelset_norm - if (eta > threshold) then - ib_markers%sf(i, j, 0) = encoded_patch_id - end if - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_model - - !> The STL patch is a 3D geometry that is imported from an STL file. - subroutine s_ib_3d_model(patch_id, ib_markers, xp, yp, zp) - - integer, intent(in) :: patch_id - type(integer_field), intent(inout) :: ib_markers - 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 :: model_id, encoded_patch_id - real(wp) :: eta, threshold - real(wp), dimension(1:3) :: offset - real(wp), dimension(1:3) :: center, xyz_local - real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation - integer :: cx, cy, cz - real(wp) :: lx(2), ly(2), lz(2) - real(wp), dimension(1:3) :: bbox_min, bbox_max, local_corner, world_corner - - center = 0._wp - 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) - inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:) - offset(:) = patch_ib(patch_id)%centroid_offset(:) - model_id = patch_ib(patch_id)%model_id - threshold = stl_models(model_id)%model_threshold - rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:) - - ! 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) - - il = -gp_layers - 1 - jl = -gp_layers - 1 - kl = -gp_layers - 1 - ir = m + gp_layers + 1 - jr = n + gp_layers + 1 - kr = p + gp_layers + 1 - - $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local, eta]', copyin='[patch_id, model_id, encoded_patch_id, center, & - & inverse_rotation, offset, threshold]', collapse=3) - do i = il, ir - do j = jl, jr - do k = kl, kr - xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)] - xyz_local = matmul(inverse_rotation, xyz_local) - xyz_local = xyz_local - offset - - eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xyz_local) - - if (eta > threshold) then - ib_markers%sf(i, j, k) = encoded_patch_id - end if - end do - end do - end do - $:END_GPU_PARALLEL_LOOP() - - end subroutine s_ib_3d_model - !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid) !! and before s_apply_ib_patches or s_apply_levelset. subroutine s_initialize_ib_airfoils() From 3cd4661d7561e0a7d6f29f3acc6ad6ad6952f27a Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 11:36:21 -0400 Subject: [PATCH 08/14] Added patch geometries to docs and fixed spelling --- docs/documentation/case.md | 3 ++- docs/module_categories.json | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index dd5f2692de..4fc1b0f98d 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -314,7 +314,8 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi | `num_ibs` | Integer | Number of immersed boundary patches | | `num_stl_models` | Integer | Number of STL/OBJ model entries in the `stl_models` array | | `num_particle_clouds` | Integer | Number of particle bed specifications to generate immersed boundary patches from | -| `ib_neighborhood_radius` | Integer | Parameter that controls the neighborhood size for IB detection. | +| `ib_neighborhood_radius` | Integer | Parameter that controls the neighborhood size for IB detection. | +| `many_ib_patch_parallelism` | Logical | Parallelize over IB patches instead of grid cells (better for many small patches). | | `geometry` | Integer | Geometry configuration of the patch.| | `x[y,z]_centroid` | Real | Centroid of the applied geometry in the [x,y,z]-direction. | | `length_x[y,z]` | Real | Length, if applicable, in the [x,y,z]-direction. | diff --git a/docs/module_categories.json b/docs/module_categories.json index 7cf0ccc995..fe38e94b3c 100644 --- a/docs/module_categories.json +++ b/docs/module_categories.json @@ -72,7 +72,8 @@ "m_checker", "m_checker_common", "m_sim_helpers", - "m_derived_variables" + "m_derived_variables", + "m_patch_geometries" ] }, { From 9b0b59d07933061442742a29dc924540c12429d1 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 11:38:28 -0400 Subject: [PATCH 09/14] Missed some spelling in my git add --- src/common/m_patch_geometries.fpp | 8 ++++---- toolchain/mfc/params/descriptions.py | 1 + 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/common/m_patch_geometries.fpp b/src/common/m_patch_geometries.fpp index 2709e89a6b..1efef96bbb 100644 --- a/src/common/m_patch_geometries.fpp +++ b/src/common/m_patch_geometries.fpp @@ -90,13 +90,13 @@ contains ! our check branches for the upper and lower half of the airfoil if (y >= 0._wp) then - ! incriment the iterator so we know where in the airfoil arrays to look + ! increment the iterator so we know where in the airfoil arrays to look k = 1 do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < x) k = k + 1 end do - ! If the values are approximately equivilent, skip the next check + ! If the values are approximately equivalent, skip the next check if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, x)) then if (y <= ib_airfoil_grids(airfoil_id)%upper(k)%y) is_inside = .true. else @@ -107,13 +107,13 @@ contains & is_inside = .true. end if else - ! incriment the iterator so we know where in the airfoil arrays to look + ! increment the iterator so we know where in the airfoil arrays to look k = 1 do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < x) k = k + 1 end do - ! If the values are approximately equivilent, skip the next check + ! If the values are approximately equivalent, skip the next check if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, x)) then if (y >= ib_airfoil_grids(airfoil_id)%lower(k)%y) is_inside = .true. else diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index c165577854..48b6f7712b 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -136,6 +136,7 @@ "num_stl_models": "Number of STL/OBJ model entries in the stl_models array", "num_particle_clouds": "Number of particle bed specifications to generate immersed boundary patches from", "ib_neighborhood_radius": "Neighborhood radius in ranks for IB awareness", + "many_ib_patch_parallelism": "Parallelize over IB patches instead of grid cells (better for many small patches)", # Acoustic sources "acoustic_source": "Enable acoustic source terms", "num_source": "Number of acoustic sources", From aeaa21e7a2e173e06a9cf7ed4d2332903f2488ed Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 12:18:10 -0400 Subject: [PATCH 10/14] Integrated with icpp patches and update documents to note the deprecation of length_y and length_z for cylinders. All tests pass on GNU compilers for full test suite. --- docs/documentation/case.md | 2 +- src/pre_process/m_icpp_patches.fpp | 37 +++++++++++++++--------------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 4fc1b0f98d..4f1687556e 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -1198,7 +1198,7 @@ Boundary is at polar angle \f$\theta = \mathrm{atan2}(y - y_{\mathrm{centroid}}, | 3 | 2D Rectangle | 2 | | 4 | 2D Airfoil | 2 | | 8 | 3D Sphere | 3 | -| 10 | 3D Cylinder | 3 | +| 10 | 3D Cylinder | 3 | `length_x` sets the axial length of the cylinder. | | 11 | 3D Airfoil | 3 | ### Acoustic Supports {#acoustic-supports} diff --git a/src/pre_process/m_icpp_patches.fpp b/src/pre_process/m_icpp_patches.fpp index 42e9332c53..f07f906a24 100644 --- a/src/pre_process/m_icpp_patches.fpp +++ b/src/pre_process/m_icpp_patches.fpp @@ -12,6 +12,7 @@ !> @brief Constructs initial condition patch geometries (lines, circles, rectangles, spheres, etc.) on the grid module m_icpp_patches + use m_patch_geometries use m_model ! Subroutine(s) related to STL files use m_derived_types ! Definitions of the derived types use m_global_parameters @@ -316,8 +317,8 @@ contains & dy)*(sqrt((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2) - radius))*(-0.5_wp) + 0.5_wp end if - if (((x_cc(i) - x_centroid)**2 + (y_cc(j) - y_centroid)**2 <= radius**2 & - & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, & + if ((f_is_inside_cylinder(x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp, radius, & + & 0._wp) .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, & & 0) == smooth_patch_id) then call s_assign_patch_primitive_variables(patch_id, i, j, 0, eta, q_prim_vf, patch_id_fp) @@ -490,8 +491,8 @@ contains & + 0.5_wp end if - if ((((x_cc(i) - x_centroid)/a)**2 + ((y_cc(j) - y_centroid)/b)**2 <= 1._wp & - & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, & + if ((f_is_inside_ellipse(x_cc(i) - x_centroid, y_cc(j) - y_centroid, [2._wp*a, 2._wp*b, & + & 0._wp]) .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) .or. patch_id_fp(i, j, & & 0) == smooth_patch_id) then call s_assign_patch_primitive_variables(patch_id, i, j, 0, eta, q_prim_vf, patch_id_fp) @@ -622,8 +623,7 @@ contains ! Assign patch vars if cell is covered and patch has write permission do j = 0, n do i = 0, m - if (x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i) .and. y_boundary%beg <= y_cc(j) & - & .and. y_boundary%end >= y_cc(j)) then + if (f_is_inside_cuboid(x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp, [length_x, length_y, 0._wp])) then if (patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) then call s_assign_patch_primitive_variables(patch_id, i, j, 0, eta, q_prim_vf, patch_id_fp) @@ -753,8 +753,8 @@ contains ! Assign patch vars if cell is covered and patch has write permission do j = 0, n do i = 0, m - if (x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i) .and. y_boundary%beg <= y_cc(j) & - & .and. y_boundary%end >= y_cc(j) .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) then + if (f_is_inside_cuboid(x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp, [length_x, length_y, & + & 0._wp]) .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) then call s_assign_patch_primitive_variables(patch_id, i, j, 0, eta, q_prim_vf, patch_id_fp) @:analytical() @@ -1003,8 +1003,8 @@ contains & - radius))*(-0.5_wp) + 0.5_wp end if - if ((((x_cc(i) - x_centroid)**2 + (cart_y - y_centroid)**2 + (cart_z - z_centroid)**2 <= radius**2) & - & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, & + if ((f_is_inside_sphere(x_cc(i) - x_centroid, cart_y - y_centroid, cart_z - z_centroid, & + & radius) .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, & & k) == smooth_patch_id) then call s_assign_patch_primitive_variables(patch_id, i, j, k, eta, q_prim_vf, patch_id_fp) @@ -1069,8 +1069,8 @@ contains cart_z = z_cc(k) end if - if (x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i) .and. y_boundary%beg <= cart_y & - & .and. y_boundary%end >= cart_y .and. z_boundary%beg <= cart_z .and. z_boundary%end >= cart_z) then + if (f_is_inside_cuboid(x_cc(i) - x_centroid, cart_y - y_centroid, cart_z - z_centroid, [length_x, length_y, & + & length_z])) then if (patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) then call s_assign_patch_primitive_variables(patch_id, i, j, k, eta, q_prim_vf, patch_id_fp) @@ -1160,12 +1160,13 @@ contains end if end if - if (((.not. f_is_default(length_x) .and. (cart_y - y_centroid)**2 + (cart_z - z_centroid)**2 <= radius**2 & - & .and. x_boundary%beg <= x_cc(i) .and. x_boundary%end >= x_cc(i)) .or. (.not. f_is_default(length_y) & - & .and. (x_cc(i) - x_centroid)**2 + (cart_z - z_centroid)**2 <= radius**2 .and. y_boundary%beg <= cart_y & - & .and. y_boundary%end >= cart_y) .or. (.not. f_is_default(length_z) .and. (x_cc(i) - x_centroid)**2 & - & + (cart_y - y_centroid)**2 <= radius**2 .and. z_boundary%beg <= cart_z .and. z_boundary%end >= cart_z) & - & .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, & + if (((.not. f_is_default(length_x) .and. f_is_inside_cylinder(cart_y - y_centroid, cart_z - z_centroid, & + & x_cc(i) - x_centroid, radius, & + & length_x)) .or. (.not. f_is_default(length_y) .and. f_is_inside_cylinder(x_cc(i) - x_centroid, & + & cart_z - z_centroid, cart_y - y_centroid, radius, & + & length_y)) .or. (.not. f_is_default(length_z) .and. f_is_inside_cylinder(x_cc(i) - x_centroid, & + & cart_y - y_centroid, cart_z - z_centroid, radius, & + & length_z)) .and. patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, k))) .or. patch_id_fp(i, j, & & k) == smooth_patch_id) then call s_assign_patch_primitive_variables(patch_id, i, j, k, eta, q_prim_vf, patch_id_fp) From 8c2ffb229b5e313194b2d2ed92264270778f847f Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 14:28:16 -0400 Subject: [PATCH 11/14] Need to hard-code index 0 in ib markers array in 2D to prevent out-of-bounds read --- src/simulation/m_ib_patches.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 2057221588..64f9018325 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -148,7 +148,7 @@ contains ! circular geometries radius = patch_ib(patch_id)%radius if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, & - & j, k) = encoded_patch_id + & j, 0) = encoded_patch_id else if (patch_ib(patch_id)%geometry == 3) then ! rectangular geometries length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] @@ -166,7 +166,7 @@ contains model_id = patch_ib(patch_id)%model_id eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local) if (eta > stl_models(model_id)%model_threshold) then - ib_markers%sf(i, j, k) = encoded_patch_id + ib_markers%sf(i, j, 0) = encoded_patch_id end if else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry From a2a3e8e027d72dc3f3e94f26d09f2e047a8a45a4 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 15:34:55 -0400 Subject: [PATCH 12/14] Updated a GPU paralellism macro that contained a removed variable --- src/simulation/m_ib_patches.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 64f9018325..59cc85f65c 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -269,7 +269,7 @@ contains do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper $:GPU_PARALLEL_LOOP(private='[xp, yp, i, il, ir, j, jl, jr, xyz_local, length, radius, & - & bounding_box_corner_distance, patch_id, airfoil_id, model_id, encoded_patch_id, center, eta]') + & patch_id, airfoil_id, model_id, encoded_patch_id, center, eta]') do patch_id = 1, num_ibs 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) From f168b40afe542a51bfd37d5a8ced99392790ad02 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 10 Jun 2026 15:37:36 -0400 Subject: [PATCH 13/14] Formatting --- src/simulation/m_ib_patches.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 59cc85f65c..aece4bbbe4 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -268,8 +268,8 @@ contains $:GPU_PARALLEL_LOOP(private='[xp, yp, patch_id, center]', collapse=3) do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper - $:GPU_PARALLEL_LOOP(private='[xp, yp, i, il, ir, j, jl, jr, xyz_local, length, radius, & - & patch_id, airfoil_id, model_id, encoded_patch_id, center, eta]') + $:GPU_PARALLEL_LOOP(private='[xp, yp, i, il, ir, j, jl, jr, xyz_local, length, radius, patch_id, airfoil_id, & + & model_id, encoded_patch_id, center, eta]') do patch_id = 1, num_ibs 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) From dac7c475d2cd8bbaf5f20e19b006f2db3059aece Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Wed, 10 Jun 2026 23:03:06 -0400 Subject: [PATCH 14/14] Passes test suite on NVHPC GPU OpenACC Parallelism --- src/common/m_model.fpp | 7 +++++-- src/simulation/m_global_parameters.fpp | 2 +- src/simulation/m_ib_patches.fpp | 5 ++++- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/common/m_model.fpp b/src/common/m_model.fpp index 5341a82ffa..76770a4a62 100644 --- a/src/common/m_model.fpp +++ b/src/common/m_model.fpp @@ -35,7 +35,7 @@ module m_model real(wp), allocatable, dimension(:,:,:,:) :: gpu_boundary_v integer, allocatable :: gpu_boundary_edge_count(:) real(wp), allocatable :: stl_bounding_boxes(:,:,:) - $:GPU_DECLARE(create='[gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count]') + $:GPU_DECLARE(create='[gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_boundary_edge_count, stl_bounding_boxes]') contains @@ -877,7 +877,7 @@ contains #endif if (num_stl_models == 0) return - allocate (stl_bounding_boxes(num_stl_models,1:3,1:3)) + @:ALLOCATE(stl_bounding_boxes(num_stl_models,1:3,1:3)) @:ALLOCATE(models(num_stl_models)) do stl_id = 1, num_stl_models @@ -950,6 +950,9 @@ contains end if end do + $:GPU_UPDATE(device='[stl_bounding_boxes]') + $:GPU_UPDATE(device='[stl_models(1:num_stl_models)]') + ! Pack and upload flat arrays for GPU (AFTER the loop) block integer :: mid, max_ntrs diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 550dd79420..ec92380f8e 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -270,7 +270,7 @@ module m_global_parameters type(ib_airfoil_grid), dimension(num_ib_airfoils_max) :: ib_airfoil_grids !< Per-airfoil computed surface grids type(ib_stl_parameters), dimension(num_stl_models_max) :: stl_models !< Per-STL model parameters (namelist) - $:GPU_DECLARE(create='[ib, num_ibs, patch_ib, ib_airfoil, ib_airfoil_grids]') + $:GPU_DECLARE(create='[ib, num_ibs, patch_ib, ib_airfoil, ib_airfoil_grids, stl_models]') $:GPU_DECLARE(create='[ib_coefficient_of_friction]') !> @} diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index aece4bbbe4..f8bb9a176b 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -265,7 +265,6 @@ contains ! get the periodicities call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper) - $:GPU_PARALLEL_LOOP(private='[xp, yp, patch_id, center]', collapse=3) do xp = xp_lower, xp_upper do yp = yp_lower, yp_upper $:GPU_PARALLEL_LOOP(private='[xp, yp, i, il, ir, j, jl, jr, xyz_local, length, radius, patch_id, airfoil_id, & @@ -548,6 +547,8 @@ contains subroutine get_indices_from_bounds(left_bound, right_bound, cell_centers, left_index, right_index) + $:GPU_ROUTINE(parallelism='[seq]') + real(wp), intent(in) :: left_bound, right_bound integer, intent(inout) :: left_index, right_index real(wp), dimension(-buff_size:), intent(in) :: cell_centers @@ -588,6 +589,8 @@ contains !> Encode the patch ID with a unique offset containing periodicity information subroutine s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id) + $:GPU_ROUTINE(parallelism='[seq]') + integer, intent(in) :: patch_id, x_periodicity, y_periodicity, z_periodicity integer, intent(out) :: encoded_patch_id integer :: temp_x_per, temp_y_per, temp_z_per, offset