Model of Inductive Heating in 1D PIC

[budjensen]

In a 2005 paper, Meige et. al. (https://doi.org/10.1063/1.1897390), introduced an approach to simulating inductive heating in 1D PIC codes. Would it be possible to add their method to WarpX?

To replicate this model, we would need an additional method of specifying a perpendicular external electric field--the user would input a time-dependent expression for a perpendicular ("inductive") current, $J_y(z,t)$. At each time step, a perpendicular electric field is applied according to the condition:

$$E_\text{y}(z_i,t) = E_\text{y}(z_\text{i},t-\Delta t) + \frac{\partial E_\text{y}(z_\text{i},t)}{\partial t}\Delta t$$

where

$$\frac{\partial E_\text{y}(z_\text{i},t)}{\partial t} = \frac{1}{\varepsilon_0} \left[ J_\text{y}(z_\text{i},t) - J_\text{cond}(z_\text{i},t) \right].$$

Would this be a reasonable addition?

[roelof-groenewald]

Hi @budjensen. I find the extensibility of WarpX through Python callbacks to be very useful for these kinds of situations where a rather particular problem is being addressed. To replicate the above model you could use the afterEsolve callback to update the Ey component of the field after the regular field update.
A rough sketch of how this would work could be:

from pywarpx import callbacks, fields

...

def update_Ey_field():
    Ey = fields.EyFPWrapper()
    J_cond = fields.JyWrapper()
    Jy = get_inductive_current()
    # update Ey field
    Ey[...] += 1/eps0 * (J_cond - Jy) * dt
...

# somewhere before calling simulation.step()
callbacks.installafterEsolve(update_Ey_field)
...

Let me know if you have any questions on how to use the callbacks or fields functionality. As a possible helping guide see https://warpx.readthedocs.io/en/latest/usage/examples.html#manipulating-fields-via-python.

[budjensen]

Hi @roelof-groenewald, thanks for the quick response. I'm working on implementing this and realizing that fields.Jx/y/zWrapper() is only deposited in electromagnetic simulations and fields.Ey/zFPWrapper() is not calculated in 1D simulations (and is reset to zero between each step).

I could implement this method exclusively in python by accessing particle data and calculating current with something like:

from pywarpx import callbacks
...
def ICP_heating():
    # Loop through particles and accumulate particle current to cells
    J_cond = get_cond_current()

    Jy = get_inductive_current()

    # update ICP field
    ICP_Ey += 1/eps0 * (J_cond - Jy) * dt

    # Push particles
   do_ICP_push(ICP_Ey)
...

# somewhere before calling simulation.step()
ICP_Ey = [0.0] * nz    # Initialize the ICP field
callbacks.installafterEsolve(ICP_heating)
...

My hesitance is that calculating the conduction current and particle push in python will be much slower than in the original C++ algorithm. Is there any places you see I could do this more efficiently within the current WarpX functionality? I'm thinking along the lines of something like a current corollary to the ParticleContainerWrapper.deposit_charge_density() function.

Thank you!

[roelof-groenewald]

Hey, sorry for the slow reply.
Having a deposit_current density() function analogous to deposit_charge_density() is something that's been on my wishlist for a while but I haven't gotten around to it.

I think making a prototype solver (including the current deposition) in Python would be a good way to confirm that things work as expected. Then afterwards functionality can be moved piece-by-piece to C++ to recover performance. Note that you wouldn't have to do the particle pushing in Python since you'll be able to write all field components during the Python field solve step. That should help to limit the performance hit.

You wouldn't have to do the particle pushing in Python since you'll be able to

[budjensen]

Hey @roelof-groenewald, thanks for the continued discussion.

I went ahead and implemented a version using python callbacks, but have ran into an amrex related issue (just submitted under #5194) while verifying the consistency of my implementation with literature. Once these are resolved I should be able to get a solid working version.

With what I've found, updating the Ex and Ey components in a 1D simulation does not lead to the particle velocities ux and uy being pushed during the Epush. Should a push happen along each of the three velocity coordinates, regardless of the simulation dimensionality?

[roelof-groenewald]

The particle velocities should definitely be updated in all three directions. A common issue that comes up when assigning new field values via a Python callback is that Python overwrites variables if you are not careful. Consider the code below to see what I mean:

from pywarpx import callbacks, fields

...

def update_Ey_field():
    J_cond = fields.JyWrapper()
    Jy = get_inductive_current()

    Ey = fields.EyFPWrapper()
    new_Ey_values = Ey[...] + 1/eps0 * (J_cond - Jy) * dt

    # correct approach: update Ey field by assigning new values to existing WarpX memory
    Ey[...] = new_Ey_values
    # incorrect approach: overwrite Ey variable without assigning new values to WarpX memory
    Ey = new_Ey_values
...

# somewhere before calling simulation.step()
callbacks.installafterEsolve(update_Ey_field)
...

Are you perhaps not writing the new field values properly into WarpX memory?