Skip to content

Use FAC data files for CR modeling #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
cjperks7 wants to merge 46 commits into
base: master
Choose a base branch
from
Open

Use FAC data files for CR modeling #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
46 commits
Select commit Hold shift + click to select a range
d7dba63
Init function to read FAC data
Dec 18, 2023
bcca908
Started formating for FAC file management
Dec 18, 2023
daf1a91
Init file loading management
Dec 18, 2023
826e4e5
Init function to read energy level FAC files
Dec 18, 2023
457100c
Init function to read Einstein coefficient FAC files
Dec 18, 2023
101ee79
Init function to read collisional excitation FAC files
Dec 18, 2023
6f60831
Init function to read radiative recombination FAC files
Dec 18, 2023
361b68e
Init function to read collisional ionization FAC files
Dec 18, 2023
66a5800
Init function to format FAC output for ColRadPy
Dec 18, 2023
283cad6
Init using FAC reading functions in populating CRM
Dec 18, 2023
044fa63
[#18] debug reading FAC data
Jan 8, 2024
090aa55
[#18] debugged level indexing
Jan 9, 2024
4503dc9
[#18] debugged expected input file Te units
Jan 9, 2024
c1f0b90
[#18] debugged having to convert excit rates into adf04 upsilon form
Jan 18, 2024
961aed8
[#18] init function to convolve cross-sections with EEDF
Jan 18, 2024
91fb794
[#18] function to calc Maxwellian dist
Jan 18, 2024
a27895b
[#18] function to preform rate coefficient integration
Jan 18, 2024
ec67969
[#18] init function to convolve EEDF for Dielectronic capture strengths
Jan 18, 2024
6f2c6be
[#18] added 2-photon decay to A-values for H-like, He-like systems
Jan 18, 2024
799fa6f
[#18] debugged convolve EEDF
Jan 19, 2024
bfcdfbc
[#18] corrected form of ionization rates to adf04 reduced form
Jan 19, 2024
d3d4a17
[#18] fixed indexing bug when loading ce data causing it to skip some…
Jan 24, 2024
a77898d
[#18] gneralized redundant code in handling mr file data
Jan 24, 2024
c2a6d0d
[#18] to shorten the script and helpful readability, I broke off a bu…
Jan 24, 2024
0d22557
[#18] made functions to use FAC excit cross-sections for convolution
Jan 24, 2024
c9248a5
[#18] added high-energy fill to the cross-section data by Bethe asymp…
Jan 24, 2024
30e31cb
[#18] added functions to convolve rad recomb and coll ioniz cross-sec…
Jan 24, 2024
2a38176
[#18] init high-energy asymptotic cross-section for rad recomb and co…
Jan 24, 2024
df62574
[#18] reduced redundancy in convolve_EEDF. Simplified high-energy cro…
Jan 25, 2024
e11cbd1
[#18] fill missing cross-section data right at threshold
Jan 25, 2024
7d04e86
[#18] reduced redundancy in loading FAC data
Jan 25, 2024
a29735b
[#18] added high-energy formula for rad recomb cross-section
Jan 25, 2024
2edfbfd
[#18] added dielectronic capture rate coefficients. NOTE: Ignored aut…
Jan 25, 2024
a7ecf7f
[#18] debugged having multiple ionized states for rad recomb cross-se…
Jan 25, 2024
86f97ae
[#18] misc debug
Feb 20, 2024
d1a289a
[#18] debug converting excitation rates to upsilon form
Feb 29, 2024
980eb31
[#18] added autoionization and various debug
Mar 7, 2024
be9b5ba
[#18] debug random nans
Mar 7, 2024
13d2993
[#18] added some convenient information in write_adf15
Apr 12, 2024
8796d43
[#18] debug for sim Kr
Apr 13, 2024
99bba99
[#18] fixed overflow error in ioniz
Apr 14, 2024
329ae16
[#18] init including excitation-autoionization to SCD rates
May 13, 2024
f24f4f4
[#18] minor debug
May 16, 2024
b27a5bc
[#18] more padding in adf15 transition table
May 16, 2024
f4cc854
[#18] added convolving Gaussian electron energy distribution function…
Jun 27, 2024
6e2f40e
[#18] various minor debug
Jun 27, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions colradpy/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,3 +7,5 @@
from .burgess_tully_rates import *
from .split_multiplet import *
from .nist_read_txt import *
from .read_FAC import *
from .convolve_EEDF import *
144 changes: 137 additions & 7 deletions colradpy/colradpy_class.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
sys.path.append('./')
from colradpy.r8necip import *
from colradpy.read_adf04_py3_class import *
from colradpy.read_FAC import *
from colradpy.ecip_rates import *
from colradpy.burgess_tully_rates import *
from colradpy.split_multiplet import *
Expand Down Expand Up @@ -160,7 +161,10 @@ def __init__(self,fil,metas=np.array([]),temp_grid=np.array([]),electron_den=np.
default_pop_norm=True,temp_dens_pair=False,rate_interp_ion = 'slinear',
rate_interp_recomb='log_slinear',rate_interp_col='log_slinear',
rate_interp_cx='log_quadratic',use_cx=False,scale_file_ioniz=False,
ne_tau = -1):
ne_tau = -1,
atomic_data_type='adf04', ele='H', nele=0, Zele=1,
EEDF='Maxwellian', atomic_physics='incl_all',
):

"""The initializing method. Sets up the nested list for data storage and starts to populate with user data
as well as reading in the adf04 file
Expand Down Expand Up @@ -238,6 +242,14 @@ def __init__(self,fil,metas=np.array([]),temp_grid=np.array([]),electron_den=np.
print('Exit here fix to run')
sys.exit()

# Atomic data file management
self.data['user']['atomic_data_type'] = atomic_data_type
self.data['user']['ele'] = ele
self.data['user']['nele'] = nele
self.data['user']['Zele'] = Zele
self.data['user']['EEDF'] = EEDF
self.data['user']['atomic_physics'] = atomic_physics

self.populate_data(fil)
self.data['atomic']['metas'] = np.asarray(metas)

Expand Down Expand Up @@ -325,8 +337,22 @@ def populate_data(self,fil):

"""
if(type(fil) == str or type(fil) == np.str_):

self.data = self.update_dict(self.data,read_adf04(fil))
# If atomic data is in adf04-format
if self.data['user']['atomic_data_type'] == 'adf04':
self.data = self.update_dict(self.data,read_adf04(fil))
# If atomic data is in FAC-format
elif self.data['user']['atomic_data_type'] == 'FAC':
self.data = self.update_dict(
self.data, read_FAC(
ele=self.data['user']['ele'],
nele=self.data['user']['nele'],
Zele=self.data['user']['Zele'],
fil=fil,
EEDF=self.data['user']['EEDF'],
reacts=self.data['user']['atomic_physics'],
Te=self.data['user']['temp_grid'], # [eV]
)
)
self.data['user']['file_loc'] = str(fil)#make this a string changed for hdfdict
else:
if(self.data):
Expand Down Expand Up @@ -860,6 +886,14 @@ def populate_cr_matrix(self):
np.sum(np.einsum('ij,j->ij',self.data['rates']['ioniz']['ionization'][i,:,:],
self.data['user']['dens_grid']),axis=0)

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['cr_matrix']['cr'][i,i,:] -= np.sum(
self.data['rates']['ioniz']['autoioniz'][i,:]
)
self.data['cr_matrix']['cr_loss'][i,i,:] -= np.sum(
self.data['rates']['ioniz']['autoioniz'][i,:]
)


self.data['cr_matrix']['cr'][i,0:len(self.data['atomic']['energy']),:] = \
Expand Down Expand Up @@ -918,6 +952,15 @@ def populate_cr_matrix(self):
0:len(self.data['atomic']['energy']),:] + \
np.einsum('ij,j->ij',self.data['rates']['ioniz']['ionization'][:,p,:],
self.data['user']['dens_grid'])

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['cr_matrix']['cr'][len(self.data['atomic']['energy'])+ p,
0:len(self.data['atomic']['energy']),:] += (
self.data['rates']['ioniz']['autoioniz'][:,p]
)[:,None]


if(self.data['user']['use_cx']):

for p in range(0,nsigmaplus_cx):
Expand Down Expand Up @@ -958,6 +1001,16 @@ def populate_cr_matrix(self):
self.data['cr_matrix']['cr_loss'][i,i,:,:] = self.data['cr_matrix']['cr'][i,i,:,:] - \
np.sum(np.einsum('ij,k->ijk',self.data['rates']['ioniz']['ionization'][i,:,:],
self.data['user']['dens_grid']),axis=0)

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['cr_matrix']['cr'][i,i,:,:] -= np.sum(
self.data['rates']['ioniz']['autoioniz'][i,:]
)
self.data['cr_matrix']['cr_loss'][i,i,:,:] -= np.sum(
self.data['rates']['ioniz']['autoioniz'][i,:]
)

#level i populating mechanisms
#these are the transition rates from higher levels into the level i
self.data['cr_matrix']['cr'][i,0:len(self.data['atomic']['energy']),:,:] = \
Expand Down Expand Up @@ -1030,6 +1083,13 @@ def populate_cr_matrix(self):
np.einsum('ij,k->ijk',self.data['rates']['ioniz']['ionization'][:,p,:],
self.data['user']['dens_grid'])

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['cr_matrix']['cr'][len(self.data['atomic']['energy'])+ p,
0:len(self.data['atomic']['energy']),:,:] += (
self.data['rates']['ioniz']['autoioniz'][:,p]
)[:,None,None]




Expand Down Expand Up @@ -1242,16 +1302,44 @@ def solve_quasi_static(self):
self.data['rates']['ioniz']['ionization'][levels_to_keep,:],
self.data['processed']['F'])

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['processed']['scd'] += np.einsum('ipl,iml->mpl',
np.einsum('ip,l->ipl',
self.data['rates']['ioniz']['autoioniz'][levels_to_keep,:],
1/self.data['user']['dens_grid']
),
self.data['processed']['F']
)

if(self.data['processed']['driving_populations_norm']):
self.data['processed']['scd'] = self.data['processed']['scd'] + \
np.einsum('ipk,ik->ipk',
self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:],
1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:],axis=0)))

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['processed']['scd'] += np.einsum('ipl,il->ipl',
np.einsum('ip,l->ipl',
self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:],
1/self.data['user']['dens_grid']
),
1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:],axis=0))
)

else:
self.data['processed']['scd'] = self.data['processed']['scd'] + \
self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:,:]

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['processed']['scd'] += (
np.einsum('ip,l->ipl',
self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:],
1/self.data['user']['dens_grid']
)
)

if(self.data['user']['use_recombination'] and self.data['user']['use_recombination_three_body']):
#this is the total recombination with three body and the rates that in included in the adf04 file
Expand Down Expand Up @@ -1373,16 +1461,45 @@ def solve_quasi_static(self):
self.data['rates']['ioniz']['ionization'][levels_to_keep,:,:],
self.data['processed']['F'])

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['processed']['scd'] += np.einsum('ipl,imkl->mpkl',
np.einsum('ip,l->ipl',
self.data['rates']['ioniz']['autoioniz'][levels_to_keep,:],
1/self.data['user']['dens_grid']
),
self.data['processed']['F']
)

if(self.data['processed']['driving_populations_norm']):
self.data['processed']['scd'] = self.data['processed']['scd'] + \
np.einsum('ipk,ikl->ipkl',
self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:,:],
1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:,:],axis=0)))

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['processed']['scd'] += np.einsum('ipl,ikl->ipkl',
np.einsum('ip,l->ipl',
self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:],
1/self.data['user']['dens_grid']
),
1/(1+np.sum(self.data['processed']['pops_no_norm'][:,self.data['atomic']['metas'],:,:],axis=0))
)

else:
self.data['processed']['scd'] = self.data['processed']['scd'] + \
self.data['rates']['ioniz']['ionization'][self.data['atomic']['metas'],:,:,None]

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['processed']['scd'] += (
np.einsum('ip,l->ipl',
self.data['rates']['ioniz']['autoioniz'][self.data['atomic']['metas'],:],
1/self.data['user']['dens_grid']
)
)[:,:,None,:]

if(self.data['user']['use_recombination'] and self.data['user']['use_recombination_three_body']):
#this is the total recombination with three body and the rates that in included in the adf04 file
recomb_coeff = np.einsum('ijk,l->ijkl',
Expand Down Expand Up @@ -1535,6 +1652,16 @@ def solve_time_dependent(self):
self.data['processed']['td']['scd'] = np.einsum('ipk,itkl->ptkl',self.data['rates']['ioniz']['ionization'],
self.data['processed']['td']['td_pop'][0:len(self.data['rates']['ioniz']['ionization']),:,:,:])

# Includes autoionization rate if given
if 'autoioniz' in self.data['rates']['ioniz'].keys():
self.data['processed']['scd'] += np.einsum('ipl,itkl->ptkl',
np.einsum('ip,l->ipl',
self.data['rates']['ioniz']['autoioniz'],
1/self.data['user']['dens_grid']
),
self.data['processed']['td']['td_pop'][0:len(self.data['rates']['ioniz']['ionization']),:,:,:]
)


def split_pec_multiplet(self):
"""This function will solve take LS resolved PECs and split them statistically among the j levels
Expand Down Expand Up @@ -2330,7 +2457,7 @@ def plot_pec_ratio_dens(self,pec1,pec2,temp=np.array([0]),meta = np.array([0]),s



def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False):
def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False, wave_lims=None):

""" This function calls the write_adf15 function.

Expand Down Expand Up @@ -2364,7 +2491,8 @@ def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False):
self.data['atomic']['charge_state'], self.data['user']['dens_grid'],
self.data['user']['temp_grid'], self.data['atomic']['metas'],
self.data['atomic']['ion_pot'],
user = self.data['user'], atomic = self.data['atomic'], num = num)
user = self.data['user'], atomic = self.data['atomic'], num = num,
wave_lims=wave_lims)

else:#write un-split PECs

Expand All @@ -2374,11 +2502,13 @@ def write_pecs_adf15(self,fil_name='', pec_inds = 0, num = 8, pecs_split=False):
len(self.data['processed']['wave_air']),dtype='int')

write_adf15(fil_name, pec_inds, self.data['processed']['wave_air']*10,
self.data['processed']['pecs'], self.data['atomic']['element'],
self.data['processed']['pecs'], self.data['processed']['pec_levels'],
self.data['atomic']['element'],
self.data['atomic']['charge_state'], self.data['user']['dens_grid'],
self.data['user']['temp_grid'], self.data['atomic']['metas'],
self.data['atomic']['ion_pot'],
user = self.data['user'], atomic = self.data['atomic'], num = num)
user = self.data['user'], atomic = self.data['atomic'], num = num,
wave_lims=wave_lims)



Expand Down
Loading