00001 subroutine iteration_poisson_bbh_2pot_test_3mpt(iter_count)
00002 use phys_constant, only : long, nnrg, nntg, nnpg, nmpt
00003 use grid_parameter
00004 use grid_parameter_interpo, only : nrg_itp, ntg_itp, npg_itp
00005 use coordinate_grav_r
00006 use coordinate_grav_phi
00007 use coordinate_grav_theta
00008 use weight_midpoint_grav
00009 use def_metric, only : psi, alph
00010 use def_metric_mpt, only : psi_, alph_
00011 use copy_array_4dto3d_mpt
00012 use make_array_2d
00013 use make_array_3d
00014 use interface_interpo_fl2gr
00015 use interface_sourceterm_poisson_solver_test
00016 use interface_sourceterm_exsurf_eqm_binary
00017 use interface_sourceterm_surface_int
00018
00019 use interface_sourceterm_2pot_test
00020 use interface_sourceterm_insurf_asympto_interpo_from_mpt
00021 use interface_sourceterm_outsurf_interpo_from_asympto_mpt
00022 use interface_poisson_solver
00023 use interface_poisson_solver_binary_bhex_homosol
00024 use interface_poisson_solver_asymptotic_patch_homosol
00025 use interface_update_grfield
00026 use interface_error_metric
00027 use interface_error_metric_type0
00028 use interface_bh_boundary_test_mpt
00029 use interface_bh_boundary_psi_test_mpt
00030 use interface_bh_boundary_alph_test_mpt
00031 use interface_bh_boundary_nh_psi_test
00032 use interface_bh_boundary_nh_alph_test
00033 use interface_interpolation_fillup_binary_mpt
00034 use interface_bh_boundary_test
00035 implicit none
00036 real(long), pointer :: sou(:,:,:)
00037 real(long), pointer :: pot(:,:,:), psi_bak(:,:,:), alph_bak(:,:,:)
00038 real(long), pointer :: sou_exsurf(:,:), dsou_exsurf(:,:)
00039 real(long), pointer :: sou_bhsurf(:,:), dsou_bhsurf(:,:)
00040 real(long), pointer :: sou_insurf(:,:), dsou_insurf(:,:)
00041 real(long), pointer :: sou_outsurf(:,:), dsou_outsurf(:,:)
00042 real(long), pointer :: fnc_itp(:,:,:)
00043 real(long) :: error_psi =0.0d0, error_alph =0.0d0, error_tmp =0.0d0, count
00044 integer :: iter_count, flag_psi = 0, flag_alph = 0, flag_tmp = 0
00045 integer :: irf, itf, ipf, irg, itg, ipg, ihy
00046 integer :: impt, impt_ex, impt_bin
00047
00048 call set_allocate_size_mpt
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(psi_bak,0,nrg,0,ntg,0,npg)
00052 call alloc_array3d(alph_bak,0,nrg,0,ntg,0,npg)
00053 call alloc_array3d(fnc_itp,0,nrg,0,ntg,0,npg)
00054 call alloc_array2d(sou_exsurf,0,ntg,0,npg)
00055 call alloc_array2d(dsou_exsurf,0,ntg,0,npg)
00056 call alloc_array2d(sou_bhsurf,0,ntg,0,npg)
00057 call alloc_array2d(dsou_bhsurf,0,ntg,0,npg)
00058 call alloc_array2d(sou_insurf,0,ntg,0,npg)
00059 call alloc_array2d(dsou_insurf,0,ntg,0,npg)
00060 call alloc_array2d(sou_outsurf,0,ntg,0,npg)
00061 call alloc_array2d(dsou_outsurf,0,ntg,0,npg)
00062
00063 iter_count = 0
00064
00065 do
00066 iter_count = iter_count + 1
00067 count = dble(iter_count)
00068
00069 do impt = 1, nmpt
00070 call copy_grid_parameter_from_mpt(impt)
00071 conv_gra = dmin1(conv0_gra,conv_ini*count)
00072 conv_den = dmin1(conv0_den,conv_ini*count)
00073 call copy_grid_parameter_to_mpt(impt)
00074
00075 call copy_from_mpatch_all_test(impt)
00076 call calc_vector_x_grav(0)
00077 call calc_vector_x_matter(0)
00078 call calc_vector_phi_grav(0)
00079 call calc_vector_phi_matter(0)
00080 call copy_to_mpatch_all_test(impt)
00081
00082
00083 sou(0:nrg,0:ntg,0:npg) = 0.0d0
00084 sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00085 sou_bhsurf(0:ntg,0:npg) = 0.0d0 ; dsou_bhsurf(0:ntg,0:npg) = 0.0d0
00086 sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00087
00088 if (impt.ne.nmpt) then
00089 if (impt.eq.1) impt_ex = 2
00090 if (impt.eq.2) impt_ex = 1
00091 call copy_from_mpatch_all_test(impt_ex)
00092 call copy_metric_and_matter_BHNS_test_from_mpt(impt_ex)
00093 call sourceterm_exsurf_eqm_binary(psi,sou_exsurf,dsou_exsurf)
00094 call copy_from_mpatch_all_test(impt)
00095 call copy_metric_and_matter_BHNS_test_from_mpt(impt)
00096
00097
00098
00099
00100
00101
00102 call copy_grid_parameter_interpo_from_mpt(nmpt)
00103 call copy_coordinate_grav_extended_interpo_from_mpt(nmpt)
00104 call copy_array4dto3d_mpt(nmpt,psi_,fnc_itp, &
00105 & 0,nrg_itp,0,ntg_itp,0,npg_itp)
00106 call copy_grid_points_binary_in_asympto_from_mpt(impt)
00107 call sourceterm_outsurf_interpo_from_asympto_mpt &
00108 & (impt,nmpt,fnc_itp,sou_outsurf,dsou_outsurf)
00109
00110
00111 call bh_boundary_psi_test_mpt(impt,'d',sou_bhsurf,dsou_bhsurf)
00112 call poisson_solver_binary_bhex_homosol('dd',sou, &
00113 & sou_exsurf,dsou_exsurf, &
00114 & sou_bhsurf,dsou_bhsurf, &
00115 & sou_outsurf,dsou_outsurf,pot)
00116 end if
00117 if (impt.eq.nmpt) then
00118 call copy_from_mpatch_all_test(impt)
00119 call copy_metric_and_matter_BHNS_test_from_mpt(impt)
00120 do impt_bin = 1, nmpt - 1
00121 call copy_grid_parameter_interpo_from_mpt(impt_bin)
00122 call copy_coordinate_grav_extended_interpo_from_mpt(impt_bin)
00123 call copy_array4dto3d_mpt(impt_bin,psi_,fnc_itp, &
00124 & 0,nrg_itp,0,ntg_itp,0,npg_itp)
00125 call copy_grid_points_asymptotic_patch_from_mpt(impt_bin)
00126 call sourceterm_insurf_asympto_interpo_from_mpt &
00127 & (impt_bin,impt,fnc_itp,sou_insurf,dsou_insurf)
00128 if (impt_bin.eq.1) then
00129 sou_bhsurf(0:ntg,0:npg) = sou_insurf(0:ntg,0:npg)
00130 dsou_bhsurf(0:ntg,0:npg) = dsou_insurf(0:ntg,0:npg)
00131 end if
00132 if (impt_bin.eq.2) then
00133 sou_insurf(0:ntg,0:npg) = &
00134 & 0.5d0*( sou_insurf(0:ntg,0:npg) + sou_bhsurf(0:ntg,0:npg))
00135 dsou_insurf(0:ntg,0:npg) = &
00136 & 0.5d0*(dsou_insurf(0:ntg,0:npg) + dsou_bhsurf(0:ntg,0:npg))
00137 end if
00138 end do
00139 call sourceterm_surface_int(psi,nrg,sou_outsurf,dsou_outsurf)
00140 call poisson_solver_asymptotic_patch_homosol('dd',sou, &
00141 & sou_insurf,dsou_insurf, &
00142 & sou_outsurf,dsou_outsurf,pot)
00143 end if
00144
00145 psi_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00146 call update_grfield(pot,psi,conv_gra)
00147
00148 if (impt.ne.nmpt) then
00149 call copy_grid_parameter_interpo_from_mpt(impt_ex)
00150 call copy_grid_parameter_binary_excision_interpo_from_mpt(impt_ex)
00151 call copy_coordinate_grav_extended_interpo_from_mpt(impt_ex)
00152 call copy_array4dto3d_mpt(impt_ex,psi_,fnc_itp, &
00153 & 0,nrg_itp,0,ntg_itp,0,npg_itp)
00154 call interpolation_fillup_binary_mpt(psi,fnc_itp)
00155 end if
00156
00157 call error_metric_type0(psi,psi_bak,error_tmp,flag_tmp,'bh')
00158 flag_psi = max ( flag_psi, flag_tmp)
00159 error_psi = dmax1(error_psi,error_tmp)
00160 call printout_error_metric(iter_count,error_psi)
00161 if (impt.eq.nmpt) call reset_bh_boundary('n')
00162 call copy_to_mpatch_all_test(impt)
00163 call copy_metric_and_matter_BHNS_test_to_mpt(impt)
00164
00165
00166 sou(0:nrg,0:ntg,0:npg) = 0.0d0
00167 sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00168 sou_bhsurf(0:ntg,0:npg) = 0.0d0 ; dsou_bhsurf(0:ntg,0:npg) = 0.0d0
00169 sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00170
00171 if (impt.ne.nmpt) then
00172 if (impt.eq.1) impt_ex = 2
00173 if (impt.eq.2) impt_ex = 1
00174 call copy_from_mpatch_all_test(impt_ex)
00175 call copy_metric_and_matter_BHNS_test_from_mpt(impt_ex)
00176 call sourceterm_exsurf_eqm_binary(alph,sou_exsurf,dsou_exsurf)
00177 call copy_from_mpatch_all_test(impt)
00178 call copy_metric_and_matter_BHNS_test_from_mpt(impt)
00179
00180
00181
00182
00183
00184
00185 call copy_grid_parameter_interpo_from_mpt(nmpt)
00186 call copy_coordinate_grav_extended_interpo_from_mpt(nmpt)
00187 call copy_array4dto3d_mpt(nmpt,alph_,fnc_itp, &
00188 & 0,nrg_itp,0,ntg_itp,0,npg_itp)
00189 call copy_grid_points_binary_in_asympto_from_mpt(impt)
00190 call sourceterm_outsurf_interpo_from_asympto_mpt &
00191 & (impt,nmpt,fnc_itp,sou_outsurf,dsou_outsurf)
00192
00193 call bh_boundary_alph_test_mpt(impt,'d',sou_bhsurf,dsou_bhsurf)
00194
00195 call sourceterm_2pot_test(sou)
00196
00197 call poisson_solver_binary_bhex_homosol('dd',sou, &
00198 & sou_exsurf,dsou_exsurf, &
00199 & sou_bhsurf,dsou_bhsurf, &
00200 & sou_outsurf,dsou_outsurf,pot)
00201 end if
00202 if (impt.eq.nmpt) then
00203 call copy_from_mpatch_all_test(impt)
00204 call copy_metric_and_matter_BHNS_test_from_mpt(impt)
00205 do impt_bin = 1, nmpt - 1
00206 call copy_grid_parameter_interpo_from_mpt(impt_bin)
00207 call copy_coordinate_grav_extended_interpo_from_mpt(impt_bin)
00208 call copy_array4dto3d_mpt(impt_bin,alph_,fnc_itp, &
00209 & 0,nrg_itp,0,ntg_itp,0,npg_itp)
00210 call copy_grid_points_asymptotic_patch_from_mpt(impt_bin)
00211 call sourceterm_insurf_asympto_interpo_from_mpt &
00212 & (impt_bin,impt,fnc_itp,sou_insurf,dsou_insurf)
00213 if (impt_bin.eq.1) then
00214 sou_bhsurf(0:ntg,0:npg) = sou_insurf(0:ntg,0:npg)
00215 dsou_bhsurf(0:ntg,0:npg) = dsou_insurf(0:ntg,0:npg)
00216 end if
00217 if (impt_bin.eq.2) then
00218 sou_insurf(0:ntg,0:npg) = &
00219 & 0.5d0*( sou_insurf(0:ntg,0:npg) + sou_bhsurf(0:ntg,0:npg))
00220 dsou_insurf(0:ntg,0:npg) = &
00221 & 0.5d0*(dsou_insurf(0:ntg,0:npg) + dsou_bhsurf(0:ntg,0:npg))
00222 end if
00223 end do
00224 call sourceterm_surface_int(alph,nrg,sou_outsurf,dsou_outsurf)
00225 sou(0:nrg,0:ntg,0:npg) = 0.0d0
00226
00227 call sourceterm_2pot_test(sou)
00228 call poisson_solver_asymptotic_patch_homosol('dd',sou, &
00229 & sou_insurf,dsou_insurf, &
00230 & sou_outsurf,dsou_outsurf,pot)
00231 end if
00232
00233 alph_bak(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)
00234 call update_grfield(pot,alph,conv_gra)
00235
00236 if (impt.ne.nmpt) then
00237 call copy_grid_parameter_interpo_from_mpt(impt_ex)
00238 call copy_grid_parameter_binary_excision_interpo_from_mpt(impt_ex)
00239 call copy_coordinate_grav_extended_interpo_from_mpt(impt_ex)
00240 call copy_array4dto3d_mpt(impt_ex,alph_,fnc_itp, &
00241 & 0,nrg_itp,0,ntg_itp,0,npg_itp)
00242 call interpolation_fillup_binary_mpt(alph,fnc_itp)
00243 end if
00244
00245 call error_metric_type0(alph,alph_bak,error_tmp,flag_tmp,'bh')
00246 flag_alph = max ( flag_alph, flag_tmp)
00247 error_alph = dmax1(error_alph,error_tmp)
00248 call printout_error_metric(iter_count,error_alph)
00249 if (impt.eq.nmpt) call reset_bh_boundary('n')
00250 call copy_to_mpatch_all_test(impt)
00251 call copy_metric_and_matter_BHNS_test_to_mpt(impt)
00252 end do
00253
00254 if ((flag_psi==0).and.(flag_alph==0)) exit
00255 if (iter_count >= iter_max) exit
00256 if (iter_count >= 150 .and. error_psi > 1.5d0) exit
00257 if (iter_count >= 150 .and. error_alph > 1.5d0) exit
00258 flag_psi = 0
00259 flag_alph = 0
00260 error_psi = 0.0d0
00261 error_alph = 0.0d0
00262 end do
00263
00264 deallocate(sou)
00265 deallocate(pot)
00266 deallocate(psi_bak)
00267 deallocate(alph_bak)
00268 deallocate(fnc_itp)
00269 deallocate(sou_exsurf)
00270 deallocate(dsou_exsurf)
00271 deallocate(sou_bhsurf)
00272 deallocate(dsou_bhsurf)
00273 deallocate(sou_insurf)
00274 deallocate(dsou_insurf)
00275 deallocate(sou_outsurf)
00276 deallocate(dsou_outsurf)
00277 end subroutine iteration_poisson_bbh_2pot_test_3mpt
