diff --git a/util/skynet/burn.py b/util/skynet/burn.py index c6154ec55e..28d8d666b7 100755 --- a/util/skynet/burn.py +++ b/util/skynet/burn.py @@ -29,7 +29,7 @@ def run_skynet(args): strongReactionLibrary = REACLIBReactionLibrary(SkyNetRoot + "/data/reaclib", ReactionType.Strong, True, LeptonMode.TreatAllAsDecayExceptLabelEC, - "Strong reactions", nuclib, opts, True) + "Strong reactions", nuclib, opts, True) symmetricFission = REACLIBReactionLibrary(SkyNetRoot + "/data/netsu_panov_symmetric_0neut", ReactionType.Strong, False, LeptonMode.TreatAllAsDecayExceptLabelEC, @@ -42,8 +42,8 @@ def run_skynet(args): # use only REACLIB weak rates weakReactionLibrary = REACLIBReactionLibrary(SkyNetRoot + "/data/reaclib", - ReactionType.Weak, False, LeptonMode.TreatAllAsDecayExceptLabelEC, - "Weak reactions", nuclib, opts, True) + ReactionType.Weak, False, LeptonMode.TreatAllAsDecayExceptLabelEC, + "Weak reactions", nuclib, opts, True) reactionLibraries = [strongReactionLibrary, symmetricFission, spontaneousFission, weakReactionLibrary] @@ -59,7 +59,7 @@ def run_skynet(args): #Y0[nuclib.NuclideIdsVsNames()["ne23"]] = 4e-4 / 23.0 # above threshold Y0[nuclib.NuclideIdsVsNames()["na23"]] = 4e-4 / 23.0 # below threshold - #XRB + #XRB #Y0[nuclib.NuclideIdsVsNames()["he4"]] = 1.0 / 4.0 #SubCh @@ -74,12 +74,12 @@ def run_skynet(args): #s = 10.0 # initial entropy in k_B / baryon #tau = 7.1 # expansion timescale in ms - time = np.linspace(0.0, 1e9) + time = np.linspace(0.0, 1e9) density_vs_time = ConstantFunction(rho) #isochoric # run through network, will save as "prefix".h5 and .log files output = net.EvolveSelfHeatingWithInitialTemperature(Y0, 0.0, 1.0E1, - T0, density_vs_time, "urca_isoc_belowthresh") + T0, density_vs_time, "urca_isoc_belowthresh") # save some additional information, not necessary as all will be # contained in h5 file @@ -87,37 +87,37 @@ def run_skynet(args): A = np.arange(len(YvsA)) np.savetxt("urca_isoc_belowthresh.txt", np.array([A, YvsA]).transpose(), - "%6i %30.20E") + "%6i %30.20E") return if __name__ == '__main__': - num_cores = multiprocessing.cpu_count() - print "Running with %i worker threads" % num_cores - pool = multiprocessing.Pool(num_cores) - num_done.value = 0 - - # use if you have trajectory files to read - #files = os.listdir(input_dir) - - # otherwise, give the initial density - #densities = [1e6, 5e6, 1e7, 5e7, 1e8, 5e8, 1e9] # multiple - densities = [1.0e9] # URCA - #densities = [2e6] #XRB Aprox13 - #densities = [1e6] #SubCh - - # number of jobs to run - size = len(densities) - args = [(r, size) for r in densities] - - # run skynet in parallel - pool.map_async(run_skynet, args) - - # done submitting jobs - pool.close() - pool.join() + num_cores = multiprocessing.cpu_count() + print("Running with %i worker threads" % num_cores) + pool = multiprocessing.Pool(num_cores) + num_done.value = 0 + + # use if you have trajectory files to read + #files = os.listdir(input_dir) + + # otherwise, give the initial density + #densities = [1e6, 5e6, 1e7, 5e7, 1e8, 5e8, 1e9] # multiple + densities = [1.0e9] # URCA + #densities = [2e6] #XRB Aprox13 + #densities = [1e6] #SubCh + + # number of jobs to run + size = len(densities) + args = [(r, size) for r in densities] + + # run skynet in parallel + pool.map_async(run_skynet, args) + + # done submitting jobs + pool.close() + pool.join() # monitor progress with # $ cd work_dir diff --git a/util/skynet/compare.py b/util/skynet/compare.py index 48e6d7684f..c728082c0b 100755 --- a/util/skynet/compare.py +++ b/util/skynet/compare.py @@ -14,151 +14,151 @@ # Get Star Killer output class starkiller_data(): - # modified from Microphysics/unit_test/burn_cell/burn_cell.py - def __init__(self, inpath, runprefix='react_urca'): - files = glob.glob(inpath + runprefix + r'_[0-9]*') - data = [] - for fn in files: - d = {} - f = open(fn, 'r') - - # Get file number - fnsplit = fn.split('_') - fnum = int(fnsplit[-1]) - d['fnum'] = fnum - - # Eat lines - flines = [] - for l in f: - flines.append(l.strip()) - f.close() - - # Fill dataset - n = len(flines) - for i, xi in enumerate(flines): - if xi[-1] == ':': - # This is a field - key = xi[:-1] - xd = [] - # Collect following lines of field data - for j in range(i+1,n): - xj = flines[j] - if xj[-1] == ':': - break - else: - xd.append(xj.split()) - - # Interpret data - if (key=='nspec' or key=='neqs' or key=='i' or key=='j' or key=='k' - or key=='n_rhs' or key=='n_jac'): - # Integer - d[key] = int(xd[0][0]) - elif key=='short_spec_names': - # Strings of nuclide abbreviations - d[key] = [xj[0] for xj in xd] - elif key=='self_heat': - # Boolean - d[key] = (xd[0][0]=='T') - elif key=='xn' or key=='ydot': - # 1-D Float Array - d[key] = np.array([float(y) for y in xd[0]]) - elif key=='jac': - # 2-D Float Array - d[key] = np.array([[float(xij) for xij in xdi] for xdi in xd]) - else: - # Float - d[key] = float(xd[0][0]) - - # Store data - data.append(d) - - ## SORT file data by file number - data_sorted = sorted(data, key=lambda d: d['fnum']) - data = data_sorted - - ## INIT VARIABLES - nspec = data[0]['nspec'] - neqs = data[0]['neqs'] - self.short_spec_names = data[0]['short_spec_names'] - - # Init time, temp, ener - fnum = [] - time = [] - temp = [] - ener = [] - - # Init abundance vectors - xn = [[] for i in range(nspec)] - - # Init ydot vectors - ydot = [[] for i in range(neqs)] - - ## DATA LOOP - # Loop through data and collect - for d in data: - fnum.append(d['fnum']) - temp.append(d['T']) - ener.append(d['e']) - time.append(d['time']) - for i, xi in enumerate(d['xn']): - xn[i].append(xi) - for i, xi in enumerate(d['ydot']): - ydot[i].append(xi) - - # Convert data to numpy arrays - for i in range(nspec): - xn[i] = np.array(xn[i]) - self.X = xn - #for i in range(neqs): - # ydot[i] = np.array(ydot[i]) - # ydot[i][0] = ydot[i][1] - self.fnum = np.array(fnum) - self.temp = np.array(temp) - self.time = np.array(time) - self.dtime = np.ediff1d(time, to_begin=0) - self.ener = np.array(ener) - denerdt = np.zeros_like(ener) - denerdt[1:] = self.ener[1:]/self.dtime[1:] - self.edot = denerdt - # for now ydot is garbage might be useful later? + # modified from Microphysics/unit_test/burn_cell/burn_cell.py + def __init__(self, inpath, runprefix='react_urca'): + files = glob.glob(inpath + runprefix + r'_[0-9]*') + data = [] + for fn in files: + d = {} + f = open(fn, 'r') + + # Get file number + fnsplit = fn.split('_') + fnum = int(fnsplit[-1]) + d['fnum'] = fnum + + # Eat lines + flines = [] + for l in f: + flines.append(l.strip()) + f.close() + + # Fill dataset + n = len(flines) + for i, xi in enumerate(flines): + if xi[-1] == ':': + # This is a field + key = xi[:-1] + xd = [] + # Collect following lines of field data + for j in range(i+1,n): + xj = flines[j] + if xj[-1] == ':': + break + else: + xd.append(xj.split()) + + # Interpret data + if (key=='nspec' or key=='neqs' or key=='i' or key=='j' or key=='k' + or key=='n_rhs' or key=='n_jac'): + # Integer + d[key] = int(xd[0][0]) + elif key=='short_spec_names': + # Strings of nuclide abbreviations + d[key] = [xj[0] for xj in xd] + elif key=='self_heat': + # Boolean + d[key] = (xd[0][0]=='T') + elif key=='xn' or key=='ydot': + # 1-D Float Array + d[key] = np.array([float(y) for y in xd[0]]) + elif key=='jac': + # 2-D Float Array + d[key] = np.array([[float(xij) for xij in xdi] for xdi in xd]) + else: + # Float + d[key] = float(xd[0][0]) + + # Store data + data.append(d) + + ## SORT file data by file number + data_sorted = sorted(data, key=lambda d: d['fnum']) + data = data_sorted + + ## INIT VARIABLES + nspec = data[0]['nspec'] + neqs = data[0]['neqs'] + self.short_spec_names = data[0]['short_spec_names'] + + # Init time, temp, ener + fnum = [] + time = [] + temp = [] + ener = [] + + # Init abundance vectors + xn = [[] for i in range(nspec)] + + # Init ydot vectors + ydot = [[] for i in range(neqs)] + + ## DATA LOOP + # Loop through data and collect + for d in data: + fnum.append(d['fnum']) + temp.append(d['T']) + ener.append(d['e']) + time.append(d['time']) + for i, xi in enumerate(d['xn']): + xn[i].append(xi) + for i, xi in enumerate(d['ydot']): + ydot[i].append(xi) + + # Convert data to numpy arrays + for i in range(nspec): + xn[i] = np.array(xn[i]) + self.X = xn + #for i in range(neqs): + # ydot[i] = np.array(ydot[i]) + # ydot[i][0] = ydot[i][1] + self.fnum = np.array(fnum) + self.temp = np.array(temp) + self.time = np.array(time) + self.dtime = np.ediff1d(time, to_begin=0) + self.ener = np.array(ener) + denerdt = np.zeros_like(ener) + denerdt[1:] = self.ener[1:]/self.dtime[1:] + self.edot = denerdt + # for now ydot is garbage might be useful later? class skydata(): - def __init__(self, inpath, pressure = False): - - out = NetworkOutput.ReadFromFile(inpath) - self.time = np.array(out.Times()) - self.temp = np.array(out.TemperatureVsTime()) * 1.0E9 - self.rho = np.array(out.DensityVsTime()) - self.edot = np.array(out.HeatingRateVsTime()) - self.Y = np.array(out.YVsTime()) - self.A = np.array(out.A()) - self.Z = np.array(out.Z()) - self.X = np.array([self.A * y for y in self.Y[:]]) - - if pressure: - jmax = len(self.temp) - - # writing intermediate thermo file for EoS - # calls the Helmholtz EoS for pressure as a.out - f = open("tmp.txt", 'w') - f.write(str(jmax) + "\n") - - for i in range(jmax): - wbar = 1.0 / np.sum(self.Y[i,:]) - abar = wbar * np.sum(np.dot(self.A,self.Y[i,:])) - zbar = wbar * np.sum(np.dot(self.Z,self.Y[i,:])) - f.write("{:.15e} {:.15e} {:.15e} {:.15e}\n".format(self.temp[i], - self.rho[i], abar, zbar)) - f.close() - - process = subprocess.Popen("./a.out", shell=True) - process.wait() - self.pressure = np.loadtxt("tmp_out.txt") - os.system("rm tmp*") - - def isoidx(self, a, z): + def __init__(self, inpath, pressure = False): + + out = NetworkOutput.ReadFromFile(inpath) + self.time = np.array(out.Times()) + self.temp = np.array(out.TemperatureVsTime()) * 1.0E9 + self.rho = np.array(out.DensityVsTime()) + self.edot = np.array(out.HeatingRateVsTime()) + self.Y = np.array(out.YVsTime()) + self.A = np.array(out.A()) + self.Z = np.array(out.Z()) + self.X = np.array([self.A * y for y in self.Y[:]]) + + if pressure: + jmax = len(self.temp) + + # writing intermediate thermo file for EoS + # calls the Helmholtz EoS for pressure as a.out + f = open("tmp.txt", 'w') + f.write(str(jmax) + "\n") + + for i in range(jmax): + wbar = 1.0 / np.sum(self.Y[i,:]) + abar = wbar * np.sum(np.dot(self.A,self.Y[i,:])) + zbar = wbar * np.sum(np.dot(self.Z,self.Y[i,:])) + f.write("{:.15e} {:.15e} {:.15e} {:.15e}\n".format(self.temp[i], + self.rho[i], abar, zbar)) + f.close() + + process = subprocess.Popen("./a.out", shell=True) + process.wait() + self.pressure = np.loadtxt("tmp_out.txt") + os.system("rm tmp*") + + def isoidx(self, a, z): ''' return index of isotope in the self.z, self.z, self.isos, self.misos, self.xisos[i,:] arrays @@ -168,7 +168,7 @@ def isoidx(self, a, z): print('isotope not found') else: return wh[0] - + #sky = skydata(inpath = "./urca_abovethresh/urca_isoc_abovethresh.h5", pressure = False) #star = starkiller_data("./urca_abovethresh/", runprefix='react_urca')