From a31f6a24f81e2c612733980f13f64a2fd32ebcb5 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Mon, 2 Dec 2024 11:30:58 -0800 Subject: [PATCH 1/4] utilities to write sesame files from opacplot --- opacplot2/constants.py | 9 + opacplot2/opg_ionmix.py | 38 ++- opacplot2/opg_sesame.py | 381 ++++++++++++++++---------- opacplot2/scripts/opac_convert.py | 440 ++++++++++++++++++++---------- 4 files changed, 573 insertions(+), 295 deletions(-) diff --git a/opacplot2/constants.py b/opacplot2/constants.py index b9dd3d7..f4b1dbc 100644 --- a/opacplot2/constants.py +++ b/opacplot2/constants.py @@ -16,9 +16,18 @@ # Convert MJ/kg to Ergs/cc: MJKG_TO_ERGCC = 1.0e+10 +# Convert Ergs/cc to GPa: +ERGCC_TO_GPA =1./GPA_TO_ERGCC + +# Convert Ergs/g to MJ/kg: +ERGG_TO_MJKG =1./MJKG_TO_ERGG + # Convert Kelvin to eV: KELVIN_TO_EV = 1.0/11604.55 +# Convert eV to Kelvin: +EV_TO_KELVIN = 11604.55 + # Avogadros number NA = 6.0221415e+23 diff --git a/opacplot2/opg_ionmix.py b/opacplot2/opg_ionmix.py index bf2be12..55863b4 100644 --- a/opacplot2/opg_ionmix.py +++ b/opacplot2/opg_ionmix.py @@ -6,6 +6,7 @@ import numpy as np import re import math +import opacplot2 as opp from .opl_grid import OplGrid from .constants import ERG_TO_JOULE @@ -197,7 +198,7 @@ def u(x): self.temps = self.get_block(self.ntemp) self.numDens = self.get_block(self.ndens) - self.dens = self.numDens * self.mpi + self.dens = self.numDens * self.mpi/opp.NA if self.verb: print(" Number of temperatures: %i" % self.ntemp) @@ -507,6 +508,40 @@ def extendToZero(self): self.temps = arr self.ntemp += 1 + def toEosDict(self, Znum=None,Xnum=None, log=None): + eos_dict = {} + if Znum is None: + raise ValueError('Znum Varray should be provided!') + + if Xnum is None: + if len(Znum) == 1: + Xnum2 = np.array([1.0]) + else: + raise ValueError('Xnum array should be provided!') + else: + Xnum2 = np.array(Xnum) + eos_dict['idens' ] = self.numDens + eos_dict['temp' ] = self.temps + eos_dict['dens' ] = self.dens + eos_dict['Zf_DT' ] = self.zbar + eos_dict['Ut_DT' ] = self.eion+self.eele#self.etot + eos_dict['Uec_DT' ] = self.eele + eos_dict['Ui_DT' ] = self.eion + eos_dict['Pi_DT' ] = self.pion + eos_dict['Pec_DT' ] = self.pele + eos_dict['Znum' ] = np.array(Znum) + eos_dict['Xnum' ] = Xnum2 + eos_dict['BulkMod'] = 1. + print(self.mpi[0], opp.NA) + eos_dict['Abar' ] = self.mpi[0] + eos_dict['Zmax' ] = np.dot(np.array(Znum), Xnum2) + # mean opacities + eos_dict['opp_int'] = self.planck_absorb[:,:,0] + eos_dict['opr_int'] = self.rosseland[:,:,0] + + return eos_dict + + def writeIonmixFile(fn, zvals, fracs, numDens, temps, @@ -697,6 +732,7 @@ def write_opac_block(var, name): write_block(pele.flatten()*ERG_TO_JOULE, 'electron pressure') write_block(dpidt.flatten()*ERG_TO_JOULE, 'D(ion pressure)_DT') write_block(dpedt.flatten()*ERG_TO_JOULE, 'D(electron pressure)_DT') + print (eion[0,0], eion[0,0]*ERG_TO_JOULE, ERG_TO_JOULE) write_block(eion.flatten()*ERG_TO_JOULE, 'ion energy') write_block(eele.flatten()*ERG_TO_JOULE, 'electron energy') write_block(cvion.flatten()*ERG_TO_JOULE, 'ion CV') diff --git a/opacplot2/opg_sesame.py b/opacplot2/opg_sesame.py index 58ae982..e516a1b 100644 --- a/opacplot2/opg_sesame.py +++ b/opacplot2/opg_sesame.py @@ -1,19 +1,28 @@ from __future__ import absolute_import from __future__ import division from __future__ import print_function -#from __future__ import unicode_literals + +# from __future__ import unicode_literals from io import open import six +import textwrap import opacplot2 as opp import opacplot2.utils import numpy as np -from .constants import KELVIN_TO_EV, GPA_TO_ERGCC, MJKG_TO_ERGCC +from .constants import ( + KELVIN_TO_EV, + GPA_TO_ERGCC, + MJKG_TO_ERGCC, + ERGCC_TO_GPA, + ERGG_TO_MJKG, +) import periodictable as ptab + class OpgSesame: """ This class is responsible for loading all SESAME formatted data files. @@ -61,38 +70,39 @@ class OpgSesame: def __init__(self, filename, precision, verbose=False): self.verbose = verbose - self.fhand = open(filename, encoding='utf-8') - if(precision == self.SINGLE): - self.entry_len = 15 - elif(precision == self.DOUBLE): + self.fhand = open(filename, encoding="utf-8") + if precision == self.SINGLE: + self.entry_len = 15 + elif precision == self.DOUBLE: self.entry_len = 22 else: raise ValueError("precision must be SINGLE or DOUBLE") - - self.fdict = { 101 : self.parseComment, - 102 : self.parseComment, - 103 : self.parseComment, - 104 : self.parseComment, - 201 : self.parseInfo, - 301 : self.parseEos, - 303 : self.parseEos, - 304 : self.parseEos, - 305 : self.parseEos, - 306 : self.parseEos, - 311 : self.parseEos, - 321 : self.parseEos, - 401 : self.parseVap, - 411 : self.parseSolid, - 412 : self.parseLiquid, - 431 : self.parseShear, - 432 : self.parseShear, - 601 : self.parseZbar, - 602 : self.parseEcond, - 603 : self.parseTcond, - 604 : self.parseTcond, - 605 : self.parseTcond, - } + self.fdict = { + 101: self.parseComment, + 102: self.parseComment, + 103: self.parseComment, + 104: self.parseComment, + 201: self.parseInfo, + 301: self.parseEos, + 303: self.parseEos, + 304: self.parseEos, + 305: self.parseEos, + 306: self.parseEos, + 311: self.parseEos, + 321: self.parseEos, + 401: self.parseVap, + 411: self.parseSolid, + 412: self.parseLiquid, + 431: self.parseShear, + 432: self.parseShear, + 504: self.parseZbar, + 601: self.parseZbar, + 602: self.parseEcond, + 603: self.parseTcond, + 604: self.parseTcond, + 605: self.parseTcond, + } self.data = {} @@ -100,35 +110,41 @@ def __init__(self, filename, precision, verbose=False): self.parse() - def parse(self): while True: header = self.fhand.readline() - if not header: break # Reached EOF - if header[:3] == " 2 ": break + print(header.split()) + if not header: + break # Reached EOF + if header[:3] == " 2 ": + break matid = int(header.split()[1]) recid = int(header.split()[2]) nentries = int(header.split()[3]) - if not matid in self.data: self.data[matid] = {} + if not matid in self.data: + self.data[matid] = {} if self.verbose and (recid > 104): - print("Material = %8i Record = %8i Entries = %8i" % (matid, recid, nentries)) + print( + "Material = %8i Record = %8i Entries = %8i" + % (matid, recid, nentries) + ) # If there are opacity tables (series 500), # we do not currently read them. - if (600>recid) and (recid>=500): - print("Warning: No reader for SESAME opacity tables. The " - "data in table number {} will not be parsed".format(recid)) - + if (600 > recid) and (recid >= 500): + print( + "Warning: No reader for SESAME opacity tables. The " + "data in table number {} will not be parsed".format(recid) + ) if not recid in self.fdict: raise ValueError("No handling function for record %d" % recid) - - self.fdict[recid](nentries,matid, recid) + self.fdict[recid](nentries, matid, recid) if matid not in self.recs.keys(): self.recs[matid] = [recid] @@ -137,7 +153,7 @@ def parse(self): def parseComment(self, nentries, matid, recid): - nlines = (nentries-1) // self.CHAR_LINE_LEN + 1 + nlines = (nentries - 1) // self.CHAR_LINE_LEN + 1 nchar = nentries + nlines return self.fhand.read(nchar) @@ -149,18 +165,15 @@ def parseInfo(self, nentries, matid, recid): self.data[matid]["bulkmod"] = words[3] self.data[matid]["excoef"] = words[4] - if self.verbose: print(" zbar = %g abar = %g rho0 = %g" % (words[0],words[1],words[2])) + if self.verbose: + print(" zbar = %g abar = %g rho0 = %g" % (words[0], words[1], words[2])) return self.data def parseEos(self, nentries, matid, recid): words = self.readEntries(nentries) - prefixes = { 301 : "total_", - 303 : "ioncc_", - 304 : "ele_", - 305 : "ion_", - 306 : "cc_" } + prefixes = {301: "total_", 303: "ioncc_", 304: "ele_", 305: "ion_", 306: "cc_"} prefix = prefixes[recid] @@ -168,13 +181,13 @@ def parseEos(self, nentries, matid, recid): ntemp = int(words[1]) start = 2 - self.data[matid][prefix+"ndens"] = ndens - self.data[matid][prefix+"ntemp"] = ntemp + self.data[matid][prefix + "ndens"] = ndens + self.data[matid][prefix + "ntemp"] = ntemp - self.data[matid][prefix+"dens"] = words[start:start+ndens] # In g/cc - start = start+ndens + self.data[matid][prefix + "dens"] = words[start : start + ndens] # In g/cc + start = start + ndens - self.data[matid][prefix+"temps"] = words[start:start+ntemp]*KELVIN_TO_EV + self.data[matid][prefix + "temps"] = words[start : start + ntemp] * KELVIN_TO_EV start = start + ntemp if self.verbose and prefix == "total_": @@ -185,14 +198,19 @@ def parseEos(self, nentries, matid, recid): print("dens[0] = %13.6e dens[-1] = %13.6e" % (dens[0], dens[-1])) print("temp[0] = %13.6e temp[-1] = %13.6e" % (temps[0], temps[-1])) - # Read pressure array (in GPa and convert to ergs/cc): - self.data[matid][prefix+"pres"] = words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()*GPA_TO_ERGCC - start = start + ntemp*ndens + self.data[matid][prefix + "pres"] = ( + words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() + * GPA_TO_ERGCC + ) + start = start + ntemp * ndens # Read specific internal energy array (in MJ/kg and convert to ergs/g): - self.data[matid][prefix+"eint"] = words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()*MJKG_TO_ERGCC - start = start + ntemp*ndens + self.data[matid][prefix + "eint"] = ( + words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() + * MJKG_TO_ERGCC + ) + start = start + ntemp * ndens if start == nentries: # There is no Helmholtz free energy data, this is the end @@ -200,8 +218,11 @@ def parseEos(self, nentries, matid, recid): return # Read specific Helmholtz free energy array (in MJ/kg and convert to ergs/g): - self.data[matid][prefix+"free"] = words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()*MJKG_TO_ERGCC - start = start + ntemp*ndens + self.data[matid][prefix + "free"] = ( + words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() + * MJKG_TO_ERGCC + ) + start = start + ntemp * ndens # exit() @@ -226,16 +247,18 @@ def parseZbar(self, nentries, matid, recid): ntemp = int(words[1]) start = 2 - self.data[matid][prefix+"_ndens"] = ndens - self.data[matid][prefix+"_ntemp"] = ntemp + self.data[matid][prefix + "_ndens"] = ndens + self.data[matid][prefix + "_ntemp"] = ntemp - self.data[matid][prefix+"_dens"] = 10**words[start:start+ndens] - start = start+ndens + self.data[matid][prefix + "_dens"] = 10 ** words[start : start + ndens] + start = start + ndens - self.data[matid][prefix+"_temps"] = 10**words[start:start+ntemp] + self.data[matid][prefix + "_temps"] = 10 ** words[start : start + ntemp] start = start + ntemp - self.data[matid][prefix] = 10**(words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()) + self.data[matid][prefix] = 10 ** ( + words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() + ) def parseEcond(self, nentries, matid, recid): words = self.readEntries(nentries) @@ -243,32 +266,42 @@ def parseEcond(self, nentries, matid, recid): def parseTcond(self, nentries, matid, recid): words = self.readEntries(nentries) - def readEntries(self,nentries): - nlines = (nentries-1) // self.WORDS_PER_LINE + 1 + def readEntries(self, nentries): + nlines = (nentries - 1) // self.WORDS_PER_LINE + 1 data = np.empty(nentries) string = "" for i in range(nlines): - string += self.fhand.readline()[:self.WORDS_PER_LINE*self.entry_len] + string += self.fhand.readline()[: self.WORDS_PER_LINE * self.entry_len] for i in range(nentries): - word = string[self.entry_len*i:self.entry_len*(i+1)] - if word[-4] == '-': word = word[:-4] + "E" + word[-4:] + word = string[self.entry_len * i : self.entry_len * (i + 1)] + if word[-4] == "-": + word = word[:-4] + "E" + word[-4:] data[i] = float(word) return data - def toEosDict(self, Znum=None, Anum=None, - Xnum=None, qeos=False, log=None, - filter_dens=0., filter_temps=0., - tabnum=None): + def toEosDict( + self, + Znum=None, + Anum=None, + Xnum=None, + qeos=False, + log=None, + filter_dens=0.0, + filter_temps=0.0, + tabnum=None, + ): # For SESAME, we need the hedp package to calculate zbar. try: from hedp import eos except ImportError: - raise ImportError('You need the hedp module. You can get it here: ' - 'https://github.com/luli/hedp.') + raise ImportError( + "You need the hedp module. You can get it here: " + "https://github.com/luli/hedp." + ) if tabnum is None: # Select the last table (newest) table available. @@ -277,7 +310,7 @@ def toEosDict(self, Znum=None, Anum=None, try: opp_ses_data = self.data[tabnum] except KeyError: - raise KeyError('Invalid table number!') + raise KeyError("Invalid table number!") # Sesame has extra data points in it, so we must merge them down. We are # merging the default grids ioncc_ and ele_. @@ -286,100 +319,110 @@ def toEosDict(self, Znum=None, Anum=None, # We must filter dens > 0. in order to avoid problems with calculating # zbar below, since eos.thomas_fermi_ionization() returns `nan` where # density is 0. - if not qeos: - opp_ses_data = opp.utils.EosMergeGrids( - opp_ses_data, - filter_dens=lambda x: (x>filter_dens), - filter_temps=lambda x: (x>filter_temps)) + # if not qeos: + # opp_ses_data = opp.utils.EosMergeGrids( + # opp_ses_data, + # filter_dens=lambda x: (x>filter_dens), + # filter_temps=lambda x: (x>filter_temps)) # If we are dealing with SESAME generated by qeos, we need to merge # ion_ and ele_ grids rather than ioncc_ and ele_. We use the same # filters as above. if qeos: opp_ses_data = opp.utils.EosMergeGrids( - opp_ses_data, intersect=['ele', 'ion'], - filter_dens=lambda x: (x>filter_dens), - filter_temps=lambda x: (x>filter_temps), - qeos=True) + opp_ses_data, + intersect=["ele", "ion"], + filter_dens=lambda x: (x > filter_dens), + filter_temps=lambda x: (x > filter_temps), + qeos=True, + ) # Converting density to ion number density. - opp_ses_data['idens'] = ((opp.NA * opp_ses_data['ele_dens']) - / opp_ses_data['abar']) - + opp_ses_data["idens"] = (opp.NA * opp_ses_data["ele_dens"]) / opp_ses_data[ + "abar" + ] # Adjust for Znum, Xnum, and Anum. if Znum is None: - if 'Znum' in opp_ses_data: - Znum = opp_ses_data['Znum'] + if "Znum" in opp_ses_data: + Znum = opp_ses_data["Znum"] else: - raise ValueError('Znum Varray should be provided!') + raise ValueError("Znum Varray should be provided!") if type(Znum) is int: Znum = [Znum] - opp_ses_data['Znum'] = np.array(Znum, dtype='int') + opp_ses_data["Znum"] = np.array(Znum, dtype="int") if Anum is None: - opp_ses_data['Anum'] = np.array([ptab.elements[el].mass for el in opp_ses_data['Znum']]) + opp_ses_data["Anum"] = np.array( + [ptab.elements[el].mass for el in opp_ses_data["Znum"]] + ) else: - opp_ses_data['Anum'] = np.array(Anum) - opp_ses_data['Zsymb'] = np.array([ptab.elements[el].symbol for el in opp_ses_data['Znum']], dtype='|S2') + opp_ses_data["Anum"] = np.array(Anum) + opp_ses_data["Zsymb"] = np.array( + [ptab.elements[el].symbol for el in opp_ses_data["Znum"]], dtype="|S2" + ) if Xnum is None: if len(Znum) == 1: - opp_ses_data['Xnum'] = np.array([1.0]) + opp_ses_data["Xnum"] = np.array([1.0]) else: - raise ValueError('Xnum array should be provided!') + raise ValueError("Xnum array should be provided!") else: - opp_ses_data['Xnum'] = np.array(Xnum) - - - # Calculate zbar using thomas_fermi_ionization. - # If there are multiple elements, it suffices to use the average - # atomic number in this calculation - JTL - dens_arr, temp_arr = np.meshgrid(opp_ses_data['ele_dens'], - opp_ses_data['ele_temps']) - zbar = eos.thomas_fermi_ionization(dens_arr, - temp_arr, - np.average(opp_ses_data['Znum'], weights=opp_ses_data['Xnum']), - opp_ses_data['abar']).T - opp_ses_data['zbar'] = zbar + opp_ses_data["Xnum"] = np.array(Xnum) + + # turn this off later. just for testing sesame writer + # ionization is possibly on a completely different grid + if not "zbar" in opp_ses_data.keys(): + # Calculate zbar using thomas_fermi_ionization. + # If there are multiple elements, it suffices to use the average + # atomic number in this calculation - JTL + dens_arr, temp_arr = np.meshgrid( + opp_ses_data["ele_dens"], opp_ses_data["ele_temps"] + ) + zbar = eos.thomas_fermi_ionization( + dens_arr, temp_arr, opp_ses_data["Znum"].mean(), opp_ses_data["abar"] + ).T + opp_ses_data["zbar"] = zbar # Translating SESAME names to common dictionary format. if qeos: # Names are slightly different for QEOS SESAME - names_dict = {'idens':'idens', - 'ele_temps':'temp', # We merged ele_ and ion_ dens & - # temp grids for qeos. - 'ele_dens':'dens', - 'zbar':'Zf_DT', - 'total_eint':'Ut_DT', # But not their energies. - 'ele_eint':'Uec_DT', - 'ion_eint':'Ui_DT', - 'ion_pres':'Pi_DT', - 'ele_pres':'Pec_DT', - 'Znum':'Znum', - 'Xnum':'Xnum', - 'bulkmod':'BulkMod', - 'abar':'Abar', - 'zmax':'Zmax' - } + names_dict = { + "idens": "idens", + "ele_temps": "temp", # We merged ele_ and ion_ dens & + # temp grids for qeos. + "ele_dens": "dens", + "zbar": "Zf_DT", + "total_eint": "Ut_DT", # But not their energies. + "ele_eint": "Uec_DT", + "ion_eint": "Ui_DT", + "ion_pres": "Pi_DT", + "ele_pres": "Pec_DT", + "Znum": "Znum", + "Xnum": "Xnum", + "bulkmod": "BulkMod", + "abar": "Abar", + "zmax": "Zmax", + } else: - names_dict = {'idens':'idens', - 'ele_temps':'temp', # We merged ele_ and ioncc_ dens & - # temp grids. - 'ele_dens':'dens', - 'zbar':'Zf_DT', - 'total_eint':'Ut_DT', # But not their energies. - 'ele_eint':'Uec_DT', - 'ioncc_eint':'Ui_DT', - 'ioncc_pres':'Pi_DT', - 'ele_pres':'Pec_DT', - 'Znum':'Znum', - 'Xnum':'Xnum', - 'bulkmod':'BulkMod', - 'abar':'Abar', - 'zmax':'Zmax' - } + names_dict = { + "idens": "idens", + "ele_temps": "temp", # We merged ele_ and ioncc_ dens & + # temp grids. + "ele_dens": "dens", + "zbar": "Zf_DT", + "total_eint": "Ut_DT", # But not their energies. + "ele_eint": "Uec_DT", + "ioncc_eint": "Ui_DT", + "ioncc_pres": "Pi_DT", + "ele_pres": "Pec_DT", + "Znum": "Znum", + "Xnum": "Xnum", + "bulkmod": "BulkMod", + "abar": "Abar", + "zmax": "Zmax", + } # Initialize dictionary. eos_dict = {} @@ -389,7 +432,7 @@ def toEosDict(self, Znum=None, Anum=None, try: eos_dict[eos_key] = opp_ses_data[ses_key] except KeyError: - print('No data for {} is being written.'.format(eos_key)) + print("No data for {} is being written.".format(eos_key)) # Handle the logarithmic data. if log is not None: @@ -397,3 +440,51 @@ def toEosDict(self, Znum=None, Anum=None, if key in log: eos_dict[key] = np.log10(eos_dict[key]) return eos_dict + + +def writeSesameFile(fn, t201, t301, t303, t304, t305, t502, t504, t505, t601): + CHAR_LINE_LEN = 80 + WORDS_PER_LINE = 5 + comment = "This Sesame table was generated using OpacPlot2 to convert an IONMIX EoS table to the Sesame format" + f = open(fn, "w") + # write comments + comment = comment + " " * (80 - len(comment) % 80) + header = " 0 9999 101 %d r 82803 22704 1" % len(comment) + header = header + (79 - len(header)) * " " + "0\n" + f.write(header) + for i in range(len(comment) // 80): + f.write(comment[i * 80 : (i + 1) * 80] + "\n") + + def theader(num): + return ( + " 1 9999 %d 240 r 82803 22704 1 1\n" + % num + ) + + def write_block(fid, num, data): + header = " 1 9999 %d %d r 82803 22704 1" % (num, len(data)) + header = header + (79 - len(header)) * " " + "1\n" + fid.write(header) + count = 0 + for n in range(len(data)): + count += 1 + # fid.write('%22.15E' % data[n]) + fid.write("%15.8E" % data[n]) + if count == 5: + count = 0 + fid.write("11111\n") + if count > 0: + m = 5 - count + fid.write(m * 22 * " " + count * "1" + m * "0" + "\n") + # fid.write(m*16*' ' + count*'1' + m*'0' + '\n') + + write_block(f, 201, t201) + write_block(f, 301, t301) + write_block(f, 303, t303) + write_block(f, 304, t304) + write_block(f, 305, t305) + write_block(f, 502, t502) + write_block(f, 504, t504) + write_block(f, 505, t505) + write_block(f, 601, t601) + f.close() diff --git a/opacplot2/scripts/opac_convert.py b/opacplot2/scripts/opac_convert.py index d65dd35..e7a414e 100644 --- a/opacplot2/scripts/opac_convert.py +++ b/opacplot2/scripts/opac_convert.py @@ -1,57 +1,76 @@ import opacplot2 as opp import argparse import os.path +import numpy as np +from opacplot2.constants import EV_TO_KELVIN, ERGCC_TO_GPA, ERGG_TO_MJKG + def get_input_data(): # Available formats. - avail_output_formats = ['ionmix'] - avail_input_formats = ['propaceos', 'multi', 'sesame', 'sesame-qeos', 'tops'] + avail_output_formats = ["ionmix", "sesame"] + avail_input_formats = ["propaceos", "multi", "sesame", "sesame-qeos", "ionmix"] # Creating the argument parser. parser = argparse.ArgumentParser( - description="This script is used to browse various" - "EoS/Opacity tables formats.") - - parser.add_argument('-v', '--verbose', - action='store_const', const=True, - help='Verbosity option.') - - parser.add_argument('--Znum', - action='store', type=str, - help='Comma separated list of Z' - 'numbers for every component.') - - parser.add_argument('--Xfracs', - action='store', type=str, - help='Comma separated list of X' - 'fractions for every component.') - - parser.add_argument('--outname', - action='store', type=str, - help='Name for output file without extension.') - - parser.add_argument('-i', '--input', - action='store', type=str, - choices=avail_input_formats, - help='Input filetype.') - - parser.add_argument('-o', '--output', - action='store', type=str, - choices=avail_output_formats, - default='ionmix', - help='Output filetype. Default: IONMIX.') - - parser.add_argument('input_file', - action='store', type=str, - help='Input file.') - - parser.add_argument('--log', - action='store', type=str, - help='Logarithmic data keys.') - - parser.add_argument('--tabnum', - action='store', type=str, - help='Specify the SESAME table number.') + description="This script is used to browse various" + "EoS/Opacity tables formats." + ) + + parser.add_argument( + "-v", "--verbose", action="store_const", const=True, help="Verbosity option." + ) + + parser.add_argument( + "--Znum", + action="store", + type=str, + help="Comma separated list of Z" "numbers for every component.", + ) + + parser.add_argument("--mpi", action="store", type=str, help="Mass per ion in grams") + + parser.add_argument( + "--Xfracs", + action="store", + type=str, + help="Comma separated list of X" "fractions for every component.", + ) + + parser.add_argument( + "--outname", + action="store", + type=str, + help="Name for output file without extension.", + ) + + parser.add_argument( + "-i", + "--input", + action="store", + type=str, + choices=avail_input_formats, + help="Input filetype.", + ) + + parser.add_argument( + "-o", + "--output", + action="store", + type=str, + choices=avail_output_formats, + default="ionmix", + help="Output filetype. Default: IONMIX.", + ) + + parser.add_argument("input_file", action="store", type=str, help="Input file.") + + parser.add_argument( + "--log", action="store", type=str, help="Logarithmic data keys." + ) + + parser.add_argument( + "--tabnum", action="store", type=str, help="Specify the SESAME table number." + ) args = parser.parse_args() @@ -69,43 +88,50 @@ def get_input_data(): else: args.outname = os.path.join(basedir, basename) - # Create lists out of the strings for Znum, Xfracs, and log if given. + # Create lists out of the strings for Znum, Xfracs, mpi, and log if given. if args.Znum is not None: - args.Znum = [int(num) for num in args.Znum.split(',')] + args.Znum = [int(num) for num in args.Znum.split(",")] if args.Xfracs is not None: - args.Xfracs = [float(num) for num in args.Xfracs.split(',')] + args.Xfracs = [float(num) for num in args.Xfracs.split(",")] + if args.mpi is not None: + args.mpi = [float(num) for num in args.mpi.split(",")] if args.log is not None: - args.log = [str(key) for key in args.log.split(',')] + args.log = [str(key) for key in args.log.split(",")] # Convert tabnum into int. if args.tabnum is not None: try: args.tabnum = int(args.tabnum) except ValueError: - raise ValueError('Please provide a valid SESAME table number.') + raise ValueError("Please provide a valid SESAME table number.") - input_data = {'args' : args, - 'basename' : basename, - 'path_in' : path_in, - 'basedir' : basedir, - 'fn_in' : fn_in} + input_data = { + "args": args, + "basename": basename, + "path_in": path_in, + "basedir": basedir, + "fn_in": fn_in, + } return input_data + def read_format_ext(args, fn_in): # Try to read from the input file extension. - ext_dict = {'.prp':'propaceos', - '.eps':'multi', - '.opp':'multi', - '.opz':'multi', - '.opr':'multi', - '.mexport':'sesame-qeos', - '.ses':'sesame', - '.html':'tops', - '.tops':'tops', - } + ext_dict = { + ".prp": "propaceos", + ".eps": "multi", + ".opp": "multi", + ".opz": "multi", + ".opr": "multi", + ".cn4": "ionmix", + ".mexport": "sesame-qeos", + ".ses": "sesame", + ".html": "tops", + ".tops": "tops", + } # If the input file is compressed, choose the next extension. - if os.path.splitext(fn_in)[1] == '.gz': + if os.path.splitext(fn_in)[1] == ".gz": _, ext = os.path.splitext(os.path.splitext(fn_in)[0]) else: _, ext = os.path.splitext(fn_in) @@ -115,14 +141,18 @@ def read_format_ext(args, fn_in): if ext in ext_dict.keys(): args.input = ext_dict[ext] else: - raise Warning('Cannot tell filetype from extension. Please specify ' - 'input file type with --input.') + raise Warning( + "Cannot tell filetype from extension. Please specify " + "input file type with --input." + ) + class Formats_toEosDict(object): """ Contains handling functions to convert different types of tables into a common EoS dictionary for IONMIX. """ + def __init__(self, args, basedir, basename, path_in): # Initialize the dictionary for handling functions. self.set_handle_dict() @@ -137,32 +167,44 @@ def __init__(self, args, basedir, basename, path_in): try: self.eos_dict = self.handle_dict[args.input]() except KeyError: - raise KeyError('Must use valid format name.') + raise KeyError("Must use valid format name.") def set_handle_dict(self): - self.handle_dict = {'propaceos' : self.propaceos_toEosDict, - 'multi' : self.multi_toEosDict, - 'sesame' : self.sesame_toEosDict, - 'sesame-qeos' : self.sesame_qeos_toEosDict, - 'tops' : self.tops_toEosDict, - } + self.handle_dict = { + "propaceos": self.propaceos_toEosDict, + "multi": self.multi_toEosDict, + "sesame": self.sesame_toEosDict, + "ionmix": self.ionmix_toEosDict, + "tops": self.tops_toEosDict, + "sesame-qeos": self.sesame_qeos_toEosDict, + } def propaceos_toEosDict(self): # If we are unable to find the correct library for opg_propaceos # we need to let the user know. try: import opacplot2.opg_propaceos - op = opp.opg_propaceos.OpgPropaceosAscii(self.path_in) + + op = opp.opg_propaceos.OpgPropaceosAscii(self.path_in, mpi=self.args.mpi) eos_dict = op.toEosDict(log=self.args.log) return eos_dict except ImportError: - raise ImportError('You do not have the opg_propaceos script.') + raise ImportError("You do not have the opg_propaceos script.") def multi_toEosDict(self): op = opp.OpgMulti.open_file(self.basedir, self.basename) - eos_dict = op.toEosDict(Znum=self.args.Znum, - Xnum=self.args.Xfracs, - log=self.args.log) + eos_dict = op.toEosDict( + Znum=self.args.Znum, Xnum=self.args.Xfracs, log=self.args.log + ) + return eos_dict + + def ionmix_toEosDict(self): + print(self.args.mpi) + op = opp.OpacIonmix(self.path_in, self.args.mpi, man=True, twot=True) + eos_dict = op.toEosDict( + Znum=self.args.Znum, Xnum=self.args.Xfracs, log=self.args.log + ) + return eos_dict def sesame_toEosDict(self): @@ -171,43 +213,52 @@ def sesame_toEosDict(self): except ValueError: op = opp.OpgSesame(self.path_in, opp.OpgSesame.DOUBLE) - if len(op.data.keys()) > 1: - raise Warning('More than one material ID found. ' - 'Use sesame-extract to create a file ' - 'with only one material first.') + raise Warning( + "More than one material ID found. " + "Use sesame-extract to create a file " + "with only one material first." + ) if self.args.tabnum is not None: - eos_dict = op.toEosDict(Znum=self.args.Znum, - Xnum=self.args.Xfracs, - log=self.args.log, - tabnum=self.args.tabnum) + eos_dict = op.toEosDict( + Znum=self.args.Znum, + Xnum=self.args.Xfracs, + log=self.args.log, + tabnum=self.args.tabnum, + ) else: - eos_dict = op.toEosDict(Znum=self.args.Znum, - Xnum=self.args.Xfracs, - log=self.args.log) + eos_dict = op.toEosDict( + Znum=self.args.Znum, Xnum=self.args.Xfracs, log=self.args.log + ) return eos_dict def sesame_qeos_toEosDict(self): - raise Warning('QEOS-SESAME is not ready yet!') + raise Warning("QEOS-SESAME is not ready yet!") try: op = opp.OpgSesame(self.path_in, opp.OpgSesame.SINGLE) except ValueError: op = opp.OpgSesame(self.path_in, opp.OpgSesame.DOUBLE) - if len(op.data.keys()) > 1: - raise Warning('More than one material ID found. ' - 'Use sesame-extract to create a file ' - 'with only one material first.') + raise Warning( + "More than one material ID found. " + "Use sesame-extract to create a file " + "with only one material first." + ) if self.args.tabnum is not None: - eos_dict = op.toEosDict(Znum=self.args.Znum, Xnum=self.args.Xfracs, - qeos=True, log=self.args.log, - tabnum=self.args.tabnum) + eos_dict = op.toEosDict( + Znum=self.args.Znum, + Xnum=self.args.Xfracs, + qeos=True, + log=self.args.log, + tabnum=self.args.tabnum, + ) else: - eos_dict = op.toEosDict(Znum=self.args.Znum, Xnum=self.args.Xfracs, - qeos=True, log=self.args.log) + eos_dict = op.toEosDict( + Znum=self.args.Znum, Xnum=self.args.Xfracs, qeos=True, log=self.args.log + ) return eos_dict def tops_toEosDict(self): @@ -215,10 +266,88 @@ def tops_toEosDict(self): eos_dict = op.toEosDict(fill_eos=True) return eos_dict + +class EosDict_toSesameFile(object): + """ + Takes a common EoS dictionary and writes it to the correct output format. + """ + + def __init__(self, args, eos_dict): + self.set_handle_dict() + self.args = args + self.eos_dict = eos_dict + + self.handle_dict[args.output]() + + def set_handle_dict(self): + self.handle_dict = {"sesame": self.eosDict_toSesame} + + def eosDict_toSesame(self): + # initialize sesame argument dictionary + ses_dict = {} + # we should need to convert units to what sesame needs + dens = np.array(self.eos_dict["dens"]) + temp = np.array(self.eos_dict["temp"]) + pele = np.array(self.eos_dict["Pec_DT"]) + pion = np.array(self.eos_dict["Pi_DT"]) + uele = np.array(self.eos_dict["Uec_DT"]) + uion = np.array(self.eos_dict["Ui_DT"]) + utot = np.array(self.eos_dict["Ut_DT"]) + dummy = np.array(self.eos_dict["Ut_DT"]) + zbar = np.array(self.eos_dict["Zf_DT"]) + plnk = np.array(self.eos_dict["opp_int"]) + rsln = np.array(self.eos_dict["opr_int"]) + ptot = pele + pion + + if len(self.eos_dict["Znum"]) > 1: + zz = self.eos_dict["Znum"] + xx = self.eos_dict["Xnum"] + znum = 0.0 + for i in range(len(self.eos_dict["Znum"])): + print(zz[i]) + znum += zz[i] * xx[i] + print(znum) + else: + znum = self.eos_dict["Znum"][0] + + ses_dict["t201"] = np.array([znum, self.eos_dict["Abar"], 1.0, 1.0, 1.0]) + ses_dict["t301"] = self.tables_toSesame(dens, temp, ptot, utot, dummy) + ses_dict["t303"] = self.tables_toSesame(dens, temp, pion, uion, dummy) + ses_dict["t304"] = self.tables_toSesame(dens, temp, pele, uele, dummy) + ses_dict["t305"] = self.tables_toSesame(dens, temp, pion, uion, dummy) + ses_dict["t502"] = self.zbar_toSesame(dens, temp, rsln) + ses_dict["t505"] = self.zbar_toSesame(dens, temp, plnk) + ses_dict["t504"] = self.zbar_toSesame(dens, temp, zbar) + ses_dict["t601"] = self.zbar_toSesame(dens, temp, zbar) + + opp.writeSesameFile(self.args.outname + ".ses", **ses_dict) + + def tables_toSesame(self, dens, temp, pres, enrg, fnrg): + # flatten (n,t) tables into sesame array for 301-305 tables + ses_tab = np.array([len(dens), len(temp)]) + ses_tab = np.append(ses_tab, dens) + ses_tab = np.append(ses_tab, temp * EV_TO_KELVIN) + ses_tab = np.append(ses_tab, np.transpose(pres).flatten() * ERGCC_TO_GPA) + ses_tab = np.append(ses_tab, np.transpose(enrg).flatten() * ERGG_TO_MJKG) + ses_tab = np.append(ses_tab, np.transpose(fnrg).flatten()) + return ses_tab + + def zbar_toSesame(self, dens, temp, data): + Ldens = np.log10(dens) + Ltemp = np.log10(temp) + Ldata = np.log10(np.transpose(data).flatten()) + ses_tab = np.array([len(Ldens), len(Ltemp)]) + ses_tab = np.append(ses_tab, Ldens) + ses_tab = np.append(ses_tab, Ltemp) + ses_tab = np.append(ses_tab, Ldata) + return ses_tab + + class EosDict_toIonmixFile(object): """ Takes a common EoS dictionary and writes it to the correct output format. """ + def __init__(self, args, eos_dict): # Initialize the handling function dictionary. self.set_handle_dict() @@ -231,23 +360,25 @@ def __init__(self, args, eos_dict): self.handle_dict[args.output]() def set_handle_dict(self): - self.handle_dict = {'ionmix' : self.eosDict_toIonmix} + self.handle_dict = {"ionmix": self.eosDict_toIonmix} def eosDict_toIonmix(self): # These are the naming conventions translated to ionmix arguments. - imx_conv = {'Znum':'zvals', - 'Xnum':'fracs', - 'idens':'numDens', - 'temp':'temps', - 'Zf_DT':'zbar', - 'Pi_DT':'pion', - 'Pec_DT':'pele', - 'Ui_DT':'eion', - 'Uec_DT':'eele', - 'groups':'opac_bounds', - 'opr_mg':'rosseland', - 'opp_mg':'planck_absorb', - 'emp_mg':'planck_emiss'} + imx_conv = { + "Znum": "zvals", + "Xnum": "fracs", + "idens": "numDens", + "temp": "temps", + "Zf_DT": "zbar", + "Pi_DT": "pion", + "Pec_DT": "pele", + "Ui_DT": "eion", + "Uec_DT": "eele", + "groups": "opac_bounds", + "opr_mg": "rosseland", + "opp_mg": "planck_absorb", + "emp_mg": "planck_emiss", + } # Initialize ionmix argument dictionary. imx_dict = {} @@ -258,63 +389,74 @@ def eosDict_toIonmix(self): imx_dict[imx_conv[key]] = self.eos_dict[key] # Set ngroups if opacity bounds are present. - if 'opac_bounds' in imx_dict: - imx_dict['ngroups'] = len(self.eos_dict['groups']) - 1 + if "opac_bounds" in imx_dict: + imx_dict["ngroups"] = len(self.eos_dict["groups"]) - 1 # Check if required FLASH EoS tables are present. - imx_req_keys = ['zbar', 'eion', 'eele', 'pion', 'pele'] + imx_req_keys = ["zbar", "eion", "eele", "pion", "pele"] if not all(key in imx_dict for key in imx_req_keys): - print("The required EoS data for FLASH is not present...\n" - "Aborting the IONMIX file creation...") + print( + "The required EoS data for FLASH is not present...\n" + "Aborting the IONMIX file creation..." + ) raise Warning("Missing EoS data for IONMIX file to run in FLASH") # For verbose flag. if self.args.verbose: - verb_conv = {'zvals':'Atomic numbers', - 'fracs':'Element fractions', - 'numDens':'Ion number densities', - 'temps':'Temperatures', - 'zbar':'Average ionizations', - 'pion':'Ion pressure', - 'pele':'Electron pressure', - 'eion':'Ion internal energy', - 'eele':'Electron internal energy', - 'opac_bounds':'Opacity bounds', - 'rosseland':'Rosseland mean opacity', - 'planck_absorb':'Absorption Planck mean opacity', - 'planck_emiss':'Emission Planck mean opacity', - 'ngroups':'Number of opacity groups'} - - verb_str = 'Wrote the following data to IONMIX file:\n' + verb_conv = { + "zvals": "Atomic numbers", + "fracs": "Element fractions", + "numDens": "Ion number densities", + "temps": "Temperatures", + "zbar": "Average ionizations", + "pion": "Ion pressure", + "pele": "Electron pressure", + "eion": "Ion internal energy", + "eele": "Electron internal energy", + "opac_bounds": "Opacity bounds", + "rosseland": "Rosseland mean opacity", + "planck_absorb": "Absorption Planck mean opacity", + "planck_emiss": "Emission Planck mean opacity", + "ngroups": "Number of opacity groups", + } + + verb_str = "Wrote the following data to IONMIX file:\n" i = 0 for key in imx_dict.keys(): i = i + 1 if i == len(imx_dict.keys()): - verb_str = verb_str + '{}. {}'.format(i, verb_conv[key]) + verb_str = verb_str + "{}. {}".format(i, verb_conv[key]) else: - verb_str = verb_str + '{}. {} \n'.format(i, verb_conv[key]) + verb_str = verb_str + "{}. {} \n".format(i, verb_conv[key]) print(verb_str) # Write the ionmix file based on what data is stored in imx_dict. - opp.writeIonmixFile(self.args.outname + '.cn4', **imx_dict) + opp.writeIonmixFile(self.args.outname + ".cn4", **imx_dict) + def convert_tables(): # Grab the input data. input_data = get_input_data() # Read the file extension if the user did not specify an input. - if input_data['args'].input is None: - read_format_ext(input_data['args'], input_data['fn_in']) - - + if input_data["args"].input is None: + read_format_ext(input_data["args"], input_data["fn_in"]) # Reading in data and converting it to the common dictionary format. - eos_dict = Formats_toEosDict(input_data['args'], - input_data['basedir'], - input_data['basename'], - input_data['path_in']).eos_dict + eos_dict = Formats_toEosDict( + input_data["args"], + input_data["basedir"], + input_data["basename"], + input_data["path_in"], + ).eos_dict + + output_type = input_data["args"].output + print(output_type) + if output_type == "sesame": + EosDict_toSesameFile(input_data["args"], eos_dict) + else: + EosDict_toIonmixFile(input_data["args"], eos_dict) - EosDict_toIonmixFile(input_data['args'], eos_dict) -if __name__=='__main__': +if __name__ == "__main__": convert_tables() From 60fcd1341af824025d4ef421bd16de61ba527b29 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Mon, 2 Dec 2024 12:47:44 -0800 Subject: [PATCH 2/4] remove some debug comments --- opacplot2/opg_ionmix.py | 475 ++++++++++++++++++------------ opacplot2/scripts/opac_convert.py | 4 - 2 files changed, 280 insertions(+), 199 deletions(-) diff --git a/opacplot2/opg_ionmix.py b/opacplot2/opg_ionmix.py index 55863b4..f193fea 100644 --- a/opacplot2/opg_ionmix.py +++ b/opacplot2/opg_ionmix.py @@ -11,6 +11,7 @@ from .opl_grid import OplGrid from .constants import ERG_TO_JOULE + class OpacIonmix: """ Class to read in IONMIX EOS and Opacity Files. @@ -136,27 +137,28 @@ class OpacIonmix: setting ``man=True`` may help to fix this. """ - joules_to_ergs = 1.0e+07 - + joules_to_ergs = 1.0e07 def __init__(self, fn, mpi, twot=False, man=False, hassele=False, verbose=False): self.fn = fn self.mpi = mpi self.twot = twot - self.man = man + self.man = man self.hassele = hassele self.verb = verbose - if verbose: print("Reading IONMIX file \"%s\"\n" % (fn)) + if verbose: + print('Reading IONMIX file "%s"\n' % (fn)) - f = open(fn,'r') + f = open(fn, "r") # Read the number of temperatures/densities: self.ntemp = int(f.read(10)) self.ndens = int(f.read(10)) # Skip the next three lines: - for i in range(3): f.readline() + for i in range(3): + f.readline() # Setup temperature/density grid: if self.man == False: @@ -167,13 +169,17 @@ def __init__(self, fn, mpi, twot=False, man=False, hassele=False, verbose=False) self.temp0_log10 = float(f.read(12)) # Compute number densities: - self.numDens = np.logspace(self.dens0_log10, - self.dens0_log10+self.ddens_log10*(self.ndens-1), - self.ndens) - - self.temps = np.logspace(self.temp0_log10, - self.temp0_log10+self.dtemp_log10*(self.ntemp-1), - self.ntemp) + self.numDens = np.logspace( + self.dens0_log10, + self.dens0_log10 + self.ddens_log10 * (self.ndens - 1), + self.ndens, + ) + + self.temps = np.logspace( + self.temp0_log10, + self.temp0_log10 + self.dtemp_log10 * (self.ntemp - 1), + self.ntemp, + ) # Read number of groups: self.ngroups = int(f.read(12)) @@ -183,12 +189,14 @@ def __init__(self, fn, mpi, twot=False, man=False, hassele=False, verbose=False) # Read the rest of the file, remove all of the white space, # and store the string in self.data: - txt = re.sub(r'\s', '', f.read()) - if sys.version < '3': + txt = re.sub(r"\s", "", f.read()) + if sys.version < "3": # converting to unicode if needed for python2 import codecs + def u(x): return codecs.unicode_escape_decode(x)[0] + txt = u(txt) self.data = StringIO(txt) @@ -198,7 +206,7 @@ def u(x): self.temps = self.get_block(self.ntemp) self.numDens = self.get_block(self.ndens) - self.dens = self.numDens * self.mpi/opp.NA + self.dens = self.numDens * self.mpi / opp.NA if self.verb: print(" Number of temperatures: %i" % self.ntemp) @@ -212,7 +220,7 @@ def u(x): self.read_eos() self.read_opac() - def get_block(self,n): + def get_block(self, n): arr = np.zeros(n) for i in range(n): arr[i] = float(self.data.read(12)) @@ -225,31 +233,31 @@ def read_eos(self): nd = self.ndens ng = self.ngroups - self.zbar = self.get_block(nd*nt).reshape(nd,nt) + self.zbar = self.get_block(nd * nt).reshape(nd, nt) if self.twot == False: # Read in e and cv, but convert from J to ergs: - self.etot = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.cvtot = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.dedn = self.get_block(nd*nt).reshape(nd,nt) + self.etot = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.cvtot = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.dedn = self.get_block(nd * nt).reshape(nd, nt) else: # Read in pressure, specific internal energies and # specific heats, but convert from J to ergs: - self.dzdt = self.get_block(nd*nt).reshape(nd,nt) - self.pion = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.pele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.dpidt = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.dpedt = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.eion = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.eele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.cvion = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.cvele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.deidn = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs - self.deedn = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.dzdt = self.get_block(nd * nt).reshape(nd, nt) + self.pion = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.pele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.dpidt = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.dpedt = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.eion = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.eele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.cvion = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.cvele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.deidn = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.deedn = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs if self.hassele: - self.sele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.sele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs def read_opac(self): # Load the opacities from the file. @@ -264,42 +272,53 @@ def read_opac(self): ng = self.ngroups # Read group bounds in eV and convert to ergs: - self.opac_bounds = self.get_block(ng+1) + self.opac_bounds = self.get_block(ng + 1) if self.verb: print("\n Number of Energy Groups: %i" % self.ngroups) - for i in range(0, self.ngroups+1): + for i in range(0, self.ngroups + 1): print("%6i%15.6e" % (i, self.opac_bounds[i])) + self.rosseland = np.empty((nd, nt, ng)) + self.planck_absorb = np.empty((nd, nt, ng)) + self.planck_emiss = np.empty((nd, nt, ng)) - self.rosseland = np.empty((nd,nt,ng)) - self.planck_absorb = np.empty((nd,nt,ng)) - self.planck_emiss = np.empty((nd,nt,ng)) - - arr_ro = self.get_block(nd*nt*ng) - arr_pa = self.get_block(nd*nt*ng) - arr_pe = self.get_block(nd*nt*ng) + arr_ro = self.get_block(nd * nt * ng) + arr_pa = self.get_block(nd * nt * ng) + arr_pe = self.get_block(nd * nt * ng) i = 0 for g in range(ng): for d in range(nd): for t in range(nt): - self.rosseland[d,t,g] = arr_ro[i] - self.planck_absorb[d,t,g] = arr_pa[i] - self.planck_emiss[d,t,g] = arr_pe[i] + self.rosseland[d, t, g] = arr_ro[i] + self.planck_absorb[d, t, g] = arr_pa[i] + self.planck_emiss[d, t, g] = arr_pe[i] i += 1 def oplAbsorb(self): - return OplGrid(self.dens, self.temps, self.opac_bounds, - lambda jd, jt: self.planck_absorb[jd,jt,:]) + return OplGrid( + self.dens, + self.temps, + self.opac_bounds, + lambda jd, jt: self.planck_absorb[jd, jt, :], + ) def oplEmiss(self): - return OplGrid(self.dens, self.temps, self.opac_bounds, - lambda jd, jt: self.planck_emiss[jd,jt,:]) + return OplGrid( + self.dens, + self.temps, + self.opac_bounds, + lambda jd, jt: self.planck_emiss[jd, jt, :], + ) def oplRosseland(self): - return OplGrid(self.dens, self.temps, self.opac_bounds, - lambda jd, jt: self.rosseland[jd,jt,:]) + return OplGrid( + self.dens, + self.temps, + self.opac_bounds, + lambda jd, jt: self.rosseland[jd, jt, :], + ) def write(self, fn, zvals, fracs, twot=None, man=None): """ @@ -327,27 +346,31 @@ def write(self, fn, zvals, fracs, twot=None, man=None): >>> op.extendToZero() # Add temperature point at zero. >>> op.write('imx-0.cn4', (13,), (1,)) """ - if twot is None: twot = self.twot + if twot is None: + twot = self.twot if twot == True and self.twot == False: raise ValueError("Error: Cannot write two-temperature data") - if man is None: man = self.man + if man is None: + man = self.man if man == False and self.man == True: raise ValueError("Error: Cannot write manual temp/dens points") # Write the header: - f = open(fn,'w') - f.write("%10i%10i\n" % (self.ntemp,self.ndens)) + f = open(fn, "w") + f.write("%10i%10i\n" % (self.ntemp, self.ndens)) f.write(" atomic #s of gases: ") - for z in zvals: f.write("%10i" % z) + for z in zvals: + f.write("%10i" % z) f.write("\n relative fractions: ") - for frac in fracs: f.write("%10.2E" % frac) + for frac in fracs: + f.write("%10.2E" % frac) f.write("\n") # Write temperature/density grid and number of groups: def convert(num): string_org = "%12.5E" % (num) - negative = (string_org[0] == "-") + negative = string_org[0] == "-" lead = "-." if negative else "0." string = lead + string_org[1] + string_org[3:8] + "E" @@ -377,7 +400,8 @@ def write_block(var): count = 0 f.write("\n") - if count != 0: f.write("\n") + if count != 0: + f.write("\n") def write_opac_block(var): count = 0 @@ -386,19 +410,25 @@ def write_opac_block(var): for jt in range(self.ntemp): count += 1 - f.write("%s" % convert(var[jd,jt,g])) + f.write("%s" % convert(var[jd, jt, g])) if count == 4: count = 0 f.write("\n") - if count != 0: f.write("\n") + if count != 0: + f.write("\n") if man == False: - f.write("%s%s%s%s" % (convert(self.ddens_log10), - convert(self.dens0_log10), - convert(self.dtemp_log10), - convert(self.temp0_log10)) ) + f.write( + "%s%s%s%s" + % ( + convert(self.ddens_log10), + convert(self.dens0_log10), + convert(self.dtemp_log10), + convert(self.temp0_log10), + ) + ) f.write("%12i\n" % self.ngroups) @@ -409,29 +439,28 @@ def write_opac_block(var): write_block(self.zbar.flatten()) if twot == False: - write_block(self.etot.flatten()/self.joules_to_ergs) - write_block(self.cvtot.flatten()/self.joules_to_ergs) + write_block(self.etot.flatten() / self.joules_to_ergs) + write_block(self.cvtot.flatten() / self.joules_to_ergs) write_block(self.enntab.flatten()) else: - write_block( self.dzdt.flatten()) - write_block( self.pion.flatten()/self.joules_to_ergs) - write_block( self.pele.flatten()/self.joules_to_ergs) - write_block(self.dpidt.flatten()/self.joules_to_ergs) - write_block(self.dpedt.flatten()/self.joules_to_ergs) - write_block(self.eion.flatten()/self.joules_to_ergs) - write_block(self.eele.flatten()/self.joules_to_ergs) - write_block(self.cvion.flatten()/self.joules_to_ergs) - write_block(self.cvele.flatten()/self.joules_to_ergs) - write_block(self.deidn.flatten()/self.joules_to_ergs) - write_block(self.deedn.flatten()/self.joules_to_ergs) + write_block(self.dzdt.flatten()) + write_block(self.pion.flatten() / self.joules_to_ergs) + write_block(self.pele.flatten() / self.joules_to_ergs) + write_block(self.dpidt.flatten() / self.joules_to_ergs) + write_block(self.dpedt.flatten() / self.joules_to_ergs) + write_block(self.eion.flatten() / self.joules_to_ergs) + write_block(self.eele.flatten() / self.joules_to_ergs) + write_block(self.cvion.flatten() / self.joules_to_ergs) + write_block(self.cvele.flatten() / self.joules_to_ergs) + write_block(self.deidn.flatten() / self.joules_to_ergs) + write_block(self.deedn.flatten() / self.joules_to_ergs) write_block(self.opac_bounds) write_opac_block(self.rosseland) write_opac_block(self.planck_absorb) write_opac_block(self.planck_emiss) - def extendToZero(self): """ This routine adds another temperature point at zero. @@ -441,116 +470,132 @@ def extendToZero(self): nt = self.ntemp ng = self.ngroups - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.dzdt[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.dzdt[:, :] self.dzdt = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.pion[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.pion[:, :] self.pion = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.pele[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.pele[:, :] self.pele = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.dpidt[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.dpidt[:, :] self.dpidt = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.dpedt[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.dpedt[:, :] self.dpedi = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.eion[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.eion[:, :] self.eion = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.eele[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.eele[:, :] self.eele = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.zbar[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.zbar[:, :] self.zbar = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.cvion[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.cvion[:, :] self.cvion = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.cvele[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.cvele[:, :] self.cvele = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.deidn[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.deidn[:, :] self.deidn = arr - arr = np.zeros((nd,nt+1)) - arr[:,1:] = self.deedn[:,:] + arr = np.zeros((nd, nt + 1)) + arr[:, 1:] = self.deedn[:, :] self.deedn = arr - arr = np.zeros((nd,nt+1,ng)) - arr[:,1:,:] = self.rosseland[:,:,:] + arr = np.zeros((nd, nt + 1, ng)) + arr[:, 1:, :] = self.rosseland[:, :, :] self.rosseland = arr - arr = np.zeros((nd,nt+1,ng)) - arr[:,1:,:] = self.planck_absorb[:,:,:] + arr = np.zeros((nd, nt + 1, ng)) + arr[:, 1:, :] = self.planck_absorb[:, :, :] self.planck_absorb = arr - arr = np.zeros((nd,nt+1,ng)) - arr[:,1:,:] = self.planck_emiss[:,:,:] + arr = np.zeros((nd, nt + 1, ng)) + arr[:, 1:, :] = self.planck_emiss[:, :, :] self.planck_emiss = arr # Reset temperatures: - arr = np.zeros((nt+1)) + arr = np.zeros((nt + 1)) arr[1:] = self.temps[:] self.temps = arr self.ntemp += 1 - def toEosDict(self, Znum=None,Xnum=None, log=None): + + def toEosDict(self, Znum=None, Xnum=None, log=None): eos_dict = {} if Znum is None: - raise ValueError('Znum Varray should be provided!') + raise ValueError("Znum Varray should be provided!") if Xnum is None: if len(Znum) == 1: Xnum2 = np.array([1.0]) else: - raise ValueError('Xnum array should be provided!') + raise ValueError("Xnum array should be provided!") else: Xnum2 = np.array(Xnum) - eos_dict['idens' ] = self.numDens - eos_dict['temp' ] = self.temps - eos_dict['dens' ] = self.dens - eos_dict['Zf_DT' ] = self.zbar - eos_dict['Ut_DT' ] = self.eion+self.eele#self.etot - eos_dict['Uec_DT' ] = self.eele - eos_dict['Ui_DT' ] = self.eion - eos_dict['Pi_DT' ] = self.pion - eos_dict['Pec_DT' ] = self.pele - eos_dict['Znum' ] = np.array(Znum) - eos_dict['Xnum' ] = Xnum2 - eos_dict['BulkMod'] = 1. - print(self.mpi[0], opp.NA) - eos_dict['Abar' ] = self.mpi[0] - eos_dict['Zmax' ] = np.dot(np.array(Znum), Xnum2) + eos_dict["idens"] = self.numDens + eos_dict["temp"] = self.temps + eos_dict["dens"] = self.dens + eos_dict["Zf_DT"] = self.zbar + eos_dict["Ut_DT"] = self.eion + self.eele # self.etot + eos_dict["Uec_DT"] = self.eele + eos_dict["Ui_DT"] = self.eion + eos_dict["Pi_DT"] = self.pion + eos_dict["Pec_DT"] = self.pele + eos_dict["Znum"] = np.array(Znum) + eos_dict["Xnum"] = Xnum2 + eos_dict["BulkMod"] = 1.0 + eos_dict["Abar"] = self.mpi[0] + eos_dict["Zmax"] = np.dot(np.array(Znum), Xnum2) # mean opacities - eos_dict['opp_int'] = self.planck_absorb[:,:,0] - eos_dict['opr_int'] = self.rosseland[:,:,0] + eos_dict["opp_int"] = self.planck_absorb[:, :, 0] + eos_dict["opr_int"] = self.rosseland[:, :, 0] return eos_dict - - -def writeIonmixFile(fn, zvals, fracs, numDens, temps, - zbar=None, dzdt=None, pion=None, pele=None, - dpidt=None, dpedt=None, eion=None, eele=None, - cvion=None, cvele=None, deidn=None, deedn=None, - ngroups=None, opac_bounds=None, - rosseland=None, planck_absorb=None, planck_emiss=None, - sele=None): +def writeIonmixFile( + fn, + zvals, + fracs, + numDens, + temps, + zbar=None, + dzdt=None, + pion=None, + pele=None, + dpidt=None, + dpedt=None, + eion=None, + eele=None, + cvion=None, + cvele=None, + deidn=None, + deedn=None, + ngroups=None, + opac_bounds=None, + rosseland=None, + planck_absorb=None, + planck_emiss=None, + sele=None, +): """ ``opacplot2.writeIonmixFile()`` provides an explicit and flexible way to write IONMIX files. @@ -623,52 +668,88 @@ def writeIonmixFile(fn, zvals, fracs, numDens, temps, ndens, ntemps = len(numDens), len(temps) - if zbar is None: zbar = np.zeros((ndens,ntemps)) - if dzdt is None: dzdt = np.zeros((ndens,ntemps)) - if pion is None: pion = np.zeros((ndens,ntemps)) - if pele is None: pele = np.zeros((ndens,ntemps)) - if dpidt is None: dpidt = np.zeros((ndens,ntemps)) - if dpedt is None: dpedt = np.zeros((ndens,ntemps)) - if eion is None: eion = np.zeros((ndens,ntemps)) - if eele is None: eele = np.zeros((ndens,ntemps)) - if cvion is None: cvion = np.zeros((ndens,ntemps)) - if cvele is None: cvele = np.zeros((ndens,ntemps)) - if deidn is None: deidn = np.zeros((ndens,ntemps)) - if deedn is None: deedn = np.zeros((ndens,ntemps)) + if zbar is None: + zbar = np.zeros((ndens, ntemps)) + if dzdt is None: + dzdt = np.zeros((ndens, ntemps)) + if pion is None: + pion = np.zeros((ndens, ntemps)) + if pele is None: + pele = np.zeros((ndens, ntemps)) + if dpidt is None: + dpidt = np.zeros((ndens, ntemps)) + if dpedt is None: + dpedt = np.zeros((ndens, ntemps)) + if eion is None: + eion = np.zeros((ndens, ntemps)) + if eele is None: + eele = np.zeros((ndens, ntemps)) + if cvion is None: + cvion = np.zeros((ndens, ntemps)) + if cvele is None: + cvele = np.zeros((ndens, ntemps)) + if deidn is None: + deidn = np.zeros((ndens, ntemps)) + if deedn is None: + deedn = np.zeros((ndens, ntemps)) if ngroups is None: ngroups = 0 - if opac_bounds is None: opac_bounds = (0.0,1.0) - if rosseland is None: rosseland = np.zeros((ndens,ntemps,ngroups)) - if planck_absorb is None: planck_absorb = np.zeros((ndens,ntemps,ngroups)) - if planck_emiss is None: planck_emiss = np.zeros((ndens,ntemps,ngroups)) - - for tab in ["zbar", "dzdt", "pion", "pele", "dpidt", "dpedt", "eion", - "eele", "cvion", "cvele", "deidn", "deedn"]: + if opac_bounds is None: + opac_bounds = (0.0, 1.0) + if rosseland is None: + rosseland = np.zeros((ndens, ntemps, ngroups)) + if planck_absorb is None: + planck_absorb = np.zeros((ndens, ntemps, ngroups)) + if planck_emiss is None: + planck_emiss = np.zeros((ndens, ntemps, ngroups)) + + for tab in [ + "zbar", + "dzdt", + "pion", + "pele", + "dpidt", + "dpedt", + "eion", + "eele", + "cvion", + "cvele", + "deidn", + "deedn", + ]: ctab = locals()[tab] if ctab.shape != (ndens, ntemps): - raise ValueError('Table {0} has shape {1}, expected {2}!'.format( - tab, str(ctab.shape), str((ndens, ntemps)))) - for tab in ['rosseland', 'planck_absorb', 'planck_emiss']: + raise ValueError( + "Table {0} has shape {1}, expected {2}!".format( + tab, str(ctab.shape), str((ndens, ntemps)) + ) + ) + for tab in ["rosseland", "planck_absorb", "planck_emiss"]: ctab = locals()[tab] if ctab.shape != (ndens, ntemps, ngroups): - raise ValueError('Table {0} has shape {1}, expected {2}!'.format( - tab, str(ctab.shape), str((ndens, ntemps, ngroups)))) + raise ValueError( + "Table {0} has shape {1}, expected {2}!".format( + tab, str(ctab.shape), str((ndens, ntemps, ngroups)) + ) + ) # Write the header: - f = open(fn,'w') - f.write("%10i%10i\n" % (ntemps,ndens)) + f = open(fn, "w") + f.write("%10i%10i\n" % (ntemps, ndens)) f.write(" atomic #s of gases: ") - for z in zvals: f.write("%10i" % z) + for z in zvals: + f.write("%10i" % z) f.write("\n relative fractions: ") - for frac in fracs: f.write("%10.2E" % frac) + for frac in fracs: + f.write("%10.2E" % frac) f.write("\n") # Write temperature/density grid and number of groups: # name argument is for error reporting purposes. def convert(num, name): string_org = "%12.5E" % (num) - negative = (string_org[0] == "-") + negative = string_org[0] == "-" lead = "-." if negative else "0." string = lead + string_org[1] + string_org[3:8] + "E" @@ -679,9 +760,11 @@ def convert(num, name): if int(string_org[1] + string_org[3:8]) == 0: return string + "+00" except ValueError: - raise ValueError('There was a problem writing the data in ' - 'the {} block to IONMIX. Try writing it in ' - 'log format.'.format(name)) + raise ValueError( + "There was a problem writing the data in " + "the {} block to IONMIX. Try writing it in " + "log format.".format(name) + ) # Not zero: expo = int(string_org[9:]) + 1 @@ -703,7 +786,8 @@ def write_block(var, name): count = 0 f.write("\n") - if count != 0: f.write("\n") + if count != 0: + f.write("\n") def write_opac_block(var, name): count = 0 @@ -712,39 +796,40 @@ def write_opac_block(var, name): for jt in range(ntemps): count += 1 - f.write("%s" % convert(var[jd,jt,g], name)) + f.write("%s" % convert(var[jd, jt, g], name)) if count == 4: count = 0 f.write("\n") - if count != 0: f.write("\n") + if count != 0: + f.write("\n") f.write("%12i\n" % ngroups) - write_block(temps, 'temperature') - write_block(numDens, 'number density') + write_block(temps, "temperature") + write_block(numDens, "number density") - write_block(zbar.flatten(), 'average ionization') + write_block(zbar.flatten(), "average ionization") - write_block(dzdt.flatten(), 'DZ_DT') - write_block(pion.flatten()*ERG_TO_JOULE, 'ion pressure') - write_block(pele.flatten()*ERG_TO_JOULE, 'electron pressure') - write_block(dpidt.flatten()*ERG_TO_JOULE, 'D(ion pressure)_DT') - write_block(dpedt.flatten()*ERG_TO_JOULE, 'D(electron pressure)_DT') - print (eion[0,0], eion[0,0]*ERG_TO_JOULE, ERG_TO_JOULE) - write_block(eion.flatten()*ERG_TO_JOULE, 'ion energy') - write_block(eele.flatten()*ERG_TO_JOULE, 'electron energy') - write_block(cvion.flatten()*ERG_TO_JOULE, 'ion CV') - write_block(cvele.flatten()*ERG_TO_JOULE, 'electron CV') - write_block(deidn.flatten()*ERG_TO_JOULE, 'D(ion energy)_DT') - write_block(deedn.flatten()*ERG_TO_JOULE, 'D(electron energy)_DT') + write_block(dzdt.flatten(), "DZ_DT") + write_block(pion.flatten() * ERG_TO_JOULE, "ion pressure") + write_block(pele.flatten() * ERG_TO_JOULE, "electron pressure") + write_block(dpidt.flatten() * ERG_TO_JOULE, "D(ion pressure)_DT") + write_block(dpedt.flatten() * ERG_TO_JOULE, "D(electron pressure)_DT") + print(eion[0, 0], eion[0, 0] * ERG_TO_JOULE, ERG_TO_JOULE) + write_block(eion.flatten() * ERG_TO_JOULE, "ion energy") + write_block(eele.flatten() * ERG_TO_JOULE, "electron energy") + write_block(cvion.flatten() * ERG_TO_JOULE, "ion CV") + write_block(cvele.flatten() * ERG_TO_JOULE, "electron CV") + write_block(deidn.flatten() * ERG_TO_JOULE, "D(ion energy)_DT") + write_block(deedn.flatten() * ERG_TO_JOULE, "D(electron energy)_DT") # Check for electron entropy (if it is there): - if sele is not None: write_block(sele.flatten()*ERG_TO_JOULE, - 'electron entropy') + if sele is not None: + write_block(sele.flatten() * ERG_TO_JOULE, "electron entropy") - write_block(opac_bounds, 'opacity bounds') - write_opac_block(rosseland, 'rosseland opacity') - write_opac_block(planck_absorb, 'planck absorption') - write_opac_block(planck_emiss, 'planck emissivity') + write_block(opac_bounds, "opacity bounds") + write_opac_block(rosseland, "rosseland opacity") + write_opac_block(planck_absorb, "planck absorption") + write_opac_block(planck_emiss, "planck emissivity") diff --git a/opacplot2/scripts/opac_convert.py b/opacplot2/scripts/opac_convert.py index e7a414e..460bebf 100644 --- a/opacplot2/scripts/opac_convert.py +++ b/opacplot2/scripts/opac_convert.py @@ -199,7 +199,6 @@ def multi_toEosDict(self): return eos_dict def ionmix_toEosDict(self): - print(self.args.mpi) op = opp.OpacIonmix(self.path_in, self.args.mpi, man=True, twot=True) eos_dict = op.toEosDict( Znum=self.args.Znum, Xnum=self.args.Xfracs, log=self.args.log @@ -304,9 +303,7 @@ def eosDict_toSesame(self): xx = self.eos_dict["Xnum"] znum = 0.0 for i in range(len(self.eos_dict["Znum"])): - print(zz[i]) znum += zz[i] * xx[i] - print(znum) else: znum = self.eos_dict["Znum"][0] @@ -451,7 +448,6 @@ def convert_tables(): ).eos_dict output_type = input_data["args"].output - print(output_type) if output_type == "sesame": EosDict_toSesameFile(input_data["args"], eos_dict) else: From aa1fb5beb82431831a8287ba169172393e299cc0 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Mon, 2 Dec 2024 13:26:31 -0800 Subject: [PATCH 3/4] reverting some auto-formatting from black --- opacplot2/opg_ionmix.py | 473 ++++++++++++------------------ opacplot2/opg_sesame.py | 366 ++++++++++------------- opacplot2/scripts/opac_convert.py | 333 +++++++++------------ 3 files changed, 499 insertions(+), 673 deletions(-) diff --git a/opacplot2/opg_ionmix.py b/opacplot2/opg_ionmix.py index f193fea..ee32ca7 100644 --- a/opacplot2/opg_ionmix.py +++ b/opacplot2/opg_ionmix.py @@ -11,7 +11,6 @@ from .opl_grid import OplGrid from .constants import ERG_TO_JOULE - class OpacIonmix: """ Class to read in IONMIX EOS and Opacity Files. @@ -137,28 +136,27 @@ class OpacIonmix: setting ``man=True`` may help to fix this. """ - joules_to_ergs = 1.0e07 + joules_to_ergs = 1.0e+07 + def __init__(self, fn, mpi, twot=False, man=False, hassele=False, verbose=False): self.fn = fn self.mpi = mpi self.twot = twot - self.man = man + self.man = man self.hassele = hassele self.verb = verbose - if verbose: - print('Reading IONMIX file "%s"\n' % (fn)) + if verbose: print("Reading IONMIX file \"%s\"\n" % (fn)) - f = open(fn, "r") + f = open(fn,'r') # Read the number of temperatures/densities: self.ntemp = int(f.read(10)) self.ndens = int(f.read(10)) # Skip the next three lines: - for i in range(3): - f.readline() + for i in range(3): f.readline() # Setup temperature/density grid: if self.man == False: @@ -169,17 +167,13 @@ def __init__(self, fn, mpi, twot=False, man=False, hassele=False, verbose=False) self.temp0_log10 = float(f.read(12)) # Compute number densities: - self.numDens = np.logspace( - self.dens0_log10, - self.dens0_log10 + self.ddens_log10 * (self.ndens - 1), - self.ndens, - ) - - self.temps = np.logspace( - self.temp0_log10, - self.temp0_log10 + self.dtemp_log10 * (self.ntemp - 1), - self.ntemp, - ) + self.numDens = np.logspace(self.dens0_log10, + self.dens0_log10+self.ddens_log10*(self.ndens-1), + self.ndens) + + self.temps = np.logspace(self.temp0_log10, + self.temp0_log10+self.dtemp_log10*(self.ntemp-1), + self.ntemp) # Read number of groups: self.ngroups = int(f.read(12)) @@ -189,14 +183,12 @@ def __init__(self, fn, mpi, twot=False, man=False, hassele=False, verbose=False) # Read the rest of the file, remove all of the white space, # and store the string in self.data: - txt = re.sub(r"\s", "", f.read()) - if sys.version < "3": + txt = re.sub(r'\s', '', f.read()) + if sys.version < '3': # converting to unicode if needed for python2 import codecs - def u(x): return codecs.unicode_escape_decode(x)[0] - txt = u(txt) self.data = StringIO(txt) @@ -206,7 +198,7 @@ def u(x): self.temps = self.get_block(self.ntemp) self.numDens = self.get_block(self.ndens) - self.dens = self.numDens * self.mpi / opp.NA + self.dens = self.numDens * self.mpi/opp.NA if self.verb: print(" Number of temperatures: %i" % self.ntemp) @@ -220,7 +212,7 @@ def u(x): self.read_eos() self.read_opac() - def get_block(self, n): + def get_block(self,n): arr = np.zeros(n) for i in range(n): arr[i] = float(self.data.read(12)) @@ -233,31 +225,31 @@ def read_eos(self): nd = self.ndens ng = self.ngroups - self.zbar = self.get_block(nd * nt).reshape(nd, nt) + self.zbar = self.get_block(nd*nt).reshape(nd,nt) if self.twot == False: # Read in e and cv, but convert from J to ergs: - self.etot = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.cvtot = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.dedn = self.get_block(nd * nt).reshape(nd, nt) + self.etot = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.cvtot = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.dedn = self.get_block(nd*nt).reshape(nd,nt) else: # Read in pressure, specific internal energies and # specific heats, but convert from J to ergs: - self.dzdt = self.get_block(nd * nt).reshape(nd, nt) - self.pion = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.pele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.dpidt = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.dpedt = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.eion = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.eele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.cvion = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.cvele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.deidn = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs - self.deedn = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.dzdt = self.get_block(nd*nt).reshape(nd,nt) + self.pion = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.pele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.dpidt = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.dpedt = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.eion = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.eele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.cvion = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.cvele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.deidn = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs + self.deedn = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs if self.hassele: - self.sele = self.get_block(nd * nt).reshape(nd, nt) * self.joules_to_ergs + self.sele = self.get_block(nd*nt).reshape(nd,nt) * self.joules_to_ergs def read_opac(self): # Load the opacities from the file. @@ -272,53 +264,42 @@ def read_opac(self): ng = self.ngroups # Read group bounds in eV and convert to ergs: - self.opac_bounds = self.get_block(ng + 1) + self.opac_bounds = self.get_block(ng+1) if self.verb: print("\n Number of Energy Groups: %i" % self.ngroups) - for i in range(0, self.ngroups + 1): + for i in range(0, self.ngroups+1): print("%6i%15.6e" % (i, self.opac_bounds[i])) - self.rosseland = np.empty((nd, nt, ng)) - self.planck_absorb = np.empty((nd, nt, ng)) - self.planck_emiss = np.empty((nd, nt, ng)) - arr_ro = self.get_block(nd * nt * ng) - arr_pa = self.get_block(nd * nt * ng) - arr_pe = self.get_block(nd * nt * ng) + self.rosseland = np.empty((nd,nt,ng)) + self.planck_absorb = np.empty((nd,nt,ng)) + self.planck_emiss = np.empty((nd,nt,ng)) + + arr_ro = self.get_block(nd*nt*ng) + arr_pa = self.get_block(nd*nt*ng) + arr_pe = self.get_block(nd*nt*ng) i = 0 for g in range(ng): for d in range(nd): for t in range(nt): - self.rosseland[d, t, g] = arr_ro[i] - self.planck_absorb[d, t, g] = arr_pa[i] - self.planck_emiss[d, t, g] = arr_pe[i] + self.rosseland[d,t,g] = arr_ro[i] + self.planck_absorb[d,t,g] = arr_pa[i] + self.planck_emiss[d,t,g] = arr_pe[i] i += 1 def oplAbsorb(self): - return OplGrid( - self.dens, - self.temps, - self.opac_bounds, - lambda jd, jt: self.planck_absorb[jd, jt, :], - ) + return OplGrid(self.dens, self.temps, self.opac_bounds, + lambda jd, jt: self.planck_absorb[jd,jt,:]) def oplEmiss(self): - return OplGrid( - self.dens, - self.temps, - self.opac_bounds, - lambda jd, jt: self.planck_emiss[jd, jt, :], - ) + return OplGrid(self.dens, self.temps, self.opac_bounds, + lambda jd, jt: self.planck_emiss[jd,jt,:]) def oplRosseland(self): - return OplGrid( - self.dens, - self.temps, - self.opac_bounds, - lambda jd, jt: self.rosseland[jd, jt, :], - ) + return OplGrid(self.dens, self.temps, self.opac_bounds, + lambda jd, jt: self.rosseland[jd,jt,:]) def write(self, fn, zvals, fracs, twot=None, man=None): """ @@ -346,31 +327,27 @@ def write(self, fn, zvals, fracs, twot=None, man=None): >>> op.extendToZero() # Add temperature point at zero. >>> op.write('imx-0.cn4', (13,), (1,)) """ - if twot is None: - twot = self.twot + if twot is None: twot = self.twot if twot == True and self.twot == False: raise ValueError("Error: Cannot write two-temperature data") - if man is None: - man = self.man + if man is None: man = self.man if man == False and self.man == True: raise ValueError("Error: Cannot write manual temp/dens points") # Write the header: - f = open(fn, "w") - f.write("%10i%10i\n" % (self.ntemp, self.ndens)) + f = open(fn,'w') + f.write("%10i%10i\n" % (self.ntemp,self.ndens)) f.write(" atomic #s of gases: ") - for z in zvals: - f.write("%10i" % z) + for z in zvals: f.write("%10i" % z) f.write("\n relative fractions: ") - for frac in fracs: - f.write("%10.2E" % frac) + for frac in fracs: f.write("%10.2E" % frac) f.write("\n") # Write temperature/density grid and number of groups: def convert(num): string_org = "%12.5E" % (num) - negative = string_org[0] == "-" + negative = (string_org[0] == "-") lead = "-." if negative else "0." string = lead + string_org[1] + string_org[3:8] + "E" @@ -400,8 +377,7 @@ def write_block(var): count = 0 f.write("\n") - if count != 0: - f.write("\n") + if count != 0: f.write("\n") def write_opac_block(var): count = 0 @@ -410,25 +386,19 @@ def write_opac_block(var): for jt in range(self.ntemp): count += 1 - f.write("%s" % convert(var[jd, jt, g])) + f.write("%s" % convert(var[jd,jt,g])) if count == 4: count = 0 f.write("\n") - if count != 0: - f.write("\n") + if count != 0: f.write("\n") if man == False: - f.write( - "%s%s%s%s" - % ( - convert(self.ddens_log10), - convert(self.dens0_log10), - convert(self.dtemp_log10), - convert(self.temp0_log10), - ) - ) + f.write("%s%s%s%s" % (convert(self.ddens_log10), + convert(self.dens0_log10), + convert(self.dtemp_log10), + convert(self.temp0_log10)) ) f.write("%12i\n" % self.ngroups) @@ -439,28 +409,29 @@ def write_opac_block(var): write_block(self.zbar.flatten()) if twot == False: - write_block(self.etot.flatten() / self.joules_to_ergs) - write_block(self.cvtot.flatten() / self.joules_to_ergs) + write_block(self.etot.flatten()/self.joules_to_ergs) + write_block(self.cvtot.flatten()/self.joules_to_ergs) write_block(self.enntab.flatten()) else: - write_block(self.dzdt.flatten()) - write_block(self.pion.flatten() / self.joules_to_ergs) - write_block(self.pele.flatten() / self.joules_to_ergs) - write_block(self.dpidt.flatten() / self.joules_to_ergs) - write_block(self.dpedt.flatten() / self.joules_to_ergs) - write_block(self.eion.flatten() / self.joules_to_ergs) - write_block(self.eele.flatten() / self.joules_to_ergs) - write_block(self.cvion.flatten() / self.joules_to_ergs) - write_block(self.cvele.flatten() / self.joules_to_ergs) - write_block(self.deidn.flatten() / self.joules_to_ergs) - write_block(self.deedn.flatten() / self.joules_to_ergs) + write_block( self.dzdt.flatten()) + write_block( self.pion.flatten()/self.joules_to_ergs) + write_block( self.pele.flatten()/self.joules_to_ergs) + write_block(self.dpidt.flatten()/self.joules_to_ergs) + write_block(self.dpedt.flatten()/self.joules_to_ergs) + write_block(self.eion.flatten()/self.joules_to_ergs) + write_block(self.eele.flatten()/self.joules_to_ergs) + write_block(self.cvion.flatten()/self.joules_to_ergs) + write_block(self.cvele.flatten()/self.joules_to_ergs) + write_block(self.deidn.flatten()/self.joules_to_ergs) + write_block(self.deedn.flatten()/self.joules_to_ergs) write_block(self.opac_bounds) write_opac_block(self.rosseland) write_opac_block(self.planck_absorb) write_opac_block(self.planck_emiss) + def extendToZero(self): """ This routine adds another temperature point at zero. @@ -470,132 +441,115 @@ def extendToZero(self): nt = self.ntemp ng = self.ngroups - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.dzdt[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.dzdt[:,:] self.dzdt = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.pion[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.pion[:,:] self.pion = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.pele[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.pele[:,:] self.pele = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.dpidt[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.dpidt[:,:] self.dpidt = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.dpedt[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.dpedt[:,:] self.dpedi = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.eion[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.eion[:,:] self.eion = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.eele[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.eele[:,:] self.eele = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.zbar[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.zbar[:,:] self.zbar = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.cvion[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.cvion[:,:] self.cvion = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.cvele[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.cvele[:,:] self.cvele = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.deidn[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.deidn[:,:] self.deidn = arr - arr = np.zeros((nd, nt + 1)) - arr[:, 1:] = self.deedn[:, :] + arr = np.zeros((nd,nt+1)) + arr[:,1:] = self.deedn[:,:] self.deedn = arr - arr = np.zeros((nd, nt + 1, ng)) - arr[:, 1:, :] = self.rosseland[:, :, :] + arr = np.zeros((nd,nt+1,ng)) + arr[:,1:,:] = self.rosseland[:,:,:] self.rosseland = arr - arr = np.zeros((nd, nt + 1, ng)) - arr[:, 1:, :] = self.planck_absorb[:, :, :] + arr = np.zeros((nd,nt+1,ng)) + arr[:,1:,:] = self.planck_absorb[:,:,:] self.planck_absorb = arr - arr = np.zeros((nd, nt + 1, ng)) - arr[:, 1:, :] = self.planck_emiss[:, :, :] + arr = np.zeros((nd,nt+1,ng)) + arr[:,1:,:] = self.planck_emiss[:,:,:] self.planck_emiss = arr # Reset temperatures: - arr = np.zeros((nt + 1)) + arr = np.zeros((nt+1)) arr[1:] = self.temps[:] self.temps = arr self.ntemp += 1 - - def toEosDict(self, Znum=None, Xnum=None, log=None): + def toEosDict(self, Znum=None,Xnum=None, log=None): eos_dict = {} if Znum is None: - raise ValueError("Znum Varray should be provided!") + raise ValueError('Znum Varray should be provided!') if Xnum is None: if len(Znum) == 1: Xnum2 = np.array([1.0]) else: - raise ValueError("Xnum array should be provided!") + raise ValueError('Xnum array should be provided!') else: Xnum2 = np.array(Xnum) - eos_dict["idens"] = self.numDens - eos_dict["temp"] = self.temps - eos_dict["dens"] = self.dens - eos_dict["Zf_DT"] = self.zbar - eos_dict["Ut_DT"] = self.eion + self.eele # self.etot - eos_dict["Uec_DT"] = self.eele - eos_dict["Ui_DT"] = self.eion - eos_dict["Pi_DT"] = self.pion - eos_dict["Pec_DT"] = self.pele - eos_dict["Znum"] = np.array(Znum) - eos_dict["Xnum"] = Xnum2 - eos_dict["BulkMod"] = 1.0 - eos_dict["Abar"] = self.mpi[0] - eos_dict["Zmax"] = np.dot(np.array(Znum), Xnum2) + eos_dict['idens' ] = self.numDens + eos_dict['temp' ] = self.temps + eos_dict['dens' ] = self.dens + eos_dict['Zf_DT' ] = self.zbar + eos_dict['Ut_DT' ] = self.eion+self.eele#self.etot + eos_dict['Uec_DT' ] = self.eele + eos_dict['Ui_DT' ] = self.eion + eos_dict['Pi_DT' ] = self.pion + eos_dict['Pec_DT' ] = self.pele + eos_dict['Znum' ] = np.array(Znum) + eos_dict['Xnum' ] = Xnum2 + eos_dict['BulkMod'] = 1. + eos_dict['Abar' ] = self.mpi[0] + eos_dict['Zmax' ] = np.dot(np.array(Znum), Xnum2) # mean opacities - eos_dict["opp_int"] = self.planck_absorb[:, :, 0] - eos_dict["opr_int"] = self.rosseland[:, :, 0] + eos_dict['opp_int'] = self.planck_absorb[:,:,0] + eos_dict['opr_int'] = self.rosseland[:,:,0] return eos_dict -def writeIonmixFile( - fn, - zvals, - fracs, - numDens, - temps, - zbar=None, - dzdt=None, - pion=None, - pele=None, - dpidt=None, - dpedt=None, - eion=None, - eele=None, - cvion=None, - cvele=None, - deidn=None, - deedn=None, - ngroups=None, - opac_bounds=None, - rosseland=None, - planck_absorb=None, - planck_emiss=None, - sele=None, -): + + +def writeIonmixFile(fn, zvals, fracs, numDens, temps, + zbar=None, dzdt=None, pion=None, pele=None, + dpidt=None, dpedt=None, eion=None, eele=None, + cvion=None, cvele=None, deidn=None, deedn=None, + ngroups=None, opac_bounds=None, + rosseland=None, planck_absorb=None, planck_emiss=None, + sele=None): """ ``opacplot2.writeIonmixFile()`` provides an explicit and flexible way to write IONMIX files. @@ -668,88 +622,52 @@ def writeIonmixFile( ndens, ntemps = len(numDens), len(temps) - if zbar is None: - zbar = np.zeros((ndens, ntemps)) - if dzdt is None: - dzdt = np.zeros((ndens, ntemps)) - if pion is None: - pion = np.zeros((ndens, ntemps)) - if pele is None: - pele = np.zeros((ndens, ntemps)) - if dpidt is None: - dpidt = np.zeros((ndens, ntemps)) - if dpedt is None: - dpedt = np.zeros((ndens, ntemps)) - if eion is None: - eion = np.zeros((ndens, ntemps)) - if eele is None: - eele = np.zeros((ndens, ntemps)) - if cvion is None: - cvion = np.zeros((ndens, ntemps)) - if cvele is None: - cvele = np.zeros((ndens, ntemps)) - if deidn is None: - deidn = np.zeros((ndens, ntemps)) - if deedn is None: - deedn = np.zeros((ndens, ntemps)) + if zbar is None: zbar = np.zeros((ndens,ntemps)) + if dzdt is None: dzdt = np.zeros((ndens,ntemps)) + if pion is None: pion = np.zeros((ndens,ntemps)) + if pele is None: pele = np.zeros((ndens,ntemps)) + if dpidt is None: dpidt = np.zeros((ndens,ntemps)) + if dpedt is None: dpedt = np.zeros((ndens,ntemps)) + if eion is None: eion = np.zeros((ndens,ntemps)) + if eele is None: eele = np.zeros((ndens,ntemps)) + if cvion is None: cvion = np.zeros((ndens,ntemps)) + if cvele is None: cvele = np.zeros((ndens,ntemps)) + if deidn is None: deidn = np.zeros((ndens,ntemps)) + if deedn is None: deedn = np.zeros((ndens,ntemps)) if ngroups is None: ngroups = 0 - if opac_bounds is None: - opac_bounds = (0.0, 1.0) - if rosseland is None: - rosseland = np.zeros((ndens, ntemps, ngroups)) - if planck_absorb is None: - planck_absorb = np.zeros((ndens, ntemps, ngroups)) - if planck_emiss is None: - planck_emiss = np.zeros((ndens, ntemps, ngroups)) - - for tab in [ - "zbar", - "dzdt", - "pion", - "pele", - "dpidt", - "dpedt", - "eion", - "eele", - "cvion", - "cvele", - "deidn", - "deedn", - ]: + if opac_bounds is None: opac_bounds = (0.0,1.0) + if rosseland is None: rosseland = np.zeros((ndens,ntemps,ngroups)) + if planck_absorb is None: planck_absorb = np.zeros((ndens,ntemps,ngroups)) + if planck_emiss is None: planck_emiss = np.zeros((ndens,ntemps,ngroups)) + + for tab in ["zbar", "dzdt", "pion", "pele", "dpidt", "dpedt", "eion", + "eele", "cvion", "cvele", "deidn", "deedn"]: ctab = locals()[tab] if ctab.shape != (ndens, ntemps): - raise ValueError( - "Table {0} has shape {1}, expected {2}!".format( - tab, str(ctab.shape), str((ndens, ntemps)) - ) - ) - for tab in ["rosseland", "planck_absorb", "planck_emiss"]: + raise ValueError('Table {0} has shape {1}, expected {2}!'.format( + tab, str(ctab.shape), str((ndens, ntemps)))) + for tab in ['rosseland', 'planck_absorb', 'planck_emiss']: ctab = locals()[tab] if ctab.shape != (ndens, ntemps, ngroups): - raise ValueError( - "Table {0} has shape {1}, expected {2}!".format( - tab, str(ctab.shape), str((ndens, ntemps, ngroups)) - ) - ) + raise ValueError('Table {0} has shape {1}, expected {2}!'.format( + tab, str(ctab.shape), str((ndens, ntemps, ngroups)))) # Write the header: - f = open(fn, "w") - f.write("%10i%10i\n" % (ntemps, ndens)) + f = open(fn,'w') + f.write("%10i%10i\n" % (ntemps,ndens)) f.write(" atomic #s of gases: ") - for z in zvals: - f.write("%10i" % z) + for z in zvals: f.write("%10i" % z) f.write("\n relative fractions: ") - for frac in fracs: - f.write("%10.2E" % frac) + for frac in fracs: f.write("%10.2E" % frac) f.write("\n") # Write temperature/density grid and number of groups: # name argument is for error reporting purposes. def convert(num, name): string_org = "%12.5E" % (num) - negative = string_org[0] == "-" + negative = (string_org[0] == "-") lead = "-." if negative else "0." string = lead + string_org[1] + string_org[3:8] + "E" @@ -760,11 +678,9 @@ def convert(num, name): if int(string_org[1] + string_org[3:8]) == 0: return string + "+00" except ValueError: - raise ValueError( - "There was a problem writing the data in " - "the {} block to IONMIX. Try writing it in " - "log format.".format(name) - ) + raise ValueError('There was a problem writing the data in ' + 'the {} block to IONMIX. Try writing it in ' + 'log format.'.format(name)) # Not zero: expo = int(string_org[9:]) + 1 @@ -786,8 +702,7 @@ def write_block(var, name): count = 0 f.write("\n") - if count != 0: - f.write("\n") + if count != 0: f.write("\n") def write_opac_block(var, name): count = 0 @@ -796,40 +711,38 @@ def write_opac_block(var, name): for jt in range(ntemps): count += 1 - f.write("%s" % convert(var[jd, jt, g], name)) + f.write("%s" % convert(var[jd,jt,g], name)) if count == 4: count = 0 f.write("\n") - if count != 0: - f.write("\n") + if count != 0: f.write("\n") f.write("%12i\n" % ngroups) - write_block(temps, "temperature") - write_block(numDens, "number density") + write_block(temps, 'temperature') + write_block(numDens, 'number density') - write_block(zbar.flatten(), "average ionization") + write_block(zbar.flatten(), 'average ionization') - write_block(dzdt.flatten(), "DZ_DT") - write_block(pion.flatten() * ERG_TO_JOULE, "ion pressure") - write_block(pele.flatten() * ERG_TO_JOULE, "electron pressure") - write_block(dpidt.flatten() * ERG_TO_JOULE, "D(ion pressure)_DT") - write_block(dpedt.flatten() * ERG_TO_JOULE, "D(electron pressure)_DT") - print(eion[0, 0], eion[0, 0] * ERG_TO_JOULE, ERG_TO_JOULE) - write_block(eion.flatten() * ERG_TO_JOULE, "ion energy") - write_block(eele.flatten() * ERG_TO_JOULE, "electron energy") - write_block(cvion.flatten() * ERG_TO_JOULE, "ion CV") - write_block(cvele.flatten() * ERG_TO_JOULE, "electron CV") - write_block(deidn.flatten() * ERG_TO_JOULE, "D(ion energy)_DT") - write_block(deedn.flatten() * ERG_TO_JOULE, "D(electron energy)_DT") + write_block(dzdt.flatten(), 'DZ_DT') + write_block(pion.flatten()*ERG_TO_JOULE, 'ion pressure') + write_block(pele.flatten()*ERG_TO_JOULE, 'electron pressure') + write_block(dpidt.flatten()*ERG_TO_JOULE, 'D(ion pressure)_DT') + write_block(dpedt.flatten()*ERG_TO_JOULE, 'D(electron pressure)_DT') + write_block(eion.flatten()*ERG_TO_JOULE, 'ion energy') + write_block(eele.flatten()*ERG_TO_JOULE, 'electron energy') + write_block(cvion.flatten()*ERG_TO_JOULE, 'ion CV') + write_block(cvele.flatten()*ERG_TO_JOULE, 'electron CV') + write_block(deidn.flatten()*ERG_TO_JOULE, 'D(ion energy)_DT') + write_block(deedn.flatten()*ERG_TO_JOULE, 'D(electron energy)_DT') # Check for electron entropy (if it is there): - if sele is not None: - write_block(sele.flatten() * ERG_TO_JOULE, "electron entropy") + if sele is not None: write_block(sele.flatten()*ERG_TO_JOULE, + 'electron entropy') - write_block(opac_bounds, "opacity bounds") - write_opac_block(rosseland, "rosseland opacity") - write_opac_block(planck_absorb, "planck absorption") - write_opac_block(planck_emiss, "planck emissivity") + write_block(opac_bounds, 'opacity bounds') + write_opac_block(rosseland, 'rosseland opacity') + write_opac_block(planck_absorb, 'planck absorption') + write_opac_block(planck_emiss, 'planck emissivity') diff --git a/opacplot2/opg_sesame.py b/opacplot2/opg_sesame.py index e516a1b..a5c821c 100644 --- a/opacplot2/opg_sesame.py +++ b/opacplot2/opg_sesame.py @@ -1,8 +1,7 @@ from __future__ import absolute_import from __future__ import division from __future__ import print_function - -# from __future__ import unicode_literals +#from __future__ import unicode_literals from io import open import six @@ -12,17 +11,10 @@ import opacplot2.utils import numpy as np -from .constants import ( - KELVIN_TO_EV, - GPA_TO_ERGCC, - MJKG_TO_ERGCC, - ERGCC_TO_GPA, - ERGG_TO_MJKG, -) +from .constants import KELVIN_TO_EV, GPA_TO_ERGCC, MJKG_TO_ERGCC, ERGCC_TO_GPA, ERGG_TO_MJKG import periodictable as ptab - class OpgSesame: """ This class is responsible for loading all SESAME formatted data files. @@ -70,39 +62,39 @@ class OpgSesame: def __init__(self, filename, precision, verbose=False): self.verbose = verbose - self.fhand = open(filename, encoding="utf-8") - if precision == self.SINGLE: - self.entry_len = 15 - elif precision == self.DOUBLE: + self.fhand = open(filename, encoding='utf-8') + if(precision == self.SINGLE): + self.entry_len = 15 + elif(precision == self.DOUBLE): self.entry_len = 22 else: raise ValueError("precision must be SINGLE or DOUBLE") - self.fdict = { - 101: self.parseComment, - 102: self.parseComment, - 103: self.parseComment, - 104: self.parseComment, - 201: self.parseInfo, - 301: self.parseEos, - 303: self.parseEos, - 304: self.parseEos, - 305: self.parseEos, - 306: self.parseEos, - 311: self.parseEos, - 321: self.parseEos, - 401: self.parseVap, - 411: self.parseSolid, - 412: self.parseLiquid, - 431: self.parseShear, - 432: self.parseShear, - 504: self.parseZbar, - 601: self.parseZbar, - 602: self.parseEcond, - 603: self.parseTcond, - 604: self.parseTcond, - 605: self.parseTcond, - } + + self.fdict = { 101 : self.parseComment, + 102 : self.parseComment, + 103 : self.parseComment, + 104 : self.parseComment, + 201 : self.parseInfo, + 301 : self.parseEos, + 303 : self.parseEos, + 304 : self.parseEos, + 305 : self.parseEos, + 306 : self.parseEos, + 311 : self.parseEos, + 321 : self.parseEos, + 401 : self.parseVap, + 411 : self.parseSolid, + 412 : self.parseLiquid, + 431 : self.parseShear, + 432 : self.parseShear, + 504 : self.parseZbar, + 601 : self.parseZbar, + 602 : self.parseEcond, + 603 : self.parseTcond, + 604 : self.parseTcond, + 605 : self.parseTcond, + } self.data = {} @@ -110,41 +102,36 @@ def __init__(self, filename, precision, verbose=False): self.parse() + def parse(self): while True: header = self.fhand.readline() - print(header.split()) - if not header: - break # Reached EOF - if header[:3] == " 2 ": - break + print (header.split()) + if not header: break # Reached EOF + if header[:3] == " 2 ": break matid = int(header.split()[1]) recid = int(header.split()[2]) nentries = int(header.split()[3]) - if not matid in self.data: - self.data[matid] = {} + if not matid in self.data: self.data[matid] = {} if self.verbose and (recid > 104): - print( - "Material = %8i Record = %8i Entries = %8i" - % (matid, recid, nentries) - ) + print("Material = %8i Record = %8i Entries = %8i" % (matid, recid, nentries)) # If there are opacity tables (series 500), # we do not currently read them. - if (600 > recid) and (recid >= 500): - print( - "Warning: No reader for SESAME opacity tables. The " - "data in table number {} will not be parsed".format(recid) - ) + if (600>recid) and (recid>=500): + print("Warning: No reader for SESAME opacity tables. The " + "data in table number {} will not be parsed".format(recid)) + if not recid in self.fdict: raise ValueError("No handling function for record %d" % recid) - self.fdict[recid](nentries, matid, recid) + + self.fdict[recid](nentries,matid, recid) if matid not in self.recs.keys(): self.recs[matid] = [recid] @@ -153,7 +140,7 @@ def parse(self): def parseComment(self, nentries, matid, recid): - nlines = (nentries - 1) // self.CHAR_LINE_LEN + 1 + nlines = (nentries-1) // self.CHAR_LINE_LEN + 1 nchar = nentries + nlines return self.fhand.read(nchar) @@ -165,15 +152,18 @@ def parseInfo(self, nentries, matid, recid): self.data[matid]["bulkmod"] = words[3] self.data[matid]["excoef"] = words[4] - if self.verbose: - print(" zbar = %g abar = %g rho0 = %g" % (words[0], words[1], words[2])) + if self.verbose: print(" zbar = %g abar = %g rho0 = %g" % (words[0],words[1],words[2])) return self.data def parseEos(self, nentries, matid, recid): words = self.readEntries(nentries) - prefixes = {301: "total_", 303: "ioncc_", 304: "ele_", 305: "ion_", 306: "cc_"} + prefixes = { 301 : "total_", + 303 : "ioncc_", + 304 : "ele_", + 305 : "ion_", + 306 : "cc_" } prefix = prefixes[recid] @@ -181,13 +171,13 @@ def parseEos(self, nentries, matid, recid): ntemp = int(words[1]) start = 2 - self.data[matid][prefix + "ndens"] = ndens - self.data[matid][prefix + "ntemp"] = ntemp + self.data[matid][prefix+"ndens"] = ndens + self.data[matid][prefix+"ntemp"] = ntemp - self.data[matid][prefix + "dens"] = words[start : start + ndens] # In g/cc - start = start + ndens + self.data[matid][prefix+"dens"] = words[start:start+ndens] # In g/cc + start = start+ndens - self.data[matid][prefix + "temps"] = words[start : start + ntemp] * KELVIN_TO_EV + self.data[matid][prefix+"temps"] = words[start:start+ntemp]*KELVIN_TO_EV start = start + ntemp if self.verbose and prefix == "total_": @@ -198,19 +188,14 @@ def parseEos(self, nentries, matid, recid): print("dens[0] = %13.6e dens[-1] = %13.6e" % (dens[0], dens[-1])) print("temp[0] = %13.6e temp[-1] = %13.6e" % (temps[0], temps[-1])) + # Read pressure array (in GPa and convert to ergs/cc): - self.data[matid][prefix + "pres"] = ( - words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() - * GPA_TO_ERGCC - ) - start = start + ntemp * ndens + self.data[matid][prefix+"pres"] = words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()*GPA_TO_ERGCC + start = start + ntemp*ndens # Read specific internal energy array (in MJ/kg and convert to ergs/g): - self.data[matid][prefix + "eint"] = ( - words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() - * MJKG_TO_ERGCC - ) - start = start + ntemp * ndens + self.data[matid][prefix+"eint"] = words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()*MJKG_TO_ERGCC + start = start + ntemp*ndens if start == nentries: # There is no Helmholtz free energy data, this is the end @@ -218,11 +203,8 @@ def parseEos(self, nentries, matid, recid): return # Read specific Helmholtz free energy array (in MJ/kg and convert to ergs/g): - self.data[matid][prefix + "free"] = ( - words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() - * MJKG_TO_ERGCC - ) - start = start + ntemp * ndens + self.data[matid][prefix+"free"] = words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()*MJKG_TO_ERGCC + start = start + ntemp*ndens # exit() @@ -247,18 +229,16 @@ def parseZbar(self, nentries, matid, recid): ntemp = int(words[1]) start = 2 - self.data[matid][prefix + "_ndens"] = ndens - self.data[matid][prefix + "_ntemp"] = ntemp + self.data[matid][prefix+"_ndens"] = ndens + self.data[matid][prefix+"_ntemp"] = ntemp - self.data[matid][prefix + "_dens"] = 10 ** words[start : start + ndens] - start = start + ndens + self.data[matid][prefix+"_dens"] = 10**words[start:start+ndens] + start = start+ndens - self.data[matid][prefix + "_temps"] = 10 ** words[start : start + ntemp] + self.data[matid][prefix+"_temps"] = 10**words[start:start+ntemp] start = start + ntemp - self.data[matid][prefix] = 10 ** ( - words[start : start + ntemp * ndens].reshape((ntemp, ndens)).transpose() - ) + self.data[matid][prefix] = 10**(words[start:start+ntemp*ndens].reshape((ntemp,ndens)).transpose()) def parseEcond(self, nentries, matid, recid): words = self.readEntries(nentries) @@ -266,42 +246,32 @@ def parseEcond(self, nentries, matid, recid): def parseTcond(self, nentries, matid, recid): words = self.readEntries(nentries) - def readEntries(self, nentries): - nlines = (nentries - 1) // self.WORDS_PER_LINE + 1 + def readEntries(self,nentries): + nlines = (nentries-1) // self.WORDS_PER_LINE + 1 data = np.empty(nentries) string = "" for i in range(nlines): - string += self.fhand.readline()[: self.WORDS_PER_LINE * self.entry_len] + string += self.fhand.readline()[:self.WORDS_PER_LINE*self.entry_len] for i in range(nentries): - word = string[self.entry_len * i : self.entry_len * (i + 1)] - if word[-4] == "-": - word = word[:-4] + "E" + word[-4:] + word = string[self.entry_len*i:self.entry_len*(i+1)] + if word[-4] == '-': word = word[:-4] + "E" + word[-4:] data[i] = float(word) return data - def toEosDict( - self, - Znum=None, - Anum=None, - Xnum=None, - qeos=False, - log=None, - filter_dens=0.0, - filter_temps=0.0, - tabnum=None, - ): + def toEosDict(self, Znum=None, Anum=None, + Xnum=None, qeos=False, log=None, + filter_dens=0., filter_temps=0., + tabnum=None): # For SESAME, we need the hedp package to calculate zbar. try: from hedp import eos except ImportError: - raise ImportError( - "You need the hedp module. You can get it here: " - "https://github.com/luli/hedp." - ) + raise ImportError('You need the hedp module. You can get it here: ' + 'https://github.com/luli/hedp.') if tabnum is None: # Select the last table (newest) table available. @@ -310,7 +280,7 @@ def toEosDict( try: opp_ses_data = self.data[tabnum] except KeyError: - raise KeyError("Invalid table number!") + raise KeyError('Invalid table number!') # Sesame has extra data points in it, so we must merge them down. We are # merging the default grids ioncc_ and ele_. @@ -319,110 +289,102 @@ def toEosDict( # We must filter dens > 0. in order to avoid problems with calculating # zbar below, since eos.thomas_fermi_ionization() returns `nan` where # density is 0. - # if not qeos: - # opp_ses_data = opp.utils.EosMergeGrids( - # opp_ses_data, - # filter_dens=lambda x: (x>filter_dens), - # filter_temps=lambda x: (x>filter_temps)) + if not qeos: + opp_ses_data = opp.utils.EosMergeGrids( + opp_ses_data, + filter_dens=lambda x: (x>filter_dens), + filter_temps=lambda x: (x>filter_temps)) # If we are dealing with SESAME generated by qeos, we need to merge # ion_ and ele_ grids rather than ioncc_ and ele_. We use the same # filters as above. if qeos: opp_ses_data = opp.utils.EosMergeGrids( - opp_ses_data, - intersect=["ele", "ion"], - filter_dens=lambda x: (x > filter_dens), - filter_temps=lambda x: (x > filter_temps), - qeos=True, - ) + opp_ses_data, intersect=['ele', 'ion'], + filter_dens=lambda x: (x>filter_dens), + filter_temps=lambda x: (x>filter_temps), + qeos=True) # Converting density to ion number density. - opp_ses_data["idens"] = (opp.NA * opp_ses_data["ele_dens"]) / opp_ses_data[ - "abar" - ] + opp_ses_data['idens'] = ((opp.NA * opp_ses_data['ele_dens']) + / opp_ses_data['abar']) + # Adjust for Znum, Xnum, and Anum. if Znum is None: - if "Znum" in opp_ses_data: - Znum = opp_ses_data["Znum"] + if 'Znum' in opp_ses_data: + Znum = opp_ses_data['Znum'] else: - raise ValueError("Znum Varray should be provided!") + raise ValueError('Znum Varray should be provided!') if type(Znum) is int: Znum = [Znum] - opp_ses_data["Znum"] = np.array(Znum, dtype="int") + opp_ses_data['Znum'] = np.array(Znum, dtype='int') if Anum is None: - opp_ses_data["Anum"] = np.array( - [ptab.elements[el].mass for el in opp_ses_data["Znum"]] - ) + opp_ses_data['Anum'] = np.array([ptab.elements[el].mass for el in opp_ses_data['Znum']]) else: - opp_ses_data["Anum"] = np.array(Anum) - opp_ses_data["Zsymb"] = np.array( - [ptab.elements[el].symbol for el in opp_ses_data["Znum"]], dtype="|S2" - ) + opp_ses_data['Anum'] = np.array(Anum) + opp_ses_data['Zsymb'] = np.array([ptab.elements[el].symbol for el in opp_ses_data['Znum']], dtype='|S2') if Xnum is None: if len(Znum) == 1: - opp_ses_data["Xnum"] = np.array([1.0]) + opp_ses_data['Xnum'] = np.array([1.0]) else: - raise ValueError("Xnum array should be provided!") + raise ValueError('Xnum array should be provided!') else: - opp_ses_data["Xnum"] = np.array(Xnum) + opp_ses_data['Xnum'] = np.array(Xnum) + - # turn this off later. just for testing sesame writer - # ionization is possibly on a completely different grid - if not "zbar" in opp_ses_data.keys(): + #ionization is possibly on a completely different grid + if (not 'zbar' in opp_ses_data.keys()): # Calculate zbar using thomas_fermi_ionization. # If there are multiple elements, it suffices to use the average # atomic number in this calculation - JTL - dens_arr, temp_arr = np.meshgrid( - opp_ses_data["ele_dens"], opp_ses_data["ele_temps"] - ) - zbar = eos.thomas_fermi_ionization( - dens_arr, temp_arr, opp_ses_data["Znum"].mean(), opp_ses_data["abar"] - ).T - opp_ses_data["zbar"] = zbar + dens_arr, temp_arr = np.meshgrid(opp_ses_data['ele_dens'], + opp_ses_data['ele_temps']) + zbar = eos.thomas_fermi_ionization(dens_arr, + temp_arr, + opp_ses_data['Znum'].mean(), + opp_ses_data['abar']).T + opp_ses_data['zbar'] = zbar # Translating SESAME names to common dictionary format. if qeos: # Names are slightly different for QEOS SESAME - names_dict = { - "idens": "idens", - "ele_temps": "temp", # We merged ele_ and ion_ dens & - # temp grids for qeos. - "ele_dens": "dens", - "zbar": "Zf_DT", - "total_eint": "Ut_DT", # But not their energies. - "ele_eint": "Uec_DT", - "ion_eint": "Ui_DT", - "ion_pres": "Pi_DT", - "ele_pres": "Pec_DT", - "Znum": "Znum", - "Xnum": "Xnum", - "bulkmod": "BulkMod", - "abar": "Abar", - "zmax": "Zmax", - } + names_dict = {'idens' :'idens' , + 'ele_temps' :'temp' , # We merged ele_ and ion_ dens & + # temp grids for qeos. + 'ele_dens' :'dens' , + 'zbar' :'Zf_DT' , + 'total_eint':'Ut_DT' , # But not their energies. + 'ele_eint' :'Uec_DT' , + 'ion_eint' :'Ui_DT' , + 'ion_pres' :'Pi_DT' , + 'ele_pres' :'Pec_DT' , + 'Znum' :'Znum' , + 'Xnum' :'Xnum' , + 'bulkmod' :'BulkMod', + 'abar' :'Abar' , + 'zmax' :'Zmax' + } else: - names_dict = { - "idens": "idens", - "ele_temps": "temp", # We merged ele_ and ioncc_ dens & - # temp grids. - "ele_dens": "dens", - "zbar": "Zf_DT", - "total_eint": "Ut_DT", # But not their energies. - "ele_eint": "Uec_DT", - "ioncc_eint": "Ui_DT", - "ioncc_pres": "Pi_DT", - "ele_pres": "Pec_DT", - "Znum": "Znum", - "Xnum": "Xnum", - "bulkmod": "BulkMod", - "abar": "Abar", - "zmax": "Zmax", - } + names_dict = {'idens':'idens', + 'ele_temps':'temp', # We merged ele_ and ioncc_ dens & + # temp grids. + 'ele_dens':'dens', + 'zbar':'Zf_DT', + 'total_eint':'Ut_DT', # But not their energies. + 'ele_eint':'Uec_DT', + 'ioncc_eint':'Ui_DT', + 'ioncc_pres':'Pi_DT', + 'ele_pres':'Pec_DT', + 'Znum':'Znum', + 'Xnum':'Xnum', + 'bulkmod':'BulkMod', + 'abar':'Abar', + 'zmax':'Zmax' + } # Initialize dictionary. eos_dict = {} @@ -432,7 +394,7 @@ def toEosDict( try: eos_dict[eos_key] = opp_ses_data[ses_key] except KeyError: - print("No data for {} is being written.".format(eos_key)) + print('No data for {} is being written.'.format(eos_key)) # Handle the logarithmic data. if log is not None: @@ -440,43 +402,35 @@ def toEosDict( if key in log: eos_dict[key] = np.log10(eos_dict[key]) return eos_dict - - def writeSesameFile(fn, t201, t301, t303, t304, t305, t502, t504, t505, t601): CHAR_LINE_LEN = 80 WORDS_PER_LINE = 5 - comment = "This Sesame table was generated using OpacPlot2 to convert an IONMIX EoS table to the Sesame format" - f = open(fn, "w") - # write comments - comment = comment + " " * (80 - len(comment) % 80) - header = " 0 9999 101 %d r 82803 22704 1" % len(comment) - header = header + (79 - len(header)) * " " + "0\n" + comment = 'This Sesame table was generated using OpacPlot2 to convert an IONMIX EoS table to the Sesame format' + f = open(fn, 'w') + #write comments + comment = comment + ' '*(80-len(comment)%80) + header = ' 0 9999 101 %d r 82803 22704 1' % len(comment) + header = header + (79-len(header))*' ' + '0\n' f.write(header) - for i in range(len(comment) // 80): - f.write(comment[i * 80 : (i + 1) * 80] + "\n") + for i in range(len(comment)//80): + f.write(comment[i*80:(i+1)*80]+'\n') def theader(num): - return ( - " 1 9999 %d 240 r 82803 22704 1 1\n" - % num - ) - + return(' 1 9999 %d 240 r 82803 22704 1 1\n' % num) def write_block(fid, num, data): - header = " 1 9999 %d %d r 82803 22704 1" % (num, len(data)) - header = header + (79 - len(header)) * " " + "1\n" + header = ' 1 9999 %d %d r 82803 22704 1' % (num,len(data)) + header = header + (79-len(header))*' ' + '1\n' fid.write(header) count = 0 for n in range(len(data)): count += 1 - # fid.write('%22.15E' % data[n]) - fid.write("%15.8E" % data[n]) + fid.write('%22.15E' % data[n]) if count == 5: count = 0 - fid.write("11111\n") + fid.write('11111\n') if count > 0: - m = 5 - count - fid.write(m * 22 * " " + count * "1" + m * "0" + "\n") - # fid.write(m*16*' ' + count*'1' + m*'0' + '\n') + m = 5-count + fid.write(m*22*' ' + count*'1' + m*'0' + '\n') write_block(f, 201, t201) write_block(f, 301, t301) diff --git a/opacplot2/scripts/opac_convert.py b/opacplot2/scripts/opac_convert.py index 460bebf..27697f8 100644 --- a/opacplot2/scripts/opac_convert.py +++ b/opacplot2/scripts/opac_convert.py @@ -8,70 +8,56 @@ def get_input_data(): # Available formats. avail_output_formats = ["ionmix", "sesame"] - avail_input_formats = ["propaceos", "multi", "sesame", "sesame-qeos", "ionmix"] + avail_input_formats = ["propaceos", "multi", "sesame", "sesame-qeos", "ionmix", "tops"] # Creating the argument parser. parser = argparse.ArgumentParser( - description="This script is used to browse various" - "EoS/Opacity tables formats." - ) - - parser.add_argument( - "-v", "--verbose", action="store_const", const=True, help="Verbosity option." - ) - - parser.add_argument( - "--Znum", - action="store", - type=str, - help="Comma separated list of Z" "numbers for every component.", - ) + description="This script is used to browse various" + "EoS/Opacity tables formats.") + + parser.add_argument('-v', '--verbose', + action='store_const', const=True, + help='Verbosity option.') + + parser.add_argument('--Znum', + action='store', type=str, + help='Comma separated list of Z' + 'numbers for every component.') + + parser.add_argument('--Xfracs', + action='store', type=str, + help='Comma separated list of X' + 'fractions for every component.') + + parser.add_argument('--outname', + action='store', type=str, + help='Name for output file without extension.') + + parser.add_argument('-i', '--input', + action='store', type=str, + choices=avail_input_formats, + help='Input filetype.') + + parser.add_argument('-o', '--output', + action='store', type=str, + choices=avail_output_formats, + default='ionmix', + help='Output filetype. Default: IONMIX.') + + parser.add_argument('input_file', + action='store', type=str, + help='Input file.') + + parser.add_argument('--log', + action='store', type=str, + help='Logarithmic data keys.') + + parser.add_argument('--tabnum', + action='store', type=str, + help='Specify the SESAME table number.') parser.add_argument("--mpi", action="store", type=str, help="Mass per ion in grams") - parser.add_argument( - "--Xfracs", - action="store", - type=str, - help="Comma separated list of X" "fractions for every component.", - ) - - parser.add_argument( - "--outname", - action="store", - type=str, - help="Name for output file without extension.", - ) - - parser.add_argument( - "-i", - "--input", - action="store", - type=str, - choices=avail_input_formats, - help="Input filetype.", - ) - - parser.add_argument( - "-o", - "--output", - action="store", - type=str, - choices=avail_output_formats, - default="ionmix", - help="Output filetype. Default: IONMIX.", - ) - - parser.add_argument("input_file", action="store", type=str, help="Input file.") - - parser.add_argument( - "--log", action="store", type=str, help="Logarithmic data keys." - ) - - parser.add_argument( - "--tabnum", action="store", type=str, help="Specify the SESAME table number." - ) - args = parser.parse_args() # Get the relevant paths and filenames. @@ -90,48 +76,44 @@ def get_input_data(): # Create lists out of the strings for Znum, Xfracs, mpi, and log if given. if args.Znum is not None: - args.Znum = [int(num) for num in args.Znum.split(",")] + args.Znum = [int(num) for num in args.Znum.split(',')] if args.Xfracs is not None: - args.Xfracs = [float(num) for num in args.Xfracs.split(",")] + args.Xfracs = [float(num) for num in args.Xfracs.split(',')] if args.mpi is not None: args.mpi = [float(num) for num in args.mpi.split(",")] if args.log is not None: - args.log = [str(key) for key in args.log.split(",")] + args.log = [str(key) for key in args.log.split(',')] # Convert tabnum into int. if args.tabnum is not None: try: args.tabnum = int(args.tabnum) except ValueError: - raise ValueError("Please provide a valid SESAME table number.") + raise ValueError('Please provide a valid SESAME table number.') - input_data = { - "args": args, - "basename": basename, - "path_in": path_in, - "basedir": basedir, - "fn_in": fn_in, - } + input_data = {'args' : args, + 'basename' : basename, + 'path_in' : path_in, + 'basedir' : basedir, + 'fn_in' : fn_in} return input_data - def read_format_ext(args, fn_in): # Try to read from the input file extension. - ext_dict = { - ".prp": "propaceos", - ".eps": "multi", - ".opp": "multi", - ".opz": "multi", - ".opr": "multi", - ".cn4": "ionmix", - ".mexport": "sesame-qeos", - ".ses": "sesame", - ".html": "tops", - ".tops": "tops", - } + ext_dict = {'.prp':'propaceos', + '.eps':'multi', + '.opp':'multi', + '.opz':'multi', + '.opr':'multi', + ".cn4": "ionmix", + '.mexport':'sesame-qeos', + '.ses':'sesame', + '.html':'tops', + '.tops':'tops', + } # If the input file is compressed, choose the next extension. - if os.path.splitext(fn_in)[1] == ".gz": + if os.path.splitext(fn_in)[1] == '.gz': _, ext = os.path.splitext(os.path.splitext(fn_in)[0]) else: _, ext = os.path.splitext(fn_in) @@ -141,18 +123,14 @@ def read_format_ext(args, fn_in): if ext in ext_dict.keys(): args.input = ext_dict[ext] else: - raise Warning( - "Cannot tell filetype from extension. Please specify " - "input file type with --input." - ) - + raise Warning('Cannot tell filetype from extension. Please specify ' + 'input file type with --input.') class Formats_toEosDict(object): """ Contains handling functions to convert different types of tables into a common EoS dictionary for IONMIX. """ - def __init__(self, args, basedir, basename, path_in): # Initialize the dictionary for handling functions. self.set_handle_dict() @@ -167,35 +145,33 @@ def __init__(self, args, basedir, basename, path_in): try: self.eos_dict = self.handle_dict[args.input]() except KeyError: - raise KeyError("Must use valid format name.") + raise KeyError('Must use valid format name.') def set_handle_dict(self): - self.handle_dict = { - "propaceos": self.propaceos_toEosDict, - "multi": self.multi_toEosDict, - "sesame": self.sesame_toEosDict, - "ionmix": self.ionmix_toEosDict, - "tops": self.tops_toEosDict, - "sesame-qeos": self.sesame_qeos_toEosDict, - } + self.handle_dict = {'propaceos' : self.propaceos_toEosDict, + 'multi' : self.multi_toEosDict, + 'sesame' : self.sesame_toEosDict, + "ionmix": self.ionmix_toEosDict, + 'sesame-qeos' : self.sesame_qeos_toEosDict, + 'tops' : self.tops_toEosDict, + } def propaceos_toEosDict(self): # If we are unable to find the correct library for opg_propaceos # we need to let the user know. try: import opacplot2.opg_propaceos - op = opp.opg_propaceos.OpgPropaceosAscii(self.path_in, mpi=self.args.mpi) eos_dict = op.toEosDict(log=self.args.log) return eos_dict except ImportError: - raise ImportError("You do not have the opg_propaceos script.") + raise ImportError('You do not have the opg_propaceos script.') def multi_toEosDict(self): op = opp.OpgMulti.open_file(self.basedir, self.basename) - eos_dict = op.toEosDict( - Znum=self.args.Znum, Xnum=self.args.Xfracs, log=self.args.log - ) + eos_dict = op.toEosDict(Znum=self.args.Znum, + Xnum=self.args.Xfracs, + log=self.args.log) return eos_dict def ionmix_toEosDict(self): @@ -212,52 +188,43 @@ def sesame_toEosDict(self): except ValueError: op = opp.OpgSesame(self.path_in, opp.OpgSesame.DOUBLE) + if len(op.data.keys()) > 1: - raise Warning( - "More than one material ID found. " - "Use sesame-extract to create a file " - "with only one material first." - ) + raise Warning('More than one material ID found. ' + 'Use sesame-extract to create a file ' + 'with only one material first.') if self.args.tabnum is not None: - eos_dict = op.toEosDict( - Znum=self.args.Znum, - Xnum=self.args.Xfracs, - log=self.args.log, - tabnum=self.args.tabnum, - ) + eos_dict = op.toEosDict(Znum=self.args.Znum, + Xnum=self.args.Xfracs, + log=self.args.log, + tabnum=self.args.tabnum) else: - eos_dict = op.toEosDict( - Znum=self.args.Znum, Xnum=self.args.Xfracs, log=self.args.log - ) + eos_dict = op.toEosDict(Znum=self.args.Znum, + Xnum=self.args.Xfracs, + log=self.args.log) return eos_dict def sesame_qeos_toEosDict(self): - raise Warning("QEOS-SESAME is not ready yet!") + raise Warning('QEOS-SESAME is not ready yet!') try: op = opp.OpgSesame(self.path_in, opp.OpgSesame.SINGLE) except ValueError: op = opp.OpgSesame(self.path_in, opp.OpgSesame.DOUBLE) + if len(op.data.keys()) > 1: - raise Warning( - "More than one material ID found. " - "Use sesame-extract to create a file " - "with only one material first." - ) + raise Warning('More than one material ID found. ' + 'Use sesame-extract to create a file ' + 'with only one material first.') if self.args.tabnum is not None: - eos_dict = op.toEosDict( - Znum=self.args.Znum, - Xnum=self.args.Xfracs, - qeos=True, - log=self.args.log, - tabnum=self.args.tabnum, - ) + eos_dict = op.toEosDict(Znum=self.args.Znum, Xnum=self.args.Xfracs, + qeos=True, log=self.args.log, + tabnum=self.args.tabnum) else: - eos_dict = op.toEosDict( - Znum=self.args.Znum, Xnum=self.args.Xfracs, qeos=True, log=self.args.log - ) + eos_dict = op.toEosDict(Znum=self.args.Znum, Xnum=self.args.Xfracs, + qeos=True, log=self.args.log) return eos_dict def tops_toEosDict(self): @@ -344,7 +311,6 @@ class EosDict_toIonmixFile(object): """ Takes a common EoS dictionary and writes it to the correct output format. """ - def __init__(self, args, eos_dict): # Initialize the handling function dictionary. self.set_handle_dict() @@ -357,25 +323,23 @@ def __init__(self, args, eos_dict): self.handle_dict[args.output]() def set_handle_dict(self): - self.handle_dict = {"ionmix": self.eosDict_toIonmix} + self.handle_dict = {'ionmix' : self.eosDict_toIonmix} def eosDict_toIonmix(self): # These are the naming conventions translated to ionmix arguments. - imx_conv = { - "Znum": "zvals", - "Xnum": "fracs", - "idens": "numDens", - "temp": "temps", - "Zf_DT": "zbar", - "Pi_DT": "pion", - "Pec_DT": "pele", - "Ui_DT": "eion", - "Uec_DT": "eele", - "groups": "opac_bounds", - "opr_mg": "rosseland", - "opp_mg": "planck_absorb", - "emp_mg": "planck_emiss", - } + imx_conv = {'Znum':'zvals', + 'Xnum':'fracs', + 'idens':'numDens', + 'temp':'temps', + 'Zf_DT':'zbar', + 'Pi_DT':'pion', + 'Pec_DT':'pele', + 'Ui_DT':'eion', + 'Uec_DT':'eele', + 'groups':'opac_bounds', + 'opr_mg':'rosseland', + 'opp_mg':'planck_absorb', + 'emp_mg':'planck_emiss'} # Initialize ionmix argument dictionary. imx_dict = {} @@ -386,66 +350,61 @@ def eosDict_toIonmix(self): imx_dict[imx_conv[key]] = self.eos_dict[key] # Set ngroups if opacity bounds are present. - if "opac_bounds" in imx_dict: - imx_dict["ngroups"] = len(self.eos_dict["groups"]) - 1 + if 'opac_bounds' in imx_dict: + imx_dict['ngroups'] = len(self.eos_dict['groups']) - 1 # Check if required FLASH EoS tables are present. - imx_req_keys = ["zbar", "eion", "eele", "pion", "pele"] + imx_req_keys = ['zbar', 'eion', 'eele', 'pion', 'pele'] if not all(key in imx_dict for key in imx_req_keys): - print( - "The required EoS data for FLASH is not present...\n" - "Aborting the IONMIX file creation..." - ) + print("The required EoS data for FLASH is not present...\n" + "Aborting the IONMIX file creation...") raise Warning("Missing EoS data for IONMIX file to run in FLASH") # For verbose flag. if self.args.verbose: - verb_conv = { - "zvals": "Atomic numbers", - "fracs": "Element fractions", - "numDens": "Ion number densities", - "temps": "Temperatures", - "zbar": "Average ionizations", - "pion": "Ion pressure", - "pele": "Electron pressure", - "eion": "Ion internal energy", - "eele": "Electron internal energy", - "opac_bounds": "Opacity bounds", - "rosseland": "Rosseland mean opacity", - "planck_absorb": "Absorption Planck mean opacity", - "planck_emiss": "Emission Planck mean opacity", - "ngroups": "Number of opacity groups", - } - - verb_str = "Wrote the following data to IONMIX file:\n" + verb_conv = {'zvals':'Atomic numbers', + 'fracs':'Element fractions', + 'numDens':'Ion number densities', + 'temps':'Temperatures', + 'zbar':'Average ionizations', + 'pion':'Ion pressure', + 'pele':'Electron pressure', + 'eion':'Ion internal energy', + 'eele':'Electron internal energy', + 'opac_bounds':'Opacity bounds', + 'rosseland':'Rosseland mean opacity', + 'planck_absorb':'Absorption Planck mean opacity', + 'planck_emiss':'Emission Planck mean opacity', + 'ngroups':'Number of opacity groups'} + + verb_str = 'Wrote the following data to IONMIX file:\n' i = 0 for key in imx_dict.keys(): i = i + 1 if i == len(imx_dict.keys()): - verb_str = verb_str + "{}. {}".format(i, verb_conv[key]) + verb_str = verb_str + '{}. {}'.format(i, verb_conv[key]) else: - verb_str = verb_str + "{}. {} \n".format(i, verb_conv[key]) + verb_str = verb_str + '{}. {} \n'.format(i, verb_conv[key]) print(verb_str) # Write the ionmix file based on what data is stored in imx_dict. - opp.writeIonmixFile(self.args.outname + ".cn4", **imx_dict) - + opp.writeIonmixFile(self.args.outname + '.cn4', **imx_dict) def convert_tables(): # Grab the input data. input_data = get_input_data() # Read the file extension if the user did not specify an input. - if input_data["args"].input is None: - read_format_ext(input_data["args"], input_data["fn_in"]) + if input_data['args'].input is None: + read_format_ext(input_data['args'], input_data['fn_in']) + + # Reading in data and converting it to the common dictionary format. - eos_dict = Formats_toEosDict( - input_data["args"], - input_data["basedir"], - input_data["basename"], - input_data["path_in"], - ).eos_dict + eos_dict = Formats_toEosDict(input_data['args'], + input_data['basedir'], + input_data['basename'], + input_data['path_in']).eos_dict output_type = input_data["args"].output if output_type == "sesame": From 1695eaac7cb591e3c9d4872261ae67506e21f9a2 Mon Sep 17 00:00:00 2001 From: adam reyes Date: Mon, 2 Dec 2024 13:33:36 -0800 Subject: [PATCH 4/4] missed a debug print --- opacplot2/opg_sesame.py | 1 - 1 file changed, 1 deletion(-) diff --git a/opacplot2/opg_sesame.py b/opacplot2/opg_sesame.py index a5c821c..7d25b3a 100644 --- a/opacplot2/opg_sesame.py +++ b/opacplot2/opg_sesame.py @@ -107,7 +107,6 @@ def parse(self): while True: header = self.fhand.readline() - print (header.split()) if not header: break # Reached EOF if header[:3] == " 2 ": break