Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine_mpatch/calc_lecc_spin_BNS_CF

00001 subroutine calc_lecc_spin_BNS_CF_mpt(total_iteration)
00002   use phys_constant, only  : long, nmpt
00003   use grid_parameter
00004   use def_quantities
00005   use def_matter_parameter
00006   use def_matter_parameter_mpt
00007   use interface_violation_midpoint_MoC_CF_peos_spin
00008   use interface_violation_gridpoint_MoC_CF_peos_spin
00009   use interface_violation_midpoint_HaC_CF_peos_spin
00010   use interface_violation_gridpoint_HaC_CF_peos_spin
00011   use interface_IO_output_3D_general
00012   use interface_IO_output_2D_general
00013   use interface_IO_output_1D_general
00014   use make_array_2d
00015   use make_array_3d
00016   use make_array_4d
00017   implicit none
00018   real(long), pointer :: pot(:,:,:), HaC_vio(:,:,:), MoC_vio(:,:,:,:)
00019   real(long) :: xa,xb,fa,fb,rm_eps,delta,f1,f2, iter_eps
00020   integer :: total_iteration, iter_count, i, icycle, Nmax=15
00021   integer :: impt
00022   character(30) :: char1, char2, char3, char4, char5
00023 
00024   rm_eps = 1.0d-06
00025   total_iteration = 0
00026   char3 = 'main_bnsphys_all_mpt.txt'
00027 !
00028   call open_directory_mpt(0)
00029   iter_eps = 1.0d0
00030   call iter_lecc_spin_BNS_CF_mpt(iter_count,1,iter_eps)     !  1 for freezing hydro in first 5 iterations
00031   total_iteration = total_iteration + iter_count
00032   write(6,'(a21,i5)') '--- total iteration =', total_iteration
00033   call calc_physical_quantities_spin_BNS_CF_mpt   
00034   call write_last_physq_BNS_mpt
00035   call printout_physq_BNS_all_mpt(char3)
00036 !  call save_solution_BNS_mpt(0)
00037   call copy_def_matter_parameter_from_mpt(1)
00038   xa = emdc
00039   call copy_def_quantities_BNS_from_mpt(1)
00040   fa = (restmass - restmass_sph)/restmass_sph
00041   write(6,*)  "xa,fa=", xa, fa,  "           restmass, restmass_sph = ", restmass, restmass_sph
00042   call chdir('../')
00043 !
00044   call open_directory_mpt(1)
00045 !  xb = xa - 0.25d0*xa
00046   xb = xa - 0.01d0*xa
00047   do impt = 1, 2
00048     call copy_def_matter_parameter_from_mpt(impt)
00049     emdc = xb
00050     call copy_def_matter_parameter_to_mpt(impt)
00051   end do
00052   write(6,'(a14,1p,2e20.12)')  "emdc COCP1,2: ", def_matter_param_real_(2,1), def_matter_param_real_(2,2)
00053   iter_eps = 1.0d0
00054   call iter_lecc_spin_BNS_CF_mpt(iter_count,2,iter_eps)     ! 2 to compute hydro from first iteration
00055   total_iteration = total_iteration + iter_count
00056   write(6,'(a21,i5)') '--- total iteration =', total_iteration
00057   call calc_physical_quantities_spin_BNS_CF_mpt
00058   call write_last_physq_BNS_mpt
00059   call printout_physq_BNS_all_mpt(char3)
00060 !  call save_solution_BNS_mpt(1)
00061   call copy_def_quantities_BNS_from_mpt(1)
00062   fb = (restmass - restmass_sph)/restmass_sph
00063   write(6,*)  "xb,fb=", xb, fb,  "           restmass, restmass_sph = ", restmass, restmass_sph
00064   call chdir('../')
00065 !
00066   write(6,*) '--------- Initial values for emdc and f(emdc)=(restmass - restmass_sph)/restmass_sph ---------------------------'
00067   write(6,'(a30,i2,1p,2e20.12)')  '(cycle,emdc,(M0-M0sph)/M0sph)=', 0,xa,fa
00068   write(6,'(a30,i2,1p,2e20.12)')  '(cycle,emdc,(M0-M0sph)/M0sph)=', 1,xb,fb
00069 !
00070   if ( dabs(fa)>dabs(fb) ) then  ! interchange xa,xb  and fa,fb
00071     f1 = xa
00072     xa = xb
00073     xb = f1
00074 
00075     f1 = fa
00076     fa = fb
00077     fb = f1
00078   endif
00079 !
00080   icycle = 1
00081   do i = 2,Nmax
00082     icycle = icycle + 1
00083     if ( dabs(fa)>dabs(fb) ) then
00084       f1 = xa
00085       xa = xb
00086       xb = f1
00087       f1 = fa
00088       fa = fb
00089       fb = f1
00090     end if
00091     delta = (xb-xa)/(fb-fa)
00092     xb = xa
00093     fb = fa
00094     delta = delta*fa
00095     if ( dabs(delta) < rm_eps ) then
00096       write(6,*) 'Convergence...'
00097       call copy_def_matter_parameter_from_mpt(1)
00098       write(6,*) '------------------------------------------------------------------'
00099       write(6,*) 'Central Emden = ',emdc
00100 
00101 !      write(6,*) "Calculating constraint violations at gridpoints of COCP-1..."
00102 !      call copy_from_mpatch_all_BNS_CF(1)
00103 !      call copy_def_metric_and_matter_from_mpt(1)
00104 !      call copy_def_matter_spin_from_mpt(1)
00105 !      call alloc_array4d(MoC_vio,0,nrg,0,ntg,0,npg,1,3)
00106 !      call alloc_array3d(HaC_vio,0,nrg,0,ntg,0,npg)
00107 !      call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00108 !      call excurve_CF('ns')             !   3rd order midpoint from ir0=-2,...
00109 !      call excurve_CF_gridpoint         !   4th order from ir0=-2,...
00110 
00111 !      MoC_vio(0:nrg,0:ntg,0:npg,1:3) = 0.0d0
00112 !      call violation_gridpoint_MoC_CF_peos_spin(MoC_vio,'ns')
00113 !      pot = 0.0d0
00114 !      pot(0:nrg,0:ntg,0:npg) = MoC_vio(0:nrg,0:ntg,0:npg,2)
00115 !      char1 = 'MoC_by_3D_mpt1.txt'
00116 !      call IO_output_3D_general(char1,'g','g',pot)
00117 !      char1 = 'MoC_by_xy_mpt1.txt'
00118 !      call IO_output_2D_general(char1,'g','g',pot,'xy')
00119 !      char1 = 'MoC_by_phi000_mpt1.txt'
00120 !      call IO_output_1D_general(char1,'g','g',pot,-1,ntg/2,0)
00121 !      char1 = 'MoC_by_phi180_mpt1.txt'
00122 !      call IO_output_1D_general(char1,'g','g',pot,-1,ntg/2,npg/2)
00123 
00124 !      HaC_vio = 0.0d0
00125 !      call violation_gridpoint_HaC_CF_peos_spin(HaC_vio,'ns')
00126 !      char1 = 'HaC_3D_mpt1.txt'
00127 !      call IO_output_3D_general(char1,'g','g',HaC_vio)
00128 !      char1 = 'HaC_xy_mpt1.txt'
00129 !      call IO_output_2D_general(char1,'g','g',HaC_vio,'xy')
00130 !      char1 = 'HaC_phi000_mpt1.txt'
00131 !      call IO_output_1D_general(char1,'g','g',HaC_vio,-1,ntg/2,0)
00132 !      char1 = 'HaC_phi180_mpt1.txt'
00133 !      call IO_output_1D_general(char1,'g','g',HaC_vio,-1,ntg/2,npg/2)
00134 
00135 !      deallocate(MoC_vio)
00136 !      deallocate(HaC_vio)
00137 !      deallocate(pot)
00138 !      call printout_physq_console_BBH
00139       return
00140     end if
00141 !    xa = xa - 0.1d0*delta
00142     if (i.lt.10) then
00143       xa = xa - i*delta/10.0d0
00144     else
00145       xa = xa - delta
00146     end if 
00147     if (i.lt.7) then
00148       iter_eps = dmin1(10.0d0**(-i), dabs(delta))
00149     else
00150       iter_eps = rm_eps
00151     end if
00152 
00153     call open_directory_mpt(icycle)
00154     do impt = 1, 2
00155       call copy_def_matter_parameter_from_mpt(impt)
00156       emdc = xa
00157       call copy_def_matter_parameter_to_mpt(impt)
00158     end do
00159     write(6,'(a14,1p,2e20.12)')  "emdc COCP1,2: ", def_matter_param_real_(2,1), def_matter_param_real_(2,2)
00160     write(6,*) '------------------------------------------------------------------'
00161 !    write (6,'(a7,i2,a20,e20.12)') 'cycle =', icycle, '        New emdc  = ',emdc
00162 
00163     call iter_lecc_spin_BNS_CF_mpt(iter_count,icycle,iter_eps)
00164     total_iteration = total_iteration + iter_count
00165     write(6,'(a21,i5)') '--- total iteration =', total_iteration
00166     call calc_physical_quantities_spin_BNS_CF_mpt
00167     call write_last_physq_BNS_mpt
00168     call printout_physq_BNS_all_mpt(char3)
00169 !    call save_solution_BNS_mpt(icycle)
00170     call copy_def_quantities_BNS_from_mpt(1)
00171     fa = (restmass - restmass_sph)/restmass_sph
00172     write(6,'(a30,i2,1p,2e20.12)')   '(cycle,emdc,(M0-M0sph)/M0sph)=',icycle,xa,fa
00173     call chdir('../')
00174 !
00175   end do
00176   write(6,*) 'No convergence after ',Nmax,' iterations.'
00177 
00178 end subroutine calc_lecc_spin_BNS_CF_mpt
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine_mpatch/calc_lecc_spin_BNS_CF_mpt.f90
