"""
Deposition process code
The Process class implements the methods necessary to support the deposition process.
"""
# Default packages
import math
import os, sys
# Core packages
import warnings
from threading import Lock
import numpy as np
import numexpr_mod as ne
# Local packages
from febid.Structure import Structure
import febid.diffusion as diffusion
import febid.heat_transfer as heat_transfer
from febid.libraries.rolling.roll import surface_temp_av
from timeit import default_timer as df
# TODO: look into k-d trees
[docs]def restrict(func):
"""
Prevent simultaneous call of the decorated methods
"""
def inner(self):
flag = self.lock.acquire()
while not flag:
flag = self.lock.acquire()
return_vals = func(self)
flag_e = self.lock.release()
return return_vals
return inner
# Deprication note:
# At some point, due to efficiency advantages, the diffusion calculation approach switched from 'rolling' to 'stencil'.
# The rolling approach explicitly requires the array of ghost cells, while stencil does not, although still relying
# on this approach. Instead of the ghost cell array, it checks the same 'precursor' array, that it gets as a base argument,
# for zero cells.
# The ghost array is still kept and maintained throughout the simulation for conceptual clearness and visualisation
[docs]class Process():
"""
Class representing the core deposition process.
It contains all necessary arrays, variables, parameters and methods to construct a continuous deposition process.
"""
### A note_ to value correspondance:
# The main reaction equation operates in absolute values of precursor density per nanometer.
# The precursor array stores values as such.
# The deposited volume is though calculated and stored as a fraction of the cell's volume.
# Thus, the precursor density has to be multiplied by the cell's base area divided by the cell volume.
def __init__(self, structure:Structure, equation_values, timings, deposition_scaling=1, temp_tracking=True, name=None):
if not name:
self.name = str(np.random.randint(000000, 999999, 1)[0])
else:
self.name = name
# Declaring necessary properties
self.structure = None
self.cell_dimension = None
self.cell_V = None
# Main arrays
# Semi-surface cells are virtual cells that can hold precursor density value, but cannot create deposit.
# Their role is to serve as pipes on steps, where regular surface cells are not in contact.
# Therefore, these cells take part in diffusion process, but are not taken into account when calculating
# other terms in the FEBID equation or deposit increment.
self.deposit = None # contains values from 0 to 1 describing the filling state of a cell
self.precursor = None # contains values from 0 to 1 describing normalized precursor density in a cell
self.surface = None # a boolean array, surface cells are True
self.semi_surface = None # a boolean array, semi-surface cells are True
self.surface_n_neighbors = None # a boolean array, surface n-nearest neighbors are True
self.ghosts = None # a boolean array, ghost cells are True
self.beam_matrix = None # contains values of the SE surface flux
self.temp = None # contains temperatures of each cell
self.surface_temp = None
# Working arrays
self.__deposit_reduced_3d = None
self.__precursor_reduced_3d = None
self.__surface_reduced_3d = None
self.__semi_surface_reduced_3d = None
self.__surface_all_reduced_3d = None
self.__ghosts_reduced_3d = None
self.__temp_reduced_3d = None
self.__ghosts_reduced_2d = None
self.__beam_matrix_reduced_2d = None
self.D_temp = None
self.tau_temp = None
# Helpers
self._surface_all = None
self.__beam_matrix_surface = None
self.__beam_matrix_effective = None
self.__deposition_index = None
self.__surface_index = None
self.__semi_surface_index = None
self._solid_index = None
self.__surface_all_index = None
self.__tau_flat = None
# Monte Carlo simulation instance
self.sim = None
# Physical variables
self.F = 0 # precursor surface flux
self.n0 = 0 # maximum absolute precursor density
self.sigma = 0 # integral cross-section
self.tau = 0 # residence time
self.V = 0 # deposit volume of a dissociated precursor molecule
self.D = 0 # surface diffusion coefficient
self.cp = 0 # heat capacity
self.heat_cond = 0 # heat conductance coefficient
self.rho = 0 # density
self.room_temp = 294 # room temperature
# Timings
self.t_diffusion = 0
self.t_dissociation = 0
self.t_desorption = 0
self.dt = 0
self.t = 0
# Utility variables
self.__neibs_sides = None # a stencil array used when surface is updated
self.__neibs_edges = None # a similar stencil
self.irradiated_area_3D = np.s_[:, :, :]
self.irradiated_area_2D = np.s_[:, :, :]
self.__area_2D_no_sub = np.s_[:, :, :]
self.deposition_scaling = deposition_scaling # multiplier of the deposit increment; used to speed up the process
self.redraw = True # flag for external functions saying that surface has been updated
self.t_prev = 0
self.vol_prev = 0
self.growth_rate = 0
self.temperature_tracking = temp_tracking
self.request_temp_recalc = temp_tracking
self.temp_step = 10000 # amount of volume to be deposited before next temperature calculation
self.temp_step_cells = 0 # number of cells to be filled before next temperature calculation
self.temp_calc_count = 0 # counting number of times temperature has been calculated
# Statistics
self.substrate_height = 0 # Thickness of the substrate
self.n_substrate_cells = 0 # the number of the cells in the substrate
self.max_neib = 0 # the number of surface nearest neighbors that could be escaped by a SE
self.max_z = 0 # maximum height of the deposited structure, cells
self.max_z_prev = 0
self.filled_cells = 0 # current number of filled cells
self.n_filled_cells = []
self.growth_rate = []
self.dep_vol = 0 # deposited volume
self.max_T = 0
self.execution_speed = 0
self.profiler = None
self.stats_freq = 1e-2 # s
self.min_precursor_covearge = 0
self.lock = Lock()
# Initialization sequence
self.__set_structure(structure)
self.__set_constants(equation_values)
self.precursor[self.surface] = self.nr
self.__get_timings(timings)
self.__expressions()
self.__get_utils()
# Initialization methods
def __set_structure(self, structure: Structure):
self.structure = structure
self.deposit = self.structure.deposit
self.precursor = self.structure.precursor
self.surface = self.structure.surface_bool
self.semi_surface = self.structure.semi_surface_bool
self.surface_n_neighbors = self.structure.surface_neighbors_bool
self._surface_all = np.logical_or(self.surface, self.semi_surface)
self.ghosts = self.structure.ghosts_bool
self.beam_matrix = np.zeros_like(structure.deposit, dtype=np.int32)
self.temp = self.structure.temperature
self.surface_temp = np.zeros_like(self.temp)
self.D_temp = np.zeros_like(self.precursor)
self.tau_temp = np.zeros_like(self.precursor)
self.cell_dimension = self.structure.cell_dimension
self.cell_V = self.cell_dimension ** 3
self.__get_max_z()
self.substrate_height = structure.substrate_height
self.irradiated_area_2D = np.s_[self.structure.substrate_height - 1:self.max_z, :, :]
self.__area_2D_no_sub = np.s_[self.substrate_height + 1:self.max_z, :, :]
self.__update_views_2d()
self.__update_views_3d()
self.__generate_surface_index()
self.__generate_deposition_index()
self._get_solid_index()
self.n_substrate_cells = self.deposit[:structure.substrate_height].size
self.temp_step_cells = self.temp_step / self.cell_V
def __set_constants(self, params):
self.kb = 0.00008617
self.F = params['F']
self.n0 = params['n0']
self.V = params['V']
self.sigma = params['sigma']
self.tau = params['tau']
self.k0 = params.get('k0', 0)
self.Ea = params.get('Ea', 0)
self.D = params['D']
self.D0 = params.get('D0', 0)
self.Ed = params.get('Ed', 0)
self.cp = params['cp']
self.heat_cond = params['heat_cond']
self.rho = params['rho']
self.deposition_scaling = params['deposition_scaling']
if self.temperature_tracking and not all([self.k0, self.Ea, self.D0, self.Ed]):
warnings.warn('Some of the temperature dependent parameters were not found! \n '
'Switch to static temperature mode? y/n')
answer: str = input().lower()
if answer in ['y', 'n']:
if answer == 'y':
self.temperature_tracking = False
print('Switching to static temperature mode.')
return
if answer == 'n':
print('Terminating.')
sys.exit('Exiting due to insufficient parameters for temperature tracking')
self.__get_surface_temp()
self.residence_time()
self.diffusion_coefficient()
def __get_timings(self, timings):
# Stability time steps
self.t_diffusion = diffusion.get_diffusion_stability_time(self.D, self.cell_dimension)
self.t_dissociation = timings['t_flux']
self.t_desorption = timings['t_desorption']
# self.get_dt()
self.dt = min(self.t_desorption, self.t_diffusion, self.t_dissociation)
self.dt = self.dt - self.dt / 10
self.t = self.dt
def get_dt(self):
if self.temperature_tracking:
D_max = self.__D_temp_reduced_2d.max()
self.t_diffusion = diffusion.get_diffusion_stability_time(D_max, self.cell_dimension)
self.t_desorption = self.__tau_temp_reduced_2d[self.__surface_reduced_2d].min()
self.t_dissociation = 1 / self.beam_matrix.max() / self.sigma
self.dt = min(self.t_desorption, self.t_diffusion, self.t_dissociation)
self.dt = self.dt - self.dt / 10
def view_dt(self, units='µs'):
m = 1E6
if units not in ['s', 'ms', 'µs', 'ns']:
print(f'Unacceptable input for time units, use one of the following: s, ms, µs, ns.')
if units == 's':
m = 1
if units == 'ms':
m = 1E3
if units == 'µs':
m = 1E6
if units == 'ns':
m = 1E9
print(f'Current time step is {self.dt * m} {units}. \n'
f'Time step is evaluated as the shortest stability time \n'
f'of the following process divided by 5: \n'
f' Diffusion: \t Dissociation: \t Desorption: \n'
f' {self.t_diffusion * m} {units} \t {self.t_dissociation * m} {units} \t {self.t_desorption * m} {units}')
# Computational methods
[docs] def check_cells_filled(self):
"""
Check if any deposit cells are fully filled
:return: bool
"""
if self.__deposit_reduced_3d.max() >= 1:
return True
return False
[docs] def update_surface(self):
"""
Updates all data arrays after a cell is filled.
:return:
"""
# What here actually done is marking the filled cell as a solid and a ghost cell and then updating surface,
# semi-surface, ghosts and precursor to describe the surface geometry around the newly filled cell.
# The approach is cell-centric, which means all the surroundings are processed
nd = (self.__deposit_reduced_3d >= 1).nonzero()
new_deposits = [(nd[0][i], nd[1][i], nd[2][i]) for i in range(nd[0].shape[0])]
self.filled_cells += len(new_deposits)
for cell in new_deposits:
surplus_deposit = self.__deposit_reduced_3d[
cell] - 1 # saving deposit overfill to distribute among the neighbors later
surplus_precursor = self.__precursor_reduced_3d[cell]
self.__deposit_reduced_3d[cell] = -1 # a fully deposited cell is always a minus unity
self.__temp_reduced_3d[cell] = self.room_temp
self.__precursor_reduced_3d[cell] = 0
self.__ghosts_reduced_3d[cell] = True # deposited cell belongs to ghost shell
self.__surface_reduced_3d[cell] = False # rising the surface one cell up (new cell)
# self.redraw = True
# Instead of using classical conditions, boolean arrays are used to select elements
# First, a condition array is created, that picks only elements that satisfy conditions
# Then this array is used as index
neibs_sides, neibs_edges = self.__neibs_sides, self.__neibs_edges
# Creating a view with the 1st nearest neighbors to the deposited cell
z_min, z_max, y_min, y_max, x_min, x_max = 0, 0, 0, 0, 0, 0
# Taking into account cases when the cell is located at the edge:
# Small note_: views should be first decreased from the end ([:2])
# and then from the beginning. Otherwise, decreasing from the end will have no effect.
if cell[0] + 2 > self.__deposit_reduced_3d.shape[0]:
z_max = self.__deposit_reduced_3d.shape[0]
neibs_sides = neibs_sides[:2, :, :]
neibs_edges = neibs_edges[:2, :, :]
else:
z_max = cell[0] + 2
if cell[0] - 1 < 0:
z_min = 0
neibs_sides = neibs_sides[1:, :, :]
neibs_edges = neibs_edges[1:, :, :]
else:
z_min = cell[0] - 1
if cell[1] + 2 > self.__deposit_reduced_3d.shape[1]:
y_max = self.__deposit_reduced_3d.shape[1]
neibs_sides = neibs_sides[:, :2, :]
neibs_edges = neibs_edges[:, :2, :]
else:
y_max = cell[1] + 2
if cell[1] - 1 < 0:
y_min = 0
neibs_sides = neibs_sides[:, 1:, :]
neibs_edges = neibs_edges[:, 1:, :]
else:
y_min = cell[1] - 1
if cell[2] + 2 > self.__deposit_reduced_3d.shape[2]:
x_max = self.__deposit_reduced_3d.shape[2]
neibs_sides = neibs_sides[:, :, :2]
neibs_edges = neibs_edges[:, :, :2]
else:
x_max = cell[2] + 2
if cell[2] - 1 < 0:
x_min = 0
neibs_sides = neibs_sides[:, :, 1:]
neibs_edges = neibs_edges[:, :, 1:]
else:
x_min = cell[2] - 1
neighbors_1st = np.s_[z_min:z_max, y_min:y_max, x_min:x_max]
# Creating a view with the 2nd nearest neighbors to the deposited cell
if cell[0] - 2 < 0:
z_min = 0
else:
z_min = cell[0] - 2
if cell[0] + 3 > self.__deposit_reduced_3d.shape[0]:
z_max = self.__deposit_reduced_3d.shape[0]
else:
z_max = cell[0] + 3
if cell[1] - 2 < 0:
y_min = 0
else:
y_min = cell[1] - 2
if cell[1] + 3 > self.__deposit_reduced_3d.shape[1]:
y_max = self.__deposit_reduced_3d.shape[1]
else:
y_max = cell[1] + 3
if cell[2] - 2 < 0:
x_min = 0
else:
x_min = cell[2] - 2
if cell[2] + 3 > self.__deposit_reduced_3d.shape[2]:
x_max = self.__deposit_reduced_3d.shape[2]
else:
x_max = cell[2] + 3
neighbors_2nd = np.s_[z_min:z_max, y_min:y_max, x_min:x_max]
# Processing cell configuration according to the cell evolution rules
deposit_kern = self.__deposit_reduced_3d[neighbors_1st]
semi_s_kern = self.__semi_surface_reduced_3d[neighbors_1st]
surf_kern = self.__surface_reduced_3d[neighbors_1st]
temp_kern = self.__temp_reduced_3d[neighbors_1st]
# Creating condition array
condition = np.logical_and(deposit_kern == 0,
neibs_sides) # True for elements that are not deposited and are side neighbors
# Updating main arrays
semi_s_kern[condition] = False
surf_kern[condition] = True
condition = np.logical_and(np.logical_and(deposit_kern == 0, surf_kern == 0),
neibs_edges) # True for elements that are not deposited, not surface cells and are edge neighbors
semi_s_kern[condition] = True
ghosts_kern = self.__ghosts_reduced_3d[neighbors_1st]
ghosts_kern[...] = False
deposit_kern[surf_kern] += surplus_deposit / np.count_nonzero(surf_kern) # distributing among the neighbors
condition = (temp_kern > self.room_temp)
if np.any(condition):
self.__temp_reduced_3d[cell] = temp_kern[condition].sum() / np.count_nonzero(condition)
surf_kern = self.__surface_reduced_3d[neighbors_2nd]
semi_s_kern = self.__semi_surface_reduced_3d[neighbors_2nd]
ghosts_kern = self.__ghosts_reduced_3d[neighbors_2nd]
condition = np.logical_and(surf_kern == 0,
semi_s_kern == 0) # True for elements that are neither surface nor semi-surface cells
ghosts_kern[condition] = True
self.__deposit_reduced_3d[cell] = -1 # a fully deposited cell is always a minus unity
self.__precursor_reduced_3d[cell] = 0 # precursor density in the deposited cell is always 0
self.__ghosts_reduced_3d[cell] = True # deposited cell belongs to ghost shell
self.__surface_reduced_3d[cell] = False # deposited cell is no longer a surface cell
precursor_kern = self.__precursor_reduced_3d[neighbors_2nd]
condition = (semi_s_kern | surf_kern) & (precursor_kern < 1e-6)
precursor_kern[condition] = surplus_precursor
# precursor_kern[condition] += 0.000001 # only for plotting purpose (to pass vtk threshold filter)
# precursor_kern[condition] += 0.000001
else:
self.__get_max_z()
if self.temperature_tracking and self.max_z - self.substrate_height - 3 > 2:
self.request_temp_recalc = self.filled_cells > self.temp_calc_count * self.temp_step_cells
if self.request_temp_recalc:
self._get_solid_index()
self.__area_2D_no_sub = np.s_[self.substrate_height + 1:self.max_z, :, :]
self.irradiated_area_2D = np.s_[self.structure.substrate_height - 1:self.max_z, :,
:] # a volume encapsulating the whole surface
self.__update_views_2d()
# Updating local surface nearest neighbors
vert_slice = (self.irradiated_area_2D[1], self.irradiated_area_3D[1], self.irradiated_area_3D[2])
deposit_red_2d = self.deposit[vert_slice]
surface_red_2d = self.surface[vert_slice]
neighbors_2d = self.surface_n_neighbors[vert_slice]
cell = (cell[0] + self.irradiated_area_3D[0].start, cell[1], cell[2])
if cell[0] - 3 < 0:
z_min = 0
else:
z_min = cell[0] - 3
if cell[0] + 4 > deposit_red_2d.shape[0]:
z_max = deposit_red_2d.shape[0]
else:
z_max = cell[0] + 4
if cell[1] - 3 < 0:
y_min = 0
else:
y_min = cell[1] - 3
if cell[1] + 4 > deposit_red_2d.shape[1]:
y_max = deposit_red_2d.shape[1]
else:
y_max = cell[1] + 4
if cell[2] - 3 < 0:
x_min = 0
else:
x_min = cell[2] - 3
if cell[2] + 4 > deposit_red_2d.shape[2]:
x_max = deposit_red_2d.shape[2]
else:
x_max = cell[2] + 4
n_3d = np.s_[z_min:z_max, y_min:y_max, x_min:x_max]
self.structure.define_surface_neighbors(self.max_neib,
deposit_red_2d[n_3d],
surface_red_2d[n_3d],
neighbors_2d[n_3d])
if self.max_z + 5 > self.structure.shape[0]:
# Here the Structure is extended in height
# and all the references to the data arrays are renewed
with self.lock: # blocks run with Lock should exclude calls of decorated functions, otherwise the thread will hang
shape_old = self.structure.shape
self.structure.resize_structure(200)
self.structure.define_surface_neighbors(self.max_neib)
beam_matrix = self.beam_matrix # taking care of the beam_matrix, because __set_structure creates it empty
self.__set_structure(self.structure)
self.beam_matrix[:shape_old[0], :shape_old[1], :shape_old[2]] = beam_matrix
self.redraw = True
# Basically, none of the slices have to be updated, because they use indexes, not references.
return True
return False
[docs] def deposition(self):
"""
Calculate an increment of a deposited volume for all irradiated cells over a time step
:return:
"""
# Instead of processing cell by cell and on the whole surface, it is implemented to process only (effectively)
# irradiated area and array-wise(thanks to Numpy)
# np.float32 — ~1E-7, produced value — ~1E-10
const = self.sigma * self.V * self.dt * 1e6 * self.deposition_scaling / self.cell_V * self.cell_dimension ** 2 # multiplying by 1e6 to preserve accuracy
self.__deposit_reduced_3d[self.__deposition_index] += self.__precursor_reduced_3d[
self.__deposition_index] * self.__beam_matrix_effective * const / 1e6
[docs] def precursor_density(self):
"""
Calculate an increment of the precursor density for every surface cell
:return:
"""
# Here, surface_all represents surface+semi_surface cells.
# It is only used in diffusion calculation, because semi_surface cells cannot take part in deposition process
precursor = self.__precursor_reduced_2d
surface_all = self.__surface_all_reduced_2d
surface = self.__surface_reduced_2d
# diffusion_matrix = self._laplace_term_stencil(precursor, surface_all) # Diffusion term is calculated separately and added in the end
diffusion_matrix = self.__rk4_diffusion(precursor, surface_all)
precursor[surface] += self.__rk4(precursor[surface],
self.__beam_matrix_surface) # An increment is calculated through Runge-Kutta method without the diffusion term
precursor[surface_all] += diffusion_matrix[surface_all] # finally adding diffusion term
[docs] def equilibrate(self, eps=1e-4, max_it=10000):
"""
Bring precursor coverage to a steady state with a given accuracy
It is advised to run this method after updating the surface in order to determine a more accurate precursor
density value for newly acquired cells
:param eps: desired accuracy
"""
start = df()
for i in range(20):
# p_prev = self.__precursor_reduced_2d.copy()
self.precursor_density()
# norm = np.linalg.norm(self.__precursor_reduced_2d - p_prev)/ np.linalg.norm(self.__precursor_reduced_2d)
# if norm < eps:
# print(f'Took {i+1} iteration(s) to equilibrate, took {df() - start}')
# return 1
else:
# acc = str(norm)[:int(3-math.log10(eps))]
# warnings.warn(f'Failed to reach {eps} accuracy in {max_it} iterations in Process.equilibrate. Acheived accuracy: {acc} \n'
# f'Terminating loop.', RuntimeWarning)
print(f'Took {i + 1} iteration(s) to equilibrate, took {df() - start}')
def __rk4(self, precursor, beam_matrix):
"""
Calculates increment of precursor density by Runge-Kutta method
:param precursor: flat precursor array
:param beam_matrix: flat surface electron flux array
:return:
"""
k1 = self.__precursor_density_increment(precursor, beam_matrix,
self.dt) # this is actually an array of k1 coefficients
k2 = self.__precursor_density_increment(precursor, beam_matrix, self.dt / 2, k1 / 2)
k3 = self.__precursor_density_increment(precursor, beam_matrix, self.dt / 2, k2 / 2)
k4 = self.__precursor_density_increment(precursor, beam_matrix, self.dt, k3)
return ne.re_evaluate("rk4", casting='same_kind')
def __rk4_diffusion(self, grid, surface):
"""
Apply Runge-Kutta 4 method to the calculation of the diffusion term.
:param grid: precursor coverage array
:param surface: surface boolean array
:return:
"""
dt = self.dt
k1 = self._diffusion(grid, surface, dt, flat=False)
k1[surface] /= 2
k2 = self._diffusion(grid, surface, dt / 2, add=k1, flat=False)
k2[surface] /= 2
k3 = self._diffusion(grid, surface, dt / 2, add=k2, flat=False)
k4 = self._diffusion(grid, surface, dt, add=k3, flat=False)
return ne.re_evaluate("rk4", casting='same_kind')
def __precursor_density_increment(self, precursor, beam_matrix, dt, addon=0.0):
"""
Calculates increment of the precursor density without a diffusion term
:param precursor: flat precursor array
:param beam_matrix: flat surface electron flux array
:param dt: time step
:param addon: Runge Kutta term
:return:
"""
if self.temperature_tracking:
tau = self.__tau_flat
else:
tau = self.tau
return ne.re_evaluate('precursor_temp',
local_dict={'F': self.F, 'dt': dt, 'n0': self.n0,
'sigma': self.sigma, 'n': precursor + addon, 'tau': tau,
'se_flux': beam_matrix}, casting='same_kind')
def _diffusion(self, grid, surface, dt=0, add=0, flat=False):
"""
Calculates diffusion term of the reaction-diffusion equation for all surface cells.
:param grid: precursor coverage array
:param surface: boolean surface array
:param add: Runge-Kutta intermediate member
:return: flat ndarray
"""
if not dt:
dt = self.dt
if self.temperature_tracking:
D_param = self.__D_temp_reduced_2d
else:
D_param = self.D
return diffusion.diffusion_ftcs(grid, surface, D_param, dt, self.cell_dimension, self.__surface_all_index,
flat=flat, add=add)
[docs] def heat_transfer(self, heating):
"""
Define heating effect on the process
:param heating: volumetric heat sources distribution
:return:
"""
# Calculating temperature profile
if self.max_z - self.substrate_height - 3 > 2:
if self.request_temp_recalc:
self.temp_calc_count += 1
slice = self.__area_2D_no_sub # using only top layer of the substrate
if type(heating) is np.ndarray:
heat = heating[slice]
else:
heat = heating
# Running solution of the heat equation
start = df()
print(f'Current max. temperature: {self.max_T} K')
print(f'Total heating power: {heat.sum() / 1e6:.3f} W/nm/K/1e6')
heat_transfer.heat_transfer_steady_sor(self.temp[slice], self.heat_cond, self.cell_dimension, heat,
1e-8)
print(f'New max. temperature {self.temp.max():.3f} K')
print(f'Temperature recalculation took {df() - start:.4f} s')
self.temp[self.substrate_height] = self.room_temp
self.__get_surface_temp() # estimating surface temperature
self.diffusion_coefficient() # calculating surface diffusion coefficients
self.residence_time() # calculating residence times
[docs] def diffusion_coefficient(self):
"""
Calculate surface diffusion coefficient for every surface cell.
:return:
"""
self.D_temp[self._surface_all] = self.diffusion_coefficient_expression(self.surface_temp[self._surface_all])
[docs] def diffusion_coefficient_expression(self, temp=294):
"""
Calculate surface diffusion coefficient at a specified temperature.
:param temp: temperature, K
:return:
"""
return self.D0 * np.exp(-self.Ed / self.kb / temp)
[docs] def residence_time(self):
"""
Calculate residence time for every surface cell.
:return:
"""
self.__tau_flat = self.residence_time_expression(self.__surface_temp_reduced_2d[self.__surface_reduced_2d])
self.__tau_temp_reduced_2d[self.__surface_reduced_2d] = self.__tau_flat
[docs] def residence_time_expression(self, temp=294):
"""
Calculate residence time at the given temperature
:param temp: temperature, K
:return:
"""
return 1 / self.k0 * np.exp(self.Ea / self.kb / temp)
# Data maintenance methods
# These methods support an optimization path that provides up to 100x speed up
# 1. By selecting and processing chunks of arrays (views) that are effectively changing
# 2. By preparing indexes for arrays
[docs] def update_helper_arrays(self):
"""
Define new views to data arrays, create axillary indexes and flatten beam_matrix array
:return:
"""
# Basically, procedure provided by this method is optional and serves as an optimisation.
self.__get_irradiated_area()
self.__generate_deposition_index()
self.__generate_surface_index()
self.__flatten_beam_matrix_effective()
self.__flatten_beam_matrix_surface()
def __update_views_3d(self):
"""
Update view-arrays in accordance with the currently irradiated area
:return:
"""
# This is a part of a helper routine destined to reduce the number of cells processed on every iteration
# All the methods operating on main arrays (deposit, precursor, etc) use a corresponding view
# on the necessary array.
# '3D view' mentioned here can be referred to as a volume that encapsulates all cells that have been irradiated
self.__deposit_reduced_3d = self.deposit[self.irradiated_area_3D]
self.__precursor_reduced_3d = self.precursor[self.irradiated_area_3D]
self.__surface_reduced_3d = self.surface[self.irradiated_area_3D]
self.__semi_surface_reduced_3d = self.semi_surface[self.irradiated_area_3D]
self.__surface_all_reduced_3d = self._surface_all[self.irradiated_area_3D]
self.__ghosts_reduced_3d = self.ghosts[self.irradiated_area_3D]
self.__temp_reduced_3d = self.temp[self.irradiated_area_3D]
self.__beam_matrix_reduced_3d = self.beam_matrix[self.irradiated_area_3D]
def __update_views_2d(self):
"""
Update view-arrays in accordance with the current max height of the structure
:return:
"""
# This is a part of a helper routine destined to reduce the number of cells processed on every iteration
# All the methods operating on main arrays (deposit, precursor, etc) use a corresponding view
# on the necessary array.
# '2D view' taken here can be referred to as a volume that encapsulates
# the whole surface of the deposited structure. This means it takes a view only along the z-axis.
self.__deposit_reduced_2d = self.deposit[self.irradiated_area_2D]
self.__precursor_reduced_2d = self.precursor[self.irradiated_area_2D]
self.__surface_reduced_2d = self.surface[self.irradiated_area_2D]
self.__semi_surface_reduced_2d = self.semi_surface[self.irradiated_area_2D]
self.__surface_all_reduced_2d = self._surface_all[self.irradiated_area_2D]
self.__ghosts_reduced_2d = self.ghosts[self.irradiated_area_2D]
self.__beam_matrix_reduced_2d = self.beam_matrix[self.irradiated_area_2D]
self.__temp_reduced_2d = self.temp[self.irradiated_area_2D]
self.__surface_temp_reduced_2d = self.surface_temp[self.irradiated_area_2D]
self.__D_temp_reduced_2d = self.D_temp[self.irradiated_area_2D]
self.__tau_temp_reduced_2d = self.tau_temp[self.irradiated_area_2D]
def __get_irradiated_area(self):
"""
Get boundaries of the irradiated area in XY plane
:return:
"""
indices = np.nonzero(self.__beam_matrix_reduced_2d)
y_start, y_end, x_start, x_end = indices[1].min(), indices[1].max() + 1, indices[2].min(), indices[2].max() + 1
self.irradiated_area_3D = np.s_[self.structure.substrate_height:self.max_z, y_start:y_end,
x_start:x_end] # a slice of the currently irradiated area
self.__update_views_3d()
def __generate_deposition_index(self):
"""
Generate a tuple of indices of the cells that are irradiated for faster indexing in 'deposition' method
:return:
"""
self.__deposition_index = self.__beam_matrix_reduced_3d.nonzero()
def __generate_surface_index(self):
"""
Generate a tuple of indices for faster indexing in 'laplace_term' method
:return:
"""
self.__surface_all_reduced_2d[:, :, :] = np.logical_or(self.__surface_reduced_2d,
self.__semi_surface_reduced_2d)
index = self.__surface_all_reduced_2d.nonzero()
self.__surface_all_index = (np.intc(index[0]), np.intc(index[1]), np.intc(index[2]))
index = self.__surface_reduced_2d.nonzero()
self.__surface_index = (np.intc(index[0]), np.intc(index[1]), np.intc(index[2]))
index = self.__semi_surface_reduced_2d.nonzero()
self.__semi_surface_index = (np.intc(index[0]), np.intc(index[1]), np.intc(index[2]))
def __flatten_beam_matrix_effective(self):
"""
Extract a flattened array of non-zero elements from beam_matrix array
:return:
"""
self.__beam_matrix_effective = self.__beam_matrix_reduced_3d[self.__deposition_index]
def __flatten_beam_matrix_surface(self):
"""
Extract a flattened array of beam_matrix array
:return:
"""
self.__beam_matrix_surface = self.__beam_matrix_reduced_2d[self.__surface_reduced_2d]
def _get_solid_index(self):
index = self.deposit[self.__area_2D_no_sub]
index = (index < 0).nonzero()
self._solid_index = (np.intc(index[0]), np.intc(index[1]), np.intc(index[2]))
return self._solid_index
@restrict
def __get_max_z(self):
"""
Get z position of the highest not empty cell in the structure
:return:
"""
self.max_z = self.deposit.nonzero()[0].max() + 3
def __get_surface_temp(self):
self.__surface_temp_reduced_2d[...] = 0
surface_temp_av(self.__surface_temp_reduced_2d, self.__temp_reduced_2d, *self.__surface_index)
surface_temp_av(self.__surface_temp_reduced_2d, self.__surface_temp_reduced_2d, *self.__semi_surface_index)
self.max_T = self.max_temperature
# Misc
def __get_utils(self):
# Kernels for choosing cells
self.__neibs_sides = np.array([[[0, 0, 0], # chooses side neighbors
[0, 1, 0],
[0, 0, 0]],
[[0, 1, 0],
[1, 0, 1],
[0, 1, 0]],
[[0, 0, 0],
[0, 1, 0],
[0, 0, 0]]])
self.__neibs_edges = np.array([[[0, 1, 0], # chooses edge neighbors
[1, 0, 1],
[0, 1, 0]],
[[1, 0, 1],
[0, 0, 0],
[1, 0, 1]],
[[0, 1, 0],
[1, 0, 1],
[0, 1, 0]]])
def __expressions(self):
"""
Prepare math expressions for faster calculations
:return:
"""
# Precompiled expressions for numexpr_mod.reevaluate function
# Creating dummy variables of necessary types
k1, k2, k3, k4, F, n, n0, tau, sigma, dt = np.arange(10, dtype=np.float64)
se_flux = np.arange(1, dtype=np.int64)
ne.cache_expression("(k1+k4)/6 +(k2+k3)/3", 'rk4')
ne.cache_expression("(F * (1 - n / n0) - n / tau - n * sigma * se_flux) * dt", 'precursor')
ne.cache_expression("F * dt * (1 - n / n0) - n * dt / tau - n * sigma * se_flux * dt", 'precursor_temp')
# Properties
@property
def kd(self):
if self.temperature_tracking:
tau = self.residence_time_expression(self.room_temp)
else:
tau = self.tau
return self.F / self.n0 + 1 / tau + self.sigma * self.beam_matrix.max()
@property
def kr(self):
if self.temperature_tracking:
tau = self.residence_time_expression(self.room_temp)
else:
tau = self.tau
return self.F / self.n0 + 1 / tau
@property
def nr(self):
"""
Calculate replenished precursor coverage
:return:
"""
return self.F / self.kr
@property
def nd(self):
"""
Calculate depleted precursor coverage
:return:
"""
return self.F / self.kd
@property
@restrict
def max_temperature(self):
"""
Get the highest current temperature of the structure.
:return:
"""
return self.__surface_temp_reduced_2d.max()
@property
@restrict
def deposited_vol(self):
"""
Get total deposited volume.
:return:
"""
return (self.filled_cells + self.__deposit_reduced_2d[self.__surface_reduced_2d].sum()) * self.cell_V
@property
@restrict
def precursor_min(self):
"""
Get the lowest precursor density at the surface.
:return:
"""
return self.__precursor_reduced_2d[self.__surface_reduced_2d].min()
if __name__ == '__main__':
print("Current script does not have an entry point.....")
input('Press Enter to exit.')