scale_atm_dyn_dgm_modalfilter.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Common / Modal filter
!!
!! @par Description
!!      Modal filter for Atmospheric dynamical process. 
!!      The modal filter suppresses the numerical instability due to the aliasing errors. 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_modalfilter
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
 
  use scale_element_base, only: elementbase
  use scale_localmesh_base, only: localmeshbase
  use scale_element_operation_base, only: elementoperationbase3d
  use scale_element_operation_gpu_driver, only: elementoperationgpudriver
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_dyn_dgm_modalfilter_apply
  public :: atm_dyn_dgm_tracer_modalfilter_apply
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
contains
 
#ifndef _OPENACC
!OCL SERIAL
  subroutine atm_dyn_dgm_modalfilter_apply(  &
    DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_,     & ! (inout)
    lmesh, elem, elem_operation,             & ! (in) 
    do_weight_gsqrt )                          ! (in)
 
    implicit none
 
    class(localmeshbase), intent(in) :: lmesh
    class(elementbase), intent(in) :: elem    
    real(rp), intent(inout)  :: ddens_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: momx_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: momy_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: momz_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: drhot_(elem%np,lmesh%nea)
    class(elementoperationbase3d), intent(in) :: elem_operation
    logical, intent(in), optional :: do_weight_gsqrt
    
    integer :: ke
    real(rp) :: tmp(elem%np,5)
    real(rp) :: work(elem%np)
    real(rp) :: tmp_out(elem%np,5)
 
    integer :: kk
    logical :: do_weight_gsqrt_
    real(rp) :: rgsqrt(elem%np)
    !------------------------------------
 
    if ( present( do_weight_gsqrt ) ) then
      do_weight_gsqrt_ = do_weight_gsqrt
    else
      do_weight_gsqrt_ = .false.
    end if
 
    if ( do_weight_gsqrt_ ) then
      !$omp parallel do private( tmp, tmp_out, work, kk, RGsqrt )
      do ke=lmesh%NeS, lmesh%NeE
        do kk=1, elem%Np
          tmp(kk,1) = lmesh%Gsqrt(kk,ke) * ddens_(kk,ke)
          tmp(kk,2) = lmesh%Gsqrt(kk,ke) * momx_(kk,ke)
          tmp(kk,3) = lmesh%Gsqrt(kk,ke) * momy_(kk,ke)
          tmp(kk,4) = lmesh%Gsqrt(kk,ke) * momz_(kk,ke)
          tmp(kk,5) = lmesh%Gsqrt(kk,ke) * drhot_(kk,ke)
        end do
 
        call elem_operation%ModalFilter_var5( tmp, work, &
          tmp_out ) ! (out)
 
        rgsqrt(:) = 1.0_rp / lmesh%Gsqrt(:,ke)
        ddens_(:,ke) = tmp_out(:,1) * rgsqrt(:)
        momx_(:,ke) = tmp_out(:,2) * rgsqrt(:)
        momy_(:,ke) = tmp_out(:,3) * rgsqrt(:)
        momz_(:,ke) = tmp_out(:,4) * rgsqrt(:)
        drhot_(:,ke) = tmp_out(:,5) * rgsqrt(:)
      end do    
 
    else
 
      !$omp parallel do private( tmp, work, kk )
      do ke=lmesh%NeS, lmesh%NeE
 
        do kk=1, elem%Np
          tmp(kk,1) = ddens_(kk,ke)
          tmp(kk,2) = momx_(kk,ke)
          tmp(kk,3) = momy_(kk,ke)
          tmp(kk,4) = momz_(kk,ke)
          tmp(kk,5) = drhot_(kk,ke)
        end do
        
        call elem_operation%ModalFilter_var5( tmp, work, &
          tmp_out )
 
        ddens_(:,ke) = tmp(:,1)
        momx_(:,ke) = tmp(:,2)
        momy_(:,ke) = tmp(:,3)
        momz_(:,ke) = tmp(:,4)
        drhot_(:,ke) = tmp(:,5)
      end do
 
    end if
 
    return
  end subroutine atm_dyn_dgm_modalfilter_apply
