IGModel-SW 1.0

src/Williamson1992_TestCase/class_TestCase3.f90

説明を見る。
00001 
00009 module class_TestCase3
00010   
00011   ! モジュール引用 ; Use statements
00012   !
00013 
00014   ! 種類型パラメータ
00015   ! Kind type parameter 
00016   !
00017   use dc_types, only: DP,  & ! 倍精度実数型. Double precision. 
00018     &                 TOKEN
00019 
00020   !
00021   !
00022   !
00023   use gtool_history
00024 
00025   ! * IGMBaseLib *
00026   ! *
00027 
00028   ! 数学
00029   ! Mathematic
00030   !
00031   use igmcore_math, only: PI
00032 
00033   ! 座標変換
00034   ! Coordinate conversion
00035   !
00036   use igmcore_coordinate_conversion, only: &
00037     & orth_to_geo_pos, geo_to_orth_pos, geo_to_orth_vec
00038 
00039   ! 正二十面体格子の管理クラス
00040   !
00041   !
00042   use IcGrid2D_FVM_Manager, only: &
00043     & IcGrid2D_FVM, &
00044     & get_EffSize_Min, get_EffSize_Max, get_glevel, get_IcRadius, check_pole, &
00045     & RC_REGIONS_NUM, NOT_POLE_FLAG, NORTH_POLE_FLAG, SOUTH_POLE_FLAG
00046 
00047   ! 正二十面体格子上に分布する物理場の管理クラス
00048   !
00049   !
00050   use Field_IcGrid2D_manager, only: &
00051     & Field_IcGrid2D, &
00052     & Field_IcGrid2D_Init, get_icgrid
00053   
00054   !
00055   !
00056   !
00057   use Field_Pattern_Builder, only: &
00058     & create_solid_rotation_field
00059 
00060   ! 正二十面体格子上の物理場における誤差ノルムチェック
00061   !
00062   !
00063   use Field_IcGrid2D_error_norm_check, only: &
00064     & numerical_error_norm1, numerical_error_norm2, &
00065     & numerical_error_norminfinity 
00066 
00067   !
00068   !
00069   !
00070   use IcGrid_ncWriter_mod, only: &
00071     & IcGrid_ncWriter, IcGrid_ncWriter_Init, &
00072     & open_ncFile, ncdef_simulation_parameter, ncdef_GridData, ncdef_FieldData, end_ncdef, &
00073     & write_GridData, write_FieldData, close_ncFile, &
00074     & increase_recorde_counter
00075 
00076 
00077   ! * IGModel-SW ****
00078   ! *
00079 
00080   ! パラメータ管理
00081   !
00082   !
00083   use param_manager, only: &
00084     & alpha, output_tick, Grav, Omega, &
00085     & integration_time, DelTime
00086 
00087   ! OpenMP
00088   !
00089   !
00090   use omp_lib
00091 
00092   ! 宣言文 ; Declaration statement
00093   !
00094   implicit none
00095   private
00096   
00097   ! 公開変数
00098   ! Public variables
00099   !
00100   
00103   real(DP) :: h_0
00104 
00107   real(DP) :: u_0
00108 
00111   real(DP) :: angular_speed
00112 
00113   ! クラス定義
00114   ! Class definition
00115   !
00116 
00117   !
00121   type, public :: TestCase3
00122 
00123     ! 宣言文 ; Declaration statement
00124     !
00125 
00126     ! 非公開メンバ変数
00127     ! Private member variable
00128     !
00129     type(Field_IcGrid2D) :: true_h
00130     type(Field_IcGrid2D) :: error_h
00131     type(Field_IcGrid2D) :: true_v
00132     type(Field_IcGrid2D) :: error_v
00133     integer :: errorh_ncVarID
00134     integer :: errorV_ncVarID
00135 
00136     ! f2003 code
00137     !type(Field_IcGrid2D), private :: true_h
00138     !type(Field_IcGrid2D), private :: error_h
00139     !type(Field_IcGrid2D), private :: true_v
00140     !type(Field_IcGrid2D), private :: error_v
00141     !integer, private :: errorh_ncVarID
00142     !integer, private :: errorV_ncVarID
00143 
00144   end type TestCase3
00145 
00146   !
00147   !
00148   !
00149   interface initialize_TestCase
00150     module procedure init_TestCase3
00151   end interface initialize_TestCase
00152 
00153   interface finalize_TestCase
00154     module procedure finalize_TestCase3
00155   end interface finalize_TestCase
00156 
00157   public :: initialize_TestCase, finalize_TestCase
00158   public :: set_initial_h, set_initial_hs, set_initial_v, timelevel_Updated
00159 
00160 
00161   ! 非公開変数
00162   ! Private variables
00163   !
00164   
00165   character(TOKEN) :: filename = 'error_norm.dat'
00166   type(IcGrid_ncWriter), save :: writer
00167 
00168 contains
00169 
00170 !
00176 subroutine init_TestCase3(self, icgrid_ref )
00177   ! 宣言文 ; Declaration statement
00178   !
00179   type(TestCase3), intent(inout) :: self
00180   type(IcGrid2D_FVM), intent(in) :: icgrid_ref
00181 
00182   ! 作業変数
00183   ! Work variables
00184   !
00185   real(DP) :: ic_radius
00186 
00187   ! 実行文 ; Executable statement
00188   !
00189 
00190   ic_radius = get_IcRadius(icgrid_ref)
00191   
00192  
00193   u_0 = 2.0d0 * PI * ic_radius / ( 12.0d0 * 24.0d0 * 3600d0 )
00194   h_0 = 2.94e04 / Grav
00195   angular_speed = 2.0d0 * PI / ( 12.0d0 * 24.0d0 * 3600d0 )
00196 
00197   !
00198   call Field_IcGrid2D_Init(self%true_h, icgrid_ref, 'true_h', 1)
00199   call Field_IcGrid2D_Init(self%error_h, icgrid_ref, 'error_h', 1)
00200 
00201   !
00202   call Field_IcGrid2D_Init(self%true_v, icgrid_ref, 'true_v', 3)
00203   call Field_IcGrid2D_Init(self%error_v, icgrid_ref, 'error_v', 3)
00204 
00205   !
00206   call create_fields(self%true_h, self%true_v)
00207 
00208   !
00209   call HistoryCreate( &
00210     & file='error_norm.nc', title='error norm check', &
00211     & source=' ', institution=' ', &
00212     & dims=(/ 't' /), dimsizes=(/ 0 /), &
00213     & longnames=(/ 'time' /), units=(/ 's' /), &
00214     & origin=real(0d0), interval=real(output_tick) )
00215   
00216   call HistoryAddVariable( &
00217     & varname='v_norm1', dims=(/ 't' /), &
00218     & longname='v_norm1', units=' ', xtype='double' )
00219 
00220   call HistoryAddVariable( &
00221     & varname='v_norm2', dims=(/ 't' /), &
00222     & longname='v_norm2', units=' ', xtype='double' )
00223 
00224   call HistoryAddVariable( &
00225     & varname='v_norm_infinity', dims=(/ 't' /), &
00226     & longname='v_norm_infinity', units=' ', xtype='double' )
00227 
00228   !
00229   call HistoryAddVariable( &
00230     & varname='h_norm1', dims=(/ 't' /), &
00231     & longname='h_norm1', units=' ', xtype='double' )
00232 
00233   call HistoryAddVariable( &
00234     & varname='h_norm2', dims=(/ 't' /), &
00235     & longname='h_norm2', units=' ', xtype='double' )
00236 
00237   call HistoryAddVariable( &
00238     & varname='h_norm_infinity', dims=(/ 't' /), &
00239     & longname='h_norm_infinity', units=' ', xtype='double' )
00240 
00241 
00242   call HistoryAddAttr('v_norm1', '+glevel', get_glevel(icgrid_ref))
00243   call HistoryAddAttr('v_norm1', '+alpha', alpha)
00244   
00245   !
00246   !
00247   call IcGrid_ncWriter_Init(writer, icgrid_ref)
00248   call open_ncFile(writer, 'h_error_check.nc')
00249   call ncdef_Simulation_Parameter(writer, integration_time, DelTime, output_tick)
00250 
00251   ! 格子点情報の定義 
00252   call ncdef_GridData(writer)
00253   self%errorh_ncVarID = ncdef_FieldData(writer, self%error_h)
00254   self%errorv_ncVarID = ncdef_FieldData(writer, self%error_v)
00255 
00256   ! 定義の終了. 
00257   call end_ncdef(writer)
00258 
00259   ! 格子情報を netCDF に書きこむ. 
00260   !
00261   call write_GridData(writer)
00262 
00263 end subroutine init_TestCase3
00264 
00265 !
00270 subroutine finalize_TestCase3(self)
00271   ! 宣言文 ; Declaration statement
00272   !
00273   type(TestCase3), intent(inout) :: self
00274 
00275   ! 実行文 ; Executable statement
00276   !
00277 
00278   !
00279   call HistoryClose()
00280   call close_ncFile(writer)
00281 
00282 end subroutine finalize_TestCase3
00283 
00284 !
00290 subroutine set_initial_v(self, init_v)
00291   ! 宣言文 ; Declaration statement
00292   !
00293   type(TestCase3), intent(inout) :: self
00294   type(Field_IcGrid2D), intent(inout) :: init_v
00295 
00296   ! 実行文 ; Executable statement
00297   !
00298 
00299   write(*,*) 'init_TestCase3_v'
00300 
00301   init_v%val(:,:,:,:) =  self%true_v%val(:,:,:,:)
00302 
00303 end subroutine set_initial_v
00304 
00305 !
00311 subroutine set_initial_h(self, init_h)
00312   ! 宣言文 ; Declaration statement
00313   !
00314   type(TestCase3), intent(inout) :: self
00315   type(Field_IcGrid2D), intent(inout) :: init_h
00316 
00317   ! 実行文 ; Executable statement
00318   !
00319   
00320   write(*,*) 'init_TestCase3_h'
00321 
00322   init_h%val(:,:,:,:) = self%true_h%val(:,:,:,:)
00323 
00324 end subroutine set_initial_h
00325 
00326 !
00332 subroutine set_initial_hs(self, init_hs)
00333 
00334   ! 宣言文 ; Declaration statement
00335   !
00336   type(TestCase3), intent(inout) :: self
00337   type(Field_IcGrid2D), intent(inout) :: init_hs
00338 
00339   ! 実行文 ; Executable statement
00340   !
00341 
00342   write(*,*) 'TestCase3: initialize height surface field'  
00343   init_hs%val(:,:,:,:) = 0.0d0
00344 
00345 end subroutine set_initial_hs
00346 
00347 !
00356 subroutine timelevel_Updated(self, tstep, dt, v_n, h_n)
00357 
00358   ! 宣言文 ; Declaration statement
00359   !
00360   type(TestCase3), intent(inout) :: self
00361   integer, intent(in) :: tstep
00362   real(DP), intent(in) :: dt
00363   type(Field_IcGrid2D), intent(inout) :: v_n
00364   type(Field_IcGrid2D), intent(inout) :: h_n
00365 
00366 
00367   ! 実行文 ; Executable statement
00368   !
00369   
00370   if ( mod(tstep * dt, output_tick) == 0 ) then
00371   !
00372     call eval_numcal_h_solution( tstep*dt, &
00373       & self%true_h, self%true_v, h_n, v_n )
00374 
00375     !
00376     self%error_h%val(:,:,:,:) = self%true_h%val(:,:,:,:) - h_n%val(:,:,:,:)
00377     call write_FieldData(writer, self%errorh_ncVarID, self%error_h)
00378     call increase_recorde_counter(writer)
00379 
00380   end if
00381 
00382 end subroutine timelevel_Updated
00383 
00384 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00385 !
00386 ! 非公開手続き
00387 ! 
00388 !!!!!!!!!!!!!!!!!!!!!!!!!
00389 
00390 !
00396 subroutine eval_numcal_h_solution( &
00397   & t, &
00398   & true_h, true_v, &
00399   & h_n, v_n &
00400   & )
00401 
00402   ! 宣言文 ; Declaration statement
00403   !
00404   real(DP), intent(in) :: t
00405   type(Field_IcGrid2D), intent(in) :: true_h
00406   type(Field_IcGrid2D), intent(in) :: true_v
00407   type(Field_IcGrid2D), intent(in) :: h_n
00408   type(Field_IcGrid2D), intent(in) :: v_n
00409   
00410   ! 作業変数
00411   ! Work variables
00412   !
00413   real(DP) norm1Val, norm2Val, normInfityVal
00414 
00415   ! 実行文 ; Executable statement
00416   !
00417   
00418   !
00419   norm1Val = numerical_error_norm1(true_h, h_n)
00420   norm2Val = numerical_error_norm2(true_h, h_n)
00421   normInfityVal = numerical_error_norminfinity(true_h, h_n, maxmin_info_flag=.false.)
00422 
00423   call HistoryPut('h_norm1', norm1Val)
00424   call HistoryPut('h_norm2', norm2Val)
00425   call HistoryPut('h_norm_infinity', normInfityVal)
00426 
00427   norm1Val = numerical_error_norm1(true_v, v_n)
00428   norm2Val = numerical_error_norm2(true_v, v_n)
00429   normInfityVal = numerical_error_norminfinity(true_v, v_n, maxmin_info_flag=.false.)
00430 
00431   call HistoryPut('v_norm1', norm1Val)
00432   call HistoryPut('v_norm2', norm2Val)
00433   call HistoryPut('v_norm_infinity', normInfityVal)
00434 
00435 end subroutine eval_numcal_h_solution
00436 
00437 subroutine create_fields(h_field, v_field)
00438   ! モジュール引用 ; Use statements
00439   !
00440 
00441   use igmcore_math, only: toDegrees
00442   use igmcore_linear_algebra, only: rotateY
00443 
00444   ! 宣言文 ; Declaration statements
00445   !
00446   type(Field_IcGrid2D), intent(inout) :: h_field
00447   type(Field_IcGrid2D), intent(inout) :: v_field
00448 
00449   ! 作業変数
00450   ! Work variables
00451   !
00452   type(IcGrid2D_FVM), pointer :: icgrid
00453   integer :: rcID, i, j
00454   integer :: EMin, EMax
00455   integer :: pole_flag
00456   integer :: ndiv
00457   
00458   real(DP) :: u, v, u_ds
00459   real(DP) :: ic_radius, s_pos(3), s_pos_dash(3)
00460   integer, parameter :: MAX_GLEVEL = 9
00461   real(DP) :: integ_dtheta
00462 
00463   !
00464  
00465   ! 実行文 ; Executable statements
00466   !
00467   
00468   !
00469   icgrid => get_icgrid(v_field)
00470   EMin = get_EffSize_Min(icgrid)
00471   EMax = get_EffSize_Max(icgrid)
00472   ic_radius = get_IcRadius(icgrid)
00473   s_pos_dash(3) = ic_radius
00474 
00475   !
00476   integ_dtheta = 2.0d0 * PI / ( 10.0d0 * 2**(MAX_GLEVEL - 1) ) / 100.0
00477 
00478   !
00479   write(*,*) 'Create the fields ..'
00480 
00481   do rcID=1, RC_REGIONS_NUM
00482     !$omp parallel do private(s_pos, s_pos_dash, u_ds, u, v, ndiv)
00483     do j=EMin, EMax
00484       do i=EMin, EMax
00485         pole_flag = check_pole(icgrid, rcID,i,j)
00486         if(  pole_flag == NORTH_POLE_FLAG ) then
00487           s_pos(:) = (/ PI, PI / 2.0d0, ic_radius /)
00488         else if( pole_flag == SOUTH_POLE_FLAG ) then
00489           s_pos(:) = (/ 0.0d0, - PI / 2.0d0, ic_radius /)
00490         else
00491           s_pos(:) = orth_to_geo_pos(icgrid%rcs_AGrid(rcID,i,j,:))
00492         end if
00493 
00494         !
00495         s_pos_dash(2) = &
00496           & dasin( &
00497           &   dsin(s_pos(2)) * dcos(alpha) + dsin(alpha) * dcos(s_pos(1)) * dcos(s_pos(2)) &
00498           & )
00499 
00500         s_pos_dash(1) = dasin( dsin(s_pos(1)) * dcos(s_pos(2)) / dcos(s_pos_dash(2)) )
00501 
00502         !
00503         if ( dcos(alpha) * dcos(s_pos(1)) * dcos(s_pos(2)) - dsin(alpha) * dsin(s_pos(2)) < 0.0d0 ) then
00504           s_pos_dash(1) = PI - s_pos_dash(1)         
00505         end if
00506 
00507         !
00508         u_ds = calc_u_dash(s_pos_dash(2)) 
00509         v = u_ds * dsin(alpha) * dsin(s_pos_dash(1)) / dcos(s_pos(2))
00510         u = v * dsin(s_pos(2)) * dtan(s_pos(1)) + u_ds * dcos(s_pos_dash(1)) / dcos(s_pos(1))
00511 
00512         if ( (abs( s_pos(2)) - PI / 2.0)**2 < 1.0e-20 ) then
00513           v_field%val(rcID,i,j,:) = geo_to_orth_vec((/ u, 0.0d0, 0.0d0 /), &
00514             &                        (/ PI, s_pos_dash(2), ic_radius /) )
00515         else
00516           !
00517           v_field%val(rcID,i,j, :) = geo_to_orth_vec((/ u, v, 0.0d0 /), s_pos)
00518         end if
00519 
00520         ndiv = int( ( s_pos_dash(2) - (- 0.5d0 * PI ) ) / integ_dtheta )
00521 
00522         !
00523         h_field%val(rcID,i,j,:) = calc_h(s_pos_dash(2), ic_radius, ndiv)
00524       end do
00525     end do
00526   end do
00527 
00528   write(*,*) 'Finish creating the fields'
00529 
00530 end subroutine create_fields
00531 
00532 function calc_h(theta_dash, ic_radius, n) result(val)
00533   real(DP), intent(in) :: theta_dash
00534   real(DP), intent(in) :: ic_radius
00535   integer, intent(in) :: n
00536   real(DP) val
00537 
00538   !
00539   !
00540   !
00541   real(DP) :: integ_ret
00542   real(DP) :: tau_lbound, tau
00543   real(DP) :: h, h2
00544   real(DP) :: e_sum, o_sum
00545   integer :: i
00546 
00547   !
00548   !
00549   if( n == 0 ) then
00550     val = h_0
00551     return
00552   end if
00553 
00554   tau_lbound = - PI / 2.0d0
00555   h = ( theta_dash - tau_lbound ) / dble(n)
00556 
00557   e_sum = 0.0d0
00558   do i=2, n-2, 2
00559     tau = tau_lbound + i * h
00560     e_sum = e_sum + h_integrand(tau, ic_radius)
00561   end do
00562 
00563   o_sum = 0.0d0
00564   do i=1, n-1, 2
00565     tau = tau_lbound + i * h
00566     o_sum = e_sum + h_integrand(tau, ic_radius)
00567   end do
00568 
00569   integ_ret =   h * ( h_integrand(tau_lbound, ic_radius) &
00570     &         + h_integrand(theta_dash, ic_radius) + 2.0d0 * e_sum + 4.0d0 * o_sum ) / 3.0d0
00571 
00572   val = h_0 - ic_radius * integ_ret / Grav
00573        
00574 end function calc_h
00575 
00576 function h_integrand(tau, a) result(val)
00577   !
00578   !
00579   real(DP), intent(in) :: tau
00580   real(DP), intent(in) :: a
00581   real(DP) val
00582 
00583   !
00584   !
00585   real(DP) :: u_ds
00586 
00587   !
00588   !
00589   u_ds = calc_u_dash(tau)
00590 
00591   val = u_ds * ( 2.0d0 * Omega * dsin(tau) + u_ds * dtan(tau) / a )
00592   
00593 end function h_integrand
00594 
00595 function calc_u_dash(theta_dash) result(val)
00596   real(DP), intent(in) :: theta_dash
00597   real(DP) val
00598 
00599   !
00600   !
00601   !
00602   real(DP), parameter :: theta_b = - PI / 6.0d0
00603   real(DP), parameter :: theta_e = PI / 2.0d0
00604   real(DP), parameter :: x_e = 0.3d0
00605   real(DP) :: x
00606 
00607   !
00608   !
00609   x = x_e * ( theta_dash - theta_b ) / ( theta_e - theta_b )
00610    
00611   val = u_0 * b(x) * b(x_e - x) * dexp(4.0d0 / x_e)
00612        
00613 end function calc_u_dash
00614 
00615 function b(x) result(val)
00616   !
00617   !
00618   real(DP), intent(in) :: x
00619   real(DP) val
00620 
00621   !
00622   !
00623 
00624   if( x <= 0.0d0 ) then
00625     val = 0.0d0
00626   else
00627     val = dexp(-1.0d0/x) 
00628   end if
00629 
00630 end function b
00631 
00632 end module class_TestCase3
 全て クラス ネームスペース ファイル 関数 変数