scale_meshfield_filter_operation_3d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Data / Filter operation
!!
!! @par Description
!!           This module provides classes to apply filter operations to MeshField data
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_meshfield_filter_operation_3d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc, only: prc_abort
 
  use scale_element_base, only: elementbase2d, elementbase3d
  use scale_element_modalfilter, only: modalfilter
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmesh_base, only: localmeshbase
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d  
 
  use scale_meshfield_base, only: &
    meshfield2d, meshfield3d
  use scale_localmeshfield_base, only: localmeshfieldbase
  use scale_mesh_rectdom2d, only: &
    meshrectdom2d
  use scale_mesh_cubedom3d, only: &
    meshcubedom3d
  use scale_meshfieldcomm_cubedom3d, only: &
    meshfieldcommcubedom3d
  use scale_meshfieldcomm_base, only: &
    meshfieldcontainer
 
  use scale_meshfield_filter_operation_base, only: &
    meshfieldfilteroperationbase, &
    meshfieldfilteroperationbase_init, meshfieldfilteroperationbase_final, &
    filter_optrtype_convfilter,                                            &
    filter_optrtype_reconstruct, filter_optrtype_reconstruct2,             &
    apply_filter1d_x => meshfieldfilteroperationbase_apply_filter1d_x,       &
    apply_reconst1d_x => meshfieldfilteroperationbase_apply_reconst1d_x,     &
    apply_reconst1d_x_2 => meshfieldfilteroperationbase_apply_reconst1d_x_2, &
    apply_filter1d_y => meshfieldfilteroperationbase_apply_filter1d_y,       &
    apply_reconst1d_y => meshfieldfilteroperationbase_apply_reconst1d_y,     &
    apply_reconst1d_y_2 => meshfieldfilteroperationbase_apply_reconst1d_y_2
  
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !> Derived type to represent filter operation for 3D mesh field
  type, extends(meshfieldfilteroperationbase), public :: meshfieldfilteroperation3d
    type(meshfieldcommcubedom3d) :: vars_comm
    integer :: vhalosize
  contains
    procedure :: init => meshfieldfilteroperation3d_init
    procedure :: final => meshfieldfilteroperation3d_final
    procedure :: apply => meshfieldfilteroperation3d_apply_filter
  end type meshfieldfilteroperation3d
 
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  !
  private :: extract_tmp3d
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
 
  !> Initialize an object to represent filter operation for 3D mesh field
!OCL SERIAL
  subroutine meshfieldfilteroperation3d_init( this, &
    FilterOptrType,                  &
    FilterShape, FilterWidthFac,     &
    Nnode_h1D_reconst, &    
    sfield_num, hvfield_num, htensorfield_num, &
    mesh3D )
    implicit none
    class(meshfieldfilteroperation3d), intent(inout) :: this
    class(meshbase3d), intent(in) :: mesh3D
    character(*), intent(in) :: FilterOptrType
    character(*), intent(in) :: FilterShape
    real(RP), intent(in) :: FilterWidthFac
    integer, intent(in) :: Nnode_h1D_reconst
    integer, intent(in) :: sfield_num
    integer, intent(in) :: hvfield_num
    integer, intent(in) :: htensorfield_num
    !------------------------------------------
 
    call meshfieldfilteroperationbase_init( this, mesh3d%refElem3D%Nnode_h1D, &
      filteroptrtype, filtershape, filterwidthfac, &
      nnode_h1d_reconst )
 
    this%vHaloSize = 1
 
    select type(mesh3d)
    class is (meshcubedom3d)
      call this%vars_comm%Init( sfield_num, hvfield_num, htensorfield_num, mesh3d, &
        halosize_h1d=this%hHaloSize, halosize_v=this%vHaloSize )
    class default
      log_info('MeshFieldFilterOperation3D_Init',*) "Unsupported mesh type is specified. Check!"
      call prc_abort
    end select
    
    return
  end subroutine meshfieldfilteroperation3d_init
 
  !> Finalize an object to represent filter operation for 3D mesh field