#else
  subroutine atm_dyn_dgm_modalfilter_apply(  &
    DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_,     & ! (inout)
    lmesh, elem, elem_operation,             & ! (in) 
    do_weight_gsqrt )                          ! (in)
 
    use scale_element_operation_gpu_driver, only: elementoperationgpu_modalfilter_var5
    implicit none
 
    class(localmeshbase), intent(in) :: lmesh
    class(elementbase), intent(in) :: elem    
    real(rp), intent(inout)  :: ddens_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: momx_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: momy_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: momz_(elem%np,lmesh%nea)
    real(rp), intent(inout)  :: drhot_(elem%np,lmesh%nea)
    class(elementoperationbase3d), intent(in) :: elem_operation
    logical, intent(in), optional :: do_weight_gsqrt
    
    integer :: ke
    real(rp) :: tmp(elem%np,lmesh%ne,5)
    real(rp) :: work(elem%np,lmesh%ne)
    real(rp) :: tmp_out(elem%np,lmesh%ne,5)
    type(elementoperationgpudriver) :: elem_operation_gpu_driver
 
    integer :: kk
    logical :: do_weight_gsqrt_
    real(rp) :: rgsqrt
    real(rp) :: gsqrt_
    !------------------------------------
 
    if ( present( do_weight_gsqrt ) ) then
      do_weight_gsqrt_ = do_weight_gsqrt
    else
      do_weight_gsqrt_ = .false.
    end if
 
 
    call elem_operation_gpu_driver%Init( elem_operation )
    !$acc data create( tmp, work, tmp_out )
 
    if ( do_weight_gsqrt_ ) then
      !$acc parallel loop collapse(2) present( DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, lmesh, elem ) async(1)
      do ke=lmesh%NeS, lmesh%NeE
      do kk=1, elem%Np
        gsqrt_ = lmesh%Gsqrt(kk,ke)
        tmp(kk,ke,1) = gsqrt_ * ddens_(kk,ke)
        tmp(kk,ke,2) = gsqrt_ * momx_(kk,ke)
        tmp(kk,ke,3) = gsqrt_ * momy_(kk,ke)
        tmp(kk,ke,4) = gsqrt_ * momz_(kk,ke)
        tmp(kk,ke,5) = gsqrt_ * drhot_(kk,ke)
      end do
      end do
    else
      !$acc parallel loop collapse(2) present( DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, lmesh, elem ) async(1)
      do ke=lmesh%NeS, lmesh%NeE
      do kk=1, elem%Np
        tmp(kk,ke,1) = ddens_(kk,ke)
        tmp(kk,ke,2) = momx_(kk,ke)
        tmp(kk,ke,3) = momy_(kk,ke)
        tmp(kk,ke,4) = momz_(kk,ke)
        tmp(kk,ke,5) = drhot_(kk,ke)
      end do
      end do      
    end if
 
    call elementoperationgpu_modalfilter_var5( elem_operation_gpu_driver, tmp, work, lmesh%Ne, &
      tmp_out ) ! (out)
 
    if ( do_weight_gsqrt_ ) then
      !$acc parallel loop collapse(2) present( DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, lmesh, elem ) async(1)
      do ke=lmesh%NeS, lmesh%NeE
      do kk=1, elem%Np
        rgsqrt = 1.0_rp / lmesh%Gsqrt(kk,ke)
        ddens_(kk,ke) = tmp_out(kk,ke,1) * rgsqrt
        momx_(kk,ke) = tmp_out(kk,ke,2) * rgsqrt
        momy_(kk,ke) = tmp_out(kk,ke,3) * rgsqrt
        momz_(kk,ke) = tmp_out(kk,ke,4) * rgsqrt
        drhot_(kk,ke) = tmp_out(kk,ke,5) * rgsqrt
      end do
      end do    
    else
      !$acc parallel loop collapse(2) present( DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, lmesh, elem ) async(1)
      do ke=lmesh%NeS, lmesh%NeE
      do kk=1, elem%Np
        ddens_(kk,ke) = tmp_out(kk,ke,1)
        momx_(kk,ke) = tmp_out(kk,ke,2)
        momy_(kk,ke) = tmp_out(kk,ke,3)
        momz_(kk,ke) = tmp_out(kk,ke,4)
        drhot_(kk,ke) = tmp_out(kk,ke,5)
      end do
      end do    
    end if
    !$acc wait(1)
    !$acc end data
 
    call elem_operation_gpu_driver%Final()
    return
  end subroutine atm_dyn_dgm_modalfilter_apply
