src/util/math/Derivate_Field_IcGrid2D_Manager.f90

Go to the documentation of this file.

module Derivate_Field_IcGrid2D_Manager

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

  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision

  ! 球面三角法
  ! Spherical trigonometry
  !
  use igmcore_spherical_trigonometry, only: &
    & spherical_tri_area


  ! 一つの格子点に対する微分演算子の提供
  ! Providing with the differential operation for a single grid-point
  !
  use FVM_HDiff_operator, only: &
    & eval_div_operator, eval_rot_operator, eval_grad_operator

  ! 正二十面体格子データの管理
  ! Manager of icosahedral grid data
  !
   use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & paste_margin_width, generate_CV5_GPindex, generate_CV6_GPindex, &
    & get_IdMin, get_EffSize_Min, get_EffSize_Max, get_IcRadius, &
    & RC_REGIONS_NUM

  ! 正二十面体格子上の物理場データの管理
  ! Manager of the physical field data on icosahedral grid
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D

  ! OpenMP
  !
  !
  use omp_lib

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開変数
  ! Public variables
  !

  !
  type, public :: Derivate_Field_IcGrid2D

    type(IcGrid2D_FVM), pointer :: icgrid
    ! f2003: type(IcGrid2D_FVM), pointer, private :: icgrid


    real(DP), pointer :: rcs_areaGPP(:,:,:,:,:)
    ! f2003: real(DP), allocatable, private :: rcs_areaGPP(:,:,:,:,:)

    integer, pointer :: rcU_CVGPindex(:,:,:,:,:)
    ! f2003: integer, allocatable, private :: rcU_CVGPindex(:,:,:,:,:)

    integer, pointer :: rcL_CVGPindex(:,:,:,:,:)
    ! integer, allocatable, private :: rcL_CVGPindex(:,:,:,:,:)

  end type Derivate_Field_IcGrid2D

  ! 公開メンバ関数
  ! Public memeber function
  !
  public :: Derivate_Field_IcGrid2D_Init, Derivate_Field_IcGrid2D_Final
  public :: divergence_op, vertical_curl_op, gradient_op

contains

!
subroutine Derivate_Field_IcGrid2D_Init( &
  & self,              &  ! (inout)
  & icgrid             &  ! (in)
  & )

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

  ! 作業変数
  ! Work variables
  !
  integer :: EMin, EMax

  ! 実行文 ; Declaration statements
  !

  self%icgrid => icgrid
  EMin = get_EffSize_Min(self%icgrid)
  EMax = get_EffSize_Max(self%icgrid)

  ! Derivate_Field_IcGrid2D クラスが作業変数として保持する配列を確保する.
  allocate(self%rcs_areaGPP(RC_REGIONS_NUM, EMin:EMax, EMin:EMax, 6, 3))
  allocate(self%rcU_CVGPindex(EMin:EMax, EMin:EMax, 6,3,2 ))
  allocate(self%rcL_CVGPindex(EMin:EMax, EMin:EMax, 6,3,2 ))

  ! コントールボリュームの頂点座標を計算するために必要な格子点のインデックスを, あらかじめ計算しておく.
  call calc_CVGPindex(self)
  ! コントールボリュームの頂点上での物理場の値を計算する際に必要となる, 面積データをあらかじめ計算しておく.
  call calc_GPP_area_weight(self)

end subroutine Derivate_Field_IcGrid2D_Init

!
subroutine Derivate_Field_IcGrid2D_Final( &
  & self &  ! (inout)
  & )

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

  ! 実行文 ; Declaration statements
  !

  deallocate(self%rcs_areaGPP)
  deallocate(self%rcU_CVGPindex)
  deallocate(self%rcL_CVGPindex)

end subroutine Derivate_Field_IcGrid2D_Final

