scale_meshfield_filter_operation_1d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Data / Filter operation 1D
!!
!! @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_1d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc, only: prc_abort
 
  use scale_element_base, only: elementbase1d
  use scale_element_modalfilter, only: modalfilter
  use scale_mesh_base1d, only: meshbase1d
  use scale_localmesh_base, only: localmeshbase
  use scale_localmesh_1d, only: localmesh1d
 
  use scale_meshfield_base, only: &
    meshfield1d
  use scale_localmeshfield_base, only: localmeshfieldbase
  use scale_meshfieldcomm_1d, only: &
    meshfieldcomm1d
  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,             &
    filter_optrtype_reconstruct2_gl, &
    filter_optrtype_interface_correction, &
    apply_filter1d_x => meshfieldfilteroperationbase_apply_filter1d_x,       &
    apply_reconst1d_x => meshfieldfilteroperationbase_apply_reconst1d_x,     &
    apply_reconst1d_x_2 => meshfieldfilteroperationbase_apply_reconst1d_x_2, &
    apply_interface_correction1d_x => meshfieldfilteroperationbase_apply_interface_correction1d_x
  
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !> Derived type to represent filter operation for 1D mesh field
  type, public, extends(meshfieldfilteroperationbase) :: meshfieldfilteroperation1d
    type(meshfieldcomm1d) :: vars_comm
 
    integer :: nnode_gl
  contains
    procedure :: init => meshfieldfilteroperation1d_init
    procedure :: final => meshfieldfilteroperation1d_final
    procedure :: apply => meshfieldfilteroperation1d_apply_filter
    procedure :: apply_gl_lc => meshfieldfilteroperation1d_apply_filter_gl_lc
  end type meshfieldfilteroperation1d
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private type & procedure
  !
  private :: extract_tmp1d
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
contains
  !> Initialize an object to represent filter operation for 1D mesh field
!OCL SERIAL
  subroutine meshfieldfilteroperation1d_init( this, &
    FilterOptrType,                  &
    FilterShape, FilterWidthFac,     &
    Nnode_1D_reconst,                &
    sfield_num,                      &
    mesh1D, Nnode_1D_GL, IF_r )
    implicit none
    class(meshfieldfilteroperation1d), intent(inout) :: this
    class(meshbase1d), intent(in) :: mesh1D
    character(*), intent(in) :: FilterOptrType
    character(*), intent(in) :: FilterShape
    real(RP), intent(in) :: FilterWidthFac
    integer, intent(in) :: Nnode_1D_reconst
    integer, intent(in) :: sfield_num
    integer, intent(in), optional :: Nnode_1D_GL
    integer, intent(in), optional :: IF_r
    !------------------------------------------
 
    call meshfieldfilteroperationbase_init( this, mesh1d%refElem1D%Np, &
      filteroptrtype, filtershape, filterwidthfac, &
      nnode_1d_reconst, nnode_1d_gl, if_r )
    
    if (present(nnode_1d_gl)) then
      this%Nnode_GL = nnode_1d_gl
    else
      this%Nnode_GL = mesh1d%refElem1D%Np
    end if
 
    select type(mesh1d)
    class is (meshbase1d)
      call this%vars_comm%Init( sfield_num, 0, mesh1d, &
        halosize_1d=this%hHaloSize )
    class default
      log_info('MeshFieldFilterOperation1D_Init',*) "Unsupported mesh type is specified. Check!"
      call prc_abort
    end select
    return
  end subroutine meshfieldfilteroperation1d_init
 
  !> Finalize an object to represent filter operation for 1D mesh field
!OCL SERIAL
  subroutine meshfieldfilteroperation1d_final( this )
    implicit none
    class(meshfieldfilteroperation1d), intent(inout) :: this
    !------------------------------------------
    call meshfieldfilteroperationbase_final( this )
    call this%vars_comm%Final()
    return
  end subroutine meshfieldfilteroperation1d_final
 
  !> Apply filter operation to 1D mesh field data
