component_perptransp_s.f90 Source File

submodule(timestep_template_m) component_perptransp_s
    !! Implementation of perpendicular (diffusive) transport model
    !! \[ \partial_t n \rightarrow \partial_t n + \nabla\cdot\left(D_n\nabla_\perp n\right) \]
    !! \[ \partial_t \Gamma_\parallel \rightarrow \partial_t \Gamma_\parallel + \nabla\cdot\left(D_{\Gamma_\parallel}\nabla_\perp \Gamma_\parallel\right) \]
    use params_template_m, only : dperp_coeff_dens, dperp_coeff_parmom, dperp_coeff_pion
    implicit none

contains

    module subroutine apply_perptransp(comm_handler)
        type(comm_handler_t), intent(in) :: comm_handler

        integer :: l, lnorth, lsouth, least, lwest
        real(GP) :: lapl_perp_dens, lapl_perp_pion, hf_inv_sqrd
        real(GP), parameter :: two_thirds = 2.0_GP / 3.0_GP

        hf_inv_sqrd = 1.0_GP / mesh_cano%get_spacing_f()**2

        !$omp parallel default(none) &
        !$omp private(l, lnorth, lsouth, least, lwest, lapl_perp_dens, lapl_perp_pion) &
        !$omp shared(np_max, hf_inv_sqrd, mesh_cano, masks_cano,  mesh_stag, masks_stag, &
        !$omp        dens, ddens, dperp_coeff_dens, parmom, dparmom, dperp_coeff_parmom, &
        !$omp        pion, dpion, dperp_coeff_pion)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_inner(l)) then
                lnorth = mesh_cano%get_index_neighbor(0, 1, l)
                lsouth = mesh_cano%get_index_neighbor(0, -1, l)
                least = mesh_cano%get_index_neighbor(1, 0, l)
                lwest = mesh_cano%get_index_neighbor(-1, 0, l)

                lapl_perp_dens = hf_inv_sqrd * (dens(lnorth) + dens(lsouth) + dens(least) + dens(lwest) - 4.0_GP * dens(l))
                lapl_perp_pion = hf_inv_sqrd * (pion(lnorth) + pion(lsouth) + pion(least) + pion(lwest) - 4.0_GP * pion(l))

                ddens(l) = ddens(l) + dperp_coeff_dens * lapl_perp_dens
                dpion(l) = dpion(l) + two_thirds * dperp_coeff_pion * lapl_perp_pion
            endif

            if (masks_stag%is_inner(l)) then
                lnorth = mesh_stag%get_index_neighbor(0, 1, l)
                lsouth = mesh_stag%get_index_neighbor(0, -1, l)
                least = mesh_stag%get_index_neighbor(1, 0, l)
                lwest = mesh_stag%get_index_neighbor(-1, 0, l)

                lapl_perp_dens = hf_inv_sqrd * (parmom(lnorth) + parmom(lsouth) + parmom(least) + parmom(lwest) - 4.0_GP * parmom(l))
                dparmom(l) = dparmom(l) + dperp_coeff_parmom * lapl_perp_dens
            endif

        enddo
        !$omp end do
        !$omp end parallel

    end subroutine

end submodule