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..ee32ca7 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,39 @@ 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. + 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, diff --git a/opacplot2/opg_sesame.py b/opacplot2/opg_sesame.py index 58ae982..7d25b3a 100644 --- a/opacplot2/opg_sesame.py +++ b/opacplot2/opg_sesame.py @@ -5,12 +5,13 @@ 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 @@ -87,6 +88,7 @@ def __init__(self, filename, precision, verbose=False): 412 : self.parseLiquid, 431 : self.parseShear, 432 : self.parseShear, + 504 : self.parseZbar, 601 : self.parseZbar, 602 : self.parseEcond, 603 : self.parseTcond, @@ -332,35 +334,37 @@ def toEosDict(self, Znum=None, Anum=None, 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 + #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: @@ -397,3 +401,43 @@ 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]) + if count == 5: + count = 0 + fid.write('11111\n') + if count > 0: + m = 5-count + fid.write(m*22*' ' + 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..27697f8 100644 --- a/opacplot2/scripts/opac_convert.py +++ b/opacplot2/scripts/opac_convert.py @@ -1,11 +1,14 @@ 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", "tops"] # Creating the argument parser. parser = argparse.ArgumentParser( @@ -53,6 +56,8 @@ def get_input_data(): action='store', type=str, help='Specify the SESAME table number.') + parser.add_argument("--mpi", action="store", type=str, help="Mass per ion in grams") + args = parser.parse_args() # Get the relevant paths and filenames. @@ -69,11 +74,13 @@ 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(',')] if args.Xfracs is not None: 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(',')] @@ -99,6 +106,7 @@ def read_format_ext(args, fn_in): '.opp':'multi', '.opz':'multi', '.opr':'multi', + ".cn4": "ionmix", '.mexport':'sesame-qeos', '.ses':'sesame', '.html':'tops', @@ -143,6 +151,7 @@ def set_handle_dict(self): 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, } @@ -152,7 +161,7 @@ def propaceos_toEosDict(self): # 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: @@ -165,6 +174,14 @@ def multi_toEosDict(self): log=self.args.log) return eos_dict + def ionmix_toEosDict(self): + 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): try: op = opp.OpgSesame(self.path_in, opp.OpgSesame.SINGLE) @@ -215,6 +232,81 @@ 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"])): + znum += zz[i] * xx[i] + 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. @@ -314,7 +406,12 @@ def convert_tables(): input_data['basename'], input_data['path_in']).eos_dict - EosDict_toIonmixFile(input_data['args'], eos_dict) + output_type = input_data["args"].output + if output_type == "sesame": + EosDict_toSesameFile(input_data["args"], eos_dict) + else: + EosDict_toIonmixFile(input_data["args"], eos_dict) + -if __name__=='__main__': +if __name__ == "__main__": convert_tables()