immersed_template_m.f90 Source File

module immersed_template_m
    !! Routines related with parallel (Sheath) boundaries treated with immersion method
    use precision_grillix_m, only : GP
    use error_handling_grillix_m, only: handle_error, error_info_t
    use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST, GRILLIX_ERR_OTHER
    use comm_handler_m, only : comm_handler_t
    use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane
    use static_data_m, only : mesh_cano, masks_cano, mesh_stag, masks_stag, map, &
        parbnd_immersed_cano, parbnd_immersed_stag
    use params_template_m, only : parbnd_type_dens, parbnd_type_parmom, parbnd_immersed_epsinv
    use state_template_m, only : np_max, dens, parmom
    implicit none

    public :: apply_immersed_sources

    real(GP), allocatable, dimension(:) :: dirind_cano, chi_cano, dirind_stag, chi_stag
    ! Penalisation function in global indexing context: Needed as a temporary workaround
contains

    subroutine apply_immersed_sources(comm_handler, ddens, dparmom, dpion)
        !! Applies immersed sources, related with parallel (Sheath) boundary conditions
        type(comm_handler_t), intent(in) :: comm_handler
        !! Communication handler
        real(GP), intent(inout), dimension(np_max) :: ddens
        !! Changerate for density
        real(GP), intent(inout), dimension(np_max) :: dparmom
        !! Changerate for parallel momentum
        real(GP), intent(inout), dimension(np_max) :: dpion
        !! Changerate for ion pressure

        logical, save :: first_call = .true.
        integer :: ki, l, nrecv

        real(GP), dimension(np_max) :: dens_target, parmom_target, wrk_fwd, wrk_bwd
        real(GP), allocatable, dimension(:) :: halo_buf

        if (first_call) then
            first_call = .false.
            allocate(dirind_cano(np_max), chi_cano(np_max), dirind_stag(np_max), chi_stag(np_max))
            ! Temporary workaround: Map penalisation fucntions from old loop convention (ki) to new loop convention
            ! Shall be resolved, when switching to PARALLAX based penalisation objects

            !$omp parallel default(none) private(l) &
            !$omp shared(np_max, dirind_cano, chi_cano, dirind_stag, chi_stag) 
            !$omp do
            do l = 1, np_max
                dirind_cano(l) = 0.0_GP
                chi_cano(l) = 0.0_GP
                dirind_stag(l) = 0.0_GP
                chi_stag(l) = 0.0_GP
            enddo
            !$omp end do
            !$omp end parallel
            do ki = 1, mesh_cano%get_n_points_inner()
                l = mesh_cano%inner_indices(ki)
                dirind_cano(l) = parbnd_immersed_cano%dirindfun(ki)
                chi_cano(l) = parbnd_immersed_cano%charfun(ki)
            enddo
            do ki = 1, mesh_stag%get_n_points_inner()
                l = mesh_stag%inner_indices(ki)
                dirind_stag(l) = parbnd_immersed_stag%dirindfun(ki)
                chi_stag(l) = parbnd_immersed_stag%charfun(ki)
            enddo
        endif

        ! Get target values for fields
        select case(parbnd_type_dens)
            case('Neumann')
                ! Get upstream values
                call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, dens, nrecv, halo_buf)
                call map%upstream_cano_from_cano_fwd(halo_buf, wrk_fwd)
                call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, dens, nrecv, halo_buf)
                call map%upstream_cano_from_cano_bwd(halo_buf, wrk_bwd)

                !$omp parallel default(none) private(l) &
                !$omp shared(np_max, masks_cano, dirind_cano, dens_target, wrk_fwd, wrk_bwd)
                !$omp do
                do l = 1, np_max
                    if (masks_cano%is_inner(l)) then
                        if (dirind_cano(l) > 0.0_GP) then
                            dens_target(l) = wrk_bwd(l)
                        else
                            dens_target(l) = wrk_fwd(l)
                        endif
                    endif
                enddo
                !$omp end do
                !$omp end parallel

            case default
                call handle_error('parbnd_type_dens not valid or available with immersed method: ' &
                    // parbnd_type_dens, GRILLIX_ERR_OTHER, __LINE__, __FILE__)
        end select

        select case(parbnd_type_parmom)
            case('Sonic')
                ! Compute parallel momentum based on sound speed
                parmom_target = 0.0_GP
                call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, dens, nrecv, halo_buf)
                call map%lift_cano_to_stag(dens, halo_buf, parmom_target)

                !$omp parallel default(none) private(l) &
                !$omp shared(np_max, masks_stag, dirind_stag, parmom_target)
                !$omp do
                do l = 1, np_max
                    if (masks_stag%is_inner(l)) then
                            parmom_target(l) = dirind_stag(l) * parmom_target(l)
                    endif
                enddo
                !$omp end do
                !$omp end parallel
         case default
                call handle_error('parbnd_type_parmom not valid or available with immersed method: ' &
                    // parbnd_type_parmom, GRILLIX_ERR_OTHER, __LINE__, __FILE__)
        end select

        !$omp parallel default(none) private(l) &
        !$omp shared(np_max, masks_cano, masks_stag, chi_cano, chi_stag, parbnd_immersed_epsinv, &
        !$omp        dens, parmom, ddens, dparmom, dens_target, parmom_target)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_inner(l)) then
                ddens(l) = (1.0_GP - chi_cano(l)) * ddens(l) &
                           + chi_cano(l) * parbnd_immersed_epsinv * (dens_target(l) - dens(l))
            endif
            if (masks_stag%is_inner(l)) then
                dparmom(l) = (1.0_GP - chi_stag(l)) * dparmom(l) &
                             + chi_stag(l) * parbnd_immersed_epsinv * (parmom_target(l) - parmom(l))
            endif
        enddo
        !$omp end do
        !$omp end parallel

        if (allocated(halo_buf)) deallocate(halo_buf)

    end subroutine

end module