mask_data_m.f90 Source File

module mask_data_m
    !! Mask arrays based on static_data set
    use precision_grillix_m, only : GP, GP_EPS
    use equilibrium_m, only : equilibrium_t
    use mesh_cart_m, only : mesh_cart_t
    use penalisation_m, only : parbnd_immersed_t => penalisation_t
    implicit none

    type, public ::  mask_data_t
        !! Mask arrays
        logical, dimension(:), allocatable :: is_inner
        !! True, if point is an inner computational point
        logical, dimension(:), allocatable :: is_boundary
        !! True, if point is a (perpendicular) boundary point
        logical, dimension(:), allocatable :: is_ghost
        !! True, if point is a (perpendicular) ghost point
        logical, dimension(:), allocatable :: is_in_vessel
        !! True if point is inside vessel
        logical, dimension(:), allocatable :: is_parbnd_immersed
        !! True if point is in region of finite immersion
    contains
        procedure, public :: build => build_mask_data
        procedure, public :: write_netcdf => write_netcdf_mask_data
        procedure, public :: read_netcdf => read_netcdf_mask_data
        final :: destructor_mask_data_t
    end type

    interface
        module subroutine write_netcdf_mask_data(self, fgid)
            !! Writes mask_data to NetCDF file
            class(mask_data_t), intent(in) :: self
            !! Instance of class
            integer, intent(in) :: fgid
            !! File or group id
        end subroutine

        module subroutine read_netcdf_mask_data(self, fgid)
            !! Reads mask_data from NetCDF file
            class(mask_data_t), intent(inout) :: self
            !! Instance of class
            integer, intent(in) :: fgid
            !! File or group id
        end subroutine
    end interface

contains

    subroutine build_mask_data(self, equi, mesh, parbnd_immersed, threshold_immersed)
        !! Builds mask arrays
        class(mask_data_t), intent(inout) :: self
        !! Instance of type
        class(equilibrium_t), intent(in) :: equi
        !! Equilibrium
        type(mesh_cart_t), intent(in) :: mesh
        !! Mesh
        type(parbnd_immersed_t), intent(in), optional :: parbnd_immersed
        !! Parallel immersion
        real(GP), intent(in), optional :: threshold_immersed
        !! Threshold above which point is considered an immersed point
        !! Default: GP_EPS

        integer :: np, l, ki
        real(GP) :: x, y, phi, thresh_im

        np = mesh%get_n_points()
        phi = mesh%get_phi()
        allocate(self%is_inner(np), source=.false.)
        allocate(self%is_boundary(np), source=.false.)
        allocate(self%is_ghost(np), source=.false.)
        allocate(self%is_in_vessel(np), source=.false.)
        allocate(self%is_parbnd_immersed(np), source=.false.)

        !$omp parallel default(none) private(l, x, y) &
        !$omp shared(self, np, equi, mesh, phi)
        !$omp do
        do l = 1, np
            if (mesh%is_inner_point(l)) then
                self%is_inner(l) = .true.
            endif
            if (mesh%is_boundary_point(l)) then
                self%is_boundary(l) = .true.
            endif
            if (mesh%is_ghost_point(l)) then
                self%is_ghost(l) = .true.
            endif
            
            x = mesh%get_x(l)
            y = mesh%get_y(l)

            self%is_in_vessel(l) = equi%in_vessel(x, y, phi)
            self%is_parbnd_immersed(l) = .false. ! first touch
        enddo
        !$omp end do
        !$omp end parallel

        if (present(parbnd_immersed)) then
            if (present(threshold_immersed)) then
                thresh_im = threshold_immersed
            else
                thresh_im = GP_EPS
            endif
            !$omp parallel default(none) private(l, ki) &
            !$omp shared(self, np, mesh, parbnd_immersed, thresh_im)
            !$omp do
            do ki = 1, mesh%get_n_points_inner()
                l = mesh%inner_indices(ki)
                if (parbnd_immersed%charfun(ki) > thresh_im) then
                    self%is_parbnd_immersed(l) = .true.
                endif
            enddo
            !$omp end do
            !$omp end parallel
        endif

    end subroutine

    subroutine destructor_mask_data_t(self)
        type(mask_data_t), intent(inout) :: self

        if (allocated(self%is_inner)) deallocate(self%is_inner)
        if (allocated(self%is_boundary)) deallocate(self%is_boundary)
        if (allocated(self%is_ghost)) deallocate(self%is_ghost)
        if (allocated(self%is_in_vessel)) deallocate(self%is_in_vessel)
        if (allocated(self%is_parbnd_immersed)) deallocate(self%is_parbnd_immersed)

    end subroutine
    
end module