refactor periodic centroid projection into a single subroutine

Author: conraddelgado

examples/2D_ibm_circle_periodic/case.py

@@ -35,7 +35,7 @@
 Nt = 5000
 t_save = 20
 
-# to fully resolve requires 40-60 cells across sphere diameter
+# to fully resolve requires 40-60 cells across circle diameter
 Nx = 511
 Ny = Nx
 Nz = 0

src/common/m_compute_levelset.fpp

@@ -24,7 +24,8 @@ module m_compute_levelset
  s_rectangle_levelset, &
  s_cuboid_levelset, &
  s_sphere_levelset, &
- s_ellipse_levelset
+ s_ellipse_levelset, &
+ s_periodic_project_center
 
 contains
 
@@ -47,21 +48,10 @@ contains
         center(1, 1) = patch_ib(ib_patch_id)%x_centroid
         center(2, 1) = patch_ib(ib_patch_id)%y_centroid
 
-        if (periodic_ibs) then ! periodically wrap spheres around domain
-            if ((center(1, 1) - domain_glb(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 2) + (center(1, 1) - domain_glb(1, 1))
-            else if ((domain_glb(1, 2) - center(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 1) - (domain_glb(1, 2) - center(1, 1))
-            else
-                center(1, 2) = center(1, 1)
-            end if
-            if ((center(2, 1) - domain_glb(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 2) + (center(2, 1) - domain_glb(2, 1))
-            else if ((domain_glb(2, 2) - center(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 1) - (domain_glb(2, 2) - center(2, 1))
-            else
-                center(2, 2) = center(2, 1)
-            end if
+        if (periodic_ibs) then ! periodically wrap circles around domain
+            do i = 1, 2
+                call s_periodic_project_center(center(i, :), domain_glb(i, :), radius)
+            end do
         end if
 
         $:GPU_PARALLEL_LOOP(private='[i,j,ix,iy,dist_vec,dist,dist_vec_temp,dist_temp]', &
@@ -71,22 +61,23 @@ contains
                 dist_vec(1) = x_cc(i) - center(1, 1)
                 dist_vec(2) = y_cc(j) - center(2, 1)
                 dist_vec(3) = 0._wp
-                dist = sqrt(sum(dist_vec**2))
+                dist = sum(dist_vec**2)
                 if (periodic_ibs) then
                     ! check all permutations of periodically projected ib to find minimum distance
                     do ix = 1, 2
                         do iy = 1, 2
                             dist_vec_temp(1) = x_cc(i) - center(1, ix)
                             dist_vec_temp(2) = y_cc(j) - center(2, iy)
                             dist_vec_temp(3) = 0._wp
-                            dist_temp = sqrt(sum(dist_vec_temp**2))
+                            dist_temp = sum(dist_vec_temp**2)
                             if (dist_temp < dist) then
                                 dist = dist_temp
                                 dist_vec = dist_vec_temp
                             end if
                         end do
                     end do
                 end if
+                dist = sqrt(dist)
                 levelset%sf(i, j, 0, ib_patch_id) = dist - radius
                 if (f_approx_equal(dist, 0._wp)) then
                     levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
@@ -558,27 +549,9 @@ contains
         center(3, 1) = patch_ib(ib_patch_id)%z_centroid
 
         if (periodic_ibs) then ! periodically wrap spheres around domain
-            if ((center(1, 1) - domain_glb(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 2) + (center(1, 1) - domain_glb(1, 1))
-            else if ((domain_glb(1, 2) - center(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 1) - (domain_glb(1, 2) - center(1, 1))
-            else
-                center(1, 2) = center(1, 1)
-            end if
-            if ((center(2, 1) - domain_glb(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 2) + (center(2, 1) - domain_glb(2, 1))
-            else if ((domain_glb(2, 2) - center(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 1) - (domain_glb(2, 2) - center(2, 1))
-            else
-                center(2, 2) = center(2, 1)
-            end if
-            if ((center(3, 1) - domain_glb(3, 1)) <= radius) then
-                center(3, 2) = domain_glb(3, 2) + (center(3, 1) - domain_glb(3, 1))
-            else if ((domain_glb(3, 2) - center(3, 1)) <= radius) then
-                center(3, 2) = domain_glb(3, 1) - (domain_glb(3, 2) - center(3, 1))
-            else
-                center(3, 2) = center(3, 1)
-            end if
+            do i = 1, 3
+                call s_periodic_project_center(center(i, :), domain_glb(i, :), radius)
+            end do
         end if
 
         $:GPU_PARALLEL_LOOP(private='[i,j,k,ix,iy,iz,dist_vec,dist,dist_vec_temp,dist_temp]', &
@@ -589,7 +562,7 @@ contains
                     dist_vec(1) = x_cc(i) - center(1, 1)
                     dist_vec(2) = y_cc(j) - center(2, 1)
                     dist_vec(3) = z_cc(k) - center(3, 1)
-                    dist = sqrt(sum(dist_vec**2))
+                    dist = sum(dist_vec**2)
                     if (periodic_ibs) then
                         ! check all permutations of periodically projected ib to find minimum distance
                         do ix = 1, 2
@@ -598,7 +571,7 @@ contains
                                     dist_vec_temp(1) = x_cc(i) - center(1, ix)
                                     dist_vec_temp(2) = y_cc(j) - center(2, iy)
                                     dist_vec_temp(3) = z_cc(k) - center(3, iz)
-                                    dist_temp = sqrt(sum(dist_vec_temp**2))
+                                    dist_temp = sum(dist_vec_temp**2)
                                     if (dist_temp < dist) then
                                         dist = dist_temp
                                         dist_vec = dist_vec_temp
@@ -607,6 +580,7 @@ contains
                             end do
                         end do
                     end if
+                    dist = sqrt(dist)
                     levelset%sf(i, j, k, ib_patch_id) = dist - radius
                     if (f_approx_equal(dist, 0._wp)) then
                         levelset_norm%sf(i, j, k, ib_patch_id, :) = (/1, 0, 0/)
@@ -700,4 +674,19 @@ contains
 
     end subroutine s_cylinder_levelset
 
+    subroutine s_periodic_project_center(center, domain_glb, radius)
+        real(wp), dimension(2), intent(inout) :: center
+        real(wp), dimension(2), intent(in) :: domain_glb
+        real(wp), intent(in) :: radius
+
+        if ((center(1) - domain_glb(1)) <= radius) then
+            center(2) = domain_glb(2) + (center(1) - domain_glb(1))
+        else if ((domain_glb(2) - center(1)) <= radius) then
+            center(2) = domain_glb(1) - (domain_glb(2) - center(1))
+        else
+            center(2) = center(1)
+        end if
+
+    end subroutine s_periodic_project_center
+
 end module m_compute_levelset

src/common/m_ib_patches.fpp

@@ -166,21 +166,10 @@ contains
         ! that cell. If both queries check out, the primitive variables of
         ! the current patch are assigned to this cell.
 
-        if (periodic_ibs) then ! periodically wrap spheres around domain
-            if ((center(1, 1) - domain_glb(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 2) + (center(1, 1) - domain_glb(1, 1))
-            else if ((domain_glb(1, 2) - center(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 1) - (domain_glb(1, 2) - center(1, 1))
-            else
-                center(1, 2) = center(1, 1)
-            end if
-            if ((center(2, 1) - domain_glb(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 2) + (center(2, 1) - domain_glb(2, 1))
-            else if ((domain_glb(2, 2) - center(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 1) - (domain_glb(2, 2) - center(2, 1))
-            else
-                center(2, 2) = center(2, 1)
-            end if
+        if (periodic_ibs) then ! periodically wrap circles around domain
+            do i = 1, 2
+                call s_periodic_project_center(center(i, :), domain_glb(i, :), radius)
+            end do
 
             $:GPU_PARALLEL_LOOP(private='[i,j,ix,iy]', copy='[ib_markers_sf]',&
                       & copyin='[patch_id,center,r2]', collapse=2)
@@ -624,27 +613,9 @@ contains
         ! the current patch are assigned to this cell.
 
         if (periodic_ibs) then ! periodically wrap spheres around domain
-            if ((center(1, 1) - domain_glb(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 2) + (center(1, 1) - domain_glb(1, 1))
-            else if ((domain_glb(1, 2) - center(1, 1)) <= radius) then
-                center(1, 2) = domain_glb(1, 1) - (domain_glb(1, 2) - center(1, 1))
-            else
-                center(1, 2) = center(1, 1)
-            end if
-            if ((center(2, 1) - domain_glb(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 2) + (center(2, 1) - domain_glb(2, 1))
-            else if ((domain_glb(2, 2) - center(2, 1)) <= radius) then
-                center(2, 2) = domain_glb(2, 1) - (domain_glb(2, 2) - center(2, 1))
-            else
-                center(2, 2) = center(2, 1)
-            end if
-            if ((center(3, 1) - domain_glb(3, 1)) <= radius) then
-                center(3, 2) = domain_glb(3, 2) + (center(3, 1) - domain_glb(3, 1))
-            else if ((domain_glb(3, 2) - center(3, 1)) <= radius) then
-                center(3, 2) = domain_glb(3, 1) - (domain_glb(3, 2) - center(3, 1))
-            else
-                center(3, 2) = center(3, 1)
-            end if
+            do i = 1, 3
+                call s_periodic_project_center(center(i, :), domain_glb(i, :), radius)
+            end do
 
             $:GPU_PARALLEL_LOOP(collapse=3, private='[i,j,k,ix,iy,iz]', copy='[ib_markers_sf]', &
                               & copyin='[patch_id,center,r2]')

src/common/m_mpi_common.fpp

@@ -1928,6 +1928,9 @@ contains
 
 #ifdef MFC_MPI
         deallocate (buff_send, buff_recv)
+#ifndef MFC_POST_PROCESS
+        @:DEALLOCATE(domain_glb)
+#endif
 #endif
 
     end subroutine s_finalize_mpi_common_module

toolchain/mfc/case_validator.py

@@ -350,6 +350,7 @@ def check_ibm(self):
         n = self.get('n', 0)
         num_ibs = self.get('num_ibs', 0)
         periodic_ibs = self.get('periodic_ibs', 'F') == 'T'
+        cyl_coord = self.get('cyl_coord', 'F') == 'T'
 
         self.prohibit(ib and n <= 0,
                      "Immersed Boundaries do not work in 1D (requires n > 0)")
@@ -358,6 +359,10 @@ def check_ibm(self):
         self.prohibit(not ib and num_ibs > 0,
                      "num_ibs is set, but ib is not enabled")
         if periodic_ibs:
+            self.prohibit(not ib,
+                          "periodic_ibs is set to true but ib is not enabled")
+            self.prohibit(cyl_coord,
+                          "periodic_ibs not compatible with cyl_coords")
             for direction in ['x', 'y', 'z']:
                 for end in ['beg', 'end']:
                     bc_val = self.get(f'bc_{direction}%{end}')