scale_element_modalfilter.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Element/ ModalFilter
!!
!! @par Description
!!      A module for modal filtering
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_element_modalfilter
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
 
  use scale_element_base, only: &
    elementbase1d, elementbase2d, elementbase3d
  
  use scale_element_line, only: lineelement
  use scale_element_quadrilateral, only: quadrilateralelement
  use scale_element_hexahedral, only: hexahedralelement
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type
  !
 
  !> Derived type representing a modal filter
  type, public :: modalfilter
    real(rp), allocatable :: filtermat(:,:)
  contains
    procedure :: init_line => modalfilter_init_line
    procedure :: init_quadrilateral => modalfilter_init_quadrilateral
    procedure :: init_hexahedral => modalfilter_init_hexahedral
    generic :: init => init_line, init_quadrilateral, init_hexahedral
    procedure :: final => modalfilter_final
  end type modalfilter
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
  private :: get_exp_filter
 
contains
  subroutine modalfilter_init_line( this,  & ! (inout)
    elem,                                  & ! (in)
    etac, alpha, ord,                      & ! (in)
    tend_flag                              ) ! (in)
 
    implicit none
    class(modalfilter), intent(inout) :: this
    class(lineelement), intent(in) :: elem
    real(RP), intent(in) :: etac
    real(RP), intent(in) :: alpha
    integer, intent(in) :: ord
    logical, intent(in), optional :: tend_flag
 
    real(RP) :: filter1D(elem%Np)
    integer :: p
    logical :: tend_flag_
    !----------------------------------------------------
    
    tend_flag_ = .false.
    if ( present(tend_flag) ) tend_flag_ =  tend_flag
 
    call get_exp_filter( filter1d,               & ! (out)
      etac, alpha, ord, elem%Np, elem%PolyOrder, & ! (in)
      tend_flag_ )                                 ! (in)
 
    allocate( this%FilterMat(elem%Np,elem%Np) )
    this%FilterMat(:,:) = 0.0_rp
    do p=1, elem%Np
      this%FilterMat(p,p) = filter1d(p)
    end do
    this%FilterMat(:,:) = matmul(this%FilterMat, elem%invV)
    this%FilterMat(:,:) = matmul(elem%V, this%FilterMat)
    !$acc enter data copyin(this%FilterMat)
    
    return
  end subroutine modalfilter_init_line
 
  subroutine modalfilter_init_quadrilateral( this,   & ! (inout)
    elem,                                            & ! (in)
    etac, alpha, ord,                                & ! (in)
    tend_flag                                        ) ! (in)
 
    implicit none
    class(modalfilter), intent(inout) :: this
    class(quadrilateralelement), intent(in) :: elem
    real(RP), intent(in) :: etac
    real(RP), intent(in) :: alpha
    integer, intent(in) :: ord
    logical, intent(in), optional :: tend_flag
 
    real(RP) :: filter1D(elem%Nfp)
    integer :: p1, p2
    integer :: l
    logical :: tend_flag_
    !----------------------------------------------------
    
    tend_flag_ = .false.
    if ( present(tend_flag) ) tend_flag_ =  tend_flag
 
    call get_exp_filter( filter1d,                & ! (out)
      etac, alpha, ord, elem%Nfp, elem%PolyOrder, & ! (in)
      tend_flag_ )                                  ! (in)
    
    allocate( this%FilterMat(elem%Np,elem%Np) )
    this%FilterMat(:,:) = 0.0_rp
    do p2=1, elem%Nfp
    do p1=1, elem%Nfp
      l = p1 + (p2-1)*elem%Nfp
      this%FilterMat(l,l) = filter1d(p1) * filter1d(p2)
    end do  
    end do
    this%FilterMat(:,:) = matmul(this%FilterMat, elem%invV)
    this%FilterMat(:,:) = matmul(elem%V, this%FilterMat)
    !$acc enter data copyin(this%FilterMat)
    return
  end subroutine modalfilter_init_quadrilateral
 
  subroutine modalfilter_init_hexahedral( this,  & ! (inout)
    elem,                                        & ! (in)
    etac_h, alpha_h, ord_h,                      & ! (in)
    etac_v, alpha_v, ord_v,                      & ! (in)
    tend_flag                                    ) ! (in)
 
    implicit none
    class(modalfilter), intent(inout) :: this
    class(hexahedralelement), intent(in) :: elem
    real(RP), intent(in) :: etac_h
    real(RP), intent(in) :: alpha_h
    integer, intent(in) :: ord_h
    real(RP), intent(in) :: etac_v
    real(RP), intent(in) :: alpha_v
    integer, intent(in) :: ord_v
    logical, intent(in), optional :: tend_flag
 
    real(RP) :: filter1D_h(elem%Nnode_h1D)
    real(RP) :: filter1D_v(elem%Nnode_v)
    integer :: p1, p2, p3
    integer :: l
    logical :: tend_flag_
    !----------------------------------------------------
    
    tend_flag_ = .false.
    if ( present(tend_flag) ) tend_flag_ =  tend_flag
 
    call get_exp_filter( filter1d_h,                            & ! (out)
      etac_h, alpha_h, ord_h, elem%Nnode_h1D, elem%PolyOrder_h, & ! (in)
      tend_flag_ )                                                ! (in)
    
    call get_exp_filter( filter1d_v,                           & ! (out)
      etac_v, alpha_v, ord_v, elem%Nnode_v, elem%PolyOrder_v,  & ! (in)
      tend_flag_ )                                               ! (in)
 
    allocate( this%FilterMat(elem%Np,elem%Np) )
    this%FilterMat(:,:) = 0.0_rp
    do p3=1, elem%Nnode_v
    do p2=1, elem%Nnode_h1D
    do p1=1, elem%Nnode_h1D
      l = p1 + (p2-1)*elem%Nnode_h1D + (p3-1)*elem%Nnode_h1D**2
      this%FilterMat(l,l) = filter1d_h(p1) * filter1d_h(p2) * filter1d_v(p3)
    end do  
    end do
    end do
    this%FilterMat(:,:) = matmul(this%FilterMat, elem%invV)
    this%FilterMat(:,:) = matmul(elem%V, this%FilterMat)
    !$acc enter data copyin(this%FilterMat)
    return
  end subroutine modalfilter_init_hexahedral
 
  subroutine modalfilter_final( this )
    implicit none
 
    class(modalfilter), intent(inout) :: this
    !--------------------------------------------
 
    if( allocated(this%FilterMat) ) then
      !$acc exit data delete(this%FilterMat)
      deallocate( this%FilterMat )
    end if
    return
  end subroutine modalfilter_final
 
!-- private --------------------------------------------------
 
  subroutine get_exp_filter( filter, &  
    etac, alpha, ord, Np, polyOrder, &
    tend_flag )
 
    implicit none
    integer, intent(in) :: Np
    real(RP), intent(out) :: filter(Np)
    real(RP), intent(in) :: etac  
    real(RP), intent(in) :: alpha
    integer , intent(in) :: ord
    integer , intent(in) :: polyOrder
    logical, intent(in) :: tend_flag
 
    integer :: p
    real(RP) :: eta
    !-----------------------------------------
 
    if ( tend_flag ) then
      filter(:) = 0.0_rp
    else
      filter(:) = 1.0_rp
    end if
 
    do p=1, np
      eta = dble(p-1)/dble(polyorder)
      if ( eta >  etac .and. p /= 1) then
        filter(p) = - alpha * ( ((eta - etac)/(1.0_rp - etac))**ord )
        if ( .not. tend_flag ) filter(p) = exp( filter(p) )
      end if
    end do
 
    return
  end subroutine get_exp_filter
 
end module scale_element_modalfilter