Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Main_code/Main_spin_BNS_CF

00001 !______________________________________________
00002 include '../Include_file/include_modulefiles_BNS_CF_3mpt.f90'
00003 include '../Include_file/include_PEOS_modulefile.f90'
00004 include '../Include_file/include_interface_modulefiles_BNS_CF_3mpt.f90'
00005 include '../Include_file/include_subroutines_BNS_CF_3mpt.f90'
00006 include '../Include_file/include_PEOS_subroutines.f90'
00007 include '../Include_file/include_functions.f90'
00008 !______________________________________________
00009 !
00010 !              Main Program
00011 !______________________________________________
00012 PROGRAM Main_spin_BNS_CF_mpt
00013 !
00014   use phys_constant, only : long, nmpt
00015   use grid_parameter, only : indata_type, outdata_type, iter_max, &
00016   &                          num_sol_seq, sw_mass_iter, NS_shape
00017   use def_matter
00018   use def_matter_parameter, only : emdc, radi, ome
00019   use def_matter_parameter_mpt
00020   use def_quantities_mpt
00021   use def_peos_parameter, only : rhoini_cgs, emdini_gcm1
00022 !  use def_bh_parameter, only : bh_soltype
00023   use grid_parameter_binary_excision
00024   use grid_points_binary_excision
00025   use grid_points_asymptotic_patch
00026   use grid_points_binary_in_asympto
00027   use weight_midpoint_binary_excision
00028   use radial_green_fn_grav
00029   use radial_green_fn_grav_bhex_nb
00030   use radial_green_fn_grav_bhex_dd
00031 !  use radial_green_fn_grav_bhex_di
00032   use radial_green_fn_grav_bhex_nd
00033   use radial_green_fn_grav_bhex_dh
00034   use radial_green_fn_grav_bhex_nh
00035   use radial_green_fn_grav_bhex_sd
00036   use interface_copy_to_hgfn_and_gfnsf
00037   use interface_calc_gradvep
00038   use interface_calc_gradvep_from_corot_id
00039   implicit none
00040 
00041 !------------------------------------------------------------------------------------------------------------!
00042 !                                                                                                            !
00043 ! numseq=1 for a single solution with given TOV initial data or given 3D initial data                        !
00044 ! numseq>1 produces a sequence of solutions with increasing emdc using routine next_solution_BNS_mpt.f90     !
00045 ! In this case the iteration starts from TOV initial data (or 3D) computes the first spinning solution       !
00046 ! and then increases emdc to move to the next spinning solution.                                             !
00047 ! When computing parallel spinning solutions starting from different TOV use numseq=1 and the script         !
00048 ! debug_BNS_CF_peos_seq_from_TOV.sh                                                                          !
00049 !                                                                                                            !
00050 ! sw_mass_iter = 'n' implements the above scenarios                                                          !
00051 ! sw_mass_iter = 'y' finds a solution for a given rest mass.                                                 !
00052 !                                                                                                            !
00053 ! print_sol_seq=0 for NOT printing solution data files                                                       !
00054 !                                                                                                            !
00055 !------------------------------------------------------------------------------------------------------------!
00056 
00057   integer :: impt, itg, print_sol_seq=1, numseq=1
00058   integer :: iseq, iter_count, total_iteration
00059   character(30) :: char1, char2, char3, char4
00060   character(100) :: dircommand
00061   real(long) :: iter_eps
00062 !
00063 ! -- Read parameters
00064   call allocate_grid_parameter_mpt
00065   call allocate_grid_parameter_binary_excision_mpt
00066   call allocate_def_matter_parameter_mpt
00067   do impt = 1, nmpt
00068     call read_parameter_mpt(impt)
00069     call read_surf_parameter_mpt(impt)
00070     call copy_grid_parameter_to_mpt(impt)
00071     call read_parameter_binary_excision_mpt(impt)
00072     call copy_grid_parameter_binary_excision_to_mpt(impt)
00073     if (impt==1 .or. impt==2) then
00074       call peos_initialize_mpt(impt)
00075       call copy_def_peos_parameter_to_mpt(impt)
00076     end if
00077     call copy_grid_parameter_to_mpt(impt)
00078   end do
00079 !
00080 ! -- Allocate arrays
00081   call set_allocate_size_mpt
00082 !
00083   call allocate_coordinate_patch_kit_grav_mpt
00084   call allocate_grid_points_binary_excision
00085   call allocate_grid_points_asymptotic_patch
00086   call allocate_grid_points_binary_in_asympto
00087   call allocate_weight_midpoint_binary_excision
00088   call allocate_hgfn_bhex
00089   call allocate_hgfn_bhex_nb
00090   call allocate_hgfn_bhex_dd
00091   call allocate_hgfn_bhex_nd
00092 !  call allocate_hgfn_bhex_di
00093   call allocate_hgfn_bhex_sd
00094   call allocate_BNS_CF
00095   call allocate_metric_on_SFC_CF
00096   call allocate_BNS_CF_spin
00097 !
00098   call allocate_mpatch_all_BNS_CF
00099   call allocate_grid_points_asymptotic_patch_mpt
00100   call allocate_grid_points_binary_in_asympto_mpt
00101   call allocate_BNS_CF_mpt
00102   call allocate_BNS_CF_spin_mpt
00103 !
00104   do impt = 1, nmpt
00105     call copy_grid_parameter_from_mpt(impt)
00106     call copy_grid_parameter_binary_excision_from_mpt(impt)
00107     call copy_def_peos_parameter_from_mpt(impt)
00108     write(6,*) "*****************************************************************************"
00109     write(6,*) "* Change integer input (1,2,3,4) ALSO in following subroutines:             *"
00110     write(6,*) "* code/Main_utility/interpolation_contour_potential_irrot_BNS_CF_3mpt.f90   *"
00111     write(6,*) "* spherical_data/initial_isotropic_coord/code/Import_isotr_peos_surf.f90    *"
00112     write(6,*) "*****************************************************************************"
00113     call coordinate_patch_kit_grav_grid_mpt(3)       ! 3:r_surf is used
00114 
00115     call calc_parameter_binary_excision
00116     call IO_printout_grid_data_mpt(impt)
00117     if (impt.ne.nmpt) then
00118       call calc_grid_points_binary_excision
00119       call calc_grid_points_binary_in_asympto(impt,nmpt)
00120       call copy_grid_points_binary_in_asympto_to_mpt(impt)
00121     end if
00122     call calc_weight_midpoint_binary_excision
00123 !      call calc_weight_midpoint_binary_excision_hybrid
00124     if (impt==1 .or. impt==2) then
00125       call calc_vector_x_grav(2)
00126       call calc_vector_phi_grav(2)
00127     else
00128       call calc_vector_x_grav(1)
00129       call calc_vector_phi_grav(1)
00130     end if
00131 !    call calc_vector_x_matter(2)
00132 !    call calc_vector_phi_matter(2)
00133     call copy_to_mpatch_all_BNS_CF(impt)
00134     call copy_def_metric_and_matter_to_mpt(impt)
00135     call copy_def_matter_spin_to_mpt(impt)
00136   end do
00137   call copy_from_mpatch_all_BNS_CF(nmpt)
00138 ! -- coordinates of asymptotic patch in central patches
00139   do impt = 1, 2
00140     call calc_grid_points_asymptotic_patch(impt,nmpt)
00141     call copy_grid_points_asymptotic_patch_to_mpt(impt)
00142   end do
00143 !
00144 ! Initialization of variables
00145   do impt = 1, nmpt
00146     call copy_from_mpatch_all_BNS_CF(impt)
00147     call copy_def_metric_and_matter_from_mpt(impt)
00148     call copy_def_matter_spin_from_mpt(impt)
00149     if (indata_type.eq.'3D') then
00150       if(NS_shape.eq.'SP')  then
00151         write(6,*) "########################### Reading 3D Spinning Initial Data ######################################"
00152         call IO_input_initial_3D_CF_spin_NS_mpt(impt)
00153         if (impt==1 .or. impt==2)    then
00154           emdc = emd(0,0,0)
00155           call calc_gradvep(vep,vepxf,vepyf,vepzf)
00156         else
00157           emdc = 0.0d0
00158         end if
00159       else if (NS_shape.eq.'IR')  then
00160         write(6,*) "########################### Reading 3D Irrotational Initial Data ####################################"
00161         call IO_input_initial_3D_CF_irrot_NS_mpt(impt)
00162         if (impt==1 .or. impt==2)    then
00163           emdc = emd(0,0,0)
00164           call calc_gradvep(vep,vepxf,vepyf,vepzf)
00165         else
00166           emdc = 0.0d0
00167         end if
00168       else if(NS_shape.eq.'CO')  then
00169         write(6,*) "########################### Reading 3D Corotating Initial Data ####################################"
00170         call IO_input_initial_3D_CF_NS_mpt(impt)
00171         if (impt==1 .or. impt==2)   then
00172           emdc = emd(0,0,0)
00173           call calc_gradvep_from_corot_id(vep,vepxf,vepyf,vepzf)
00174         else
00175           emdc = 0.0d0
00176         end if
00177       else
00178         write(6,*) "**************** indata_type=3D and NS_shape is neither IR nor CO or SP....exiting"
00179         stop
00180       end if
00181     else                                     !  1D initial data 
00182       if (impt==1 .or. impt==2)  then
00183         call  IO_input_initial_1D_CF_NS_mpt(impt) 
00184         emdc = emd(0,0,0)
00185         call initial_velocity_potential_NS_mpt(impt)
00186         if(NS_shape.eq.'SP')   call IO_input_initial_vrot_NS_mpt(impt) 
00187       else
00188         emdc = 0.0d0
00189       end if
00190     end if
00191     write (6,'(a8,i2,a16,1p,e20.12)')  "--Patch=",impt,"           emdc=",emdc
00192 !
00193     call copy_def_metric_and_matter_to_mpt(impt)
00194     call copy_def_matter_spin_to_mpt(impt)
00195     call copy_to_mpatch_all_BNS_CF(impt)
00196   end do  ! =========>>> end initialization of all patches
00197 !
00198   call copy_grid_parameter_from_mpt(1)  
00199   if ( sw_mass_iter=='y' ) then
00200     call calc_spin_BNS_CF_mpt(total_iteration)
00201     call calc_physical_quantities_spin_BNS_CF_mpt
00202     char3 = 'main_bnsphys_all_mpt.txt'
00203     call printout_physq_BNS_all_mpt(char3)
00204     call write_last_physq_BNS_mpt
00205     do impt = 1, nmpt
00206       call copy_from_mpatch_all_BNS_CF(impt)
00207       call copy_def_metric_and_matter_from_mpt(impt)
00208       call copy_def_matter_spin_from_mpt(impt)
00209       if (outdata_type.eq.'3D')  call IO_output_solution_3D_CF_spin_NS_mpt(impt,'.las')
00210       if (impt==1 .or. impt==2)  call printout_NS_shape_mpt(impt)
00211     end do
00212   else
00213     do iseq = 1, numseq   !num_sol_seq
00214       write(char1, '(i5)') iseq
00215       char2 = adjustl(char1)
00216       if (iseq < 10)  then
00217         char3 = 'iseq0' // trim(char2)
00218       else
00219         char3 = 'iseq'  // trim(char2)
00220       end if
00221       dircommand = 'mkdir ' // char3
00222       call system(dircommand)
00223 !      dircommand = 'cp -f peos_parameter_mpt* ./' // char3
00224 !      call system(dircommand)
00225       call chdir(char3)
00226 !
00227       write(6,'(a50,i2,a32)') '============================== Solution sequence #',iseq,' ============================== '
00228       write(6,'(a14,1p,2e20.12)')  "emdc COCP1,2: ", def_matter_param_real_(2,1), def_matter_param_real_(2,2)
00229 !     Change the argument iseq to a number bigger than 1 (ex 2) if you want to  compute hydro from first iteration. 
00230 !     Usually this is the case when the initial data are 3D.
00231       iter_eps = 1.0d-02
00232       call iter_spin_BNS_CF_mpt(iter_count,iseq, iter_eps)
00233       call calc_physical_quantities_spin_BNS_CF_mpt
00234       char3 = 'main_bnsphys_all_mpt.txt'
00235       call printout_physq_BNS_all_mpt(char3)
00236       call write_last_physq_BNS_mpt
00237       if (print_sol_seq==1) then
00238         do impt = 1, nmpt
00239           call copy_from_mpatch_all_BNS_CF(impt)
00240           call copy_def_metric_and_matter_from_mpt(impt)
00241           call copy_def_matter_spin_from_mpt(impt)
00242           if (outdata_type.eq.'3D') call IO_output_solution_3D_CF_spin_NS_mpt(impt,'.las')
00243           if (impt==1 .or. impt==2)  call printout_NS_shape_mpt(impt)
00244         end do
00245       end if
00246       call chdir('../')   ! now we are in work_area_BNS
00247       if ( iseq < numseq )   call next_solution_BNS_mpt(iseq)
00248     end do
00249   end if 
00250 !
00251 END PROGRAM Main_spin_BNS_CF_mpt
/home/uryu/Cocal_svn_checkinout/Cocal/code/Main_code/Main_spin_BNS_CF_mpt.f90
