Skip to content

Commit

Permalink
Merge pull request #39 from aidancrilly/better_mat_name_convention
Browse files Browse the repository at this point in the history
  • Loading branch information
@aidancrilly
aidancrilly authored Jul 7, 2024
2 parents 88cb4cc + 2454417 commit 81acd6a
Show file tree
Hide file tree
Showing 6 changed files with 212 additions and 226 deletions.
220 changes: 98 additions & 122 deletions example/NeSST Guide.ipynb
View file

Large diffs are not rendered by default.

10 changes: 7 additions & 3 deletions src/NeSST/constants.py
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,9 @@
''' Contains some physical constants used throughout the analysis '''
import scipy.constants as sc
import numpy as np
import numpy.typing as npt

''' Contains some physical constants used throughout the analysis '''
# Scipy constants uses CODATA2018 database

# Need to swap from MeV to eV
amu = sc.value('atomic mass constant energy equivalent in MeV')*1e6
c = sc.c
Expand All @@ -30,4 +29,9 @@
import os
package_directory = os.path.dirname(os.path.abspath(__file__))
data_dir = os.path.join(package_directory,"./data/")
ENDF_dir = os.path.join(data_dir,"./ENDF/")
ENDF_dir = os.path.join(data_dir,"./ENDF/")

# Materials
default_mat_list = ['H','D','T','C12','Be9']
available_materials = []
mat_dict = {}
106 changes: 67 additions & 39 deletions src/NeSST/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,14 +7,15 @@
import numpy.typing as npt
import warnings
import numpy as np
from glob import glob
from os.path import basename
# NeSST libraries
from NeSST.constants import *
import NeSST.collisions as col
import NeSST.spectral_model as sm

# Global variable defaults
col.classical_collisions = False
available_materials = list(sm.available_materials_dict.keys())
# Atomic fraction of D and T in scattering medium and source
frac_D_default = 0.5
frac_T_default = 0.5
Expand All @@ -23,6 +24,33 @@
# Energies, temperatures eV
# Velocities m/s

############################
# Default material load in #
############################
def initialise_material_data(label):
# Aliases
if(label == 'H'):
json = 'H1.json'
elif(label == 'D'):
json = 'H2.json'
elif(label == 'T'):
json = 'H3.json'
# Parse name
else:
mat_jsons = glob(data_dir+'*.json')
mats = [basename(f).split('.')[0] for f in mat_jsons]
if(label in mats):
json = label+'.json'
else:
print("Material label '"+label+"' not recognised")
return
mat_data = sm.material_data(label,json)
mat_dict[label] = mat_data
available_materials.append(label)

for mat in default_mat_list:
initialise_material_data(mat)

