Possible energy leak with PEC boundaries

[snicks11]

Hi all, sorry to potentially re-open this issue. I'd first like to say that I indeed see that PEC boundaries are indeed working better. However, with further testing I am also seeing an energy leak into the simulation for both fields and particles. (My input deck is essentially the same as before, but for the purpose of reproduction I have attached it below.) This energy leak seems to be related to the formation of a sheath (Ez) at the ends of the domain. Of course, PEC boundaries allow a nonzero Ez, but with reflecting particles, this sheath is unexpected. If I use particles of equal mass, the sheath and energy leak are suppressed. Also using a large number of particles per cell seems to help. I see the same sheath formation as well in 2D.

Energy in transverse field components:

Energy in particles:

Field traces at beginning and end of domain:

Input deck: inputs_1d_picmi.zip

Originally posted by @snicks11 in https://github.com/ECP-WarpX/WarpX/issues/3664#issuecomment-1421350444

[ax3l]

cc @RevathiJambunathan @prkkumar can you potentially follow up on this? :)

[roelof-groenewald]

I've been wondering if this is related to an issue @clarkse noticed a couple of weeks ago (which I will describe here, but Eric please feel free to correct any errors). When using PEC boundaries, the charge density (and current density, I believe) is reduced in the cells close to the boundaries, presumably due to the particle shape extended the "particle volume" over the boundary. An example of this can be seen in the attached plot, Eric made.

I'm wondering if this could be causing an asymmetric reduction in density between electrons and ions (when they have unequal masses), which is what leads to the improper formation of a sheath as @snicks11 observed.

Asking around a bit, I understand that in some PIC codes during particle deposition an appropriate mirror charge is considered for particles that extend over the PEC boundary. The mirror charge consequently extends into the simulation domain, cancelling the density reduction.

[snicks11]

@RevathiJambunathan @prkkumar has there been any update on this issue?

[snicks11]

@RevathiJambunathan please use this updated input deck. Thank you! PICMI_inputs_1d.zip

[RevathiJambunathan]

Thank you @snicks11

[snicks11]

Energy conservation becomes worse when using linear particles instead of cubic:

[snicks11]

When testing two oppositely charged particles in a small domain, the particle motion seems fine. However, field energy is obviously very low in this case. This suggests to me that the issue is arising due to the fields when large numbers of particles are present.

[RevathiJambunathan]

Looks like PR #3711 needs to be merged first, from discussions at the WarpX meeting

[snicks11]

@RevathiJambunathan as an update on this issue, I see the same phenomenon in 2D RZ. I am using a uniform quiescent (no physical instabilities) plasma with magnetic field directed along z. For particles, I use reflecting boundaries on the radial edges and periodic boundaries on the axial edges. For fields, the outer radial edge has a PEC boundary while the inner boundary has "none". I see sheath formation, most severely for the inner radial edge.

Electric fields vs. radius near the midplane

2D perturbed Bz field

Evidently a sheath is forming at the radial edges, causing a radial electric field, which then causes an ExB drift that generates an axial perturbed magnetic field. (If only this could be harnessed in a real FRC.) Energy conservation is consequently extremely bad, mostly due to fields.

Perturbed energy breakdown (particles, fields, and total.)

Input deck: inputs_2D_RZ_picmi.zip

@dpgrote if you have any input into the best way to use boundaries in RZ, it would be appreciated.

(Apologies that the plots may not show up well on dark mode) Also @roelof-groenewald

[snicks11]

Update: this phenomenon still occurs in RZ with an absorbing particle boundary at r = 0.

[ax3l]

We revisited this this week and @jlvay recommended to try the spectral solver to overcome the current show-stopper (noise is too large on axis).

We can also check with @dpgrote and @RemiLehe (will be back in August from vacation).

[dpgrote]

I'm working on an implicit PIC solver for WarpX and to get exact energy conservation for this in RZ, I've had to make a few changes to the RZ solver, specifically fixing the fields on the axis. These changes carry over to the explicit PIC solver. It would be interesting if you could try running your case using my branch https://github.com/dpgrote/WarpX/tree/implicit_picard and see if it behaves any better.

[snicks11]

Thanks @ax3l and @dpgrote. I will try to run with the spectral solver and the mentioned implicit branch to see if these options help.

[snicks11]

Strategy 1: The spectral solver seems to be currently incompatible with PEC boundaries, and PML boundaries are (as far as I know) not currently supported in RZ. The chief problem is the outer radial boundary. Are there field boundary options available for the spectral solver in RZ other than periodic?

Strategy 2: @dpgrote I do not see a substantive difference in the numerical noise on the axis when compiling with the implicit_picard branch.

[snicks11]

@dpgrote do you happen to have and RZ example scripts made for your fork?

[dpgrote]

That branch doesn't make any other changes to the explicit version (beyond the handling of the boundaries). There wouldn't be any differences in the input file.

[roelof-groenewald]

@snicks11 Sorry to just add more requested tests without giving ideas of possible causes, but I was wondering this morning: do you know if this issue only shows up when you have a finite B-field applied? This might be why the issue is not widely seen - not many people are running WarpX with a strong background magnetic field (as far as I know). Also if the issue can be made much more noticeable by applying a large external B-field it might help in the debugging process.

[snicks11]

@roelof-groenewald the axis field noise appears regardless of whether I have an external magnetic field. However, I have not quantitatively checked whether it is worse with a magnetic field. (The perturbed fields look different because there are ExB forces that do not exist in the case of no external field.)

[snicks11]

@roelof-groenewald @dpgrote @ax3l Looking at the density, I am seeing an instability that seems to be driven by the strongly nonuniform pseudoparticles/volume density (all particles treated with unity weight) with respect to radius. Below are plots of the pseudoparticle density as a function of radius for the beginning and ending timesteps. Most of the transition in fact occurs in about 1/4 as many timesteps. The initial trend matches the expected 1/r dependence (linear on a log-log plot) but dramatically flattens.

The corresponding real density (particles treated with their real weights) is plotted below. The initial distribution essentially shows a uniform radial density, as is expected. During the development of the instability, strong density fluctuations occur at the axis before the distribution settles into that shown below.

Upon reflection, this phenomenon is caused by (or at least related to) the simple diffusion of hot pseudo-particles when they are initialized in a non-uniform distribution. This raises several questions, some of which are

1. Is RZ mostly intended to be used with cold (foil) plasmas?
2. Is it possible to change the pseudo-particles per cell as a function of radius

Input deck for this run: PICMI_inputs_2D_RZ.zip

[snicks11]

When macroparticles are loaded such that macroparticle density is constant (via a callback function), the field noise at the origin persists. The electric field profile near the midplane is shown below.

Below are the pseudo(macro)particle density and real particle density. Clearly there is some spurious activity at the axis related to the field noise.

Input deck: PICMI_inputs_2D_RZ.zip