scale_meshfield_filter_operation_2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Data / Filter operation 2D
!!
!! @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_2d
  !-----------------------------------------------------------------------------
  !
  !++ 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_localmesh_base, only: localmeshbase
  use scale_localmesh_2d, only: localmesh2d
 
  use scale_meshfield_base, only: &
    meshfield2d
  use scale_localmeshfield_base, only: localmeshfieldbase
  use scale_mesh_rectdom2d, only: &
    meshrectdom2d
  use scale_meshfieldcomm_rectdom2d, only: &
    meshfieldcommrectdom2d
  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 2D mesh field
  type, public, extends(meshfieldfilteroperationbase) :: meshfieldfilteroperation2d
    type(meshfieldcommrectdom2d) :: vars_comm
  contains
    procedure :: init => meshfieldfilteroperation2d_init
    procedure :: final => meshfieldfilteroperation2d_final
    procedure :: apply => meshfieldfilteroperation2d_apply_filter
  end type meshfieldfilteroperation2d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  !
  private :: extract_tmp2d
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
  !> Initialize an object to represent filter operation for 2D mesh field
!OCL SERIAL
  subroutine meshfieldfilteroperation2d_init( this, &
    FilterOptrType,                  &
    FilterShape, FilterWidthFac,     &
    Nnode_h1D_reconst,               &
    sfield_num, hvfield_num, htensorfield_num, &
    mesh2D )
    implicit none
    class(meshfieldfilteroperation2d), intent(inout) :: this
    class(meshbase2d), intent(in) :: mesh2D
    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, mesh2d%refElem2D%Nfp, &
      filteroptrtype, filtershape, filterwidthfac, &
      nnode_h1d_reconst )
 
    select type(mesh2d)
    class is (meshrectdom2d)
      call this%vars_comm%Init( sfield_num, hvfield_num, htensorfield_num, mesh2d, &
        halosize_1d=this%hHaloSize )
    class default
      log_info('MeshFieldFilterOperation2D_Init',*) "Unsupported mesh type is specified. Check!"
      call prc_abort
    end select
    return
  end subroutine meshfieldfilteroperation2d_init
 
  !> Finalize an object to represent filter operation for 2D mesh field
!OCL SERIAL
  subroutine meshfieldfilteroperation2d_final( this )
    implicit none
    class(meshfieldfilteroperation2d), intent(inout) :: this
    !------------------------------------------
    call meshfieldfilteroperationbase_final( this )
    call this%vars_comm%Final()
    return
  end subroutine meshfieldfilteroperation2d_final
 
  !> Apply filter operation to 2D mesh field data