#endif
 
 
!> Apply a modal filtering to tracer variables 
!OCL SERIAL
  subroutine atm_dyn_dgm_tracer_modalfilter_apply(  &
    QTRC_,                                   & ! (inout)
    dens_hyd_, ddens0_, ddens_,              & ! (in)
    lmesh, elem, elem_operation              ) ! (in)
 
    implicit none
 
    class(localmeshbase), intent(in) :: lmesh
    class(elementbase), intent(in) :: elem    
    real(rp), intent(inout)  :: qtrc_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: dens_hyd_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: ddens0_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: ddens_ (elem%np,lmesh%nea)
    class(elementoperationbase3d), intent(in) :: elem_operation
 
    integer :: ke
    real(rp) :: tmp(elem%np)
    real(rp) :: work(elem%np)
    real(rp) :: tmp_out(elem%np)
    real(rp) :: weight(elem%np)
 
    real(rp) :: qtrc_avg
    real(rp) :: gsqrtrhoq_ref(elem%np)
    integer :: kk
    !------------------------------------
 
    !$omp parallel do private( tmp, tmp_out, work, kk, weight, qtrc_avg, GsqrtRHOQ_ref )
    do ke=lmesh%NeS, lmesh%NeE
 
      !    ( dens_hyd + ddens0 ) * ( qtrc_av + dqtrc )
      ! = ( dens_hyd + ddens ) * qtrc_avg + dens_ * dqtrc
      ! 
      !  dens_hyd * qtrc_avg + qtrc_avg [ MF * ddens ] 
      ! + MF_q [ dqtrc ]
      qtrc_avg = 0.125_rp * sum(elem%IntWeight_lgl(:) * qtrc_(:,ke))
      gsqrtrhoq_ref(:) = lmesh%Gsqrt(:,ke) * qtrc_avg * ( &
        dens_hyd_(:,ke) + ddens_(:,ke) - ddens0_(:,ke) )
 
      do kk=1, elem%Np
        weight(kk) = lmesh%Gsqrt(kk,ke) &
                   * ( dens_hyd_(kk,ke) + ddens0_(kk,ke) ) 
!        tmp(kk) = weight(kk) * QTRC_(kk,ke)
        tmp(kk) = weight(kk) * qtrc_(kk,ke) &
                - gsqrtrhoq_ref(kk)
      end do
 
      ! del(RHO Q) = RHOQ - RHO_hyd Q_avg
      !            = (RHO_hyd + dRHO) (Q_avg + dQ) - RHO_hyd Q_avg
      !            = RHO_hyd dQ + dRHO Q_avg + dR
      ! RHOQ = RHO_hyd Q_avg + del(RHO Q)
      call elem_operation%ModalFilter_tracer( tmp, work, &
        tmp_out )
 
      qtrc_(:,ke) = ( gsqrtrhoq_ref(:) + tmp_out(:) ) / weight(:)
    end do    
 
    return
  end subroutine atm_dyn_dgm_tracer_modalfilter_apply
 
! !OCL SERIAL
!   subroutine atm_dyn_dgm_tracer_modalfilter_apply(  &
!     QTRC_,                                   & ! (inout)
!     DENS_hyd_, DDENS_,                       & ! (in)
!     lmesh, elem, elem_operation              ) ! (in)
 
!     implicit none
 
!     class(LocalMeshBase), intent(in) :: lmesh
!     class(ElementBase), intent(in) :: elem    
!     real(RP), intent(inout)  :: QTRC_(elem%Np,lmesh%NeA)
!     real(RP), intent(in)  :: DENS_hyd_(elem%Np,lmesh%NeA)
!     real(RP), intent(in)  :: DDENS_(elem%Np,lmesh%NeA)
!     class(ElementOperationBase3D), intent(in) :: elem_operation
 
!     integer :: ke
!     real(RP) :: tmp(elem%Np)
!     real(RP) :: work(elem%Np)
!     real(RP) :: tmp_out(elem%Np)
!     real(RP) :: weight(elem%Np)
!     integer :: kk
!     !------------------------------------
 
!     !$omp parallel do private( tmp, tmp_out, work, kk, weight )
!     do ke=lmesh%NeS, lmesh%NeE
 
!       do kk=1, elem%Np
!         weight(kk) = lmesh%Gsqrt(kk,ke) &
!                    * ( DENS_hyd_(kk,ke) + DDENS_(kk,ke) ) 
!         tmp(kk) = weight(kk) * QTRC_(kk,ke)
!       end do
 
!       call elem_operation%ModalFilter_tracer( tmp, work, &
!         tmp_out )
 
!       QTRC_(:,ke) = tmp_out(:) / weight(:)
!     end do    
 
!     return
!   end subroutine atm_dyn_dgm_tracer_modalfilter_apply
 
end module scale_atm_dyn_dgm_modalfilter