scale_element_operation_general.F90

Go to the documentation of this file.

 
!-------------------------------------------------------------------------------
!> module FElib / Element / Operation with arbitary elements
!!
!! @par Description
!!           A module for providing mathematical operations with arbitary elements using a module for SpMV
!!
!! @author Yuta Kawai, Xuanzhengbo Ren, and Team SCALE
!!
!<
#include "scaleFElib.h"
module scale_element_operation_general
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
 
  use scale_element_base, only: &
    elementbase3d, &
    elementbase3d_init, elementbase3d_final
 
  use scale_sparsemat, only: &
    sparsemat, sparsemat_matmul
  use scale_element_modalfilter, only: modalfilter  
 
  use scale_element_operation_base, only: elementoperationbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !  
 
  !> Derived type for elementwise operations with arbitary elements
  type, public, extends(elementoperationbase3d) :: elementoperationgeneral
    type(sparsemat), pointer :: dx_sm
    type(sparsemat), pointer :: dy_sm
    type(sparsemat), pointer :: dz_sm
    type(sparsemat), pointer :: lift_sm
    real(rp), allocatable :: intrpmat_vpordm1(:,:)
 
    type(modalfilter) :: mfilter
    type(modalfilter) :: mfilter_tracer
  contains
    procedure, public :: init => element_operation_general_init
    procedure, public :: final => element_operation_general_final
    procedure, public :: dx => element_operation_general_dx
    procedure, public :: dy => element_operation_general_dy
    procedure, public :: dz => element_operation_general_dz
    procedure, public :: lift => element_operation_general_lift
    procedure, public :: dxdydzlift => element_operation_general_dxdydzlift
    procedure, public :: div => element_operation_general_div
    procedure, public :: div_var5 => element_operation_general_div_var5
    procedure, public :: div_var5_2 => element_operation_general_div_var5_2
    procedure, public :: lift_var5 => element_operation_general_lift_var5
    procedure, public :: vfilterpm1 => element_operation_general_vfilterpm1
    !-
    procedure, public :: setup_modalfilter => element_operation_general_setup_modalfilter
    procedure, public :: setup_modalfilter_tracer => element_operation_general_setup_modalfilter_tracer
    procedure, public :: modalfilter_tracer => element_operation_general_modalfilter_tracer
    procedure, public :: modalfilter_var5 => element_operation_general_modalfilter_var5
  end type elementoperationgeneral
  
  interface elementoperationgeneral_generate_vpordm1
    module procedure element_operation_general_generate_vpordm1
  end interface
  public :: elementoperationgeneral_generate_vpordm1
  
contains
 
  !> Initialization
  !!
!OCL SERIAL
  subroutine element_operation_general_init( this, elem3D, &
      Dx, Dy, Dz, Lift )
    use scale_element_operation_base, only: element_operator_type_general
    implicit none
    class(ElementOperationGeneral), intent(inout) :: this
    class(ElementBase3D), intent(in), target :: elem3D
    type(SparseMat), intent(in), target :: Dx
    type(SparseMat), intent(in), target :: Dy
    type(SparseMat), intent(in), target :: Dz
    type(SparseMat), intent(in), target :: Lift
    !----------------------------------------------------------
 
    this%operator_type = element_operator_type_general
    
    this%elem3D => elem3d
    this%Dx_sm => dx
    this%Dy_sm => dy
    this%Dz_sm => dz
    this%Lift_sm => lift
 
    !--
    allocate( this%IntrpMat_VPOrdM1(elem3d%Np,elem3d%Np) )
    call element_operation_general_generate_vpordm1( this%IntrpMat_VPOrdM1, &
      elem3d )
 
    return
  end subroutine element_operation_general_init
 
  !> Generate a vertical filter matrix to remove the highest mode
  !!
