|
IGModel-SW 1.0
|
00001 00013 module sw_equation_solver 00014 00015 ! モジュール引用 ; Use statements 00016 ! 00017 00018 !* gtool ** 00019 !* 00020 00021 ! 種類型パラメータ 00022 ! Kind type parameter 00023 ! 00024 use dc_types, only: DP ! 倍精度実数型. Double precision. 00025 00026 !* IGMBaseLib ** 00027 !* 00028 00029 ! 正二十面体格子の管理 00030 ! Manager of icosahedral grid 00031 ! 00032 use IcGrid2D_FVM_Manager, only: & 00033 & IcGrid2D_FVM, & 00034 & get_EffSize_Min, get_EffSize_Max, get_idMin 00035 00036 ! 正二十面体格子上の物理場データの管理 00037 ! Manager of the physical field data on icosahedral grid 00038 ! 00039 use Field_IcGrid2D_Manager, only: & 00040 & Field_IcGrid2D, & 00041 & Field_IcGrid2D_Init, Field_IcGrid2D_Final, paste_field2D_margin 00042 00043 00044 !* IGModel-SW ** 00045 !* 00046 00047 ! 物理場の管理 00048 ! Manager of the physical fields 00049 ! 00050 use field_manager, only: & 00051 & xy_Vel_TL_list, xy_Height_TL_list, xy_Htopo, & 00052 & DVelDt_TL_list, DHeightDt_TL_list, & 00053 & TL_N, TL_Nplus1, & 00054 & TL_DDT_N, TL_DDT_Nminus1, TL_DDT_Nminus2, & 00055 & xy_Coli, diff_eval 00056 00057 ! 運動方程式の評価 00058 ! Evaluate motion equation 00059 ! 00060 use motion_equation, only: & 00061 & motion_equation_Init, & 00062 & calc_motion_eq_dvdt 00063 00064 ! 連続の式の評価 00065 ! Evaluate continious equation 00066 ! 00067 use continious_equation, only: & 00068 & continious_equation_Init, & 00069 & calc_continious_eq_dhdt 00070 00071 ! シミュレーションパラメータの管理 00072 ! Manager of simulation parameter 00073 ! 00074 use param_manager, only: motionEq_flag 00075 00076 ! 宣言文 ; Declaration statement 00077 ! 00078 implicit none 00079 private 00080 00081 ! 公開手続き 00082 ! Public procedure 00083 ! 00084 public :: init_sw_equation_solver, temporal_integration 00085 00086 ! 非公開変数 00087 ! Private variables 00088 ! 00089 00090 contains 00091 00092 ! 00103 subroutine init_sw_equation_solver( & 00104 & icgrid & ! (inout) 00105 & ) 00106 ! 宣言文 ; Declaration statements 00107 ! 00108 type(IcGrid2D_FVM), intent(inout) :: icgrid 00109 00110 ! 実行文 ; Executable statements 00111 ! 00112 00113 !write(*,*) 'init sw_equation_solver' 00114 00115 ! 運動方程式・連続の式の評価用モジュールの初期化. 00116 call motion_equation_Init(icgrid) 00117 call continious_equation_Init(icgrid) 00118 00119 end subroutine init_sw_equation_solver 00120 00121 00122 ! 00137 subroutine temporal_integration( & 00138 & tstep, dt, & ! (in) 00139 & icgrid & ! (inout) 00140 ) 00141 00142 ! 宣言文 ; Declaration statement 00143 ! 00144 integer, intent(in) :: tstep 00145 real(DP), intent(in) :: dt 00146 type(IcGrid2D_FVM), intent(in) :: icgrid 00147 00148 ! 作業変数 00149 ! Work variables 00150 ! 00151 integer :: idMin, EMin, EMax 00152 00153 ! 実行文 ; Executable statement 00154 ! 00155 00156 00157 ! 最初の 2 Step は, 4 次の Runge=Kutta 法を使って時間積分する. 00158 ! それ以降は, 3 次の Adams=Bashforth 法を使って時間積分する. 00159 ! 00160 if( tstep <= 2 ) then 00161 idMin = get_IdMin(icgrid) 00162 00163 ! 4 次の Runge=Kutta 法により, 浅水方程式を時間積分する. 00164 ! 各時間スッテプにおいて, dvdt, dhdt を記憶しておく. 00165 call RungeKutta_fourth_order( & 00166 & xy_Vel_TL_list(TL_Nplus1)%val, xy_Vel_TL_list(TL_N)%val, DVelDt_TL_list(TL_DDT_N)%val, & 00167 & xy_Height_TL_list(TL_Nplus1)%val, xy_Height_TL_list(TL_N)%val, DHeightDt_TL_list(TL_DDT_N)%val, & 00168 & dt, icgrid, idMin ) 00169 00170 else 00171 00172 ! 3 次の Adams=Bashforth 法により, 浅水方程式を時間積分する. 00173 ! 運動方程式を時間積分して, タイムレベル n+1 の速度場を求める. 00174 ! motionEq_flag が偽のときは, 運動方程式の時間積分は行われない. 00175 ! 00176 00177 EMin = get_EffSize_Min(icgrid) 00178 EMax = get_EffSize_Max(icgrid) 00179 00180 if ( motionEq_flag ) then 00181 ! 速度場の時間変化率を計算する. 00182 call calc_motion_eq_dvdt( & 00183 & DVelDtN = DVelDt_TL_list(TL_DDT_N), & 00184 & xy_VelN = xy_Vel_TL_list(TL_N), & 00185 & xy_HeightN = xy_Height_TL_list(TL_N), & 00186 & xy_Coli = xy_Coli, & 00187 & diff_eval = diff_eval ) 00188 end if 00189 00190 ! 高度場の時間変化率を計算する. 00191 call calc_continious_eq_dhdt( & 00192 & DHeightDtN = DHeightDt_TL_list(TL_DDT_N), & 00193 & xy_VelN = xy_Vel_TL_list(TL_N), & 00194 & xy_HeightN = xy_Height_TL_list(TL_N), & 00195 & xy_Htopo = xy_Htopo, & 00196 & diff_eval=diff_eval ) 00197 00198 !write(*,*) tstep, ':dh=', dhdt_TL_list(TL_DDT_N)%val(10,EMax,EMax,1)*dt 00199 00200 !$omp sections 00201 00202 !$omp section 00203 if ( motionEq_flag ) then 00204 ! 速度場に対する時間積分を行う. 00205 call Adams_Bashforth_third_order( & 00206 & U_nplus1 = xy_Vel_TL_list(TL_Nplus1)%val(:,EMin:EMax,EMin:EMax,:), & 00207 & U_n = xy_Vel_TL_list(TL_N)%val(:,EMin:EMax,EMin:EMax,:), & 00208 & dUdt_n = DVelDt_TL_list(TL_DDT_N)%val(:,EMin:EMax,EMin:EMax,:), & 00209 & dUdt_nminus1 = DVelDt_TL_list(TL_DDT_Nminus1)%val(:,EMin:EMax,EMin:EMax,:), & 00210 & dUdt_nminus2 = DVelDt_TL_list(TL_DDT_Nminus2)%val(:,EMin:EMax,EMin:EMax,:), & 00211 & dt = dt ) 00212 00213 end if 00214 00215 !$omp section 00216 ! 高度場に対する時間積分を行う. 00217 call Adams_Bashforth_third_order(& 00218 & U_nplus1 = xy_Height_TL_list(TL_Nplus1)%val(:,EMin:EMax,EMin:EMax,:), & 00219 & U_n = xy_Height_TL_list(TL_N)%val(:,EMin:EMax,EMin:EMax,:), & 00220 & dUdt_n = DHeightDt_TL_list(TL_DDT_N)%val(:,EMin:EMax,EMin:EMax,:), & 00221 & dUdt_nminus1 = DHeightDt_TL_list(TL_DDT_Nminus1)%val(:,EMin:EMax,EMin:EMax,:), & 00222 & dUdt_nminus2 = DHeightDt_TL_list(TL_DDT_Nminus2)%val(:,EMin:EMax,EMin:EMax,:), & 00223 & dt = dt ) 00224 00225 !$omp end sections 00226 00227 end if 00228 00229 end subroutine temporal_integration 00230 00231 00232 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00233 ! 00234 ! 非公開手続き 00235 ! Private procedures 00236 ! 00237 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00238 00239 ! 00299 subroutine RungeKutta_fourth_order( & 00300 & xy_VelA, & ! (inout) 00301 & xy_VelN, DVelDtN, & ! (in) 00302 & xy_HeightA, & ! (inout) 00303 & xy_HeightN, DHeightDtN, & ! (in) 00304 & dt, icgrid, idMin & ! (in) 00305 & ) 00306 00307 ! 宣言文 ; Declaration statements 00308 ! 00309 integer, intent(in) :: idMin 00310 real(DP), intent(inout) :: xy_VelA(:,idMin:,idMin:,:) 00311 real(DP), intent(in) :: xy_VelN(:,idMin:,idMin:,:) 00312 real(DP), intent(inout) :: DVelDtN(:,idMin:,idMin:,:) 00313 real(DP), intent(inout) :: xy_HeightA(:,idMin:,idMin:,:) 00314 real(DP), intent(in) :: xy_HeightN(:,idMin:,idMin:,:) 00315 real(DP), intent(inout) :: DHeightDtN(:,idMin:,idMin:,:) 00316 real(DP), intent(in) :: dt 00317 type(IcGrid2D_FVM), intent(in) :: icgrid ! IcGrid2D_FVM クラスのインスタンスの参照. 00318 00319 ! 作業変数 00320 ! Work variables 00321 ! 00322 type(Field_IcGrid2D) :: k_v(4) 00323 type(Field_IcGrid2D) :: v_predictor 00324 type(Field_IcGrid2D) :: k_h(4) 00325 type(Field_IcGrid2D) :: h_predictor 00326 00327 integer, parameter :: K_NUM = 4 00328 real(DP), parameter :: k_coef(3) = (/ 0.5d0, 0.5d0, 1.0d0 /) 00329 integer :: m,n 00330 00331 ! 実行文 ; Executable statements 00332 ! 00333 00334 ! 作業配列を初期化する. 00335 ! Initializes the work variables. 00336 00337 do m=1, K_NUM 00338 call Field_IcGrid2D_Init( & 00339 & k_v(m), & 00340 & icgrid=icgrid, name='k_v', rank=3 ) 00341 00342 call Field_IcGrid2D_Init( & 00343 & k_h(m), & 00344 & icgrid=icgrid, name='k_h', rank=1 ) 00345 00346 end do 00347 00348 call Field_IcGrid2D_Init( & 00349 & self=v_predictor, & 00350 & icgrid=icgrid, name='v_predictor', rank=3 ) 00351 00352 call Field_IcGrid2D_Init( & 00353 & self=h_predictor, & 00354 & icgrid=icgrid, name='h_predictor', rank=1 ) 00355 00356 ! k_v(1), k_h(1) をまず求めるために, 半離散式 L の右辺の関数の引数 v_predictor, h_predictor として, 00357 ! v_n, h_n を与える. 00358 v_predictor%val(:,:,:,:) = xy_VelN(:,:,:,:) 00359 h_predictor%val(:,:,:,:) = xy_HeightN(:,:,:,:) 00360 00361 do m=1, K_NUM 00362 call paste_field2D_margin(v_predictor) 00363 call paste_field2D_margin(h_predictor) 00364 00365 ! (t_n, v_predictor, h_predictor) を引数とする運動方程式の半離散式を評価し, 00366 ! k_v(m) を求める. 00367 00368 call calc_motion_eq_dvdt( & 00369 & DVelDtN=k_v(m), & 00370 & xy_VelN=v_predictor, xy_HeightN=h_predictor, xy_Coli=xy_Coli, & 00371 & diff_eval=diff_eval ) 00372 00373 ! (t_n, v_predictor, h_predictor, h_s) を引数とする運動方程式の半離散式を評価し, 00374 ! k_h(m) を求める. 00375 call calc_continious_eq_dhdt( & 00376 & DHeightDtN=k_h(m), & 00377 & xy_VelN=v_predictor, xy_HeightN=h_predictor, xy_Htopo=xy_Htopo, & 00378 & diff_eval=diff_eval ) 00379 00380 00381 ! k_v(m+1), k_h(m+1) (ただし, 2 <= m < 4) を求めるために, 00382 ! 半離散式 L の右辺の関数の引数 v_predictor, h_predictor を更新する. 00383 00384 if( m < K_NUM )then 00385 00386 ! 速度場に対して 00387 if ( motionEq_flag ) then 00388 v_predictor%val(:,:,:,:) = xy_VelN(:,:,:,:) + k_coef(m) * dt * k_v(m)%val(:,:,:,:) 00389 else 00390 ! 運動方程式を使用しない場合は, 予測子にタイムレベル n 速度場を代入する. 00391 v_predictor%val(:,:,:,:) = xy_VelN(:,:,:,:) 00392 end if 00393 00394 ! 表面高度場に対して 00395 h_predictor%val(:,:,:,:) = xy_HeightN(:,:,:,:) + k_coef(m) * dt * k_h(m)%val(:,:,:,:) 00396 00397 end if 00398 00399 end do 00400 00401 ! 求めた k_v(1:4), k_h(1:4) を使って, タイムレベル n の局所時間変化率 dh/dt, dv/dt を求める. 00402 ! 00403 00404 DVelDtN(:,:,:,:) = ( k_v(1)%val(:,:,:,:) + 2.0d0 * k_v(2)%val(:,:,:,:) & 00405 & + 2.0d0 * k_v(3)%val(:,:,:,:) + k_v(4)%val(:,:,:,:) & 00406 & ) / 6.0d0 00407 00408 DHeightDtN(:,:,:,:) = ( k_h(1)%val(:,:,:,:) + 2.0d0 * k_h(2)%val(:,:,:,:) & 00409 & + 2.0d0 * k_h(3)%val(:,:,:,:) + k_h(4)%val(:,:,:,:) & 00410 & ) / 6.0d0 00411 00412 ! タイムレベル n+1 の速度場, 高度場を求める. 00413 ! 00414 00415 xy_VelA(:,:,:,:) = xy_VelN(:,:,:,:) + DVelDtN(:,:,:,:) * dt 00416 xy_HeightA(:,:,:,:) = xy_HeightN(:,:,:,:) + DHeightDtN(:,:,:,:) * dt 00417 00418 00419 ! 作業配列を破棄する. 00420 ! 00421 00422 do m=1, K_NUM 00423 call Field_IcGrid2D_Final(k_v(m)) 00424 call Field_IcGrid2D_Final(k_h(m)) 00425 end do 00426 00427 end subroutine RungeKutta_fourth_order 00428 00429 ! 00469 subroutine Adams_Bashforth_third_order( & 00470 & U_nplus1, & ! (inout) 00471 & U_n, dUdt_n, dUdt_nminus1, dUdt_nminus2, & ! (in) 00472 & dt ) ! (in) 00473 00474 ! 宣言文 ; Declaration statement 00475 ! 00476 real(DP), intent(inout) :: U_nplus1(:,:,:,:) 00477 real(DP), intent(in) :: U_n(:,:,:,:) 00478 real(DP), intent(in) :: dUdt_n(:,:,:,:) 00479 real(DP), intent(in) :: dUdt_nminus1(:,:,:,:) 00480 real(DP), intent(in) :: dUdt_nminus2(:,:,:,:) 00481 real(DP), intent(in) :: dt 00482 00483 ! 実行文 ; Executable statements 00484 ! 00485 00486 U_nplus1(:,:,:,:) = & 00487 & U_n(:,:,:,:) & 00488 & + dt * ( 23.0d0 * dUdt_n(:,:,:,:) - 16.0d0 * dUdt_nminus1(:,:,:,:) + 5.0d0 * dUdt_nminus2(:,:,:,:) ) / 12.0d0 00489 00490 end subroutine Adams_Bashforth_third_order 00491 00492 end module sw_equation_solver
1.7.3 
