add new reactions between desired particles

[mojgan-zare]

Hi,
I will provide you with the code I'm using, which is similar to the capacitive discharge example, with changes to particles, reactions, and cross-section files, and the error I faced when running this code. I hope you can guide me through debugging the code. I think the problem could be related to energy or my desire for cross-section files, which I'm using. Or maybe something else should be done on the reference code. Could you please help me? I'm interested in using the warpx for this simulation. code: #!/usr/bin/env python3

import numpy as np
from scipy.sparse import csc_matrix
from scipy.sparse import linalg as sla
from pywarpx import callbacks, fields, picmi

constants = picmi.constants

##########################

Physics Parameters
##########################

D_CA = 0.067 # m
N_INERT = 9.64e20 # m^-3
T_INERT = 300.0 # K
FREQ = 13.56e6 # Hz
VOLTAGE = 450.0
M_HMINUS = 1.00784 * constants.m_p # kg, H- mass in kg
PLASMA_DENSITY = 2.56e14 # m^-3
T_ELEC = 30000.0 # K
DT = 1.0 / (400 * FREQ)

##########################

Numerics Parameters
##########################

max_steps = 50
diagnostic_intervals = "::50"

nx, ny = 128, 8
xmin, ymin = 0.0, 0.0
xmax, ymax = D_CA, D_CA / nx * ny
number_per_cell_each_dim = [32, 16]

##########################

Cross-Section Directory
##########################

cross_sec_dir = "warpx-data/MCC_cross_sections/H/"

files = {
"elastic": "e_H2_ELASTIC_modified.dat",
"excitation": "e_H2_EXCITATION_modified.dat",
"ionization": "e_H2_IONIZATION_modified.dat",
"mutual_neutralization": "MN_modified.dat",
"associative_detachment": "AD_modified.dat",
"non_associative_detachment": "NAD_modified.dat",
"charge_exchange": "CX_modified.dat",
"electron_detachment": "ED_modified.dat",
"elastic_collision": "eEC_modified.dat",
}

#############################

Specialized Poisson Solver
#############################

class PoissonSolverPseudo1D(picmi.ElectrostaticSolver):
def init(self, grid, **kwargs):
super(PoissonSolverPseudo1D, self).init(
grid=grid,
method=kwargs.pop("method", "Multigrid"),
required_precision=1,
**kwargs,
)
self.rho_wrapper = None
self.phi_wrapper = None

def solver_initialize_inputs(self):
self.right_voltage = self.grid.potential_xmax
self.grid.potential_xmin = None
self.grid.potential_xmax = None
self.grid.potential_ymin = None
self.grid.potential_ymax = None
self.grid.potential_zmin = None
self.grid.potential_zmax = None

super(PoissonSolverPseudo1D, self).solver_initialize_inputs()

self.nx = self.grid.number_of_cells[0]
self.nz = self.grid.number_of_cells[1]
self.dx = (self.grid.upper_bound[0] - self.grid.lower_bound[0]) / self.nx
self.dz = (self.grid.upper_bound[1] - self.grid.lower_bound[1]) / self.nz

if not np.isclose(self.dx, self.dz):
    raise RuntimeError("Direct solver requires dx = dz.")

self.nxguardrho = 2
self.nzguardrho = 2
self.nxguardphi = 1
self.nzguardphi = 1

self.phi = np.zeros(
    (self.nx + 1 + 2 * self.nxguardphi, self.nz + 1 + 2 * self.nzguardphi)
)

self.decompose_matrix()
callbacks.installpoissonsolver(self._run_solve)

def decompose_matrix(self):
self.nxsolve = self.nx + 1
self.nzsolve = self.nz + 3