!OCL SERIAL
  subroutine element_operation_general_generate_vpordm1( IntrpMat_VPOrdM1, &
    elem3D )
    implicit none
    class(elementbase3d), intent(in), target :: elem3D
    real(RP), intent(out) :: IntrpMat_VPOrdM1(elem3D%Np,elem3D%Np)
    !----------------------------------------------------------
    call elem3d%Generate_ModalTruncationMat( elem3d%PolyOrder_h, elem3d%PolyOrder_v-1, & ! (in)
      intrpmat_vpordm1 ) ! (out)
    return
  end subroutine element_operation_general_generate_vpordm1
 
  !> Setup modal filter
  !!
!OCL SERIAL
  subroutine element_operation_general_setup_modalfilter( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
 
    implicit none
    class(elementoperationgeneral), intent(inout) :: this
    real(RP), intent(in) :: MF_ETAC_h
    real(RP), intent(in) :: MF_ALPHA_h
    integer, intent(in) :: MF_ORDER_h
    real(RP), intent(in) :: MF_ETAC_v
    real(RP), intent(in) :: MF_ALPHA_v
    integer, intent(in) :: MF_ORDER_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter, &
      mf_etac_h, mf_alpha_h, mf_order_h,         &
      mf_etac_v, mf_alpha_v, mf_order_v,         &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v  )
  
    return
  end subroutine element_operation_general_setup_modalfilter  
 
  !> Setup modal filter for tracer
  !!
!OCL SERIAL
  subroutine element_operation_general_setup_modalfilter_tracer( this, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v )
    implicit none
    class(elementoperationgeneral), intent(inout) :: this
    real(RP), intent(in) :: MF_ETAC_h
    real(RP), intent(in) :: MF_ALPHA_h
    integer, intent(in) :: MF_ORDER_h
    real(RP), intent(in) :: MF_ETAC_v
    real(RP), intent(in) :: MF_ALPHA_v
    integer, intent(in) :: MF_ORDER_v
    !--------------------------------------------------------
 
    call setup_modalfilter( this%MFilter_tracer, &
      mf_etac_h, mf_alpha_h, mf_order_h,         &
      mf_etac_v, mf_alpha_v, mf_order_v,         &
      this%elem3D%PolyOrder_h, this%elem3D%PolyOrder_v  )
    
    return
  end subroutine element_operation_general_setup_modalfilter_tracer
  
 
  !> Finalization
  !!
  !OCL SERIAL
  subroutine element_operation_general_final( this )
    implicit none
    class(elementoperationgeneral), intent(inout) :: this
    !----------------------------------------------------------
 
    nullify( this%elem3D )
    nullify( this%Dx_sm, this%Dy_sm, this%Dz_sm, this%Lift_sm )
 
    deallocate( this%IntrpMat_VPOrdM1 )
 
    return
  end subroutine element_operation_general_final
 
!> Calculate the differential in x-direction
!!
!OCL SERIAL
  subroutine element_operation_general_dx( this, vec_in, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np)
    real(RP), intent(out) :: vec_out(this%elem3D%Np)
    !----------------------------------------------------------
    call sparsemat_matmul( this%Dx_sm, vec_in, vec_out )
    return
  end subroutine element_operation_general_dx
 
