Current Continuity in 1D PIC

[budjensen]

When I run the conditions for Turner's benchmark (https://doi.org/10.1063/1.4775084, Examples/Physics_applications/capacitive_discharge), I can match time averaged ion densities, ionization rate, and power to ions from the benchmark:

but am unable to match power to electrons:

My guess is that the root of this problem is that total current is not continuous across the device, as should be the case from the steady state continuity equations: $\nabla\cdot(J_\text{i} + J_\text{e}) = 0 $. When I looked into this and plotted time averaged ion and electron currents, I found:

Do you know what is leading to this non-continuous current in WarpX? I have ran Turner's benchmark on another code in the past and did not have issues with current continuity.

[budjensen]

Some possibly helpful information I've found since my initial post:

Turner's benchmark utilizes a time step such that $\Delta{t}/\omega_\text{pi} = 0.12$. This yields results as shown above. When I decrease the WarpX time step to where $\Delta{t}/\omega_\text{pi} = 0.03$, I get closer agreement to power put into electrons and a somewhat improved electron current. Any thoughts what this may indicate?

Thank you!

[budjensen]

Perhaps a bit clearer of a visualization, here are the time averaged fluxes:

[roelof-groenewald]

Hi @budjensen. Thanks for bringing this up!

I'm wondering if the time-centering of the currents might play a role here. How are you getting the individual J's? In WarpX at time $t_n$ the $x_n$ positions and $u_{n-1/2}$ velocities are known, so to get $J_{t=n}$ the velocities should be pushed forward by a half time step. An issue of this type would show smaller error with smaller $\Delta t$ as you observed and be more noticeable for the electrons since their velocities change more drastically over a single step.

[budjensen]

Good thinking, thanks @roelof-groenewald.

I've written a python callback which uses the particle velocities, weights, and positions to calculate the current. I install this callback in the beforestep position.

Just to ensure I'm implementing/tracking values correctly, at the beforestep location we know the following quantities: $E_{n}$, $x_{n}$, and $v_{n-1/2}$ with $\frac{t}{\Delta{t}}=n$. From there I can do $v_{n-1/2} \rightarrow {v}_{n}$ with $E_n$ to calculate $J_n$. Is this correct?

[roelof-groenewald]

Yes, that is correct. There is actually functionality in WarpX to sync the particle velocities and positions (see here). So if you wanted to just quickly check if this really is the problem, you could modify WarpXEvolve.cpp line 126 as follows:

Synchronize();
ExecutePythonCallback("beforestep");

This should not break anything else since after this there is a check if the positions and velocities are synchronized and if they are, the velocity is moved back.

[budjensen]

I made that exact edit and things are looking exactly as I'd expect (see below). Thank you!

I moved Synchronize() to callbacks.cpp line 34 like

#include "WarpX.H"

...

// Execute Python callbacks of the type given by the input string
void ExecutePythonCallback ( const std::string& name )
{
    if ( IsPythonCallbackInstalled(name) ) {
        WARPX_PROFILE("warpx_py_" + name);

        // Sychronize before the "beforestep" callbacks
        if (name == "beforestep") {
            WarpX::GetInstance().Synchronize();
        }

        try {
            warpx_callback_py_map[name]();

...

to minimize function calls. Would it be possible to open Synchronize() up to python?

[roelof-groenewald]

That's great! Thanks for following up about this.

Yes, we can expose Synchronize() to Python. Basically, it should work by adding the following snippet to line 256 of Python/WarpX.cpp (https://github.com/ECP-WarpX/WarpX/blob/development/Source/Python/WarpX.cpp):

        .def_static("synchronize",
            [] () { WarpX::Synchronize(); },
            "Synchronize particle velocities and positions."
        )

You'd then be able to execute it using

simulation.extension.warpx.synchronize()

Would you like to open a PR for that?

[budjensen]

Sure! I'm not directly familiar with pybind11, but after a few compilation errors I got this snippet to work (and run correctly when I used it in warpx via python):

        .def("synchronize",
            [] (WarpX& wx) { wx.Synchronize(); },
            "Synchronize particle velocities and positions."
        )

I'll go ahead and create the PR, and keep an eye out here in case this implementation should be modified. Thank you for the help!