00001 subroutine iteration_BBH_CF_AH(iter_count)
00002 use phys_constant, only : long
00003 use grid_parameter
00004 use coordinate_grav_r
00005 use coordinate_grav_phi
00006 use coordinate_grav_theta
00007 use weight_midpoint_grav
00008 use make_array_2d
00009 use make_array_3d
00010 use make_array_4d
00011 use def_metric
00012 use def_matter
00013 use def_vector_bh
00014 use interface_poisson_solver_binary_bhex_homosol
00015 use interface_sourceterm_exsurf_eqm_binary
00016 use interface_sourceterm_surface_int
00017 use interface_compute_shift
00018 use interface_compute_alps2alph
00019 use interface_update_grfield
00020 use interface_update_parameter
00021 use interface_error_metric_type0
00022 use interface_error_metric_type2
00023 use interface_interpolation_fillup_binary
00024
00025 use interface_bh_boundary_d_psi
00026 use interface_bh_boundary_nh_psi_test
00027 use interface_bh_boundary_d_alps
00028 use interface_bh_boundary_AH
00029 use interface_outer_boundary_d_psi
00030 use interface_outer_boundary_d_alps
00031 use interface_outer_boundary_d_bvxd
00032 use interface_outer_boundary_d_bvyd
00033 use interface_outer_boundary_d_bvzd
00034 use interface_outer_boundary_d_Bfun
00035 use interface_outer_boundary_n_Bfun
00036 use interface_outer_boundary_d_potx
00037 use interface_outer_boundary_d_poty
00038 use interface_outer_boundary_d_potz
00039 use interface_sourceterm_HaC_CF
00040 use interface_sourceterm_trG_CF
00041 use interface_sourceterm_MoC_CF
00042 use interface_sourceterm_MoC_CF_type1_bhex
00043 use interface_sourceterm_Bfun
00044 use interface_interpolation_fillup_binary_parity
00045 use interface_sourceterm_exsurf_eqm_binary_parity
00046 use interface_compute_shift_v2
00047 use interface_compute_dBfun
00048 use gnufor2
00049
00050 implicit none
00051 real(long), pointer :: sou(:,:,:), pot(:,:,:)
00052 real(long), pointer :: Bfun(:,:,:), potx(:,:,:), poty(:,:,:), potz(:,:,:)
00053 real(long), pointer :: pot_bak(:,:,:), work(:,:,:)
00054 real(long), pointer :: potx_bak(:,:,:), poty_bak(:,:,:), potz_bak(:,:,:)
00055 real(long), pointer :: dBfundx(:,:,:), dBfundy(:,:,:), dBfundz(:,:,:)
00056 real(long), pointer :: souvec(:,:,:,:)
00057 real(long), pointer :: sou_exsurf(:,:), dsou_exsurf(:,:)
00058 real(long), pointer :: sou_bhsurf(:,:), dsou_bhsurf(:,:)
00059 real(long), pointer :: sou_outsurf(:,:), dsou_outsurf(:,:)
00060 real(long) :: count, error_psi,error_alph, error_bvxd, error_bvyd, error_bvzd
00061 real(long) :: error_potx, error_poty, error_potz, error_Bfun
00062 real(long) :: pari
00063 integer :: iter_count, flag_psi=0, flag_alph=0, flag_bvxd=0, flag_bvyd=0, flag_bvzd=0
00064 integer :: flag_potx=0, flag_poty=0, flag_potz=0, flag_Bfun=0
00065 integer :: irg, itg, ipg, ihy, nn, ii
00066 character(30) :: char1, char2, char3
00067
00068 nn = 100
00069 call alloc_array4d(souvec,0,nrg,0,ntg,0,npg,1,3)
00070 call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00071
00072
00073
00074
00075 call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00076 call alloc_array3d(potx,0,nrg,0,ntg,0,npg)
00077 call alloc_array3d(poty,0,nrg,0,ntg,0,npg)
00078 call alloc_array3d(potz,0,nrg,0,ntg,0,npg)
00079 call alloc_array3d(pot_bak,0,nrg,0,ntg,0,npg)
00080 call alloc_array3d(potx_bak,0,nrg,0,ntg,0,npg)
00081 call alloc_array3d(poty_bak,0,nrg,0,ntg,0,npg)
00082 call alloc_array3d(potz_bak,0,nrg,0,ntg,0,npg)
00083 call alloc_array3d(work,0,nrg,0,ntg,0,npg)
00084 call alloc_array2d(sou_exsurf,0,ntg,0,npg)
00085 call alloc_array2d(dsou_exsurf,0,ntg,0,npg)
00086 call alloc_array2d(sou_bhsurf,0,ntg,0,npg)
00087 call alloc_array2d(dsou_bhsurf,0,ntg,0,npg)
00088 call alloc_array2d(sou_outsurf,0,ntg,0,npg)
00089 call alloc_array2d(dsou_outsurf,0,ntg,0,npg)
00090
00091 iter_count = 0
00092 psi(0:nrg,0:ntg,0:npg) =1.0d0
00093 alph(0:nrg,0:ntg,0:npg)=1.0d0
00094 alps(0:nrg,0:ntg,0:npg)=1.0d0
00095 bvxd(0:nrg,0:ntg,0:npg)=0.0d0
00096 bvyd(0:nrg,0:ntg,0:npg)=0.0d0
00097 bvzd(0:nrg,0:ntg,0:npg)=0.0d0
00098
00099
00100 write(6 ,'(a12,6a15)') 'Errors:','psi ','alph','bvxd','bvyd','bvzd'
00101
00102
00103
00104
00105
00106 do
00107 iter_count = iter_count + 1
00108 count = dble(iter_count)
00109 conv_gra = dmin1(conv0_gra,conv_ini*count)
00110 conv_den = dmin1(conv0_den,conv_ini*count)
00111
00112
00113
00114
00115
00116 call calc_vector_bh(2)
00117 call excurve
00118 call excurve_CF_gridpoint_bhex
00119
00120 sou(0:nrg,0:ntg,0:npg) = 0.0d0
00121 call sourceterm_HaC_CF(sou)
00122 call sourceterm_exsurf_eqm_binary(psi,sou_exsurf,dsou_exsurf)
00123 call sourceterm_surface_int(psi,0,sou_bhsurf,dsou_bhsurf)
00124 call sourceterm_surface_int(psi,nrg,sou_outsurf,dsou_outsurf)
00125 call bh_boundary_AH(sou_bhsurf,dsou_bhsurf,'psi ')
00126 call outer_boundary_d_psi(sou_outsurf)
00127 call poisson_solver_binary_bhex_homosol('dd',sou, &
00128 & sou_exsurf,dsou_exsurf, &
00129 & sou_bhsurf,dsou_bhsurf, &
00130 & sou_outsurf,dsou_outsurf,pot)
00131 pot_bak = psi
00132 call update_grfield(pot,psi,conv_gra)
00133 call interpolation_fillup_binary(psi)
00134 call error_metric_type0(psi,pot_bak,error_psi,flag_psi,'bh')
00135
00136 sou(0:nrg,0:ntg,0:npg) = 0.0d0
00137 alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00138 call sourceterm_trG_CF(sou)
00139 call sourceterm_exsurf_eqm_binary(alps,sou_exsurf,dsou_exsurf)
00140 call sourceterm_surface_int(alps,0,sou_bhsurf,dsou_bhsurf)
00141 call sourceterm_surface_int(alps,nrg,sou_outsurf,dsou_outsurf)
00142 call bh_boundary_AH(sou_bhsurf,dsou_bhsurf,'alps')
00143 call outer_boundary_d_alps(sou_outsurf)
00144 call poisson_solver_binary_bhex_homosol('dd',sou, &
00145 & sou_exsurf,dsou_exsurf, &
00146 & sou_bhsurf,dsou_bhsurf, &
00147 & sou_outsurf,dsou_outsurf,pot)
00148 call compute_alps2alph(pot,psi)
00149 pot_bak = alph
00150 call update_grfield(pot,alph,conv_gra)
00151 call interpolation_fillup_binary(alph)
00152 call error_metric_type0(alph,pot_bak,error_alph,flag_alph,'bh')
00153
00154 call sourceterm_MoC_CF_type1_bhex(souvec)
00155 do ii = 1, 3
00156 sou(0:nrg,0:ntg,0:npg) = souvec(0:nrg,0:ntg,0:npg,ii)
00157 if (ii.eq.1) then
00158 work(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00159 pari = -1.0d0
00160 else if (ii.eq.2) then
00161 work(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00162 pari = -1.0d0
00163 else if (ii.eq.3) then
00164 work(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00165 pari = 1.0d0
00166 end if
00167 call sourceterm_exsurf_eqm_binary_parity(work,sou_exsurf,dsou_exsurf,&
00168 & pari)
00169 call sourceterm_surface_int(work,0,sou_bhsurf,dsou_bhsurf)
00170 call sourceterm_surface_int(work,nrg,sou_outsurf,dsou_outsurf)
00171 if (ii.eq.1) then
00172 call bh_boundary_AH(sou_bhsurf,dsou_bhsurf,'bvxd')
00173 call outer_boundary_d_bvxd(sou_outsurf)
00174 else if (ii.eq.2) then
00175 call bh_boundary_AH(sou_bhsurf,dsou_bhsurf,'bvyd')
00176 call outer_boundary_d_bvyd(sou_outsurf)
00177 else if (ii.eq.3) then
00178 call bh_boundary_AH(sou_bhsurf,dsou_bhsurf,'bvzd')
00179 call outer_boundary_d_bvzd(sou_outsurf)
00180 end if
00181 call poisson_solver_binary_bhex_homosol('dd',sou, &
00182 & sou_exsurf,dsou_exsurf, &
00183 & sou_bhsurf,dsou_bhsurf, &
00184 & sou_outsurf,dsou_outsurf,pot)
00185 if (ii.eq.1) potx(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00186 if (ii.eq.2) poty(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00187 if (ii.eq.3) potz(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00188 end do
00189
00190 potx_bak = bvxd
00191 poty_bak = bvyd
00192 potz_bak = bvzd
00193 call update_grfield(potx,bvxd,conv_gra)
00194 call update_grfield(poty,bvyd,conv_gra)
00195 call update_grfield(potz,bvzd,conv_gra)
00196
00197 call interpolation_fillup_binary_parity(bvxd,-1.0d0)
00198 call interpolation_fillup_binary_parity(bvyd,-1.0d0)
00199 call interpolation_fillup_binary_parity(bvzd,1.0d0)
00200
00201 call error_metric_type2(bvxd,potx_bak,error_bvxd,flag_bvxd,'bh')
00202 call error_metric_type2(bvyd,poty_bak,error_bvyd,flag_bvyd,'bh')
00203 call error_metric_type2(bvzd,potz_bak,error_bvzd,flag_bvzd,'bh')
00204
00205 call printout_error_metric_combined(iter_count,error_psi,error_alph,&
00206 & error_bvxd,error_bvyd,error_bvzd)
00207
00208 if (mod(iter_count,5)==0.or.iter_count==1) then
00209 write(char1, '(i5)') iter_count
00210 char2 = adjustl(char1)
00211 char3 = 'iteration' // trim(char2) // '.txt'
00212
00213 open(15,file='./'//char3,status='unknown')
00214 do irg = nrg, 0, -1
00215 write(15,'(1p,9e20.12)') -rg(irg), psi(irg,ntgeq,npgxzm) &
00216 & , alph(irg,ntgeq,npgxzm) &
00217 & , bvxd(irg,ntgeq,npgxzm) &
00218 & , bvyd(irg,ntgeq,npgxzm) &
00219 & , bvzd(irg,ntgeq,npgxzm) &
00220 & , bvxd(irg,ntg*5/8,npg*9/16) &
00221 & , bvyd(irg,ntg*5/8,npg*9/16) &
00222 & , bvzd(irg,ntg*5/8,npg*9/16)
00223 end do
00224 write(15,*) ' '
00225 do irg = 0, nrg
00226 write(15,'(1p,9e20.12)') rg(irg), psi(irg,ntgeq,npgxzp) &
00227 & , alph(irg,ntgeq,npgxzp) &
00228 & , bvxd(irg,ntgeq,npgxzp) &
00229 & , bvyd(irg,ntgeq,npgxzp) &
00230 & , bvzd(irg,ntgeq,npgxzp) &
00231 & , bvxd(irg,ntg*3/8,npg/16) &
00232 & , bvyd(irg,ntg*3/8,npg/16) &
00233 & , bvzd(irg,ntg*3/8,npg/16)
00234 end do
00235 close(15)
00236 endif
00237
00238 if ((flag_psi==0).and.(flag_alph==0).and.(flag_bvxd==0).and.(flag_bvyd==0).and.(flag_bvzd==0)) exit
00239 if (iter_count >= iter_max) exit
00240
00241
00242
00243
00244
00245 flag_psi = 0
00246 flag_alph = 0
00247 flag_bvxd = 0
00248 flag_bvyd = 0
00249 flag_bvzd = 0
00250 flag_Bfun = 0
00251 end do
00252
00253
00254
00255
00256
00257
00258
00259 deallocate(souvec)
00260 deallocate(sou)
00261
00262
00263
00264
00265
00266
00267
00268 deallocate(pot)
00269 deallocate(potx)
00270 deallocate(poty)
00271 deallocate(potz)
00272 deallocate(pot_bak)
00273 deallocate(potx_bak)
00274 deallocate(poty_bak)
00275 deallocate(potz_bak)
00276 deallocate(work)
00277 deallocate(sou_exsurf)
00278 deallocate(dsou_exsurf)
00279 deallocate(sou_bhsurf)
00280 deallocate(dsou_bhsurf)
00281 deallocate(sou_outsurf)
00282 deallocate(dsou_outsurf)
00283 end subroutine iteration_BBH_CF_AH
