Skip to content

Particle Source

ecrouse3 edited this page May 12, 2023 · 17 revisions

A particle source can be created in the setup.cfg file or with an outside code that can write a particleSource.nc file that the setup.cfg can read from. The script that creates the particleSource.nc file is also responsible for finding the sputtered flux and impurity flux of whatever surface and plasma are being analyzed. In the following example, the interaction surface is made up of tungsten (W).

Methods

.cfg method

This method requires adding a section to the impurityParticleSource section of the setup.cfg file. The initial conditions for the particles is included under impurityParticleSource.initialConditions as shown in the following block.

impurityParticleSource = 
{
    nP = <number of particles>;
    sourceStrength = <source strength>;
    Z = <number of protons>;
    source_material_Z = <number of protons of surface>;
    source_material_SurfaceBindingEnergy = <Binding energy of surface>;

	initialConditions = 
	{
		x_start = <x starting location in meters>;
		y_start = <y starting location>;
		z_start = <z starting location>;
		energy_eV_x_start = <x starting energy in eV>;
		energy_eV_y_start = <y starting energy>;
		energy_eV_z_start = <z starting energy>;
		impurity_amu = <atomic mass of impurity>;
		impurity_Z = <number of protons of impurity>;
	        charge = <initial charge in e>;
                energy_eV = <initial energy>;
                phi = <initial phi angle in degrees>;
                theta = <initial theta angle in degrees>;
    };
    ionization = {remains the same};
    recombination = {reamains the same};
};

This method is good for testing a system quickly.

.nc method

When using a particleSource.nc file to store the information for a particle source, it must be filled in using Python. The following is a simple point source code that will yield the same results as using the setup.cfg method.

def point_source(nP):
    x = <x starting location in meters>*np.ones(nP)
    y = <y starting location>*np.ones(nP)
    z = <z starting location>*np.ones(nP)
    vx = <x starting velocity in m/s>*np.ones(nP)
    vy = <y starting velocity>*np.ones(nP)
    vz = <z starting velocity>*np.ones(nP)

    #########################################
    #make NetCDF Particle Source file
    #########################################

    rootgrp = netCDF4.Dataset("particleSource.nc", "w", format="NETCDF4")
    npp = rootgrp.createDimension("nP", nP)
    xxx = rootgrp.createVariable("x","f8",("nP"))
    yyy = rootgrp.createVariable("y","f8",("nP"))
    zzz = rootgrp.createVariable("z","f8",("nP"))
    vxx = rootgrp.createVariable("vx","f8",("nP"))
    vyy = rootgrp.createVariable("vy","f8",("nP"))
    vzz = rootgrp.createVariable("vz","f8",("nP"))
    xxx[:] = x
    yyy[:] = y
    zzz[:] = z
    vxx[:] = vx
    vyy[:] = vy
    vzz[:] = vz
    rootgrp.close()

The particleSource.nc can be implemented in the setup.cfg by adding the following block. The inclusion of method = 1 ensures the particleSource.nc file is used.

particleSource =
{
    ncFileString = "particleSource.nc";
};
impurityParticleSource =
{   method=1;
    nP = <number of particles>;
    sourceStrength = <source strength>;
    Z = <Z value>;
    source_material_Z = <Z value of surface>;
	initialConditions = 
	{
		impurity_amu = <atomic mass of impurity>;
		impurity_Z = <Z value of impurity>;
	        charge = <charge in e>;
                energy_eV = <initial energy>;
                phi = <initial phi angle in degrees>;
                theta = <initial theta angle in degrees>;
    };
    ionization = {remains the same};
    recombination = {reamains the same};
};

Flux and Energy

To determine the yields of sputtering, the particle properties and the surface material must be defined. This starts by defining the interaction material's indices out of the whole geometry. This could be the indices where the divertor is or wherever particle interactions take place.

These values must be user determined as they will be called later in the script.

nP = <number of particles> 
surfW = <array of tungsten surface incdices>
geom = <geometry .cfg file>
profiles_file = <plasma profiles .nc file>
gitr_rz = <rz .txt file>
rmrs_fine_file = <r minus r sep .txt file>
W_fine_file = <tungsten fine indices .txt file>
ftFile = <ftridynBackground for desired particle>
configuration = <midpoint,uniform, or random>