!OCL serial
  subroutine meshfieldfilteroperation1d_apply_filter( this, q_list, &
    mesh1D )
    implicit none
    class(meshfieldfilteroperation1d), intent(inout) :: this
    type(meshfield1d), intent(inout), target :: q_list(this%vars_comm%field_num_tot)
    class(meshbase1d), intent(in), target :: mesh1D
 
    class(localmesh1d), pointer :: lmesh1D
    class(elementbase1d), pointer :: elem1D
    real(RP), allocatable :: tmp1D(:,:,:)
 
    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)%field1D => 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, mesh1d%LOCAL_MESH_NUM
      lmesh1d => mesh1d%lcmesh_list(n)
      elem1d => lmesh1d%refElem1D
      allocate( tmp1d(-this%hHaloSize+1:elem1d%Np+this%hHaloSize,0:2,lmesh1d%Ne) )
 
      do iv=1, this%vars_comm%field_num_tot
        call extract_tmp1d( tmp1d, &
          q_list(iv)%local(n)%val, q_list(iv)%local(n)%val, lmesh1d, elem1d, lmesh1d%VMapP, this%hHaloSize )
        
        select case (this%operator_type)
        case (filter_optrtype_convfilter)
          call apply_filter1d_x( q_list(iv)%local(n)%val, &
            tmp1d, this%FilterMat_h1D, elem1d%Np, 1, 1, lmesh1d%Ne, lmesh1d%NeA, elem1d%Np )
        case (filter_optrtype_reconstruct)
          call apply_reconst1d_x( q_list(iv)%local(n)%val, &
            tmp1d, this%Minv_Ml_tr, this%Minv_Mc_tr, this%Minv_Mr_tr, this%IntrpMat, elem1d%Np, 1, 1, lmesh1d%Ne, this%Nnode_h1D_reconst, lmesh1d )
        case (filter_optrtype_reconstruct2)
          call apply_reconst1d_x_2( q_list(iv)%local(n)%val, &
            tmp1d, this%Ml_tr, this%Mc_tr, this%Mr_tr, elem1d%Np, 1, 1, lmesh1d%Ne, this%Nnode_h1D_reconst, lmesh1d )
        case (filter_optrtype_interface_correction)
          call apply_interface_correction1d_x( q_list(iv)%local(n)%val, &
            tmp1d, this%IF_gL, this%IF_gR, elem1d%Np, 1, 1, lmesh1d%Ne, lmesh1d )
        end select
 
      end do
      deallocate( tmp1d )
    end do    
    return
  end subroutine meshfieldfilteroperation1d_apply_filter
 
 
  !> Apply filter operation to 1D mesh field data
!OCL serial
  subroutine meshfieldfilteroperation1d_apply_filter_gl_lc( this, &
    q_GL,   &
    q_list, &
    lmesh, elem1D  )
    implicit none
    class(meshfieldfilteroperation1d), intent(inout) :: this
    class(localmesh1d), intent(in) :: lmesh
    class(elementbase1d), intent(in) :: elem1D
    real(RP), intent(out) :: q_GL(this%Nnode_GL,lmesh%Ne,this%vars_comm%field_num_tot)
    type(meshfield1d), intent(inout), target :: q_list(this%vars_comm%field_num_tot)
 
    real(RP), allocatable :: tmp1D(:,:,:)
 
    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)%field1D => 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)
 
    allocate( tmp1d(-this%hHaloSize+1:elem1d%Np+this%hHaloSize,0:2,lmesh%Ne) )
 
    n = lmesh%lcdomID
    do iv=1, this%vars_comm%field_num_tot
      call extract_tmp1d( tmp1d, &
        q_list(iv)%local(n)%val, q_list(iv)%local(n)%val, lmesh, elem1d, lmesh%VMapP, this%hHaloSize )
      
      select case (this%operator_type)
      case (filter_optrtype_reconstruct2_gl)
        call apply_reconst1d_x_2_gl( q_gl(:,:,iv), &
          tmp1d, this%Ml_tr, this%Mc_tr, this%Mr_tr, elem1d%Np, this%Nnode_GL, lmesh%Ne, this%Nnode_h1D_reconst, lmesh )
      end select
 
    end do
    return
  end subroutine meshfieldfilteroperation1d_apply_filter_gl_lc
  
