component_pardiff_s.f90 Source File

submodule(timestep_template_m) component_pardiff_s
    !! Implementation of paralllel diffusion terms
    !! \[ \partial_t n \rightarrow \partial_t n + \nabla\cdot\left(D_n \mathbf{b}\nabla_parallel n\right) \]
    !! \[ \partial_t \Gamma_\parallel \rightarrow \partial_t \Gamma_\parallel + \nabla\cdot\left(D_{\Gamma_{\parallel}}\mathbf{b}\nabla_parallel \Gamma_\parallel\right) \]
    use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane
    use params_template_m, only : dpar_coeff_dens, dpar_coeff_parmom, dpar_coeff_pion
    implicit none

contains

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

        integer :: nrecv, l
        real(GP), allocatable, dimension(:) :: halo_buf
        real(GP), dimension(np_max) :: parflux, pdiff_dens, pdiff_parmom, pdiff_pion
    
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, dens, nrecv, halo_buf)
        call map%par_grad(dens, halo_buf, parflux)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, parflux, nrecv, halo_buf)
        call map%par_divb(parflux, halo_buf, pdiff_dens)

        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, pion, nrecv, halo_buf)
        call map%par_grad(pion, halo_buf, parflux)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, parflux, nrecv, halo_buf)
        call map%par_divb(parflux, halo_buf, pdiff_pion)

        call getdata_fwdbwdplane(comm_handler%get_comm_planes(),-1, np_max, parmom, nrecv, halo_buf)
        call map%par_divb(parmom, halo_buf, parflux)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, parflux, nrecv, halo_buf)
        call map%par_grad(parflux, halo_buf, pdiff_parmom)
        
        !$omp parallel default(none) private(l) &
        !$omp shared(np_max, masks_cano, masks_stag, dpar_coeff_dens, dpar_coeff_parmom, dpar_coeff_pion, &
        !$omp        ddens, dparmom, dpion, pdiff_dens, pdiff_parmom, pdiff_pion)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_inner(l)) then
                ddens(l) = ddens(l) + dpar_coeff_dens * pdiff_dens(l)
                dpion(l) = dpion(l) + dpar_coeff_pion * pdiff_pion(l)
            endif
            if (masks_stag%is_inner(l)) then
                dparmom(l) = dparmom(l) + dpar_coeff_parmom * pdiff_parmom(l)
            endif
        enddo
        !$omp end do
        !$omp end parallel

        if (allocated(halo_buf)) deallocate(halo_buf)

    end subroutine

end submodule