parallel_target_flux_build_s.f90 Source File

submodule (parallel_target_flux_m) parallel_target_flux_build_s
    !! Routines for building parallel target markers
    implicit none

contains

    module subroutine build_target_markers(self, comm_handler, mesh, map, penalisation, staggered, &
                                           pen_threshold, dbgout)
        class(parallel_target_flux_t), intent(inout) :: self
        type(comm_handler_t), intent(in) :: comm_handler
        type(mesh_cart_t), intent(in) :: mesh
        type(parallel_map_t), intent(in) :: map
        type(penalisation_t), intent(in) :: penalisation
        logical, intent(in) :: staggered
        real(GP), intent(in) :: pen_threshold
        integer, intent(in), optional :: dbgout

        integer :: i, l, ki, kb, n_fwd, n_bwd, outi, nrecv
        real(GP) :: pen_val_fwd, pen_val_bwd
        real(GP), dimension(mesh%get_n_points()) :: pen_val_mask
        real(GP), dimension(mesh%get_n_points()) :: pen_vals, &
                                buffer_inds_fwd, buffer_weights_fwd, &
                                buffer_inds_bwd, buffer_weights_bwd
        integer, dimension(mesh%get_n_points_boundary()) :: bnd_descrs_nmn
        real(GP), allocatable, dimension(:) :: pen_val_halo_fwd, pen_val_halo_bwd
        
        ! Set debug output level
        outi = 0
        if ( present(dbgout) ) then
            if (is_rank_info_writer) then
                outi = dbgout
            endif
            if ( dbgout >= 3 ) then
                outi = dbgout ! Every rank prints debug info
            endif
        endif

        ! Set grid type
        self%is_staggered = staggered

        ! Extrapolate penalisation values to full mesh
        !$omp parallel private(l, ki, kb)
        !$omp do
        do ki = 1, mesh%get_n_points_inner()
            l = mesh%inner_indices(ki)
            pen_vals(l) = penalisation%get_charfun_val(ki)
        enddo
        !$omp end do
        !$omp do
        do kb = 1, mesh%get_n_points_boundary()
            bnd_descrs_nmn(kb) = BND_TYPE_NEUMANN
        enddo
        !$omp end do
        call set_perpbnds(mesh, bnd_descrs_nmn, pen_vals)
        call extrapolate_ghost_points(mesh, pen_vals)

        ! Fill masked array
        !$omp do
        do l = 1, mesh%get_n_points()
            if ( pen_vals(l) > GP_EPS ) then
                pen_val_mask(l) = 1.0_GP
            else
                pen_val_mask(l) = 0.0_GP
            endif
        enddo
        !$omp end do
        !$omp end parallel
        
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), +1,  mesh%get_n_points(), &
                                 pen_val_mask, nrecv, pen_val_halo_fwd)
        
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1,  mesh%get_n_points(), &
                                 pen_val_mask, nrecv, pen_val_halo_bwd)
        

        ! Running tallies of points found
        n_fwd = 0
        n_bwd = 0
        
        do ki = 1, mesh%get_n_points_inner()
            l = mesh%inner_indices(ki)
            ! Assert that origin point is not in the penalisation volume
            if ( penalisation%get_charfun_val(ki) < GP_EPS ) then
                ! Map from forward plane
                if (.not. self%is_staggered) then
                    pen_val_fwd = map%upstream_cano_from_cano_fwd_single(&
                                        pen_val_halo_fwd, l)
                else
                    pen_val_fwd = map%upstream_stag_from_stag_fwd_single(&
                                        pen_val_halo_fwd, l)
                endif
                if ( pen_val_fwd > pen_threshold ) then
                    n_fwd = n_fwd + 1
                    buffer_inds_fwd(n_fwd) = ki
                    buffer_weights_fwd(n_fwd) = pen_val_fwd
                endif
                ! Map from backward plane
                if (.not. self%is_staggered) then
                    pen_val_bwd = map%upstream_cano_from_cano_bwd_single(&
                                        pen_val_halo_bwd, l)
                else
                    pen_val_bwd = map%upstream_stag_from_stag_bwd_single(&
                                        pen_val_halo_bwd, l)
                endif
                if ( pen_val_bwd > pen_threshold ) then
                    n_bwd = n_bwd + 1
                    buffer_inds_bwd(n_bwd) = ki
                    buffer_weights_bwd(n_bwd) = pen_val_bwd
                endif
            endif
        enddo
        
        deallocate(pen_val_halo_bwd)
        deallocate(pen_val_halo_fwd)

        self%n_inds_fwd = n_fwd
        self%n_inds_bwd = n_bwd

        if ( outi > 0 ) then
            write(get_stdout(),*)"Parallel target flux markers built:"
            write(get_stdout(),*)"  On staggered grid:     ", self%is_staggered
            write(get_stdout(),*)"  Forward plane points:  ", n_fwd
            write(get_stdout(),*)"  Backward plane points: ", n_bwd
            write(get_stdout(),*)
        endif

        ! Allocate and fill final arrays
        if (.not. allocated(self%inds_fwd)) allocate(self%inds_fwd(n_fwd))
        if (.not. allocated(self%inds_bwd)) allocate(self%inds_bwd(n_bwd))
        if (.not. allocated(self%weights_fwd)) allocate(self%weights_fwd(n_fwd))
        if (.not. allocated(self%weights_bwd)) allocate(self%weights_bwd(n_bwd))

        ! Forward plane
        !$omp parallel private(i)
        !$omp do
        do i = 1, n_fwd
            self%inds_fwd(i) = buffer_inds_fwd(i) 
            self%weights_fwd(i) = buffer_weights_fwd(i) 
        enddo
        !$omp end do
        ! Backward plane
        !$omp do
        do i = 1, n_bwd
            self%inds_bwd(i) = buffer_inds_bwd(i) 
            self%weights_bwd(i) = buffer_weights_bwd(i) 
        enddo
        !$omp end do
        !$omp end parallel

    end subroutine

    module subroutine destructor(self)
        type(parallel_target_flux_t), intent(inout) :: self

        if (allocated(self%inds_fwd)) deallocate(self%inds_fwd)
        if (allocated(self%inds_bwd)) deallocate(self%inds_bwd)
        if (allocated(self%weights_fwd)) deallocate(self%weights_fwd)
        if (allocated(self%weights_bwd)) deallocate(self%weights_bwd)
       
    end subroutine

end submodule