#import wall geometry to plot over
with open(gitr_rz, 'r') as file:
    wall = file.readlines()

#import refined rmrs at the W surface
with open(rmrs_fine_file, 'r') as file:
    rmrs_fine = file.readlines()   
rmrsFine = np.array(rmrs_fine,dtype='float')

#import W surface indices
with open(W_fine_file, 'r') as file:
    W_fine = file.readlines()
W_fine = np.array(W_fine,dtype='int')

#import coarse rmrs at W surface
profiles = netCDF4.Dataset(profiles_file)
rmrsCoarse = profiles.variables['rmrs_inner_target'][surfW]

The r and z location of each point on the interaction and noninteraction surface must be listed.

R = np.zeros(len(wall))
Z = np.zeros(len(wall))
for i,line in enumerate(wall):
    point = line.split()
    R[i] = float(point[0])
    Z[i] = float(point[1])

R = R[W_fine]
Z = Z[W_fine]

r1 = R[:-1]
r2 = R[1:]
z1 = Z[:-1]
z2 = Z[1:]

The following code can be used to estimate the energy and angle of a particle as it reaches a surface. This code uses plasmaProfiles.nc saved as profiles to determine the electron and ion temperature at each index of the interacted surface. The sputtering yields can also be found by comparing energies with fractal tridyn tables.

def get_incoming_IEADs(q, profiles, surfW, rmrsCoarse, rmrsFine):
    #extract plasma parameters at the wall indices
    teCoarse = profiles['te_target'][surfW]
    tiCoarse = profiles['ti_target'][surfW]
    
    #get temp as a function of rmrsCoarse
    fte = scii.interp1d(rmrsCoarse, teCoarse, fill_value='extrapolate')
    fti = scii.interp1d(rmrsCoarse, tiCoarse, fill_value='extrapolate')
    
    #interpolate temp at points in rmrsFine
    te = fte(rmrsFine)
    ti = fti(rmrsFine)

    SimpleEnergyEst = 2*ti+3*te*q

    AngleEst = 70*np.ones(len(SimpleEnergyEst))

    return SimpleEnergyEst, AngleEst

def get_ft_spyld(S, surfE, surfA, ftBFile):
    #import sputtering yield tables for incident ions on W
    ftB = netCDF4.Dataset(ftBFile, "r", format="NETCDF4")
    spyld = ftB.variables['spyld'][S][:]
    ftE = ftB.variables['E'][:]
    ftA = ftB.variables['A'][:]
    
    surfY = scii.interpn((ftE,ftA), spyld, (surfE,surfA))
    
    return surfY

#get sputtering yields for particles from fractal tridyn tables
spyld = get_ft_spyld(<Ion Charge>, <energy of particle>, <angle of particle>, <fractal tridyn tables>)

The analytical sputtering yield can be used instead of the fractal tridyn tables. It uses the Eckstein formula to obtain the sputtering yield.

def get_analytic_spyld(surfE, surfA, Z1=6, M1=12, Z2=74, M2=183.84, \
                       FitParam='N', Eth=45.3362, lam=0.0921, q=1.4389, mu=2.0225):
    #this entire function is simply the Eckstein formula
    #defaults are for C on W but with fitting parameters, Eth, and Esp for N on W
    #energies are in eV and angles are in degrees
    #M1, Z1 are for the projectile
    #M2, Z2 are for the target
    #Eth is the threshold energy for any sputtering to occur
    #Esp is the SBE for self-bombardment
    #lam, q, mu, f, b, c are fitting parameters from Eckstein tables

    e = 14.399651

    a_L = 0.0529177*((9*np.pi**2/128)**(1/3)) + \
        (Z1**(2/3) + Z2**(2/3))**-0.5

    epsilon_L = surfE * (M2/(M1+M2)) * (a_L/(Z1*Z2*e**2))

    omega = epsilon_L + 0.1728*np.sqrt(epsilon_L) + \
        0.008*epsilon_L**0.1504

    snKrC = 0.5*np.log(1+1.2288*epsilon_L)/omega
    
    #spyld(E, angle = 0) normal incidence
    Y_0 = np.zeros(len(surfE))
    for i,v in enumerate(surfE):
        if v>=Eth:
            Y_0[i] = q*snKrC[i]*(v/Eth-1)**mu / \
                (lam/omega[i] + (v/Eth-1)**mu)
    
    #choose basis for fitting parameters for different energies
    Esp, f, b, c, ThetaMax = FittingParameters_NonW(surfE)

    theta_0 = np.pi - np.arccos(np.sqrt(1/(1+surfE/Esp)))

    #spyld(E,A)    
    Y = Y_0 * (np.cos((np.pi*surfA/(2*theta_0))**c))**(-f) * \
        np.exp(b*(1-1/np.cos((np.pi*surfA/(2*theta_0))**c)))

    return Y_0

