src/core/grid/IcGrid2D_FVM_Manager.f90

Go to the documentation of this file.

module IcGrid2D_FVM_Manager

  ! モジュール引用 ; Use statements 
  !

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP ! 倍精度実数型 ; Double recision
    
  ! 数学定数
  ! Mathematic
  !
  use igmcore_math, only: PI     ! 円周率. PI.

  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: vec_length, &
    &                       rotateZ

  ! 座標変換
  ! Coordinate conversion
  !
  use igmcore_coordinate_conversion, only: geo_to_orth_pos

  ! 宣言文 ; Declaration statemenets
  !

  implicit none
  private
  
  ! 公開変数
  ! Public variables 
  !

  integer, parameter, public :: RC_REGIONS_NUM = 10     

  integer, parameter, public :: ICOSAHEDRON_VERTEX_NUM = 12     

  integer, parameter, public :: CELL_POINTS_NUM = 6     

  integer, parameter, public :: NOT_POLE_FLAG  =  0

  integer, parameter, public :: NORTH_POLE_FLAG = 1

  integer, parameter, public :: SOUTH_POLE_FLAG = 2

  ! 格子タイプを表す定数
  !

  integer, parameter, public :: GTYPE_STDGRID = 1

  integer, parameter, public :: GTYPE_STDGCGRID = 2

  integer, parameter, public :: GTYPE_SPRGCGRID = 3


  !
  type, public :: IcGrid2D_FVM

    ! 公開メンバ変数 
    ! Public memeber variables
    !

    real(DP), pointer :: rcs_AGrid(:,:,:,:) => null()
    ! f2003: real(DP), public, allocatable :: rcs_AGrid(:,:,:,:)

    real(DP), pointer :: rcs_CV(:,:,:,:,:) => null()
    ! f2003: real(DP), public, allocatable :: rcs_CV(:,:,:,:,:)

    real(DP), pointer :: rcs_CVArea(:,:,:) => null()
    ! f2003: real(DP), public, allocatable :: rcs_CVArea(:,:,:)

    ! 非公開メンバ変数
    ! Private member variables
    !

    integer :: h_glevel
    ! f2003: integer, private :: h_glevel

    real(DP) :: icRadius
    ! f2003: real(DP), private :: icRadius


    integer :: grid_EffSize
    ! f2003: integer, private :: grid_EffSize

    integer :: grid_margine = 1
    ! f2003: integer, private :: grid_margine = 1

    integer :: idMin
    ! f2003: integer, private :: idMin

    integer :: idMax
    ! f2003: integer, private :: idMax

    integer :: idEfMin = 1
    ! f2003: integer, private :: idEfMin = 1

    integer :: idEfMax
    ! f2003: integer, private :: idEfMax
    
  end type IcGrid2D_FVM
  

  ! 公開関数
  ! Public procedures
  !
  public :: IcGrid2D_FVM_Init, IcGrid2D_FVM_Final

  public :: malloc_GP_array, malloc_CV_array
  public :: generate_CV5_GPindex, generate_CV6_GPindex, check_pole

  public :: set_rcregion_bounds, paste_margin_width
  public :: get_EffSize_Min, get_EffSize_Max, get_IdMin, get_IdMax
  public :: get_glevel, get_IcRadius

  public :: calc_icosahedron_vertex

contains

!
subroutine IcGrid2D_FVM_Init( &
  & self,              &  ! (inout)
  & glevel, ic_radius  &  ! (in)
  )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self
  integer, intent(in) :: glevel
  real(DP), intent(in) :: ic_radius
  
  ! 実行文 ; Executable statements
  !
  self%h_glevel = glevel
  self%icRadius = ic_radius

  ! 指定された水平格子の分割レベルに対応する要素数を計算する. 
  self%grid_EffSize = 2**self%h_glevel + 1 
  
  ! IcGrid2D_FVM クラスの配列変数のインデックスののりしろを含めた下限・上限,
  ! 物理的に意味のある下限・上限を求める. 
  self%idEfMax = self%idEfMin + self%grid_EffSize - 1
  self%idMin = self%idEfMin - self%grid_margine
  self%idMax = self%idEfMax + self%grid_margine
  
end subroutine IcGrid2D_FVM_Init