!
subroutine divergence_op( &
  & self,                &  ! (inout)
  & vector_field,        &  ! (in)
  & ret_scalar_field     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: vector_field
  type(Field_IcGrid2D), intent(inout) :: ret_scalar_field

  ! 作業変数
  ! Work variables
  !
  integer :: idMin

  ! 実行文 ; Executable statements
  !

  idMin = get_IdMin(self%icgrid)

  call eval_diffOptr( &
    & self=self, &
    & diffOptr_fnc=eval_div_operator, rcs_GP_val=vector_field%val, val_dim=3, &  ! (in)
    & ret_rcs_val=ret_scalar_field%val,                        &       ! (inout)
    & ret_dim=1, idMin=idMin )                                         ! (in)

end subroutine divergence_op

!
subroutine vertical_curl_op( &
  & self,                &  ! (inout)
  & vector_field,        &  ! (in)
  & ret_scalar_field     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: vector_field
  type(Field_IcGrid2D), intent(inout) :: ret_scalar_field

  ! 作業変数
  ! Work variables
  !
  integer :: idMin

  ! 実行文 ; Executable statements
  !

  idMin = get_IdMin(self%icgrid)

  call eval_diffOptr( &
    & self=self, &
    & diffOptr_fnc=eval_rot_operator, rcs_GP_val=vector_field%val, val_dim=3, &  ! (in)
    & ret_rcs_val=ret_scalar_field%val,                        &  ! (inout)
    & ret_dim=1, idMin=idMin )                                    ! (in)

end subroutine vertical_curl_op

!
subroutine gradient_op( &
  & self,                &  ! (inout)
  & scalar_field,        &  ! (in)
  & ret_vector_field     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  type(Field_IcGrid2D), intent(in) :: scalar_field
  type(Field_IcGrid2D), intent(inout) :: ret_vector_field

  ! 作業変数
  ! Work variables
  !
  integer :: idMin

  ! 実行文 ; Executable statements
  !

  idMin = get_IdMin(self%icgrid)

  call eval_diffOptr( &
    & self=self, &
    & diffOptr_fnc=eval_grad_operator, rcs_GP_val=scalar_field%val, val_dim=1, &  ! (in)
    & ret_rcs_val=ret_vector_field%val,                        &   ! (inout)
    & ret_dim=3, idMin=idMin )                                     ! (in)

end subroutine gradient_op

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
!
! モジュール外非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
subroutine eval_diffOptr( &
  & self,                           &  ! (inout)
  & diffOptr_fnc, rcs_GP_val, val_dim,  &  ! (in)
  & ret_rcs_val,                    &  ! (inout)
  & ret_dim, idMin                  &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(Derivate_Field_IcGrid2D), intent(inout), target :: self
  integer, intent(in) :: idMin
  integer, intent(in) :: val_dim
  integer, intent(in) :: ret_dim

  interface
    !!
    !import  DP

    !
    function diffOptr_fnc( &
      & val_GP0, val_CV, GP0, CV,                   &  ! (in)
      & CVArea, ic_radius, CV_num, val_dim, ret_dim &  ! (in)
      & ) result(val)

      use dc_types, only: DP
      implicit none

      ! 宣言文 ; Declaration statements
      !

      integer, intent(in) :: val_dim
      integer, intent(in) :: ret_dim
      integer, intent(in) :: CV_num
      real(DP), intent(in) :: val_GP0(val_dim)
      real(DP), intent(in) :: GP0(3)
      real(DP), intent(in) :: val_CV(CV_num, val_dim)
      real(DP), intent(in) :: CV(CV_num, val_dim)
      real(DP), intent(in) :: CVArea
      real(DP), intent(in) :: ic_radius
      real(DP) val(ret_dim)

    end function diffOptr_fnc
  end interface

  real(DP), intent(in) :: rcs_GP_val(:,idMin:,idMin:,:)
  real(DP), intent(inout) :: ret_rcs_val(:,idMin:,idMin:,:)

  ! 作業変数
  ! Work variables
  !
  real(DP) :: ic_radius
  real(DP), allocatable :: CV_val(:,:,:,:,:)
  integer :: rcID, i,j
  integer :: EMin, EMax
  integer :: CV_num

  ! 実行文 ; Executable statements
  !

  EMin = get_EffSize_Min(self%icgrid)
  EMax = get_EffSize_Max(self%icgrid)
  ic_radius = get_IcRadius(self%icgrid)

  ! コントロールボリュームの頂点上の物理場の値を保持するための配列のメモリを確保する.
  allocate(CV_Val(RC_REGIONS_NUM, EMin:EMax, EMin:EMax, 6, 3))

  !
  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do private(CV_num)
    do j=EMin, EMax
      do i=EMin,EMax

        if ( ( i==EMin .and. j==EMin ) .or. ( i==EMax .and. j==EMin ) .or. &
             ( i==EMin .and. j==EMax ) .or. ( i==EMax .and. j==EMax ) ) then
          ! 特異格子点のコントールボリュームの頂点数は 5 個
          CV_num = 5
        else
          ! 通常ののコントールボリュームの頂点数は 6 個
          CV_num = 6
        end if

        ! 格子点(rcID,i,j) のコントールボリュームの頂点上における物理場の値を計算する.
        call calc_CV_val( &
          & self, &
          & CV_val(rcID, i, j, 1:CV_num, 1:val_dim), &      ! (inout)
          & rcs_GP_val(:,:,:,1:val_dim), &                  ! (in)
          & rcID, i,j, val_dim, CV_num, idMin, EMin, EMax ) ! (in)

        ! 動的に接続された微分演算オブジェクトに対して微分演算を実行し, 格子点(rcID,i,j)上の微分値を計算する.
        ! また, その結果を保存する.
        ret_rcs_val(rcID, i,j, 1:ret_dim) = diffOptr_fnc( &
          & val_GP0=rcs_GP_val(rcID,i,j,1:val_dim), val_CV=CV_val(rcID,i,j, 1:CV_num,1:val_dim), & ! (in)
          & GP0=self%icgrid%rcs_AGrid(rcID,i,j,:), CV=self%icgrid%rcs_CV(rcID,i,j,1:CV_num,:), &   ! (in)
          & CVArea=self%icgrid%rcs_CVArea(rcID,i,j), &                                             ! (in)
          & ic_radius=ic_radius, CV_num=CV_num, val_dim=val_dim, ret_dim=ret_dim )                 ! (in)

      end do
    end do
  end do

end subroutine eval_diffOptr

!
subroutine calc_CV_val( self, CV_val, rcs_GP_val, rcID, gp_i, gp_j, val_dim, CV_num, idMin, EMin, EMax)

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self
  integer, intent(in) :: val_dim, CV_num, idMin, EMin, EMax
  real(DP), intent(inout) :: CV_val(CV_num, val_dim)
  real(DP), intent(in) :: rcs_GP_val(:,idMin:,idMin:,:)
  integer, intent(in) :: rcID, gp_i, gp_j

  ! 作業変数
  ! Work variable
  !
  integer :: cvID, k
  real(DP) :: GP_val3(3,val_dim)
  real(DP) :: area3(3)

  ! 実行文 ; Executable statement
  !

  do cvID=1, CV_num

    ! 矩形領域 rcID 内の格子点 (gp_i, gp_j) に付随するコントロールボリュームの一頂点 cvID と近傍の格子点が作る 3 つの三角形の
    ! 面積を取得する.
    area3(:) = self%rcs_areaGPP(rcID,gp_i,gp_j,cvID, :)

    !
    if( rcID <= 5 ) then
      if ( gp_i==EMin .and. gp_j==EMin )then
        GP_val3(1,:) = rcs_GP_val(cvID, EMin, EMin, :); GP_val3(2,:) = rcs_GP_val(cvID, EMin+1, EMin, :);
        GP_val3(3,:) = rcs_GP_val(cvID, EMin, EMin+1, :);
      else
        do k=1,3
          GP_val3(k,:) = rcs_GP_val(rcID, self%rcU_CVGPindex(gp_i,gp_j,cvID,k,1), self%rcU_CVGPindex(gp_i,gp_j,cvID,k,2), :)
        end do
      end if

    else
      if ( gp_i==EMax .and. gp_j==EMax )then
        GP_val3(1,:) = rcs_GP_val(11-cvID, EMax, EMax, :); GP_val3(2,:) = rcs_GP_val(11-cvID, EMax-1, EMax, :);
        GP_val3(3,:) = rcs_GP_val(11-cvID, EMax, EMax-1, :);
      else
        do k=1,3
          GP_val3(k,:) = rcs_GP_val(rcID, self%rcL_CVGPindex(gp_i,gp_j,cvID,k,1), self%rcL_CVGPindex(gp_i,gp_j,cvID,k,2), :)
        end do
      end if

    end if

    ! コントロールボリュームの各頂点(点 G)上の物理場の値を, 近傍の格子点(点 P0,P1,P2)上の物理場の値から補間する.
    ! 各格子点の物理場の値は, コントロールボリュームの頂点と格子点が作る三角形(GP0P1, GP1P2, GP2P0)の面積によって重みを付けられる.
    CV_val(cvID,1:val_dim) = ( &
      &  area3(1) *  GP_val3(1,:) + area3(2) * GP_val3(2,:) + area3(3) * GP_val3(3,:) &
      &  ) / ( area3(1) + area3(2) + area3(3) )

  end do

end subroutine calc_CV_val

!
subroutine calc_CVGPindex(self)

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self

  ! 作業変数
  ! Work variable
  !
  integer :: i, j
  integer :: EMin, EMax, idMin

  ! 実行文 ; Executable statement
  !

  EMin = get_EffSize_Min(self%icgrid)
  EMax = get_EffSize_Max(self%icgrid)
  idMin = get_IdMin(self%icgrid)

  do j=EMin,EMax
    do i=EMin,EMax
      if ( ( i==EMax .and. j==EMin ) .or. ( i==EMin .and. j==EMax )  ) then
        !
        !
        self%rcU_CVGPindex(i,j,1:5,:,:) = &
          & generate_CV5_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=1 )
        self%rcL_CVGPindex(i,j,1:5,:,:) = &
          & generate_CV5_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=6 )
      else

        !
        self%rcU_CVGPindex(i,j,1:6,:,:) = &
          & generate_CV6_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=1 )
        !
        self%rcL_CVGPindex(i,j,1:6,:,:) = &
          & generate_CV6_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=6 )

      end if
    end do
  end do

  self%rcU_CVGPindex(EMax,EMax,1:5,:,:) = &
    & generate_CV5_GPindex(self%icgrid, GP_i=EMax, GP_j=EMax, rcID=1 )
  self%rcL_CVGPindex(EMin,EMin,1:5,:,:) = &
    & generate_CV5_GPindex(self%icgrid, GP_i=EMin, GP_j=EMin, rcID=6 )