A = np.zeros((self.nzsolve * self.nxsolve, self.nzsolve * self.nxsolve))
kk = 0
for ii in range(self.nxsolve):
    for jj in range(self.nzsolve):
        temp = np.zeros((self.nxsolve, self.nzsolve))
        if ii == 0 or ii == self.nxsolve - 1:
            temp[ii, jj] = 1.0
        elif ii == 1:
            temp[ii, jj] = -2.0
            temp[ii - 1, jj] = 1.0
            temp[ii + 1, jj] = 1.0
        elif ii == self.nxsolve - 2:
            temp[ii, jj] = -2.0
            temp[ii + 1, jj] = 1.0
            temp[ii - 1, jj] = 1.0
        elif jj == 0:
            temp[ii, jj] = 1.0
            temp[ii, -3] = -1.0
        elif jj == self.nzsolve - 1:
            temp[ii, jj] = 1.0
            temp[ii, 2] = -1.0
        else:
            temp[ii, jj] = -4.0
            temp[ii, jj + 1] = 1.0
            temp[ii, jj - 1] = 1.0
            temp[ii - 1, jj] = 1.0
            temp[ii + 1, jj] = 1.0

        A[kk] = temp.flatten()
        kk += 1

A = csc_matrix(A, dtype=np.float32)
self.lu = sla.splu(A)

def _run_solve(self):
if self.rho_wrapper is None:
self.rho_wrapper = fields.RhoFPWrapper(0, True)
self.rho_data = self.rho_wrapper[Ellipsis]

self.solve()

if self.phi_wrapper is None:
    self.phi_wrapper = fields.PhiFPWrapper(0, True)
self.phi_wrapper[Ellipsis] = self.phi

def solve(self):
right_voltage_expr = str(self.right_voltage) if isinstance(self.right_voltage, str) else "0.0"
right_voltage = eval(
right_voltage_expr,
{"t": sim.extension.warpx.gett_new(0), "sin": np.sin, "pi": np.pi},
)
left_voltage = 0.0

rho = (
    -self.rho_data[self.nxguardrho : -self.nxguardrho, self.nzguardrho : -self.nzguardrho]
    / constants.ep0
)

nx, nz = np.shape(rho)
source = np.zeros((nx, nz + 2), dtype=np.float32)
source[:, 1:-1] = rho * self.dx**2

source[0] = left_voltage
source[-1] = right_voltage

b = source.flatten()

flat_phi = self.lu.solve(b)
self.phi[self.nxguardphi : -self.nxguardphi] = flat_phi.reshape(np.shape(source))

self.phi[: self.nxguardphi] = left_voltage
self.phi[-self.nxguardphi :] = right_voltage

self.phi[:, : self.nzguardphi] = 0
self.phi[:, -self.nzguardphi :] = 0

##########################

Physics Components
##########################

v_rms_elec = np.sqrt(constants.kb * T_ELEC / constants.m_e)
v_rms_hminus = np.sqrt(constants.kb * T_INERT / M_HMINUS)

Uniform plasma distributions
uniform_plasma_elec = picmi.UniformDistribution(
density=PLASMA_DENSITY,
upper_bound=[None] * 3,
rms_velocity=[v_rms_elec] * 3,
directed_velocity=[0.0] * 3,
)

uniform_plasma_hminus = picmi.UniformDistribution(
density=PLASMA_DENSITY,
upper_bound=[None] * 3,
rms_velocity=[v_rms_hminus] * 3,
directed_velocity=[0.0] * 3,
)

Species definitions
electrons = picmi.Species(
particle_type="electron", name="electrons", initial_distribution=uniform_plasma_elec
)

h_minus = picmi.Species(
name="h_minus",
charge=-constants.q_e,
mass=M_HMINUS,
initial_distribution=uniform_plasma_hminus,
)

MCC collisions
mcc_electrons = picmi.MCCCollisions(
name="coll_elec",
species=electrons,
background_density=N_INERT,
background_temperature=T_INERT,
background_mass=h_minus.mass,
scattering_processes={
"elastic": {"cross_section": cross_sec_dir + "e_H2_ELASTIC.dat", "energy": 0.10},
"excitation1": {"cross_section": cross_sec_dir + files["excitation"], "energy": 0.04},
"ionization": {"cross_section": cross_sec_dir + files["ionization"], "energy": 15.40, "species": h_minus},
},
)