The ion charge can be changed to simulate different types of particles. FittingParameters_NonW is a function that estimates the angle of approach for nitrogen on tungsten.

def FittingParameters_NonW(E):
    Eth=45.3362
    Esp = np.ones(len(E))
    f = np.zeros(len(E))
    b = np.zeros(len(E))
    c = np.zeros(len(E))
    ThetaMax = np.zeros(len(E))
    for i,v in enumerate(E):
        if v<Eth:
            continue
        elif v>=Eth and v<49:
            f[i],b[i],c[i] = 2.9557, 5.8879, 0.9465
        elif v>=49 and v<51:
            f[i],b[i],c[i] = 1.7735, 4.3144, 0.9468
        elif v>=51 and v<53.5:
            f[i],b[i],c[i] = 1.2707, 3.6458, 0.8840
        elif v>=53.5 and v<57.5:
            f[i],b[i],c[i] = 1.1002, 3.3751, 0.9604
        elif v>=57.5 and v<65:
            f[i],b[i],c[i] = 0.4622, 2.5095, 1.0118
        elif v>=65 and v<75:
            f[i],b[i],c[i] = 0.4622, 2.5095, 1.0118
        elif v>=75 and v<85:
            f[i],b[i],c[i] = 2.7960, 3.4029, 0.8841
        elif v>=85 and v<95:
            f[i],b[i],c[i] = 2.1152, 2.6541, 0.9226
        elif v>=95 and v<110:
            f[i],b[i],c[i] = 1.7312, 2.1735, 0.9489
        elif v>=110 and v<130:
            f[i],b[i],c[i] = 1.6230, 1.6737, 1.0004
        elif v>=130 and v<170:
            f[i],b[i],c[i] = 1.7195, 1.5092, 1.0176
            ThetaMax[i] = 30.54
        elif v>=170 and v<250:
            f[i],b[i],c[i] = 2.0138, 1.3460, 1.0316
            ThetaMax[i] = 51.98
        elif v>=250 and v<400:
            f[i],b[i],c[i] = 2.2531, 1.2151, 1.0310
            ThetaMax[i] = 59.47
        elif v>=400 and v<750:
            f[i],b[i],c[i] = 2.4324, 1.1313, 1.0171
            ThetaMax[i] = 62.20
        elif v>=750 and v<2000:
            f[i],b[i],c[i] = 2.4383, 0.9940, 0.9936
            ThetaMax[i] = 66.00
        else:
            print('WARNING: TABULAR SPUTTERING FITTING PARAMETERS MISSING')
    return Esp, f, b, c, ThetaMax

#get coarse flux profile from background D, C and refine to rmrsFine
fluxCoarse = np.abs(profiles.variables['flux_inner_target'][<Ion Charge>][<surface of interaction indices>])
fflux = scii.interp1d(rmrsCoarse,fluxCoarse,fill_value='extrapolate')
flux = fflux(rmrsFine)

#multiply incoming ion flux by Y_s to get sputtered W flux by each species
sputt_flux = spyld*flux

# find the total sputtered flux
sputt_flux_total = <Sum of calculated flux>

#multiply by area to get the outgoing particles per second
pps = np.multiply(sputt_flux_total,area)
pps_weights = nP*pps/np.sum(pps)