end subroutine calc_CVGPindex

!
subroutine calc_GPP_area_weight( &
  & self &   ! (inout)
  & )

  ! 宣言文 ; Declaration statement
  !
  type(Derivate_Field_IcGrid2D), intent(inout) :: self

  ! 作業変数
  ! Work variables
  !
  integer :: rcID, i, j, cvID
  integer :: EMin, EMax, idMin
  real(DP) :: ic_radius
  integer :: CV5_GPindx(5,3,2), CV6_GPindx(6,3,2)

  ! 実行文 ; Executable statement
  !

  EMin = get_EffSize_Min(self%icgrid)
  EMax = get_EffSize_Max(self%icgrid)
  idMin = get_IdMin(self%icgrid)
  ic_radius = get_IcRadius(self%icgrid)

  do rcID=1, RC_REGIONS_NUM
    do j=EMin,EMax
      do i=EMin,EMax
        if ( ( rcID > 5 .and. i==EMin .and. j==EMin ) .or. ( i==EMax .and. j==EMin ) .or. &
             ( i==EMin .and. j==EMax ) .or. ( rcID <= 5 .and. i==EMax .and. j==EMax ) ) then
          !
          CV5_GPindx(:,:,:) = generate_CV5_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=rcID )

          do cvID=1,5
            self%rcs_areaGPP(rcID,i,j,cvID, 1:3) = calc_usual_area3( CV_GPindx=CV5_GPindx(cvID, :,:), &
              & g=self%icgrid%rcs_CV(rcID,i,j,cvID,:), &
              & rc_AGrid=self%icgrid%rcs_AGrid(rcID,:,:,:), ic_radius=ic_radius, idMin=idMin )
          end do
        else if( rcID <= 5 .and. i==EMin .and. j==EMin ) then
          !
          do cvID=1,5
            self%rcs_areaGPP(rcID,i,j,cvID, 1:3) = calc_area3_from_4pts( ic_radius=ic_radius, &
              & p1=self%icgrid%rcs_AGrid(cvID,EMin,EMin,:), &
              & p2=self%icgrid%rcs_AGrid(cvID,EMin+1,EMin,:), &
              & p3=self%icgrid%rcs_AGrid(cvID,EMin,EMin+1,:), &
              & g=self%icgrid%rcs_CV(rcID,EMin,EMin,cvID,:) )
          end do

        else if( rcID > 5 .and. i==EMax .and. j==EMax ) then
          !
          do cvID=1,5
            self%rcs_areaGPP(rcID,i,j,cvID, 1:3) = calc_area3_from_4pts( ic_radius=ic_radius, &
              & p1=self%icgrid%rcs_AGrid(11-cvID,EMax,EMax,:), &
              & p2=self%icgrid%rcs_AGrid(11-cvID,EMax-1,EMax,:), &
              & p3=self%icgrid%rcs_AGrid(11-cvID,EMax,EMax-1,:), &
              & g=self%icgrid%rcs_CV(rcID,EMax,EMax,cvID,:) )
          end do

        else
          !
          CV6_GPindx(:,:,:) = generate_CV6_GPindex(self%icgrid, GP_i=i, GP_j=j, rcID=rcID )

          do cvID=1,6
            self%rcs_areaGPP(rcID,i,j,cvID, 1:3) = calc_usual_area3( CV_GPindx=CV6_GPindx(cvID, :,:), &
              &  g=self%icgrid%rcs_CV(rcID,i,j,cvID,:), &
              & rc_AGrid=self%icgrid%rcs_AGrid(rcID,:,:,:), ic_radius=ic_radius, idMin=idMin )
          end do

        end if
      end do
    end do
  end do

