diff --git a/src/common/m_boundary_common.fpp b/src/common/m_boundary_common.fpp index 8689680370..a772feb931 100644 --- a/src/common/m_boundary_common.fpp +++ b/src/common/m_boundary_common.fpp @@ -14,14 +14,14 @@ module m_boundary_common use m_constants use m_boundary_primitives use m_boundary_io +#ifndef MFC_PRE_PROCESS + use m_mpi_common, only: s_mpi_sendrecv_grid_variables_buffers +#endif implicit none - private; public :: s_initialize_boundary_common_module, s_populate_variables_buffers, s_create_mpi_types, & - & s_populate_capillary_buffers, s_populate_F_igr_buffers, s_write_serial_boundary_condition_files, & - & s_write_parallel_boundary_condition_files, s_read_serial_boundary_condition_files, & - & s_read_parallel_boundary_condition_files, s_assign_default_bc_type, s_populate_grid_variables_buffers, & - & s_finalize_boundary_common_module + private; public :: s_initialize_boundary_common_module, s_populate_variables_buffers, s_populate_capillary_buffers, & + & s_populate_F_igr_buffers, s_populate_grid_variables_buffers, s_finalize_boundary_common_module public :: bc_buffers @@ -167,6 +167,230 @@ contains end subroutine s_populate_bc_direction + !> Populate ghost cell buffers for the color function and its divergence used in capillary surface tension. + impure subroutine s_populate_capillary_buffers(c_divs, bc_type, bc) + + type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs + type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type + type(bc_xyz_info), intent(in) :: bc + + call s_populate_capillary_bc_direction(1, -1, bc%x, bc_type(1, 1), c_divs) + call s_populate_capillary_bc_direction(1, 1, bc%x, bc_type(1, 2), c_divs) + + if (n == 0) return + + #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 + call s_populate_capillary_bc_direction(2, -1, bc%y, bc_type(2, 1), c_divs) + call s_populate_capillary_bc_direction(2, 1, bc%y, bc_type(2, 2), c_divs) + #:endif + + if (p == 0) return + + #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 + call s_populate_capillary_bc_direction(3, -1, bc%z, bc_type(3, 1), c_divs) + call s_populate_capillary_bc_direction(3, 1, bc%z, bc_type(3, 2), c_divs) + #:endif + + end subroutine s_populate_capillary_buffers + + !> Populate ghost cell buffers for one capillary BC direction and location, via MPI exchange for processor boundaries or by + !! dispatching the per-cell capillary BC routines over the boundary face. + impure subroutine s_populate_capillary_bc_direction(bc_dir, bc_loc, bc_bounds, bc_type_edge, c_divs) + + integer, intent(in) :: bc_dir, bc_loc + type(int_bounds_info), intent(in) :: bc_bounds + type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs + type(integer_field), intent(in) :: bc_type_edge + integer :: bc_edge, k_beg, k_end, l_beg, l_end, k, l, bc_code + + if (bc_loc == -1) then + bc_edge = bc_bounds%beg + else + bc_edge = bc_bounds%end + end if + + if (bc_edge >= 0) then + call s_mpi_sendrecv_variables_buffers(c_divs, bc_dir, bc_loc, num_dims + 1) + return + end if + + if (bc_dir == 1) then + k_beg = 0; k_end = n; l_beg = 0; l_end = p + else if (bc_dir == 2) then + k_beg = -buff_size; k_end = m + buff_size; l_beg = 0; l_end = p + else + k_beg = -buff_size; k_end = m + buff_size; l_beg = -buff_size; l_end = n + buff_size + end if + + $:GPU_PARALLEL_LOOP(private='[l, k, bc_code]', collapse=2) + do l = l_beg, l_end + do k = k_beg, k_end + if (bc_dir == 1) then + bc_code = int(bc_type_edge%sf(0, k, l)) + else if (bc_dir == 2) then + bc_code = int(bc_type_edge%sf(k, 0, l)) + else + bc_code = int(bc_type_edge%sf(k, l, 0)) + end if + + select case (bc_code) + case (BC_PERIODIC) + call s_color_function_periodic(c_divs, bc_dir, bc_loc, k, l) + case (BC_REFLECTIVE) + call s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l) + case default + call s_color_function_ghost_cell_extrapolation(c_divs, bc_dir, bc_loc, k, l) + end select + end do + end do + $:END_GPU_PARALLEL_LOOP() + + end subroutine s_populate_capillary_bc_direction + + !> Populate ghost cell buffers for the Jacobian scalar field used in the IGR elliptic solver. + impure subroutine s_populate_F_igr_buffers(bc_type, jac_sf) + + type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type + type(scalar_field), dimension(1:), intent(inout) :: jac_sf + + call s_populate_F_igr_bc_direction(1, -1, bc_x, bc_type(1, 1), jac_sf) + call s_populate_F_igr_bc_direction(1, 1, bc_x, bc_type(1, 2), jac_sf) + + if (n == 0) return + + #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 + call s_populate_F_igr_bc_direction(2, -1, bc_y, bc_type(2, 1), jac_sf) + call s_populate_F_igr_bc_direction(2, 1, bc_y, bc_type(2, 2), jac_sf) + #:endif + + if (p == 0) return + + #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 + call s_populate_F_igr_bc_direction(3, -1, bc_z, bc_type(3, 1), jac_sf) + call s_populate_F_igr_bc_direction(3, 1, bc_z, bc_type(3, 2), jac_sf) + #:endif + + end subroutine s_populate_F_igr_buffers + + !> Populate ghost cell buffers for one IGR Jacobian BC direction and location, via MPI exchange for processor boundaries or by + !! dispatching the per-cell IGR Jacobian BC routines over the boundary face. + impure subroutine s_populate_F_igr_bc_direction(bc_dir, bc_loc, bc_bounds, bc_type_edge, jac_sf) + + integer, intent(in) :: bc_dir, bc_loc + type(int_bounds_info), intent(in) :: bc_bounds + type(integer_field), intent(in) :: bc_type_edge + type(scalar_field), dimension(1:), intent(inout) :: jac_sf + integer :: bc_edge, k_beg, k_end, l_beg, l_end, k, l, j, bc_code + + if (bc_loc == -1) then + bc_edge = bc_bounds%beg + else + bc_edge = bc_bounds%end + end if + + if (bc_edge >= 0) then + call s_mpi_sendrecv_variables_buffers(jac_sf, bc_dir, bc_loc, 1) + return + end if + + if (bc_dir == 1) then + k_beg = 0; k_end = n; l_beg = 0; l_end = p + else if (bc_dir == 2) then + k_beg = idwbuff(1)%beg; k_end = idwbuff(1)%end; l_beg = 0; l_end = p + else + k_beg = idwbuff(1)%beg; k_end = idwbuff(1)%end; l_beg = idwbuff(2)%beg; l_end = idwbuff(2)%end + end if + + $:GPU_PARALLEL_LOOP(private='[l, k, bc_code]', collapse=2) + do l = l_beg, l_end + do k = k_beg, k_end + if (bc_dir == 1) then + bc_code = int(bc_type_edge%sf(0, k, l)) + else if (bc_dir == 2) then + bc_code = int(bc_type_edge%sf(k, 0, l)) + else + bc_code = int(bc_type_edge%sf(k, l, 0)) + end if + + select case (bc_code) + case (BC_PERIODIC) + call s_F_igr_periodic(jac_sf, bc_dir, bc_loc, k, l) + case (BC_REFLECTIVE) + call s_F_igr_reflective(jac_sf, bc_dir, bc_loc, k, l) + case default + call s_F_igr_ghost_cell_extrapolation(jac_sf, bc_dir, bc_loc, k, l) + end select + end do + end do + $:END_GPU_PARALLEL_LOOP() + + end subroutine s_populate_F_igr_bc_direction + + !> Populate the buffers of the grid variables, which are constituted of the cell-boundary locations and cell-width + !! distributions, based on the boundary conditions. + subroutine s_populate_grid_variables_buffers + + integer :: i + type(int_bounds_info) :: offset_x, offset_y, offset_z + + offset_x%beg = buff_size; offset_x%end = buff_size + offset_y%beg = buff_size; offset_y%end = buff_size + offset_z%beg = buff_size; offset_z%end = buff_size + +#ifndef MFC_PRE_PROCESS + call s_populate_grid_bc_direction(1, -1, bc_x, offset_x) + call s_populate_grid_bc_direction(1, 1, bc_x, offset_x) + + if (n == 0) return + + #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 + call s_populate_grid_bc_direction(2, -1, bc_y, offset_y) + call s_populate_grid_bc_direction(2, 1, bc_y, offset_y) + #:endif + + if (p == 0) return + + #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 + call s_populate_grid_bc_direction(3, -1, bc_z, offset_z) + call s_populate_grid_bc_direction(3, 1, bc_z, offset_z) + #:endif +#endif + + end subroutine s_populate_grid_variables_buffers + +#ifndef MFC_PRE_PROCESS + !> Populate grid variable buffers (cell widths and centers) for one direction and location. + subroutine s_populate_grid_bc_direction(bc_dir, bc_loc, bc_bounds, offset_dir) + + integer, intent(in) :: bc_dir, bc_loc + type(int_bounds_info), intent(in) :: bc_bounds, offset_dir + integer :: bc_edge + + if (bc_loc == -1) then + bc_edge = bc_bounds%beg + else + bc_edge = bc_bounds%end + end if + + if (bc_edge >= 0) then + call s_mpi_sendrecv_grid_variables_buffers(bc_dir, bc_loc) + return + end if + + select case (bc_edge) + case (BC_PERIODIC) + call s_grid_periodic_bc(bc_dir, bc_loc, offset_dir) + case (BC_REFLECTIVE) + call s_grid_reflective_bc(bc_dir, bc_loc, offset_dir) + case (BC_AXIS) + if (bc_dir == 2) call s_grid_axis_bc(bc_loc, offset_dir) + case default + call s_grid_ghost_cell_extrapolation_bc(bc_dir, bc_loc, offset_dir) + end select + + end subroutine s_populate_grid_bc_direction +#endif + !> Deallocate boundary condition buffer arrays allocated during module initialization. subroutine s_finalize_boundary_common_module() diff --git a/src/common/m_boundary_io.fpp b/src/common/m_boundary_io.fpp index c7d0f9ffef..e926c39570 100644 --- a/src/common/m_boundary_io.fpp +++ b/src/common/m_boundary_io.fpp @@ -25,498 +25,6 @@ module m_boundary_io contains - !> Populate ghost cell buffers for the color function and its divergence used in capillary surface tension. - impure subroutine s_populate_capillary_buffers(c_divs, bc_type, bc) - - type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs - type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type - type(bc_xyz_info), intent(in) :: bc - integer :: k, l - - !> x-direction - - if (bc%x%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(c_divs, 1, -1, num_dims + 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = 0, n - select case (bc_type(1, 1)%sf(0, k, l)) - case (BC_PERIODIC) - call s_color_function_periodic(c_divs, 1, -1, k, l) - case (BC_REFLECTIVE) - call s_color_function_reflective(c_divs, 1, -1, k, l) - case default - call s_color_function_ghost_cell_extrapolation(c_divs, 1, -1, k, l) - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - if (bc%x%end >= 0) then - call s_mpi_sendrecv_variables_buffers(c_divs, 1, 1, num_dims + 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = 0, n - select case (bc_type(1, 2)%sf(0, k, l)) - case (BC_PERIODIC) - call s_color_function_periodic(c_divs, 1, 1, k, l) - case (BC_REFLECTIVE) - call s_color_function_reflective(c_divs, 1, 1, k, l) - case default - call s_color_function_ghost_cell_extrapolation(c_divs, 1, 1, k, l) - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - if (n == 0) return - - #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 - !> y-direction - if (bc%y%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(c_divs, 2, -1, num_dims + 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = -buff_size, m + buff_size - select case (bc_type(2, 1)%sf(k, 0, l)) - case (BC_PERIODIC) - call s_color_function_periodic(c_divs, 2, -1, k, l) - case (BC_REFLECTIVE) - call s_color_function_reflective(c_divs, 2, -1, k, l) - case default - call s_color_function_ghost_cell_extrapolation(c_divs, 2, -1, k, l) - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - if (bc%y%end >= 0) then - call s_mpi_sendrecv_variables_buffers(c_divs, 2, 1, num_dims + 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = -buff_size, m + buff_size - select case (bc_type(2, 2)%sf(k, 0, l)) - case (BC_PERIODIC) - call s_color_function_periodic(c_divs, 2, 1, k, l) - case (BC_REFLECTIVE) - call s_color_function_reflective(c_divs, 2, 1, k, l) - case default - call s_color_function_ghost_cell_extrapolation(c_divs, 2, 1, k, l) - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - #:endif - - if (p == 0) return - - #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 - !> z-direction - if (bc%z%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(c_divs, 3, -1, num_dims + 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = -buff_size, n + buff_size - do k = -buff_size, m + buff_size - select case (bc_type(3, 1)%sf(k, l, 0)) - case (BC_PERIODIC) - call s_color_function_periodic(c_divs, 3, -1, k, l) - case (BC_REFLECTIVE) - call s_color_function_reflective(c_divs, 3, -1, k, l) - case default - call s_color_function_ghost_cell_extrapolation(c_divs, 3, -1, k, l) - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - if (bc%z%end >= 0) then - call s_mpi_sendrecv_variables_buffers(c_divs, 3, 1, num_dims + 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = -buff_size, n + buff_size - do k = -buff_size, m + buff_size - select case (bc_type(3, 2)%sf(k, l, 0)) - case (BC_PERIODIC) - call s_color_function_periodic(c_divs, 3, 1, k, l) - case (BC_REFLECTIVE) - call s_color_function_reflective(c_divs, 3, 1, k, l) - case default - call s_color_function_ghost_cell_extrapolation(c_divs, 3, 1, k, l) - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - #:endif - - end subroutine s_populate_capillary_buffers - - !> Apply periodic boundary conditions to the color function and its divergence fields. - subroutine s_color_function_periodic(c_divs, bc_dir, bc_loc, k, l) - - $:GPU_ROUTINE(function_name='s_color_function_periodic', parallelism='[seq]', cray_inline=True) - type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs - integer, intent(in) :: bc_dir, bc_loc - integer, intent(in) :: k, l - integer :: j, i - - if (bc_dir == 1) then !< x-direction - if (bc_loc == -1) then ! bc_x%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(-j, k, l) = c_divs(i)%sf(m - (j - 1), k, l) - end do - end do - else !< bc_x%end - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(m + j, k, l) = c_divs(i)%sf(j - 1, k, l) - end do - end do - end if - else if (bc_dir == 2) then !< y-direction - if (bc_loc == -1) then !< bc_y%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, -j, l) = c_divs(i)%sf(k, n - (j - 1), l) - end do - end do - else !< bc_y%end - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, n + j, l) = c_divs(i)%sf(k, j - 1, l) - end do - end do - end if - else if (bc_dir == 3) then !< z-direction - if (bc_loc == -1) then !< bc_z%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, l, -j) = c_divs(i)%sf(k, l, p - (j - 1)) - end do - end do - else !< bc_z%end - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, l, p + j) = c_divs(i)%sf(k, l, j - 1) - end do - end do - end if - end if - - end subroutine s_color_function_periodic - - !> Apply reflective boundary conditions to the color function and its divergence fields. - subroutine s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l) - - $:GPU_ROUTINE(function_name='s_color_function_reflective', parallelism='[seq]', cray_inline=True) - type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs - integer, intent(in) :: bc_dir, bc_loc - integer, intent(in) :: k, l - integer :: j, i - - if (bc_dir == 1) then !< x-direction - if (bc_loc == -1) then ! bc_x%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - if (i == bc_dir) then - c_divs(i)%sf(-j, k, l) = -c_divs(i)%sf(j - 1, k, l) - else - c_divs(i)%sf(-j, k, l) = c_divs(i)%sf(j - 1, k, l) - end if - end do - end do - else !< bc_x%end - do i = 1, num_dims + 1 - do j = 1, buff_size - if (i == bc_dir) then - c_divs(i)%sf(m + j, k, l) = -c_divs(i)%sf(m - (j - 1), k, l) - else - c_divs(i)%sf(m + j, k, l) = c_divs(i)%sf(m - (j - 1), k, l) - end if - end do - end do - end if - else if (bc_dir == 2) then !< y-direction - if (bc_loc == -1) then !< bc_y%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - if (i == bc_dir) then - c_divs(i)%sf(k, -j, l) = -c_divs(i)%sf(k, j - 1, l) - else - c_divs(i)%sf(k, -j, l) = c_divs(i)%sf(k, j - 1, l) - end if - end do - end do - else !< bc_y%end - do i = 1, num_dims + 1 - do j = 1, buff_size - if (i == bc_dir) then - c_divs(i)%sf(k, n + j, l) = -c_divs(i)%sf(k, n - (j - 1), l) - else - c_divs(i)%sf(k, n + j, l) = c_divs(i)%sf(k, n - (j - 1), l) - end if - end do - end do - end if - else if (bc_dir == 3) then !< z-direction - if (bc_loc == -1) then !< bc_z%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - if (i == bc_dir) then - c_divs(i)%sf(k, l, -j) = -c_divs(i)%sf(k, l, j - 1) - else - c_divs(i)%sf(k, l, -j) = c_divs(i)%sf(k, l, j - 1) - end if - end do - end do - else !< bc_z%end - do i = 1, num_dims + 1 - do j = 1, buff_size - if (i == bc_dir) then - c_divs(i)%sf(k, l, p + j) = -c_divs(i)%sf(k, l, p - (j - 1)) - else - c_divs(i)%sf(k, l, p + j) = c_divs(i)%sf(k, l, p - (j - 1)) - end if - end do - end do - end if - end if - - end subroutine s_color_function_reflective - - !> Extrapolate the color function and its divergence into ghost cells by copying boundary values. - subroutine s_color_function_ghost_cell_extrapolation(c_divs, bc_dir, bc_loc, k, l) - - $:GPU_ROUTINE(function_name='s_color_function_ghost_cell_extrapolation', parallelism='[seq]', cray_inline=True) - type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs - integer, intent(in) :: bc_dir, bc_loc - integer, intent(in) :: k, l - integer :: j, i - - if (bc_dir == 1) then !< x-direction - if (bc_loc == -1) then ! bc_x%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(-j, k, l) = c_divs(i)%sf(0, k, l) - end do - end do - else !< bc_x%end - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(m + j, k, l) = c_divs(i)%sf(m, k, l) - end do - end do - end if - else if (bc_dir == 2) then !< y-direction - if (bc_loc == -1) then !< bc_y%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, -j, l) = c_divs(i)%sf(k, 0, l) - end do - end do - else !< bc_y%end - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, n + j, l) = c_divs(i)%sf(k, n, l) - end do - end do - end if - else if (bc_dir == 3) then !< z-direction - if (bc_loc == -1) then !< bc_z%beg - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, l, -j) = c_divs(i)%sf(k, l, 0) - end do - end do - else !< bc_z%end - do i = 1, num_dims + 1 - do j = 1, buff_size - c_divs(i)%sf(k, l, p + j) = c_divs(i)%sf(k, l, p) - end do - end do - end if - end if - - end subroutine s_color_function_ghost_cell_extrapolation - - !> Populate ghost cell buffers for the Jacobian scalar field used in the IGR elliptic solver. - impure subroutine s_populate_F_igr_buffers(bc_type, jac_sf) - - type(integer_field), dimension(1:num_dims,1:2), intent(in) :: bc_type - type(scalar_field), dimension(1:), intent(inout) :: jac_sf - integer :: j, k, l - - if (bc_x%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(jac_sf, 1, -1, 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = 0, n - select case (bc_type(1, 1)%sf(0, k, l)) - case (BC_PERIODIC) - do j = 1, buff_size - jac_sf(1)%sf(-j, k, l) = jac_sf(1)%sf(m - j + 1, k, l) - end do - case (BC_REFLECTIVE) - do j = 1, buff_size - jac_sf(1)%sf(-j, k, l) = jac_sf(1)%sf(j - 1, k, l) - end do - case default - do j = 1, buff_size - jac_sf(1)%sf(-j, k, l) = jac_sf(1)%sf(0, k, l) - end do - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - if (bc_x%end >= 0) then - call s_mpi_sendrecv_variables_buffers(jac_sf, 1, 1, 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = 0, n - select case (bc_type(1, 2)%sf(0, k, l)) - case (BC_PERIODIC) - do j = 1, buff_size - jac_sf(1)%sf(m + j, k, l) = jac_sf(1)%sf(j - 1, k, l) - end do - case (BC_REFLECTIVE) - do j = 1, buff_size - jac_sf(1)%sf(m + j, k, l) = jac_sf(1)%sf(m - (j - 1), k, l) - end do - case default - do j = 1, buff_size - jac_sf(1)%sf(m + j, k, l) = jac_sf(1)%sf(m, k, l) - end do - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - #:if not MFC_CASE_OPTIMIZATION or num_dims > 1 - if (n == 0) then - return - else if (bc_y%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(jac_sf, 2, -1, 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = idwbuff(1)%beg, idwbuff(1)%end - select case (bc_type(2, 1)%sf(k, 0, l)) - case (BC_PERIODIC) - do j = 1, buff_size - jac_sf(1)%sf(k, -j, l) = jac_sf(1)%sf(k, n - j + 1, l) - end do - case (BC_REFLECTIVE) - do j = 1, buff_size - jac_sf(1)%sf(k, -j, l) = jac_sf(1)%sf(k, j - 1, l) - end do - case default - do j = 1, buff_size - jac_sf(1)%sf(k, -j, l) = jac_sf(1)%sf(k, 0, l) - end do - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - if (bc_y%end >= 0) then - call s_mpi_sendrecv_variables_buffers(jac_sf, 2, 1, 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = 0, p - do k = idwbuff(1)%beg, idwbuff(1)%end - select case (bc_type(2, 2)%sf(k, 0, l)) - case (BC_PERIODIC) - do j = 1, buff_size - jac_sf(1)%sf(k, n + j, l) = jac_sf(1)%sf(k, j - 1, l) - end do - case (BC_REFLECTIVE) - do j = 1, buff_size - jac_sf(1)%sf(k, n + j, l) = jac_sf(1)%sf(k, n - (j - 1), l) - end do - case default - do j = 1, buff_size - jac_sf(1)%sf(k, n + j, l) = jac_sf(1)%sf(k, n, l) - end do - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - #:endif - - #:if not MFC_CASE_OPTIMIZATION or num_dims > 2 - if (p == 0) then - return - else if (bc_z%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(jac_sf, 3, -1, 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = idwbuff(2)%beg, idwbuff(2)%end - do k = idwbuff(1)%beg, idwbuff(1)%end - select case (bc_type(3, 1)%sf(k, l, 0)) - case (BC_PERIODIC) - do j = 1, buff_size - jac_sf(1)%sf(k, l, -j) = jac_sf(1)%sf(k, l, p - j + 1) - end do - case (BC_REFLECTIVE) - do j = 1, buff_size - jac_sf(1)%sf(k, l, -j) = jac_sf(1)%sf(k, l, j - 1) - end do - case default - do j = 1, buff_size - jac_sf(1)%sf(k, l, -j) = jac_sf(1)%sf(k, l, 0) - end do - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - - if (bc_z%end >= 0) then - call s_mpi_sendrecv_variables_buffers(jac_sf, 3, 1, 1) - else - $:GPU_PARALLEL_LOOP(private='[l, k]', collapse=2) - do l = idwbuff(2)%beg, idwbuff(2)%end - do k = idwbuff(1)%beg, idwbuff(1)%end - select case (bc_type(3, 2)%sf(k, l, 0)) - case (BC_PERIODIC) - do j = 1, buff_size - jac_sf(1)%sf(k, l, p + j) = jac_sf(1)%sf(k, l, j - 1) - end do - case (BC_REFLECTIVE) - do j = 1, buff_size - jac_sf(1)%sf(k, l, p + j) = jac_sf(1)%sf(k, l, p - (j - 1)) - end do - case default - do j = 1, buff_size - jac_sf(1)%sf(k, l, p + j) = jac_sf(1)%sf(k, l, p) - end do - end select - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - #:endif - - end subroutine s_populate_F_igr_buffers - !> Create MPI derived datatypes for boundary condition type arrays and buffer arrays used in parallel I/O. impure subroutine s_create_mpi_types(bc_type) @@ -847,189 +355,4 @@ contains end subroutine s_assign_default_bc_type - !> Populate the buffers of the grid variables, which are constituted of the cell-boundary locations and cell-width - !! distributions, based on the boundary conditions. - subroutine s_populate_grid_variables_buffers - - integer :: i - -#ifdef MFC_SIMULATION - ! Required for compatibility between codes - type(int_bounds_info) :: offset_x, offset_y, offset_z - - offset_x%beg = buff_size; offset_x%end = buff_size - offset_y%beg = buff_size; offset_y%end = buff_size - offset_z%beg = buff_size; offset_z%end = buff_size -#endif - -#ifndef MFC_PRE_PROCESS - ! Population of Buffers in x-direction - - ! Populating cell-width distribution buffer at bc_x%beg - if (bc_x%beg >= 0) then - call s_mpi_sendrecv_grid_variables_buffers(1, -1) - else if (bc_x%beg <= BC_GHOST_EXTRAP) then - do i = 1, buff_size - dx(-i) = dx(0) - end do - else if (bc_x%beg == BC_REFLECTIVE) then - do i = 1, buff_size - dx(-i) = dx(i - 1) - end do - else if (bc_x%beg == BC_PERIODIC) then - do i = 1, buff_size - dx(-i) = dx(m - (i - 1)) - end do - end if - - ! Computing the cell-boundary and center locations buffer at bc_x%beg - do i = 1, offset_x%beg - x_cb(-1 - i) = x_cb(-i) - dx(-i) - end do - - do i = 1, buff_size - x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp - end do - - ! Populating the cell-width distribution buffer at bc_x%end - if (bc_x%end >= 0) then - call s_mpi_sendrecv_grid_variables_buffers(1, 1) - else if (bc_x%end <= BC_GHOST_EXTRAP) then - do i = 1, buff_size - dx(m + i) = dx(m) - end do - else if (bc_x%end == BC_REFLECTIVE) then - do i = 1, buff_size - dx(m + i) = dx(m - (i - 1)) - end do - else if (bc_x%end == BC_PERIODIC) then - do i = 1, buff_size - dx(m + i) = dx(i - 1) - end do - end if - - ! Populating the cell-boundary and center locations buffer at bc_x%end - do i = 1, offset_x%end - x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) - end do - - do i = 1, buff_size - x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp - end do - - ! Population of Buffers in y-direction - - ! Populating cell-width distribution buffer at bc_y%beg - if (n == 0) then - return - else if (bc_y%beg >= 0) then - call s_mpi_sendrecv_grid_variables_buffers(2, -1) - else if (bc_y%beg <= BC_GHOST_EXTRAP .and. bc_y%beg /= BC_AXIS) then - do i = 1, buff_size - dy(-i) = dy(0) - end do - else if (bc_y%beg == BC_REFLECTIVE .or. bc_y%beg == BC_AXIS) then - do i = 1, buff_size - dy(-i) = dy(i - 1) - end do - else if (bc_y%beg == BC_PERIODIC) then - do i = 1, buff_size - dy(-i) = dy(n - (i - 1)) - end do - end if - - ! Computing the cell-boundary and center locations buffer at bc_y%beg - do i = 1, offset_y%beg - y_cb(-1 - i) = y_cb(-i) - dy(-i) - end do - - do i = 1, buff_size - y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp - end do - - ! Populating the cell-width distribution buffer at bc_y%end - if (bc_y%end >= 0) then - call s_mpi_sendrecv_grid_variables_buffers(2, 1) - else if (bc_y%end <= BC_GHOST_EXTRAP) then - do i = 1, buff_size - dy(n + i) = dy(n) - end do - else if (bc_y%end == BC_REFLECTIVE) then - do i = 1, buff_size - dy(n + i) = dy(n - (i - 1)) - end do - else if (bc_y%end == BC_PERIODIC) then - do i = 1, buff_size - dy(n + i) = dy(i - 1) - end do - end if - - ! Populating the cell-boundary and center locations buffer at bc_y%end - do i = 1, offset_y%end - y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) - end do - - do i = 1, buff_size - y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp - end do - - ! Population of Buffers in z-direction - - ! Populating cell-width distribution buffer at bc_z%beg - if (p == 0) then - return - else if (Bc_z%beg >= 0) then - call s_mpi_sendrecv_grid_variables_buffers(3, -1) - else if (bc_z%beg <= BC_GHOST_EXTRAP) then - do i = 1, buff_size - dz(-i) = dz(0) - end do - else if (bc_z%beg == BC_REFLECTIVE) then - do i = 1, buff_size - dz(-i) = dz(i - 1) - end do - else if (bc_z%beg == BC_PERIODIC) then - do i = 1, buff_size - dz(-i) = dz(p - (i - 1)) - end do - end if - - ! Computing the cell-boundary and center locations buffer at bc_z%beg - do i = 1, offset_z%beg - z_cb(-1 - i) = z_cb(-i) - dz(-i) - end do - - do i = 1, buff_size - z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp - end do - - ! Populating the cell-width distribution buffer at bc_z%end - if (bc_z%end >= 0) then - call s_mpi_sendrecv_grid_variables_buffers(3, 1) - else if (bc_z%end <= BC_GHOST_EXTRAP) then - do i = 1, buff_size - dz(p + i) = dz(p) - end do - else if (bc_z%end == BC_REFLECTIVE) then - do i = 1, buff_size - dz(p + i) = dz(p - (i - 1)) - end do - else if (bc_z%end == BC_PERIODIC) then - do i = 1, buff_size - dz(p + i) = dz(i - 1) - end do - end if - - ! Populating the cell-boundary and center locations buffer at bc_z%end - do i = 1, offset_z%end - z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) - end do - - do i = 1, buff_size - z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp - end do -#endif - - end subroutine s_populate_grid_variables_buffers - end module m_boundary_io diff --git a/src/common/m_boundary_primitives.fpp b/src/common/m_boundary_primitives.fpp index 41916ac115..ecc66ee8f5 100644 --- a/src/common/m_boundary_primitives.fpp +++ b/src/common/m_boundary_primitives.fpp @@ -1060,4 +1060,575 @@ contains end subroutine s_qbmm_extrapolation + !> Apply periodic boundary conditions to the color function and its divergence fields. + subroutine s_color_function_periodic(c_divs, bc_dir, bc_loc, k, l) + + $:GPU_ROUTINE(function_name='s_color_function_periodic', parallelism='[seq]', cray_inline=True) + type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs + integer, intent(in) :: bc_dir, bc_loc + integer, intent(in) :: k, l + integer :: j, i + + if (bc_dir == 1) then !< x-direction + if (bc_loc == -1) then ! bc_x%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(-j, k, l) = c_divs(i)%sf(m - (j - 1), k, l) + end do + end do + else !< bc_x%end + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(m + j, k, l) = c_divs(i)%sf(j - 1, k, l) + end do + end do + end if + else if (bc_dir == 2) then !< y-direction + if (bc_loc == -1) then !< bc_y%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, -j, l) = c_divs(i)%sf(k, n - (j - 1), l) + end do + end do + else !< bc_y%end + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, n + j, l) = c_divs(i)%sf(k, j - 1, l) + end do + end do + end if + else if (bc_dir == 3) then !< z-direction + if (bc_loc == -1) then !< bc_z%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, l, -j) = c_divs(i)%sf(k, l, p - (j - 1)) + end do + end do + else !< bc_z%end + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, l, p + j) = c_divs(i)%sf(k, l, j - 1) + end do + end do + end if + end if + + end subroutine s_color_function_periodic + + !> Apply reflective boundary conditions to the color function and its divergence fields. + subroutine s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l) + + $:GPU_ROUTINE(function_name='s_color_function_reflective', parallelism='[seq]', cray_inline=True) + type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs + integer, intent(in) :: bc_dir, bc_loc + integer, intent(in) :: k, l + integer :: j, i + + if (bc_dir == 1) then !< x-direction + if (bc_loc == -1) then ! bc_x%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + if (i == bc_dir) then + c_divs(i)%sf(-j, k, l) = -c_divs(i)%sf(j - 1, k, l) + else + c_divs(i)%sf(-j, k, l) = c_divs(i)%sf(j - 1, k, l) + end if + end do + end do + else !< bc_x%end + do i = 1, num_dims + 1 + do j = 1, buff_size + if (i == bc_dir) then + c_divs(i)%sf(m + j, k, l) = -c_divs(i)%sf(m - (j - 1), k, l) + else + c_divs(i)%sf(m + j, k, l) = c_divs(i)%sf(m - (j - 1), k, l) + end if + end do + end do + end if + else if (bc_dir == 2) then !< y-direction + if (bc_loc == -1) then !< bc_y%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + if (i == bc_dir) then + c_divs(i)%sf(k, -j, l) = -c_divs(i)%sf(k, j - 1, l) + else + c_divs(i)%sf(k, -j, l) = c_divs(i)%sf(k, j - 1, l) + end if + end do + end do + else !< bc_y%end + do i = 1, num_dims + 1 + do j = 1, buff_size + if (i == bc_dir) then + c_divs(i)%sf(k, n + j, l) = -c_divs(i)%sf(k, n - (j - 1), l) + else + c_divs(i)%sf(k, n + j, l) = c_divs(i)%sf(k, n - (j - 1), l) + end if + end do + end do + end if + else if (bc_dir == 3) then !< z-direction + if (bc_loc == -1) then !< bc_z%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + if (i == bc_dir) then + c_divs(i)%sf(k, l, -j) = -c_divs(i)%sf(k, l, j - 1) + else + c_divs(i)%sf(k, l, -j) = c_divs(i)%sf(k, l, j - 1) + end if + end do + end do + else !< bc_z%end + do i = 1, num_dims + 1 + do j = 1, buff_size + if (i == bc_dir) then + c_divs(i)%sf(k, l, p + j) = -c_divs(i)%sf(k, l, p - (j - 1)) + else + c_divs(i)%sf(k, l, p + j) = c_divs(i)%sf(k, l, p - (j - 1)) + end if + end do + end do + end if + end if + + end subroutine s_color_function_reflective + + !> Extrapolate the color function and its divergence into ghost cells by copying boundary values. + subroutine s_color_function_ghost_cell_extrapolation(c_divs, bc_dir, bc_loc, k, l) + + $:GPU_ROUTINE(function_name='s_color_function_ghost_cell_extrapolation', parallelism='[seq]', cray_inline=True) + type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs + integer, intent(in) :: bc_dir, bc_loc + integer, intent(in) :: k, l + integer :: j, i + + if (bc_dir == 1) then !< x-direction + if (bc_loc == -1) then ! bc_x%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(-j, k, l) = c_divs(i)%sf(0, k, l) + end do + end do + else !< bc_x%end + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(m + j, k, l) = c_divs(i)%sf(m, k, l) + end do + end do + end if + else if (bc_dir == 2) then !< y-direction + if (bc_loc == -1) then !< bc_y%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, -j, l) = c_divs(i)%sf(k, 0, l) + end do + end do + else !< bc_y%end + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, n + j, l) = c_divs(i)%sf(k, n, l) + end do + end do + end if + else if (bc_dir == 3) then !< z-direction + if (bc_loc == -1) then !< bc_z%beg + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, l, -j) = c_divs(i)%sf(k, l, 0) + end do + end do + else !< bc_z%end + do i = 1, num_dims + 1 + do j = 1, buff_size + c_divs(i)%sf(k, l, p + j) = c_divs(i)%sf(k, l, p) + end do + end do + end if + end if + + end subroutine s_color_function_ghost_cell_extrapolation + + !> Apply periodic boundary conditions to the IGR Jacobian field by copying values from the opposite domain boundary. + subroutine s_F_igr_periodic(jac_sf, bc_dir, bc_loc, k, l) + + $:GPU_ROUTINE(function_name='s_F_igr_periodic', parallelism='[seq]', cray_inline=True) + type(scalar_field), dimension(1:), intent(inout) :: jac_sf + integer, intent(in) :: bc_dir, bc_loc + integer, intent(in) :: k, l + integer :: j + + if (bc_dir == 1) then + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(-j, k, l) = jac_sf(1)%sf(m - j + 1, k, l) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(m + j, k, l) = jac_sf(1)%sf(j - 1, k, l) + end do + end if + else if (bc_dir == 2) then + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(k, -j, l) = jac_sf(1)%sf(k, n - j + 1, l) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(k, n + j, l) = jac_sf(1)%sf(k, j - 1, l) + end do + end if + else + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(k, l, -j) = jac_sf(1)%sf(k, l, p - j + 1) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(k, l, p + j) = jac_sf(1)%sf(k, l, j - 1) + end do + end if + end if + + end subroutine s_F_igr_periodic + + !> Apply reflective boundary conditions to the IGR Jacobian field by mirroring values across the boundary. + subroutine s_F_igr_reflective(jac_sf, bc_dir, bc_loc, k, l) + + $:GPU_ROUTINE(function_name='s_F_igr_reflective', parallelism='[seq]', cray_inline=True) + type(scalar_field), dimension(1:), intent(inout) :: jac_sf + integer, intent(in) :: bc_dir, bc_loc + integer, intent(in) :: k, l + integer :: j + + if (bc_dir == 1) then + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(-j, k, l) = jac_sf(1)%sf(j - 1, k, l) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(m + j, k, l) = jac_sf(1)%sf(m - (j - 1), k, l) + end do + end if + else if (bc_dir == 2) then + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(k, -j, l) = jac_sf(1)%sf(k, j - 1, l) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(k, n + j, l) = jac_sf(1)%sf(k, n - (j - 1), l) + end do + end if + else + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(k, l, -j) = jac_sf(1)%sf(k, l, j - 1) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(k, l, p + j) = jac_sf(1)%sf(k, l, p - (j - 1)) + end do + end if + end if + + end subroutine s_F_igr_reflective + + !> Extrapolate the IGR Jacobian field into ghost cells by copying boundary values. + subroutine s_F_igr_ghost_cell_extrapolation(jac_sf, bc_dir, bc_loc, k, l) + + $:GPU_ROUTINE(function_name='s_F_igr_ghost_cell_extrapolation', parallelism='[seq]', cray_inline=True) + type(scalar_field), dimension(1:), intent(inout) :: jac_sf + integer, intent(in) :: bc_dir, bc_loc + integer, intent(in) :: k, l + integer :: j + + if (bc_dir == 1) then + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(-j, k, l) = jac_sf(1)%sf(0, k, l) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(m + j, k, l) = jac_sf(1)%sf(m, k, l) + end do + end if + else if (bc_dir == 2) then + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(k, -j, l) = jac_sf(1)%sf(k, 0, l) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(k, n + j, l) = jac_sf(1)%sf(k, n, l) + end do + end if + else + if (bc_loc == -1) then + do j = 1, buff_size + jac_sf(1)%sf(k, l, -j) = jac_sf(1)%sf(k, l, 0) + end do + else + do j = 1, buff_size + jac_sf(1)%sf(k, l, p + j) = jac_sf(1)%sf(k, l, p) + end do + end if + end if + + end subroutine s_F_igr_ghost_cell_extrapolation + +#ifndef MFC_PRE_PROCESS + !> Apply periodic boundary conditions to grid variables by copying cell widths from opposite domain boundary. + subroutine s_grid_periodic_bc(bc_dir, bc_loc, offset_dir) + + integer, intent(in) :: bc_dir, bc_loc + type(int_bounds_info), intent(in) :: offset_dir + integer :: i + + if (bc_dir == 1) then + if (bc_loc == -1) then + do i = 1, buff_size + dx(-i) = dx(m - (i - 1)) + end do + do i = 1, offset_dir%beg + x_cb(-1 - i) = x_cb(-i) - dx(-i) + end do + do i = 1, buff_size + x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp + end do + else + do i = 1, buff_size + dx(m + i) = dx(i - 1) + end do + do i = 1, offset_dir%end + x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) + end do + do i = 1, buff_size + x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp + end do + end if + else if (bc_dir == 2) then + if (bc_loc == -1) then + do i = 1, buff_size + dy(-i) = dy(n - (i - 1)) + end do + do i = 1, offset_dir%beg + y_cb(-1 - i) = y_cb(-i) - dy(-i) + end do + do i = 1, buff_size + y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp + end do + else + do i = 1, buff_size + dy(n + i) = dy(i - 1) + end do + do i = 1, offset_dir%end + y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) + end do + do i = 1, buff_size + y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp + end do + end if + else + if (bc_loc == -1) then + do i = 1, buff_size + dz(-i) = dz(p - (i - 1)) + end do + do i = 1, offset_dir%beg + z_cb(-1 - i) = z_cb(-i) - dz(-i) + end do + do i = 1, buff_size + z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp + end do + else + do i = 1, buff_size + dz(p + i) = dz(i - 1) + end do + do i = 1, offset_dir%end + z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) + end do + do i = 1, buff_size + z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp + end do + end if + end if + + end subroutine s_grid_periodic_bc + + !> Apply reflective boundary conditions to grid variables by mirroring cell widths across the boundary. + subroutine s_grid_reflective_bc(bc_dir, bc_loc, offset_dir) + + integer, intent(in) :: bc_dir, bc_loc + type(int_bounds_info), intent(in) :: offset_dir + integer :: i + + if (bc_dir == 1) then + if (bc_loc == -1) then + do i = 1, buff_size + dx(-i) = dx(i - 1) + end do + do i = 1, offset_dir%beg + x_cb(-1 - i) = x_cb(-i) - dx(-i) + end do + do i = 1, buff_size + x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp + end do + else + do i = 1, buff_size + dx(m + i) = dx(m - (i - 1)) + end do + do i = 1, offset_dir%end + x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) + end do + do i = 1, buff_size + x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp + end do + end if + else if (bc_dir == 2) then + if (bc_loc == -1) then + do i = 1, buff_size + dy(-i) = dy(i - 1) + end do + do i = 1, offset_dir%beg + y_cb(-1 - i) = y_cb(-i) - dy(-i) + end do + do i = 1, buff_size + y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp + end do + else + do i = 1, buff_size + dy(n + i) = dy(n - (i - 1)) + end do + do i = 1, offset_dir%end + y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) + end do + do i = 1, buff_size + y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp + end do + end if + else + if (bc_loc == -1) then + do i = 1, buff_size + dz(-i) = dz(i - 1) + end do + do i = 1, offset_dir%beg + z_cb(-1 - i) = z_cb(-i) - dz(-i) + end do + do i = 1, buff_size + z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp + end do + else + do i = 1, buff_size + dz(p + i) = dz(p - (i - 1)) + end do + do i = 1, offset_dir%end + z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) + end do + do i = 1, buff_size + z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp + end do + end if + end if + + end subroutine s_grid_reflective_bc + + !> Extrapolate grid variables by copying boundary cell width into ghost cells. + subroutine s_grid_ghost_cell_extrapolation_bc(bc_dir, bc_loc, offset_dir) + + integer, intent(in) :: bc_dir, bc_loc + type(int_bounds_info), intent(in) :: offset_dir + integer :: i + + if (bc_dir == 1) then + if (bc_loc == -1) then + do i = 1, buff_size + dx(-i) = dx(0) + end do + do i = 1, offset_dir%beg + x_cb(-1 - i) = x_cb(-i) - dx(-i) + end do + do i = 1, buff_size + x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp + end do + else + do i = 1, buff_size + dx(m + i) = dx(m) + end do + do i = 1, offset_dir%end + x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) + end do + do i = 1, buff_size + x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp + end do + end if + else if (bc_dir == 2) then + if (bc_loc == -1) then + do i = 1, buff_size + dy(-i) = dy(0) + end do + do i = 1, offset_dir%beg + y_cb(-1 - i) = y_cb(-i) - dy(-i) + end do + do i = 1, buff_size + y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp + end do + else + do i = 1, buff_size + dy(n + i) = dy(n) + end do + do i = 1, offset_dir%end + y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) + end do + do i = 1, buff_size + y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp + end do + end if + else + if (bc_loc == -1) then + do i = 1, buff_size + dz(-i) = dz(0) + end do + do i = 1, offset_dir%beg + z_cb(-1 - i) = z_cb(-i) - dz(-i) + end do + do i = 1, buff_size + z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp + end do + else + do i = 1, buff_size + dz(p + i) = dz(p) + end do + do i = 1, offset_dir%end + z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) + end do + do i = 1, buff_size + z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp + end do + end if + end if + + end subroutine s_grid_ghost_cell_extrapolation_bc + + !> Apply axis boundary conditions to grid variables for cylindrical coordinates. + subroutine s_grid_axis_bc(bc_loc, offset_dir) + + integer, intent(in) :: bc_loc + type(int_bounds_info), intent(in) :: offset_dir + integer :: i + + if (bc_loc == -1) then + do i = 1, buff_size + dy(-i) = dy(i - 1) + end do + do i = 1, offset_dir%beg + y_cb(-1 - i) = y_cb(-i) - dy(-i) + end do + do i = 1, buff_size + y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp + end do + end if + + end subroutine s_grid_axis_bc +#endif end module m_boundary_primitives diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index a1c1498f11..dfdcff7678 100644 --- a/src/common/m_checker_common.fpp +++ b/src/common/m_checker_common.fpp @@ -15,7 +15,7 @@ module m_checker_common implicit none - private; public :: s_check_inputs_common, wp + private; public :: s_check_inputs_common contains diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index d0fdb0dd74..c0eb912744 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -1439,6 +1439,7 @@ contains #ifdef MFC_MPI integer :: ierr !< Generic flag used to identify and report MPI errors + integer :: i if (mpi_dir == 1) then if (pbc_loc == -1) then ! PBC at the beginning @@ -1449,6 +1450,10 @@ contains call MPI_SENDRECV(dx(0), buff_size, mpi_p, bc_x%beg, 1, dx(-buff_size), buff_size, mpi_p, bc_x%beg, 0, & & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) end if + do i = 1, buff_size + x_cb(-1 - i) = x_cb(-i) - dx(-i) + x_cc(-i) = x_cc(1 - i) - (dx(1 - i) + dx(-i))/2._wp + end do else ! PBC at the end if (bc_x%beg >= 0) then ! PBC at the end and beginning call MPI_SENDRECV(dx(0), buff_size, mpi_p, bc_x%beg, 1, dx(m + 1), buff_size, mpi_p, bc_x%end, 1, & @@ -1457,6 +1462,10 @@ contains call MPI_SENDRECV(dx(m - buff_size + 1), buff_size, mpi_p, bc_x%end, 0, dx(m + 1), buff_size, mpi_p, & & bc_x%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) end if + do i = 1, buff_size + x_cb(m + i) = x_cb(m + (i - 1)) + dx(m + i) + x_cc(m + i) = x_cc(m + (i - 1)) + (dx(m + (i - 1)) + dx(m + i))/2._wp + end do end if else if (mpi_dir == 2) then if (pbc_loc == -1) then ! PBC at the beginning @@ -1467,6 +1476,10 @@ contains call MPI_SENDRECV(dy(0), buff_size, mpi_p, bc_y%beg, 1, dy(-buff_size), buff_size, mpi_p, bc_y%beg, 0, & & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) end if + do i = 1, buff_size + y_cb(-1 - i) = y_cb(-i) - dy(-i) + y_cc(-i) = y_cc(1 - i) - (dy(1 - i) + dy(-i))/2._wp + end do else ! PBC at the end if (bc_y%beg >= 0) then ! PBC at the end and beginning call MPI_SENDRECV(dy(0), buff_size, mpi_p, bc_y%beg, 1, dy(n + 1), buff_size, mpi_p, bc_y%end, 1, & @@ -1475,6 +1488,10 @@ contains call MPI_SENDRECV(dy(n - buff_size + 1), buff_size, mpi_p, bc_y%end, 0, dy(n + 1), buff_size, mpi_p, & & bc_y%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) end if + do i = 1, buff_size + y_cb(n + i) = y_cb(n + (i - 1)) + dy(n + i) + y_cc(n + i) = y_cc(n + (i - 1)) + (dy(n + (i - 1)) + dy(n + i))/2._wp + end do end if else if (pbc_loc == -1) then ! PBC at the beginning @@ -1485,6 +1502,10 @@ contains call MPI_SENDRECV(dz(0), buff_size, mpi_p, bc_z%beg, 1, dz(-buff_size), buff_size, mpi_p, bc_z%beg, 0, & & MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) end if + do i = 1, buff_size + z_cb(-1 - i) = z_cb(-i) - dz(-i) + z_cc(-i) = z_cc(1 - i) - (dz(1 - i) + dz(-i))/2._wp + end do else ! PBC at the end if (bc_z%beg >= 0) then ! PBC at the end and beginning call MPI_SENDRECV(dz(0), buff_size, mpi_p, bc_z%beg, 1, dz(p + 1), buff_size, mpi_p, bc_z%end, 1, & @@ -1493,6 +1514,10 @@ contains call MPI_SENDRECV(dz(p - buff_size + 1), buff_size, mpi_p, bc_z%end, 0, dz(p + 1), buff_size, mpi_p, & & bc_z%end, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) end if + do i = 1, buff_size + z_cb(p + i) = z_cb(p + (i - 1)) + dz(p + i) + z_cc(p + i) = z_cc(p + (i - 1)) + (dz(p + (i - 1)) + dz(p + i))/2._wp + end do end if end if #endif diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index 287a8a505c..b9d6799024 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -15,6 +15,7 @@ module m_data_input use m_mpi_common use m_compile_specific use m_boundary_common + use m_boundary_io use m_helper implicit none diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 65e3fedfd6..fd5a1e15a0 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -15,6 +15,7 @@ module m_start_up use m_mpi_proxy use m_mpi_common use m_boundary_common + use m_boundary_io use m_variables_conversion use m_data_input use m_data_output diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index c0dd4c906a..0993f1ba76 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -19,6 +19,7 @@ module m_data_output use m_delay_file_access use m_boundary_common use m_boundary_conditions + use m_boundary_io use m_thermochem, only: species_names use m_helper use m_constants, only: model_eqns_5eq, precision_single diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 6bb59e7a64..506eff1939 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -19,6 +19,7 @@ module m_start_up use m_riemann_solvers use m_cbc use m_boundary_common + use m_boundary_io use m_acoustic_src use m_rhs use m_chemistry