Numerical instability problem using azimuthal decomposition and prescribed fields

[CernuschiA]

Hi everyone,
I'm trying to run a 2D simulation with azimuthal mode decomposition, using some prescribed fields (a homogeneous magnetic field increasing over time and an electric field modelled as plane waves at 2.45 GHz. However, I'm clearly seeing sign of numerical instabilities:



Here are some solutions I tried to mitigate the problem:

• I tried to simplify the simulation as much as possible, removing any extra features except the ones strictly needed for the simulations, consisting in the plasma (e- and He2+) and the prescribed field (a homogeneous magnetic field increasing in time "Bl_m_mode_0" and an electric field modeled with plane waves "Er_mode_1" and "Et_mode_1"). As boundary condition I used [ ["silver-muller","silver-muller"],["reflective","buneman"]] for the EM (reflective in the radial direction along the cavity axis) and [["remove", "remove"],["reflective", "remove"],] for the particles. In this configuration I tried reducing as much as possible the cell size and the time-step: the cell size has been reduced up to 1/4 of the estimated Debye length (dz=l_D/4 and dr=l_D/2), while the time-step has been reduced to 1/2 of the estimated CFL condition for azimuthal mode decomposition. I'm having difficulties to reduce them further than this due to the limited computational resources on which I'm performing my test.
• I also tried adding further features which should mitigate numerical instabilities: a current filter, B-TIS3 interpolation method and PML boundary condition. Nevertheless, nothing is working.

I'm a little bit out of ideas, if anyone has some further suggestion or can try to take a look at the namelist it would be great.
In my simulations I actually implemented the prescribed field using cython, which helped me to speed up the use of prescribed field by quite a margin. I'm uploading a version of the namelist which uses direct python functions, since the cython implementation requires to know the specific version of cython you are using. If anyone need the faster version for some testing I would upload the complete namelist with cython functions.
test.zip

Thanks to everyone in advance!

[mccoys]

I did not follow the discussion, but is the problem present without prescribed fields? I wonder whether there is some particular numerical instability at r=0 with these nonphysical fields.

[CernuschiA]

I just made a test without the prescribed field and indeed I don't observe instabilities. However, the use of prescribed field is mandatory for the phenomenon I want to observe. I also made some previous tests with a more realistic magnetic field (generated with COMSOL and implemented through field maps), but I still observed numerical instabilities. I'm pretty sure the fields are suitable for the phenomenon I want to observe (I made some Monte Carlo simulation with single-electron approximation and indeed the electron behaviour reproduce the phenomenon I want to observe). Do you have any suggestion to debug the problem? Maybe put some boundary condition in the prescribed field functions to treat the case with r=0?

[mccoys]

Why not use external fields instead?

[CernuschiA]

The fields are both time-dependent (both magnetic and electric field) and space-dependent (electric field). In the full version of the simulation, the magnetic field will also be space dependent. From the explanation in the namelist, it would seem to me that the external fields are not suitable for my simulation. If I'm incorrect, please let me know.

[Z10Frank]

I am not sure to understand which fields you want to initialise in AM.
What are their analytical expressions in cylindrical coordinates and components?

Also, the fields evolved through the Maxwell solver explode from either one or more points of the grids where a "strange" field value appears, or a current density accumulates from the macro-particles. Do you see an anomalous current density J when the instability starts?

[CernuschiA]

The field I'm initializing in AM coordinates are:

• A constant magnetic field increasing in time according to a sinusoidal law (only axial component with perfect cylindrical symmetry, hence Bl_m_mode_0)
  B=Bmax*(1.-(Bmax-B0)/Bmax*sin(omega_pul*t+pi/2))

• An electric field implemented through plane waves (Ex=E0*cos(omega_hf*t-k*x) and Ey=E0*sin(omega_hf*t-k*x)) converted to cylindrical coordinates (dependent on cos(theta) and sin(theta), hence with components Er_mode_1 and Et_mode_1)
  Er=E0*cos(omega_hf*t-k*x)+E0*sin(omega_hf*t-k*x)*j
Et=-E0*sin(omega_hf*t-k*x)+E0*cos(omega_hf*t-k*x)*j
  One thing that I noted and that I forgot to mention is that at the beginning of the simulation (when the code is solving the poisson equation before starting the time evolution), the phase omega_hf*t-k*x of Er_mode_1 and Et_mode_1 seems to be non-finite and hence I have to add a line like
  if not np.isfinite(phase): phase = 0.
  to let the simulation work properly.

Regarding anomalies in the current density, I'm currently running a new simulation and I will let you know as soon as it finishes.

[Z10Frank]

Wait, if by x you mean the propagation direction as in the code convention (see figure in https://smileipic.github.io/Smilei/Understand/azimuthal_modes_decomposition.html), then your Er and Et are only proportional to Ey. The Ex would be the El component.

I don't know if this is the origin of the instability, but better to debug a namelist that reflects what you want to simulate. If the field configuration that is forced through the Prescribed field gives a strange behaviour (e.g. an accumulation of particles in a region), everything is possible as instability.

[CernuschiA]

I'm sorry for the confusion, I'm used to another nomenclature and I wrote the message a little bit too fast. With Ex and Ey I meant the transversal direction to the cavity axis (hence Ey and Ez with Smilei nomenclature), I'm pretty sure I used the correct convention in the namelist.
Indeed, the phenomenon I want to observe would lead to an accumulation of the electron population in a bunch, which can cause this strange behavior (from my diagnostics I haven't still clearly observed the phenomenon). However, it would seem to me that the problem I'm encountering it's a numerical instability, and not a physical one (maybe the two can be related, I'm quite new to PIC codes).
By the way, can it be that the source of instabilities is a low number of particles per cell? Currently, I'm using 32 for each species (2 in the radial and axial direction, 16 in the azimuthal one), since I have limited computational resources.
In any case I added some further diagnostics in the new simulation I'm running, as soon as I have some results I will try to observe the current density as you suggested.

[Z10Frank]

I obtain for mode 1

Er_1=Ey+1j*Ez
Et_1=-1j*Ey+Ez,

and Bl_0=B_x for mode 0, maybe we just have a different sign?

I wonder if we need to explicitly add also functions for the missing fields (a complex zero for all x,r,t), i.e. El_0,Er_0,Et_0, El_1, Br_0, Bt_0,Bl_1,Br_1,Bt_1, can you try?