Moving Levelset Calculation to Pre-Process

src/common/m_derived_types.fpp

@@ -34,6 +34,16 @@ module m_derived_types
         integer, pointer, dimension(:, :, :) :: sf => null()
     end type integer_field
 
+    !> Derived type for levelset
+    type levelset_field
+        real(kind(0d0)), pointer, dimension(:, :, :, :) :: sf => null()
+    end type levelset_field
+
+    !> Derived type for levelset norm
+    type levelset_norm_field
+        real(kind(0d0)), pointer, dimension(:, :, :, :, :) :: sf => null()
+    end type levelset_norm_field
+
     type mpi_io_var
         integer, allocatable, dimension(:) :: view
         type(scalar_field), allocatable, dimension(:) :: var
@@ -44,6 +54,16 @@ module m_derived_types
         type(integer_field) :: var
     end type mpi_io_ib_var
 
+    type mpi_io_levelset_var
+        integer :: view
+        type(levelset_field) :: var
+    end type mpi_io_levelset_var
+
+    type mpi_io_levelset_norm_var
+        integer :: view
+        type(levelset_norm_field) :: var
+    end type mpi_io_levelset_norm_var
+
     !> Derived type annexing a vector field (VF)
     type vector_field
         type(scalar_field), allocatable, dimension(:) :: vf !< Vector field

src/common/m_mpi_common.fpp

@@ -62,7 +62,9 @@
 
     !! @param q_cons_vf Conservative variables
     !! @param ib_markers track if a cell is within the immersed boundary
-    subroutine s_initialize_mpi_data(q_cons_vf, ib_markers)
+    !! @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
+    subroutine s_initialize_mpi_data(q_cons_vf, ib_markers, levelset, levelset_norm)
 
         type(scalar_field), &
             dimension(sys_size), &
@@ -72,6 +74,14 @@
             optional, &
             intent(in) :: ib_markers
 
+        type(levelset_field), &
+            optional, &
+            intent(IN) :: levelset
+
+        type(levelset_norm_field), &
+            optional, &
+            intent(IN) :: levelset_norm
+
         integer, dimension(num_dims) :: sizes_glb, sizes_loc
         integer, dimension(1) :: airfoil_glb, airfoil_loc, airfoil_start
 
@@ -137,13 +147,30 @@
 
 #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_DOUBLE_PRECISION, MPI_IO_levelset_DATA%view, ierr)
+            call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
+                                          MPI_ORDER_FORTRAN, MPI_DOUBLE_PRECISION, 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

src/post_process/m_global_parameters.fpp

@@ -155,10 +155,14 @@ module m_global_parameters
     integer, allocatable, dimension(:) :: start_idx !<
     !! Starting cell-center index of local processor in global grid
 
+    integer :: num_ibs  !< Number of immersed boundaries
+
 #ifdef MFC_MPI
 
     type(mpi_io_var), public :: MPI_IO_DATA
     type(mpi_io_ib_var), public :: MPI_IO_IB_DATA
+    type(mpi_io_levelset_var), public :: MPI_IO_levelset_DATA
+    type(mpi_io_levelset_norm_var), public :: MPI_IO_levelsetnorm_DATA
 
 #endif
 
@@ -401,6 +405,9 @@ contains
         surface_tension = .false.
         adv_n = .false.
 
+        ! IBM
+        num_ibs = dflt_int
+
     end subroutine s_assign_default_values_to_user_inputs
 
     !>  Computation of parameters, allocation procedures, and/or

src/post_process/m_mpi_proxy.fpp

@@ -159,7 +159,7 @@ contains
             & 'model_eqns', 'num_fluids', 'bc_x%beg', 'bc_x%end', 'bc_y%beg',  &
             & 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'flux_lim', 'format',        &
             & 'precision', 'fd_order', 'thermal', 'nb', 'relax_model',         &
-            & 'n_start' ]
+            & 'n_start', 'num_ibs' ]
             call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
         #:endfor
 

src/post_process/m_start_up.f90

@@ -80,7 +80,7 @@ subroutine s_read_input_file
             parallel_io, rhoref, pref, bubbles, qbmm, sigR, &
             R0ref, nb, polytropic, thermal, Ca, Web, Re_inv, &
             polydisperse, poly_sigma, file_per_process, relax, &
-            relax_model, cf_wrt, sigma, adv_n, ib, &
+            relax_model, cf_wrt, sigma, adv_n, ib, num_ibs, &
             cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, &
             cfl_target, surface_tension
 

src/simulation/m_compute_levelset.fpp → src/pre_process/m_compute_levelset.fpp