!
subroutine malloc_GP_array( &
  & self                    &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self  ! IcGrid2D_FVM クラスのインスタンスの参照.
  
  ! 実行文 ; Executable statements
  !
  
  !write(*,*) 'allocate grid point array'
  allocate(self%rcs_AGrid(RC_REGIONS_NUM, self%idMin:self%idMax, self%idMin:self%idMax, 3))  

end subroutine malloc_GP_array

!
subroutine malloc_CV_array( &
  & self                    &  ! (inout)
  & )
  
  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self
  
  ! 実行文 ; Executable statements
  !

  !write(*,*) 'allocate controll volume array'
  allocate( self%rcs_CV(RC_REGIONS_NUM, self%idEfMin:self%idEfMax, self%idEfMin:self%idEfMax, CELL_POINTS_NUM, 3) )
  allocate( self%rcs_CVArea(RC_REGIONS_NUM, self%idEfMin:self%idEfMax, self%idEfMin:self%idEfMax) )
  
end subroutine malloc_CV_array
 
!
subroutine IcGrid2D_FVM_Final( &
  & self             &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self

  ! 実行文 ; Executable statements
  !

  ! 本モジュールにおいて, 動的に確保した配列を解放する.
  deallocate(self%rcs_AGrid)
  deallocate(self%rcs_CV)
  deallocate(self%rcs_CVArea)

end subroutine IcGrid2D_FVM_Final

!
function get_EffSize_Min( &
  & self                  &  ! (in)
  & ) result(gridEffMin)
 
  type(IcGrid2D_FVM), intent(in) :: self
  integer gridEffMin
  
  gridEffMin = self%idEfMin

end function get_EffSize_Min

!
function get_EffSize_Max( &
  & self                  &  ! (in)
  & ) result(gridEffMax)

  type(IcGrid2D_FVM), intent(in) :: self
  integer gridEffMax
  
  gridEffMax = self%idEfMax

end function get_EffSize_Max

!
function get_IdMax( &
  & self            &  ! (in)
  & ) result(id_Max)

  type(IcGrid2D_FVM), intent(in) :: self
  integer id_Max
  
  id_Max = self%idMax

end function get_IdMax

!
function get_IdMin( &
  & self            &  ! (in)
  & ) result(id_Min)

  type(IcGrid2D_FVM), intent(in) :: self
  integer id_Min
  
  id_Min = self%idMin

end function get_IdMin

!
function get_IcRadius( &
  & self               &  ! (in)
  & ) result(ic_radius)
 
  type(IcGrid2D_FVM), intent(in) :: self
  real(DP) ic_radius
  
  ic_radius = self%icRadius
end function get_IcRadius

!
function get_glevel( &
  & self             &  ! (in)
  & ) result(glevel)
 
  type(IcGrid2D_FVM), intent(in) :: self  ! IcGrid2D_FVM クラスのインスタンスの参照
  integer glevel
  
  glevel = self%h_glevel
end function get_glevel

!
subroutine set_rcregion_bounds(    &
  & self,                          &  ! (inout)
  & rcID, top, left, bottom, right &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self
  integer, intent(in) :: rcID
  real(DP), intent(in) :: top(3)
  real(DP), intent(in) :: left(3)
  real(DP), intent(in) :: bottom(3)
  real(DP), intent(in) :: right(3)
  
  self%rcs_AGrid(rcID, self%idEfMin, self%idEfMin, :) = top
  self%rcs_AGrid(rcID, self%idEfMax, self%idEfMin, :) = left
  self%rcs_AGrid(rcID, self%idEfMax, self%idEfMax, :) = bottom
  self%rcs_AGrid(rcID, self%idEfMin, self%idEfMax, :) = right

end subroutine set_rcregion_bounds

!
subroutine paste_margin_width( &
  & self                       &  ! (inout) 
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self

  ! 作業変数
  ! Work variables
  !
  integer :: idlist1(7) = (/ 5, 1, 2, 3, 4, 5, 1 /)
  integer :: idlist2(7) = (/ 10, 6, 7, 8, 9, 10, 6 /)
  integer :: EMin, EMax, marginw
  integer i, k

  ! 実行文 ; Executable statements
  !

  EMin = get_EffSize_Min(self)
  EMax = get_EffSize_Max(self)
  marginw = EMin - get_IdMin(self)

  do i=1, RC_REGIONS_NUM/2
    do k=1, marginw
      ! 北半球側の矩形領域ののりしろに値をコピーする.
      self%rcs_AGrid(i, EMin:EMax, EMin-k, :) = self%rcs_AGrid(idlist1(i), EMin+k, EMin:EMax, :)
      self%rcs_AGrid(i, EMax+k, EMin:EMax, :) = self%rcs_AGrid(idlist2(i), EMin+k, EMin:EMax, :)
      self%rcs_AGrid(i, EMin:EMax, EMax+k, :) = self%rcs_AGrid(idlist2(i+1), EMin:EMax, EMin+k, :)
      self%rcs_AGrid(i, EMin-k, EMin:EMax, :) = self%rcs_AGrid(idlist1(i+2), EMin:EMax, EMin+k, :)

      ! 南半球側の矩形領域ののりしろに値をコピーする.
      self%rcs_AGrid(i+5, EMin:EMax, EMin-k, :) = self%rcs_AGrid(idlist1(i+1), EMin:EMax, EMax-k, :)
      self%rcs_AGrid(i+5, EMax+k, EMin:EMax, :) = self%rcs_AGrid(idlist2(i), EMin:EMax, EMax-k, :)
      self%rcs_AGrid(i+5, EMin:EMax, EMax+k, :) = self%rcs_AGrid(idlist2(i+2), EMax-k, EMin:EMax, :)
      self%rcs_AGrid(i+5, EMin-k, EMin:EMax, :) = self%rcs_AGrid(idlist1(i+2), EMax-k, EMin:EMax, :)
    end do
  end do

end subroutine paste_margin_width

!
function generate_CV5_GPindex( &
  & self,                      &  ! (inout)
  & GP_i, GP_j, rcID           &  ! (in)
  & ) result(CV_GPId)
  
  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self
  integer, intent(in) :: GP_i
  integer, intent(in) :: GP_j
  integer, intent(in) :: rcID
  integer CV_GPId(5,3,2)  

  ! 作業変数
  ! Work variables
  !
  integer :: EMin,  EMax
  integer i,j 
  integer :: CV_GPIdShape(2) = (/ 3, 2 /)

  ! 実行文 ; Executable statements
  !
  
  EMin = get_EffSize_Min(self)
  EMax = get_EffSize_Max(self)

  i = GP_i; j = GP_j
  if( rcID <= 5 )then
    ! 北半球に対して
    if( i==EMax .and. j == EMax ) then
      ! 下端特異点
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j+1 /), CV_GPIdShape ) 
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j+1,j /), CV_GPIdShape ) 
      CV_GPId(3,:,:) = reshape( (/ i,i-1,i, j,j,j-1 /), CV_GPIdShape ) 
      CV_GPId(4,:,:) = reshape( (/ i,i,i+1, j,j-1,j-1 /), CV_GPIdShape ) 
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i+1, j,j-1,j /), CV_GPIdShape ) 
    else if( j == EMin ) then
      ! 左端特異点
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j+1 /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j+1,j /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i-1,i-1, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i-1,i, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i, j,j,j+1 /), CV_GPIdShape )
      return
    else if( j == EMax ) then
      ! 右端特異点
      CV_GPId(1,:,:) = reshape( (/ i,i-1,i-1, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i,i+1, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i+1,i+1, j,j-1,j /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i, j,j,j+1 /), CV_GPIdShape )
      return
    end if
  else
    ! 南半球に対して
    if ( i==EMin .and. j==EMin ) then
      ! 上端特異点
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j+1 /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j+1,j /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i,i+1, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i+1,i+1, j,j-1,j /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i, j,j,j+1 /), CV_GPIdShape ) 
      return
    else if( j == EMin ) then
      ! 左端特異点
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j+1 /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j+1,j /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i-1,i, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i+1,i+1, j,j,j+1 /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i, j,j+1,j+1 /), CV_GPIdShape )
      return
    else if( j == EMax ) then
      ! 右端特異点
      CV_GPId(1,:,:) = reshape( (/ i,i-1,i, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i,i+1, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i+1,i+1, j,j-1,j /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i+1,i+1, j,j,j+1 /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i, j,j+1,j+1 /), CV_GPIdShape )
      return
    end if
  end if

