Add sesame tables as an output

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
 

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,

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()

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()