Question regarding de-excitation process

[zhongerting]

Hello WarpX Community,
I'm using the two_product_reaction method in the new WarpX DSMC module to model multi-step ionization (ground state -> excited state -> [higher excited state] -> ion).

After checking the documentation, I haven't found a process for the de-excitation from a higher excited state back to the ground state or a lower excited state.

How is this physical phenomenon handled in the current implementation?

Thanks for your insights!

[zhongerting]

In other words, does WarpX support setting an attribute called "lifetime" (so to speak) for excited state particles, combined with random number calculations for radiative de-excitation within given time steps?

[RemiLehe]

Hi @zhongerting
This is a great suggestion! We currently don't have a way to handle de-excitation in WarpX, but it would probably be useful to a number of users, and have a broad range of application.
If you need this urgently, as a temporary solution, there is always the possibilty to use callback functions to try to implement this yourself in Python: https://warpx.readthedocs.io/en/latest/usage/workflows/python_extend.html

[zhongerting]

Hi @RemiLehe
Thanks for your guidence and support. I see the ParticleContainerWrapper in the High-Level Particle Wrapper chapter. However, I can't find the member function for removing particles.
So is it necessary to start building from the particles removal function? I am still confused. And are there any similar examples that can be referred to?

[roelof-groenewald]

Hi @zhongerting. To implement a de-excitation function in Python, I would think you would want to move particles from an "excited_species" container to a "ground_state" container. Here is an example Python function that I've used in the past to move particles between species containers:

    def _move_ions(self):
        """Function to move particles from thermal ion to fast ion containers."""

        if cp is not None:
            xp = cp
        else:
            xp = np

        if not hasattr(self, 'pc'):
            # get thermal ion container
            multi_pc = twxrun.sim_ext.warpx.multi_particle_container()
            self.pc = multi_pc.get_particle_container_from_name(self.ions.name)

            self.x_idx = self.pc.get_real_comp_index('x')
            self.y_idx = self.pc.get_real_comp_index('y')
            self.z_idx = self.pc.get_real_comp_index('z')
            self.ux_idx = self.pc.get_real_comp_index('ux')
            self.uy_idx = self.pc.get_real_comp_index('uy')
            self.uz_idx = self.pc.get_real_comp_index('uz')
            self.w_idx = self.pc.get_real_comp_index('w')

        # get the local number of ions
        nps = int(self.pc.total_number_of_particles(True, True) * 0.01)

        # allocate arrays for new fast ions
        x = np.zeros(nps)
        y = np.zeros(nps)
        z = np.zeros(nps)
        ux = np.zeros(nps)
        uy = np.zeros(nps)
        uz = np.zeros(nps)
        w = np.zeros(nps)
        start_idx = 0
        end_idx = 0

        for lvl in range(self.pc.finest_level + 1):
            # get every local chunk of particles
            for pti in self.pc.iterator(self.pc, level=lvl):
                # attributes in SoA format
                soa = pti.soa()

                # get x and z coordinate
                px = xp.array(soa.get_real_data(self.x_idx), copy=False)
                py = xp.array(soa.get_real_data(self.y_idx), copy=False)
                pz = xp.array(soa.get_real_data(self.z_idx), copy=False)

                # get velocities
                pux = xp.array(soa.get_real_data(self.ux_idx), copy=False)
                puy = xp.array(soa.get_real_data(self.uy_idx), copy=False)
                puz = xp.array(soa.get_real_data(self.uz_idx), copy=False)

                pw = xp.array(soa.get_real_data(self.w_idx), copy=False)

                idx = [insert some logic to choose which particles to move]

                end_idx = start_idx + len(idx)
                if end_idx > len(x):
                    print("Overflow on particles to move", flush=True)
                    end_idx = len(x)
                    idx = idx[:end_idx - start_idx]

                if xp == np:
                    x[start_idx:end_idx] = px[idx]
                    y[start_idx:end_idx] = py[idx]
                    z[start_idx:end_idx] = pz[idx]
                    ux[start_idx:end_idx] = pux[idx]
                    uy[start_idx:end_idx] = puy[idx]
                    uz[start_idx:end_idx] = puz[idx]
                    w[start_idx:end_idx] = pw[idx]
                else:
                    x[start_idx:end_idx] = px[idx].get()
                    y[start_idx:end_idx] = py[idx].get()
                    z[start_idx:end_idx] = pz[idx].get()
                    ux[start_idx:end_idx] = pux[idx].get()
                    uy[start_idx:end_idx] = puy[idx].get()
                    uz[start_idx:end_idx] = puz[idx].get()
                    w[start_idx:end_idx] = pw[idx].get()

                start_idx = end_idx

                # delete particles that will be moved
                cpuid = xp.array(soa.get_idcpu_data(), copy=False)
                cpuid[idx] = 19999999

        # trim the new particle arrays
        x  = x[:end_idx]
        y  = y[:end_idx]
        z  = z[:end_idx]
        ux = ux[:end_idx]
        uy = uy[:end_idx]
        uz = uz[:end_idx]
        w  = w[:end_idx]

        # create fast ion particles from the removed thermal ions
        pc = particle_containers.ParticleContainerWrapper(
            self.fast_ions.name
        )
        pc.add_particles(x=x, y=y, z=z, ux=ux, uy=uy, uz=uz, w=w)  

Hopefully you'll be able to work from this example for your purposes.