00001 subroutine iteration_BH_WL(iter_count)
00002 use phys_constant, only : long, nnrg, nntg, nnpg
00003 use def_formulation, only : chgra, chope
00004 use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd
00005 use grid_parameter
00006 use coordinate_grav_r
00007 use coordinate_grav_phi
00008 use coordinate_grav_theta
00009 use weight_midpoint_grav
00010 use def_metric
00011 use make_array_2d
00012 use make_array_3d
00013 use make_array_4d
00014 use interface_interpo_fl2gr
00015 use interface_poisson_solver_1bh_homosol
00016 use interface_bh_boundary_analytic
00017 use interface_compute_shift
00018 use interface_compute_alps2alph
00019 use interface_update_grfield
00020
00021 use interface_source_HaC_WL_BH
00022 use interface_source_MoC_WL_BH
00023 use interface_source_trfreeG_WL_BH
00024 use interface_source_trG_WL_BH
00025 use interface_compute_fnc_multiple
00026 use interface_error_metric_type3
00027 implicit none
00028 real(long), pointer :: sou(:,:,:), pot(:,:,:)
00029 real(long), pointer :: sou1(:,:,:), sou2(:,:,:), sou3(:,:,:)
00030 real(long), pointer :: sou4(:,:,:), sou5(:,:,:), sou6(:,:,:)
00031 real(long), pointer :: pot_bak(:,:,:)
00032 real(long), pointer :: potx(:,:,:), poty(:,:,:), potz(:,:,:)
00033 real(long), pointer :: potx_bak(:,:,:), poty_bak(:,:,:), potz_bak(:,:,:)
00034 real(long), pointer :: souvec(:,:,:,:), souten(:,:,:,:)
00035 real(long), pointer :: pot1(:,:,:), pot2(:,:,:), pot3(:,:,:)
00036 real(long), pointer :: pot4(:,:,:), pot5(:,:,:), pot6(:,:,:)
00037 real(long), pointer :: bk1(:,:,:), bk2(:,:,:), bk3(:,:,:)
00038 real(long), pointer :: bk4(:,:,:), bk5(:,:,:), bk6(:,:,:)
00039 real(long), pointer :: work(:,:,:)
00040 real(long), pointer :: sou_bhsurf(:,:), dsou_bhsurf(:,:)
00041 real(long), pointer :: sou_outsurf(:,:), dsou_outsurf(:,:)
00042 real(long) :: count
00043 integer :: iter_count, flag = 0
00044 integer :: irf, itf, ipf, irg, itg, ipg, ihy
00045
00046 real(long) :: error_all(20)
00047 integer :: flag_all(20) = (/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/)
00048
00049 call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00050 call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00051 call alloc_array3d(potx,0,nrg,0,ntg,0,npg)
00052 call alloc_array3d(poty,0,nrg,0,ntg,0,npg)
00053 call alloc_array3d(potz,0,nrg,0,ntg,0,npg)
00054 call alloc_array4d(souvec,0,nrg,0,ntg,0,npg,1,3)
00055 call alloc_array4d(souten,0,nrg,0,ntg,0,npg,1,6)
00056 call alloc_array3d(sou1,0,nrg,0,ntg,0,npg)
00057 call alloc_array3d(sou2,0,nrg,0,ntg,0,npg)
00058 call alloc_array3d(sou3,0,nrg,0,ntg,0,npg)
00059 call alloc_array3d(sou4,0,nrg,0,ntg,0,npg)
00060 call alloc_array3d(sou5,0,nrg,0,ntg,0,npg)
00061 call alloc_array3d(sou6,0,nrg,0,ntg,0,npg)
00062 call alloc_array3d(pot1,0,nrg,0,ntg,0,npg)
00063 call alloc_array3d(pot2,0,nrg,0,ntg,0,npg)
00064 call alloc_array3d(pot3,0,nrg,0,ntg,0,npg)
00065 call alloc_array3d(pot4,0,nrg,0,ntg,0,npg)
00066 call alloc_array3d(pot5,0,nrg,0,ntg,0,npg)
00067 call alloc_array3d(pot6,0,nrg,0,ntg,0,npg)
00068 call alloc_array3d(pot_bak,0,nrg,0,ntg,0,npg)
00069 call alloc_array3d(potx_bak,0,nrg,0,ntg,0,npg)
00070 call alloc_array3d(poty_bak,0,nrg,0,ntg,0,npg)
00071 call alloc_array3d(potz_bak,0,nrg,0,ntg,0,npg)
00072 call alloc_array3d(bk1,0,nrg,0,ntg,0,npg)
00073 call alloc_array3d(bk2,0,nrg,0,ntg,0,npg)
00074 call alloc_array3d(bk3,0,nrg,0,ntg,0,npg)
00075 call alloc_array3d(bk4,0,nrg,0,ntg,0,npg)
00076 call alloc_array3d(bk5,0,nrg,0,ntg,0,npg)
00077 call alloc_array3d(bk6,0,nrg,0,ntg,0,npg)
00078 call alloc_array3d(work,0,nrg,0,ntg,0,npg)
00079 call alloc_array2d(sou_bhsurf,0,ntg,0,npg)
00080 call alloc_array2d(dsou_bhsurf,0,ntg,0,npg)
00081 call alloc_array2d(sou_outsurf,0,ntg,0,npg)
00082 call alloc_array2d(dsou_outsurf,0,ntg,0,npg)
00083
00084 iter_count = 0
00085 do
00086 iter_count = iter_count + 1
00087 write(6,'(a14,i5)') ' Iteration # =', iter_count
00088 count = dble(iter_count)
00089 conv_gra = dmin1(conv0_gra,conv_ini*count)
00090 conv_den = dmin1(conv0_den,conv_ini*count)
00091
00092 call calc_vector_x_grav(1)
00093 call calc_vector_phi_grav(1)
00094 call calc_vector_bh(1)
00095 call calc_shift_down2up
00096 call calc_shift2rotshift
00097
00098
00099 if (chgra == 'i') call cleargeometry
00100 if (chgra /= 'i') then
00101 call cristoffel_midpoint
00102 call cristoffel_gridpoint
00103 call transverse_part
00104
00105 call riccitensor_midpoint
00106 end if
00107 if (chgra == 'h'.or.chgra == 'c'.or.chgra == 'C' &
00108 &.or.chgra == 'H'.or.chgra == 'W') then
00109 call liegmab_midpoint
00110 call liegmab_gridpoint
00111 end if
00112
00113 call excurve_WL_midpoint
00114 call excurve_WL_gridpoint
00115
00116 call excurve_traceK_WL
00117 call ap2alps
00118
00119
00120
00121
00122
00123 call source_HaC_WL_BH(sou)
00124 call bh_boundary_analytic('psi ',0,sou_bhsurf,dsou_bhsurf)
00125 call bh_boundary_analytic('psi ',nrg,sou_outsurf,dsou_outsurf)
00126 call poisson_solver_1bh_homosol('dd',sou, &
00127 & sou_bhsurf,dsou_bhsurf, &
00128 & sou_outsurf,dsou_outsurf,pot)
00129 pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00130
00131
00132
00133
00134
00135
00136
00137 call update_grfield(pot,psi,conv_gra)
00138 call error_metric_type3(' psi',psi,pot_bak,error_all(1),flag_all(1),'bh')
00139
00140
00141
00142 call source_trG_WL_BH(sou)
00143 call bh_boundary_analytic('alps',0,sou_bhsurf,dsou_bhsurf)
00144 call bh_boundary_analytic('alps',nrg,sou_outsurf,dsou_outsurf)
00145 call poisson_solver_1bh_homosol('dd',sou, &
00146 & sou_bhsurf,dsou_bhsurf, &
00147 & sou_outsurf,dsou_outsurf,pot)
00148 pot_bak(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)
00149
00150
00151
00152
00153
00154
00155
00156 call compute_alps2alph(pot,psi)
00157 call update_grfield(pot,alph,conv_gra)
00158 call error_metric_type3('alph',alph,pot_bak,error_all(2),flag_all(2),'bh')
00159
00160
00161
00162 call source_MoC_WL_BH(souvec)
00163 sou1(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 1)
00164 sou2(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 2)
00165 sou3(1:nrg, 1:ntg, 1:npg) = souvec(1:nrg, 1:ntg, 1:npg, 3)
00166 call bh_boundary_analytic('bvxd',0,sou_bhsurf,dsou_bhsurf)
00167 call bh_boundary_analytic('bvxd',nrg,sou_outsurf,dsou_outsurf)
00168 call poisson_solver_1bh_homosol('dd',sou1, &
00169 & sou_bhsurf,dsou_bhsurf, &
00170 & sou_outsurf,dsou_outsurf,potx)
00171 call bh_boundary_analytic('bvyd',0,sou_bhsurf,dsou_bhsurf)
00172 call bh_boundary_analytic('bvyd',nrg,sou_outsurf,dsou_outsurf)
00173 call poisson_solver_1bh_homosol('dd',sou2, &
00174 & sou_bhsurf,dsou_bhsurf, &
00175 & sou_outsurf,dsou_outsurf,poty)
00176 call bh_boundary_analytic('bvzd',0,sou_bhsurf,dsou_bhsurf)
00177 call bh_boundary_analytic('bvzd',nrg,sou_outsurf,dsou_outsurf)
00178 call poisson_solver_1bh_homosol('dd',sou3, &
00179 & sou_bhsurf,dsou_bhsurf, &
00180 & sou_outsurf,dsou_outsurf,potz)
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 potx_bak(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00193 poty_bak(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00194 potz_bak(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00195 bvxd(0:nrg,0:ntg,0:npg) = potx(0:nrg,0:ntg,0:npg)
00196 bvyd(0:nrg,0:ntg,0:npg) = poty(0:nrg,0:ntg,0:npg)
00197 bvzd(0:nrg,0:ntg,0:npg) = potz(0:nrg,0:ntg,0:npg)
00198
00199 call calc_shift_down2up
00200 call gauge_shift
00201 call calc_shift_up2down
00202
00203
00204 potx(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00205 poty(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00206 potz(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00207 bvxd(0:nrg,0:ntg,0:npg) = potx_bak(0:nrg,0:ntg,0:npg)
00208 bvyd(0:nrg,0:ntg,0:npg) = poty_bak(0:nrg,0:ntg,0:npg)
00209 bvzd(0:nrg,0:ntg,0:npg) = potz_bak(0:nrg,0:ntg,0:npg)
00210
00211
00212
00213
00214 call update_grfield(potx,bvxd,conv_gra)
00215 call update_grfield(poty,bvyd,conv_gra)
00216 call update_grfield(potz,bvzd,conv_gra)
00217 call error_metric_type3('bvxd',bvxd,potx_bak,error_all(3),flag_all(3),'bh')
00218 call error_metric_type3('bvyd',bvyd,poty_bak,error_all(4),flag_all(4),'bh')
00219 call error_metric_type3('bvzd',bvzd,potz_bak,error_all(5),flag_all(5),'bh')
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 call source_trfreeG_WL_BH(souten)
00230 itg = ntg/2-2; ipg = 2
00231 open(16,file='test_vec_sou',status='unknown')
00232 do irg = 1, nrg
00233 write(16,'(1p,12e20.12)')hrg(irg), souvec(irg,itg,ipg,1) &
00234 & , souvec(irg,itg,ipg,2) &
00235 & , souvec(irg,itg,ipg,3) &
00236 & , souten(irg,itg,ipg,1) &
00237 & , souten(irg,itg,ipg,2) &
00238 & , souten(irg,itg,ipg,3) &
00239 & , souten(irg,itg,ipg,4) &
00240 & , souten(irg,itg,ipg,5) &
00241 & , souten(irg,itg,ipg,6)
00242 end do
00243 close(16)
00244
00245
00246
00247 if (chope == 'L') then
00248
00249 sou1(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,1)
00250 sou2(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,2)
00251 sou3(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,3)
00252 sou4(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,4)
00253 sou5(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,5)
00254 sou6(0:nrg,0:ntg,0:npg) = souten(0:nrg,0:ntg,0:npg,6)
00255
00256 call bh_boundary_analytic('hxxd',0,sou_bhsurf,dsou_bhsurf)
00257 call bh_boundary_analytic('hxxd',nrg,sou_outsurf,dsou_outsurf)
00258 call poisson_solver_1bh_homosol('dd',sou1, &
00259 & sou_bhsurf,dsou_bhsurf, &
00260 & sou_outsurf,dsou_outsurf,pot1)
00261 call bh_boundary_analytic('hxyd',0,sou_bhsurf,dsou_bhsurf)
00262 call bh_boundary_analytic('hxyd',nrg,sou_outsurf,dsou_outsurf)
00263 call poisson_solver_1bh_homosol('dd',sou2, &
00264 & sou_bhsurf,dsou_bhsurf, &
00265 & sou_outsurf,dsou_outsurf,pot2)
00266 call bh_boundary_analytic('hxzd',0,sou_bhsurf,dsou_bhsurf)
00267 call bh_boundary_analytic('hxzd',nrg,sou_outsurf,dsou_outsurf)
00268 call poisson_solver_1bh_homosol('dd',sou3, &
00269 & sou_bhsurf,dsou_bhsurf, &
00270 & sou_outsurf,dsou_outsurf,pot3)
00271 call bh_boundary_analytic('hyyd',0,sou_bhsurf,dsou_bhsurf)
00272 call bh_boundary_analytic('hyyd',nrg,sou_outsurf,dsou_outsurf)
00273 call poisson_solver_1bh_homosol('dd',sou4, &
00274 & sou_bhsurf,dsou_bhsurf, &
00275 & sou_outsurf,dsou_outsurf,pot4)
00276 call bh_boundary_analytic('hyzd',0,sou_bhsurf,dsou_bhsurf)
00277 call bh_boundary_analytic('hyzd',nrg,sou_outsurf,dsou_outsurf)
00278 call poisson_solver_1bh_homosol('dd',sou5, &
00279 & sou_bhsurf,dsou_bhsurf, &
00280 & sou_outsurf,dsou_outsurf,pot5)
00281 call bh_boundary_analytic('hzzd',0,sou_bhsurf,dsou_bhsurf)
00282 call bh_boundary_analytic('hzzd',nrg,sou_outsurf,dsou_outsurf)
00283 call poisson_solver_1bh_homosol('dd',sou6, &
00284 & sou_bhsurf,dsou_bhsurf, &
00285 & sou_outsurf,dsou_outsurf,pot6)
00286
00287 end if
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310 pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00311 bk1(0:nrg,0:ntg,0:npg) = hxxd(0:nrg,0:ntg,0:npg)
00312 bk2(0:nrg,0:ntg,0:npg) = hxyd(0:nrg,0:ntg,0:npg)
00313 bk3(0:nrg,0:ntg,0:npg) = hxzd(0:nrg,0:ntg,0:npg)
00314 bk4(0:nrg,0:ntg,0:npg) = hyyd(0:nrg,0:ntg,0:npg)
00315 bk5(0:nrg,0:ntg,0:npg) = hyzd(0:nrg,0:ntg,0:npg)
00316 bk6(0:nrg,0:ntg,0:npg) = hzzd(0:nrg,0:ntg,0:npg)
00317
00318 hxxd(0:nrg,0:ntg,0:npg) = pot1(0:nrg,0:ntg,0:npg)
00319 hxyd(0:nrg,0:ntg,0:npg) = pot2(0:nrg,0:ntg,0:npg)
00320 hxzd(0:nrg,0:ntg,0:npg) = pot3(0:nrg,0:ntg,0:npg)
00321 hyyd(0:nrg,0:ntg,0:npg) = pot4(0:nrg,0:ntg,0:npg)
00322 hyzd(0:nrg,0:ntg,0:npg) = pot5(0:nrg,0:ntg,0:npg)
00323 hzzd(0:nrg,0:ntg,0:npg) = pot6(0:nrg,0:ntg,0:npg)
00324
00325
00326
00327
00328 call invhij
00329
00330 call gauge_hiju(conv_gra)
00331 call invhij_up2down
00332 call adjusthij
00333
00334
00335
00336 pot(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00337 pot1(0:nrg,0:ntg,0:npg) = hxxd(0:nrg,0:ntg,0:npg)
00338 pot2(0:nrg,0:ntg,0:npg) = hxyd(0:nrg,0:ntg,0:npg)
00339 pot3(0:nrg,0:ntg,0:npg) = hxzd(0:nrg,0:ntg,0:npg)
00340 pot4(0:nrg,0:ntg,0:npg) = hyyd(0:nrg,0:ntg,0:npg)
00341 pot5(0:nrg,0:ntg,0:npg) = hyzd(0:nrg,0:ntg,0:npg)
00342 pot6(0:nrg,0:ntg,0:npg) = hzzd(0:nrg,0:ntg,0:npg)
00343
00344 psi(0:nrg,0:ntg,0:npg) = pot_bak(0:nrg,0:ntg,0:npg)
00345 hxxd(0:nrg,0:ntg,0:npg) = bk1(0:nrg,0:ntg,0:npg)
00346 hxyd(0:nrg,0:ntg,0:npg) = bk2(0:nrg,0:ntg,0:npg)
00347 hxzd(0:nrg,0:ntg,0:npg) = bk3(0:nrg,0:ntg,0:npg)
00348 hyyd(0:nrg,0:ntg,0:npg) = bk4(0:nrg,0:ntg,0:npg)
00349 hyzd(0:nrg,0:ntg,0:npg) = bk5(0:nrg,0:ntg,0:npg)
00350 hzzd(0:nrg,0:ntg,0:npg) = bk6(0:nrg,0:ntg,0:npg)
00351
00352 call update_grfield(pot_bak,psi,conv_gra)
00353
00354 call update_grfield(pot1,hxxd,conv_gra)
00355 call update_grfield(pot2,hxyd,conv_gra)
00356 call update_grfield(pot3,hxzd,conv_gra)
00357 call update_grfield(pot4,hyyd,conv_gra)
00358 call update_grfield(pot5,hyzd,conv_gra)
00359 call update_grfield(pot6,hzzd,conv_gra)
00360
00361
00362
00363
00364
00365
00366
00367 call invhij
00368 call compute_fnc_multiple(alph,psi,alps)
00369
00370 call error_metric_type3('hxxd',hxxd,bk1,error_all(6),flag_all(6),'bh')
00371 call error_metric_type3('hxyd',hxyd,bk2,error_all(7),flag_all(7),'bh')
00372 call error_metric_type3('hxzd',hxzd,bk3,error_all(8),flag_all(8),'bh')
00373 call error_metric_type3('hyyd',hyyd,bk4,error_all(9),flag_all(9),'bh')
00374 call error_metric_type3('hyzd',hyzd,bk5,error_all(10),flag_all(10),'bh')
00375 call error_metric_type3('hzzd',hzzd,bk6,error_all(11),flag_all(11),'bh')
00376
00377
00378
00379 itg = ntg/2-2; ipg = npg/2-2
00380 open(16,file='test_vec2',status='unknown')
00381 do irg = 0, nrg
00382 write(16,'(1p,12e20.12)') rg(irg), psi(irg,itg,ipg) &
00383 & , alph(irg,itg,ipg) &
00384 & , bvxd(irg,itg,ipg) &
00385 & , bvyd(irg,itg,ipg) &
00386 & , bvzd(irg,itg,ipg) &
00387 & , hxxd(irg,itg,ipg) &
00388 & , hxyd(irg,itg,ipg) &
00389 & , hxzd(irg,itg,ipg) &
00390 & , hyyd(irg,itg,ipg) &
00391 & , hyzd(irg,itg,ipg) &
00392 & , hzzd(irg,itg,ipg)
00393 end do
00394 close(16)
00395
00396
00397
00398 if (flag_all(1)==0.and.flag_all(2)==0.and.flag_all(3)==0.and.&
00399 & flag_all(4)==0.and.flag_all(5)==0.and.flag_all(6)==0.and.&
00400 & flag_all(7)==0.and.flag_all(8)==0.and.flag_all(9)==0.and.&
00401 & flag_all(10)==0.and.flag_all(11)==0.and.flag_all(12)==0.and.&
00402 & flag_all(13)==0.and.flag_all(14)==0.and.flag_all(15)==0.and.&
00403 & flag_all(16)==0.and.flag_all(17)==0.and.flag_all(18)==0.and.&
00404 & flag_all(19)==0.and.flag_all(20)==0) exit
00405 if (iter_count >= iter_max) exit
00406 flag = 0
00407 flag_all(1:20) = 0
00408
00409 end do
00410
00411 deallocate(sou)
00412 deallocate(pot)
00413 deallocate(potx)
00414 deallocate(poty)
00415 deallocate(potz)
00416 deallocate(souvec)
00417 deallocate(souten)
00418 deallocate(sou1)
00419 deallocate(sou2)
00420 deallocate(sou3)
00421 deallocate(sou4)
00422 deallocate(sou5)
00423 deallocate(sou6)
00424 deallocate(pot1)
00425 deallocate(pot2)
00426 deallocate(pot3)
00427 deallocate(pot4)
00428 deallocate(pot5)
00429 deallocate(pot6)
00430 deallocate(pot_bak)
00431 deallocate(potx_bak)
00432 deallocate(poty_bak)
00433 deallocate(potz_bak)
00434 deallocate(bk1)
00435 deallocate(bk2)
00436 deallocate(bk3)
00437 deallocate(bk4)
00438 deallocate(bk5)
00439 deallocate(bk6)
00440 deallocate(work)
00441 deallocate( sou_bhsurf)
00442 deallocate(dsou_bhsurf)
00443 deallocate( sou_outsurf)
00444 deallocate(dsou_outsurf)
00445 end subroutine iteration_BH_WL
