00001 subroutine iter_lecc_spin_BNS_CF_mpt(iter_count, iseq, iter_eps)
00002 use phys_constant, only : long, nmpt
00003 use grid_parameter
00004 use grid_parameter_interpo, only : nrg_itp, ntg_itp, npg_itp
00005 use def_metric, only : psi, alph, alps, bvxd, bvyd, bvzd
00006 use coordinate_grav_r
00007 use grid_points_binary_excision, only : irg_exin, irg_exout
00008
00009 use def_matter
00010 use def_velocity_rot
00011 use def_velocity_potential
00012 use def_matter_parameter
00013 use interface_interpo_fl2gr
00014 use interface_sourceterm_HaC_peos_spin
00015 use interface_sourceterm_HaC_CF
00016 use interface_sourceterm_trG_peos_spin
00017 use interface_sourceterm_trG_CF
00018 use interface_sourceterm_MoC_CF_with_divshift_peos_spin
00019 use interface_sourceterm_MoC_CF_with_divshift_r3rd
00020
00021
00022
00023
00024 use interface_hydrostatic_eq_CF_peos_spin
00025 use interface_calc_gradvep
00026 use interface_hydro_spbns_vep_CF_peos
00027 use interface_calc_surface
00028 use interface_update_matter
00029 use interface_update_parameter_BNS_spin
00030 use interface_update_surface
00031 use interface_error_matter
00032 use interface_error_matter_type2
00033 use interface_error_velocity_potential
00034
00035
00036
00037 use def_binary_parameter, only : dis
00038 use copy_array_4dto3d_mpt
00039 use make_array_2d
00040 use make_array_3d
00041 use make_array_4d
00042 use interface_compute_alps2alph
00043 use interface_sourceterm_surface_int
00044 use interface_sourceterm_exsurf_binary_COCP
00045 use interface_sourceterm_outsurf_COCP_from_ARCP
00046 use interface_sourceterm_insurf_ARCP_from_COCP
00047 use interface_outer_boundary_d_BBH_CF
00048 use interface_outer_boundary_d_psi
00049 use interface_outer_boundary_d_alps
00050 use interface_outer_boundary_d_bvxd
00051 use interface_outer_boundary_d_bvyd
00052 use interface_outer_boundary_d_bvzd
00053
00054 use interface_poisson_solver_binary_star_homosol
00055 use interface_poisson_solver_asymptotic_patch_homosol
00056 use interface_interpolation_fillup_binary_COCP
00057 use interface_update_grfield
00058 use interface_error_metric
00059 use interface_error_metric_type0
00060 use interface_error_metric_type2
00061 use interface_error_metric_type2_mpt
00062 implicit none
00063 real(long), pointer :: sou(:,:,:), pot(:,:,:)
00064 real(long), pointer :: potx(:,:,:), poty(:,:,:), potz(:,:,:)
00065 real(long), pointer :: pot_bak(:,:,:), work(:,:,:)
00066 real(long), pointer :: potx_bak(:,:,:), poty_bak(:,:,:), potz_bak(:,:,:)
00067 real(long), pointer :: souvec(:,:,:,:), MoC_vio(:,:,:,:)
00068 real(long), pointer :: sou_exsurf(:,:), dsou_exsurf(:,:)
00069
00070 real(long), pointer :: sou_insurf(:,:), dsou_insurf(:,:)
00071 real(long), pointer :: sou_outsurf(:,:), dsou_outsurf(:,:)
00072 real(long), pointer :: fnc_itp(:,:,:)
00073 real(long), pointer :: potf(:,:,:), emd_bak(:,:,:), vep_bak(:,:,:)
00074 real(long), pointer :: potrs(:,:)
00075 real(long), pointer :: potut(:,:,:), potux(:,:,:), potuy(:,:,:), potuz(:,:,:)
00076
00077 real(long) :: count, iter_eps
00078 real(long) :: error_psi = 0.0d0, error_alph = 0.0d0, error_tmp = 0.0d0,
00079 error_bvxd = 0.0d0, error_bvyd = 0.0d0, error_bvzd = 0.0d0, &
00080 error_ome = 2.0d0, error_emd = 0.0d0, error_metr = 0.0d0, &
00081 error_vep = 0.0d0, error_emd_last, error_psi_last, error_alph_last, &
00082 error_vep_last
00083 real(long) :: pari, ome_bh_new, work_shift, error_metric_mpt_mid(10,10)
00084 integer :: iter_count, iter_extra = 0, istep_niq = -1,
00085 flag_tmp = 0, flag_param = 99, flag_emd = 0, &
00086 flag_psi = 0, flag_alph = 0, flag_vep = 0, &
00087 flag_bvxd = 0, flag_bvyd = 0, flag_bvzd = 0
00088 integer :: irg, itg, ipg, ii, iseq, jj, sumjj, je
00089 integer :: impt, impt_ex, impt_bin
00090 character(len=2) :: chgreen, chpa, chpB, cocp
00091 character(len=4) :: chbe
00092 integer :: flag = 0, hydro_iter = 4, ishift, shift_iter=1,
00093 flag_mid, flag_mpt_mid(10)
00094 integer :: irf, itf, ipf, ihy, npg_l, npg_r, icount
00095 character(30) :: char1, char2, char3, char4, char5
00096 character(len=1) :: char_1
00097
00098 call set_allocate_size_mpt
00099 call alloc_array4d(souvec ,0,nrg,0,ntg,0,npg,1,3)
00100 call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00101 call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00102 call alloc_array3d(potx,0,nrg,0,ntg,0,npg)
00103 call alloc_array3d(poty,0,nrg,0,ntg,0,npg)
00104 call alloc_array3d(potz,0,nrg,0,ntg,0,npg)
00105 call alloc_array3d(pot_bak,0,nrg,0,ntg,0,npg)
00106 call alloc_array3d(potx_bak,0,nrg,0,ntg,0,npg)
00107 call alloc_array3d(poty_bak,0,nrg,0,ntg,0,npg)
00108 call alloc_array3d(potz_bak,0,nrg,0,ntg,0,npg)
00109 call alloc_array3d(work,0,nrg,0,ntg,0,npg)
00110 call alloc_array3d(fnc_itp,0,nrg,0,ntg,0,npg)
00111 call alloc_array2d(sou_exsurf,0,ntg,0,npg)
00112 call alloc_array2d(dsou_exsurf,0,ntg,0,npg)
00113 call alloc_array2d(sou_insurf,0,ntg,0,npg)
00114 call alloc_array2d(dsou_insurf,0,ntg,0,npg)
00115 call alloc_array2d(sou_outsurf,0,ntg,0,npg)
00116 call alloc_array2d(dsou_outsurf,0,ntg,0,npg)
00117 call alloc_array3d( potf,0,nrf,0,ntf,0,npf)
00118 call alloc_array3d(potut,0,nrf,0,ntf,0,npf)
00119 call alloc_array3d(potux,0,nrf,0,ntf,0,npf)
00120 call alloc_array3d(potuy,0,nrf,0,ntf,0,npf)
00121 call alloc_array3d(potuz,0,nrf,0,ntf,0,npf)
00122 call alloc_array3d(emd_bak,0,nrf,0,ntf,0,npf)
00123 call alloc_array3d(vep_bak,0,nrf,0,ntf,0,npf)
00124 call alloc_array2d(potrs,0,ntf,0,npf)
00125
00126 iter_count = 0
00127 error_emd_last =1.0d0
00128 error_psi_last =1.0d0
00129 error_alph_last=1.0d0
00130 flag_mid=0
00131
00132 do
00133 iter_count = iter_count + 1
00134
00135 count = dble(iter_count)
00136 flag_mpt_mid(1:10) =0
00137 write(6,'(a10,i4,a108)') "Iteration:",iter_count, &
00138 & " ---psi------------alph------------bvxd------------bvyd------------bvzd------------emd------------vep-------"
00139
00140 do impt = 1, nmpt
00141 write(6,'(a8,i2,a1)', ADVANCE = "NO") " Patch #", impt, ""
00142 write(6,*) " "
00143 call copy_grid_parameter_from_mpt(impt)
00144 conv_gra = dmin1(conv0_gra,conv_ini*count)
00145 conv_den = dmin1(conv0_den,conv_ini*count)
00146 conv_vep = dmin1(conv0_vep,conv_ini*count)
00147 call copy_grid_parameter_to_mpt(impt)
00148
00149 call copy_from_mpatch_all_BNS_CF(impt)
00150 call copy_def_metric_and_matter_from_mpt(impt)
00151 call copy_def_matter_spin_from_mpt(impt)
00152
00153
00154
00155
00156
00157 if (impt==1 .or. impt==2) then
00158 cocp = 'ns'
00159 call interpo_fl2gr(emd, emdg)
00160 call interpo_fl2gr(vepxf, vepxg)
00161 call interpo_fl2gr(vepyf, vepyg)
00162 call interpo_fl2gr(vepzf, vepzg)
00163 call interpo_fl2gr(wxspf, wxspg)
00164 call interpo_fl2gr(wyspf, wyspg)
00165 call interpo_fl2gr(wzspf, wzspg)
00166
00167 if (mod(iter_count,5)==0.and.impt==10) then
00168 write(char1, '(i5)') iter_count
00169 char2 = adjustl(char1)
00170 write(char4, '(i5)') impt
00171 char5 = adjustl(char4)
00172
00173 char3 = 'vep_iteration' // trim(char2) // '_mpt' // trim(char5) // '.txt'
00174 open(12,file=char3,status='unknown')
00175 itg=12; ipg=6
00176 write(12,'(a1,a2,4a16,a30)') '#', 'ir', 'rg*rs', 'vepxf', 'vepyf', 'vepzf', 'plot using every :::0::0'
00177 do irg=0,nrf
00178 write(12,'(i3,1p,4e16.6)') irg, rg(irg)*rs(itg,ipg), vepxf(irg,itg,ipg), vepyf(irg,itg,ipg), vepzf(irg,itg,ipg)
00179 end do
00180 write(12,*) '#------------------------------------------------------------------------------------------------------'
00181 write(12,*) ""
00182 write(12,'(a1,a2,4a16,a30)') '#', 'ir', 'rg', 'vepxg', 'vepyg', 'vepzg', 'plot using every :::1::1'
00183 do irg=0,nrg
00184 write(12,'(i3,1p,4e16.6)') irg, rg(irg), vepxg(irg,itg,ipg), vepyg(irg,itg,ipg), vepzg(irg,itg,ipg)
00185 end do
00186 write(12,*) ""
00187 close(12)
00188 endif
00189
00190
00191 call calc_vector_x_matter(2)
00192 call calc_vector_phi_matter(2)
00193 call excurve_CF(cocp)
00194 call excurve_CF_gridpoint
00195 call copy_def_metric_and_matter_to_mpt(impt)
00196 call copy_def_matter_spin_to_mpt(impt)
00197 if (iter_count==1 .and. impt==1) then
00198 open(2,file='iter_data.txt',status='unknown')
00199 end if
00200 write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00201 else
00202 cocp = 'bh'
00203 call excurve_CF(cocp)
00204 call excurve_CF_gridpoint_bhex
00205 call copy_def_metric_and_matter_to_mpt(impt)
00206 call copy_def_matter_spin_to_mpt(impt)
00207 write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00208 write(2,*) '-----------------------------------------------------------------------'
00209 end if
00210
00211 if (impt.eq.nmpt) call calc_mass_asympto('ns')
00212 call copy_to_mpatch_all_BNS_CF(impt)
00213
00214
00215 if(error_psi_last.le.1.0d-13) then
00216 write(6,'(a16)', ADVANCE = "NO") ' ## SKIP psi ##'
00217 write(6,*) " "
00218 else
00219 sou(0:nrg,0:ntg,0:npg) = 0.0d0
00220 sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00221 sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00222 if (impt==1 .or. impt==2) then
00223 call sourceterm_HaC_peos_spin(sou)
00224 else
00225 call sourceterm_HaC_CF(sou)
00226 end if
00227
00228 if (impt==1 .or. impt==2) then
00229 call sourceterm_exsurf_binary_COCP(impt,'psi ', 'ev', sou_exsurf, dsou_exsurf)
00230 call sourceterm_outsurf_COCP_from_ARCP(impt,'psi ', 'ev', sou_outsurf, dsou_outsurf)
00231 call poisson_solver_binary_star_homosol('sd', sou, sou_exsurf, dsou_exsurf, &
00232 & sou_outsurf,dsou_outsurf,pot)
00233 else if (impt.eq.nmpt) then
00234 call sourceterm_insurf_ARCP_from_COCP(impt,'psi ', 'ev', sou_insurf, dsou_insurf)
00235 call outer_boundary_d_psi(sou_outsurf)
00236
00237 call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf, dsou_insurf, &
00238 & sou_outsurf,dsou_outsurf,pot)
00239 end if
00240
00241 pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00242 call update_grfield(pot,psi,conv_gra)
00243
00244 if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00245 & call update_parameter_BNS_spin(conv_gra)
00246
00247 if (impt.eq.nmpt) psi(nrg,0:ntg,0:npg) = 1.0d0
00248
00249 if (impt.ne.nmpt) call interpolation_fillup_binary_COCP(impt,'psi ','ev',psi)
00250
00251 call error_metric_type2_mpt(psi,pot_bak,error_tmp,flag_tmp,cocp,impt)
00252 flag_psi = max ( flag_psi, flag_tmp)
00253 error_psi = dmax1(error_psi,error_tmp)
00254
00255 error_psi_last = error_psi
00256
00257 call copy_to_mpatch_all_BNS_CF(impt)
00258 call copy_def_metric_and_matter_to_mpt(impt)
00259 call copy_def_matter_spin_to_mpt(impt)
00260 end if
00261
00262
00263 if(error_alph_last.le.1.0d-13) then
00264 write(6,'(a16)', ADVANCE = "NO") ' ## SKIP alph ##'
00265 write(6,*) " "
00266 else
00267 sou(0:nrg,0:ntg,0:npg) = 0.0d0
00268 sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00269 sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00270 alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00271 if (impt==1 .or. impt==2) then
00272 call sourceterm_trG_peos_spin(sou)
00273 else
00274 call sourceterm_trG_CF(sou)
00275 end if
00276 if (impt==1 .or. impt==2) then
00277 call sourceterm_exsurf_binary_COCP(impt,'alps','ev',sou_exsurf, dsou_exsurf)
00278 call sourceterm_outsurf_COCP_from_ARCP(impt,'alps','ev',sou_outsurf, dsou_outsurf)
00279 call poisson_solver_binary_star_homosol('sd',sou, sou_exsurf, dsou_exsurf, &
00280 & sou_outsurf,dsou_outsurf,pot)
00281 else if (impt.eq.nmpt) then
00282 call sourceterm_insurf_ARCP_from_COCP(impt,'alps','ev',sou_insurf, dsou_insurf)
00283 call outer_boundary_d_alps(sou_outsurf)
00284
00285 call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf,dsou_insurf, &
00286 & sou_outsurf,dsou_outsurf,pot)
00287 end if
00288 call compute_alps2alph(pot,psi)
00289 pot_bak(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)
00290 call update_grfield(pot,alph,conv_gra)
00291
00292 if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00293 & call update_parameter_BNS_spin(conv_gra)
00294
00295 if (impt.eq.nmpt) alph(nrg,0:ntg,0:npg) = 1.0d0
00296
00297 if (impt.ne.nmpt) call interpolation_fillup_binary_COCP(impt,'alph','ev',alph)
00298
00299 call error_metric_type2_mpt(alph,pot_bak,error_tmp,flag_tmp,cocp,impt)
00300 flag_alph = max ( flag_alph, flag_tmp)
00301 error_alph = dmax1(error_alph,error_tmp)
00302
00303 error_alph_last = error_alph
00304
00305 alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00306 call copy_to_mpatch_all_BNS_CF(impt)
00307 call copy_def_metric_and_matter_to_mpt(impt)
00308 call copy_def_matter_spin_to_mpt(impt)
00309 end if
00310
00311
00312 do ishift=1,shift_iter
00313 souvec(0:nrg,0:ntg,0:npg,1:3) = 0.0d0
00314 if (impt==1 .or. impt==2) then
00315 call sourceterm_MoC_CF_with_divshift_peos_spin(souvec)
00316 else
00317 call sourceterm_MoC_CF_with_divshift_r3rd(souvec)
00318 end if
00319
00320 do ii = 1, 3
00321 sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00322 sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00323 sou(0:nrg,0:ntg,0:npg) = souvec(0:nrg,0:ntg,0:npg,ii)
00324 if (ii.eq.1) work(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00325 if (ii.eq.2) work(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00326 if (ii.eq.3) work(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00327 if (ii.eq.1) then; chbe = 'bvxd'; chpa = 'od'; end if
00328 if (ii.eq.2) then; chbe = 'bvyd'; chpa = 'od'; end if
00329 if (ii.eq.3) then; chbe = 'bvzd'; chpa = 'ev'; end if
00330 if (impt==1 .or. impt==2) then
00331 chpB = 'ev'
00332 if (impt.eq.2.and.ii.eq.1) chpB = 'od'
00333 if (impt.eq.2.and.ii.eq.2) chpB = 'od'
00334 call sourceterm_exsurf_binary_COCP(impt,chbe,chpa,sou_exsurf, dsou_exsurf)
00335 call sourceterm_outsurf_COCP_from_ARCP(impt,chbe,chpB,sou_outsurf, dsou_outsurf)
00336 call poisson_solver_binary_star_homosol('sd',sou, sou_exsurf,dsou_exsurf, &
00337 & sou_outsurf,dsou_outsurf,pot)
00338 else if (impt.eq.nmpt) then
00339 call sourceterm_insurf_ARCP_from_COCP(impt,chbe,chpa,sou_insurf, dsou_insurf)
00340 if (ii.eq.1) call outer_boundary_d_bvxd(sou_outsurf)
00341 if (ii.eq.2) call outer_boundary_d_bvyd(sou_outsurf)
00342 if (ii.eq.3) call outer_boundary_d_bvzd(sou_outsurf)
00343
00344
00345
00346 call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf,dsou_insurf, &
00347 & sou_outsurf,dsou_outsurf,pot)
00348 end if
00349 if (ii.eq.1) potx(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00350 if (ii.eq.2) poty(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00351 if (ii.eq.3) potz(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00352 end do
00353
00354 potx_bak(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00355 poty_bak(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00356 potz_bak(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00357 call update_grfield(potx,bvxd,conv_gra)
00358 call update_grfield(poty,bvyd,conv_gra)
00359 call update_grfield(potz,bvzd,conv_gra)
00360
00361 if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00362 & call update_parameter_BNS_spin(conv_gra)
00363
00364 if (impt.ne.nmpt) then
00365 call interpolation_fillup_binary_COCP(impt,'bvxd','od',bvxd)
00366 call interpolation_fillup_binary_COCP(impt,'bvyd','od',bvyd)
00367 call interpolation_fillup_binary_COCP(impt,'bvzd','ev',bvzd)
00368 end if
00369
00370 call error_metric_type2_mpt(bvxd,potx_bak,error_tmp,flag_tmp,cocp,impt)
00371 flag_bvxd = max ( flag_bvxd, flag_tmp)
00372 error_bvxd = dmax1(error_bvxd,error_tmp)
00373 call error_metric_type2_mpt(bvyd,poty_bak,error_tmp,flag_tmp,cocp,impt)
00374 flag_bvyd = max ( flag_bvyd, flag_tmp)
00375 error_bvyd = dmax1(error_bvyd,error_tmp)
00376 call error_metric_type2_mpt(bvzd,potz_bak,error_tmp,flag_tmp,cocp,impt)
00377 flag_bvzd = max ( flag_bvzd, flag_tmp)
00378 error_bvzd = dmax1(error_bvzd,error_tmp)
00379
00380
00381
00382 end do
00383
00384
00385
00386
00387 call copy_to_mpatch_all_BNS_CF(impt)
00388 call copy_def_metric_and_matter_to_mpt(impt)
00389 call copy_def_matter_spin_to_mpt(impt)
00390
00391 error_metr = dmax1(error_psi,error_alph,error_bvxd,error_bvyd,error_bvzd)
00392
00393 if(iter_count<5 .and. iseq<=1 .and. indata_type.eq.'1D') then
00394 write(6,'(a14)') '#SKIP HYDRO a#'
00395 error_emd_last = 1.0d0
00396 cycle
00397 end if
00398
00399
00400
00401 if (error_emd_last.le.eps*0.0001d0 .and. error_vep_last.le.eps*0.0001d0 ) then
00402 write(6,'(a14)') '#SKIP HYDRO b#'
00403 else
00404 if (impt==1 .or. impt==2) then
00405 do ihy = 1, hydro_iter
00406 call interpo_gr2fl_metric_CF
00407 call hydrostatic_eq_CF_peos_spin(potf,potut,potux,potuy,potuz)
00408 call calc_surface(potrs, potf)
00409 emd_bak = emd
00410 call update_matter( potf,emd,conv_den)
00411
00412
00413
00414
00415 call update_surface(potrs,rs,conv_den)
00416
00417 call reset_fluid_BNS
00418 call update_parameter_BNS_spin(conv_den)
00419 call calc_vector_x_matter(2)
00420 call calc_vector_phi_matter(2)
00421
00422 call hydro_spbns_vep_CF_peos(potf,iter_count,impt,10)
00423 vep_bak = vep
00424 call update_matter(potf,vep,conv_vep)
00425 call reset_fluid_vep
00426 call calc_gradvep(vep,vepxf,vepyf,vepzf)
00427 call reset_fluid_gradvep
00428 call update_parameter_BNS_spin(conv_vep)
00429
00430 if (ihy.eq.hydro_iter) then
00431
00432 call error_matter_type2(emd,emd_bak,error_tmp,flag_tmp)
00433 flag_emd = max ( flag_emd, flag_tmp)
00434 error_emd = dmax1(error_emd,error_tmp)
00435 error_emd_last = error_emd
00436
00437
00438 call error_velocity_potential(vep,vep_bak,error_tmp,flag_tmp)
00439 flag_vep = max ( flag_vep, flag_tmp)
00440 error_vep = dmax1(error_vep,error_tmp)
00441 error_vep_last = error_vep
00442
00443 end if
00444 end do
00445 call copy_to_mpatch_all_BNS_CF(impt)
00446 call copy_def_metric_and_matter_to_mpt(impt)
00447 call copy_def_matter_spin_to_mpt(impt)
00448 else
00449 write(6,*) ' ## NO HYDRO ##'
00450 end if
00451 end if
00452
00453
00454 if(impt>5) then
00455 if (mod(iter_count,5)==0.or.iter_count==1) then
00456 write(char1, '(i5)') iter_count
00457 char2 = adjustl(char1)
00458 write(char4, '(i5)') impt
00459 char5 = adjustl(char4)
00460 char3 = 'iteration' // trim(char2) // '_mpt' // trim(char5) // '.txt'
00461 call printout_axis_spin_BNS_mpt(impt,char3)
00462 endif
00463 end if
00464
00465
00466 if (flag_mid==0) then
00467 if(error_metr < iter_eps ) then
00468 flag_mpt_mid(impt) = 1
00469 error_metric_mpt_mid(impt,1) = error_psi
00470 error_metric_mpt_mid(impt,2) = error_alph
00471 error_metric_mpt_mid(impt,3) = error_bvxd
00472 error_metric_mpt_mid(impt,4) = error_bvyd
00473 error_metric_mpt_mid(impt,5) = error_bvzd
00474 end if
00475 sumjj=0
00476 do jj=1,nmpt
00477 sumjj = sumjj + flag_mpt_mid(jj)
00478 end do
00479 if (sumjj==nmpt) then
00480 open(3,file='mid_error_data.txt',status='unknown')
00481 write(3,'(a11,i5)') 'Iteration =', iter_count
00482 write(3,'(a8,1p,5e17.6)') 'COCP #=1', (error_metric_mpt_mid(1,je), je=1,5)
00483 write(3,'(a8,1p,5e17.6)') 'COCP #=2', (error_metric_mpt_mid(2,je), je=1,5)
00484 write(3,'(a8,1p,5e17.6)') 'ARCP # ', (error_metric_mpt_mid(3,je), je=1,5)
00485 close(3)
00486
00487
00488
00489
00490
00491
00492
00493 flag_mid = 1
00494 end if
00495 end if
00496 end do
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510 if (flag_mid==1) then
00511 do impt=1, nmpt
00512 call copy_from_mpatch_all_BNS_CF(impt)
00513 write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00514 end do
00515 write(2,*) '-----------------------------------------------------------------------'
00516 close(2)
00517 exit
00518 end if
00519
00520 if ( (flag_vep==0).and.(flag_emd==0).and.(flag_psi==0).and.(flag_alph==0).and. &
00521 & (flag_bvxd==0).and.(flag_bvyd==0).and.(flag_bvzd==0)) then
00522 do impt=1, nmpt
00523 call copy_from_mpatch_all_BNS_CF(impt)
00524 write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00525 end do
00526 write(2,*) '-----------------------------------------------------------------------'
00527 close(2)
00528 exit
00529 endif
00530
00531 if (iter_count >= iter_max) then
00532 write(6,*) "*************** iter_count==iter_max ....exiting *******************"
00533 stop
00534
00535 end if
00536
00537
00538 flag_psi = 0
00539 flag_alph = 0
00540 flag_bvxd = 0
00541 flag_bvyd = 0
00542 flag_bvzd = 0
00543 flag_emd = 0
00544 flag_vep = 0
00545 error_psi = 0.0d0
00546 error_alph = 0.0d0
00547 error_bvxd = 0.0d0
00548 error_bvyd = 0.0d0
00549 error_bvzd = 0.0d0
00550 error_emd = 0.0d0
00551 error_vep = 0.0d0
00552 end do
00553
00554 deallocate(souvec)
00555 deallocate(sou)
00556 deallocate(pot)
00557 deallocate(potx)
00558 deallocate(poty)
00559 deallocate(potz)
00560 deallocate(pot_bak)
00561 deallocate(potx_bak)
00562 deallocate(poty_bak)
00563 deallocate(potz_bak)
00564 deallocate(work)
00565 deallocate(fnc_itp)
00566 deallocate(sou_exsurf)
00567 deallocate(dsou_exsurf)
00568 deallocate(sou_insurf)
00569 deallocate(dsou_insurf)
00570 deallocate(sou_outsurf)
00571 deallocate(dsou_outsurf)
00572 deallocate(potf)
00573 deallocate(potut)
00574 deallocate(potux)
00575 deallocate(potuy)
00576 deallocate(potuz)
00577 deallocate(emd_bak)
00578 deallocate(vep_bak)
00579 deallocate(potrs)
00580
00581 end subroutine iter_lecc_spin_BNS_CF_mpt
