Levelset refactor

examples/2D_forward_facing_step/case.py

@@ -70,10 +70,10 @@
             "patch_icpp(1)%alpha(1)": 1.0,
             # Patch: No slip rectangle
             "patch_ib(1)%geometry": 3,
-            "patch_ib(1)%x_centroid": 11.5 * h,
-            "patch_ib(1)%y_centroid": 0 * h,
-            "patch_ib(1)%length_x": 17 * h,
-            "patch_ib(1)%length_y": 2 * h,
+            "patch_ib(1)%x_centroid": 9 * h,
+            "patch_ib(1)%y_centroid": 0.5 * h,
+            "patch_ib(1)%length_x": 12 * h,
+            "patch_ib(1)%length_y": h,
             "patch_ib(1)%slip": "T",
             # Fluids Physical Parameters
             "fluid_pp(1)%gamma": 1.0 / (gam_a - 1.0),

src/common/m_derived_types.fpp

@@ -180,8 +180,18 @@ module m_derived_types
     type :: t_model
         integer :: ntrs   ! Number of triangles
         type(t_triangle), allocatable :: trs(:) ! Triangles
+
     end type t_model
 
+    type :: t_model_array
+        type(t_model), allocatable :: model
+        real(wp), allocatable, dimension(:, :, :) :: boundary_v
+        real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v
+        integer :: boundary_edge_count
+        integer :: total_vertices
+        logical :: interpolate
+    end type t_model_array
+
     !> Derived type adding initial condition (ic) patch parameters as attributes
     !! NOTE: The requirements for the specification of the above parameters
     !! are strongly dependent on both the choice of the multicomponent flow
@@ -437,6 +447,8 @@ module m_derived_types
         integer, dimension(3) :: ip_grid !< Top left grid point of IP
         real(wp), dimension(2, 2, 2) :: interp_coeffs !< Interpolation Coefficients of image point
         integer :: ib_patch_id !< ID of the IB Patch the ghost point is part of
+        real(wp) :: levelset
+        real(wp), dimension(1:3) :: levelset_norm
         logical :: slip
         integer, dimension(3) :: DB
     end type ghost_point

src/common/m_model.fpp

@@ -487,6 +487,8 @@ contains
     !! @return True if the point is inside the octree, false otherwise.
     impure function f_model_is_inside(model, point, spacing, spc) result(fraction)
 
+        ! $:GPU_ROUTINE(parallelism='[seq]')
+
         type(t_model), intent(in) :: model
         real(wp), dimension(1:3), intent(in) :: point
         real(wp), dimension(1:3), intent(in) :: spacing
@@ -499,6 +501,15 @@ contains
 
         real(wp), dimension(1:spc, 1:3) :: ray_origins, ray_dirs
 
+        ! TODO :: The random number generation prohibits GPU compute due to the subroutine not being able to be called in kernels
+        ! This should be swapped out with something that allows GPU compute. I recommend the fibonacci sphere:
+        ! do i = 1, spc
+        !   phi = acos(1.0 - 2.0*(i-1.0)/(spc-1.0))
+        !   theta = pi * (1.0 + sqrt(5.0)) * (i-1.0)
+        !   ray_dirs(i,:) = [cos(theta)*sin(phi), sin(theta)*sin(phi), cos(phi)]
+        !   ray_origins(i,:) = point
+        ! end do
+
         do i = 1, spc
             call random_number(ray_origins(i, :))
             ray_origins(i, :) = point + (ray_origins(i, :) - 0.5_wp)*spacing(:)
@@ -739,7 +750,6 @@ contains
         interpolate = .false.
 
         do j = 1, boundary_edge_count
-
             l1 = sqrt((boundary_v(j, 2, 1) - boundary_v(j, 1, 1))**2 + &
                       (boundary_v(j, 2, 2) - boundary_v(j, 1, 2))**2)
 
@@ -1049,6 +1059,8 @@ contains
     !! @param distance     The output levelset distance
     subroutine f_distance_normals_3D(model, point, normals, distance)
 
+        $:GPU_ROUTINE(parallelism='[seq]')
+
         type(t_model), intent(IN) :: model
         real(wp), dimension(1:3), intent(in) :: point
         real(wp), dimension(1:3), intent(out) :: normals
@@ -1112,8 +1124,11 @@ contains
     !! @return                             Distance which the levelset distance without interpolation
     function f_distance(boundary_v, boundary_edge_count, point) result(distance)
 
+        $:GPU_ROUTINE(parallelism='[seq]')
+
         integer, intent(in) :: boundary_edge_count
-        real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v
+        real(wp), intent(in), dimension(:, :, :) :: boundary_v
+        ! real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v
         real(wp), dimension(1:3), intent(in) :: point
 
         integer :: i
@@ -1143,8 +1158,10 @@ contains
     !! @param normals                      Output levelset normals without interpolation
     subroutine f_normals(boundary_v, boundary_edge_count, point, normals)
 
+        $:GPU_ROUTINE(parallelism='[seq]')
+
         integer, intent(in) :: boundary_edge_count
-        real(wp), intent(in), dimension(1:boundary_edge_count, 1:3, 1:2) :: boundary_v
+        real(wp), intent(in), dimension(:, :, :) :: boundary_v
         real(wp), dimension(1:3), intent(in) :: point
         real(wp), dimension(1:3), intent(out) :: normals
 