!OCL SERIAL
  subroutine meshfieldfilteroperation3d_final( this )
    implicit none
    class(meshfieldfilteroperation3d), intent(inout) :: this
    !------------------------------------------    
    call meshfieldfilteroperationbase_final( this )
    call this%vars_comm%Final()
    return
  end subroutine meshfieldfilteroperation3d_final
 
  !> Apply filter operation to 3D mesh field data
!OCL serial
  subroutine meshfieldfilteroperation3d_apply_filter( this, q_list, &
    mesh3D )
    implicit none
    class(meshfieldfilteroperation3d), intent(inout) :: this
    type(meshfield3d), intent(inout), target :: q_list(this%vars_comm%field_num_tot)
    class(meshbase3d), intent(in), target :: mesh3D
 
    class(localmesh3d), pointer :: lmesh3D
    class(elementbase3d), pointer :: elem3D
    real(RP), allocatable :: tmp3D(:,:,:,:)
 
    integer :: iv
    integer :: n
    type(meshfieldcontainer) :: comm_vars_list(this%vars_comm%field_num_tot)
    !-----------------------------------------------------
 
    do iv=1, this%vars_comm%field_num_tot
      comm_vars_list(iv)%field3d => q_list(iv)
    end do
    call this%vars_comm%Put(comm_vars_list, 1)
    call this%vars_comm%Exchange()
    call this%vars_comm%Get(comm_vars_list, 1)
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lmesh3d => mesh3d%lcmesh_list(n)
      elem3d => lmesh3d%refElem3D
      allocate( tmp3d(-this%hHaloSize+1:elem3d%Nnode_h1D+this%hHaloSize,-this%hHaloSize+1:elem3d%Nnode_h1D+this%hHaloSize,elem3d%Nnode_v,lmesh3d%Ne) )
 
      do iv=1, this%vars_comm%field_num_tot
        call extract_tmp3d( tmp3d, &
          q_list(iv)%local(n)%val, q_list(iv)%local(n)%val, lmesh3d, elem3d, lmesh3d%VMapP, this%hHaloSize, 0 )
        
        select case (this%operator_type)
        case (filter_optrtype_convfilter)
          call apply_filter1d_x( q_list(iv)%local(n)%val, &
            tmp3d, this%FilterMat_h1D, elem3d%Nnode_h1D, elem3d%Nnode_h1D, elem3d%Nnode_v, lmesh3d%Ne, lmesh3d%NeA, elem3d%Nnode_h1D )
        case (filter_optrtype_reconstruct)
          call apply_reconst1d_x( q_list(iv)%local(n)%val, &
            tmp3d, this%Minv_Ml_tr, this%Minv_Mc_tr, this%Minv_Mr_tr, this%IntrpMat, elem3d%Nnode_h1D, elem3d%Nnode_h1D, elem3d%Nnode_v, lmesh3d%Ne, this%Nnode_h1D_reconst, lmesh3d )
        case (filter_optrtype_reconstruct2)
          call apply_reconst1d_x_2( q_list(iv)%local(n)%val, &
            tmp3d, this%Ml_tr, this%Mc_tr, this%Mr_tr, elem3d%Nnode_h1D, elem3d%Nnode_h1D, elem3d%Nnode_v, lmesh3d%Ne, this%Nnode_h1D_reconst, lmesh3d )
        end select
 
      end do
      deallocate( tmp3d )
    end do
 
    call this%vars_comm%Put(comm_vars_list, 1)
    call this%vars_comm%Exchange()
    call this%vars_comm%Get(comm_vars_list, 1)
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lmesh3d => mesh3d%lcmesh_list(n)
      elem3d => lmesh3d%refElem3D
      allocate( tmp3d(-this%hHaloSize+1:elem3d%Nnode_h1D+this%hHaloSize,-this%hHaloSize+1:elem3d%Nnode_h1D+this%hHaloSize,elem3d%Nnode_v,lmesh3d%Ne) )
 
      do iv=1, this%vars_comm%field_num_tot
        call extract_tmp3d( tmp3d, &
          q_list(iv)%local(n)%val, q_list(iv)%local(n)%val, lmesh3d, elem3d, lmesh3d%VMapP, this%hHaloSize, 1 )
 
        select case (this%operator_type)
        case (filter_optrtype_convfilter)
          call apply_filter1d_y( q_list(iv)%local(n)%val, &
            tmp3d, this%FilterMat_h1D, elem3d%Nnode_h1D, elem3d%Nnode_h1D, elem3d%Nnode_v,lmesh3d%Ne, lmesh3d%NeA, elem3d%Nnode_h1D )
        case (filter_optrtype_reconstruct)
          call apply_reconst1d_y( q_list(iv)%local(n)%val, &
            tmp3d, this%Minv_Ml_tr, this%Minv_Mc_tr, this%Minv_Mr_tr, this%IntrpMat, elem3d%Nnode_h1D, elem3d%Nnode_h1D, elem3d%Nnode_v, lmesh3d%Ne, this%Nnode_h1D_reconst, lmesh3d )
        case (filter_optrtype_reconstruct2)
          call apply_reconst1d_y_2( q_list(iv)%local(n)%val, &
            tmp3d, this%Ml_tr, this%Mc_tr, this%Mr_tr, elem3d%Nnode_h1D, elem3d%Nnode_h1D, elem3d%Nnode_v, lmesh3d%Ne, this%Nnode_h1D_reconst, lmesh3d )
        end select
        
      end do
      deallocate( tmp3d )
    end do
    
    return
  end subroutine meshfieldfilteroperation3d_apply_filter
  