The sputtering yield, incoming flux, and sputtered flux can be determined for each particle type by repeating this process with different ion properties.


The particles must have an initial starting location and starting velocity. Since the location of particles along the wall can vary, there is a portion of the program that can decide whether the distribution of the surface particles in between each line segment is random, uniform, or they all start in the midpoint of the segment. The functions for these distribution settings are below.

def midpoints(nP,pps_weights,adj,slope,Beta, r1,z1,r2,z2):
    #get midpoints of coords
    r_mid = np.zeros(len(r1)-1)
    z_mid = np.zeros(len(z1)-1)
    for i in range(len(r1)):
        r_mid[i] = np.average(np.array([r1[i],r2[i]]))
        z_mid[i] = np.average(np.array([z1[i],z2[i]]))
        
    x = np.zeros(nP)
    y = np.zeros(nP)
    z = np.zeros(nP)
    counter = 0
    for i in range(len(pps_weights)):
        x[counter:counter+pps_weights[i]] = r_mid[i] - adj*np.abs(np.cos(Beta[i]))
        z[counter:counter+pps_weights[i]] = z_mid[i] + np.sign(slope[i]) * adj*np.abs(np.sin(Beta[i]))
        counter += pps_weights[i]

    return x,y,z

def uniform(nP,pps_weights,adj,slope,Beta, r1,z1,r2,z2):
    x = np.zeros(nP)
    y = np.zeros(nP)
    z = np.zeros(nP)
    tally = 0
    for i in range(len(pps_weights)):
        if pps_weights[i]>0:
            dr = (r2[i]-r1[i])/pps_weights[i]
            dz = (z2[i]-z1[i])/pps_weights[i]
            r_segment = np.linspace(r1[i]+dr/2, r2[i]-dr/2, pps_weights[i])
            z_segment = np.linspace(z1[i]+dz/2, z2[i]-dz/2, pps_weights[i])
            
            x[tally:tally+pps_weights[i]] = r_segment - adj*np.abs(np.cos(Beta[i]))
            z[tally:tally+pps_weights[i]] = z_segment + np.sign(slope[i]) * adj*np.abs(np.sin(Beta[i]))
            tally += pps_weights[i]

    return x,y,z

def random(nP,pps_weights,adj,slope,Beta, r1,z1,r2,z2):
    x = np.zeros(nP)
    y = np.zeros(nP)
    z = np.zeros(nP)
    counter = 0
    for i in range(len(pps_weights)):
        for j in range(pps_weights[i]):
            chi = np.random.rand(1)
            x[counter+j] = r1[i]+chi*(r2[i]-r1[i]) - adj*np.abs(np.cos(Beta[i]))
            z[counter+j] = z1[i]+chi*(z2[i]-z1[i]) + np.sign(slope[i]) * adj*np.abs(np.sin(Beta[i]))
        counter += pps_weights[i]

    return x,y,z
#########################################
#get x,y,z distributions for sputtered W
#########################################

#confirm nP stays constant
for i in range(len(pps_weights)): pps_weights[i] = round(pps_weights[i])
nP_diff = int(nP-np.sum(pps_weights))

if nP_diff > 0:
    for i in range(abs(nP_diff)):
        rand_index = np.random.choice(len(pps_weights))
        pps_weights[rand_index] += 1
elif nP_diff < 0:
    while nP_diff<0:
        rand_index = np.random.choice(len(pps_weights))
        if pps_weights[rand_index] - 1 > 0:
            pps_weights[rand_index] -= 1
            nP_diff+=1

pps_weights = np.round(pps_weights)
int_weights = np.array(pps_weights,dtype='int')
nP_diff = nP-np.sum(int_weights)

print('total nP', nP)
print('pps over nP', np.sum(pps)/nP)
print('nP(r_mid):', int_weights)
print('nP_diff should be 0: ', nP_diff)


#define adjustment into the sheath because particles can't start exactly on the wall
adj = 1e-7