@@ -38,12 +38,11 @@
 
 contains
 
-    !>  Initialize IBM module
     subroutine s_compute_circle_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
-        integer, intent(in) :: ib_patch_id
+        type(levelset_field), intent(INOUT) :: levelset
+        type(levelset_norm_field), intent(INOUT) :: levelset_norm
+        integer, intent(IN) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist
         real(kind(0d0)) :: x_centroid, y_centroid
@@ -62,11 +61,11 @@
                 dist_vec(2) = y_cc(j) - y_centroid
                 dist_vec(3) = 0
                 dist = dsqrt(sum(dist_vec**2))
-                levelset(i, j, 0, ib_patch_id) = dist - radius
+                levelset%sf(i, j, 0, ib_patch_id) = dist - radius
                 if (dist == 0) then
-                    levelset_norm(i, j, 0, ib_patch_id, :) = 0
+                    levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
                 else
-                    levelset_norm(i, j, 0, ib_patch_id, :) = &
+                    levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
                         dist_vec(:)/dist
                 end if
 
@@ -77,9 +76,9 @@
 
     subroutine s_compute_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
-        integer, intent(in) :: ib_patch_id
+        type(levelset_field), intent(INOUT) :: levelset
+        type(levelset_norm_field), intent(INOUT) :: levelset_norm
+        integer, intent(IN) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist, global_dist
         integer :: global_id
@@ -145,11 +144,11 @@
                     dist = global_dist
                 end if
 
-                levelset(i, j, 0, ib_patch_id) = dist
+                levelset%sf(i, j, 0, ib_patch_id) = dist
                 if (dist == 0) then
-                    levelset_norm(i, j, 0, ib_patch_id, :) = 0
+                    levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
                 else
-                    levelset_norm(i, j, 0, ib_patch_id, :) = &
+                    levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
                         dist_vec(:)/dist
                 end if
 
@@ -160,9 +159,9 @@
 
     subroutine s_compute_3D_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
-        integer, intent(in) :: ib_patch_id
+        type(levelset_field), intent(INOUT) :: levelset
+        type(levelset_norm_field), intent(INOUT) :: levelset_norm
+        integer, intent(IN) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist, dist_surf, dist_side, global_dist
         integer :: global_id
@@ -237,18 +236,18 @@
                     dist_side = min(abs(z_cc(l) - z_min), abs(z_max - z_cc(l)))
 
                     if (dist_side < dist_surf) then
-                        levelset(i, j, l, ib_patch_id) = dist_side
+                        levelset%sf(i, j, l, ib_patch_id) = dist_side
                         if (dist_side == abs(z_cc(l) - z_min)) then
-                            levelset_norm(i, j, l, ib_patch_id, :) = (/0, 0, -1/)
+                            levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, -1/)
                         else
-                            levelset_norm(i, j, l, ib_patch_id, :) = (/0, 0, 1/)
+                            levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, 1/)
                         end if
                     else
-                        levelset(i, j, l, ib_patch_id) = dist_surf
+                        levelset%sf(i, j, l, ib_patch_id) = dist_surf
                         if (dist == 0) then
-                            levelset_norm(i, j, l, ib_patch_id, :) = 0
+                            levelset_norm%sf(i, j, l, ib_patch_id, :) = 0
                         else
-                            levelset_norm(i, j, l, ib_patch_id, :) = &
+                            levelset_norm%sf(i, j, l, ib_patch_id, :) = &
                                 dist_vec(:)/dist_surf
                         end if
                     end if
@@ -262,10 +261,10 @@
     !>  Initialize IBM module
     subroutine s_compute_rectangle_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
-        integer, intent(in) :: ib_patch_id
+        type(levelset_field), intent(INOUT) :: levelset
+        type(levelset_norm_field), intent(INOUT) :: levelset_norm
 
+        integer :: ib_patch_id
         real(kind(0d0)) :: top_right(2), bottom_left(2)
         real(kind(0d0)) :: x, y, min_dist
         real(kind(0d0)) :: side_dists(4)
@@ -303,37 +302,37 @@
                     end if
 
                     if (min_dist == abs(side_dists(1))) then
-                        levelset(i, j, 0, ib_patch_id) = side_dists(1)
+                        levelset%sf(i, j, 0, ib_patch_id) = side_dists(1)
                         if (side_dists(1) == 0) then
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
                         else
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(1)/ &
-                                                                     abs(side_dists(1))
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(1)/ &
+                                                                        abs(side_dists(1))
                         end if
 
                     else if (min_dist == abs(side_dists(2))) then
