diff_coeff_neutrals_m.f90 Source File

module diff_coeff_neutrals_m
!! Module for calculating the diffusion coefficient D_N
    use precision_grillix_m, only : GP, GP_NAN, GP_EPS
    use parallel_map_m, only : parallel_map_t
    use variable_m, only: variable_t
    ! From PARALLAX
    use mesh_cart_m, only: mesh_cart_t
    use equilibrium_m, only : equilibrium_t
    ! Parameters
    use params_brag_model_m, only : &
        tratio
    use params_neut_model_m, only : &
        s_cx, mach_limit, diff_tn_min
    use params_neut_floors_m, only : &
        floor_temp
    use params_neut_switches_m, only : &
        swn_dens_diff_tor
    implicit none
    
    public :: eval_diff_coeff
 
contains

    subroutine eval_diff_coeff(equi, mesh_cano, map, ne, ti, neutrals_temp, neutrals_pressure, diff_coeff)
        !! Evaluate the normalized pressure diffusion coefficient
        class(equilibrium_t), intent(in) :: equi
        !! Equilibrium
        type(mesh_cart_t), intent(in) :: mesh_cano
        !! Mesh (canonical)
        type(parallel_map_t), intent(in) :: map
        !! Parallel map
        type(variable_t), intent(in) :: ne
        !! Plasma density
        type(variable_t), intent(in) :: ti
        !! Ion temperature
        type(variable_t), intent(in) :: neutrals_temp
        !! Neutrals temperature
        type(variable_t), intent(in) :: neutrals_pressure
        !! Neutrals pressure N * TN
        real(GP), intent(out), dimension(mesh_cano%get_n_points()) :: diff_coeff
        !! Diffusion coefficient values

        ! Intermediary variable for | grad( N * TN ) |
        real(GP) :: absgrad
        ! Intermediary variable for ( N * c_s,N / absgrad )
        real(GP) :: diff_limit
        ! Mesh indices
        integer :: l, ki, kb, kg, lwest, least, lsouth, lnorth
        ! Derivatives of neutrals pressure ( N*TN ) in x, y, z direction
        real(GP) :: ddx_neutrals_pressure, ddy_neutrals_pressure, ddz_neutrals_pressure
        ! Mesh spacing and coordinates
        real(GP) :: hf, prefac, phi, x, y

        hf = mesh_cano%get_spacing_f()
        prefac = 1.0_GP/(2.0_GP*hf)
        phi = mesh_cano%get_phi()

        !$omp parallel default(none) &
        !$omp          private(l, ki, kb, kg, x, y, lwest, least, lsouth, lnorth, absgrad, diff_limit, &
        !$omp                  ddx_neutrals_pressure, ddy_neutrals_pressure, ddz_neutrals_pressure) &
        !$omp          shared(equi, mesh_cano, map, ne, ti, neutrals_temp, neutrals_pressure, diff_coeff, &
        !$omp                 diff_tn_min, hf, prefac, phi, swn_dens_diff_tor, mach_limit, tratio, s_cx)
        !$omp do
        do ki = 1, mesh_cano%get_n_points_inner()
            l = mesh_cano%inner_indices(ki)
            ! Compute unlimited diffusion coefficient
            diff_coeff(l) = max(1.0_GP, diff_tn_min/neutrals_temp%vals(l)) * sqrt(tratio) &
                            / ( 2.93_GP * s_cx * ne%vals(l) * sqrt(ti%vals(l)) )

            ! Limit diffusivity, can only apply on inner points (required for gradient)
            !   Can disable limit by setting parameter mach_limit <= 0.0
            if ( mach_limit > GP_EPS ) then

                x = mesh_cano%get_x(l)
                y = mesh_cano%get_y(l)

                ! Compute neighbor indices
                lwest   = mesh_cano%get_index_neighbor(-1,  0, l)
                least   = mesh_cano%get_index_neighbor( 1,  0, l)
                lsouth  = mesh_cano%get_index_neighbor( 0, -1, l)
                lnorth  = mesh_cano%get_index_neighbor( 0,  1, l)

                ! Compute derivatives
                ddx_neutrals_pressure   = prefac * (neutrals_pressure%vals(least) - neutrals_pressure%vals(lwest))
                ddy_neutrals_pressure   = prefac * (neutrals_pressure%vals(lnorth) - neutrals_pressure%vals(lsouth))
                if (swn_dens_diff_tor > GP_EPS) then
                    ddz_neutrals_pressure = (neutrals_pressure%hfwd(l) - neutrals_pressure%hbwd(l)) &
                                    / (2.0_GP * equi%jacobian(x, y, phi) * map%get_dphi())
                else
                    ddz_neutrals_pressure = 0.0_GP
                end if

                absgrad = sqrt( ddx_neutrals_pressure**2 + ddy_neutrals_pressure**2 + ddz_neutrals_pressure**2 )

                if ( absgrad > GP_EPS ) then
                    ! Compute diffusion coefficient limit
                    diff_limit = mach_limit * neutrals_pressure%vals(l) * sqrt( tratio / neutrals_temp%vals(l)) / absgrad

                    diff_coeff(l) = min( diff_coeff(l), diff_limit )
                endif
            endif
        enddo
        !$omp end do
        !$omp do
        do kb = 1, mesh_cano%get_n_points_boundary()
            l = mesh_cano%boundary_indices(kb)
            diff_coeff(l) = max(1.0_GP, diff_tn_min/neutrals_temp%vals(l)) * sqrt(tratio) &
                            / ( 2.93_GP * s_cx * ne%vals(l) * sqrt(ti%vals(l)) )
        enddo
        !$omp end do
        !$omp do
        do kg = 1, mesh_cano%get_n_points_ghost()
            l = mesh_cano%ghost_indices(kg)
            diff_coeff(l) = max(1.0_GP, diff_tn_min/neutrals_temp%vals(l)) * sqrt(tratio) &
                            / ( 2.93_GP * s_cx * ne%vals(l) * sqrt(ti%vals(l)) )
        enddo
        !$omp end do
        !$omp end parallel

    end subroutine

end module