end function generate_CV5_GPindex

!
function generate_CV6_GPindex( &
  & self,                      &  ! (inout)
  & GP_i, GP_j, rcID           &  ! (in)
  & ) result(CV_GPId)
 
  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self
  integer, intent(in) :: GP_i
  integer, intent(in) :: GP_j
  integer, intent(in) :: rcID
  integer :: CV_GPId(6,3,2)
  
  ! 作業変数
  ! Work variables
  !
  integer :: i,j, EMin, EMax 
  integer :: CV_GPIdShape(2) = (/ 3, 2 /)

  ! 実行文 ; Executable statements
  ! 

  EMin = get_EffSize_Min(self)
  EMax = get_EffSize_Max(self)
  
  i = GP_i; j = GP_j

  ! 先に, インデックス i,j が矩形領域の境界に当たる場合を処理する. 
  if( rcID <= 5 )then
    ! 北半球の矩形領域に対して
    if(  i == 1 ) then
      ! 右上端
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i-1,i, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i,i+1, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i+1, j,j-1,j /), CV_GPIdShape )
      CV_GPId(6,:,:) = reshape( (/ i,i+1,i, j,j,j+1  /), CV_GPIdShape ) 
      return
    else if( j == 1 ) then
      ! 左上端
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j+1 /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j+1,j /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i-1,i-1, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i-1,i, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i,i+1, j,j-1,j /), CV_GPIdShape )
      CV_GPId(6,:,:) = reshape( (/ i,i+1,i, j,j,j+1 /), CV_GPIdShape )           
      return
    end if

  else
    ! 南半球の矩形領域に対して
    if( i == EMax ) then
      ! 左下端
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j+1 /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j+1,j /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i-1,i, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i,i+1, j,j-1,j /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i+1, j,j,j+1 /), CV_GPIdShape )
      CV_GPId(6,:,:) = reshape( (/ i,i+1,i, j,j+1,j+1 /), CV_GPIdShape ) 
      return
    else if( j == EMax ) then
      ! 右下端
      CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j /), CV_GPIdShape )
      CV_GPId(2,:,:) = reshape( (/ i,i-1,i, j,j,j-1 /), CV_GPIdShape )
      CV_GPId(3,:,:) = reshape( (/ i,i,i+1, j,j-1,j-1 /), CV_GPIdShape )
      CV_GPId(4,:,:) = reshape( (/ i,i+1,i+1, j,j-1,j /), CV_GPIdShape )
      CV_GPId(5,:,:) = reshape( (/ i,i+1,i+1, j,j,j+1 /), CV_GPIdShape )
      CV_GPId(6,:,:) = reshape( (/ i,i,i+1, j,j+1,j+1 /), CV_GPIdShape )          
      return
    end if
  end if
  
  ! インデックス i,j が矩形領域の内部(境界以外)に対応する場合
  CV_GPId(1,:,:) = reshape( (/ i,i,i-1, j,j+1,j+1 /), CV_GPIdShape )
  CV_GPId(2,:,:) = reshape( (/ i,i-1,i-1, j,j+1,j /), CV_GPIdShape )
  CV_GPId(3,:,:) = reshape( (/ i,i-1,i, j,j,j-1 /), CV_GPIdShape )
  CV_GPId(4,:,:) = reshape( (/ i,i,i+1, j,j-1,j-1 /), CV_GPIdShape )
  CV_GPId(5,:,:) = reshape( (/ i,i+1,i+1, j,j-1,j /), CV_GPIdShape )
  CV_GPId(6,:,:) = reshape( (/ i,i+1,i, j,j,j+1 /), CV_GPIdShape ) 

