diff --git a/Makefile b/Makefile index 16bac61..f8004dc 100644 --- a/Makefile +++ b/Makefile @@ -16,7 +16,7 @@ vpath %.f90 $(SRC) INCFFTW= /usr/include LIBFFTW= /lib/x86_64-linux-gnu #========================================================================== -FC = ifort#gfortran# +FC = gfortran#ifort# FFLAGS= -g -fmax-errors=5 -O2 diff --git a/setup.py b/setup.py index 6ba43e2..c6d5fda 100755 --- a/setup.py +++ b/setup.py @@ -16,11 +16,11 @@ sim.black_hole_spin= round(0.7*sim.black_hole_mass,4) sim.compactification_length= float(1) #============================================================================= -sim.evolve_time= float(150) ## units of black hole mass -sim.num_saved_times= int(500) +sim.evolve_time= float(20) ## units of black hole mass +sim.num_saved_times= int(40) #============================================================================= -sim.nx= 112 ## num radial pts -sim.nl= 28 ## num angular values +sim.nx= 64 ## num radial pts +sim.nl= 16 ## num angular values #============================================================================= ## evolution and write: take boolean values @@ -33,7 +33,7 @@ sim.write_scd_order_source= True sim.write_coefs= False #============================================================================= -sim.computer= 'feynman'#'home'# +sim.computer= 'home'#'feynman'# sim.feyn_out_stem= '/mnt/grtheory/tf-out/' ## for feynman cluster/slurm script @@ -56,49 +56,63 @@ sim.start_multiple= float(1.0) #============================================================================= +## Initial data +#============================================================================= +## p(m)n: +(-) m angular number +## l_ang: initial data is a particular swal function +## initial_data_direction: which way pulse is approximately "heading" +## amp_re(im)_pm: initial amplitude of real/imaginary parts of psi4 +## rl(ru)_pm_0: lower(upper) bounds of initial data as a multiple +## of the black hole horizon +#============================================================================= ## initial data for mode m1 #============================================================================= -sim.pm1_ang = int(2) ## m_ang is preserved by time evolution +sim.pm1_ang = int(2) #----------------------------------------------------------------------------- -sim.l_ang_pm1= int(2) ## support over single spin weighted spherical harmonic +sim.l_ang_pm1= int(2) sim.initial_data_direction_pm1= "ingoing"#"outgoing"#"time_symmetric"# -sim.amp_pm1= float( 0.1) ## amplitude of the initial perturbation +sim.amp_re_pm1= float(0.1) +sim.amp_im_pm1= float(0.1) -sim.rl_pm1_0= float(-2.0) ## lower r value as a multiple of horizon -sim.ru_pm1_0= float( 2.0) ## upper r value as a multiple of horizon +sim.rl_pm1_0= float(-2.0) +sim.ru_pm1_0= float( 2.0) #----------------------------------------------------------------------------- -sim.l_ang_nm1= int(2) ## support over single spin weighted spherical harmonic +sim.l_ang_nm1= int(2) sim.initial_data_direction_nm1= "ingoing"#"time_symmetric"#"outgoing"# -sim.amp_nm1= float( 0.0) ## amplitude of the initial perturbation +sim.amp_re_nm1= float( 0.0) +sim.amp_im_nm1= float( 0.0) -sim.rl_nm1_0= float( 1.1) ## lower r value as multiple of horizon -sim.ru_nm1_0= float( 3.0) ## upper r value as multiple of horizon +sim.rl_nm1_0= float( 1.1) +sim.ru_nm1_0= float( 3.0) #============================================================================= ## initial data for mode m2 #============================================================================= sim.pm2_ang = int(3) ## m_ang is preserved by time evolution #----------------------------------------------------------------------------- -sim.l_ang_pm2= int(3) ## support over single spin weighted spherical harmonic +sim.l_ang_pm2= int(3) sim.initial_data_direction_pm2= "ingoing"#"time_symmetric"#"outgoing"# -sim.amp_pm2= float(0.05) ## amplitude of the initial perturbation +sim.amp_re_pm2= float(0.05) +sim.amp_im_pm2= float(0.05) -sim.rl_pm2_0= float(-1.5) ## lower r value as multiple of horizon -sim.ru_pm2_0= float( 1.5) ## upper r value as multiple of horizon +sim.rl_pm2_0= float(-1.5) +sim.ru_pm2_0= float( 1.5) #----------------------------------------------------------------------------- -sim.l_ang_nm2= int(3) ## support over single spin weighted spherical harmonic +sim.l_ang_nm2= int(3) +sim.l_ang_nm2= int(3) sim.initial_data_direction_nm2= "ingoing"#"time_symmetric"#"outgoing"# -sim.amp_nm2= float(0.0) ## amplitude of the initial perturbation +sim.amp_re_nm2= float(0.0) +sim.amp_im_nm2= float(0.0) -sim.rl_nm2_0= float(-1.5) ## lower r value as multiple of horizon -sim.ru_nm2_0= float( 1.5) ## upper r value as multiple of horizon +sim.rl_nm2_0= float(-1.5) +sim.ru_nm2_0= float( 1.5) #============================================================================= ## which m angular values to evolve diff --git a/src/mod_initial_data.f90 b/src/mod_initial_data.f90 index 5bb156b..bc5b239 100644 --- a/src/mod_initial_data.f90 +++ b/src/mod_initial_data.f90 @@ -16,14 +16,18 @@ module mod_initial_data pm1_ang, & initial_data_direction_nm1, & initial_data_direction_pm1, & - amp_nm1, rl_nm1, ru_nm1, l_ang_nm1, & - amp_pm1, rl_pm1, ru_pm1, l_ang_pm1, & + amp_re_nm1, rl_nm1, ru_nm1, l_ang_nm1, & + amp_im_nm1, & + amp_re_pm1, rl_pm1, ru_pm1, l_ang_pm1, & + amp_im_pm1, & pm2_ang, & initial_data_direction_nm2, & initial_data_direction_pm2, & - amp_nm2, rl_nm2, ru_nm2, l_ang_nm2, & - amp_pm2, rl_pm2, ru_pm2, l_ang_pm2 + amp_re_nm2, rl_nm2, ru_nm2, l_ang_nm2, & + amp_im_nm2, & + amp_re_pm2, rl_pm2, ru_pm2, l_ang_pm2, & + amp_im_pm2 !============================================================================= implicit none private @@ -41,33 +45,35 @@ subroutine set_initial_data(m_ang,p,q,f) integer(ip) :: i, j integer(ip) :: l_ang real(rp) :: max_val, bump, r, y - real(rp) :: width, amp, rl, ru + real(rp) :: width, rl, ru + + complex(rp) :: amp character(:), allocatable :: initial_data_direction if (m_ang==pm1_ang) then - amp = amp_pm1 + amp = amp_re_pm1 + ZI*amp_im_pm1 ru = ru_pm1 rl = rl_pm1 l_ang = l_ang_pm1 initial_data_direction = initial_data_direction_pm1 else if (m_ang==-pm1_ang) then - amp = amp_nm1 + amp = amp_re_nm1 + ZI*amp_im_nm1 ru = ru_nm1 rl = rl_nm1 l_ang = l_ang_nm1 initial_data_direction = initial_data_direction_nm1 else if (m_ang==pm2_ang) then - amp = amp_pm2 + amp = amp_re_pm2 + ZI*amp_im_pm2 ru = ru_pm2 rl = rl_pm2 l_ang = l_ang_pm2 initial_data_direction = initial_data_direction_pm2 else if (m_ang==-pm2_ang) then - amp = amp_nm2 + amp = amp_re_nm2 + ZI*amp_im_nm2 ru = ru_nm2 rl = rl_nm2 l_ang = l_ang_nm2 diff --git a/src/mod_params.f90 b/src/mod_params.f90 index 37ec212..02c34f0 100644 --- a/src/mod_params.f90 +++ b/src/mod_params.f90 @@ -5,15 +5,15 @@ module mod_params use mod_prec implicit none character(*), parameter :: home_dir = '/home/jripley/teuk-fortran' - character(*), parameter :: run_type = 'multiple_runs' + character(*), parameter :: run_type = 'basic_run' logical, parameter :: debug = .false. real(rp), parameter :: black_hole_mass = 0.5_rp - real(rp), parameter :: black_hole_spin = 0.499_rp + real(rp), parameter :: black_hole_spin = 0.35_rp real(rp), parameter :: compactification_length = 1.0_rp - real(rp), parameter :: evolve_time = 150.0_rp - integer(ip), parameter :: num_saved_times = 500_ip - integer(ip), parameter :: nx = 192_ip - integer(ip), parameter :: nl = 36_ip + real(rp), parameter :: evolve_time = 20.0_rp + integer(ip), parameter :: num_saved_times = 40_ip + integer(ip), parameter :: nx = 64_ip + integer(ip), parameter :: nl = 16_ip logical, parameter :: metric_recon = .true. logical, parameter :: scd_order = .true. logical, parameter :: constrained_evo = .true. @@ -21,53 +21,57 @@ module mod_params logical, parameter :: write_metric_recon_fields = .false. logical, parameter :: write_scd_order_source = .true. logical, parameter :: write_coefs = .false. - character(*), parameter :: computer = 'feynman' + character(*), parameter :: computer = 'home' character(*), parameter :: feyn_out_stem = '/mnt/grtheory/tf-out/' character(*), parameter :: walltime = '168:00:00' character(*), parameter :: memory = '1024' character(*), parameter :: email = 'lloydripley@gmail.com' integer(ip), parameter :: psi_spin = -2_ip integer(ip), parameter :: psi_boost = -2_ip - integer(ip), parameter :: start_multiple = 3_ip + real(rp), parameter :: start_multiple = 1.0_rp integer(ip), parameter :: pm1_ang = 2_ip integer(ip), parameter :: l_ang_pm1 = 2_ip character(*), parameter :: initial_data_direction_pm1 = 'ingoing' - real(rp), parameter :: amp_pm1 = 0.1_rp + real(rp), parameter :: amp_re_pm1 = 0.1_rp + real(rp), parameter :: amp_im_pm1 = 0.1_rp real(rp), parameter :: rl_pm1_0 = -2.0_rp real(rp), parameter :: ru_pm1_0 = 2.0_rp integer(ip), parameter :: l_ang_nm1 = 2_ip character(*), parameter :: initial_data_direction_nm1 = 'ingoing' - real(rp), parameter :: amp_nm1 = 0.1_rp + real(rp), parameter :: amp_re_nm1 = 0.0_rp + real(rp), parameter :: amp_im_nm1 = 0.0_rp real(rp), parameter :: rl_nm1_0 = 1.1_rp real(rp), parameter :: ru_nm1_0 = 3.0_rp integer(ip), parameter :: pm2_ang = 3_ip integer(ip), parameter :: l_ang_pm2 = 3_ip character(*), parameter :: initial_data_direction_pm2 = 'ingoing' - real(rp), parameter :: amp_pm2 = 0.05_rp + real(rp), parameter :: amp_re_pm2 = 0.05_rp + real(rp), parameter :: amp_im_pm2 = 0.05_rp real(rp), parameter :: rl_pm2_0 = -1.5_rp real(rp), parameter :: ru_pm2_0 = 1.5_rp integer(ip), parameter :: l_ang_nm2 = 3_ip character(*), parameter :: initial_data_direction_nm2 = 'ingoing' - real(rp), parameter :: amp_nm2 = 0.0_rp + real(rp), parameter :: amp_re_nm2 = 0.0_rp + real(rp), parameter :: amp_im_nm2 = 0.0_rp real(rp), parameter :: rl_nm2_0 = -1.5_rp real(rp), parameter :: ru_nm2_0 = 1.5_rp integer(ip), dimension(2), parameter :: lin_m = [2_ip,-2_ip] integer(ip), dimension(3), parameter :: scd_m = [0_ip,-4_ip,4_ip] integer(ip), dimension(2), parameter :: lin_write_m = [2_ip,-2_ip] integer(ip), dimension(3), parameter :: scd_write_m = [0_ip,-4_ip,4_ip] - integer(ip), parameter :: max_l = 35_ip - real(rp), parameter :: horizon = 0.5315753542972996_rp - real(rp), parameter :: R_max = 1.8812008343048945_rp - real(rp), parameter :: rl_pm1 = -1.0631507085945993_rp - real(rp), parameter :: ru_pm1 = 1.0631507085945993_rp - real(rp), parameter :: rl_nm1 = 0.5847328897270296_rp - real(rp), parameter :: ru_nm1 = 1.5947260628918989_rp - real(rp), parameter :: rl_pm2 = -0.7973630314459494_rp - real(rp), parameter :: ru_pm2 = 0.7973630314459494_rp - real(rp), parameter :: rl_nm2 = -0.7973630314459494_rp - real(rp), parameter :: ru_nm2 = 0.7973630314459494_rp - real(rp), parameter :: constraint_damping = 447.21348369659154_rp - real(rp), parameter :: scd_order_start_time = 25.94115596715251_rp + integer(ip), parameter :: max_l = 15_ip + real(rp), parameter :: horizon = 0.8567143500057153_rp + real(rp), parameter :: R_max = 1.167250204217227_rp + real(rp), parameter :: rl_pm1 = -1.7134287000114305_rp + real(rp), parameter :: ru_pm1 = 1.7134287000114305_rp + real(rp), parameter :: rl_nm1 = 0.9423857850062869_rp + real(rp), parameter :: ru_nm1 = 2.570143050017146_rp + real(rp), parameter :: rl_pm2 = -1.285071525008573_rp + real(rp), parameter :: ru_pm2 = 1.285071525008573_rp + real(rp), parameter :: rl_nm2 = -1.285071525008573_rp + real(rp), parameter :: ru_nm2 = 1.285071525008573_rp + real(rp), parameter :: constraint_damping = 36.51483710615301_rp + real(rp), parameter :: scd_order_start_time = 11.248163954718162_rp integer(ip), dimension(1), parameter :: lin_pos_m = [2_ip] integer(ip), parameter :: len_lin_pos_m = 1_ip integer(ip), parameter :: len_lin_m = 2_ip @@ -75,17 +79,16 @@ module mod_params integer(ip), parameter :: len_lin_write_m = 2_ip integer(ip), parameter :: len_scd_write_m = 3_ip integer(ip), parameter :: num_threads = 1_ip - integer(ip), parameter :: ny = 50_ip - real(rp), parameter :: dt = 0.000244140625_rp - integer(ip), parameter :: nt = 307200_ip - integer(ip), parameter :: t_step_save = 614_ip + integer(ip), parameter :: ny = 30_ip + real(rp), parameter :: dt = 0.002197265625_rp + integer(ip), parameter :: nt = 4551_ip + integer(ip), parameter :: t_step_save = 113_ip integer(ip), parameter :: max_m = 4_ip integer(ip), parameter :: min_m = -4_ip integer(ip), parameter :: max_s = 3_ip integer(ip), parameter :: min_s = -3_ip - character(*), parameter :: output_stem = 'Sep_12_18_10_39_m0.5_s0.499_nx192_nl36' - character(*), parameter :: output_dir = '/mnt/grtheory/tf-out/Sep_12_18_10_39_m0.5_s0.499_nx192_nl36' - character(*), parameter :: bin_name = 'Sep_12_18_10_39_m0.5_s0.499_nx192_nl36.run' - character(*), parameter :: tables_dir = '/mnt/grtheory/tf-out/Sep_12_18_10_39_m0.5_s0.499_nx192_nl36/tables' - character(*), parameter :: output_file = '/mnt/grtheory/tf-out/Sep_12_18_09_44_m0.5_s0.499_nx192_nl36/output.txt' + character(*), parameter :: output_stem = 'Sep_12_15_38_40_m0.5_s0.35_nx64_nl16' + character(*), parameter :: output_dir = 'output/Sep_12_15_38_40_m0.5_s0.35_nx64_nl16' + character(*), parameter :: bin_name = 'Sep_12_15_38_40_m0.5_s0.35_nx64_nl16.run' + character(*), parameter :: tables_dir = 'output/Sep_12_15_38_40_m0.5_s0.35_nx64_nl16/tables' end module mod_params