Uniform plasma in external current and magnetic field

[cuphyalex]

Hi WarpX community,

I am trying to inject background ions using NUniformPerCell and add an external current using Jx_external_grid_function(x,y,z,t) and an external magnetic field in hybrid solver. However, the simulation doesn't run after exporting the first step of diagnostics.

3d_rotating_jet.txt

[RemiLehe]

Thanks for reaching out!

I think that the reason the run stalls right after the first diagnostic is that your external B-field and external current expressions are singular at the origin.

Look at the radial factor 1/sqrt(x**2+y**2) that appears in all four of them:

warpx.Bx_external_grid_function(x,y,z) =  y/(x**2+y**2)**(1/2)*B0/2
warpx.By_external_grid_function(x,y,z) = -x/(x**2+y**2)**(1/2)*B0/2
hybrid_pic_model.Jx_external_grid_function(x,y,z,t) = -J0*x/sqrt(x**2+y**2)
hybrid_pic_model.Jy_external_grid_function(x,y,z,t) = -J0*y/sqrt(x**2+y**2)

Your grid is collocated, n_cell = 8 8 8, and the domain is centered on 0 (prob_lo = -0.5*L, prob_hi = 0.5*L). That means there is a grid node sitting exactly at x = y = 0, where these expressions evaluate to 0/0 = NaN. WarpX writes the initial (step-0) diagnostic and then, as soon as the hybrid solver uses these fields, the NaN propagates through the Ohm's-law E-field and the B-field sub-stepping, so the simulation appears to "hang" / stop progressing.

A few ways to fix it on your side:

• Add a small regularization term under the square root so the denominator is never zero:

    warpx.Bx_external_grid_function(x,y,z) =  y/sqrt(x**2+y**2+eps**2)*B0/2
    warpx.By_external_grid_function(x,y,z) = -x/sqrt(x**2+y**2+eps**2)*B0/2
    hybrid_pic_model.Jx_external_grid_function(x,y,z,t) = -J0*x/sqrt(x**2+y**2+eps**2)
    hybrid_pic_model.Jy_external_grid_function(x,y,z,t) = -J0*y/sqrt(x**2+y**2+eps**2)

with e.g. `my_constants.eps = 1e-3*dx` (small compared to a cell, so it only affects the single on-axis point). At the origin the numerator is also 0, so the field smoothly goes to 0 there.

• Guard the on-axis point explicitly, e.g.

hybrid_pic_model.Jx_external_grid_function(x,y,z,t) = if(x**2+y**2>0, -J0*x/sqrt(x**2+y**2), 0)

Physically, these two solutions say "the azimuthal field/current is zero on the axis", which is the natural value for these ~1/r profiles.