module taylor_template_m !! Routines related with parallel (Sheath) boundaries treated with Taylor's method use precision_grillix_m, only : GP use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use mpi_mapping_auxiliaries_m, only : getdata_fwdbwdplane use comm_handler_m, only : comm_handler_t use params_template_m, only : parbnd_taylor_order, parbnd_type_dens, parbnd_type_parmom use static_data_m, only : equi, mesh_cano, masks_cano, map, parbnd_taylor_cano, parbnd_taylor_stag use state_template_m, only : np_max, dens, parmom implicit none public :: apply_taylor_parbnds integer, private :: ndepth_max !! Maximum depth of guard points contains subroutine apply_taylor_parbnds(comm_handler) !! Applies parallel (Sheath) boundary conditions with Taylor's method type(comm_handler_t), intent(in) :: comm_handler !! Communication handler logical, save :: first_call = .true. integer :: k, nrecv, l real(GP) :: depth_max_cano, depth_max_stag real(GP), allocatable, dimension(:) :: halo_buf real(GP), dimension(np_max) :: wrk_zeros, wrk_bwd2, wrk_bwd1, wrk_fwd1, wrk_fwd2, ux, & parmom_cano, parmom_sonic, parmom_extrapolate, psec_sonic, psec_extrapolate ! Determine ndepth_max if (first_call) then first_call = .false. depth_max_cano = maxval(abs(parbnd_taylor_cano%depth_in_bnd)) depth_max_stag = maxval(abs(parbnd_taylor_stag%depth_in_bnd)) ndepth_max = ceiling(max(depth_max_cano, depth_max_stag)) endif ! Set parallel boundary conditions on density do k = 1, ndepth_max call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -2, np_max, dens, nrecv, halo_buf) call map%upstream_cano_from_cano_bwd2(halo_buf, wrk_bwd2) call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, dens, nrecv, halo_buf) call map%upstream_cano_from_cano_bwd(halo_buf, wrk_bwd1) call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, dens, nrecv, halo_buf) call map%upstream_cano_from_cano_fwd(halo_buf, wrk_fwd1) call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 2, np_max, dens, nrecv, halo_buf) call map%upstream_cano_from_cano_fwd2(halo_buf, wrk_fwd2) select case(parbnd_type_dens) case('Neumann') wrk_zeros = 0.0_GP call parbnd_taylor_cano%set_guard_points(equi, mesh_cano, 0.0_GP, 1.0_GP, wrk_zeros, & dens, wrk_bwd2, wrk_bwd1, wrk_fwd1, wrk_fwd2, parbnd_taylor_order, ux) case('Extrapolate') call parbnd_taylor_cano%extrapolate_guard_points(equi, mesh_cano, & dens, wrk_bwd2, wrk_bwd1, wrk_fwd1, wrk_fwd2, ux) case default call handle_error('parbnd_type_dens not valid or available with Taylor method: ' & // parbnd_type_dens, GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select enddo ! Target values for parallel momentum !$omp parallel default(none) private(l) shared(np_max, parbnd_taylor_cano, ux) !$omp do do l = 1, np_max ux(l) = ux(l) * parbnd_taylor_cano%b_bnd_direction(l) enddo !$omp end do !$omp end parallel ! Map parallel momentum to canonical mesh call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, parmom, nrecv, halo_buf) call map%flat_stag_to_cano(parmom, halo_buf, parmom_cano) ! Set boundary conditions on parallel momentum on canonical mesh do k = 1, ndepth_max call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -2, np_max, parmom_cano, nrecv, halo_buf) call map%upstream_cano_from_cano_bwd2(halo_buf, wrk_bwd2) call getdata_fwdbwdplane(comm_handler%get_comm_planes(), -1, np_max, parmom_cano, nrecv, halo_buf) call map%upstream_cano_from_cano_bwd(halo_buf, wrk_bwd1) call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, parmom_cano, nrecv, halo_buf) call map%upstream_cano_from_cano_fwd(halo_buf, wrk_fwd1) call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 2, np_max, parmom_cano, nrecv, halo_buf) call map%upstream_cano_from_cano_fwd2(halo_buf, wrk_fwd2) ! parmom_sonic = n*cs parmom_sonic = parmom_cano call parbnd_taylor_cano%set_guard_points(equi, mesh_cano, 1.0_GP, 0.0_GP, ux, & parmom_sonic, wrk_bwd2, wrk_bwd1, wrk_fwd1, wrk_fwd2, parbnd_taylor_order, psec_sonic) ! Extrapolation of parmom parmom_extrapolate = parmom_cano call parbnd_taylor_cano%extrapolate_guard_points(equi, mesh_cano, & parmom_extrapolate, wrk_bwd2, wrk_bwd1, wrk_fwd1, wrk_fwd2, psec_extrapolate) select case(parbnd_type_parmom) case('Sonic') ! parmom = n*c_s !$omp parallel default(none) private(l) shared(np_max, masks_cano, parmom_cano, parmom_sonic) !$omp do do l = 1, np_max if (.not.masks_cano%is_in_vessel(l)) then parmom_cano(l) = parmom_sonic(l) endif enddo !$omp end do !$omp end parallel case('Bohm-Chodura') ! parmom > n*c_s !$omp parallel default(none) private(l) & !$omp shared(np_max, masks_cano, parmom_cano, parmom_extrapolate, parmom_sonic, & !$omp psec_extrapolate, psec_sonic) !$omp do do l = 1, np_max if (.not.masks_cano%is_in_vessel(l)) then if (abs(psec_extrapolate(l)) > abs(psec_sonic(l))) then parmom_cano(l) = parmom_extrapolate(l) else parmom_cano(l) = parmom_sonic(l) endif endif enddo !$omp end do !$omp end parallel case default call handle_error('parbnd_type_parmom not valid or available with Taylor method: ' & // parbnd_type_parmom, GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select enddo ! Map back to staggered mesh call getdata_fwdbwdplane(comm_handler%get_comm_planes(), 1, np_max, parmom_cano, nrecv, halo_buf) call map%lift_cano_to_stag(parmom_cano, halo_buf, parmom) if (allocated(halo_buf)) deallocate(halo_buf) end subroutine end module