##########################################
# Primary spectral shapes & reactivities #
##########################################
Expand Down Expand Up @@ -237,10 +265,10 @@ def init_DT_scatter(Eout : npt.NDArray, Ein : npt.NDArray):
Eout (numpy.array): the array on outgoing neutron energies
"""
sm.mat_D.init_energy_grids(Eout,Ein)
sm.mat_T.init_energy_grids(Eout,Ein)
sm.mat_D.init_station_scatter_matrices()
sm.mat_T.init_station_scatter_matrices()
mat_dict['D'].init_energy_grids(Eout,Ein)
mat_dict['T'].init_energy_grids(Eout,Ein)
mat_dict['D'].init_station_scatter_matrices()
mat_dict['T'].init_station_scatter_matrices()

def init_DT_ionkin_scatter(varr : npt.NDArray, nT: bool = False, nD: bool = False):
"""Initialise the scattering matrices including the effect of ion
Expand All @@ -255,31 +283,31 @@ def init_DT_ionkin_scatter(varr : npt.NDArray, nT: bool = False, nD: bool = Fals
"""
if(nT):
if(sm.mat_T.Ein is None):
if(mat_dict['T'].Ein is None):
print("nT - Needed to initialise energy grids - see init_DT_scatter")
else:
sm.mat_T.full_scattering_matrix_create(varr)
mat_dict['T'].full_scattering_matrix_create(varr)
if(nD):
if(sm.mat_D.Ein is None):
if(mat_dict['D'].Ein is None):
print("nD - Needed to initialise energy grids - see init_DT_scatter")
else:
sm.mat_D.full_scattering_matrix_create(varr)
mat_dict['D'].full_scattering_matrix_create(varr)

def calc_DT_ionkin_primspec_rhoL_integral(I_E : npt.NDArray, rhoL_func=None, nT: bool = False, nD: bool = False):
if(nT):
if(sm.mat_T.vvec is None):
if(mat_dict['T'].vvec is None):
print("nT - Needed to initialise velocity grid - see init_DT_ionkin_scatter")
else:
if(rhoL_func is not None):
sm.mat_T.scattering_matrix_apply_rhoLfunc(rhoL_func)
sm.mat_T.matrix_primspec_int(I_E)
mat_dict['T'].scattering_matrix_apply_rhoLfunc(rhoL_func)
mat_dict['T'].matrix_primspec_int(I_E)
if(nD):
if(sm.mat_D.vvec is None):
if(mat_dict['D'].vvec is None):
print("nD - Needed to initialise velocity grid - see init_DT_ionkin_scatter")
else:
if(rhoL_func is not None):
sm.mat_D.scattering_matrix_apply_rhoLfunc(rhoL_func)
sm.mat_D.matrix_primspec_int(I_E)
mat_dict['D'].scattering_matrix_apply_rhoLfunc(rhoL_func)
mat_dict['D'].matrix_primspec_int(I_E)


###################################
Expand All @@ -297,7 +325,7 @@ def init_mat_scatter(Eout : npt.NDArray, Ein : npt.NDArray, mat_label : str):
"""
mat = sm.available_materials_dict[mat_label]
mat = mat_dict[mat_label]
mat.init_energy_grids(Eout,Ein)
mat.init_station_scatter_matrices()
return mat
Expand Down Expand Up @@ -327,12 +355,12 @@ def DT_sym_scatter_spec(I_E : npt.NDArray
"""
rhoL_func = lambda x : np.ones_like(x)
sm.mat_D.calc_dNdEs(I_E,rhoL_func)
sm.mat_T.calc_dNdEs(I_E,rhoL_func)
nD = frac_D*sm.mat_D.elastic_dNdE
nT = frac_T*sm.mat_T.elastic_dNdE
Dn2n = frac_D*sm.mat_D.n2n_dNdE
Tn2n = frac_T*sm.mat_T.n2n_dNdE
mat_dict['D'].calc_dNdEs(I_E,rhoL_func)
mat_dict['T'].calc_dNdEs(I_E,rhoL_func)
nD = frac_D*mat_dict['D'].elastic_dNdE
nT = frac_T*mat_dict['T'].elastic_dNdE
Dn2n = frac_D*mat_dict['D'].n2n_dNdE
Tn2n = frac_T*mat_dict['T'].n2n_dNdE
total = nD+nT+Dn2n+Tn2n
return total,(nD,nT,Dn2n,Tn2n)

Expand Down Expand Up @@ -363,12 +391,12 @@ def DT_asym_scatter_spec(I_E : npt.NDArray, rhoL_func : callable,
(nD,nT,Dn2n,Tn2n)
"""
sm.mat_D.calc_dNdEs(I_E,rhoL_func)
sm.mat_T.calc_dNdEs(I_E,rhoL_func)
nD = frac_D*sm.mat_D.elastic_dNdE
nT = frac_T*sm.mat_T.elastic_dNdE
Dn2n = frac_D*sm.mat_D.n2n_dNdE
Tn2n = frac_T*sm.mat_T.n2n_dNdE
mat_dict['D'].calc_dNdEs(I_E,rhoL_func)
mat_dict['T'].calc_dNdEs(I_E,rhoL_func)
nD = frac_D*mat_dict['D'].elastic_dNdE
nT = frac_T*mat_dict['T'].elastic_dNdE
Dn2n = frac_D*mat_dict['D'].n2n_dNdE
Tn2n = frac_T*mat_dict['T'].n2n_dNdE
total = nD+nT+Dn2n+Tn2n
return total,(nD,nT,Dn2n,Tn2n)

Expand Down Expand Up @@ -403,20 +431,20 @@ def DT_scatter_spec_w_ionkin(I_E : npt.NDArray, vbar : float, dv : float, rhoL_f
"""
rhoL_func = lambda x : np.ones_like(x)
sm.mat_D.calc_dNdEs(I_E,rhoL_func)
sm.mat_T.calc_dNdEs(I_E,rhoL_func)
if(sm.mat_D.vvec is None):
dNdE_nD = sm.mat_D.elastic_dNdE
mat_dict['D'].calc_dNdEs(I_E,rhoL_func)
mat_dict['T'].calc_dNdEs(I_E,rhoL_func)
if(mat_dict['D'].vvec is None):
dNdE_nD = mat_dict['D'].elastic_dNdE
else:
dNdE_nD = sm.mat_D.matrix_interpolate_gaussian(sm.mat_D.Eout,vbar,dv)
if(sm.mat_T.vvec is None):
dNdE_nT = sm.mat_T.elastic_dNdE
dNdE_nD = mat_dict['D'].matrix_interpolate_gaussian(mat_dict['D'].Eout,vbar,dv)
if(mat_dict['T'].vvec is None):
dNdE_nT = mat_dict['T'].elastic_dNdE
else:
dNdE_nT = sm.mat_T.matrix_interpolate_gaussian(sm.mat_T.Eout,vbar,dv)
dNdE_nT = mat_dict['T'].matrix_interpolate_gaussian(mat_dict['T'].Eout,vbar,dv)
nD = frac_D*dNdE_nD
nT = frac_T*dNdE_nT
Dn2n = frac_D*sm.mat_D.n2n_dNdE
Tn2n = frac_T*sm.mat_T.n2n_dNdE
Dn2n = frac_D*mat_dict['D'].n2n_dNdE
Tn2n = frac_T*mat_dict['T'].n2n_dNdE
total = nD+nT+Dn2n+Tn2n
return total,(nD,nT,Dn2n,Tn2n)

Expand Down Expand Up @@ -467,7 +495,7 @@ def calc_DT_sigmabar(Ein : npt.NDArray, I_E : npt.NDArray,
Returns:
float : the spectrally averaged total DT cross section
"""
sigmabar = frac_D*np.trapz(sm.mat_D.sigma_tot(Ein)*I_E,Ein)+frac_T*np.trapz(sm.mat_T.sigma_tot(Ein)*I_E,Ein)
sigmabar = frac_D*np.trapz(mat_dict['D'].sigma_tot(Ein)*I_E,Ein)+frac_T*np.trapz(mat_dict['T'].sigma_tot(Ein)*I_E,Ein)
return sigmabar

def rhoR_2_A1s(rhoR : typing.Union[float,npt.NDArray],
Expand Down
68 changes: 34 additions & 34 deletions src/NeSST/fitting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,25 +60,25 @@ def init_symmetric_model(self):

if(self.ion_kinematics):
calc_DT_ionkin_primspec_rhoL_integral(self.dNdE_DT,nT=True,nD=True)
sm.mat_D.calc_n2n_dNdE(self.dNdE_DT,rhoL_func)
sm.mat_T.calc_n2n_dNdE(self.dNdE_DT,rhoL_func)
mat_dict['D'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func)
mat_dict['T'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func)
else:
sm.mat_D.calc_station_elastic_dNdE(self.dNdE_DT,rhoL_func)
sm.mat_T.calc_station_elastic_dNdE(self.dNdE_DT,rhoL_func)
sm.mat_D.calc_n2n_dNdE(self.dNdE_DT,rhoL_func)
sm.mat_T.calc_n2n_dNdE(self.dNdE_DT,rhoL_func)
mat_dict['D'].calc_station_elastic_dNdE(self.dNdE_DT,rhoL_func)
mat_dict['T'].calc_station_elastic_dNdE(self.dNdE_DT,rhoL_func)
mat_dict['D'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func)
mat_dict['T'].calc_n2n_dNdE(self.dNdE_DT,rhoL_func)

dNdE_Dn2n = interpolate_1d(self.E_sspec,sm.mat_D.n2n_dNdE,fill_value=0.0,bounds_error=False)
dNdE_Tn2n = interpolate_1d(self.E_sspec,sm.mat_T.n2n_dNdE,fill_value=0.0,bounds_error=False)
dNdE_Dn2n = interpolate_1d(self.E_sspec,mat_dict['D'].n2n_dNdE,fill_value=0.0,bounds_error=False)
dNdE_Tn2n = interpolate_1d(self.E_sspec,mat_dict['T'].n2n_dNdE,fill_value=0.0,bounds_error=False)

if(self.ion_kinematics):
def model(E,rhoL,vbar,dv,fT,fD,Yn):
"""
Symmetric areal density model with scattering ion velocity distribution with mean and std dev, vbar and dv in m/s
"""
A_1S = rhoR_2_A1s(rhoL,frac_D=fD,frac_T=fT)
dNdE_nT = sm.mat_T.matrix_interpolate_gaussian(E,vbar,dv)
dNdE_nD = sm.mat_D.matrix_interpolate_gaussian(E,vbar,dv)
dNdE_nT = mat_dict['T'].matrix_interpolate_gaussian(E,vbar,dv)
dNdE_nD = mat_dict['D'].matrix_interpolate_gaussian(E,vbar,dv)
dNdE_tot = A_1S*(fT*dNdE_nT+fD*dNdE_nD+fD*dNdE_Dn2n(E)+fT*dNdE_Tn2n(E))
return Yn*(dNdE_tot+(fD/fT)*(frac_T_default/frac_D_default)*self.I_DD(E)+(fT/fD)*(frac_D_default/frac_T_default)*self.I_TT(E))
else:
Expand All @@ -88,8 +88,8 @@ def model(E,rhoL,Ts,fT,fD,Yn):
Symmetric areal density model with scattering temperature Ts, in keV
"""
A_1S = rhoR_2_A1s(rhoL,frac_D=fD,frac_T=fT)
dNdE_nT = sm.mat_T.elastic_dNdE.copy()
dNdE_nD = sm.mat_D.elastic_dNdE.copy()
dNdE_nT = mat_dict['T'].elastic_dNdE.copy()
dNdE_nD = mat_dict['D'].elastic_dNdE.copy()
if(Ts > 0.1):
T_MeV = Ts/1e3
E_nT0 = ((sm.A_T-1.0)/(sm.A_T+1.0))**2*self.DTmean
Expand All @@ -114,44 +114,44 @@ def init_modeone_model(self,P1_arr):
self.P1_arr = P1_arr

# T(n,2n)
sm.mat_T.n2n_ddx.rgrid_IE = np.trapz(sm.mat_T.n2n_ddx.rgrid*self.dNdE_DT[:,None,None],self.E_DTspec,axis=0)
sm.mat_T.n2n_dNdE_mode1 = np.trapz(sm.mat_T.n2n_ddx.rgrid_IE[:,:,None]*
(1.0+self.P1_arr[None,None,:]*sm.mat_T.n2n_mu[:,None,None])
,sm.mat_T.n2n_mu,axis=0)
mat_dict['T'].n2n_ddx.rgrid_IE = np.trapz(mat_dict['T'].n2n_ddx.rgrid*self.dNdE_DT[:,None,None],self.E_DTspec,axis=0)
mat_dict['T'].n2n_dNdE_mode1 = np.trapz(mat_dict['T'].n2n_ddx.rgrid_IE[:,:,None]*
(1.0+self.P1_arr[None,None,:]*mat_dict['T'].n2n_mu[:,None,None])
,mat_dict['T'].n2n_mu,axis=0)

sm.mat_T.n2n_dNdE_mode1 = interpolate_2d(self.E_sspec,self.P1_arr,sm.mat_T.n2n_dNdE_mode1,bounds_error=False)
mat_dict['T'].n2n_dNdE_mode1 = interpolate_2d(self.E_sspec,self.P1_arr,mat_dict['T'].n2n_dNdE_mode1,bounds_error=False)

# D(n,2n)
sm.mat_D.n2n_ddx.rgrid_IE = np.trapz(sm.mat_D.n2n_ddx.rgrid*self.dNdE_DT[:,None,None],self.E_DTspec,axis=0)
sm.mat_D.n2n_dNdE_mode1 = np.trapz(sm.mat_D.n2n_ddx.rgrid_IE[:,:,None]*
(1.0+self.P1_arr[None,None,:]*sm.mat_D.n2n_mu[:,None,None])
,sm.mat_D.n2n_mu,axis=0)
mat_dict['D'].n2n_ddx.rgrid_IE = np.trapz(mat_dict['D'].n2n_ddx.rgrid*self.dNdE_DT[:,None,None],self.E_DTspec,axis=0)
mat_dict['D'].n2n_dNdE_mode1 = np.trapz(mat_dict['D'].n2n_ddx.rgrid_IE[:,:,None]*
(1.0+self.P1_arr[None,None,:]*mat_dict['D'].n2n_mu[:,None,None])
,mat_dict['D'].n2n_mu,axis=0)

sm.mat_D.n2n_dNdE_mode1 = interpolate_2d(self.E_sspec,self.P1_arr,sm.mat_D.n2n_dNdE_mode1,bounds_error=False)
mat_dict['D'].n2n_dNdE_mode1 = interpolate_2d(self.E_sspec,self.P1_arr,mat_dict['D'].n2n_dNdE_mode1,bounds_error=False)

# nT
M_mode1 = np.trapz((1.0+self.P1_arr[None,None,None,:]*sm.mat_T.full_scattering_mu[:,:,:,None])
*sm.mat_T.full_scattering_M[:,:,:,None]
M_mode1 = np.trapz((1.0+self.P1_arr[None,None,None,:]*mat_dict['T'].full_scattering_mu[:,:,:,None])
*mat_dict['T'].full_scattering_M[:,:,:,None]
*self.dNdE_DT[None,None,:,None],self.E_DTspec,axis=2)

sm.mat_T.M_mode1_interp = interpolate_1d(self.P1_arr,M_mode1,axis=-1,bounds_error=False)
mat_dict['T'].M_mode1_interp = interpolate_1d(self.P1_arr,M_mode1,axis=-1,bounds_error=False)

# nD
M_mode1 = np.trapz((1.0+self.P1_arr[None,None,None,:]*sm.mat_D.full_scattering_mu[:,:,:,None])
*sm.mat_D.full_scattering_M[:,:,:,None]
M_mode1 = np.trapz((1.0+self.P1_arr[None,None,None,:]*mat_dict['D'].full_scattering_mu[:,:,:,None])
*mat_dict['D'].full_scattering_M[:,:,:,None]
*self.dNdE_DT[None,None,:,None],self.E_DTspec,axis=2)

sm.mat_D.M_mode1_interp = interpolate_1d(self.P1_arr,M_mode1,axis=-1,bounds_error=False)
mat_dict['D'].M_mode1_interp = interpolate_1d(self.P1_arr,M_mode1,axis=-1,bounds_error=False)

if(self.ion_kinematics):
def model(E,rhoL,P1,vbar,dv,fT,fD,Yn):
A_1S = rhoR_2_A1s(rhoL,frac_D=fD,frac_T=fT)
sm.mat_T.M_prim = sm.mat_T.M_mode1_interp(P1)
sm.mat_D.M_prim = sm.mat_D.M_mode1_interp(P1)
dNdE_nT = sm.mat_T.matrix_interpolate_gaussian(E,vbar,dv)
dNdE_nD = sm.mat_D.matrix_interpolate_gaussian(E,vbar,dv)
dNdE_Tn2n = sm.mat_T.n2n_dNdE_mode1(E,P1)
dNdE_Dn2n = sm.mat_D.n2n_dNdE_mode1(E,P1)
mat_dict['T'].M_prim = mat_dict['T'].M_mode1_interp(P1)
mat_dict['D'].M_prim = mat_dict['D'].M_mode1_interp(P1)
dNdE_nT = mat_dict['T'].matrix_interpolate_gaussian(E,vbar,dv)
dNdE_nD = mat_dict['D'].matrix_interpolate_gaussian(E,vbar,dv)
dNdE_Tn2n = mat_dict['T'].n2n_dNdE_mode1(E,P1)
dNdE_Dn2n = mat_dict['D'].n2n_dNdE_mode1(E,P1)
dNdE_tot = A_1S*(fT*dNdE_nT+fD*dNdE_nD+fD*dNdE_Dn2n+fT*dNdE_Tn2n)
# Primary
dNdE_DD = (fD/fT)*(frac_T_default/frac_D_default)*self.I_DD(E)
Expand Down
25 changes: 2 additions & 23 deletions src/NeSST/spectral_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,21 +24,10 @@ def unity(x):

class material_data:

def __init__(self,label):
def __init__(self,label,json):
self.label = label
if(self.label == 'H'):
self.json = 'H1.json'
elif(self.label == 'D'):
self.json = 'H2.json'
elif(self.label == 'T'):
self.json = 'H3.json'
elif(self.label == '12C'):
self.json = 'C12.json'
elif(self.label == '9Be'):
self.json = 'Be9.json'
else:
print("Material label "+self.label+" not recognised")

self.json = json
ENDF_data = retrieve_ENDF_data(self.json)

self.A = ENDF_data['A']
Expand Down Expand Up @@ -244,16 +233,6 @@ def matrix_interpolate_gaussian(self,E,vbar,dv):
interp = interpolate_1d(self.Eout,M_v,method='linear',bounds_error=False)
return interp(E)


# Default load
mat_H = material_data('H')
mat_D = material_data('D')
mat_T = material_data('T')
mat_9Be = material_data('9Be')
mat_12C = material_data('12C')

available_materials_dict = {"H" : mat_H, "D" : mat_D, "T" : mat_T, "12C" : mat_12C, "9Be" : mat_9Be}

# Load in TT spectrum
# Based on Appelbe, stationary emitter, temperature range between 1 and 10 keV
# https://www.sciencedirect.com/science/article/pii/S1574181816300295
Expand Down
Loading

0 comments on commit 81acd6a

Please sign in to comment.