Commit 6fcda20a authored by Roland Haas's avatar Roland Haas
Browse files

POWER: make usable as python module

parent 7b20453a
Pipeline #667 failed with stages
...@@ -371,193 +371,194 @@ def get_angular_momentum(python_strain): ...@@ -371,193 +371,194 @@ def get_angular_momentum(python_strain):
#-----Main-----# #-----Main-----#
#Initialize simulation data if __name__ == "__main__":
if(len(sys.argv) < 2): #Initialize simulation data
print("Pass in the number n of the n innermost detector radii to be used in the extrapolation (optional, default=all) and the simulation folders (e.g., ./power.py 6 ./simulations/J0040_N40 /path/to/my_simulation_folder).") if(len(sys.argv) < 2):
sys.exit() print("Pass in the number n of the n innermost detector radii to be used in the extrapolation (optional, default=all) and the simulation folders (e.g., ./power.py 6 ./simulations/J0040_N40 /path/to/my_simulation_folder).")
elif(os.path.isdir(sys.argv[1])): sys.exit()
radiiUsedForExtrapolation = 7 #use the first n radii available elif(os.path.isdir(sys.argv[1])):
paths = sys.argv[1:] radiiUsedForExtrapolation = 7 #use the first n radii available
elif(not os.path.isdir(sys.argv[1])): paths = sys.argv[1:]
radiiUsedForExtrapolation = int(sys.argv[1]) #use the first n radii available elif(not os.path.isdir(sys.argv[1])):
if(radiiUsedForExtrapolation < 1 or radiiUsedForExtrapolation > 7): radiiUsedForExtrapolation = int(sys.argv[1]) #use the first n radii available
print("Invalid specified radii number") if(radiiUsedForExtrapolation < 1 or radiiUsedForExtrapolation > 7):
sys.exit() print("Invalid specified radii number")
paths = sys.argv[2:] sys.exit()
paths = sys.argv[2:]
for sim_path in paths:
main_dir = sim_path for sim_path in paths:
sim = sim_path.split("/")[-1] main_dir = sim_path
sim = sim_path.split("/")[-1]
simdirs = main_dir+"/output-????/%s/" % (sim)
f0 = getCutoffFrequency(sim) simdirs = main_dir+"/output-????/%s/" % (sim)
f0 = getCutoffFrequency(sim)
#Check if necessary files exist
par_file = main_dir+"/output-0000/%s.par" % (sim) #Check if necessary files exist
two_punctures_file = main_dir+"/output-0000/%s/TwoPunctures.bbh" % (sim) par_file = main_dir+"/output-0000/%s.par" % (sim)
if(not os.path.isfile(par_file) or not os.path.isfile(two_punctures_file)): two_punctures_file = main_dir+"/output-0000/%s/TwoPunctures.bbh" % (sim)
continue if(not os.path.isfile(par_file) or not os.path.isfile(two_punctures_file)):
continue
#Create data directories
main_directory = "Extrapolated_Strain" #Create data directories
sim_dir = main_directory+"/"+sim main_directory = "Extrapolated_Strain"
if not os.path.exists(main_directory): sim_dir = main_directory+"/"+sim
os.makedirs(main_directory) if not os.path.exists(main_directory):
if not os.path.exists(sim_dir): os.makedirs(main_directory)
os.makedirs(sim_dir) if not os.path.exists(sim_dir):
os.makedirs(sim_dir)
#Get ADMMass
ADMMass = -1 #Get ADMMass
filename = main_dir+"/output-0000/%s/TwoPunctures.bbh" % (sim) ADMMass = -1
with open(filename) as file: filename = main_dir+"/output-0000/%s/TwoPunctures.bbh" % (sim)
contents = file.readlines() with open(filename) as file:
for line in contents: contents = file.readlines()
line_elems = line.split(" ") for line in contents:
if(line_elems[0] == "initial-ADM-energy"): line_elems = line.split(" ")
ADMMass = float(line_elems[-1]) if(line_elems[0] == "initial-ADM-energy"):
ADMMass = float(line_elems[-1])
# TODO: fix this. It will fail if output-0000 does not contain any mp
# output and also will open the output files multiple times # TODO: fix this. It will fail if output-0000 does not contain any mp
fn = sorted(glob.glob(simdirs+"mp_psi4.h5"))[0] # output and also will open the output files multiple times
with h5py.File(fn, "r") as fh: fn = sorted(glob.glob(simdirs+"mp_psi4.h5"))[0]
# get all radii with h5py.File(fn, "r") as fh:
radii = set() # get all radii
modes = set() radii = set()
dsets = dict() modes = set()
for dset in fh: dsets = dict()
# TODO: extend Multipole to save the radii as attributes and/or for dset in fh:
# use a group structure in the hdf5 file # TODO: extend Multipole to save the radii as attributes and/or
m = re.match(r'l(\d*)_m(-?\d*)_r(\d*\.\d)', dset) # use a group structure in the hdf5 file
if m: m = re.match(r'l(\d*)_m(-?\d*)_r(\d*\.\d)', dset)
radius = float(m.group(3)) if m:
mode = (int(m.group(1)), int(m.group(2))) radius = float(m.group(3))
modes.add(mode) mode = (int(m.group(1)), int(m.group(2)))
radii.add(radius) modes.add(mode)
dsets[(radius, mode)] = dset radii.add(radius)
modes = sorted(modes) dsets[(radius, mode)] = dset
radii = sorted(radii) modes = sorted(modes)
radii = sorted(radii)
#Get Psi4
for (l,m) in modes: #Get Psi4
mp_psi4_vars = [] for (l,m) in modes:
for radius in radii: mp_psi4_vars = []
psi4dsetname = dsets[(radius, (l,m))] for radius in radii:
mp_psi4 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname) psi4dsetname = dsets[(radius, (l,m))]
mp_psi4_vars.append(mp_psi4) mp_psi4 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname)
mp_psi4_vars.append(mp_psi4)
#Get Tortoise Coordinate
tortoise = [] #Get Tortoise Coordinate
for radius in radii: tortoise = []
tortoise.append(-RadialToTortoise(radius, ADMMass)) for radius in radii:
tortoise.append(-RadialToTortoise(radius, ADMMass))
strain = []
phase = [] strain = []
amp = [] phase = []
for i in range(len(radii)): amp = []
#Get modified Psi4 (Multiply real and imaginary psi4 columns by radii and add the tortoise coordinate to the time colum) for i in range(len(radii)):
mp_psi4_vars[i][:, 0] += tortoise[i] #Get modified Psi4 (Multiply real and imaginary psi4 columns by radii and add the tortoise coordinate to the time colum)
mp_psi4_vars[i][:, 1] *= radii[i] mp_psi4_vars[i][:, 0] += tortoise[i]
mp_psi4_vars[i][:, 2] *= radii[i] mp_psi4_vars[i][:, 1] *= radii[i]
mp_psi4_vars[i][:, 2] *= radii[i]
#Check for psi4 amplitude going to zero
cur_psi4_amp = (mp_psi4_vars[i][0, 1]**2 + mp_psi4_vars[i][0, 2]**2)**(1/2) #Check for psi4 amplitude going to zero
min_psi4_amp = cur_psi4_amp cur_psi4_amp = (mp_psi4_vars[i][0, 1]**2 + mp_psi4_vars[i][0, 2]**2)**(1/2)
# TODO: use array notatino for this since it finds the minimum amplitude min_psi4_amp = cur_psi4_amp
for j in range(0, len(mp_psi4_vars[i][:, 0])): # TODO: use array notatino for this since it finds the minimum amplitude
cur_psi4_amp = (mp_psi4_vars[i][j, 1]**2 + mp_psi4_vars[i][j, 2]**2)**(1/2) for j in range(0, len(mp_psi4_vars[i][:, 0])):
if(cur_psi4_amp < min_psi4_amp): cur_psi4_amp = (mp_psi4_vars[i][j, 1]**2 + mp_psi4_vars[i][j, 2]**2)**(1/2)
min_psi4_amp = cur_psi4_amp if(cur_psi4_amp < min_psi4_amp):
if(min_psi4_amp < np.finfo(float).eps and l >= 2): min_psi4_amp = cur_psi4_amp
print("The psi4 amplitude is near zero. The phase is ill-defined.") if(min_psi4_amp < np.finfo(float).eps and l >= 2):
print("The psi4 amplitude is near zero. The phase is ill-defined.")
#Fixed-frequency integration twice to get strain
hTable = psi4ToStrain(mp_psi4_vars[i], f0) #Fixed-frequency integration twice to get strain
time = hTable[:, 0] hTable = psi4ToStrain(mp_psi4_vars[i], f0)
h = hTable[:, 1] time = hTable[:, 0]
hplus = h.real h = hTable[:, 1]
hcross = h.imag hplus = h.real
newhTable = np.column_stack((time, hplus, hcross)) hcross = h.imag
warnings.filterwarnings('ignore') newhTable = np.column_stack((time, hplus, hcross))
finalhTable = newhTable.astype(float) warnings.filterwarnings('ignore')
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_strain_at_"+str(radii[i])+"_l"+str(l)+"_m"+str(m)+".dat", finalhTable) finalhTable = newhTable.astype(float)
strain.append(finalhTable) np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_strain_at_"+str(radii[i])+"_l"+str(l)+"_m"+str(m)+".dat", finalhTable)
strain.append(finalhTable)
#Get phase and amplitude of strain
h_phase = np.unwrap(np.angle(h)) #Get phase and amplitude of strain
while(h_phase[0] < 0): h_phase = np.unwrap(np.angle(h))
h_phase[:] += 2*math.pi while(h_phase[0] < 0):
angleTable = np.column_stack((time, h_phase)) h_phase[:] += 2*math.pi
angleTable = angleTable.astype(float) angleTable = np.column_stack((time, h_phase))
phase.append(angleTable) angleTable = angleTable.astype(float)
h_amp = np.absolute(h) phase.append(angleTable)
ampTable = np.column_stack((time, h_amp)) h_amp = np.absolute(h)
ampTable = ampTable.astype(float) ampTable = np.column_stack((time, h_amp))
amp.append(ampTable) ampTable = ampTable.astype(float)
amp.append(ampTable)
#Interpolate phase and amplitude
t = phase[0][:, 0] #Interpolate phase and amplitude
last_t = phase[radiiUsedForExtrapolation - 1][-1, 0] t = phase[0][:, 0]
last_index = 0; last_t = phase[radiiUsedForExtrapolation - 1][-1, 0]
# TODO: use array notation for this (this is a boolean last_index = 0;
# plus a first_of or so) # TODO: use array notation for this (this is a boolean
for i in range(0, len(phase[0][:, 0])): # plus a first_of or so)
if(t[i] > last_t): for i in range(0, len(phase[0][:, 0])):
last_index = i if(t[i] > last_t):
break last_index = i
last_index = last_index-1 break
t = phase[0][0:last_index, 0] last_index = last_index-1
dts = t[1:] - t[:-1] t = phase[0][0:last_index, 0]
dt = float(np.amin(dts)) dts = t[1:] - t[:-1]
t = np.arange(phase[0][0, 0], phase[0][last_index, 0], dt) dt = float(np.amin(dts))
interpolation_order = 9 t = np.arange(phase[0][0, 0], phase[0][last_index, 0], dt)
for i in range(0, radiiUsedForExtrapolation): interpolation_order = 9
interp_function = scipy.interpolate.interp1d(phase[i][:, 0], phase[i][:, 1], kind=interpolation_order) for i in range(0, radiiUsedForExtrapolation):
resampled_phase_vals = interp_function(t) interp_function = scipy.interpolate.interp1d(phase[i][:, 0], phase[i][:, 1], kind=interpolation_order)
while(resampled_phase_vals[0] < 0): resampled_phase_vals = interp_function(t)
resampled_phase_vals[:] += 2*math.pi while(resampled_phase_vals[0] < 0):
phase[i] = np.column_stack((t, resampled_phase_vals)) resampled_phase_vals[:] += 2*math.pi
interp_function = scipy.interpolate.interp1d(amp[i][:, 0], amp[i][:, 1], kind=interpolation_order) phase[i] = np.column_stack((t, resampled_phase_vals))
resampled_amp_vals = interp_function(t) interp_function = scipy.interpolate.interp1d(amp[i][:, 0], amp[i][:, 1], kind=interpolation_order)
amp[i] = np.column_stack((t, resampled_amp_vals)) resampled_amp_vals = interp_function(t)
amp[i] = np.column_stack((t, resampled_amp_vals))
#Extrapolate
phase_extrapolation_order = 1 #Extrapolate
amp_extrapolation_order = 2 phase_extrapolation_order = 1
radii = np.asarray(radii, dtype=float) amp_extrapolation_order = 2
radii = radii[0:radiiUsedForExtrapolation] radii = np.asarray(radii, dtype=float)
# TODO: replace by np.ones (which is all it does anyway) radii = radii[0:radiiUsedForExtrapolation]
A_phase = np.power(radii, 0) # TODO: replace by np.ones (which is all it does anyway)
A_amp = np.power(radii, 0) A_phase = np.power(radii, 0)
for i in range(1, phase_extrapolation_order+1): A_amp = np.power(radii, 0)
A_phase = np.column_stack((A_phase, np.power(radii, -1*i))) for i in range(1, phase_extrapolation_order+1):
A_phase = np.column_stack((A_phase, np.power(radii, -1*i)))
for i in range(1, amp_extrapolation_order+1):
A_amp = np.column_stack((A_amp, np.power(radii, -1*i))) for i in range(1, amp_extrapolation_order+1):
A_amp = np.column_stack((A_amp, np.power(radii, -1*i)))
radially_extrapolated_phase = np.empty(0)
radially_extrapolated_amp = np.empty(0) radially_extrapolated_phase = np.empty(0)
for i in range(0, len(t)): radially_extrapolated_amp = np.empty(0)
b_phase = np.empty(0) for i in range(0, len(t)):
for j in range(0, radiiUsedForExtrapolation): b_phase = np.empty(0)
b_phase = np.append(b_phase, phase[j][i, 1]) for j in range(0, radiiUsedForExtrapolation):
x_phase = np.linalg.lstsq(A_phase, b_phase)[0] b_phase = np.append(b_phase, phase[j][i, 1])
radially_extrapolated_phase = np.append(radially_extrapolated_phase, x_phase[0]) x_phase = np.linalg.lstsq(A_phase, b_phase)[0]
radially_extrapolated_phase = np.append(radially_extrapolated_phase, x_phase[0])
b_amp = np.empty(0)
for j in range(0, radiiUsedForExtrapolation): b_amp = np.empty(0)
b_amp = np.append(b_amp, amp[j][i, 1]) for j in range(0, radiiUsedForExtrapolation):
x_amp = np.linalg.lstsq(A_amp, b_amp)[0] b_amp = np.append(b_amp, amp[j][i, 1])
radially_extrapolated_amp = np.append(radially_extrapolated_amp, x_amp[0]) x_amp = np.linalg.lstsq(A_amp, b_amp)[0]
radially_extrapolated_amp = np.append(radially_extrapolated_amp, x_amp[0])
radially_extrapolated_h_plus = np.empty(0)
radially_extrapolated_h_cross = np.empty(0) radially_extrapolated_h_plus = np.empty(0)
for i in range(0, len(radially_extrapolated_amp)): radially_extrapolated_h_cross = np.empty(0)
radially_extrapolated_h_plus = np.append(radially_extrapolated_h_plus, radially_extrapolated_amp[i] * math.cos(radially_extrapolated_phase[i])) for i in range(0, len(radially_extrapolated_amp)):
radially_extrapolated_h_cross = np.append(radially_extrapolated_h_cross, radially_extrapolated_amp[i] * math.sin(radially_extrapolated_phase[i])) radially_extrapolated_h_plus = np.append(radially_extrapolated_h_plus, radially_extrapolated_amp[i] * math.cos(radially_extrapolated_phase[i]))
radially_extrapolated_h_cross = np.append(radially_extrapolated_h_cross, radially_extrapolated_amp[i] * math.sin(radially_extrapolated_phase[i]))
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_strain_l"+str(l)+"_m"+str(m)+".dat", np.column_stack((t, radially_extrapolated_h_plus, radially_extrapolated_h_cross)))
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_amplitude_l"+str(l)+"_m"+str(m)+".dat", np.column_stack((t, radially_extrapolated_amp))) np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_strain_l"+str(l)+"_m"+str(m)+".dat", np.column_stack((t, radially_extrapolated_h_plus, radially_extrapolated_h_cross)))
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_phase_l"+str(l)+"_m"+str(m)+".dat", np.column_stack((t, radially_extrapolated_phase[:]))) np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_amplitude_l"+str(l)+"_m"+str(m)+".dat", np.column_stack((t, radially_extrapolated_amp)))
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_phase_l"+str(l)+"_m"+str(m)+".dat", np.column_stack((t, radially_extrapolated_phase[:])))
get_energy(sim)
get_angular_momentum(sim) get_energy(sim)
get_angular_momentum(sim)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment