RustBCA Coupling
Problem Description
RustBCA is a plasma-surface interactions code which hPIC2 optionally couples to through its boundary conditions. It models the sputtering, reflection, and implantation processes resulting from incident ion impacts.
In this problem, we consider a plasma made of the following species:
Two Boris-Buneman uniform beam ion species,
"H"and"B", with equal temperatures \(T_i\)One Boltzmann electron population
"e-"with temperature \(T_e\)
Both boundaries are hydrogen-boron RustBCA walls. Incident particle interactions are simulated in RustBCA, where sputtered, reflected, and deposited particles are subsequently added or removed from the hPIC2 particle list. The two ion populations interact electrostatically with each other and with the electrons, in addition to through electron impact ionization events.
Theory
The binary-collision approximation assumes that when charged particles impact a material surface, they undergo a cascade of elastic binary collisions with the surface atoms while experiencing inelastic electronic stopping. Each cascade continues until the particle has either been reflected out of the simulation domain or become implanted in the surface by dropping below a defined cutoff energy. Additionally, collisions can dislodge surface atoms which can sputter out from the surface.
Simulation Setup
One-dimensional domain
Uniform mesh
Hydrogen-boron RustBCA wall boundaries
Plasma density \(n_0 = 10^{18} \; [m^{-3}]\)
Number of computational particles = 100,000 hydrogen particles/population, 10,000 boron particles/population
Plasma temperature: \(T_i = T_e = 10 \; [eV]\)
Input File (TOML)
[input_mode]
hpic_mode = "pic"
simulation_tag = "rustbca_coupling"
units = "si"
[mesh]
type = "uniform"
[mesh.type_specification]
x1_points = [0.0, 0.01]
x1_elem_size = 1.0e-4 # 100 elements
[time]
num_time_steps = 1000
dt = 2.0e-9
[species]
[species."e-"]
type = "boltzmann"
[species."e-".type_params]
temperature = 116045.0
[species."H"]
mass = 1.67262192369e-27
type = "boris_buneman"
[species."H".type_params]
atomic_number = 1
initial_condition = "uniform_beam"
[species."H".type_params.initial_condition_params]
num_particles = 100000
charge_states = [ { charge_number = 1, density = 1.0e18 } ]
flow_velocity_1 = 0.0
flow_velocity_2 = 0.0
flow_velocity_3 = 0.0
temperature = 116045.0
[[species."H".type_params.boundary_conditions]]
boundary = "west"
type = "rust_bca"
[species."H".type_params.boundary_conditions.type_params]
incident_species_cutoff_energy = 1.60217662e-19
incident_species_surface_binding_energy = 1.60217662e-19
[[species."H".type_params.boundary_conditions.type_params.targets]]
species = "B"
wall_density = 1.309e29
cutoff_energy = 8.01088e-19
surface_binding_energy = 9.22854e-19
bulk_binding_energy = 0.0
[[species."H".type_params.boundary_conditions.type_params.targets]]
species = "H"
wall_density = 1.309e28
cutoff_energy = 1.52207e-19
surface_binding_energy = 2.40326e-19
bulk_binding_energy = 0.0
[[species."H".type_params.boundary_conditions]]
boundary = "east"
type = "rust_bca"
[species."H".type_params.boundary_conditions.type_params]
incident_species_cutoff_energy = 1.60217662e-19
incident_species_surface_binding_energy = 1.60217662e-19
[[species."H".type_params.boundary_conditions.type_params.targets]]
species = "B"
wall_density = 1.309e29
cutoff_energy = 8.01088e-19
surface_binding_energy = 9.22854e-19
bulk_binding_energy = 0.0
[[species."H".type_params.boundary_conditions.type_params.targets]]
species = "H"
wall_density = 1.309e28
cutoff_energy = 1.52207e-19
surface_binding_energy = 2.40326e-19
bulk_binding_energy = 0.0
[[species."H".type_params.volumetric_sources]]
type = "minimum_mass"
temperature = 116045.0
[species."B"]
mass = 1.79521e-26
type = "boris_buneman"
[species."B".type_params]
atomic_number = 5
initial_condition = "uniform_beam"
[species."B".type_params.initial_condition_params]
num_particles = 10000
charge_states = [ { charge_number = 1, density = 1.0e12 },
{ charge_number = 2, density = 2.0e13 },
{ charge_number = 3, density = 4.0e13 },
{ charge_number = 4, density = 9.0e13 },
{ charge_number = 5, density = 2.0e13 } ]
flow_velocity_1 = 0.0
flow_velocity_2 = 0.0
flow_velocity_3 = 0.0
temperature = 116045.0
[[species."B".type_params.boundary_conditions]]
boundary = "west"
type = "rust_bca"
[species."B".type_params.boundary_conditions.type_params]
incident_species_cutoff_energy = 8.01088e-19
incident_species_surface_binding_energy = 9.22854e-19
[[species."B".type_params.boundary_conditions.type_params.targets]]
species = "B"
wall_density = 1.309e29
cutoff_energy = 8.01088e-19
surface_binding_energy = 9.22854e-19
bulk_binding_energy = 0.0
[[species."B".type_params.boundary_conditions.type_params.targets]]
species = "H"
wall_density = 1.309e28
cutoff_energy = 1.52207e-19
surface_binding_energy = 2.40326e-19
bulk_binding_energy = 0.0
[[species."B".type_params.boundary_conditions]]
boundary = "east"
type = "rust_bca"
[species."B".type_params.boundary_conditions.type_params]
incident_species_cutoff_energy = 8.01088e-19
incident_species_surface_binding_energy = 9.22854e-19
[[species."B".type_params.boundary_conditions.type_params.targets]]
species = "B"
wall_density = 1.309e29
cutoff_energy = 8.01088e-19
surface_binding_energy = 9.22854e-19
bulk_binding_energy = 0.0
[[species."B".type_params.boundary_conditions.type_params.targets]]
species = "H"
wall_density = 1.309e28
cutoff_energy = 1.52207e-19
surface_binding_energy = 2.40326e-19
bulk_binding_energy = 0.0
[[species."B".type_params.volumetric_sources]]
type = "minimum_mass"
temperature = 116045.0
[interactions.electron_impact_ionization."B"]
electron_species = "e-"
A = [2.83984e-55, 2.32824e-55, 1.93601e-55, 2.04331e-55, 1.02679e-55] # [1.1063, 0.9070, 0.7542, 0.7960, 0.4000] x 10^-13 eV^2 cm^2
B = [[-2.74512e-55, -2.25637e-56], # [-1.0694, -0.0879] x 10^-13 eV^2 cm^2
[-1.22444e-55, 5.05693e-56], # [-0.4770, 0.1970] x 10^-13 eV^2 cm^2
[-4.85157e-57, -7.60285e-55, 1.9299e-54, -2.19299e-54, 7.98533e-55], # [-0.0189, -2.9618, 7.5182, -8.5431, 3.1108] x 10^-13 eV^2 cm^2
[-1.28451e-55, 2.26817e-55], # [-0.5004, 0.8836] x 10^-13 eV^2 cm^2
[]]
I = [1.32981e-18, 4.0294742e-18, 6.07706e-18, 4.15557e-17, 5.45093e-17] # [8.30, 25.15, 37.93, 259.37, 340.22] eV
[interactions.electron_impact_ionization."H"]
electron_species = "e-"
A = [4.73606e-56] # 0.1845e-13 eV^2 cm^2
B = [[-4.77456e-57, 3.15994e-56, -4.87981e-56, 2.44555e-55]] # [-0.0186, 0.1231, -0.1901, 0.9527] x 10^-13 eV^2 cm^2
I = [2.17896e-18] # 13.6 eV
[magnetic_field]
type = "uniform"
[magnetic_field.type_params]
b1 = 0.0
b2 = 0.0
b3 = 0.0
[electric_potential]
poisson_solver = "tridiag"
[[electric_potential.boundary_conditions]]
boundary = "west"
type = "dirichlet"
function = "constant"
[electric_potential.boundary_conditions.function_params]
value = 0.0
[[electric_potential.boundary_conditions]]
boundary = "east"
type = "dirichlet"
function = "constant"
[electric_potential.boundary_conditions.function_params]
value = 0.0
[output_diagnostics]
output_dir = "rustbca_coupling_example_out"
[output_diagnostics.logging]
log_level = "info"
timing_log_enabled = false
[output_diagnostics.particle_output]
stride = 50
source_counters = true
species = ["B", "H"]
[output_diagnostics.field_output]
stride = 50
The mesh used in this example is type = "uniform".
The domain size \(L\) was set to \(0.01 \; [m]\).
\(dx\) is set such that there are 100 total elements,
and is represented in the toml file as x1_elem_size = 1.0e-4.
Time step is set to \(2e-9 \; [s]\), for 1000 total steps.
Three species are specified under [species]: two ion populations,
"H" and "B"; with an electron background "e+". "H" ions are set to a charge state of +1 with a density of
\(1e18 \; [m^{-3}]\) and 100,000 computational particles. "B" ions are set to charge states +1 through +5 with
densities batween \(1e12 \; [m^{-3}]\) and \(9e13 \; [m^{-3}]\), split between 10,000 computational particles.
These parameters are listed in the "charge_states" and "num_particles" blocks, respectively, under each [species]
table.
mass is the species mass.
type is the way which the species is represented or simulated.
Two RustBCA wall boundaries are specified under the [[boundary_conditions]] subtables of each [species] table.
Each boundary consists of both "H" and "B", and requires the cutoff and surface binding energies of the
incident species as well as the density, cutoff energy, surface binding energy, and bulk binding energy of each wall species.
The [interactions.electron_impact_ionization] tables define coefficients for electron impact ionizations between each
"boris_buneman" ion population and the "boltzmann" electron population.
Running
To run the simulation, use the following command:
$ mpiexec -np 1 ./hpic2 -i rustbca_coupling.toml
The simulation will run for \(2e-6 \; [s]\) over 1000 time steps,
and all outputs will be tagged with the rustbca_coupling prefix.
The standard output will show the simulation progress:
$ mpiexec -np 1 ~/hpic2_dev/omp_opt_build/hpic2 -i rustbca_coupling.toml
2023-05-18 14:23:45.282732 UTC-05:00 info [main] [mpi rank 0] starting simulation.
2023-05-18 14:23:45.310319 UTC-05:00 info [ElectricFieldTridiag] [mpi rank 0] saving potential
2023-05-18 14:23:45.310503 UTC-05:00 info [ElectricFieldTridiag] [mpi rank 0] saving E field
2023-05-18 14:23:45.310627 UTC-05:00 info [ElectricFieldTridiag] [mpi rank 0] saving B field
2023-05-18 14:23:45.310728 UTC-05:00 info [ElectricFieldTridiag] [mpi rank 0] saving electrostatic energy to rustbca_coupling_example_out/rustbca_coupling_ENERGY_ELECTROSTATIC.dat
2023-05-18 14:23:45.314217 UTC-05:00 info [State] [mpi rank 0] At timestep 1
2023-05-18 14:23:45.433586 UTC-05:00 info [State] [mpi rank 0] At timestep 2
2023-05-18 14:23:45.555824 UTC-05:00 info [State] [mpi rank 0] At timestep 3
2023-05-18 14:23:45.669392 UTC-05:00 info [State] [mpi rank 0] At timestep 4
2023-05-18 14:23:45.800298 UTC-05:00 info [State] [mpi rank 0] At timestep 5
2023-05-18 14:23:45.922554 UTC-05:00 info [State] [mpi rank 0] At timestep 6
2023-05-18 14:23:46.054694 UTC-05:00 info [State] [mpi rank 0] At timestep 7
2023-05-18 14:23:46.188037 UTC-05:00 info [State] [mpi rank 0] At timestep 8
2023-05-18 14:23:46.321542 UTC-05:00 info [State] [mpi rank 0] At timestep 9
2023-05-18 14:23:46.457601 UTC-05:00 info [State] [mpi rank 0] At timestep 10
...
Plotting the results
The simulation results are stored in the rustbca_coupling_example_out folder.
The output files are split between a particle file for each ion species, a general fields file, and a
total electrostatic energy file. Their structures and
naming conventions are described here.
You can visualize the results in multiple ways, for example
using Python and the h5py package (follow instructions, etc.).
import h5py
import matplotlib.pyplot as plt
import numpy as np
output_prefix = 'rustbca_coupling_example_out/rustbca_coupling_'
# Reading field data
fields_file = h5py.File(output_prefix+'fields.hdf','r')
points = fields_file["VTKHDF/Points"][:,0]
nsteps = fields_file["VTKHDF/Steps"].attrs["NSteps"]
Efield = fields_file["VTKHDF/PointData/field_E"]
Ex = np.reshape(Efield[:,0], (nsteps,-1))
phi = fields_file["VTKHDF/PointData/field_phi"]
phi = np.reshape(phi,(nsteps,-1))
# Reading particle data
H_particle_file = h5py.File(output_prefix+'H.h5part', 'r')
B_particle_file = h5py.File(output_prefix+'B.h5part', 'r')
Hx = H_particle_file["Step#20/x"]
Bx = B_particle_file["Step#20/x"]
# Plot the final electric field
fig, ax = plt.subplots()
ax.plot(points, Ex[-1,:], label = 'E_x')
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('E')
fig.savefig('rustbca_coupling_E.png')
# Plot the final electric potential
fig, ax = plt.subplots()
ax.plot(points, phi[-1,:], label = 'phi')
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('phi')
fig.savefig('rustbca_coupling_phi.png')
# Plot the final particle positions
fig, ax = plt.subplots()
ax.scatter(Hx, np.zeros_like(Hx), label = 'H')
ax.scatter(Bx, np.zeros_like(Bx), label = 'B')
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.savefig('rustbca_coupling_particles.png')
fields_file.close()
H_particle_file.close()
B_particle_file.close()