m_compute_levelset.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!>
!!@file
!! @brief Contains module m_compute_levelset
 
# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
 
# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 207 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 232 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 243 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 245 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 283 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 293 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 303 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 312 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 329 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 346 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 358 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 364 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 370 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 376 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
! New line at end of file is required for FYPP
# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
 
# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 213 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 241 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 256 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 266 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 275 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 291 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 301 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 308 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
! New line at end of file is required for FYPP
# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
 
# 21 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 50 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 104 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 119 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 171 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 182 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 204 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 214 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 225 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 236 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 246 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 252 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 258 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 264 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 270 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 272 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
# 273 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
 
# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Caution:
! This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI rank.
! That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0.
! For an example see misc/nvidia_uvm/bind.sh.
# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 167 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
! New line at end of file is required for FYPP
# 6 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp" 2
 
!> @brief Computes signed-distance level-set fields and surface normals for immersed-boundary patch geometries
module m_compute_levelset
 
    use m_ib_patches           !< the ib patch parameters
 
    use m_model                !< subroutine(s) related to stl files
 
    use m_derived_types        !< definitions of the derived types
 
    use m_global_parameters    !< definitions of the global parameters
 
    use m_mpi_proxy            !< message passing interface (mpi) module proxy
 
    use m_helper_basic         !< functions to compare floating point numbers
 
    implicit none
 
    private; public :: s_apply_levelset
 
contains
 
    !> @brief Dispatches level-set distance and normal computations for all ghost points based on their patch geometry type.
    impure subroutine s_apply_levelset(gps, num_gps)
 
        type(ghost_point), dimension(:), intent(inout) :: gps
        integer, intent(in) :: num_gps
 
        integer :: i, patch_id, patch_geometry
 
        !  3D Patch Geometries
        if (p > 0) then
 
 
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if defined(MFC_OpenACC)
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc parallel loop gang vector default(present) private(i, patch_id, patch_geometry) copy(gps) copyin(patch_ib(1:num_ibs), Np) 
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif defined(MFC_OpenMP)
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, patch_id, patch_geometry) map(tofrom:gps) map(to:patch_ib(1:num_ibs), Np) 
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
# 39 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
            do i = 1, num_gps
 
                patch_id = gps(i)%ib_patch_id
                patch_geometry = patch_ib(patch_id)%geometry
 
                if (patch_geometry == 8) then
                    call s_sphere_levelset(gps(i))
                elseif (patch_geometry == 9) then
                    call s_cuboid_levelset(gps(i))
                elseif (patch_geometry == 10) then
                    call s_cylinder_levelset(gps(i))
                elseif (patch_geometry == 11) then
                    call s_3d_airfoil_levelset(gps(i))
                elseif (patch_geometry == 12) then
                    call s_model_levelset(gps(i))
                end if
            end do
 
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if defined(MFC_OpenACC)
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc end parallel loop
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif defined(MFC_OpenMP)
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp end target teams loop
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
# 57 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
 
            ! 2D Patch Geometries
        elseif (n > 0) then
 
 
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if defined(MFC_OpenACC)
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc parallel loop gang vector default(present) private(i, patch_id, patch_geometry) copy(gps) copyin(Np, patch_ib(1:num_ibs)) 
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif defined(MFC_OpenMP)
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, patch_id, patch_geometry) map(tofrom:gps) map(to:Np, patch_ib(1:num_ibs)) 
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
# 62 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
            do i = 1, num_gps
 
                patch_id = gps(i)%ib_patch_id
                patch_geometry = patch_ib(patch_id)%geometry
 
                if (patch_geometry == 2) then
                    call s_circle_levelset(gps(i))
                elseif (patch_geometry == 3) then
                    call s_rectangle_levelset(gps(i))
                elseif (patch_geometry == 4) then
                    call s_airfoil_levelset(gps(i))
                elseif (patch_geometry == 5) then
                    call s_model_levelset(gps(i))
                elseif (patch_geometry == 6) then
                    call s_ellipse_levelset(gps(i))
                end if
            end do
 
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if defined(MFC_OpenACC)
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc end parallel loop
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif defined(MFC_OpenMP)
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp end target teams loop
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
# 80 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
 
        end if
 
    end subroutine s_apply_levelset
 
    !> @brief Computes the signed distance and outward normal from a ghost point to a circular immersed boundary.
    subroutine s_circle_levelset(gp)
 
 
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 89 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: radius, dist
        real(wp), dimension(2) :: center
        real(wp), dimension(3) :: dist_vec
 
        integer :: i, j, ib_patch_id !< Loop index variables
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
 
        radius = patch_ib(ib_patch_id)%radius
 
        dist_vec(1) = x_cc(i) - patch_ib(ib_patch_id)%x_centroid - real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        dist_vec(2) = y_cc(j) - patch_ib(ib_patch_id)%y_centroid - real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        dist_vec(3) = 0._wp
        dist = sqrt(sum(dist_vec**2))
 
        gp%levelset = dist - radius
        if (f_approx_equal(dist, 0._wp)) then
            gp%levelset_norm = 0._wp
        else
            gp%levelset_norm = dist_vec(:)/dist
        end if
 
    end subroutine s_circle_levelset
 
    !> @brief Computes the signed distance and outward normal from a ghost point to a 2D NACA airfoil surface.
    subroutine s_airfoil_levelset(gp)
 
 
