module nonaligned_operators_m !! Implementation of parallel operators with non-aligned methods, !! These methods are mostly for testing purpose or exploring numerical techniques, !! i.e. there is no focus on performance optimisation and no unit tests were created for these routines !! - Guenter's scheme second order: (see S. Guenter et al., Journal of Computational Physics 209 (2005) 354–370 !! - Guenter's scheme fourth order: (see S. Guenter et al., Journal of Computational Physics 226 (2007) 2306-2316 use precision_grillix_m, only : GP use error_handling_grillix_m, only : handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use comm_handler_m, only : comm_handler_t use parallel_map_m, only : parallel_map_t use variable_m, only: variable_t use communication_planes_m, only : start_comm_planes, & finalize_comm_planes ! From PARALLAX use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t implicit none public :: gradient_3D_staggered public :: parallel_gradient_nonaligned public :: parallel_divergence_nonaligned private :: dvdx_2nd_order private :: dvdy_2nd_order private :: dvdz_2nd_order private :: dvdx_4th_order private :: dvdy_4th_order private :: dvdz_4th_order contains subroutine gradient_3D_staggered(comm_handler, equi, mesh, map, u, gradu, order) !! Computes gradient (dudx, dudy, dudphi = dudz) at 3D staggered positions !! Gradient is computed according to symmetric formulation of Guenter's scheme !! Staggered point is at northeast-forward position w.r.t. its considered grid point on full grid !! Second and fourth order are implemented type(comm_handler_t), intent(in) :: comm_handler !! Communicator class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(variable_t) :: u !! Variable, must be defined on canonical grid real(GP), dimension(mesh%get_n_points(), 3), intent(out) :: gradu !! Values of gradient integer, intent(in), optional :: order !! Order of scheme default: 2 integer :: gorder, l, i, j, in, jn, ln real(GP) :: hf, dphi, dudx, dudy, dudz real(GP), allocatable, dimension(:,:,:) :: v real(GP), dimension(:), pointer :: ufwd2 => null() if (present(order)) then gorder = order else gorder = 2 endif if ((gorder /= 2).and.(gorder /= 4)) then call handle_error('Order not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Order: ',[gorder])) endif hf = mesh%get_spacing_f() dphi = map%get_dphi() ! v stores values of stencil to be used allocate(v(gorder, gorder, gorder)) ufwd2 => u%get_halo(comm_handler, 2) !$OMP PARALLEL PRIVATE(l, i, j, in, jn, ln, v, dudx, dudy, dudz) !$OMP DO do l = 1, mesh%get_n_points() do i = 1, order do j = 1, order in = i - gorder/2 jn = j - gorder/2 ln = mesh%get_index_neighbor(in, jn, l) if (ln == 0) then v(i,j,:) = 0.0_GP else if (gorder == 2) then v(i,j,1) = u%vals(ln) v(i,j,2) = u%hfwd(ln) elseif (gorder == 4) then v(i,j,1) = u%hbwd(ln) v(i,j,2) = u%vals(ln) v(i,j,3) = u%hfwd(ln) v(i,j,4) = ufwd2(ln) endif endif enddo enddo if (gorder == 2) then dudx = dvdx_2nd_order(hf, v) dudy = dvdy_2nd_order(hf, v) dudz = dvdz_2nd_order(dphi, v) elseif (gorder == 4) then dudx = dvdx_4th_order(hf, v) dudy = dvdy_4th_order(hf, v) dudz = dvdz_4th_order(dphi, v) endif gradu(l,1) = dudx gradu(l,2) = dudy gradu(l,3) = dudz enddo !$OMP END DO !$OMP END PARALLEL ufwd2 => null() deallocate(v) end subroutine subroutine parallel_gradient_nonaligned(comm_handler, equi, mesh, map, u, pargradu, order) !! Computes parallel gradient at 3D staggered positions !! Based on gradient_3D_staggered type(comm_handler_t), intent(in) :: comm_handler !! Communicator class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(variable_t), intent(in) :: u !! Variable, must be defined on canonical grid real(GP), dimension(mesh%get_n_points()), intent(out) :: pargradu !! Values of parallel gradient integer, intent(in), optional :: order !! Order of scheme default: 2 integer :: gorder, l real(GP) :: phi, hf, xs, ys, bxs, bys, bzs, absbs real(GP), dimension(mesh%get_n_points(), 3) :: gradu if (present(order)) then gorder = order else gorder = 2 endif if ((gorder /= 2).and.(gorder /= 4)) then call handle_error('Order not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Order: ',[gorder])) endif ! Compute 3D gradient call gradient_3D_staggered(comm_handler, equi, mesh, map, u, gradu, gorder) hf = mesh%get_spacing_f() phi = mesh%get_phi() !$OMP PARALLEL PRIVATE(l, xs, ys, bxs, bys, bzs, absbs) !$OMP DO do l = 1, mesh%get_n_points() ! Coordinates of staggered positions xs = mesh%get_x(l) + hf / 2.0_GP ys = mesh%get_y(l) + hf / 2.0_GP bxs = equi%bx(xs, ys, phi) bys = equi%by(xs, ys, phi) bzs = equi%btor(xs,ys, phi) absbs = equi%absb(xs,ys, phi) bxs = bxs / absbs bys = bys / absbs bzs = bzs / absbs pargradu(l) = bxs*gradu(l,1) + bys*gradu(l,2) + bzs*gradu(l,3) enddo !$OMP END DO !$OMP END PARALLEL end subroutine subroutine parallel_divergence_nonaligned(comm, equi, mesh, map, uflx, pardivu, order) !! Computes parallel divergence from parallel flux defined at 3D staggered positions !! Result (pardivu) is defined on full grid !! Based on symmetric formulation of Guenter's scheme (2nd and 4th order available) integer, intent(in) :: comm !! MPI communicator class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(inout) :: mesh !! Mesh within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map real(GP), dimension(mesh%get_n_points()), intent(in) :: uflx !! Parallel flux real(GP), dimension(mesh%get_n_points_inner()), intent(out) :: pardivu !! Parallel divergence of uflx integer, intent(in), optional :: order !! Order of scheme default: 2 integer :: gorder, ndim, ki, l, i, j, in, jn, ln real(GP) :: phi, hf, dphi, xs, ys, absbs, bxs, bys, bzs, dudx, dudy, dudz real(GP), allocatable, dimension(:,:,:) :: qparx, qpary, qparz, v integer, dimension(2) :: req_bwd, req_bwd2, req_fwd real(GP), dimension(mesh%get_n_points()) :: uflx_bwd2, uflx_bwd, uflx_fwd phi = mesh%get_phi() if (present(order)) then gorder = order else gorder = 2 endif if ((gorder /= 2).and.(gorder /= 4)) then call handle_error('Order not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Order: ',[gorder])) endif hf = mesh%get_spacing_f() dphi = map%get_dphi() ! Communicate fluxes ndim = mesh%get_n_points() call start_comm_planes(comm, ndim, ndim, uflx, -1, uflx_bwd, req_bwd) call finalize_comm_planes(req_bwd) if (gorder == 4) then call start_comm_planes(comm, ndim, ndim, uflx, -2, uflx_bwd2, req_bwd2) call finalize_comm_planes(req_bwd2) call start_comm_planes(comm, ndim, ndim, uflx, +1, uflx_fwd, req_fwd) call finalize_comm_planes(req_fwd) endif ! qparx, qpary, qparz, v stores values of stencil to be used allocate(qparx(gorder, gorder, gorder)) allocate(qpary(gorder, gorder, gorder)) allocate(qparz(gorder, gorder, gorder)) allocate(v(gorder, gorder, gorder)) !$OMP PARALLEL PRIVATE(ki, l, i, j, in, jn, ln, xs, ys, absbs, bxs, bys, bzs, qparx, qpary, qparz, & !$OMP v, dudx, dudy, dudz) !$OMP DO do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) do i = 1, gorder do j = 1, gorder in = i - gorder/2-1 jn = j - gorder/2-1 ln = mesh%get_index_neighbor(in, jn, l) if (ln == 0) then call handle_error('Neighbour not found', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('l, in, jn: ',[l,in,jn])) else xs = mesh%get_x(ln) + hf/2.0_GP ys = mesh%get_y(ln) + hf/2.0_GP absbs = equi%absb(xs, ys, phi) bxs = equi%bx(xs, ys, phi) / absbs bys = equi%by(xs, ys, phi) / absbs bzs = equi%btor(xs,ys, phi) / absbs if (gorder == 2) then qparx(i,j,1) = bxs * uflx_bwd(ln) qparx(i,j,2) = bxs * uflx(ln) qpary(i,j,1) = bys * uflx_bwd(ln) qpary(i,j,2) = bys * uflx(ln) qparz(i,j,1) = bzs * uflx_bwd(ln) qparz(i,j,2) = bzs * uflx(ln) elseif (gorder == 4) then qparx(i,j,1) = bxs * uflx_bwd2(ln) qparx(i,j,2) = bxs * uflx_bwd(ln) qparx(i,j,3) = bxs * uflx(ln) qparx(i,j,4) = bxs * uflx_fwd(ln) qpary(i,j,1) = bys * uflx_bwd2(ln) qpary(i,j,2) = bys * uflx_bwd(ln) qpary(i,j,3) = bys * uflx(ln) qpary(i,j,4) = bys * uflx_fwd(ln) qparz(i,j,1) = bzs * uflx_bwd2(ln) qparz(i,j,2) = bzs * uflx_bwd(ln) qparz(i,j,3) = bzs * uflx(ln) qparz(i,j,4) = bzs * uflx_fwd(ln) endif endif enddo enddo v = qparx if (gorder == 2) then dudx = dvdx_2nd_order(hf, v) elseif (gorder == 4) then dudx = dvdx_4th_order(hf, v) endif pardivu(ki) = dudx v = qpary if (gorder == 2) then dudy = dvdy_2nd_order(hf, v) elseif (gorder == 4) then dudy = dvdy_4th_order(hf, v) endif pardivu(ki) = pardivu(ki) + dudy v = qparz if (gorder == 2) then dudz = dvdz_2nd_order(dphi, v) elseif (gorder == 4) then dudz = dvdz_4th_order(dphi, v) endif pardivu(ki) = pardivu(ki) + dudz enddo !$OMP END DO !$OMP END PARALLEL deallocate(v) deallocate(qparz) deallocate(qpary) deallocate(qparx) end subroutine real(GP) function dvdx_2nd_order(dx, v) !! Computes dvdx according to second order symmetric finite difference real(GP), intent(in) :: dx !! Grid spacing in x-direction real(GP), dimension(2,2,2) :: v !! Values of v at stencil dvdx_2nd_order = -(v(1,1,1) + v(1,1,2) + v(1,2,1) + v(1,2,2) - v(2,1,1) - v(2,1,2) - v(2,2,1) - v(2,2,2))/(4*dx) end function real(GP) function dvdy_2nd_order(dy, v) !! Computes dvdy according to second order symmetric finite difference real(GP), intent(in) :: dy !! Grid spacing in y-direction real(GP), dimension(2,2,2) :: v !! Values of v at stencil dvdy_2nd_order = -(v(1,1,1) + v(1,1,2) - v(1,2,1) - v(1,2,2) + v(2,1,1) + v(2,1,2) - v(2,2,1) - v(2,2,2))/(4*dy) end function real(GP) function dvdz_2nd_order(dz, v) !! Computes dvdz according to second order symmetric finite difference real(GP), intent(in) :: dz !! Grid spacing in z-direction real(GP), dimension(2,2,2) :: v !! Values of v at stencil dvdz_2nd_order = -(v(1,1,1) - v(1,1,2) + v(1,2,1) - v(1,2,2) + v(2,1,1) - v(2,1,2) + v(2,2,1) - v(2,2,2))/(4*dz) end function real(GP) function dvdx_4th_order(dx, v) !! Computes dvdx according to fourth order symmetric finite difference real(GP), intent(in) :: dx !! Grid spacing in x-direction real(GP), dimension(4,4,4) :: v !! Values of v at stencil dvdx_4th_order = & (v(1,1,1) - 9*v(1,1,2) - 9*v(1,1,3) + v(1,1,4) - 9*v(1,2,1) + 81*v(1,2,2) + 81*v(1,2,3) - 9*v(1,2,4) - 9*v(1,3,1) & + 81*v(1,3,2) + 81*v(1,3,3) - 9*v(1,3,4) + v(1,4,1) - 9*v(1,4,2) - 9*v(1,4,3) + v(1,4,4) - 27*v(2,1,1) + 243*v(2,1,2) + & 243*v(2,1,3) - 27*v(2,1,4) + 243*v(2,2,1) - 2187*v(2,2,2) - 2187*v(2,2,3) + 243*v(2,2,4) + 243*v(2,3,1) - 2187*v(2,3,2) - & 2187*v(2,3,3) + 243*v(2,3,4) - 27*v(2,4,1) + 243*v(2,4,2) + 243*v(2,4,3) - 27*v(2,4,4) + 27*v(3,1,1) - 243*v(3,1,2) - & 243*v(3,1,3) + 27*v(3,1,4) - 243*v(3,2,1) + 2187*v(3,2,2) + 2187*v(3,2,3) - 243*v(3,2,4) - 243*v(3,3,1) + 2187*v(3,3,2) + & 2187*v(3,3,3) - 243*v(3,3,4) + 27*v(3,4,1) - 243*v(3,4,2) - 243*v(3,4,3) + 27*v(3,4,4) - v(4,1,1) + 9*v(4,1,2) + 9*v(4,1,3)& - v(4,1,4) + 9*v(4,2,1) - 81*v(4,2,2) - 81*v(4,2,3) + 9*v(4,2,4) + 9*v(4,3,1) - 81*v(4,3,2) - 81*v(4,3,3) + 9*v(4,3,4) - & v(4,4,1) + 9*v(4,4,2) + 9*v(4,4,3) - v(4,4,4))/(6144*dx) end function real(GP) function dvdy_4th_order(dy, v) !! Computes dvdy according to fourth order symmetric finite difference real(GP), intent(in) :: dy !! Grid spacing in x-direction real(GP), dimension(4,4,4) :: v !! Values of v at stencil dvdy_4th_order = & (v(1,1,1) - 9*v(1,1,2) - 9*v(1,1,3) + v(1,1,4) - 27*v(1,2,1) + 243*v(1,2,2) + 243*v(1,2,3) - 27*v(1,2,4) + & 27*v(1,3,1) - 243*v(1,3,2) - 243*v(1,3,3) + 27*v(1,3,4) - v(1,4,1) + 9*v(1,4,2) + 9*v(1,4,3) - v(1,4,4) - 9*v(2,1,1) + & 81*v(2,1,2) + 81*v(2,1,3) - 9*v(2,1,4) + 243*v(2,2,1) - 2187*v(2,2,2) - 2187*v(2,2,3) + 243*v(2,2,4) - 243*v(2,3,1) + & 2187*v(2,3,2) + 2187*v(2,3,3) - 243*v(2,3,4) + 9*v(2,4,1) - 81*v(2,4,2) - 81*v(2,4,3) + 9*v(2,4,4) - 9*v(3,1,1) & + 81*v(3,1,2) & + 81*v(3,1,3) - 9*v(3,1,4) + 243*v(3,2,1) - 2187*v(3,2,2) - 2187*v(3,2,3) + 243*v(3,2,4) - 243*v(3,3,1) + 2187*v(3,3,2) + & 2187*v(3,3,3) - 243*v(3,3,4) + 9*v(3,4,1) - 81*v(3,4,2) - 81*v(3,4,3) + 9*v(3,4,4) + v(4,1,1) - 9*v(4,1,2) - 9*v(4,1,3) + & v(4,1,4) - 27*v(4,2,1) + 243*v(4,2,2) + 243*v(4,2,3) - 27*v(4,2,4) + 27*v(4,3,1) - 243*v(4,3,2) - 243*v(4,3,3) & + 27*v(4,3,4) - v(4,4,1) + 9*v(4,4,2) + 9*v(4,4,3) - v(4,4,4))/(6144*dy) end function real(GP) function dvdz_4th_order(dz, v) !! Computes dvdz according to fourth order symmetric finite difference real(GP), intent(in) :: dz !! Grid spacing in x-direction real(GP), dimension(4,4,4) :: v !! Values of v at stencil dvdz_4th_order = & (v(1,1,1) - 27*v(1,1,2) + 27*v(1,1,3) - v(1,1,4) - 9*v(1,2,1) + 243*v(1,2,2) - 243*v(1,2,3) + 9*v(1,2,4) - 9*v(1,3,1) & + 243*v(1,3,2) - 243*v(1,3,3) + 9*v(1,3,4) + v(1,4,1) - 27*v(1,4,2) + 27*v(1,4,3) - v(1,4,4) - 9*v(2,1,1) + 243*v(2,1,2) - & 243*v(2,1,3) + 9*v(2,1,4) + 81*v(2,2,1) - 2187*v(2,2,2) + 2187*v(2,2,3) - 81*v(2,2,4) + 81*v(2,3,1) - 2187*v(2,3,2) + & 2187*v(2,3,3) - 81*v(2,3,4) - 9*v(2,4,1) + 243*v(2,4,2) - 243*v(2,4,3) + 9*v(2,4,4) - 9*v(3,1,1) + 243*v(3,1,2) - & 243*v(3,1,3) + 9*v(3,1,4) + 81*v(3,2,1) - 2187*v(3,2,2) + 2187*v(3,2,3) - 81*v(3,2,4) + 81*v(3,3,1) - 2187*v(3,3,2) + & 2187*v(3,3,3) - 81*v(3,3,4) - 9*v(3,4,1) + 243*v(3,4,2) - 243*v(3,4,3) + 9*v(3,4,4) + v(4,1,1) - 27*v(4,1,2) + & 27*v(4,1,3) - & v(4,1,4) - 9*v(4,2,1) + 243*v(4,2,2) - 243*v(4,2,3) + 9*v(4,2,4) - 9*v(4,3,1) + 243*v(4,3,2) - 243*v(4,3,3) + 9*v(4,3,4) + & v(4,4,1) - 27*v(4,4,2) + 27*v(4,4,3) - v(4,4,4))/(6144*dz) end function end module