parallel_target_flux_compute_s.f90 Source File

submodule (parallel_target_flux_m) parallel_target_flux_compute_s
    !! Routines for computing parallel fluxes into target
    implicit none

contains

    module subroutine compute_target_flux(self, comm_handler, mesh, penalisation, equi_storage, &
                                          g, flux_fwd, flux_bwd, mode)
        class(parallel_target_flux_t), intent(in) :: self
        type(comm_handler_t), intent(in) :: comm_handler
        type(mesh_cart_t), intent(in) :: mesh
        type(penalisation_t), intent(in) :: penalisation
        type(equilibrium_storage_t), intent(in) :: equi_storage
        real(GP), dimension(mesh%get_n_points()), intent(in) :: g
        real(GP), intent(out) :: flux_fwd
        real(GP), intent(out) :: flux_bwd
        character(*), intent(in), optional :: mode

        integer i, ki, l, ierr
        real(GP) :: flux_fwd_plane, flux_bwd_plane

        flux_fwd_plane = 0.0_GP
        flux_bwd_plane = 0.0_GP

        if ( trim(mode) == 'DEFAULT' .or. (.not. present(mode))) then
            !$omp parallel do private(i, l, ki) reduction(+: flux_fwd_plane)
            do i = 1, self%n_inds_fwd
                ki = self%inds_fwd(i)
                l = mesh%inner_indices(ki)
                flux_fwd_plane = flux_fwd_plane + g(l) &
                                 * penalisation%get_dirindfun_val(ki) &
                                 * equi_storage%btor(l) / equi_storage%absb(l)
            enddo
            !$omp end parallel do
            !$omp parallel do private(i, l, ki) reduction(+: flux_bwd_plane)
            do i = 1, self%n_inds_bwd
                ki = self%inds_bwd(i)
                l = mesh%inner_indices(ki)
                flux_bwd_plane = flux_bwd_plane + g(l) &
                                 * penalisation%get_dirindfun_val(ki) &
                                 * equi_storage%btor(l) / equi_storage%absb(l)
            enddo
            !$omp end parallel do
        else if ( trim(mode) == 'WEIGHTED' ) then
            !$omp parallel do private(i, l, ki) reduction(+: flux_fwd_plane)
            do i = 1, self%n_inds_fwd
                ki = self%inds_fwd(i)
                l = mesh%inner_indices(ki)
                flux_fwd_plane = flux_fwd_plane + g(l) * self%weights_fwd(i) &
                                 * penalisation%get_dirindfun_val(ki) &
                                 * equi_storage%btor(l) / equi_storage%absb(l)
            enddo
            !$omp end parallel do
            !$omp parallel do private(i, l, ki) reduction(+: flux_bwd_plane)
            do i = 1, self%n_inds_bwd
                ki = self%inds_bwd(i)
                l = mesh%inner_indices(ki)
                flux_bwd_plane = flux_bwd_plane + g(l) * self%weights_bwd(i) &
                                 * penalisation%get_dirindfun_val(ki) &
                                 * equi_storage%btor(l) / equi_storage%absb(l)
            enddo
            !$omp end parallel do
        else
            call handle_error("Invalid mode, must be ['DEFAULT','WEIGHTED']", &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                              error_info_t('mode: '//trim(mode)))
        endif

        flux_fwd_plane = flux_fwd_plane * mesh%get_spacing_f()**2 
        flux_bwd_plane = flux_bwd_plane * mesh%get_spacing_f()**2

        call MPI_allreduce(flux_fwd_plane, flux_fwd, 1, MPI_GP, MPI_SUM, &
                           comm_handler%get_comm_planes(), ierr)
        call MPI_allreduce(flux_bwd_plane, flux_bwd, 1, MPI_GP, MPI_SUM, &
                           comm_handler%get_comm_planes(), ierr)

    end subroutine

    module function compute_target_flux_single(self, mesh, penalisation, equi_storage, &
                                               i_marker, g, dirind, mode) result(flux_out)
        class(parallel_target_flux_t), intent(in) :: self
        type(mesh_cart_t), intent(in) :: mesh
        type(penalisation_t), intent(in) :: penalisation
        type(equilibrium_storage_t), intent(in) :: equi_storage
        integer, intent(in) :: i_marker
        real(GP), dimension(mesh%get_n_points()), intent(in) :: g
        character(*), intent(in) :: dirind
        character(*), intent(in), optional :: mode

        real(GP) :: flux_out, weight
        integer :: ki, l

        if (trim(dirind) == 'fwd') then
            ki = self%inds_fwd(i_marker)
            weight = self%weights_fwd(i_marker)
        else if (trim(dirind) == 'bwd') then
            ki = self%inds_bwd(i_marker)
            weight = self%weights_bwd(i_marker)
        else
            call handle_error("Invalid dirind, must be ['fwd','bwd']", &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                              error_info_t('dirind: '//trim(dirind)))
        endif

        l = mesh%inner_indices(ki)
        flux_out = 0.0_GP

        if ( trim(mode) == 'DEFAULT' .or. (.not. present(mode))) then
            flux_out = g(l) * penalisation%get_dirindfun_val(ki) &
                       * equi_storage%btor(l)/equi_storage%absb(l)

        else if ( trim(mode) == 'WEIGHTED' ) then
            flux_out = g(l) * weight * penalisation%get_dirindfun_val(ki)  &
                       * equi_storage%btor(l)/equi_storage%absb(l)
        else
            call handle_error("Invalid mode, must be ['DEFAULT','WEIGHTED']", &
                              GRILLIX_ERR_OTHER, __LINE__, __FILE__, &
                              error_info_t('mode: '//trim(mode)))
        endif

        flux_out = flux_out * mesh%get_spacing_f()**2 

    end function

end submodule