diff --git a/docs/documentation/case.md b/docs/documentation/case.md index b62b3adcbd..1e62ecf520 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -320,7 +320,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. | @@ -1207,7 +1208,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/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" ] }, { diff --git a/examples/2D_mibm_particle_cloud/case.py b/examples/2D_mibm_particle_cloud/case.py index e8214d7d82..97c14b81d0 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..1efef96bbb --- /dev/null +++ b/src/common/m_patch_geometries.fpp @@ -0,0 +1,143 @@ +!> +!! @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_sphere, f_is_inside_cylinder, f_is_inside_cuboid, f_is_inside_airfoil, f_is_inside_ellipse + +contains + + !> 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(x, y, z, radius) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: radius, x, y, z + logical :: is_inside + + 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(polar_x, polar_y, height, radius, length) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + 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 <= 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 + + 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(x, y, z, length) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + 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*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*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(x, y, z, length, airfoil_id) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: x, y, z, length + integer, intent(in) :: airfoil_id + logical :: is_inside + integer :: k + real(wp) :: f + + 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*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 + ! 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 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 + ! 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 + ! 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 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 + ! 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(x, y, length) result(is_inside) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: x, y + real(wp), dimension(3), intent(in) :: length + logical :: is_inside + + ! 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 + + 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 82869f7709..e11fbcbfd4 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -101,12 +101,16 @@ 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) + !> 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 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 aa6e5a482b..2c117967d4 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -96,7 +96,9 @@ module m_global_parameters !> @{ 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) + !> 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/pre_process/m_icpp_patches.fpp b/src/pre_process/m_icpp_patches.fpp index f526c42672..9775d15bc7 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 @@ -323,8 +324,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) @@ -497,8 +498,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) @@ -629,8 +630,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) @@ -760,8 +760,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() @@ -1010,8 +1010,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) @@ -1076,8 +1076,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) @@ -1167,12 +1167,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) 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 77eb90a31d..550dd79420 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -535,6 +535,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_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 3f01f4df8e..aece4bbbe4 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 @@ -26,6 +27,309 @@ module m_ib_patches contains + !> 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 + + 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, 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, eta + + ! 3D Patch Geometries + + if (num_dims == 3) then + 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 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) + + ! 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 + 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, eta]', & + & 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(k) - 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(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(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(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(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 + 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 + end if + end if + end do + end do + end do + $:END_GPU_PARALLEL_LOOP() + end do + end do + end do + end do + + ! 2D Patch Geometries + else if (num_dims == 2) then + 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 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 + + ! 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 + call get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) + + $: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 + ! 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) + + ! perform the interior check for the patch geometry of this IB + if (patch_ib(patch_id)%geometry == 2) then + ! 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, 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] + 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(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 + 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, 0) = 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] + 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 + $:END_GPU_PARALLEL_LOOP() + end do + end do + end do + end if + + end subroutine s_apply_ib_patches_grid_cell_parallelism + + 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) :: radius, eta + + 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, & + & 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) + 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 + 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 + ! get coordinate frame centered on IB + 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) + + ! 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(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(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(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(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 + 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 + end if + 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 + 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, & + & 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) + center(3) = 0._wp + + ! 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 + ! 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 + ! 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) + + ! perform the interior check for the patch geometry of this IB + if (patch_ib(patch_id)%geometry == 2) then + ! 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 + 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(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(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 + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset + 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 + 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] + 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 + end do + $:END_GPU_PARALLEL_LOOP() + end do + end do + end if + + end subroutine s_apply_ib_patches_ib_parallelism + !> 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() @@ -98,730 +402,6 @@ contains 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 - integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds - - ! 3D Patch Geometries - - if (p > 0) 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 - end do - end do - end do - end do - !> @} - - ! 2D Patch Geometries - else if (n > 0) 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 - 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 - 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 - - 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() - - end subroutine s_ib_ellipse - - !> 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 - - ! 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 - 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(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 - - ! 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 - 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(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 - !> Compute a rotation matrix for converting to the rotating frame of the boundary subroutine s_update_ib_rotation_matrix(patch_id) @@ -869,7 +449,104 @@ contains end subroutine s_update_ib_rotation_matrix - subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index) + subroutine get_bounding_indices(center, patch_id, il, ir, jl, jr, kl, kr) + + $:GPU_ROUTINE(parallelism='[seq]') + + 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 .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 + 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 + 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 + 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 + + 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_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 @@ -906,7 +583,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) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 99aeb1e382..23f70b140e 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 48ce9d74e7..3d283aa822 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -643,6 +643,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"]: @@ -1276,6 +1277,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", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index f9e7e13863..ceab47a585 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",