From b6dff9e88c87702d1ffafe82c5d8ba110daa4be7 Mon Sep 17 00:00:00 2001 From: Matt Kriete Date: Tue, 2 Jul 2024 18:23:12 -0400 Subject: [PATCH 1/2] Can now pull charge exchange rate coefficients from ADF11 files for performing ionization balance --- colradpy/ionization_balance_class.py | 28 ++++++++++++++++------------ colradpy/read_adf11.py | 16 ++++++++-------- 2 files changed, 24 insertions(+), 20 deletions(-) diff --git a/colradpy/ionization_balance_class.py b/colradpy/ionization_balance_class.py index 5fd2b40..696b527 100644 --- a/colradpy/ionization_balance_class.py +++ b/colradpy/ionization_balance_class.py @@ -35,16 +35,18 @@ class ionization_balance(): :param metas: List of arrays for the metastable levels in a charge state :type metas: list - :param temp_grid: Array of the temperatures in (eV) + :param temp_grid: Array of electron temperatures in (eV) :type temp_grid: float array - :param dens_grid: Array of the densities in (cm-3) + :param dens_grid: Array of electron densities in (cm-3) :type dens_grid: float array - :param htemp_grid: Temperature grid of thermal CX hydrogen (eV) + :param htemp_grid: Temperature grid of thermal CX hydrogen (eV). Must be + aligned with the electron temperature grid. :type htemp_grid: float array - :param hdens_grid: Density grid of the thermal CX hydrogen densities in (cm-3) + :param hdens_grid: Density grid of the thermal CX hydrogen densities in + (cm^-3). Must be aligned with the electron density grid. :type hdens_grid: float array :param soln_times: Times to calculate the solution for the time dependent solutions (s) @@ -120,7 +122,9 @@ def __init__(self,fils,metas,temp_grid, dens_grid, htemp_grid=np.array([]), hden self.data['cr_data']['gcrs'] = {} self.data['user'] = {} self.data['user']['temp_grid'] = np.asarray(temp_grid) #eV - self.data['user']['dens_grid'] = np.asarray(dens_grid)#cm-3 + self.data['user']['dens_grid'] = np.asarray(dens_grid) #cm^-3 + self.data['user']['htemp_grid'] = np.asarray(htemp_grid) #eV + self.data['user']['hdens_grid'] = np.asarray(hdens_grid) # cm^-3 self.data['user']['fils'] = np.asarray(fils) self.data['user']['init_abund'] = np.asarray(init_abund) self.data['user']['soln_times'] = np.asarray(soln_times) @@ -147,9 +151,9 @@ def __init__(self,fils,metas,temp_grid, dens_grid, htemp_grid=np.array([]), hden if(type(use_cx) == bool): - self.data['user']['use_cx'] = np.ones(len(adf11['input_file']['metas'])-1,dtype=bool) - self.data['user']['use_cx'][:] = False - + self.data['user']['use_cx'] = np.ones(len(adf11['input_file']['metas'])-1, dtype=bool) + self.data['user']['use_cx'][:] = use_cx + for j in range(0,len(adf11['input_file']['metas'])-1): if( str(j) not in self.data['cr_data']['gcrs']): @@ -168,20 +172,20 @@ def __init__(self,fils,metas,temp_grid, dens_grid, htemp_grid=np.array([]), hden self.data['input_file']['acd'] = adf11['input_file'] if( 'qcd' in fils[i]): - self.data['cr_data']['gcrs'][str(j)]['qcd']= interp_rates_adf11(adf11['input_file']['temp_grid'], + self.data['cr_data']['gcrs'][str(j)]['qcd'] = interp_rates_adf11(adf11['input_file']['temp_grid'], adf11['input_file']['dens_grid'], temp_grid,dens_grid,adf11['input_file'][str(j)]) self.data['input_file']['qcd'] = adf11['input_file'] if( 'xcd' in fils[i]): - self.data['cr_data']['gcrs'][str(j)]['xcd']= interp_rates_adf11(adf11['input_file']['temp_grid'], + self.data['cr_data']['gcrs'][str(j)]['xcd'] = interp_rates_adf11(adf11['input_file']['temp_grid'], adf11['input_file']['dens_grid'], temp_grid,dens_grid,adf11['input_file'][str(j)]) self.data['input_file']['xcd'] = adf11['input_file'] if( 'ccd' in fils[i]): - self.data['cr_data']['gcrs'][str(j)]['ccd']= interp_rates_adf11(adf11['input_file']['temp_grid'], + self.data['cr_data']['gcrs'][str(j)]['ccd'] = interp_rates_adf11(adf11['input_file']['temp_grid'], adf11['input_file']['dens_grid'], - temp_grid,dens_grid,adf11['input_file'][str(j)]) + htemp_grid,hdens_grid,adf11['input_file'][str(j)]) self.data['input_file']['ccd'] = adf11['input_file'] if('qcd' not in self.data['cr_data']['gcrs'][str(j)]): diff --git a/colradpy/read_adf11.py b/colradpy/read_adf11.py index b49c7c5..5b781a9 100644 --- a/colradpy/read_adf11.py +++ b/colradpy/read_adf11.py @@ -25,8 +25,8 @@ def read_adf11(fil): adf11['input_file']['charge_max'] = int(tmp[4]) f.readline() #reading '-------------' - if('r' in re.split('_', os.path.split(fil)[-1])[0]): - + file_type = re.split('_', os.path.split(fil)[-1])[0] + if('r' in file_type): adf11['input_file']['metas'] = np.array(list( #metastables map(int,re.findall('(\d+)',f.readline())))) f.readline() #reading '---------------' @@ -49,32 +49,32 @@ def read_adf11(fil): num_stages = len(adf11['input_file']['metas']) for i in range(0,num_stages): adf11['input_file'][str(i)] = {} - if( 'scd' in fil or 'acd' in fil): + if any(x in file_type for x in ['scd', 'acd', 'ccd']): adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i], adf11['input_file']['metas'][i+1], len(adf11['input_file']['temp_grid']), len(adf11['input_file']['dens_grid']))) - if( 'qcd' in fil ): + if( 'qcd' in file_type ): adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i], adf11['input_file']['metas'][i], len(adf11['input_file']['temp_grid']), len(adf11['input_file']['dens_grid']))) - if('xcd' in fil): + if('xcd' in file_type): adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i+1], adf11['input_file']['metas'][i+1], len(adf11['input_file']['temp_grid']), len(adf11['input_file']['dens_grid']))) - if('plt' in fil): + if('plt' in file_type): adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i], adf11['input_file']['metas'][i], len(adf11['input_file']['temp_grid']), len(adf11['input_file']['dens_grid']))) - if('prb' in fil): + if('prb' in file_type): adf11['input_file'][str(i)] = np.zeros((adf11['input_file']['metas'][i], adf11['input_file']['metas'][i], len(adf11['input_file']['temp_grid']), len(adf11['input_file']['dens_grid']))) - + #Reading the GCR value portion gcr_line = f.readline() ii = 0 From 2e8643d5428fc3e4c27304563b0cf02522155c9d Mon Sep 17 00:00:00 2001 From: Matt Kriete Date: Wed, 3 Jul 2024 13:06:09 -0400 Subject: [PATCH 2/2] Energy balance can now include CX in the calculation of the underlying ionization balance --- colradpy/energy_balance.py | 32 +++++++++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 3 deletions(-) diff --git a/colradpy/energy_balance.py b/colradpy/energy_balance.py index d460cb0..b4bd5ab 100644 --- a/colradpy/energy_balance.py +++ b/colradpy/energy_balance.py @@ -24,6 +24,11 @@ class energy_balance(object): densities/temperatures at a time (it is not vectorized across a grid of temperatures/densities like the ionization balance). + Currently, if charge exchange effects are turned on, they only directly + affect the ionization balance. The energy balance is then indirectly + affected through the ionization balance. The direct effects of charge + exchange on the energy balance are not modelled. + :Args: :param fils: Array of the input adf04 or adf11 files :type fils: string array @@ -54,7 +59,16 @@ class energy_balance(object): :param ion_charge_number: Main ion species charge number :type ion_charge_number: int - + + :param use_cx: Include thermal charge exchange in the ionization balance + :param use_cx: bool + + :param hydrogen_temp: Neutral hydrogen temperature in (eV) + :type hydrogen_temp: float + + :param hydrogen_dens: Neutral hydrogen density in (cm^-3) + :type hydrogen_dens: float + :param init_temp: Initial temperatures in (eV) of the charge states for the modelled species :type init_temp: float array @@ -63,12 +77,19 @@ class energy_balance(object): def __init__(self, fils, metas, mass, polarizability, electron_temp, electron_dens, ion_temp, ion_dens, ion_mass, ion_charge_number, + use_cx=False, hydrogen_temp=None, hydrogen_dens=None, init_temp=np.array([]), **kwargs): # Set up ionization balance self.ion_balance = ionization_balance( - fils, metas, temp_grid=np.array([electron_temp]), - dens_grid=np.array([electron_dens]), **kwargs, + fils, + metas, + temp_grid=np.array([electron_temp]), + dens_grid=np.array([electron_dens]), + htemp_grid=np.array([hydrogen_temp]), + hdens_grid=np.array([hydrogen_dens]), + use_cx=use_cx, + **kwargs, ) # Initialize data @@ -83,6 +104,10 @@ def __init__(self, fils, metas, mass, polarizability, electron_temp, self.data["user"]["ion_mass"] = ion_mass self.data["user"]["ion_charge_number"] = ion_charge_number self.data["user"]["init_temp"] = init_temp + self.data["user"]["use_cx"] = use_cx + if use_cx: + self.data["user"]["hydrogen_temp"] = hydrogen_temp + self.data["user"]["hydrogen_dens"] = hydrogen_dens def solve(self, n0=np.array([]), T0=np.array([]), td_t=np.array([])): @@ -169,6 +194,7 @@ def _evolve_energy(self, time, energies, pops_fun): Returns a vector of the time derivative of the energy in each charge state """ + # TODO: Implement effects of charge exchange reactions on energy balance # Get fractional abundance for this time and calculate state temperatures pops = pops_fun(np.array([time]))[:, 0, 0, 0] # Handle special cases where populations are zero (typically at beginning)