Non-Physical temporal magnetic jump on Smilei

[Gaetan-Gauth]

To resume our problem, we simulate 2D bump-on-tail configuration (i.e. an electron Maxwellian beam at temperature Te, at speed v=vb ~0.02 c and at density nb in [0.005, 0.0005]n0, a bulk of Maxwellian electron at temperature Te, at speed =0 and at density ne=n0-nb and ion neutralized background at temperature Ti=Te/10, at spped =0 and at density n0. We used Smilei 4.3 version. On a simulation, we observe non-physical temporal magnetic jump on Bz field on scalar and field diagnostics. We have obtain these figures using happi (+ python load text for scalar diagnostics) .

1. 2D - 512x512 cells, dt = 0.006 and nb = 0.005 nb :

• Figures with scalar diagnostic : We have got magnetic field energy for Bx, By, Bz save every 10 times step during 50 000 dt thus 5 000 saved : with a zoom during first time step :

2. 2D - 1024x1024 cells, dt = 0.01 and nb = 0.0005 nb :

• Figure with scalar diagnostic : We have got magnetic field energy for Bx, By, Bz save every 10 times step during 60 000 dt thus 6 000 saved. A magnetic jump appears on Bz at t = 120, 180, 400 :

• Figure with field diagnostic : This Bz jump appears on fields diagnostic save every 20 time step during 60 000 dt thus 6 000 saved. We have recalculated magnetic Bz energy using this diagnostic (sum on x and y of |Bz(x,y,t)|^2 in Smilei units) :
  This jump does not appear on Bx and By. Due to the difference between beam's density (10 times lower in second case), this instability grows up 10 times later. This two simulations use vectorization (=on), load balancing (every 150, cell load = 1, particle_load =0.1), without external field, interpolation order 2 (we have tested with order 4 without magnetic jump).

[mccoys]

Is it possible for you to obtain precisely the timesteps where this occured ?

[Gaetan-Gauth]

yes, this is occurred at 12580, 17780, 18140, 38830 timesteps.

[MickaelGrech]

@Gaetan-Gauth (1) could you confirm that these timesteps do no correspond to checkpoints? (2) could you add your simulation input file (change the .py extension by a .txt one)?

[Gaetan-Gauth]

@MickaelGrech (1) These timesteps are not correspond to restart point and checkpoints. (2) Find attach my input files input_Smilei2DGG.txt

Thanks !

[Z10Frank]

@Gaetan-Gauth To have some clues about where this anomalous field could come from, please could you show us some snapshots of Bx in the window, in time steps when the anomalous values start to grow? Also, saturating the colormap could show if these values come e.g. from the borders and start diffusing in the window.

Do not hesitate to navigate the time steps using the slide() function of the DiagFields and DiagProbes to spot the generation of anomalous points in the grid.

[Gaetan-Gauth]

@Z10Frank - This is some snapshots of Bz(t) around magnetic jump at t = 120 separated of dt * every saved field (look at top right) :

[Gaetan-Gauth]

Here are some new snapshots of Bz around the second and third magnetic jump. We observe localized instabilities:

[Z10Frank]

Thank you!

We will investigate, your plots surely are giving us some information.

Meanwhile, are you able to craft a smaller version of your simulation (less grid points, less particles, ...) where you can observe the same phenomenon?

[mccoys]

Let me comment on something unrelated to the bug. In your input you have

    if (x<=NumberOfCell_inX):

This compares a distance to a number of cells, which is incorrect.

Currently I am trying to find a smaller case to reproduce the bug.

[iltommi]

@Gaetan-Gauth sorry, I deleted the last comment by mistake, could you please re-post it?

[Gaetan-Gauth]

We found a solution ; when we reduce the patch size (16x16 to 8x8 with same simulation case i.e. 1800 part./cell and 2D 1024x1024), the non-physical magnetic jump disappears (see figures).

• 16x16 case

• 8x8 case

Thank a lot !

[mccoys]

@Gaetan-Gauth Could you please check whether the jump occurs on Ex, Ey or Ez ?

[Gaetan-Gauth]

@mccoys I can not see jumps in Ex field. Find attached Ex snapshots around first and second magnetic jump. But, the electric field amplitude is 100 times bigger than magnetic field amplitude.

• Energy (using scalar diagnostic)

• Ex field (using fields diagnostic)

[mccoys]

Do you have the scalars Uelm_Ex, Uelm_Ey and Uelm_Ez ?

[Gaetan-Gauth]

@mccoys Yes, I have just add U_elm_Ej (j =x,y,z) in previous message.

[mccoys]

@jderouillat, I have been told by @Gaetan-Gauth that the bug disappears when vectorisation is turned off. Is this something you can investigate ?

[beck-llr]

In order for us to investigate we would need to know which compiler you used and which compilation flag on which system ? Is there any chance we can access the data of these runs ? If it is on Irene for instance ?

[Gaetan-Gauth]

I work on Occigen with Smilei4.4 (+ this phenomenon appears with Smilei4.3 too) installed using make machine=occigen. My module list is 1) c/intel/.17.0 5) idb/17.0 9) hdf5/1.8.14 2) c++/intel/.17.0 6) intel/17.0 10) qt/4.8.6 3) fortran/intel/.17.0 7) intelmpi/5.1.3.258 11) python/2.7.13 4) mkl/17.0 8) szip/2.1 Where can I drop the interesting 2D fields/scalar files (About 10 Mo for scalars data and about 100 Go for fields) ?

