immersed_m.f90 Source File

module immersed_m
    !! Abstract immersed boundary module to set boundary conditions at targets
    use precision_m, only : FP, FP_EPS
    use equilibrium_m, only : equilibrium_t
    use mesh_cart_m, only : mesh_cart_t
    use list_operations_m, only : unique_tuples
    implicit none
    
    type, public, abstract :: immersed_t
        !! Abstract immersed boundary type
        real(FP), public, allocatable, dimension(:) :: charfun
        !! Characteristic function of immersed boundary,
        !! commonly referred chi
        real(FP), public, allocatable, dimension(:) :: dirindfun
        !! Function indicating direction of magnetic field 
        !! (towards/away from target), commonly referred xi
        integer, public, allocatable, dimension(:) :: inds
        !! Mesh indices which lie in immersed boundary region
        integer, public :: n_inds
        !! Size of inds
        integer, public, allocatable, dimension(:) :: adj_inds
        !! Mesh indices, which are adjacent to the immersed boundary region
        integer, public :: n_adj_inds
        !! Size of adj_inds
    contains
        procedure(init), public, deferred :: init
        procedure(display), public, deferred :: display
        procedure, public :: build_inds
        procedure, public :: build_adj_inds
        procedure, public :: write_netcdf => write_netcdf_immersed
        procedure, public :: read_netcdf => read_netcdf_immersed
    end type

    abstract interface
    
        subroutine init(self, equi, mesh, filename)
            !! Initialises the immersed boundary fields charfun and dirind,
            !! depending on specialised implementation
            import :: immersed_t
            import :: equilibrium_t
            import :: mesh_cart_t
            class(immersed_t), intent(inout) :: self
            !! Instance of type
            class(equilibrium_t), intent(inout) :: equi
            !! Equilibrium
            type(mesh_cart_t), intent(in) :: mesh
            !! Mesh
            character(len=*), intent(in), optional :: filename
            !! Filename, where to read parameters from, 
            !! if not provided default parameters will be used
        end subroutine

        subroutine display(self)
            !! Writes basic information of immersed boundary to stdout
            import :: immersed_t
            class(immersed_t), intent(in) :: self
            !! Instance of type
        end subroutine

    end interface

    interface

        module subroutine write_netcdf_immersed(self, fgid)
            !! Writes immersed boundary data to NetCDF id
            class(immersed_t), intent(in) :: self
            !! Instance of class
            integer, intent(in) :: fgid
            !! NetCDF file or group id
        end subroutine

        module subroutine read_netcdf_immersed(self, fgid)
            !! Reads immersed boundary data from NetCDF id
            class(immersed_t), intent(inout) :: self
            !! Instance of class
            integer, intent(in) :: fgid
            !! NetCDF file or group id
        end subroutine

    end interface

contains

        subroutine build_inds(self, mesh)
            !! Computes inds.
            !! The construction is common to all derived immersed types
            class(immersed_t), intent(inout) :: self
            !! Instance of type
            type(mesh_cart_t), intent(in) :: mesh
            !! Mesh

            integer :: l, ic

            ! Find mesh indices that have a finite immersed value
            self%n_inds = count(self%charfun > FP_EPS)
            allocate(self%inds(self%n_inds))
            ic  = 0
            do l = 1, mesh%get_n_points()
                if (self%charfun(l) > FP_EPS) then
                    ic = ic + 1
                    self%inds(ic) = l
                endif
            enddo

        end subroutine

        subroutine build_adj_inds(self, mesh)
            !! Computes adj_inds.
            !! The construction is common to all derived immersed types
            class(immersed_t), intent(inout) :: self
            !! Instance of type
            type(mesh_cart_t), intent(in) :: mesh
            !! Mesh

            integer :: l, ic, k, ln, fill_it
            integer, allocatable, dimension(:) :: wrk

            integer, parameter, dimension(4) :: ing = [0, -1, 1, 0]
            integer, parameter, dimension(4) :: jng = [-1, 0, 0, 1]
            
            ! Find mesh points, that have a direct neighbour to the immersed
            ! region
            do fill_it = 1, 2
                ! fill_it = 1: Determines size and allocates wrk
                ! fill_it = 2: Fills wrk
                ic = 0
                mloop: do l = 1, mesh%get_n_points()
                    if (self%charfun(l) > FP_EPS) cycle mloop
                    do k = 1, 4
                        ln = mesh%get_index_neighbor(ing(k), jng(k), l)
                        if (ln == 0) cycle mloop
                        if (self%charfun(ln) > FP_EPS) then
                            ic = ic + 1
                            if (fill_it == 2) then
                                wrk(ic) = l
                            endif
                            cycle mloop
                        endif
                    enddo
                enddo mloop
                if (fill_it == 1) then
                    self%n_adj_inds = ic
                    allocate(wrk(self%n_adj_inds))
                endif
            enddo

            ! wrk may contain duplicate entries, unique them into adj_inds
            call unique_tuples(1, self%n_adj_inds, wrk)
            allocate(self%adj_inds(self%n_adj_inds))
            do ic = 1, self%n_adj_inds
                self%adj_inds(ic) = wrk(ic)
            enddo
            deallocate(wrk)

        end subroutine

end module