pnpn_res_cpu.f90

Go to the documentation of this file.

 
module pnpn_res_cpu
  use gather_scatter, only : gs_t, gs_op_add
  use operators, only : opgrad, curl, cdtp, rotate_cyc
  use field, only : field_t
  use ax_product, only : ax_t
  use coefs, only : coef_t
  use facet_normal, only : facet_normal_t
  use pnpn_residual, only : pnpn_prs_res_t, pnpn_vel_res_t
  use scratch_registry, only : neko_scratch_registry
  use mesh, only : mesh_t
  use num_types, only : rp
  use space, only : space_t
  use math, only : copy, cmult2, invers2, rzero
  use, intrinsic :: iso_c_binding, only : c_ptr
  implicit none
  private
 
  type, public, extends(pnpn_prs_res_t) :: pnpn_prs_res_cpu_t
   contains
     procedure, nopass :: compute => pnpn_prs_res_cpu_compute
  end type pnpn_prs_res_cpu_t
 
  type, public, extends(pnpn_vel_res_t) :: pnpn_vel_res_cpu_t
   contains
     procedure, nopass :: compute => pnpn_vel_res_cpu_compute
  end type pnpn_vel_res_cpu_t
 
contains
 
  subroutine pnpn_prs_res_cpu_compute(p, p_res, u, v, w, u_e, v_e, w_e, f_x, &
       f_y, f_z, c_Xh, gs_Xh, bc_prs_surface, bc_sym_surface, Ax, bd, dt, mu, &
       rho, event)
    type(field_t), intent(inout) :: p, u, v, w
    type(field_t), intent(in) :: u_e, v_e, w_e
    type(field_t), intent(inout) :: p_res
    type(field_t), intent(in) :: f_x, f_y, f_z
    type(coef_t), intent(inout) :: c_Xh
    type(gs_t), intent(inout) :: gs_Xh
    type(facet_normal_t), intent(in) :: bc_prs_surface
    type(facet_normal_t), intent(in) :: bc_sym_surface
    class(ax_t), intent(inout) :: Ax
    real(kind=rp), intent(in) :: bd
    real(kind=rp), intent(in) :: dt
    type(field_t), intent(in) :: mu
    type(field_t), intent(in) :: rho
    type(c_ptr), intent(inout) :: event
    real(kind=rp) :: dtbd, rho_inv, mu_rho
    integer :: n
    integer :: i
    type(field_t), pointer :: ta1, ta2, ta3, wa1, wa2, wa3, work1, work2
    integer :: temp_indices(8)
 
    call neko_scratch_registry%request_field(ta1, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(ta2, temp_indices(2), .false.)
    call neko_scratch_registry%request_field(ta3, temp_indices(3), .false.)
    call neko_scratch_registry%request_field(wa1, temp_indices(4), .false.)
    call neko_scratch_registry%request_field(wa2, temp_indices(5), .false.)
    call neko_scratch_registry%request_field(wa3, temp_indices(6), .false.)
    call neko_scratch_registry%request_field(work1, temp_indices(7), .false.)
    call neko_scratch_registry%request_field(work2, temp_indices(8), .false.)
 
    n = c_xh%dof%size()
 
    ! We assume the material properties are constant, and only their
    ! reciprocals are needed here.
    rho_inv = 1.0_rp / rho%x(1,1,1,1)
    mu_rho = mu%x(1,1,1,1) / rho%x(1,1,1,1)
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1, n
       c_xh%h1(i,1,1,1) = rho_inv
       c_xh%h2(i,1,1,1) = 0.0_rp
    end do
    !$omp end parallel do
    c_xh%ifh2 = .false.
 
    call curl(ta1, ta2, ta3, u_e, v_e, w_e, work1, work2, c_xh)
    call curl(wa1, wa2, wa3, ta1, ta2, ta3, work1, work2, c_xh)
 
    ! ta = f / rho - wa * mu / rho * B
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1,n
       ta1%x(i,1,1,1) = f_x%x(i,1,1,1) * rho_inv &
            - ((wa1%x(i,1,1,1) * mu_rho) * c_xh%B(i,1,1,1))
       ta2%x(i,1,1,1) = f_y%x(i,1,1,1) * rho_inv &
            - ((wa2%x(i,1,1,1) * mu_rho) * c_xh%B(i,1,1,1))
       ta3%x(i,1,1,1) = f_z%x(i,1,1,1) * rho_inv &
            - ((wa3%x(i,1,1,1) * mu_rho) * c_xh%B(i,1,1,1))
    end do
    !$omp end parallel do
 
    call rotate_cyc(ta1%x, ta2%x, ta3%x, 1, c_xh)
    call gs_xh%op(ta1%x, ta2%x, ta3%x, n, gs_op_add)
    call rotate_cyc(ta1%x, ta2%x, ta3%x, 0, c_xh)
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1,n
       ta1%x(i,1,1,1) = ta1%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
       ta2%x(i,1,1,1) = ta2%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
       ta3%x(i,1,1,1) = ta3%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
    end do
    !$omp end parallel do
 
    call cdtp(wa1%x, ta1%x, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
    call cdtp(wa2%x, ta2%x, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
    call cdtp(wa3%x, ta3%x, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
 
    call ax%compute(p_res%x, p%x, c_xh, p%msh, p%Xh)
 
    dtbd = bd / dt
 
    !$omp parallel private(i)
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp do
    do i = 1,n
       p_res%x(i,1,1,1) = (-p_res%x(i,1,1,1)) &
            + wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1)
    end do
    !$omp end do
 
    !
    ! Surface velocity terms
    !
    call bc_prs_surface%apply_surfvec_sub(p_res%x, u%x, v%x, w%x, dtbd, n)
    call bc_sym_surface%apply_surfvec_sub(p_res%x, ta1%x, ta2%x, ta3%x, &
         1.0_rp, n)
    !$omp end parallel
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine pnpn_prs_res_cpu_compute
 
  subroutine pnpn_vel_res_cpu_compute(Ax, u, v, w, u_res, v_res, w_res, &
       p, f_x, f_y, f_z, c_Xh, msh, Xh, mu, rho, bd, dt, n)
    class(ax_t), intent(in) :: Ax
    type(mesh_t), intent(inout) :: msh
    type(space_t), intent(inout) :: Xh
    type(field_t), intent(inout) :: p, u, v, w
    type(field_t), intent(inout) :: u_res, v_res, w_res
    type(field_t), intent(in) :: f_x, f_y, f_z
    type(coef_t), intent(inout) :: c_Xh
    type(field_t), intent(in) :: mu
    type(field_t), intent(in) :: rho
    real(kind=rp), intent(in) :: bd
    real(kind=rp), intent(in) :: dt
    real(kind=rp) :: rho_val, mu_val
    integer :: temp_indices(3)
    type(field_t), pointer :: ta1, ta2, ta3
    integer, intent(in) :: n
    integer :: i
 
    ! We assume the material properties are constant
    rho_val = rho%x(1,1,1,1)
    mu_val = mu%x(1,1,1,1)
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1, n
       c_xh%h1(i,1,1,1) = mu_val
       c_xh%h2(i,1,1,1) = rho_val * bd / dt
    end do
    !$omp end parallel do
    c_xh%ifh2 = .true.
 
    ! One fused pass: streams the geometric factors and h1/h2 once instead
    ! of three times, and opens one parallel region instead of three.
    call ax%compute_vector(u_res%x, v_res%x, w_res%x, u%x, v%x, w%x, &
         c_xh, msh, xh)
    call neko_scratch_registry%request_field(ta1, temp_indices(1), .false.)
    call neko_scratch_registry%request_field(ta2, temp_indices(2), .false.)
    call neko_scratch_registry%request_field(ta3, temp_indices(3), .false.)
 
    call opgrad(ta1%x, ta2%x, ta3%x, p%x, c_xh)
 
    !OCL NORECURRENCE, NOVREC, NOALIAS
    !DIR$ CONCURRENT
    !DIR$ IVDEP
    !GCC$ ivdep
    !$omp parallel do
    do i = 1, n
       u_res%x(i,1,1,1) = (-u_res%x(i,1,1,1)) - ta1%x(i,1,1,1) + f_x%x(i,1,1,1)
       v_res%x(i,1,1,1) = (-v_res%x(i,1,1,1)) - ta2%x(i,1,1,1) + f_y%x(i,1,1,1)
       w_res%x(i,1,1,1) = (-w_res%x(i,1,1,1)) - ta3%x(i,1,1,1) + f_z%x(i,1,1,1)
    end do
    !$omp end parallel do
 
    call neko_scratch_registry%relinquish_field(temp_indices)
 
  end subroutine pnpn_vel_res_cpu_compute
 
end module pnpn_res_cpu