[jderouillat]

Hi Gaetan, To be able to understand the origin of your problem, I'm trying to reproduce it on a downsized test case extracted from the input file provided, and unfortunately on a supercomputer which is not Occigen.
To be sure that we are not looking for a chimera, could you tell me if the bug is reproducible on the initial configuration (with vecto, patch size : 16 x 16) ? And the absence too (either without vecto or patch size : 8 x 8) ? It's not clear in your description Regards

[jderouillat]

Could you try your simulation commenting the line 755 of src/Projector/Projector2D2OrderV.cppand lines 766 to 777 of the same file to force to use the routine currents and not currentsAndDensity which is suspect as shown below ?

Rho (and per species Rho) will then not be available in diagnostics but the simulation will still be ok.

    // If no field diagnostics this timestep, then the projection is done directly on the total arrays
//    if( !diag_flag ) {
        if( !is_spectral ) {
            double *b_Jx =  &( *EMfields->Jx_ )( 0 );
            double *b_Jy =  &( *EMfields->Jy_ )( 0 );
            double *b_Jz =  &( *EMfields->Jz_ )( 0 );
            currents( b_Jx, b_Jy, b_Jz, particles,  istart, iend, invgf, iold, &( *delta )[0], ipart_ref );
        } else {
            ERROR( "TO DO with rho" );
        }

        // Otherwise, the projection may apply to the species-specific arrays
//    } else {
//        double *b_Jx  = EMfields->Jx_s [ispec] ? &( *EMfields->Jx_s [ispec] )( 0 ) : &( *EMfields->Jx_ )( 0 ) ;
//        double *b_Jy  = EMfields->Jy_s [ispec] ? &( *EMfields->Jy_s [ispec] )( 0 ) : &( *EMfields->Jy_ )( 0 ) ;
//        double *b_Jz  = EMfields->Jz_s [ispec] ? &( *EMfields->Jz_s [ispec] )( 0 ) : &( *EMfields->Jz_ )( 0 ) ;
//        double *b_rho = EMfields->rho_s[ispec] ? &( *EMfields->rho_s[ispec] )( 0 ) : &( *EMfields->rho_ )( 0 ) ;
//        for( int ipart=istart ; ipart<iend; ipart++ ) {
//            //Do not use cells sorting for now : f(ipart) for now, f(istart) laterfor now,
//            //(*iold)[ipart       ] = round( particles.position(0, ipart)* dx_inv_ - dt*particles.momentum(0, ipart)*(*invgf)[ipart] * dx_inv_ ) - i_domain_begin ;
//            //(*iold)[ipart+nparts] = round( particles.position(1, ipart)* dy_inv_ - dt*particles.momentum(1, ipart)*(*invgf)[ipart] * dy_inv_ ) - j_domain_begin ;
//            currentsAndDensity( b_Jx, b_Jy, b_Jz, b_rho, particles,  ipart, ( *invgf )[ipart-ipart_ref], iold, &( *delta )[ipart-ipart_ref], invgf->size() );
//        }
//    }

[Gaetan-Gauth]

Hi Julien, This bug is reproducible (I got it 3 times. For example when I try to remove load_balancing with initial conditions, I obtain again this bug) and the absence too. To sum up :

• Size_patch = 16x16, 1024x1024, 1800 part/cell, vecto and load_balancing : JUMP
• Size_patch = 16x16, 1024x1024, 1800 part/cell, vecto and NO load_balancing : JUMP
• Size_patch = 16x16, 1024x1024, 10 part/cell, vecto and load_balancing : NO JUMP
• Size_patch = 8x8, 1024x1024, 1800 part/cell, vecto and load_balancing : JUMP (later in time)
• Size_patch = 16x16, 1024x1024, 1800 part/cell, NO vecto and load_balancing : NO JUMP I hope to be clearer now. Regards,

[jderouillat]

Thank you for these details, the intersection of requirements (vectorization x patch size x number of particles per cell, especially the patch size) to produce the bug is very surprising !

Concerning the additional test that I asked you, the vectorized currentsAndDensity is a cheap adaptation of the routine used without vectorization (I didn't think that the vectorization was used in 2D ...) which was done just to validate the main vectorized algorithm.
Replacing it by currents make some noise disappear on the reproducer that I extracted from your data (I didn't reproduce your jump but a very light noise in vectorized mode). So I would be sure that it's a good path to solve your problem.
It's certainly expensive in terms of simulation time but very useful to the understanding.

[Gaetan-Gauth]

@jderouillat I have corrected my previous message, the case 8x8, 1024x1024, vecto and load balancing. Magnetic jump appears later in time (we have seen that now because tests did not go beyond 180 w_p t). The requested test is in queue now. Thanks !

[mccoys]

@Gaetan-Gauth Any news ?

[Gaetan-Gauth]

@mccoys -- sorry about late reply. When I have comment lines like jderouillat asked, I have results but nothing really append (see below magnetic figure and distribution function). But no magnetic jump ! For example, magnetic energy values are too low.

[mccoys]

@Gaetan-Gauth I am not sure I understand your results. Why is the first plot so low (10⁻³³) but your older plots had 10⁻⁶ or so ? I dont understand why you use the distribution function. Does it tell you the magnetic jump is gone ? Basically I don't understand if you are saying that the problem has disappeared or not.