slope = np.zeros(len(r1))
Alpha = np.zeros(len(r1))
for i in range(len(r1)):
    if (r2[i]-r1[i])!=0:
        slope[i] = (z2[i]-z1[i])/(r2[i]-r1[i])
        Alpha[i] = np.abs(np.arctan((z2[i]-z1[i]) / (r2[i]-r1[i])))
    elif (r2[i]-r1[i])==0:
        slope[i] = 100
        Alpha[i] = 89.999*np.pi/180

Alpha = np.abs(Alpha) 
Beta = np.abs(np.pi/2 - Alpha)

#populate x,y,z with r_mid,0,z_mid
if configuration == 'random': 
    x,y,z = random(nP,int_weights,adj,slope,Beta, r1,z1,r2,z2)
elif configuration == 'uniform': 
    x,y,z = uniform(nP,int_weights,adj,slope,Beta, r1,z1,r2,z2)
elif configuration == 'midpoint': 
    x,y,z = midpoints(nP,int_weights, adj,slope,Beta, r1,z1,r2,z2)
else:
    print('(x,y,z) configuration not set')

#########################################
#get vx,vy,vz from IEADs
#########################################

# set up vx,vy,vz,E,theta,psi distributions
PartDist = ParticleDistribution.ParticleDistribution(nP, ListAttr=['vx','vy','vz'])

vx = np.zeros(1)
vy = np.zeros(1)
vz = np.zeros(1)

for i in range(len(int_weights)):

    weight = int(int_weights[i])

    if weight>0:
        m = np.sign(slope[i])
        #get IEADs for sputtered W
        E = PartDist.Generate(weight, 'Thomson')
        PolAng = PartDist.Generate(weight, 'SinCos', x=np.linspace(0,np.pi/2,10*weight))
        AziAng = PartDist.Generate(weight, 'Uniform', x=np.linspace(0,2*np.pi,10*weight))
        
        #convert IEADs to vx,vy,vz unit vectors in particle frame of ref
        PartDist.SetAttr('vx', np.multiply(np.cos(PolAng), np.cos(AziAng)))
        PartDist.SetAttr('vy', np.multiply(np.cos(PolAng), np.sin(AziAng)))
        PartDist.SetAttr('vz', m*np.sin(PolAng))
        
        vx_prime = PartDist.Particles['vx']
        vy_prime = PartDist.Particles['vy']
        vz_prime = PartDist.Particles['vz']

        #rotate vx,vy,vz from particle frame to lab frame
        PartDist.RotateAngle('v', -m*Alpha[i],0, Degree=False)
        
        vx_lab = PartDist.Particles['vx']
        vy_lab = PartDist.Particles['vy']
        vz_lab = PartDist.Particles['vz']

        #convert unit vectors to vx,vy,vz
        W_kg = 183.84 * 1.6605e-27 #mass of W in kg
        vtot = np.sqrt(2*E*1.6022e-19/W_kg) #convert eV to m/s
        
        vx = np.append(vx, vtot*vx_lab)
        vy = np.append(vy, vtot*vy_lab)
        vz = np.append(vz, vtot*vz_lab)

vx = np.delete(vx,0)
vy = np.delete(vy,0)
vz = np.delete(vz,0)

This block contains a call to the file and class ParticleDistribution which can be downloaded here: Particles.zip . The author of this code is Jerome Guterl from General Atomics.


Last, the particleSource.nc file must be created and its variables must be filled.

#########################################
#make NetCDF Particle Source file
#########################################

rootgrp = netCDF4.Dataset("particleSource.nc", "w", format="NETCDF4")
npp = rootgrp.createDimension("nP", nP)
xxx = rootgrp.createVariable("x","f8",("nP"))
yyy = rootgrp.createVariable("y","f8",("nP"))
zzz = rootgrp.createVariable("z","f8",("nP"))
vxx = rootgrp.createVariable("vx","f8",("nP"))
vyy = rootgrp.createVariable("vy","f8",("nP"))
vzz = rootgrp.createVariable("vz","f8",("nP"))
xxx[:] = x
yyy[:] = y
zzz[:] = z
vxx[:] = vx
vyy[:] = vy
vzz[:] = vz
rootgrp.close()