...
 
Commits (5)
......@@ -228,6 +228,16 @@ def angular_momentum(x, q, m, chi1, chi2, LInitNR):
* x**3.))
return l - LInitNR
# Convert Cacus truth values to python's
def to_bool(s):
s = tolower(s)
if s in ["true", "yes", "1"]:
return True
elif s in ["false", "no", "0"]:
return False
else:
raise(ValueError("Not a boolean values: %s" % s))
#Get cutoff frequency
def getCutoffFrequency(sim_name):
"""
......@@ -236,11 +246,20 @@ def getCutoffFrequency(sim_name):
sim_name = string of simulation
return = cutoff frequency
"""
center_offset = 0.
par_b = 1.0
pyp = 0.
pym = 0.
target_m1 = 1.0
target_m2 = 1.0
give_bare_mass = True
S1 = 0.
S2 = 0.
filename = main_dir+"/output-0000/%s.par" % (sim_name)
with open(filename) as file:
contents = file.readlines()
for line in contents:
line_elems = line.split(" ")
line_elems = re.sub("#.*", "", line).split(" ")
if(line_elems[0] == "TwoPunctures::par_b"):
par_b = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::center_offset[0]"):
......@@ -250,21 +269,55 @@ def getCutoffFrequency(sim_name):
if(line_elems[0] == "TwoPunctures::par_P_minus[1]"):
pym = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::target_M_plus"):
m1 = float(line_elems[-1])
target_m1 = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::target_M_minus"):
m2 = float(line_elems[-1])
target_m2 = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::par_S_plus[2]"):
S1 = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::par_S_minus[2]"):
S2 = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::give_bare_mass"):
give_bare_mass = to_bool(line_elems[-1])
adm_m1 = None
adm_m2 = None
two_punctures_files = sorted(glob.glob(main_dir+"/output-????/%s/TwoPunctures.bbh" % (sim)))
out_files = sorted(glob.glob(main_dir+"/output-????/%s.out" % (sim)))
par_files = sorted(glob.glob(main_dir+"/output-????/%s.par" % (sim)))
if(two_punctures_files):
two_punctures_file = two_punctures_files[0]
with open(two_punctures_file) as file:
contents = file.readlines()
for line in contents:
line_elems = re.sub("#.*", "", line).split(" ")
if(line_elems[0] == "initial-bh-puncture-adm-mass1"):
adm_m1 = float(line_elems[-1])
if(line_elems[0] == "initial-bh-puncture-adm-mass2"):
adm_m1 = float(line_elems[-1])
elif(out_files):
out_file = out_files[0]
with open(out_file) as file:
contents = file.readlines()
for line in contents:
m = re.match("INFO \(TwoPunctures\): Puncture 1 ADM mass is (.*)", line)
if(m):
adm_m1 = float(m.group(1))
m = re.match("INFO \(TwoPunctures\): Puncture 2 ADM mass is (.*)", line)
if(m):
adm_m2 = float(m.group(1))
elif(not give_bare_mass):
adm_m1 = target_m1
adm_m2 = target_m2
else:
print("Cannot determine ADM masses of punctures")
raise ValueError
xp = par_b + center_offset
xm = -1*par_b + center_offset
LInitNR = xp*pyp + xm*pym
M = m1+m2
q = m1/m2
chi1 = S1/m1**2
chi2 = S2/m2**2
M = adm_m1+adm_m2
q = adm_m1/adm_m2
chi1 = S1/adm_m1**2
chi2 = S2/adm_m2**2
# .014 is the initial guess for cutoff frequency
omOrbPN = scipy.optimize.fsolve(angular_momentum, .014, (q, M, chi1, chi2, LInitNR))[0]
omOrbPN = omOrbPN**(3./2.)
......@@ -432,7 +485,7 @@ if __name__ == "__main__":
# TODO: fix this. It will fail if output-0000 does not contain any mp
# output and also will open the output files multiple times
fn = sorted(glob.glob(simdirs+"mp_psi4.h5"))[0]
fn = sorted(glob.glob(simdirs+"mp_[Pp]si4.h5"))[0]
with h5py.File(fn, "r") as fh:
# get all radii
radii = set()
......@@ -456,7 +509,7 @@ if __name__ == "__main__":
mp_psi4_vars = []
for radius in radii:
psi4dsetname = dsets[(radius, (l,m))]
mp_psi4 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname)
mp_psi4 = loadHDF5Series(simdirs+"mp_[Pp]si4.h5", psi4dsetname)
mp_psi4_vars.append(mp_psi4)
#Get Tortoise Coordinate
......