# 122 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 122 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 122 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 122 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 122 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 122 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 122 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: dist, global_dist
        integer :: global_id
        real(wp), dimension(3) :: dist_vec
 
        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
 
        integer :: i, j, k, ib_patch_id !< Loop index variables
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
 
        center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        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(:)
 
        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 - 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)
            do k = 1, np
                dist_vec(1) = xy_local(1) - airfoil_grid_u(k)%x
                dist_vec(2) = xy_local(2) - airfoil_grid_u(k)%y
                dist_vec(3) = 0._wp
                dist = sqrt(sum(dist_vec**2))
                if (k == 1) then
                    global_dist = dist
                    global_id = k
                else
                    if (dist < global_dist) then
                        global_dist = dist
                        global_id = k
                    end if
                end if
            end do
            dist_vec(1) = xy_local(1) - airfoil_grid_u(global_id)%x
            dist_vec(2) = xy_local(2) - airfoil_grid_u(global_id)%y
            dist_vec(3) = 0
            dist = global_dist
        else
            do k = 1, np
                dist_vec(1) = xy_local(1) - airfoil_grid_l(k)%x
                dist_vec(2) = xy_local(2) - airfoil_grid_l(k)%y
                dist_vec(3) = 0
                dist = sqrt(sum(dist_vec**2))
                if (k == 1) then
                    global_dist = dist
                    global_id = k
                else
                    if (dist < global_dist) then
                        global_dist = dist
                        global_id = k
                    end if
                end if
            end do
            dist_vec(1) = xy_local(1) - airfoil_grid_l(global_id)%x
            dist_vec(2) = xy_local(2) - airfoil_grid_l(global_id)%y
            dist_vec(3) = 0
            dist = global_dist
        end if
 
        gp%levelset = dist
        if (f_approx_equal(dist, 0._wp)) then
            gp%levelset_norm = 0._wp
        else
            gp%levelset_norm = matmul(rotation, dist_vec(:))/dist ! convert the normal vector back to global grid coordinates
        end if
 
    end subroutine s_airfoil_levelset
 
    !> @brief Computes the signed distance and outward normal from a ghost point to a 3D extruded airfoil surface including spanwise end caps.
    subroutine s_3d_airfoil_levelset(gp)
 
 
