Merge branch 'levelset-refactor' of github.com:danieljvickers/MFC into levelset-refactor

danieljvickers · danieljvickers · commit 3e097d721ec7 · 2026-01-28T20:50:42.000-05:00
diff --git a/src/simulation/m_compute_levelset.fpp b/src/simulation/m_compute_levelset.fpp
@@ -77,7 +77,7 @@ contains
         integer :: global_id
         real(wp), dimension(3) :: dist_vec
 
-        real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
+        real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
         real(wp), dimension(1:2) :: center
         real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
 
@@ -87,14 +87,15 @@ contains
         center(2) = patch_ib(ib_patch_id)%y_centroid
         inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
         rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
+        offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)
 
         $:GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,dist_vec,dist,global_dist,global_id]', &
-                  & copyin='[ib_patch_id,center,rotation,inverse_rotation,airfoil_grid_u,airfoil_grid_l]', collapse=2)
+                  & copyin='[ib_patch_id,center,rotation,inverse_rotation,offset,airfoil_grid_u,airfoil_grid_l]', collapse=2)
         do i = 0, m
             do j = 0, n
                 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
                 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
-                xy_local = xy_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
+                xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
 
                 if (xy_local(2) >= 0._wp) then
                     ! finds the location on the airfoil grid with the minimum distance (closest)
@@ -165,7 +166,7 @@ contains
         real(wp) :: lz, z_max, z_min
         real(wp), dimension(3) :: dist_vec
 
-        real(wp), dimension(1:3) :: xyz_local, center !< x, y, z coordinates in local IB frame
+        real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame
         real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
 
         real(wp) :: length_z
@@ -178,19 +179,20 @@ contains
         lz = patch_ib(ib_patch_id)%length_z
         inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
         rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
+        offset(:) = patch_ib(ib_patch_id)%centroid_offset(:)
 
         z_max = center(3) + lz/2
         z_min = center(3) - lz/2
 
         $:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,dist_vec,dist,global_dist,global_id,dist_side,dist_surf]', &
-                  & copyin='[ib_patch_id,center,rotation,inverse_rotation,airfoil_grid_u,airfoil_grid_l,z_min,z_max]', collapse=3)
+                  & copyin='[ib_patch_id,center,rotation,inverse_rotation,offset,airfoil_grid_u,airfoil_grid_l,z_min,z_max]', collapse=3)
         do l = 0, p
             do j = 0, n
                 do i = 0, m
 
                     xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
                     xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
-                    xyz_local = xyz_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
+                    xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
 
                     if (xyz_local(2) >= center(2)) then
                         do k = 1, Np
@@ -239,12 +241,12 @@ contains
                     if (dist_side < dist_surf) then
                         levelset%sf(i, j, l, ib_patch_id) = dist_side
                         if (f_approx_equal(dist_side, abs(z_cc(l) - z_min))) then
-                            levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, -1/)
+                            levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0._wp, 0._wp, -1._wp/)
                         else
-                            levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, 1/)
+                            levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0._wp, 0._wp, 1._wp/)
                         end if
                         levelset_norm%sf(i, j, l, ib_patch_id, :) = &
-                            matmul(rotation, levelset_norm%sf(i, j, l, ib_patch_id, :)/dist_surf)
+                            matmul(rotation, levelset_norm%sf(i, j, l, ib_patch_id, :))
                     else
                         levelset%sf(i, j, l, ib_patch_id) = dist_surf
                         if (f_approx_equal(dist_surf, 0._wp)) then
diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp
@@ -192,7 +192,7 @@ contains
         integer :: i, j, k
         integer :: Np1, Np2
 
-        real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
+        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
         real(wp), dimension(1:3, 1:3) :: inverse_rotation
 
@@ -203,6 +203,7 @@ contains
         ma = patch_ib(patch_id)%m
         ta = patch_ib(patch_id)%t
         inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
+        offset(:) = patch_ib(patch_id)%centroid_offset(:)
 
         ! rank(dx) is not consistent between pre_process and simulation. This IFDEF prevents compilation errors
 #ifdef MFC_PRE_PROCESS
@@ -274,12 +275,12 @@ contains
         end if
 
         $:GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,f]', copy='[ib_markers_sf]',&
-                  & copyin='[patch_id,center,inverse_rotation,ma,ca_in,airfoil_grid_u,airfoil_grid_l]', collapse=2)
+                  & copyin='[patch_id,center,inverse_rotation,offset,ma,ca_in,airfoil_grid_u,airfoil_grid_l]', collapse=2)
         do j = 0, n
             do i = 0, m
                 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
                 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates
-                xy_local = xy_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
+                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
                     xa = xy_local(1)/ca_in
@@ -345,7 +346,7 @@ contains
         integer :: i, j, k, l
         integer :: Np1, Np2
 
-        real(wp), dimension(1:3) :: xyz_local, center !< x, y, z coordinates in local IB frame
+        real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame
         real(wp), dimension(1:3, 1:3) :: inverse_rotation
 
         center(1) = patch_ib(patch_id)%x_centroid
@@ -357,6 +358,7 @@ contains
         ma = patch_ib(patch_id)%m
         ta = patch_ib(patch_id)%t
         inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
+        offset(:) = patch_ib(patch_id)%centroid_offset(:)
 
         ! rank(dx) is not consistent between pre_process and simulation. This IFDEF prevents compilation errors
 #ifdef MFC_PRE_PROCESS
@@ -428,13 +430,13 @@ contains
         end if
 
         $:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,f]', copy='[ib_markers_sf]',&
-                  & copyin='[patch_id,center,inverse_rotation,ma,ca_in,airfoil_grid_u,airfoil_grid_l]', collapse=3)
+                  & copyin='[patch_id,center,inverse_rotation,offset,ma,ca_in,airfoil_grid_u,airfoil_grid_l]', collapse=3)
         do l = 0, p
             do j = 0, n
                 do i = 0, m
                     xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
                     xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
-                    xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
+                    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
 
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
@@ -505,9 +505,15 @@ contains
                                .or. temp_loc > s_cc(index + 1)))
                         index = index + dir
                         if (index < -buff_size .or. index > bound) then
-                            print *, "temp_loc=", temp_loc, " s_cc(index)=", s_cc(index), " s_cc(index+1)=", s_cc(index + 1)
-                            print *, "Increase buff_size further in m_helper_basic (currently set to a minimum of 10)"
-                            error stop "Increase buff_size"
+                            print *, "A required image point is not located in this computational domain."
+                            print *, "Ghost Point is located at ", [x_cc(i), y_cc(j), z_cc(k)], " while moving in dimension ", dim
+                            print *, "We are searching for image point at ", ghost_points_in(q)%ip_loc(:)
+                            print *, "We can only support points located inside the box from ", [x_cc(-buff_size), y_cc(-buff_size), z_cc(-buff_size)]
+                            print *, "To ", [x_cc( m + buff_size - 1), y_cc(n + buff_size - 1), z_cc(p + buff_size - 1)]
+                            print *, "Image point is located approximately ", (ghost_points_in(q)%loc(dim) - ghost_points_in(q)%ip_loc(dim)) / (s_cc(1) - s_cc(0)), " grid cells away"
+                            print *, "Levelset ", dist, " and Norm: ", norm(:)
+                            print *, "A short term fix may include increasing buff_size further in m_helper_basic (currently set to a minimum of 10)"
+                            error stop "Ghost Point and Image Point on Different Processors"
                         end if
                     end do
                     ghost_points_in(q)%ip_grid(dim) = index