mcc_h_minus = picmi.MCCCollisions(
name="coll_hminus",
species=h_minus,
background_density=N_INERT,
background_temperature=T_INERT,
scattering_processes={
"mutual_neutralization": {"cross_section": cross_sec_dir + files["mutual_neutralization"]},
"associative_detachment": {"cross_section": cross_sec_dir + files["associative_detachment"]},
"non_associative_detachment": {"cross_section": cross_sec_dir + files["non_associative_detachment"]},
"charge_exchange": {"cross_section": cross_sec_dir + files["charge_exchange"]},
"electron_detachment": {"cross_section": cross_sec_dir + files["electron_detachment"]},
"elastic_collision": {"cross_section": cross_sec_dir + files["elastic_collision"]},
},
)

Magnetic field
uniform_magnetic_field = picmi.AnalyticInitialField(
Bx_expression="0.0",
By_expression="0.0",
Bz_expression="0.1", # Tesla
)

##########################

Numerics Components
##########################

grid = picmi.Cartesian2DGrid(
number_of_cells=[nx, ny],
warpx_max_grid_size=128,
lower_bound=[xmin, ymin],
upper_bound=[xmax, ymax],
bc_xmin="dirichlet",
bc_xmax="dirichlet",
bc_ymin="dirichlet",
bc_ymax="dirichlet",
warpx_potential_hi_x="%.1fsin(2pi*%.5e*t)" % (VOLTAGE, FREQ),
lower_boundary_conditions_particles=["absorbing", "absorbing"],
upper_boundary_conditions_particles=["absorbing", "absorbing"],
)

solver = PoissonSolverPseudo1D(grid=grid)

##########################

Diagnostics
##########################

particle_diag = picmi.ParticleDiagnostic(
name="diag1", period=diagnostic_intervals,
)

field_diag = picmi.FieldDiagnostic(
name="diag1", grid=grid, period=diagnostic_intervals, data_list=["rho_electrons", "rho_h_minus"],
)

##########################

Simulation Setup
##########################

sim = picmi.Simulation(
solver=solver,
time_step_size=DT,
max_steps=max_steps,
warpx_collisions=[mcc_electrons, mcc_h_minus],
)

sim.add_species(
electrons,
layout=picmi.GriddedLayout(
n_macroparticle_per_cell=number_per_cell_each_dim, grid=grid
),
)

sim.add_species(
h_minus,
layout=picmi.GriddedLayout(
n_macroparticle_per_cell=number_per_cell_each_dim, grid=grid
),
)

sim.add_applied_field(uniform_magnetic_field)
sim.add_diagnostic(particle_diag)
sim.add_diagnostic(field_diag)

##########################

Simulation Run
##########################

sim.step(max_steps)

Confirm that the external solver was run
assert hasattr(solver, "phi")

[mojgan-zare]

and i will provide you the error:
Initializing AMReX (24.10)...
OMP initialized with 32 OMP threads
AMReX (24.10) initialized
0::Assertion `(std::abs(energies[i] - energies[i-1] - dE) < dE / 100.0)' failed, file "/home/conda/feedstock_root/build_artifacts/warpx_1727944891349/work/Source/Particles/Collision/ScatteringProcess.cpp", line 122, Msg:

ERROR : Energy grid not evenly spaced.
!!!
SIGABRT
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
addr2line: DWARF error: can't find .debug_ranges section.
See Backtrace.0.0 file for details

[mojgan-zare]

Based on the error, I don't understand which of the energy spaces must be uniform. The energies in the first column of my cross-section files have uniform space! Or maybe it's a misunderstanding of the energies used in the MCC part of the code. Could you please clarify it for me? I hope you answer me as soon as possible. Thank you.