# Copyright (c) 2020 ETH Zurich
"""
Module for the MonteCarlo class.
The MonteCarlo class implements all random-number based methods to simulate the motion
of electrons. The simulation time-step is determined with the null collision technique
(calculate_max_coll_freq creates a lookup table for the choice of the
trial collision frequency, determine_timestep calculates the time-step based on the
current maximum electron energy and acceleration). The collision processes are randomly
chosen based on the collision frequency of each process. The scattering angles are
randomly chosen with either an isotropic or an anisotropic model.
"""
# Import Packages
import numpy as np
import numpy.linalg as linalg
import scipy.constants as csts
import scipy.interpolate
from typing import Tuple
# Import modules
import pyMETHES.utils as utils
from pyMETHES.__about__ import __version__
from pyMETHES.config import Config
from pyMETHES.gas_mixture import GasMixture
np.seterr(all='raise')
[docs]class MonteCarlo:
"""
Class implementing all random-number based simulation methods
Attributes:
config (Config): configuration of the simulation
trial_coll_freq (float): trial collision frequency for the null collision
technique
max_coll_freq (interp1d): cumulative maximum of the collision frequency as a
function of the electron energy
max_coll_period (ndarray): array inversely proportional to the cumulative
maximum of the collision frequency
max_coll_period_squared (ndarray): array inversely proportional to the square
of the cumulative maximum of the collision frequency
collision_by_electron (ndarray): array of collision index for each electron
(starting at 0), an index equal to the number of cross sections indicates
a null collision
"""
version = f"pyMETHES version {__version__}\n"
[docs] def __init__(self, cfg: Config):
"""
Instantiates the MonteCarlo class.
Args:
cfg (Config): configuration of the simulation
"""
self.config = cfg
# trial collision frequency
self.trial_coll_freq = None
# data for the calculation of the trial collision frequency
self.max_coll_freq = None
self.max_coll_period = None
self.max_coll_period_squared = None
# vector of which collision which electron undergoes
self.collision_by_electron = None
[docs] def calculate_max_coll_freq(self, gas_mixture: GasMixture):
"""
Calculates the maximum collision frequency in the given gas mixture.
Args:
gas_mixture (GasMixture): gas mixture
"""
gas_density = self.config.gas_number_density
velocity = utils.velocity_from_energy(gas_mixture.energy_vector)
freq = gas_density * velocity * gas_mixture.total_cross_section
# replace zero value to avoid divide by zero:
freq[0] = freq[1]
# take cumulative maximum of the collision frequency (make it monotonous)
freq = np.maximum.accumulate(freq)
self.max_coll_freq = scipy.interpolate.interp1d(gas_mixture.energy_vector,
freq,
kind='linear')
self.max_coll_period = \
0.5 * csts.electron_mass / csts.elementary_charge / freq
self.max_coll_period_squared = \
0.5 * csts.electron_mass / csts.elementary_charge / freq ** 2
[docs] def determine_timestep(self, max_velocity: float, max_acceleration: float) -> float:
"""
Determine the duration of the next time-step in the simulation with the
null-collision technique.
Args:
max_velocity (float): current maximum electron velocity
max_acceleration (float): current maximum electron acceleration
Returns: time-step duration (s)
"""
# Draw random number
s = - np.log(np.random.random())
# Calculate energy after step (worst case: max_velocity + max_acceleration * dt)
max_energy = utils.energy_from_velocity(max_velocity)
freq = self.max_coll_freq.y
de = 2 * max_velocity * max_acceleration * s * self.max_coll_period \
+ max_acceleration ** 2 * s ** 2 * self.max_coll_period_squared
# Clip energy after step to maximum energy
e_end = np.clip(max_energy + de, a_min=None,
a_max=self.config.max_cross_section_energy)
# Find trial collision frequency that fulfills the condition
ok_frequencies = freq[self.max_coll_freq(e_end) <= freq]
# Take the first (= smallest) frequency that fulfills the condition
self.trial_coll_freq = 1.01 * ok_frequencies[0]
# Calculate the time-step
dt = s/self.trial_coll_freq
return dt
[docs] def determine_collisions(self, gas_mixture: GasMixture,
velocity_norm: np.ndarray,
energy: np.ndarray) -> None:
"""
Calculates the collision frequencies for all electrons with all cross sections,
and chooses a collision type via a random number.
Args:
gas_mixture (GasMixture): cross section data
velocity_norm (ndarray): norm of the velocity of each electron
energy (ndarray): energy of each electron
"""
# calculate the coll frequency for each cross section and each electron
collision_matrix = np.array([
self.config.gas_number_density * interp(energy) * velocity_norm
/ self.trial_coll_freq for interp in gas_mixture.cross_sections])
rand_vector = np.random.rand(1, energy.size)
rand_matrix = np.repeat(rand_vector, collision_matrix.shape[0], axis=0)
coll_matrix_cum = np.cumsum(collision_matrix, axis=0)
self.collision_by_electron = np.sum(rand_matrix > coll_matrix_cum, axis=0)
[docs] @staticmethod
def unit_scattered_velocity(energy: np.ndarray, velocity: np.ndarray,
iso: bool) -> Tuple[np.ndarray, np.ndarray]:
"""
Calculates the new direction of the velocity vector after scattering.
Args:
energy (ndarray): energy of each electron (eV)
velocity (ndarray): velocity of each electron in (x,y,z) directions before
the collision (m.s-1)
iso (bool): isotropic scattering (True) or anisotropic scattering (False)
Returns: normed velocities after collisions, cosine of polar scattering angle
"""
cos_chi, sin_chi, cos_phi, sin_phi = MonteCarlo.scattering_angles(energy, iso)
# Normalize the velocity
v_hat = velocity / linalg.norm(velocity, axis=0, keepdims=True)
# unit vector in x direction
e_x = np.zeros_like(velocity)
e_x[0, :] = 1
# calculate theta (angle between v_hat and e_x)
theta = np.arccos(v_hat[0, :])
sin_theta = np.sin(theta)
# calculate the direction of the new velocity, see Vahedi 1995 equation (11)
v_new_dir = \
v_hat * cos_chi \
+ np.cross(v_hat, e_x, axis=0) * sin_chi * sin_phi / sin_theta \
+ np.cross(v_hat, np.cross(e_x, v_hat, axis=0), axis=0) \
* sin_chi * cos_phi / sin_theta
return v_new_dir, cos_chi
[docs] @staticmethod
def scattering_angles(energy: np.ndarray, iso: bool
) -> (np.ndarray, np.ndarray, np.ndarray, np.ndarray):
"""
Generates values for the polar (chi) and azimuthal (phi)
isotropic or anisotropic scattering angles according to Vahedi (1995).
Args:
energy: array of electron energies
iso: isotropic scattering or not
Returns: 4 arrays cos(chi), sin(chi), cos(phi), sin(phi)
"""
# choose scattering angle phi
phi = 2 * np.pi * np.random.random(energy.size)
sin_phi = np.sin(phi)
cos_phi = np.cos(phi)
# choose scattering angle chi
if iso:
cos_chi = 1 - 2 * np.random.random(energy.size)
else:
cos_chi = (2 + energy
- 2 * (1 + energy) ** np.random.random(energy.size)) / energy
sin_chi = np.sqrt(1 - cos_chi ** 2)
return cos_chi, sin_chi, cos_phi, sin_phi