1454-projection-optimization

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. |
@@ -1248,7 +1249,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}

docs/module_categories.json

@@ -81,7 +81,8 @@
             "m_checker",
             "m_checker_common",
             "m_sim_helpers",
-            "m_derived_variables"
+            "m_derived_variables",
+            "m_patch_geometries"
         ]
     },
     {

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,

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

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

src/post_process/m_global_parameters.fpp

@@ -97,6 +97,10 @@ module m_global_parameters
 
     ! shear_num, shear_indices, shear_BC_flip_num, shear_BC_flip_indices: in m_global_parameters_common
     ! proc_coords, start_idx, mpiiofs, mpi_info_int: in m_global_parameters_common
+    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

src/pre_process/m_global_parameters.fpp

@@ -70,7 +70,12 @@ module m_global_parameters
     integer                             :: nmom  !< Number of carried moments
     !> @}
 
+    !> @name Immersed Boundaries
+    !> @{
     ! patch_ib, ib_airfoil, stl_models: auto-generated in generated_decls.fpp
+    !> Per-airfoil computed surface grids (unused in pre_process)
+    type(ib_airfoil_grid), allocatable, dimension(:) :: ib_airfoil_grids
+    !> @}
 
     !> @name Non-polytropic bubble gas compression
     !> @{

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)
 

src/simulation/m_checker.fpp

@@ -38,6 +38,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

src/simulation/m_global_parameters.fpp

@@ -426,6 +426,7 @@ contains
         collision_time = dflt_real
         ib_coefficient_of_friction = dflt_real
         ib_state_wrt = .false.
+        many_ib_patch_parallelism = .false.
 
         ! Bubble modeling (sim-specific)
         bubble_model = 1