!- Private -----------------------
 
!OCL SERIAL
  subroutine extract_tmp3d( tmp3D, q0, q0_, lmesh3D, elem3D, vmapP, hHaloSize, x_or_y )
    implicit none
    class(localmesh3d), intent(in) :: lmesh3D
    class(elementbase3d), intent(in) :: elem3D
    integer, intent(in) :: hHaloSize
    real(RP), intent(out) :: tmp3D(-hHaloSize+1:elem3D%Nnode_h1D+hHaloSize,-hHaloSize+1:elem3D%Nnode_h1D+hHaloSize,elem3D%Nnode_v,lmesh3D%Ne)
    real(RP), intent(in) :: q0(elem3D%Np*lmesh3D%NeA) 
    real(RP), intent(in) :: q0_(elem3D%Nnode_h1D,elem3D%Nnode_h1D,elem3D%Nnode_v,lmesh3D%NeA) 
    integer, intent(in) :: vmapP(elem3D%NfpTot,lmesh3D%Ne) 
    integer, intent(in) :: x_or_y !< x: 0, y: 1
 
    integer :: ke, ke_h, ke_x, ke_y, ke_z, px, py, ph, pz, i, f
    integer :: fso, fs, fe
    integer :: iP(elem3D%Nnode_h1D)
    real(RP) :: halo_h(elem3D%Nnode_h1D,elem3D%Nnode_v,hHaloSize,max(lmesh3D%NeX,lmesh3D%NeY),lmesh3D%NeZ,4)
    integer :: Neh_f(4)
    integer :: NehxNeZ_os(4)
    !------------------------------
 
    neh_f(:) = (/ lmesh3d%NeX, lmesh3d%NeY, lmesh3d%NeX, lmesh3d%NeY /)
    nehxnez_os(1) = 0
    do f=2, 4
      nehxnez_os(f) = nehxnez_os(f-1) + neh_f(f-1) * lmesh3d%NeZ 
    end do
    
    !$omp parallel private(ke, ke_x, ke_y, ke_h, ke_z, i, px, py, ph, pz, fso, fs, fe, iP, f)
 
    !$omp do collapse(2)
    do f=1, 4
      do ke_z=1, lmesh3d%NeZ
      do ke_h=1, neh_f(f)
        fso = lmesh3d%Ne * elem3d%Np &
            + elem3d%Nfp_h*hhalosize * ( nehxnez_os(f) + (ke_h-1) + (ke_z-1)*neh_f(f) )
 
        do ph=1, hhalosize
        do pz=1, elem3d%Nnode_v
          fs = fso + 1 + (pz-1)*elem3d%Nnode_h1D + (ph-1)*elem3d%Nfp_h
          fe = fso + elem3d%Nnode_h1D - 1
          halo_h(:,pz,ph,ke_h,ke_z,f) = q0(fs:fe)
        end do
        end do
      end do
      end do
    end do
 
    !$omp do collapse(3)
    do ke_z=1, lmesh3d%NeZ
    do ke_y=1, lmesh3d%NeY
    do ke_x=1, lmesh3d%NeX
      ke = ke_x + (ke_y-1) * lmesh3d%NeX + (ke_z-1) * lmesh3d%NeX * lmesh3d%NeY
 
      do pz=1, elem3d%Nnode_v
      do py=1, elem3d%Nnode_h1D
      do px=1, elem3d%Nnode_h1D
        tmp3d(px,py,pz,ke) = q0_(px,py,pz,ke)
      end do
      end do
      end do
 
      if ( x_or_y == 1 ) then ! y-direction
        ! Face 1
        fso = 0
        do pz=1, elem3d%Nnode_v
          fs = fso + 1 + (pz-1)*elem3d%Nnode_h1D; fe = fs + elem3d%Nnode_h1D - 1
          ip(:) = vmapp(fs:fe,ke)
          if ( ip(1) <= lmesh3d%Ne * elem3d%Np ) then
            do ph=1, hhalosize
              tmp3d(1:elem3d%Nnode_h1D,-ph+1,pz,ke) = q0(ip(:) - (ph-1)*elem3d%Nnode_h1D)
            end do
          else
            do ph=1, hhalosize
              tmp3d(1:elem3d%Nnode_h1D,-ph+1,pz,ke) = halo_h(:,pz,ph,ke_x,ke_z,1)
            end do
          end if
        enddo
 
        ! Face 3
        fso = 2 * elem3d%Nfp_h
        do pz=1, elem3d%Nnode_v
          fs = fso + 1 + (pz-1)*elem3d%Nnode_h1D; fe = fs + elem3d%Nnode_h1D - 1
          ip(:) = vmapp(fs:fe,ke)  
          if ( ip(1) <= lmesh3d%Ne * elem3d%Np ) then
            do ph=1, hhalosize
              tmp3d(1:elem3d%Nnode_h1D,ph+elem3d%Nnode_h1D,pz,ke) = q0(ip(:) + (ph-1)*elem3d%Nnode_h1D)
            end do
          else
            do ph=1, hhalosize
              tmp3d(1:elem3d%Nnode_h1D,ph+elem3d%Nnode_h1D,pz,ke) = halo_h(:,pz,ph,ke_x,ke_z,3)
            end do
          end if          
        enddo
      else if ( x_or_y == 0 ) then ! x-direction
        ! Face 2
        fso = elem3d%Nfp_h
        do pz=1, elem3d%Nnode_v
          fs = fso + 1 + (pz-1)*elem3d%Nnode_h1D; fe = fs + elem3d%Nnode_h1D - 1
          ip(:) = vmapp(fs:fe,ke)  
          if ( ip(1) <= lmesh3d%Ne * elem3d%Np ) then
            do ph=1, hhalosize
              tmp3d(ph+elem3d%Nnode_h1D,1:elem3d%Nnode_h1D,pz,ke) = q0(ip(:) + (ph-1))
            end do
          else
            do ph=1, hhalosize
              tmp3d(ph+elem3d%Nnode_h1D,1:elem3d%Nnode_h1D,pz,ke) = halo_h(:,pz,ph,ke_y,ke_z,2)
            end do
          end if           
        enddo
 
        ! Face 4
        fso = 3 * elem3d%Nfp_h
        do pz=1, elem3d%Nnode_v
          fs = fso + 1 + (pz-1)*elem3d%Nnode_h1D; fe = fs + elem3d%Nnode_h1D - 1
          ip(:) = vmapp(fs:fe,ke)  
          if ( ip(1) <= lmesh3d%Ne * elem3d%Np ) then
            do ph=1, hhalosize
              tmp3d(-ph+1,1:elem3d%Nnode_h1D,pz,ke) = q0(ip(:) - (ph-1))
            end do
          else
            do ph=1, hhalosize
              tmp3d(-ph+1,1:elem3d%Nnode_h1D,pz,ke) = halo_h(:,pz,ph,ke_y,ke_z,4)
            end do
          end if           
        enddo
      end if
    
    end do
    end do
    end do
 
    !$omp end parallel
    return
  end subroutine extract_tmp3d
  
end module scale_meshfield_filter_operation_3d