component_mms_sources_s.f90 Source File

submodule(timestep_template_m) component_mms_sources_s
    !! Application of MMS sources and boundary conditions
    !! \[ \partial_t n \rightarrow \partial_t n + S^{MMS}_{n}
    !! \[ \partial_t \Gamma_\parallel \rightarrow \partial_t \Gamma_\parallel +  S^{MMS}_{\Gamma_{\parallel}}) \]
    use mms_template_m, only: mms_sol_dens_template, mms_source_dens_template, mms_sol_parmom_template, mms_source_parmom_template, &
        mms_sol_pion_template, mms_source_pion_template
    implicit none

contains

    module subroutine apply_mms_sources(comm_handler, tau)
        type(comm_handler_t), intent(in) :: comm_handler
        real(GP), intent(in) :: tau

        integer :: l
        real(GP) :: x, y, phi_cano, phi_stag

        phi_cano = mesh_cano%get_phi()
        phi_stag = mesh_stag%get_phi()

        !$omp parallel default(none) private(l, x, y) &
        !$omp shared(np_max, tau, mesh_cano, masks_cano, phi_cano, mesh_stag, masks_stag, phi_stag, ddens, dparmom, &
        !$omp        dpion, mms_source_dens_template, mms_source_parmom_template, mms_source_pion_template)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_inner(l)) then
                x = mesh_cano%get_x(l)
                y = mesh_cano%get_y(l)
                ddens(l) = ddens(l) + mms_source_dens_template(x, y, phi_cano, tau)
                dpion(l) = dpion(l) + mms_source_pion_template(x, y, phi_cano, tau)
            endif
            if (masks_stag%is_inner(l)) then
                x = mesh_stag%get_x(l)
                y = mesh_stag%get_y(l)
                dparmom(l) = dparmom(l) + mms_source_parmom_template(x, y, phi_stag, tau)
            endif
        enddo
        !$omp end do
        !$omp end parallel

    end subroutine

    module subroutine apply_mms_boundaries(comm_handler, tau)
        type(comm_handler_t), intent(in) :: comm_handler
        real(GP), intent(in) :: tau

        integer :: l
        real(GP) :: tau_adv, x, y, phi_cano, phi_stag

        phi_cano = mesh_cano%get_phi()
        phi_stag = mesh_stag%get_phi()
        
        !$omp parallel default(none) private(l, x, y) &
        !$omp shared(np_max, mesh_cano, masks_cano, phi_cano, mesh_stag, masks_stag, phi_stag, &
        !$omp        tau, dens, parmom, pion, &
        !$omp        mms_sol_dens_template, mms_sol_parmom_template, mms_sol_pion_template)
        !$omp do
        do l = 1, np_max
            if (masks_cano%is_boundary(l) .or. masks_cano%is_ghost(l)) then
                x = mesh_cano%get_x(l)
                y = mesh_cano%get_y(l)
                dens(l) = mms_sol_dens_template(x, y, phi_cano, tau)
                pion(l) = mms_sol_pion_template(x, y, phi_cano, tau)
            endif
            if (masks_stag%is_boundary(l) .or. masks_stag%is_ghost(l)) then
                x = mesh_stag%get_x(l)
                y = mesh_stag%get_y(l)
                parmom(l) = mms_sol_parmom_template(x, y, phi_stag, tau)
            endif
        enddo
        !$omp end do
        !$omp end parallel

    end subroutine

end submodule