src/core/grid/SPRGCGrid_Builder.f90

Go to the documentation of this file.

module SPRGCGrid_Builder

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

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

  ! 数学定数
  ! Mathematic
  !
  use igmcore_math, only: PI   ! 円周率

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

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

  ! 正二十面体格子データの管理
  ! Manager of Icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & paste_margin_width, generate_CV5_GPindex, &
    & get_idMin, get_idMax, get_EffSize_Min, get_EffSize_Max, &
    & get_IcRadius, get_glevel, &
    & check_pole, NOT_POLE_FLAG, &
    & RC_REGIONS_NUM


  ! STD-GC-grid の構築
  ! STD-GC-grid construction
  !
  use STDGCGrid_Builder, only: &
    & calc_STDGrid_GP_vertex, calc_CV_vertex, calc_CV_Area, &
    & move_GP_into_CVGravPt, project_on_sphere

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: construct_SPRGCGrid, set_iteration_limiter

  ! 非公開変数
  ! Private variables
  !

  real(DP) :: eq_d

  real(DP) :: radius

  real(DP) :: beta = 0.4d0

  ! ばねの振動の減衰のふるまいがおよそ臨界減衰となるように,
  ! 摩擦定数とバネ定数を選ぶ.

  real(DP) :: alpha = 2000.0d0

  real(DP) :: k = 1.0d06

  real(DP) :: dt = 0.001d0

  integer :: end_tstep = 20000

  logical :: tstep_limiter = .false.


  integer :: conv_check_step = 100

  integer, allocatable :: GP_IJ(:,:,:,:,:)

contains

!
subroutine construct_SPRGCGrid( &
  & icgrid,     &  ! (inout)
  & tstep_limit &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid
  integer, intent(in), optional :: tstep_limit

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

  ! 実行文 ; Executable statements
  !

  radius = get_IcRadius(icgrid)
  glevel = get_glevel(icgrid)

  !
  if ( present(tstep_limit) ) then
    tstep_limiter = .false.
    end_tstep = tstep_limit
  end if

  ! つり合いの状態における, 近傍の格子点までの距離の推定値を定める.
  ! ( Tomita,etal(2001)の式(23) )
  eq_d = beta * 2.0d0 * PI * radius / ( 10.0d0 * 2.0d0**( glevel-1 ) )

  ! SPR-GC-grid を生成する.
  !

  ! 最初に STD-grid の格子点座標を求める.
  call calc_STDGrid_GP_vertex(icgrid)
  ! バネ力学を用いて, STD-grid の格子の空間分布を準一様にする.
  call modify_STDGrid_with_spring_dyn(icgrid)

  ! 改良した格子でコントールボリュームの格子位置, 面積を求める.
  call calc_CV_vertex(icgrid)
  call calc_CV_Area(icgrid)

  ! 格子位置を, コントールボリュームの重心位置に移動させる.
  call move_GP_into_CVGravPt(icgrid)

end subroutine construct_SPRGCGrid

!
subroutine set_iteration_limiter( &
  & limiter_flag &  ! (in)
  & )

  ! 宣言文 ; Declaration statement
  !
  logical, intent(in)  :: limiter_flag

  ! 実行文 ; Executable statement
  !
  tstep_limiter = limiter_flag

end subroutine set_iteration_limiter

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

!
subroutine modify_STDGrid_with_spring_dyn( &
  & icgrid )                                    ! (inout)

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

  ! 作業変数
  ! Work variables
  !
  integer :: idMin,  idMax
  integer :: EMin,  EMax
  integer :: rcID, i, j
  integer :: pt_num

  ! 実行文 ; Executable statements
  !

  !write(*,*) 'modify STD-GC-grid'
  ! 配列の境界インデックスを取得
  idMin = get_idMin(icgrid)
  idMax = get_idMax(icgrid)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  ! バネ力学を利用して, STD-grid の格子位置を改良する.
  !

  allocate(GP_IJ(RC_REGIONS_NUM, EMin:EMax, EMin:EMax, 6, 2))

  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==EMin .and. j==EMax ) .or. &
            & ( i==EMax .and. j==EMin ) .or. ( rcID <= 5 .and. i==EMax .and. j==EMax )     &
            &       ) then
            pt_num = 5
            GP_IJ(rcID,i,j, 1:pt_num,:) = generate_near_GP5index(icgrid, rcID,i,j)
          else
            pt_num = 6
            GP_IJ(rcID,i,j, 1:pt_num,:) = generate_near_GP6index(rcID,i,j, EMin, EMax)
          end if
      end do
    end do
  end do

  call move_GP_with_spring_dyn(icgrid, idMin, idMax, EMin, EMax)
  deallocate(GP_IJ)