-                        levelset(i, j, 0, ib_patch_id) = side_dists(2)
+                        levelset%sf(i, j, 0, ib_patch_id) = side_dists(2)
                         if (side_dists(2) == 0) then
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
                         else
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(2)/ &
-                                                                     abs(side_dists(2))
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(2)/ &
+                                                                        abs(side_dists(2))
                         end if
 
                     else if (min_dist == abs(side_dists(3))) then
                         if (side_dists(3) == 0) then
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
                         else
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(3)/ &
-                                                                     abs(side_dists(3))
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(3)/ &
+                                                                        abs(side_dists(3))
                         end if
 
                     else if (min_dist == abs(side_dists(4))) then
                         if (side_dists(4) == 0) then
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = 0d0
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0d0
                         else
-                            levelset_norm(i, j, 0, ib_patch_id, 1) = side_dists(4)/ &
-                                                                     abs(side_dists(4))
+                            levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(4)/ &
+                                                                        abs(side_dists(4))
                         end if
 
                     end if
@@ -347,9 +346,9 @@
 
     subroutine s_compute_sphere_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
-        integer, intent(in) :: ib_patch_id
+        type(levelset_field), intent(INOUT) :: levelset
+        type(levelset_norm_field), intent(INOUT) :: levelset_norm
+        integer, intent(IN) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist
         real(kind(0d0)) :: x_centroid, y_centroid, z_centroid
@@ -369,11 +368,11 @@
                     dist_vec(2) = y_cc(j) - y_centroid
                     dist_vec(3) = z_cc(k) - z_centroid
                     dist = dsqrt(sum(dist_vec**2))
-                    levelset(i, j, k, ib_patch_id) = dist - radius
+                    levelset%sf(i, j, k, ib_patch_id) = dist - radius
                     if (dist == 0) then
-                        levelset_norm(i, j, k, ib_patch_id, :) = (/1, 0, 0/)
+                        levelset_norm%sf(i, j, k, ib_patch_id, :) = (/1, 0, 0/)
                     else
-                        levelset_norm(i, j, k, ib_patch_id, :) = &
+                        levelset_norm%sf(i, j, k, ib_patch_id, :) = &
                             dist_vec(:)/dist
                     end if
                 end do
@@ -384,9 +383,9 @@
 
     subroutine s_compute_cylinder_levelset(levelset, levelset_norm, ib_patch_id)
 
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs), intent(inout) :: levelset
-        real(kind(0d0)), dimension(0:m, 0:n, 0:p, num_ibs, 3), intent(inout) :: levelset_norm
-        integer, intent(in) :: ib_patch_id
+        type(levelset_field), intent(INOUT) :: levelset
+        type(levelset_norm_field), intent(INOUT) :: levelset_norm
+        integer, intent(IN) :: ib_patch_id
 
         real(kind(0d0)) :: radius, dist
         real(kind(0d0)) :: x_centroid, y_centroid, z_centroid
@@ -434,20 +433,20 @@
                                    - radius
 
                     if (dist_side < abs(dist_surface)) then
-                        levelset(i, j, k, ib_patch_id) = -dist_side
+                        levelset%sf(i, j, k, ib_patch_id) = -dist_side
                         if (dist_side == abs(side_pos - boundary%beg)) then
-                            levelset_norm(i, j, k, ib_patch_id, :) = -dist_sides_vec
+                            levelset_norm%sf(i, j, k, ib_patch_id, :) = -dist_sides_vec
                         else
-                            levelset_norm(i, j, k, ib_patch_id, :) = dist_sides_vec
+                            levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_sides_vec
                         end if
                     else
-                        levelset(i, j, k, ib_patch_id) = dist_surface
+                        levelset%sf(i, j, k, ib_patch_id) = dist_surface
 
-                        levelset_norm(i, j, k, ib_patch_id, :) = &
+                        levelset_norm%sf(i, j, k, ib_patch_id, :) = &
                             (pos_vec - centroid_vec)*dist_surface_vec
-                        levelset_norm(i, j, k, ib_patch_id, :) = &
-                            levelset_norm(i, j, k, ib_patch_id, :)/ &
-                            norm2(levelset_norm(i, j, k, ib_patch_id, :))
+                        levelset_norm%sf(i, j, k, ib_patch_id, :) = &
+                            levelset_norm%sf(i, j, k, ib_patch_id, :)/ &
+                            norm2(levelset_norm%sf(i, j, k, ib_patch_id, :))
                     end if
                 end do
             end do