diff --git a/util/skynet/README.md b/util/skynet/README.md deleted file mode 100644 index 47d208837e..0000000000 --- a/util/skynet/README.md +++ /dev/null @@ -1,35 +0,0 @@ -This directory contains: - -"compare.py" - loads the output of skynet and starkiller networks - for comparison tests -"burn.py" - runs skynet in parallel on problems stated therein - -NOTE: The "pressure" option in the skydata class in compare.py can be - used if the helmoholtz eos executable for finding pressure is - located in the directory. The helmholtz eos is freely available - from: https://cococubed.com/code_pages/eos.shtml - -Quick Start: - 1. Make sure Skynet is installed and working. Be sure to set your - SKYNET_HOME variable, and to prepend SKYNET_HOME to your - PYTHON_PATH. The above depend on importing built-in scripts from - Skynet, so we need to tell python where to find them. - - 2. Change desired problem parameters in burn.py as needed. - - >> python burn.py - - 3. Load results using the sky class in "compare.py" - - >> import compare as com - >> data = com.skydata(/path/to/h5/file/) - >> edot = data.edot - - - - - - - - - diff --git a/util/skynet/burn.py b/util/skynet/burn.py deleted file mode 100755 index e48f04f63a..0000000000 --- a/util/skynet/burn.py +++ /dev/null @@ -1,126 +0,0 @@ -#!/usr/bin/env python - -from SkyNet import * -import numpy as np -import shutil -import multiprocessing - -num_done = multiprocessing.Value('i', 0) -output_dir = 'isochoric_tests/' - -def run_skynet(args): - global num_done - # unpack arguments - rho, total_size = args - - # set up network and options - nuclib = NuclideLibrary.CreateFromWebnucleoXML(SkyNetRoot - + "/data/webnucleo_nuc_v2.0.xml") - - opts = NetworkOptions() - opts.ConvergenceCriterion = NetworkConvergenceCriterion.Mass - opts.MassDeviationThreshold = 1.0E-10 - opts.IsSelfHeating = True - opts.EnableScreening = True - opts.DisableStdoutOutput = False - - screen = SkyNetScreening(nuclib) - helm = HelmholtzEOS(SkyNetRoot + "/data/helm_table.dat") - - strongReactionLibrary = REACLIBReactionLibrary(SkyNetRoot + "/data/reaclib", - ReactionType.Strong, True, LeptonMode.TreatAllAsDecayExceptLabelEC, - "Strong reactions", nuclib, opts, True) - symmetricFission = REACLIBReactionLibrary(SkyNetRoot - + "/data/netsu_panov_symmetric_0neut", ReactionType.Strong, False, - LeptonMode.TreatAllAsDecayExceptLabelEC, - "Symmetric neutron induced fission with 0 neutrons emitted", nuclib, opts, - False) - spontaneousFission = REACLIBReactionLibrary(SkyNetRoot + - "/data/netsu_sfis_Roberts2010rates", ReactionType.Strong, False, - LeptonMode.TreatAllAsDecayExceptLabelEC, "Spontaneous fission", nuclib, opts, - False) - - # use only REACLIB weak rates - weakReactionLibrary = REACLIBReactionLibrary(SkyNetRoot + "/data/reaclib", - ReactionType.Weak, False, LeptonMode.TreatAllAsDecayExceptLabelEC, - "Weak reactions", nuclib, opts, True) - reactionLibraries = [strongReactionLibrary, symmetricFission, - spontaneousFission, weakReactionLibrary] - - net = ReactionNetwork(nuclib, reactionLibraries, helm, screen, opts) - nuclib = net.GetNuclideLibrary() - - # set up initial conditions - Y0 = np.zeros(nuclib.NumNuclides()) - - #URCA - Y0[nuclib.NuclideIdsVsNames()["c12"]] = 0.4998 / 12.0 - Y0[nuclib.NuclideIdsVsNames()["o16"]] = 0.4998 / 16.0 - #Y0[nuclib.NuclideIdsVsNames()["ne23"]] = 4e-4 / 23.0 # above threshold - Y0[nuclib.NuclideIdsVsNames()["na23"]] = 4e-4 / 23.0 # below threshold - - #XRB - #Y0[nuclib.NuclideIdsVsNames()["he4"]] = 1.0 / 4.0 - - #SubCh - #Y0[nuclib.NuclideIdsVsNames()["c12"]] = 0.08 / 12.0 - #Y0[nuclib.NuclideIdsVsNames()["o16"]] = 0.08 / 16.0 - #Y0[nuclib.NuclideIdsVsNames()["n14"]] = 0.01 / 14.0 - #Y0[nuclib.NuclideIdsVsNames()["he4"]] = 0.83 / 4.0 - - # T0 = 0.55 URCA, 0.2 XRB, 0.3 SubCh - T0 = 0.55 # initial temperature in GK - Ye = 0.5 # initial Ye - #s = 10.0 # initial entropy in k_B / baryon - #tau = 7.1 # expansion timescale in ms - - 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") - - # save some additional information, not necessary as all will be - # contained in h5 file - YvsA = np.array(output.FinalYVsA()) - A = np.arange(len(YvsA)) - - np.savetxt("urca_isoc_belowthresh.txt", np.array([A, YvsA]).transpose(), - "%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() - -# monitor progress with -# $ cd work_dir -# $ watch "bash -c 'tail -n 3 *.log'" - - diff --git a/util/skynet/compare.py b/util/skynet/compare.py deleted file mode 100755 index c728082c0b..0000000000 --- a/util/skynet/compare.py +++ /dev/null @@ -1,174 +0,0 @@ -#!/usr/bin/env python -from __future__ import print_function -from SkyNet import * -import SkyNet -import sys -import os -import subprocess -import argparse -import glob -import numpy as np -from cycler import cycler -import matplotlib.pyplot as plt - -# 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? - - -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): - ''' - return index of isotope in the self.z, self.z, self.isos, self.misos, - self.xisos[i,:] arrays - ''' - wh = np.where((self.A == a) & (self.Z == z))[0] - if len(wh) == 0: - 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')