add new reactions between desired particles at capacitive discharge 2d Example

[mojgan-zare]

Hello,
I need to know how can I add new particles and new reactions witch are not in the excitation, ionization, and elastic collisions category to the 2d capacitive discharge example?
it is important for me.
Thank you.

[mojgan-zare]

Can I simulate negative H plasma ion source by warpx ?

[mojgan-zare]

Is there any way that I send my code ro you to see , why I faced error ?

[RemiLehe]

Thanks for your interest in WarpX.

I will answer your first question (how can I add new particles and new reactions). For the two other questions: could you open a separate Github issue (https://github.com/ECP-WarpX/WarpX/issues), so that we have one issue per independent question? (It makes it easier to track and answer.)

Regarding the new reaction: which type of reaction would you like to add?

[mojgan-zare]

Hi. Thanks for your answer.
some reactions like electron detachment process of H minus and mutual neutralization of H minus and charge exchange of H minus ?

[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.