end subroutine modify_STDGrid_with_spring_dyn

!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
subroutine move_GP_with_spring_dyn( &
  & icgrid,                           & ! (inout)
  & idMin, idMax, EMin, EMax )        ! (in)

  ! 宣言文 ; Declaration statements
  !

  type(IcGrid2D_FVM), intent(inout) :: icgrid
  integer, intent(in) :: idMin,  idMax
  integer, intent(in) :: EMin,  EMax


  ! 作業変数
  ! Work variables
  !
  integer :: tstep
  integer rcID, i,j, k, convrecID
  real(DP) :: GP(6,3)
  integer :: pt_num
  real(DP) :: w0(RC_REGIONS_NUM, idMin:idMax, idMin:idMax, 1:3)

  integer, parameter :: DCONV_REC_NUM = 20
  integer, parameter :: DCONV_PLUS_ULIMIT = 10
  real(DP) :: conv_max, conv, conv_rec, conv_total
  real(DP) :: dconv_recs(DCONV_REC_NUM);
  integer :: dconv_plus_counter = 0
  logical :: conv_flag = .false.

  ! 実行文 ; Executable statements
  !

  !
  tstep = 1
  ! 初期速度の設定
  w0(:,:,:,:) = 0.0d0

  do
    do rcID=1, RC_REGIONS_NUM
      !$omp parallel do private(i, k, GP, pt_num)
      do j=EMin,EMax
        do i=EMin,EMax
          if( ( rcID > 5 .and. i==EMin .and. j==EMin ) .or. ( i==EMin .and. j==EMax ) .or. &
            & ( i==EMax .and. j==EMin ) .or. ( rcID <= 5 .and. i==EMax .and. j==EMax )     &
            &       ) then
             ! 特異点
             ! 格子点(rcID,i,j)の近傍にある格子点 5 個の格子点を求める.
            pt_num = 5
          else
             ! 格子点(rcID,i,j)の近傍にある格子点 6 個の格子点を求める.
            pt_num = 6
          end if

          if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
            do k=1,pt_num
              GP(k,1:3) = icgrid%rcs_AGrid(rcID, GP_IJ(rcID,i,j,k,1), GP_IJ(rcID,i,j,k,2), 1:3)
            end do

             ! 常微分方程式系を一タイムスッテプ分時間積分して, 一タイムスッテプ後の格子点位置を計算する.
            call integ_ode_1step_spring_dynamics( &
                 & r0=icgrid%rcs_AGrid(rcID,i,j,:), w0=w0(rcID,i,j,:), &
                 & GP=GP(1:pt_num,:), pt_num=pt_num &
                 & )

            call project_on_sphere( radius, icgrid%rcs_AGrid(rcID,i,j,:))
          end if

        end do ! for i
      end do ! for j
    end do ! for rcID

    ! 次のタイムステップに移る前に, のりしろを貼り合わせる.
    call paste_margin_width(icgrid)

    ! 収束判定
    if ( mod(tstep, conv_check_step) == 0 ) then

      !
      call check_global_convergence(icgrid, conv_total, conv_max, EMin, EMax)

      convrecID = mod(int(tstep / conv_check_step), DCONV_REC_NUM)+1
      dconv_recs(convrecID) = ( conv_total - conv_rec ) / abs(conv_max)
      conv_rec = conv_total

      if ( int(tstep / conv_check_step) > DCONV_REC_NUM ) then
        if ( sum(dconv_recs) > 0.0d0 ) then
           dconv_plus_counter = dconv_plus_counter + 1
             if ( dconv_plus_counter > DCONV_PLUS_ULIMIT ) then
               conv_flag = .true.
             end if
        end if
      end if

    end if

    if ( mod(tstep, 500) == 0 ) then
      write(*,*) '-- tstep:', tstep, '/', end_tstep
      write(*,*) '* conv_max=', conv_max, ', dconv=', sum(dconv_recs) / DCONV_REC_NUM
      write(*,*) '* conv_total=', conv_total,  ', dconv_plus_counter=', dconv_plus_counter
     ! write(*,*) conv_total / ( RC_REGIONS_NUM * ( (EMax - EMin + 1)**2 -4 ) )
    end if

    if ( ( tstep_limiter .and. ( tstep > end_tstep ) ) .or. conv_flag ) then
      write(*,*) tstep, tstep_limiter, end_tstep
      exit
    end if

    tstep = tstep + 1
  end do ! 時間積分ループ

  ! 反復回数と収束の程度を表す値を出力する.
  !

  write(*,*) 'SPR-GC-grid generation: the number of interation=', tstep
  write(*,*) 'The value for checking convergence:', conv_max