# 205 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 205 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 205 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 205 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 205 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 205 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 205 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: dist, dist_surf, dist_side, global_dist
        integer :: global_id
        real(wp) :: lz, z_max, z_min
        real(wp), dimension(3) :: dist_vec
 
        real(wp), dimension(1:3) :: xyz_local, center, offset, normal !< x, y, z coordinates in local IB frame
        real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
 
        real(wp) :: length_z
 
        integer :: i, j, k, l, ib_patch_id !< Loop index variables
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
        l = gp%loc(3)
 
        center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
        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 = lz/2
        z_min = -lz/2
 
        xyz_local = [x_cc(i), y_cc(j), z_cc(l)] - center
        xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
        xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
 
        if (xyz_local(2) >= 0._wp) then
            do k = 1, np
                dist_vec(1) = xyz_local(1) - airfoil_grid_u(k)%x
                dist_vec(2) = xyz_local(2) - airfoil_grid_u(k)%y
                dist_vec(3) = 0._wp
                dist_surf = sqrt(sum(dist_vec**2))
                if (k == 1) then
                    global_dist = dist_surf
                    global_id = k
                else
                    if (dist_surf < global_dist) then
                        global_dist = dist_surf
                        global_id = k
                    end if
                end if
            end do
            dist_vec(1) = xyz_local(1) - airfoil_grid_u(global_id)%x
            dist_vec(2) = xyz_local(2) - airfoil_grid_u(global_id)%y
            dist_vec(3) = 0._wp
            dist_surf = global_dist
        else
            do k = 1, np
                dist_vec(1) = xyz_local(1) - airfoil_grid_l(k)%x
                dist_vec(2) = xyz_local(2) - airfoil_grid_l(k)%y
                dist_vec(3) = 0
                dist_surf = sqrt(sum(dist_vec**2))
                if (k == 1) then
                    global_dist = dist_surf
                    global_id = k
                else
                    if (dist_surf < global_dist) then
                        global_dist = dist_surf
                        global_id = k
                    end if
                end if
            end do
            dist_vec(1) = xyz_local(1) - airfoil_grid_l(global_id)%x
            dist_vec(2) = xyz_local(2) - airfoil_grid_l(global_id)%y
            dist_vec(3) = 0._wp
            dist_surf = global_dist
        end if
 
        dist_side = min(abs(xyz_local(3) - z_min), abs(z_max - xyz_local(3)))
 
        if (dist_side < dist_surf) then
            gp%levelset = dist_side
            normal = 0._wp
            if (f_approx_equal(dist_side, abs(xyz_local(3) - z_min))) then
                normal(3) = -1._wp
            else
                normal(3) = 1._wp
            end if
            gp%levelset_norm = matmul(rotation, normal)
        else
            gp%levelset = dist_surf
            if (f_approx_equal(dist_surf, 0._wp)) then
                gp%levelset_norm = 0._wp
            else
                gp%levelset_norm = matmul(rotation, dist_vec(:)/dist_surf)
            end if
        end if
 
    end subroutine s_3d_airfoil_levelset
 
    !>  Subroutine for computing the levelset values at a ghost point belonging to the rectangle IB
    subroutine s_rectangle_levelset(gp)
 
 
# 308 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 308 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 308 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 308 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 308 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 308 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 308 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: top_right(2), bottom_left(2)
        real(wp) :: min_dist
        real(wp) :: side_dists(4)
 
        real(wp) :: length_x, length_y
        real(wp), dimension(1:3) :: xy_local, dist_vec !< x and y coordinates in local IB frame
        real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
        real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
 
        integer :: i, j, k !< Loop index variables
        integer :: idx !< Shortest path direction indicator
        integer :: ib_patch_id !< patch ID
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
 
        length_x = patch_ib(ib_patch_id)%length_x
        length_y = patch_ib(ib_patch_id)%length_y
        center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
        rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
 
        top_right(1) = length_x/2
        top_right(2) = length_y/2
        bottom_left(1) = -length_x/2
        bottom_left(2) = -length_y/2
 
        ! convert grid to local coordinates
        xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
        xy_local = matmul(inverse_rotation, xy_local)
 
        side_dists(1) = bottom_left(1) - xy_local(1)
        side_dists(2) = top_right(1) - xy_local(1)
        side_dists(3) = bottom_left(2) - xy_local(2)
        side_dists(4) = top_right(2) - xy_local(2)
        min_dist = side_dists(1)
        idx = 1
 
        do k = 2, 4
            if (abs(side_dists(k)) < abs(min_dist)) then
                idx = k
                min_dist = side_dists(idx)
            end if
        end do
 
        gp%levelset = side_dists(idx)
        dist_vec = 0._wp
        if (.not. f_approx_equal(side_dists(idx), 0._wp)) then
            if (idx == 1 .or. idx == 2) then
                ! vector points along the x axis
                dist_vec(1) = side_dists(idx)/abs(side_dists(idx))
            else
                ! vector points along the y axis
                dist_vec(2) = side_dists(idx)/abs(side_dists(idx))
            end if
            ! convert the normal vector back into the global coordinate system
            gp%levelset_norm = matmul(rotation, dist_vec)
        else
            gp%levelset_norm = 0._wp
        end if
 
    end subroutine s_rectangle_levelset
 
    !> @brief Computes the signed distance and outward normal from a ghost point to an elliptical immersed boundary via a quadratic projection.
    subroutine s_ellipse_levelset(gp)
 
 