end subroutine calc_GPP_area_weight

!
function calc_usual_area3( &
  & CV_GPindx, g, rc_AGrid, ic_radius, idMin &  ! (in)
  & ) result(area3)

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: idMin
  integer, intent(in) :: CV_GPindx(3,2)
  real(DP), intent(in) :: g(3)
  real(DP), intent(in) :: rc_AGrid(idMin:,idMin:,:)
  real(DP), intent(in) :: ic_radius

  ! 作業変数
  ! Work variable
  !
  real(DP) area3(3)

  ! 実行文 ; Executable statement
  !

  area3(:) = calc_area3_from_4pts( &
    & p1=rc_AGrid(CV_GPindx(1,1), CV_GPindx(1,2), :), &
    & p2=rc_AGrid(CV_GPindx(2,1), CV_GPindx(2,2), :), &
    & p3=rc_AGrid(CV_GPindx(3,1), CV_GPindx(3,2), :), &
    & g=g(:), &
    & ic_radius=ic_radius )

end function calc_usual_area3

!
function calc_area3_from_4pts(p1, p2, p3, g, ic_radius) result(area3)

  ! 宣言文 ; Declaration statement
  !
  real(DP), intent(in) :: p1(3), p2(3), p3(3)
  real(DP), intent(in) :: g(3)
  real(DP), intent(in) :: ic_radius
  real(DP) area3(3)

  ! 実行文 ; Executable statement
  !

  area3(1) = spherical_tri_area(ic_radius, g, p2, p3)
  area3(2) = spherical_tri_area(ic_radius, g, p3, p1)
  area3(3) = spherical_tri_area(ic_radius, g, p1, p2)

end function calc_area3_from_4pts

end module Derivate_Field_IcGrid2D_Manager