!> Calculate the differential in y-direction
!!
!OCL SERIAL
  subroutine element_operation_general_dy( this, vec_in, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np)
    real(RP), intent(out) :: vec_out(this%elem3D%Np)
    !----------------------------------------------------------
    call sparsemat_matmul( this%Dy_sm, vec_in, vec_out )
    return
  end subroutine element_operation_general_dy
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_general_dz( this, vec_in, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np)
    real(RP), intent(out) :: vec_out(this%elem3D%Np)
    !----------------------------------------------------------
    call sparsemat_matmul( this%Dz_sm, vec_in, vec_out )
    return
  end subroutine element_operation_general_dz  
 
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_general_lift( this, vec_in, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%NfpTot)
    real(RP), intent(out) :: vec_out(this%elem3D%Np)
    !----------------------------------------------------------
    call sparsemat_matmul( this%Lift_sm, vec_in, vec_out )
    return
  end subroutine element_operation_general_lift 
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_general_dxdydzlift( this, vec_in, vec_in_lift, vec_out_dx, vec_out_dy, vec_out_dz, vec_out_lift )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np)
    real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
    real(RP), intent(out) :: vec_out_dx(this%elem3D%Np)
    real(RP), intent(out) :: vec_out_dy(this%elem3D%Np)
    real(RP), intent(out) :: vec_out_dz(this%elem3D%Np)
    real(RP), intent(out) :: vec_out_lift(this%elem3D%Np)
    !----------------------------------------------------------
 
    call sparsemat_matmul( this%Dx_sm, vec_in, vec_out_dx )
    call sparsemat_matmul( this%Dy_sm, vec_in, vec_out_dy )
    call sparsemat_matmul( this%Dz_sm, vec_in, vec_out_dz )
    call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out_lift )
    return
  end subroutine element_operation_general_dxdydzlift
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_general_div( this, vec_in, vec_in_lift, &
    vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np,3)
    real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot)
    real(RP), intent(out) :: vec_out(this%elem3D%Np,4)
    !---------------------------------------------------------------
 
    call sparsemat_matmul( this%Dx_sm, vec_in(:,1), vec_out(:,1) )
    call sparsemat_matmul( this%Dy_sm, vec_in(:,2), vec_out(:,2) )
    call sparsemat_matmul( this%Dz_sm, vec_in(:,3), vec_out(:,3) )
    call sparsemat_matmul( this%Lift_sm, vec_in_lift, vec_out(:,4) )
    return
  end subroutine element_operation_general_div    
 
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_general_div_var5( this, vec_in, vec_in_lift, &
    vec_out_d )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np,3,5)
    real(RP), intent(in) :: vec_in_lift(this%elem3D%NfpTot,5)
    real(RP), intent(out) :: vec_out_d(this%elem3D%Np,4,5)
 
    integer :: iv
    !---------------------------------------------------------------
 
    do iv=1, 5
      call sparsemat_matmul( this%Dx_sm, vec_in(:,1,iv), vec_out_d(:,1,iv) )
      call sparsemat_matmul( this%Dy_sm, vec_in(:,2,iv), vec_out_d(:,2,iv) )
      call sparsemat_matmul( this%Dz_sm, vec_in(:,3,iv), vec_out_d(:,3,iv) )
    end do
    do iv=1, 5
      call sparsemat_matmul( this%Lift_sm, vec_in_lift(:,iv), vec_out_d(:,4,iv) )
    end do
    return
  end subroutine element_operation_general_div_var5
 
!> Calculate the 3D gradient
!!
!OCL SERIAL
  subroutine element_operation_general_div_var5_2( this, vec_in,  &
    vec_out_d )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np,3,5)
    real(RP), intent(out) :: vec_out_d(this%elem3D%Np,3,5)
 
    integer :: iv
    !---------------------------------------------------------------
 
    do iv=1, 5
      call sparsemat_matmul( this%Dx_sm, vec_in(:,1,iv), vec_out_d(:,1,iv) )
      call sparsemat_matmul( this%Dy_sm, vec_in(:,2,iv), vec_out_d(:,2,iv) )
      call sparsemat_matmul( this%Dz_sm, vec_in(:,3,iv), vec_out_d(:,3,iv) )
    end do
    return
  end subroutine element_operation_general_div_var5_2
!> Calculate the differential in z-direction
!!
!OCL SERIAL
  subroutine element_operation_general_lift_var5( this, vec_in, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%NfpTot,5)
    real(RP), intent(out) :: vec_out(this%elem3D%Np,5)
 
    integer :: iv
    !----------------------------------------------------------
    do iv=1,5
      call sparsemat_matmul( this%Lift_sm, vec_in(:,iv), vec_out(:,iv) )
    end do
    return
  end subroutine element_operation_general_lift_var5
 
!OCL SERIAL
  subroutine element_operation_general_vfilterpm1( this, vec_in, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np)
    real(RP), intent(out) :: vec_out(this%elem3D%Np)
    !---------------------------------------------------------------
 
    call matmul_( this%IntrpMat_VPOrdM1, vec_in, this%elem3D%Np, &
      vec_out )    
    return
  end subroutine element_operation_general_vfilterpm1 