# 380 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 380 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 380 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 380 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 380 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 380 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 380 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: ellipse_coeffs(2) ! a and b in the ellipse equation
        real(wp) :: quadratic_coeffs(3) ! A, B, C in the quadratic equation to compute levelset
 
        real(wp) :: length_x, length_y
        real(wp), dimension(1:3) :: xy_local, normal_vector !< x and y coordinates in local IB frame
        real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
        real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
 
        integer :: i, j, k !< Loop index variables
        integer :: idx !< Shortest path direction indicator
        integer :: ib_patch_id !< patch ID
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
 
        length_x = patch_ib(ib_patch_id)%length_x
        length_y = patch_ib(ib_patch_id)%length_y
        center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
        rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
 
        ellipse_coeffs(1) = 0.5_wp*length_x
        ellipse_coeffs(2) = 0.5_wp*length_y
 
        xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
        xy_local = matmul(inverse_rotation, xy_local)
 
        normal_vector = xy_local
        normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transformation method
        normal_vector = normal_vector/sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
        gp%levelset_norm = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame
 
        ! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3
        quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2
        quadratic_coeffs(2) = 2._wp*((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2)))
        quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp
 
        ! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
        gp%levelset = -0.5_wp*(-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp*quadratic_coeffs(1)*quadratic_coeffs(3)))/quadratic_coeffs(1)
 
    end subroutine s_ellipse_levelset
 
    !> @brief Computes the signed distance and outward normal from a ghost point to the nearest face of a cuboid immersed boundary.
    subroutine s_cuboid_levelset(gp)
 
 
# 431 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 431 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 431 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 431 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 431 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 431 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 431 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: Right, Left, Bottom, Top, Front, Back
        real(wp) :: min_dist
        real(wp) :: dist_left, dist_right, dist_bottom, dist_top, dist_back, dist_front
 
        real(wp), dimension(3) :: center
        real(wp) :: length_x, length_y, length_z
        real(wp), dimension(1:3) :: xyz_local, dist_vec !< x and y coordinates in local IB frame
        real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
 
        integer :: i, j, k !< Loop index variables
        integer :: ib_patch_id !< patch ID
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
        k = gp%loc(3)
 
        length_x = patch_ib(ib_patch_id)%length_x
        length_y = patch_ib(ib_patch_id)%length_y
        length_z = patch_ib(ib_patch_id)%length_z
 
        center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
 
        inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
        rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)
 
        right = length_x/2
        left = -length_x/2
        top = length_y/2
        bottom = -length_y/2
        front = length_z/2
        back = -length_z/2
 
        xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
        xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinate
 
        dist_left = left - xyz_local(1)
        dist_right = xyz_local(1) - right
        dist_bottom = bottom - xyz_local(2)
        dist_top = xyz_local(2) - top
        dist_back = back - xyz_local(3)
        dist_front = xyz_local(3) - front
 
        min_dist = min(abs(dist_left), abs(dist_right), abs(dist_bottom), &
                       abs(dist_top), abs(dist_back), abs(dist_front))
        dist_vec = 0._wp
 
        if (f_approx_equal(min_dist, abs(dist_left))) then
            gp%levelset = dist_left
            if (.not. f_approx_equal(dist_left, 0._wp)) then
                dist_vec(1) = dist_left/abs(dist_left)
            end if
        else if (f_approx_equal(min_dist, abs(dist_right))) then
            gp%levelset = dist_right
            if (.not. f_approx_equal(dist_right, 0._wp)) then
                dist_vec(1) = -dist_right/abs(dist_right)
            end if
        else if (f_approx_equal(min_dist, abs(dist_bottom))) then
            gp%levelset = dist_bottom
            if (.not. f_approx_equal(dist_bottom, 0._wp)) then
                dist_vec(2) = dist_bottom/abs(dist_bottom)
            end if
        else if (f_approx_equal(min_dist, abs(dist_top))) then
            gp%levelset = dist_top
            if (.not. f_approx_equal(dist_top, 0._wp)) then
                dist_vec(2) = -dist_top/abs(dist_top)
            end if
        else if (f_approx_equal(min_dist, abs(dist_back))) then
            gp%levelset = dist_back
            if (.not. f_approx_equal(dist_back, 0._wp)) then
                dist_vec(3) = dist_back/abs(dist_back)
            end if
        else if (f_approx_equal(min_dist, abs(dist_front))) then
            gp%levelset = dist_front
            if (.not. f_approx_equal(dist_front, 0._wp)) then
                dist_vec(3) = -dist_front/abs(dist_front)
            end if
        end if
 
        gp%levelset_norm = matmul(rotation, dist_vec)
 
    end subroutine s_cuboid_levelset
 
    !> @brief Computes the signed distance and outward normal from a ghost point to a spherical immersed boundary.
    subroutine s_sphere_levelset(gp)
 
 
