diff --git a/README.md b/README.md index 7a6ca47..69eb053 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ A Fortran ('08) and python code that solves Teukolsky equation for the linearly perturbed Newman-Penrose scalar Psi4 about a Kerr black hole. -The code also reconstructs +The code also directly reconstructs the linear spacetime metric in outgoing raditation gauge from the linearized Newman-Penrose scalar Psi\_4, and then solves the equations of motion for the second order Psi\_4. @@ -21,6 +21,9 @@ Runtime parameters are configured in the `setup.py` file. * Fastest Fourier Transform in the West (FFTW): http://fftw.org +* OpenMP (this is optional, and can be deactivated in the Makefile): + https://www.openmp.org/ + I have successfully compiled the code with gfortran (version 9) and ifort (version 17). diff --git a/setup.py b/setup.py index c6d5fda..e446db9 100755 --- a/setup.py +++ b/setup.py @@ -13,14 +13,14 @@ sim= Sim(args) #============================================================================= sim.black_hole_mass= float(0.5) -sim.black_hole_spin= round(0.7*sim.black_hole_mass,4) +sim.black_hole_spin= round(0.998*sim.black_hole_mass,4) sim.compactification_length= float(1) #============================================================================= -sim.evolve_time= float(20) ## units of black hole mass -sim.num_saved_times= int(40) +sim.evolve_time= float(150) ## units of black hole mass +sim.num_saved_times= int(500) #============================================================================= -sim.nx= 64 ## num radial pts -sim.nl= 16 ## num angular values +sim.nx= 192 ## num radial pts +sim.nl= 36 ## 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= 'home'#'feynman'# +sim.computer= 'feynman'#'home'# sim.feyn_out_stem= '/mnt/grtheory/tf-out/' ## for feynman cluster/slurm script @@ -74,10 +74,10 @@ sim.initial_data_direction_pm1= "ingoing"#"outgoing"#"time_symmetric"# sim.amp_re_pm1= float(0.1) -sim.amp_im_pm1= float(0.1) +sim.amp_im_pm1= float(0.0) -sim.rl_pm1_0= float(-2.0) -sim.ru_pm1_0= float( 2.0) +sim.rl_pm1_0= float( 1.1) +sim.ru_pm1_0= float( 2.5) #----------------------------------------------------------------------------- sim.l_ang_nm1= int(2) @@ -87,7 +87,7 @@ sim.amp_im_nm1= float( 0.0) sim.rl_nm1_0= float( 1.1) -sim.ru_nm1_0= float( 3.0) +sim.ru_nm1_0= float( 2.5) #============================================================================= ## initial data for mode m2 #============================================================================= @@ -97,8 +97,8 @@ sim.initial_data_direction_pm2= "ingoing"#"time_symmetric"#"outgoing"# -sim.amp_re_pm2= float(0.05) -sim.amp_im_pm2= float(0.05) +sim.amp_re_pm2= float(0.0) +sim.amp_im_pm2= float(0.0) sim.rl_pm2_0= float(-1.5) sim.ru_pm2_0= float( 1.5) @@ -145,28 +145,28 @@ #============================================================================= elif (sim.run_type == "multiple_runs"): #----------------------------------------------------------------------------- - default_amp = 0.1 default_sm = 1 default_nx_07 = 176 default_nl_07 = 32 - - default_nx_099 = 132 - default_nl_099 = 28 +# default_nx_07 = 144 +# default_nl_07 = 24 default_nx_0998 = 192 default_nl_0998 = 36 +# default_nx_0998 = 160 +# default_nl_0998 = 28 #----------------------------------------------------------------------------- #----------------------------------------------------------------------------- sim.black_hole_spin= round(0.7*sim.black_hole_mass,4) sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp sim.nx = default_nx_07 sim.nl = default_nl_07 nxs = [160, 176, 192] nls = [ 28, 32, 36] +# nxs = [128, 144, 160] +# nls = [ 20, 24, 28] for i in range(len(nxs)): sim.nx = nxs[i] sim.nl = nls[i] @@ -174,21 +174,6 @@ #----------------------------------------------------------------------------- sim.black_hole_spin= round(0.7*sim.black_hole_mass,4) sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp - sim.nx = default_nx_07 - sim.nl = default_nl_07 - - amps = [0.3, 0.5] - for amp in amps: - sim.amp_pm1= amp - sim.amp_nm1= amp - sim.launch_run() -#----------------------------------------------------------------------------- - sim.black_hole_spin= round(0.7*sim.black_hole_mass,4) - sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp sim.nx = default_nx_07 sim.nl = default_nl_07 @@ -200,13 +185,13 @@ #----------------------------------------------------------------------------- sim.black_hole_spin= round(0.998*sim.black_hole_mass,4) sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp sim.nx = default_nx_0998 sim.nl = default_nl_0998 nxs = [176, 192, 208] nls = [32, 36, 40] +# nxs = [144, 160, 176] +# nls = [ 24, 28, 32] for i in range(len(nxs)): sim.nx = nxs[i] sim.nl = nls[i] @@ -214,21 +199,6 @@ #----------------------------------------------------------------------------- sim.black_hole_spin= round(0.998*sim.black_hole_mass,4) sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp - sim.nx = default_nx_0998 - sim.nl = default_nl_0998 - - amps = [0.3, 0.5] - for amp in amps: - sim.amp_pm1= amp - sim.amp_nm1= amp - sim.launch_run() -#----------------------------------------------------------------------------- - sim.black_hole_spin= round(0.998*sim.black_hole_mass,4) - sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp sim.nx = default_nx_0998 sim.nl = default_nl_0998 @@ -236,55 +206,6 @@ for sm in sms: sim.start_multiple= sm sim.launch_run() - - sys.exit() -#----------------------------------------------------------------------------- -#----------------------------------------------------------------------------- - sim.black_hole_spin= round(0.99*sim.black_hole_mass,4) - sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp - sim.nx = default_nx_099 - sim.nl = default_nl_099 - - nxs = [96, 112, 128] - nls = [20, 24, 28] - for i in range(len(nxs)): - sim.nx = nxs[i] - sim.nl = nls[i] - sim.launch_run() -#----------------------------------------------------------------------------- - sim.black_hole_spin= round(0.99*sim.black_hole_mass,4) - sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp - sim.nx = default_nx_099 - sim.nl = default_nl_099 - - amps = [0.3, 0.5] - for amp in amps: - sim.amp_pm1= amp - sim.amp_nm1= amp - sim.launch_run() -#----------------------------------------------------------------------------- - sim.black_hole_spin= round(0.99*sim.black_hole_mass,4) - sim.start_multiple = default_sm - sim.amp_pm1 = default_amp - sim.amp_nm1 = default_amp - sim.nx = default_nx_099 - sim.nl = default_nl_099 - - sms = [2,3] - for sm in sms: - sim.start_multiple= sm - sim.launch_run() -#============================================================================= -elif (sim.run_type == "start_times"): - sms = [1,2,3] - - for sm in sms: - sim.start_multiple= sm - sim.launch_run() #============================================================================= elif (sim.run_type == "spin_ramp"): for bhs in [0,0.01,0.02,0.04,0.08,0.12,0.16,0.2,0.24,0.28,0.32]: diff --git a/src/mod_params.f90 b/src/mod_params.f90 index 02c34f0..dede0c1 100644 --- a/src/mod_params.f90 +++ b/src/mod_params.f90 @@ -8,12 +8,12 @@ module mod_params 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.35_rp + real(rp), parameter :: black_hole_spin = 0.499_rp real(rp), parameter :: compactification_length = 1.0_rp - 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 + 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 logical, parameter :: metric_recon = .true. logical, parameter :: scd_order = .true. logical, parameter :: constrained_evo = .true. @@ -21,7 +21,7 @@ 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 = 'home' + character(*), parameter :: computer = 'feynman' character(*), parameter :: feyn_out_stem = '/mnt/grtheory/tf-out/' character(*), parameter :: walltime = '168:00:00' character(*), parameter :: memory = '1024' @@ -33,20 +33,20 @@ module mod_params integer(ip), parameter :: l_ang_pm1 = 2_ip character(*), parameter :: initial_data_direction_pm1 = 'ingoing' 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 + real(rp), parameter :: amp_im_pm1 = 0.0_rp + real(rp), parameter :: rl_pm1_0 = 1.1_rp + real(rp), parameter :: ru_pm1_0 = 2.5_rp integer(ip), parameter :: l_ang_nm1 = 2_ip character(*), parameter :: initial_data_direction_nm1 = 'ingoing' 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 + real(rp), parameter :: ru_nm1_0 = 2.5_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_re_pm2 = 0.05_rp - real(rp), parameter :: amp_im_pm2 = 0.05_rp + real(rp), parameter :: amp_re_pm2 = 0.0_rp + real(rp), parameter :: amp_im_pm2 = 0.0_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 @@ -59,19 +59,19 @@ module mod_params 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 = 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), 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 = 0.5847328897270296_rp + real(rp), parameter :: ru_pm1 = 1.328938385743249_rp + real(rp), parameter :: rl_nm1 = 0.5847328897270296_rp + real(rp), parameter :: ru_nm1 = 1.328938385743249_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 = 6.854615053280418_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 @@ -79,16 +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 = 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 :: 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 :: 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_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' + character(*), parameter :: output_stem = 'Sep_15_19_22_37_m0.5_s0.499_nx192_nl36' + character(*), parameter :: output_dir = '/mnt/grtheory/tf-out/Sep_15_19_22_37_m0.5_s0.499_nx192_nl36' + character(*), parameter :: bin_name = 'Sep_15_19_22_37_m0.5_s0.499_nx192_nl36.run' + character(*), parameter :: tables_dir = '/mnt/grtheory/tf-out/Sep_15_19_22_37_m0.5_s0.499_nx192_nl36/tables' end module mod_params