end function generate_CV6_GPindex

!
function check_pole(self, rcID, i, j) result (pole_flag)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: self
  integer, intent(in) :: rcID
  integer, intent(in) :: i
  integer, intent(in) :: j
  integer pole_flag

  if( rcID <= 5 .and. i==self%idEfMin .and. j== self%idEfMin ) then
    pole_flag = NORTH_POLE_FLAG
  else if( rcID > 5 .and. i==self%idEfMax .and. j== self%idEfMax ) then
    pole_flag = SOUTH_POLE_FLAG
  else
    pole_flag = NOT_POLE_FLAG
  end if

end function check_pole

!
subroutine calc_icosahedron_vertex(orth_icvertex)
  
  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(inout) :: orth_icvertex(ICOSAHEDRON_VERTEX_NUM, 3)   

  ! 作業変数
  ! Work variable
  !
  real(DP) :: unit_radius = 1.0d0
  real(DP) :: del_lambda = 2.0d0 * PI / 5.0d0
  real(DP) :: pos2_theta
  real(DP) :: pos2_lambda 
  integer i

  ! 実行文 ; Executable statements
  ! 
  
  !　頂点1, 12 を z 軸上に配置する. 
  orth_icvertex(1, :) = (/ 0.0d0, 0.0d0, unit_radius /)
  orth_icvertex(12, :) = (/ 0.0d0, 0.0d0, - unit_radius /)
   
  ! 頂点 2 の緯度・経度を設定する. 
  pos2_lambda = 0.0d0 
  pos2_theta = 2.0d0 * dasin( 1.0d0 / (2.0d0 * dsin(0.2d0 * PI)) ) - PI/2.0d0
  
  ! 地理座標系における頂点 2, 7 の座標をデカルト座標系における座標へ変換する. 
  ! なお, 頂点 7 は頂点 2 と xy 平面を挟んで反対側に存在する. 
  orth_icvertex(2, :) = geo_to_orth_pos( (/ pos2_lambda, pos2_theta, unit_radius /) )
  orth_icvertex(7, :) = geo_to_orth_pos( (/ pos2_lambda + del_lambda/2.0d0, - pos2_theta, unit_radius /) )

  ! z 軸を軸にして, 座標が既知である頂点を del_lambda づつ回転させて, 残りの頂点座標を計算する.  
  do i=3, 6
    orth_icvertex(i, :) = rotateZ( orth_icvertex(i-1, :), del_lambda )
    orth_icvertex(i+5, :) = rotateZ( orth_icvertex( (i+5)-1, : ), del_lambda )
  end do

end subroutine calc_icosahedron_vertex

end module IcGrid2D_FVM_Manager