# 523 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 523 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 523 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 523 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 523 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 523 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 523 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: radius, dist
        real(wp), dimension(3) :: dist_vec, center, periodicity
 
        integer :: i, j, k, ib_patch_id !< Loop index variables
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
        k = gp%loc(3)
 
        radius = patch_ib(ib_patch_id)%radius
        periodicity(1) = real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        periodicity(2) = real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        periodicity(3) = real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
        center(1) = patch_ib(ib_patch_id)%x_centroid
        center(2) = patch_ib(ib_patch_id)%y_centroid
        center(3) = patch_ib(ib_patch_id)%z_centroid
        center = center + periodicity
 
        dist_vec(1) = x_cc(i) - center(1)
        dist_vec(2) = y_cc(j) - center(2)
        dist_vec(3) = z_cc(k) - center(3)
        dist = sqrt(sum(dist_vec**2))
        gp%levelset = dist - radius
        if (f_approx_equal(dist, 0._wp)) then
            gp%levelset_norm = (/1, 0, 0/)
        else
            gp%levelset_norm = dist_vec(:)/dist
        end if
 
    end subroutine s_sphere_levelset
 
    !> @brief Computes the signed distance and outward normal from a ghost point to a cylindrical immersed boundary surface and end caps.
    subroutine s_cylinder_levelset(gp)
 
 