@@ -1203,8 +1220,10 @@ contains
     !! @return                             Distance which the levelset distance without interpolation
     function f_interpolated_distance(interpolated_boundary_v, total_vertices, point) result(distance)
 
+        $:GPU_ROUTINE(parallelism='[seq]')
+
         integer, intent(in) :: total_vertices
-        real(wp), intent(in), dimension(1:total_vertices, 1:3) :: interpolated_boundary_v
+        real(wp), intent(in), dimension(:, :) :: interpolated_boundary_v
         real(wp), dimension(1:3), intent(in) :: point
 
         integer :: i !< Loop iterator

src/common/m_mpi_common.fpp

@@ -124,15 +124,11 @@ contains
 
     !! @param q_cons_vf Conservative variables
     !! @param ib_markers track if a cell is within the immersed boundary
-    !! @param levelset closest distance from every cell to the IB
-    !! @param levelset_norm normalized vector from every cell to the closest point to the IB
     !! @param beta Eulerian void fraction from lagrangian bubbles
-    impure subroutine s_initialize_mpi_data(q_cons_vf, ib_markers, levelset, levelset_norm, beta)
+    impure subroutine s_initialize_mpi_data(q_cons_vf, ib_markers, beta)
 
         type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
         type(integer_field), optional, intent(in) :: ib_markers
-        type(levelset_field), optional, intent(IN) :: levelset
-        type(levelset_norm_field), optional, intent(IN) :: levelset_norm
         type(scalar_field), intent(in), optional :: beta
 
         integer, dimension(num_dims) :: sizes_glb, sizes_loc
@@ -203,73 +199,13 @@ contains
         end if
 #endif
 
+#ifndef MFC_PRE_PROCESS
         if (present(ib_markers)) then
-
-#ifdef MFC_PRE_PROCESS
-            MPI_IO_IB_DATA%var%sf => ib_markers%sf
-            MPI_IO_levelset_DATA%var%sf => levelset%sf
-            MPI_IO_levelsetnorm_DATA%var%sf => levelset_norm%sf
-#else
             MPI_IO_IB_DATA%var%sf => ib_markers%sf(0:m, 0:n, 0:p)
 
-#ifndef MFC_POST_PROCESS
-            MPI_IO_levelset_DATA%var%sf => levelset%sf(0:m, 0:n, 0:p, 1:num_ibs)
-            MPI_IO_levelsetnorm_DATA%var%sf => levelset_norm%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3)
-#endif
-
-#endif
             call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
                                           MPI_ORDER_FORTRAN, MPI_INTEGER, MPI_IO_IB_DATA%view, ierr)
             call MPI_TYPE_COMMIT(MPI_IO_IB_DATA%view, ierr)
-
-#ifndef MFC_POST_PROCESS
-            call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
-                                          MPI_ORDER_FORTRAN, mpi_p, MPI_IO_levelset_DATA%view, ierr)
-            call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
-                                          MPI_ORDER_FORTRAN, mpi_p, MPI_IO_levelsetnorm_DATA%view, ierr)
-
-            call MPI_TYPE_COMMIT(MPI_IO_levelset_DATA%view, ierr)
-            call MPI_TYPE_COMMIT(MPI_IO_levelsetnorm_DATA%view, ierr)
-#endif
-        end if
-
-#ifndef MFC_POST_PROCESS
-        if (present(ib_markers)) then
-            do j = 1, num_ibs
-                if (patch_ib(j)%c > 0) then
-
-#ifdef MFC_PRE_PROCESS
-                    allocate (MPI_IO_airfoil_IB_DATA%var(1:2*Np))
-#endif
-
-                    airfoil_glb(1) = 3*Np*num_procs
-                    airfoil_loc(1) = 3*Np
-                    airfoil_start(1) = 3*proc_rank*Np
-
-#ifdef MFC_PRE_PROCESS
-                    do i = 1, Np
-                        MPI_IO_airfoil_IB_DATA%var(i)%x = airfoil_grid_l(i)%x
-                        MPI_IO_airfoil_IB_DATA%var(i)%y = airfoil_grid_l(i)%y
-                    end do
-#endif
-
-                    call MPI_TYPE_CREATE_SUBARRAY(1, airfoil_glb, airfoil_loc, airfoil_start, &
-                                                  MPI_ORDER_FORTRAN, mpi_p, MPI_IO_airfoil_IB_DATA%view(1), ierr)
-                    call MPI_TYPE_COMMIT(MPI_IO_airfoil_IB_DATA%view(1), ierr)
-
-#ifdef MFC_PRE_PROCESS
-                    do i = 1, Np
-                        MPI_IO_airfoil_IB_DATA%var(Np + i)%x = airfoil_grid_u(i)%x
-                        MPI_IO_airfoil_IB_DATA%var(Np + i)%y = airfoil_grid_u(i)%y
-                    end do
-#endif
-                    call MPI_TYPE_CREATE_SUBARRAY(1, airfoil_glb, airfoil_loc, airfoil_start, &
-                                                  MPI_ORDER_FORTRAN, mpi_p, MPI_IO_airfoil_IB_DATA%view(2), ierr)
-                    call MPI_TYPE_COMMIT(MPI_IO_airfoil_IB_DATA%view(2), ierr)
-
-                end if
-            end do
-
         end if
 #endif