!OCL serial
  subroutine meshfieldfilteroperation2d_apply_filter( this, q_list, &
    mesh2D )
    implicit none
    class(meshfieldfilteroperation2d), intent(inout) :: this
    type(meshfield2d), intent(inout), target :: q_list(this%vars_comm%field_num_tot)
    class(meshbase2d), intent(in), target :: mesh2D
 
    class(localmesh2d), pointer :: lmesh2D
    class(elementbase2d), pointer :: elem2D
    real(RP), allocatable :: tmp2D(:,:,:)
 
    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)%field2d => 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, mesh2d%LOCAL_MESH_NUM
      lmesh2d => mesh2d%lcmesh_list(n)
      elem2d => lmesh2d%refElem2D
      allocate( tmp2d(-this%hHaloSize+1:elem2d%Nfp+this%hHaloSize,-this%hHaloSize+1:elem2d%Nfp+this%hHaloSize,lmesh2d%Ne) )
 
      do iv=1, this%vars_comm%field_num_tot
        call extract_tmp2d( tmp2d, &
          q_list(iv)%local(n)%val, q_list(iv)%local(n)%val, lmesh2d, elem2d, lmesh2d%VMapP, this%hHaloSize, 0 )
        
        select case (this%operator_type)
        case (filter_optrtype_convfilter)
          call apply_filter1d_x( q_list(iv)%local(n)%val, &
            tmp2d, this%FilterMat_h1D, elem2d%Nfp, elem2d%Nfp, 1, lmesh2d%Ne, lmesh2d%NeA, elem2d%Nfp )
        case (filter_optrtype_reconstruct)
          call apply_reconst1d_x( q_list(iv)%local(n)%val, &
            tmp2d, this%Minv_Ml_tr, this%Minv_Mc_tr, this%Minv_Mr_tr, this%IntrpMat, elem2d%Nfp, elem2d%Nfp, 1, lmesh2d%Ne, this%Nnode_h1D_reconst, lmesh2d )
        case (filter_optrtype_reconstruct2)
          call apply_reconst1d_x_2( q_list(iv)%local(n)%val, &
            tmp2d, this%Ml_tr, this%Mc_tr, this%Mr_tr, elem2d%Nfp, elem2d%Nfp, 1, lmesh2d%Ne, this%Nnode_h1D_reconst, lmesh2d )
        end select
 
      end do
      deallocate( tmp2d )
    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, mesh2d%LOCAL_MESH_NUM
      lmesh2d => mesh2d%lcmesh_list(n)
      elem2d => lmesh2d%refElem2D
      allocate( tmp2d(-this%hHaloSize+1:elem2d%Nfp+this%hHaloSize,-this%hHaloSize+1:elem2d%Nfp+this%hHaloSize,lmesh2d%Ne) )
 
      do iv=1, this%vars_comm%field_num_tot
        call extract_tmp2d( tmp2d, &
          q_list(iv)%local(n)%val, q_list(iv)%local(n)%val, lmesh2d, elem2d, lmesh2d%VMapP, this%hHaloSize, 1 )
 
        select case (this%operator_type)
        case (filter_optrtype_convfilter)
          call apply_filter1d_y( q_list(iv)%local(n)%val, &
            tmp2d, this%FilterMat_h1D, elem2d%Nfp, elem2d%Nfp, 1, lmesh2d%Ne, lmesh2d%NeA, elem2d%Nfp )
        case (filter_optrtype_reconstruct)
          call apply_reconst1d_y( q_list(iv)%local(n)%val, &
            tmp2d, this%Minv_Ml_tr, this%Minv_Mc_tr, this%Minv_Mr_tr, this%IntrpMat, elem2d%Nfp, elem2d%Nfp, 1, lmesh2d%Ne, this%Nnode_h1D_reconst, lmesh2d )
        case (filter_optrtype_reconstruct2)
          call apply_reconst1d_y_2( q_list(iv)%local(n)%val, &
            tmp2d, this%Ml_tr, this%Mc_tr, this%Mr_tr, elem2d%Nfp, elem2d%Nfp, 1, lmesh2d%Ne, this%Nnode_h1D_reconst, lmesh2d )
        end select
        
      end do
      deallocate( tmp2d )
    end do
    
    return
  end subroutine meshfieldfilteroperation2d_apply_filter
  
!- Private -----------------------
 
