00001 subroutine iteration_helmholtz_solver_test(iter_count)
00002 use phys_constant, only : long, nnrg, nntg, nnpg
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 def_metric, only : psi
00009 use def_matter, only : emd, emdg
00010 use make_array_2d
00011 use make_array_3d
00012 use interface_helmholtz_solver
00013 use interface_update_grfield
00014 use interface_error_metric
00015 use interface_sourceterm_helmholtz_solver_test
00016 use interface_helmholtz_solver
00017 implicit none
00018 real(long), pointer :: sou(:,:,:), pot(:,:,:), psi_bak(:,:,:)
00019 real(long), pointer :: sou_exsurf(:,:), dsou_exsurf(:,:)
00020 real(long) :: error_psi, count
00021 integer :: iter_count, flag = 0
00022 integer :: irf, itf, ipf, irg, itg, ipg, ihy
00023
00024 call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00025 call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00026 call alloc_array3d(psi_bak,0,nrg,0,ntg,0,npg)
00027 call alloc_array2d(sou_exsurf,0,ntg,0,npg)
00028 call alloc_array2d(dsou_exsurf,0,ntg,0,npg)
00029
00030 iter_count = 0
00031 do
00032 iter_count = iter_count + 1
00033 count = dble(iter_count)
00034 conv_gra = dmin1(conv0_gra,conv_ini*count)
00035 conv_den = dmin1(conv0_den,conv_ini*count)
00036
00037 call calc_vector_x_grav(1)
00038 call calc_vector_x_matter(1)
00039 call calc_vector_phi_grav(1)
00040 call calc_vector_phi_matter(1)
00041
00042 call sourceterm_helmholtz_solver_test(sou)
00043 call helmholtz_solver(sou,pot)
00044 psi_bak = psi
00045 call update_grfield(pot,psi,conv_gra)
00046 call error_metric(psi,psi_bak,error_psi,flag)
00047 call printout_error_metric(iter_count,error_psi)
00048
00049 if (flag == 0) exit
00050 if (iter_count >= iter_max) exit
00051 flag = 0
00052 end do
00053
00054 deallocate(sou)
00055 deallocate(pot)
00056 deallocate(sou_exsurf)
00057 deallocate(dsou_exsurf)
00058 deallocate(psi_bak)
00059 end subroutine iteration_helmholtz_solver_test
