component_partransp_s.f90 Source File

submodule(timestep_template_m) component_partransp_s
    !! Implementation of paralllel transport model
    !! \[ \partial_t n \rightarrow \partial_t n - \nabla\cdot\left(\mathbf{b}\Gamma_\parallel\right) \]
    !! \[ \partial_t \Gamma_\parallel \rightarrow \partial_t \Gamma_\parallel - \nabla\cdot\left(\mathbf{b}\Gamma_\parallel u_\parallel\right) - \nabla_\parallel p_i \]
    !! \[ \partial_t p_i \rightarrow \partial_t p_i - \nabla\cdot\left(\mathbf{b}p_i u_\parallel\right) - \frac23p_i\nabla_\parallel u_\parallel  + \frac23\nabla\cdot\left(\mathbf{b}T_i^{5/2}\nabla_\parallel T_i\right)\]
    use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane 
    use params_template_m, only : ion_heat_cond_coeff
    implicit none

contains

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

        integer :: l, nrecv
        real(GP), parameter :: two_thirds = 2.0_GP / 3.0_GP
        real(GP), allocatable, dimension(:) :: halo_buf
        real(GP), dimension(np_max) :: ddt_wrk1, ddt_wrk2, ddt_wrk3, ddt_wrk4, pion_stag, tion_stag, parflux, &
        parpionflux, pariheatcond, divb_parflux, divb_upar, gradb_tion

        ! ddt_wrk1 = divb(parmom)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, parmom, nrecv, halo_buf)
        call map%par_divb(parmom, halo_buf, ddt_wrk1)

        ! divb_upar = divb(upar)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, upar, nrecv, halo_buf)
        call map%par_divb(upar, halo_buf, divb_upar)

        ! ddt_wrk2 = gradb(pion)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, pion, nrecv, halo_buf)
        call map%par_grad(pion, halo_buf, ddt_wrk2)
        ! map pion to staggered mesh -> pion_stag
        call map%lift_cano_to_stag(pion, halo_buf, pion_stag)

        ! gradb_tion = gradb(tion)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, tion, nrecv, halo_buf)
        call map%par_grad(tion, halo_buf, gradb_tion)
        ! map tion to staggered mesh -> tion_stag
        call map%lift_cano_to_stag(tion, halo_buf, tion_stag)

        ! ddens += - ddt_wrk1
        ! dparmom += - ddt_wrk2
        ! dpion += - 2/3 * pion * divb_upar
        ! parflux = parmom * upar
        ! parpionflux = pion_stag * upar
        ! pariheatcond = tion_stag**(5/2) * gradb_tion
        !$omp parallel default(none) private(l) &
        !$omp shared(np_max, masks_cano, masks_stag, ddens, dparmom, dpion, ddt_wrk1, ddt_wrk2, divb_upar, &
        !$omp        parflux, parpionflux, pariheatcond, parmom, upar, pion, pion_stag, tion_stag, gradb_tion)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_inner(l)) then
                ddens(l) = ddens(l) - ddt_wrk1(l)
                dpion(l) = dpion(l) - two_thirds * pion(l) * divb_upar(l)
            endif
            if (masks_stag%is_inner(l)) then
                dparmom(l) = dparmom(l) - ddt_wrk2(l)
            endif
            if (masks_stag%is_inner(l) .or. masks_stag%is_boundary(l) .or. masks_stag%is_ghost(l)) then
                parflux(l) = parmom(l) * upar(l)
                parpionflux(l) = pion_stag(l) * upar(l)
                pariheatcond(l) = sqrt(tion_stag(l)**5) * gradb_tion(l)
            endif
        enddo
        !$omp end do
        !$omp end parallel

        ! divb(parflux)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, parflux, nrecv, halo_buf)
        call map%par_divb(parflux, halo_buf, divb_parflux)
        ! map this to staggered mesh -> ddt_wrk1
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, divb_parflux, nrecv, halo_buf)
        call map%lift_cano_to_stag(divb_parflux, halo_buf, ddt_wrk1)

        ! ddt_wrk3 = divb(parpionflux)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, parpionflux, nrecv, halo_buf)
        call map%par_divb(parpionflux, halo_buf, ddt_wrk3)

        ! ddt_wrk4 = divb(pariheatcond)
        call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, pariheatcond, nrecv, halo_buf)
        call map%par_divb(pariheatcond, halo_buf, ddt_wrk4)

        ! dpion += - ddt_wrk3 + 2/3 * ion_heat_cond_coeff * ddt_wrk4
        ! dparmom += -ddt_wrk1
        
        !$omp parallel default(none) private(l) shared(np_max, masks_cano, masks_stag, dpion, dparmom, &
        !$omp                                          ddt_wrk1, ddt_wrk3, ddt_wrk4, ion_heat_cond_coeff)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_inner(l)) then
                dpion(l) = dpion(l) - ddt_wrk3(l) + two_thirds * ion_heat_cond_coeff * ddt_wrk4(l)
            end if
            if (masks_stag%is_inner(l)) then
                dparmom(l) = dparmom(l) - ddt_wrk1(l)
            endif
        enddo
       !$omp end do
       !$omp end parallel

    end subroutine

end submodule