!- Private -----------------------
 
!OCL SERIAL
  subroutine apply_reconst1d_x_2_gl( q, q0, Minv_Ml_tr, Minv_Mc_tr, Minv_Mr_tr, &
    Npx, NpxGL, Ne, Npx_reconst, lmesh )
    implicit none
    class(localmeshbase), intent(in) :: lmesh
    integer, intent(in) :: Npx, NpxGL, Ne
    integer, intent(in) :: Npx_reconst
    real(RP), intent(out) :: q(NpxGL,lmesh%Ne)
    real(RP), intent(in) :: q0(-Npx+1:2*Npx,0:2,Ne)
    real(RP), intent(in) :: Minv_Ml_tr(Npx,NpxGL)
    real(RP), intent(in) :: Minv_Mc_tr(Npx,NpxGL)
    real(RP), intent(in) :: Minv_Mr_tr(Npx,NpxGL)
 
    integer :: ke
 
    integer :: i,j
 
    real(RP) :: tmp(NpxGL)
    real(RP) :: s
    !-------------------------------------
 
    !$omp parallel do private(ke,i,j, tmp,s)
    do ke=1, ne
      do i=1, npxgl
        s = 0.0_rp
        do j=1, npx
          s = s &
            + minv_ml_tr(j,i) * q0(j-npx,1,ke) &
            + minv_mc_tr(j,i) * q0(    j,1,ke) &
            + minv_mr_tr(j,i) * q0(j+npx,1,ke)
        end do
        tmp(i) = s
      end do
      q(:,ke) = tmp(:)
    end do
    return
  end subroutine apply_reconst1d_x_2_gl
 
!OCL SERIAL
  subroutine extract_tmp1d( tmp1D, q0, q0_, lmesh1D, elem1D, vmapP, hHaloSize )
    implicit none
    class(localmesh1d), intent(in) :: lmesh1D
    class(elementbase1d), intent(in) :: elem1D
    integer, intent(in) :: hHaloSize
    real(RP), intent(out) :: tmp1D(-hHaloSize+1:elem1D%Np+hHaloSize,0:2,lmesh1D%Ne)
    real(RP), intent(in) :: q0(elem1D%Np*lmesh1D%NeA) 
    real(RP), intent(in) :: q0_(elem1D%Np,lmesh1D%NeA) 
    integer, intent(in) :: vmapP(elem1D%NfpTot,lmesh1D%Ne) 
 
    integer :: ke, i, f
    integer :: fso
    integer :: ph
    integer :: iP
    real(RP) :: halo_h(hHaloSize,2)
    !------------------------------
 
    do f=1, 2
      fso = lmesh1d%Ne * elem1d%Np + hhalosize * ( f-1 )
      do ph=1, hhalosize
        halo_h(ph,f) = q0(fso+ph)
      end do
    end do
    
    !$omp parallel private(ke, i, ph, iP)
    !$omp do
    do ke=1, lmesh1d%Ne
      tmp1d(1:elem1d%Np,1,ke) = q0_(:,ke)
 
      ! Face 2
      ip = vmapp(2,ke)  
      if ( ip <= lmesh1d%Ne * elem1d%Np ) then
        do ph=1, hhalosize
          tmp1d(ph+elem1d%Np,1,ke) = q0(ip + (ph-1))
        end do
      else
        do ph=1, hhalosize
          tmp1d(ph+elem1d%Np,1,ke) = halo_h(ph,2)
        end do
      end if           
 
      ! Face 1
      ip = vmapp(1,ke)  
      if ( ip <= lmesh1d%Ne * elem1d%Np ) then
        do ph=1, hhalosize
          tmp1d(-ph+1,1,ke) = q0(ip - (ph-1))
        end do
      else
        do ph=1, hhalosize
          tmp1d(-ph+1,1,ke) = halo_h(ph,1)
        end do
      end if         
    end do
 
    !$omp end parallel
    return
  end subroutine extract_tmp1d
end module scale_meshfield_filter_operation_1d