# 562 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 562 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 562 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 562 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 562 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 562 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 562 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        real(wp) :: radius
        real(wp), dimension(3) :: dist_sides_vec, dist_surface_vec, length
        real(wp), dimension(2) :: boundary
        real(wp) :: dist_side, dist_surface, side_pos
        integer :: i, j, k !< Loop index variables
        integer :: ib_patch_id !< patch ID
 
        real(wp), dimension(1:3) :: xyz_local, center !< x and y coordinates in local IB frame
        real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation
 
        ib_patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
        k = gp%loc(3)
 
        radius = patch_ib(ib_patch_id)%radius
        center(1) = patch_ib(ib_patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        center(2) = patch_ib(ib_patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        center(3) = patch_ib(ib_patch_id)%z_centroid + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
        length(1) = patch_ib(ib_patch_id)%length_x
        length(2) = patch_ib(ib_patch_id)%length_y
        length(3) = 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(:, :)
 
        if (.not. f_approx_equal(length(1), 0._wp)) then
            boundary(1) = -0.5_wp*length(1)
            boundary(2) = 0.5_wp*length(1)
            dist_sides_vec = (/1, 0, 0/)
            dist_surface_vec = (/0, 1, 1/)
        else if (.not. f_approx_equal(length(2), 0._wp)) then
            boundary(1) = -0.5_wp*length(2)
            boundary(2) = 0.5_wp*length(2)
            dist_sides_vec = (/0, 1, 0/)
            dist_surface_vec = (/1, 0, 1/)
        else if (.not. f_approx_equal(length(3), 0._wp)) then
            boundary(1) = -0.5_wp*length(3)
            boundary(2) = 0.5_wp*length(3)
            dist_sides_vec = (/0, 0, 1/)
            dist_surface_vec = (/1, 1, 0/)
        end if
 
        xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
        xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
 
        ! get distance to flat edge of cylinder
        side_pos = dot_product(xyz_local, dist_sides_vec)
        dist_side = min(abs(side_pos - boundary(1)), &
                        abs(boundary(2) - side_pos))
        ! get distance to curved side of cylinder
        dist_surface = norm2(xyz_local*dist_surface_vec) &
                       - radius
 
        if (dist_side < abs(dist_surface)) then
            ! if the closest edge is flat
            gp%levelset = -dist_side
            if (f_approx_equal(dist_side, abs(side_pos - boundary(1)))) then
                gp%levelset_norm = matmul(rotation, -dist_sides_vec)
            else
                gp%levelset_norm = matmul(rotation, dist_sides_vec)
            end if
        else
            gp%levelset = dist_surface
            xyz_local = xyz_local*dist_surface_vec
            xyz_local = xyz_local/max(norm2(xyz_local), sgm_eps)
            gp%levelset_norm = matmul(rotation, xyz_local)
        end if
 
    end subroutine s_cylinder_levelset
 
    !> The STL patch is a 2/3D geometry that is imported from an STL file.
    !! @param gp Ghost point to compute levelset for
    subroutine s_model_levelset(gp)
 
 
# 641 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#if MFC_OpenACC
# 641 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$acc routine seq 
# 641 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#elif MFC_OpenMP
# 641 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
    
# 641 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
 
# 641 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
!$omp declare target device_type(any) 
# 641 "/home/runner/work/MFC/MFC/src/simulation/m_compute_levelset.fpp"
#endif
 
        type(ghost_point), intent(inout) :: gp
 
        integer :: i, j, k, patch_id, boundary_edge_count, total_vertices
        real(wp), dimension(1:3) :: center, xyz_local
        real(wp) :: normals(1:3) !< Boundary normal buffer
        real(wp) :: distance
        real(wp), dimension(1:3, 1:3) :: inverse_rotation, rotation
 
        patch_id = gp%ib_patch_id
        i = gp%loc(1)
        j = gp%loc(2)
        k = gp%loc(3)
 
        ! load in model values
        boundary_edge_count = gpu_boundary_edge_count(patch_id)
        total_vertices = gpu_total_vertices(patch_id)
 
        center = 0._wp
        if (.not. f_is_default(patch_ib(patch_id)%x_centroid)) center(1) = patch_ib(patch_id)%x_centroid + real(gp%x_periodicity, wp)*(x_domain%end - x_domain%beg)
        if (.not. f_is_default(patch_ib(patch_id)%y_centroid)) center(2) = patch_ib(patch_id)%y_centroid + real(gp%y_periodicity, wp)*(y_domain%end - y_domain%beg)
        if (p > 0) then
            if (.not. f_is_default(patch_ib(patch_id)%z_centroid)) center(3) = patch_ib(patch_id)%z_centroid + real(gp%z_periodicity, wp)*(z_domain%end - z_domain%beg)
        end if
 
        inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
        rotation(:, :) = patch_ib(patch_id)%rotation_matrix(:, :)
 
        ! determine where we are located in space
        xyz_local = (/x_cc(i) - center(1), y_cc(j) - center(2), 0._wp/)
        if (p > 0) then
            xyz_local(3) = z_cc(k) - center(3)
        end if
        xyz_local = matmul(inverse_rotation, xyz_local)
 
        ! 3D models
        if (p > 0) then
            ! Get the boundary normals and shortest distance between the cell center and the model boundary
            call s_distance_normals_3d(gpu_ntrs(patch_id), patch_id, xyz_local, normals, distance)
 
            ! Get the shortest distance between the cell center and the model boundary
            gp%levelset = distance
            gp%levelset = -abs(gp%levelset)
 
            ! Assign the levelset_norm
            gp%levelset_norm = matmul(rotation, normals(1:3))
        else
            ! 2D models
            call s_distance_normals_2d(patch_id, &
                                       boundary_edge_count, &
                                       xyz_local, normals, distance)
            gp%levelset = -abs(distance)
            gp%levelset_norm = matmul(rotation, normals(1:3))
        end if
 
    end subroutine s_model_levelset
 
end module m_compute_levelset