!--
!OCL SERIAL
  subroutine matmul_( IntrpMat_VPOrdM1, vec_in_, Np, vec_out_ )
    implicit none
    integer, intent(in) :: Np
    real(RP), intent(in) :: IntrpMat_VPOrdM1(Np,Np)
    real(RP), intent(in) :: vec_in_(Np)
    real(RP), intent(out) :: vec_out_(Np)
    !-------------------------------------------
    vec_out_(:) = matmul( intrpmat_vpordm1(:,:), vec_in_(:) )
    return
  end subroutine matmul_
 
!OCL SERIAL
  subroutine element_operation_general_modalfilter_tracer( this, vec_in, vec_work, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np)
    real(RP), intent(out) :: vec_work(this%elem3D%Np)
    real(RP), intent(out) :: vec_out(this%elem3D%Np)
 
    integer :: ii, kk, Np
    real(RP) :: Mik
    !---------------------------------------------
 
    np = this%elem3D%Np 
    vec_out(:) = 0.0_rp
 
    do ii=1, np
    do kk=1, np
      mik = this%MFilter_tracer%FilterMat(ii,kk)
      vec_out(ii) = vec_out(ii) + mik * vec_in(kk)
    end do
    end do
 
    return
  end subroutine element_operation_general_modalfilter_tracer
 
!OCL SERIAL
  subroutine element_operation_general_modalfilter_var5( this, vec_in, vec_work, vec_out )
    implicit none
    class(elementoperationgeneral), intent(in) :: this
    real(RP), intent(in) :: vec_in(this%elem3D%Np,5)
    real(RP), intent(out) :: vec_work(this%elem3D%Np)
    real(RP), intent(out) :: vec_out(this%elem3D%Np,5)
 
    integer :: ii, kk, Np
    real(RP) :: Mik
    !---------------------------------------------
 
    np = this%elem3D%Np 
    vec_out(:,:) = 0.0_rp
 
    do ii=1, np
    do kk=1, np
      mik = this%MFilter%FilterMat(ii,kk)
 
      vec_out(ii,1) = vec_out(ii,1) + mik * vec_in(kk,1)
      vec_out(ii,2) = vec_out(ii,2) + mik * vec_in(kk,2)
      vec_out(ii,3) = vec_out(ii,3) + mik * vec_in(kk,3)
      vec_out(ii,4) = vec_out(ii,4) + mik * vec_in(kk,4)
      vec_out(ii,5) = vec_out(ii,5) + mik * vec_in(kk,5)
    end do
    end do
 
    return
  end subroutine element_operation_general_modalfilter_var5
 
!- private -
 
!OCL SERIAL
  subroutine setup_modalfilter( MFilter, &
    MF_ETAC_h, MF_ALPHA_h, MF_ORDER_h, &
    MF_ETAC_v, MF_ALPHA_v, MF_ORDER_v, &
    PolyOrder_h, PolyOrder_v )
 
    use scale_element_hexahedral, only: hexahedralelement
    implicit none
 
    class(modalfilter), intent(inout) :: MFilter
    real(RP), intent(in) :: MF_ETAC_h
    real(RP), intent(in) :: MF_ALPHA_h
    integer, intent(in) :: MF_ORDER_h
    real(RP), intent(in) :: MF_ETAC_v
    real(RP), intent(in) :: MF_ALPHA_v
    integer, intent(in) :: MF_ORDER_v
    integer, intent(in) :: PolyOrder_h
    integer, intent(in) :: PolyOrder_v
 
    type(hexahedralelement) :: elem3D
    !--------------------------------------------------------
 
    call elem3d%Init( polyorder_h, polyorder_v, .false. )
 
    call mfilter%Init( &
      elem3d,                              & ! (in)
      mf_etac_h, mf_alpha_h, mf_order_h,   & ! (in)
      mf_etac_v, mf_alpha_v, mf_order_v    ) ! (in)
    
    call elem3d%Final()
    return
  end subroutine setup_modalfilter
 
end module scale_element_operation_general