search

ECP-WarpX / WarpX

WarpX is an advanced electromagnetic & electrostatic Particle-In-Cell code.
https://ecp-warpx.github.io
Other
309 stars 196 forks source link

add new reactions between desired particles #5478

Open mojgan-zare opened 2 hours ago

mojgan-zare commented 2 hours ago

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%.5et)" % (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 commented 2 hours ago

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 commented 2 hours ago

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.