!OCL SERIAL
  subroutine extract_tmp2d( tmp2D, q0, q0_, lmesh2D, elem2D, vmapP, hHaloSize, x_or_y )
    implicit none
    class(localmesh2d), intent(in) :: lmesh2D
    class(elementbase2d), intent(in) :: elem2D
    integer, intent(in) :: hHaloSize
    real(RP), intent(out) :: tmp2D(-hHaloSize+1:elem2D%Nfp+hHaloSize,-hHaloSize+1:elem2D%Nfp+hHaloSize,lmesh2D%Ne)
    real(RP), intent(in) :: q0(elem2D%Np*lmesh2D%NeA) 
    real(RP), intent(in) :: q0_(elem2D%Nfp,elem2D%Nfp,lmesh2D%NeA) 
    integer, intent(in) :: vmapP(elem2D%NfpTot,lmesh2D%Ne) 
    integer, intent(in) :: x_or_y !< x: 0, y: 1
 
    integer :: ke, ke_h, ke_x, ke_y, px, py, ph, i, f
    integer :: fso, fs, fe
    integer :: iP(elem2D%Nfp)
    real(RP) :: halo_h(elem2D%Nfp,hHaloSize,max(lmesh2D%NeX,lmesh2D%NeY),4)
    integer :: Neh_f(4)
    integer :: Neh_os(4)
    !------------------------------
 
    neh_f(:) = (/ lmesh2d%NeX, lmesh2d%NeY, lmesh2d%NeX, lmesh2d%NeY /)
    neh_os(1) = 0
    do f=2, 4
      neh_os(f) = neh_os(f-1) + neh_f(f-1)
    end do
    
    !$omp parallel private(ke, ke_x, ke_y, ke_h, i, px, py, ph, fso, fs, fe, iP, f)
 
    !$omp do
    do f=1, 4
      do ke_h=1, neh_f(f)
        fso = lmesh2d%Ne * elem2d%Np &
            + elem2d%Nfp*hhalosize * ( neh_os(f) + (ke_h-1) )
 
        do ph=1, hhalosize
          fs = fso + 1 + (ph-1)*elem2d%Nfp
          fe = fso + elem2d%Nfp - 1
          halo_h(:,ph,ke_h,f) = q0(fs:fe)
        end do
      end do
    end do
 
    !$omp do collapse(2)
    do ke_y=1, lmesh2d%NeY
    do ke_x=1, lmesh2d%NeX
      ke = ke_x + (ke_y-1) * lmesh2d%NeX
 
      do py=1, elem2d%Nfp
      do px=1, elem2d%Nfp
        tmp2d(px,py,ke) = q0_(px,py,ke)
      end do
      end do
 
      if ( x_or_y == 1 ) then ! y-direction
        ! Face 1
        fso = 0
        fs = fso + 1; fe = fs + elem2d%Nfp - 1
        ip(:) = vmapp(fs:fe,ke)
        if ( ip(1) <= lmesh2d%Ne * elem2d%Np ) then
          do ph=1, hhalosize
            tmp2d(1:elem2d%Nfp,-ph+1,ke) = q0(ip(:) - (ph-1)*elem2d%Nfp)
          end do
        else
          do ph=1, hhalosize
            tmp2d(1:elem2d%Nfp,-ph+1,ke) = halo_h(:,ph,ke_x,1)
          end do
        end if
 
        ! Face 3
        fso = 2 * elem2d%Nfp
        fs = fso + 1; fe = fs + elem2d%Nfp - 1
        ip(:) = vmapp(fs:fe,ke)  
        if ( ip(1) <= lmesh2d%Ne * elem2d%Np ) then
          do ph=1, hhalosize
            tmp2d(1:elem2d%Nfp,ph+elem2d%Nfp,ke) = q0(ip(:) + (ph-1)*elem2d%Nfp)
          end do
        else
          do ph=1, hhalosize
            tmp2d(1:elem2d%Nfp,ph+elem2d%Nfp,ke) = halo_h(:,ph,ke_x,3)
          end do
        end if          
      else if ( x_or_y == 0 ) then ! x-direction
        ! Face 2
        fso = elem2d%Nfp
        fs = fso + 1; fe = fs + elem2d%Nfp - 1
        ip(:) = vmapp(fs:fe,ke)  
        if ( ip(1) <= lmesh2d%Ne * elem2d%Np ) then
          do ph=1, hhalosize
            tmp2d(ph+elem2d%Nfp,1:elem2d%Nfp,ke) = q0(ip(:) + (ph-1))
          end do
        else
          do ph=1, hhalosize
            tmp2d(ph+elem2d%Nfp,1:elem2d%Nfp,ke) = halo_h(:,ph,ke_y,2)
          end do
        end if           
 
        ! Face 4
        fso = 3 * elem2d%Nfp
        fs = fso + 1; fe = fs + elem2d%Nfp - 1
        ip(:) = vmapp(fs:fe,ke)  
        if ( ip(1) <= lmesh2d%Ne * elem2d%Np ) then
          do ph=1, hhalosize
            tmp2d(-ph+1,1:elem2d%Nfp,ke) = q0(ip(:) - (ph-1))
          end do
        else
          do ph=1, hhalosize
            tmp2d(-ph+1,1:elem2d%Nfp,ke) = halo_h(:,ph,ke_y,4)
          end do
        end if         
      end if
    end do
    end do
 
    !$omp end parallel
    return
  end subroutine extract_tmp2d
end module scale_meshfield_filter_operation_2d