Forked from
Eliu Huerta Escudero / Waveform_Extractor
116 commits behind the upstream repository.
-
Daniel Johnson authoredDaniel Johnson authored
power.py 15.43 KiB
#!/usr/bin/env python
# Copyright (c) 2017 The Board of Trustees of the University of Illinois
# All rights reserved.
#
# Developed by: Daniel Johnson, E. A. Huerta, Roland Haas
# NCSA Gravity Group
# National Center for Supercomputing Applications
# University of Illinois at Urbana-Champaign
# http://gravity.ncsa.illinois.edu/
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal with the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# Redistributions of source code must retain the above copyright notice, this
# list of conditions and the following disclaimers.
#
# Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimers in the documentation
# and/or other materials provided with the distribution.
#
# Neither the names of the National Center for Supercomputing Applications,
# University of Illinois at Urbana-Champaign, nor the names of its
# contributors may be used to endorse or promote products derived from this
# Software without specific prior written permission.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# WITH THE SOFTWARE.
# Based off of SimulationTools Mathematica Package
# http://www.simulationtools.org/
import numpy as np
import glob
import os
import h5py
import string
import math
import sys
import warnings
import scipy.optimize
import scipy.interpolate
#-----Function Definitions-----#
#Function used in getting psi4 from simulation
def joinDsets(dsets):
"""joints multiple datasets which each have a
time like first column, eg iteration number of
time. Removes overlapping segments, keeping the
last segment.
dsets = iterable of 2d array like objects with data"""
# joins multiple datasets of which the first column is assumed to be "time"
if(not dsets):
return None
length = 0
for d in dsets:
length += len(d)
newshape = list(dsets[0].shape)
newshape[0] = length
dset = np.empty(shape=newshape, dtype=dsets[0].dtype)
usedlength = 0
for d in dsets:
insertpointidx = np.where(dset[0:usedlength,0] >= d[0,0])
if(insertpointidx[0].size):
insertpoint = insertpointidx[0][0]
else:
insertpoint = usedlength
newlength = insertpoint+len(d)
dset[insertpoint:newlength] = d
usedlength = newlength
return dset[0:usedlength]
#Function used in getting psi4 from simulation
def loadHDF5Series(nameglob, series):
"""load HDF5 timeseries data and concatenate the content of multiple files
nameglob = a shell glob that matches all files to be loaded,
files are sorted alphabetically
series = HDF5 dataset name of dataset to load from files"""
dsets = list()
for fn in sorted(glob.glob(nameglob)):
fh = h5py.File(fn)
dsets.append(fh[series])
return joinDsets(dsets)
#Convert radial to tortoise coordinates
def RadialToTortoise(r, M):
"""
Convert the radial coordinate to the tortoise coordinate
r = radial coordinate
M = ADMMass used to convert coordinate
return = tortoise coordinate value
"""
return r + 2. * M * math.log( r / (2. * M) - 1.)
#Convert modified psi4 to strain
def psi4ToStrain(mp_psi4, f0):
"""
Convert the input mp_psi4 data to the strain of the gravitational wave
mp_psi4 = Weyl scalar result from simulation
f0 = cutoff frequency
return = strain (h) of the gravitational wave
"""
#TODO: Check for uniform spacing in time
t0 = mp_psi4[:, 0]
list_len = len(t0)
complexPsi = np.zeros(list_len, dtype=np.complex_)
complexPsi = mp_psi4[:, 1]+1.j*mp_psi4[:, 2]
freq, psif = myFourierTransform(t0, complexPsi)
dhf = ffi(freq, psif, f0)
hf = ffi(freq, dhf, f0)
time, h = myFourierTransformInverse(freq, hf, t0[0])
hTable = np.column_stack((time, h))
return hTable
#Fixed frequency integration
# See https://arxiv.org/abs/1508.07250 for method
def ffi(freq, data, f0):
"""
Integrates the data according to the input frequency and cutoff frequency
freq = fourier transform frequency
data = input on which ffi is performed
f0 = cutoff frequency
"""
f1 = f0/(2*math.pi)
fs = freq
gs = data
mask1 = (np.sign((fs/f1) - 1) + 1)/2
mask2 = (np.sign((-fs/f1) - 1) + 1)/2
mask = 1 - (1 - mask1) * (1 - mask2)
fs2 = mask * fs + (1-mask) * f1 * np.sign(fs - np.finfo(float).eps)
new_gs = gs/(2*math.pi*1.j*fs2)
return new_gs
#Fourier Transform
def myFourierTransform(t0, complexPsi):
"""
Transforms the complexPsi data to frequency space
t0 = time data points
complexPsi = data points of Psi to be transformed
"""
psif = np.fft.fft(complexPsi, norm="ortho")
l = len(complexPsi)
n = int(math.floor(l/2))
newpsif = psif[l-n:]
newpsif = np.append(newpsif, psif[:l-n])
T = np.amin(np.diff(t0))*l
freq = range(-n, l-n)/T
return freq, newpsif
#Inverse Fourier Transform
def myFourierTransformInverse(freq, hf, t0):
l = len(hf)
n = int(math.floor(l/2))
newhf = hf[n:]
newhf = np.append(newhf, hf[:n])
amp = np.fft.ifft(newhf, norm="ortho")
df = np.amin(np.diff(freq))
time = t0 + range(0, l)/(df*l)
return time, amp
def angular_momentum(x, q, m, chi1, chi2, LInitNR):
eta = q/(1.+q)**2.
m1 = (1.+(1.-4.*eta)**0.5)/2.
m2 = m - m1
S1 = m1**2. * chi1
S2 = m2**2. * chi2
Sl = S1+S2
Sigmal = S2/m2 - S1/m1
DeltaM = m1 - m2
mu = eta
nu = eta
GammaE = 0.5772156649;
e4 = -(123671./5760.)+(9037.* math.pi**2.)/1536.+(896.*GammaE)/15.+(-(498449./3456.)+(3157.*math.pi**2.)/576.)*nu+(301. * nu**2.)/1728.+(77.*nu**3.)/31104.+(1792. *math.log(2.))/15.
e5 = -55.13
j4 = -(5./7.)*e4+64./35.
j5 = -(2./3.)*e5-4988./945.-656./135. * eta;
a1 = -2.18522;
a2 = 1.05185;
a3 = -2.43395;
a4 = 0.400665;
a5 = -5.9991;
CapitalDelta = (1.-4.*eta)**0.5
l = (eta/x**(1./2.)*(
1. +
x*(3./2. + 1./6.*eta) +
x**2. *(27./8. - 19./8.*eta + 1./24.*eta**2.) +
x**3. *(135./16. + (-6889./144. + 41./24. * math.pi**2.)*eta + 31./24.*eta**2. + 7./1296.*eta**3.) +
x**4. *((2835./128.) + eta*j4 - (64.*eta*math.log(x)/3.))+
x**5. *((15309./256.) + eta*j5 + ((9976./105.) + (1312.*eta/15.))*eta*math.log(x))+
x**(3./2.)*(-(35./6.)*Sl - 5./2.*DeltaM* Sigmal) +
x**(5./2.)*((-(77./8.) + 427./72.*eta)*Sl + DeltaM* (-(21./8.) + 35./12.*eta)*Sigmal) +
x**(7./2.)*((-(405./16.) + 1101./16.*eta - 29./16.*eta**2.)*Sl + DeltaM*(-(81./16.) + 117./4.*eta - 15./16.*eta**2.)*Sigmal) +
(1./2. + (m1 - m2)/2. - eta)* chi1**2. * x**2. +
(1./2. + (m2 - m1)/2. - eta)* chi2**2. * x**2. +
2.*eta*chi1*chi2*x**2. +
((13.*chi1**2.)/9. +
(13.*CapitalDelta*chi1**2.)/9. -
(55.*nu*chi1**2.)/9. -
29./9.*CapitalDelta*nu*chi1**2. +
(14.*nu**2. *chi1**2.)/9. +
(7.*nu*chi1*chi2)/3. +
17./18.* nu**2. * chi1 * chi2 +
(13.* chi2**2.)/9. -
(13.*CapitalDelta*chi2**2.)/9. -
(55.*nu*chi2**2.)/9. +
29./9.*CapitalDelta*nu*chi2**2. +
(14.*nu**2. * chi2**2.)/9.)
* x**3.))
return l - LInitNR
#Get cutoff frequency
def getCutoffFrequency(sim_name):
"""
Determine cutoff frequency of simulation
sim_name = string of simulation
return = cutoff frequency
"""
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(" ")
if(line_elems[0] == "TwoPunctures::par_b"):
par_b = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::center_offset[0]"):
center_offset = float(line_elems[-1])
if(line_elems[0] == "TwoPunctures::par_P_plus[1]"):
pyp = float(line_elems[-1])
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])
if(line_elems[0] == "TwoPunctures::target_M_minus"):
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])
xp = par_b + center_offset
xm = -1*par_b + center_offset
LInitNR = xp*pyp + xm*pym
M = m1+m2
q = m1/m2
chi1 = S1/math.pow(m1, 2.)
chi2 = S2/math.pow(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.)
omGWPN = 2. * omOrbPN
omCutoff = 0.75 * omGWPN
return omCutoff
#-----Main-----#
#Initialize simulation data
if(len(sys.argv) < 2):
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).")
sys.exit()
elif(os.path.isdir(sys.argv[1])):
radiiUsedForExtrapolation = 7 #use the first n radii available
paths = sys.argv[1:]
elif(not os.path.isdir(sys.argv[1])):
radiiUsedForExtrapolation = int(sys.argv[1]) #use the first n radii available
if(radiiUsedForExtrapolation < 1 or radiiUsedForExtrapolation > 7):
print("Invalid specified radii number")
sys.exit()
paths = sys.argv[2:]
for sim_path in paths:
main_dir = sim_path
sim = sim_path.split("/")[-1]
simdirs = main_dir+"/output-????/%s/" % (sim)
f0 = getCutoffFrequency(sim)
#Check if necessary files exist
par_file = main_dir+"/output-0000/%s.par" % (sim)
two_punctures_file = main_dir+"/output-0000/%s/TwoPunctures.bbh" % (sim)
if(not os.path.isfile(par_file) or not os.path.isfile(two_punctures_file)):
continue
#Create data directories
main_directory = "Extrapolated_Strain"
sim_dir = main_directory+"/"+sim
if not os.path.exists(main_directory):
os.makedirs(main_directory)
if not os.path.exists(sim_dir):
os.makedirs(sim_dir)
#Get ADMMass
ADMMass = -1
filename = main_dir+"/output-0000/%s/TwoPunctures.bbh" % (sim)
with open(filename) as file:
contents = file.readlines()
for line in contents:
line_elems = line.split(" ")
if(line_elems[0] == "initial-ADM-energy"):
ADMMass = float(line_elems[-1])
#Get Psi4
radii = [100, 115, 136, 167, 214, 300, 500]
psi4dsetname_100 = 'l2_m2_r'+str(radii[0])+'.00'
psi4dsetname_115 = 'l2_m2_r'+str(radii[1])+'.00'
psi4dsetname_136 = 'l2_m2_r'+str(radii[2])+'.00'
psi4dsetname_167 = 'l2_m2_r'+str(radii[3])+'.00'
psi4dsetname_214 = 'l2_m2_r'+str(radii[4])+'.00'
psi4dsetname_300 = 'l2_m2_r'+str(radii[5])+'.00'
psi4dsetname_500 = 'l2_m2_r'+str(radii[6])+'.00'
mp_psi4_100 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname_100)
mp_psi4_115 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname_115)
mp_psi4_136 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname_136)
mp_psi4_167 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname_167)
mp_psi4_214 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname_214)
mp_psi4_300 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname_300)
mp_psi4_500 = loadHDF5Series(simdirs+"mp_psi4.h5", psi4dsetname_500)
mp_psi4_vars = [mp_psi4_100, mp_psi4_115, mp_psi4_136, mp_psi4_167, mp_psi4_214, mp_psi4_300, mp_psi4_500]
#Get Tortoise Coordinate
tortoise = [None] * 7
for i in range(0, 7):
tortoise[i] = -RadialToTortoise(radii[i], ADMMass)
strain = [None]*7
phase = [None]*7
amp = [None]*7
for i in range(0, 7):
#Get modified Psi4 (Multiply real and imaginary psi4 columns by radii and add the tortoise coordinate to the time colum)
mp_psi4_vars[i][:, 0] += tortoise[i]
mp_psi4_vars[i][:, 1] *= radii[i]
mp_psi4_vars[i][:, 2] *= radii[i]
#Fixed-frequency integration twice to get strain
hTable = psi4ToStrain(mp_psi4_vars[i], f0)
time = hTable[:, 0]
h = hTable[:, 1]
hplus = h.real
hcross = h.imag
newhTable = np.column_stack((time, hplus, hcross))
warnings.filterwarnings('ignore')
finalhTable = newhTable.astype(float)
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_strain_at_"+str(radii[i])+".dat", finalhTable)
strain[i] = finalhTable
#Get phase and amplitude of strain
h_phase = np.unwrap(np.angle(h))
while(h_phase[0] < 0):
h_phase[:] += 2*math.pi
angleTable = np.column_stack((time, h_phase))
angleTable = angleTable.astype(float)
phase[i] = angleTable
h_amp = np.absolute(h)
ampTable = np.column_stack((time, h_amp))
ampTable = ampTable.astype(float)
amp[i] = ampTable
#Interpolate phase and amplitude
t = phase[0][:, 0]
last_t = phase[radiiUsedForExtrapolation - 1][-1, 0]
last_index = 0;
for i in range(0, len(phase[0][:, 0])):
if(t[i] > last_t):
last_index = i
break
last_index = last_index-1
t = phase[0][0:last_index, 0]
dts = t[1:] - t[:-1]
dt = float(np.amin(dts))
t = np.arange(phase[0][0, 0], phase[0][last_index, 0], dt)
interpolation_order = 9
for i in range(0, radiiUsedForExtrapolation):
interp_function = scipy.interpolate.interp1d(phase[i][:, 0], phase[i][:, 1], kind=interpolation_order)
resampled_phase_vals = interp_function(t)
while(resampled_phase_vals[0] < 0):
resampled_phase_vals[:] += 2*math.pi
phase[i] = np.column_stack((t, resampled_phase_vals))
interp_function = scipy.interpolate.interp1d(amp[i][:, 0], amp[i][:, 1], kind=interpolation_order)
resampled_amp_vals = interp_function(t)
amp[i] = np.column_stack((t, resampled_amp_vals))
#Extrapolate
phase_extrapolation_order = 1
amp_extrapolation_order = 2
radii = np.asarray(radii, dtype=float)
radii = radii[0:radiiUsedForExtrapolation]
A_phase = np.power(radii, 0)
A_amp = np.power(radii, 0)
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)))
radially_extrapolated_phase = np.empty(0)
radially_extrapolated_amp = np.empty(0)
for i in range(0, len(t)):
b_phase = np.empty(0)
for j in range(0, radiiUsedForExtrapolation):
b_phase = np.append(b_phase, phase[j][i, 1])
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.append(b_amp, amp[j][i, 1])
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)
for i in range(0, len(radially_extrapolated_amp)):
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.dat", np.column_stack((t, radially_extrapolated_h_plus, radially_extrapolated_h_cross)))
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_amplitude.dat", np.column_stack((t, radially_extrapolated_amp)))
np.savetxt("./Extrapolated_Strain/"+sim+"/"+sim+"_radially_extrapolated_phase.dat", np.column_stack((t, radially_extrapolated_phase)))