end subroutine move_GP_with_spring_dyn

!
subroutine check_global_convergence( &
  & icgrid,                   &  ! (in)
  & conv_total_ave, conv_max, &  ! (out)
  & EMin, EMax                &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(in) :: icgrid
  real(DP), intent(out) :: conv_total_ave
  real(DP), intent(out) :: conv_max
  integer, intent(in) :: EMin
  integer, intent(in) :: EMax

  ! 作業変数
  ! Work variables
  !
  integer rcID, i, j, k
  real(DP) :: GP(6,3)
  real(DP) :: conv
  integer :: pt_num

  ! 実行文 ; Executable statements
  !

  conv_max = 0.0d0
  conv_total_ave = 0.0d0

  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do private(i, k, GP, conv, pt_num) reduction(+:conv_total_ave)
    do j=EMin,EMax
      do i=EMin, EMax
        if(  .not. ( ( i==EMin .and. j==EMin ) .or. ( i==EMin .and. j==EMax ) .or. &
          &         ( i==EMax .and. j==EMin ) .or. ( i==EMax .and. j==EMax ) )    &
          &         ) then

          pt_num = 6
          do k=1,pt_num
            GP(k,1:3) = icgrid%rcs_AGrid(rcID, GP_IJ(rcID,i,j,k,1), GP_IJ(rcID,i,j,k,2), 1:3)
          end do

          conv = check_convergence(icgrid%rcs_AGrid(rcID,i,j,:), GP(1:pt_num,:), pt_num)
          conv_total_ave = conv_total_ave + conv

          if( conv  > conv_max ) then
            conv_max = conv
          end if
        end if

      end do
    end  do
  end do

  conv_total_ave = conv_total_ave / ( RC_REGIONS_NUM * ( (EMax - EMin + 1)**2 -4 ) )

end subroutine check_global_convergence

!
function check_convergence( &
  & r0, pts, pt_num &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statement
  !
  real(8), intent(in) :: r0(3)    !
  integer, intent(in) :: pt_num
  real(8), intent(in) :: pts(pt_num, 3) !
  real(8) val

  ! 作業変数
  ! Work variables
  !
  real(8) vec(3)
  integer i

  ! 実行文 ; Executable statement
  !
  vec = 0.0d0

  do i=1, pt_num
    vec(:) = vec(:) + &
       ( geodesic_arc_length(radius, r0, pts(i,:) ) - eq_d ) * calc_ei(r0(1:3), pts(i,1:3))
  end do

  val = vec_length(vec)

end function check_convergence

!
function generate_near_GP6index( &
  & rcID, GP0_i, GP0_j, EMin, EMax &  ! (in)
  & ) result(GP6Id)

  ! 宣言文 ; Declaration statements
  !

  integer, intent(in) :: rcID
  integer, intent(in) :: GP0_i
  integer, intent(in) :: GP0_j
  integer :: EMin,  EMax
  integer GP6Id(6, 2)

  ! 作業変数
  ! Work variables
  !
  integer i,j
  integer :: GP6IdShape(2) = (/ 6, 2 /)

  ! 実行文 ; Executable statements
  !

  i = GP0_i; j = GP0_j
  if ( rcID <= 5 .and. j==EMin ) then
    GP6Id(:,:) = reshape( (/ i,i-1,i-1,i-1,i,i+1, j-1,j-1,j,j+1,j+1,j /), GP6IdShape)

  else if ( rcID <= 5 .and. i==EMin ) then
    GP6Id(:,:) = reshape( (/ i,i-1,i-1,i,i+1,i+1, j-1,j-1,j,j+1,j,j-1 /), GP6IdShape)

  else if ( rcID > 5 .and. i==EMax ) then
    GP6Id(:,:) = reshape( (/ i+1,i+1,i,i-1,i-1,i, j,j+1,j+1,j+1,j,j-1 /), GP6IdShape)

  else if ( rcID > 5 .and. j==EMax ) then
    GP6Id(:,:) = reshape( (/ i+1,i+1,i+1,i,i-1,i, j-1,j,j+1,j+1,j,j-1 /), GP6IdShape)

  else
    GP6Id(:,:) = reshape( (/ i,i-1,i-1,i,i+1,i+1, j-1,j,j+1,j+1,j,j-1 /), GP6IdShape)
  end if

end function generate_near_GP6index

!
function generate_near_GP5index( &
  & icgrid,            &  ! (inout)
  & rcID, GP0_i, GP0_j &  ! (in)
  & ) result(GP5Id)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid
  integer, intent(in) :: rcID
  integer, intent(in) :: GP0_i
  integer, intent(in) :: GP0_j
  integer GP5Id(5, 2)

  ! 作業変数
  ! Work variables
  !
  integer :: tmp_GP5Id(5,3,2)

  ! 実行文 ; Executable statements
  !

  tmp_GP5Id(:,:,:) = generate_CV5_GPindex(icgrid, GP0_i, GP0_j, rcID)
  GP5Id(1:5, 1:2) = tmp_GP5Id(1:5, 2, 1:2)

end function generate_near_GP5index

!
subroutine integ_ode_1step_spring_dynamics( &
  & r0, w0,    &  ! (inout)
  & GP, pt_num &  ! (in)
  & )

  !
  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(inout) :: r0(3)
  real(DP), intent(inout) :: w0(3)
  integer, intent(in) :: pt_num
  real(DP), intent(in) :: GP(pt_num,3)

  ! 作業変数
  ! Work variables
  !
  real(DP) r0_(3), w0_(3)
  real(DP) k1(2,3), k2(2,3), k3(2,3)

  ! 実行文 ; Executable statements
  !

  !
  ! 連立常微分方程式(Tomita etal, (2001) の式(20), (21) )を, 3 次の Runge=Kutta 法を使って解く.
  !

  ! k1 を求める.
  k1(1, 1:3) = f1(r0(:), w0(:))
  k1(2, 1:3) = f2(r0(:), w0(:), GP(:,:), pt_num, alpha, k, eq_d, radius)

  ! k2 を求める.
  r0_(:) = r0(:) + 0.5d0 * dt * k1(1,:)
  w0_(:) = w0(:) + 0.5d0 * dt * k1(2,:)
  k2(1, 1:3) = f1(r0_(:), w0_(:))
  k2(2, 1:3) = f2(r0_(:), w0_(:), GP(:,:), pt_num, alpha, k, eq_d, radius)

  ! k3 を求める.
  r0_(:) = r0(:) - dt * k1(1,:) + 2.0d0 * dt * k2(1,:)
  w0_(:) = w0(:) - dt * k1(2,:) + 2.0d0 * dt * k2(2,:)
  k3(1, 1:3) = f1(r0_(:), w0_(:))
  k3(2, 1:3) = f2(r0_(:), w0_(:), GP(:,:), pt_num, alpha, k, eq_d, radius)

  ! 質点の位置ベクトルと速度ベクトルを次のタイムレベルに更新する.
  r0(:) = r0(:) + ( k1(1,:) + 4.0d0 * k2(1,:) + k3(1,:) ) * dt / 6.0d0
  w0(:) = w0(:) + ( k1(2,:) + 4.0d0 * k2(2,:) + k3(2,:) ) * dt / 6.0d0

end subroutine integ_ode_1step_spring_dynamics

!
function f1( &
  & r0, w0 &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: r0(3)
  real(DP), intent(in) :: w0(3)
  real(DP) val(3)

  ! 実行文 ; Executable statements
  !
  val(1:3) = w0(1:3)

end function f1

!
function f2( &
  & r0, w0, pts, pt_num, alpha, k, eq_d, radius &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: r0(3)
  real(DP), intent(in) :: w0(3)
  integer, intent(in) :: pt_num
  real(DP), intent(in) :: pts(pt_num,3)
  real(DP), intent(in) :: alpha
  real(DP), intent(in) :: k
  real(DP), intent(in) :: eq_d
  real(DP), intent(in) :: radius
  real(8) val(1:3)

  ! 作業変数
  ! Work variables
  !
  real(8) tmpv(3)
  integer i

  ! 実行文 ; Executable statements
  !

  tmpv(:) = 0.0d0
  do i=1, pt_num
    tmpv(:) = tmpv(:) + &
          ( geodesic_arc_length(radius, r0(1:3), pts(i,1:3)) - eq_d ) * calc_ei(r0(:), pts(i,:))
  end do

  ! M=1
  val(1:3) = k * tmpv(1:3) - alpha * w0(1:3)

end function f2

!
function calc_ei( &
  & r0, ri &  ! (in)
  & ) result(e)

  ! 宣言文 ; Declaration statements
  !
  real(8), intent(in) :: r0(3)
  real(8), intent(in) :: ri(3)
  real(8) e(3)

  ! 実行文 ; Executable statements
  !
  e(:) = cross( cross(r0(:), ri(:) ), r0)
  call vec_normarize(e(:))

end function calc_ei

end module SPRGCGrid_Builder