External field

[Tissot11]

Description

I'm trying to initialise external magnetic field using trapezoidal profile. However, at the left boundary of the simulation, I see a ramp in magnetic field strength despite I have set xslope1 and xslope2 to be zero. On the right boundary I don't see any ramping but only on the right side. I attach a snapshot and

paste some lines of namelist below. I'm also injecting particles but the issue using the random initialisation as suggested #271.

EM_boundary_conditions = [ ['silver-muller','silver-muller'] ] ,

ExternalField(
    field = "Bx",
    profile = trapezoidal(0.03548, xvacuum=0.0,xplateau=400.*l0,xslope1=0.0,xslope2=0.0)
)
ExternalField(
    field = "Bz",
    profile = trapezoidal(0.03548, xvacuum=0.0,xplateau=400.*l0,xslope1=0.0,xslope2=0.0)
)
ExternalField(
    field = "Bx",
    profile = trapezoidal(0.0357,xvacuum=400.*l0, xplateau=Lsim,xslope1=0.0,xslope2=0.0,)
)
ExternalField(
    field = "Bz",
    profile = trapezoidal(0.1401,xvacuum=400.*l0, xplateau=Lsim,xslope1=0.0,xslope2=0.0,)
)

Steps to reproduce the problem

If relevant, provide a step-by-step guide

Parameters

• Smilei 4.3.2
• Input file ?
• Commands or script used to run the faulty simulation ?
• Computing resources
  • type of computer
  • OS
  • C++ compiler (g++ 4.8.5)
  • MPI version (g++ 6.2.0)
  • HDF5 version 1.8.21
  • Python version (2.7.5)
  • ...

[mccoys]

You see a ramp between the first 2 points of the grid. This is the limit of your resolution.

[Tissot11]

Thanks for the quick reply! Yes, I use a resolution of dx=0.1. But I would have thought that B should be constant at the boundary also. Actually it occurs in the middle of the simulation too where I define different value of the same field. I can make the resolution smaller but I can't go really small otherwise computation will take longer. What is the best way to deal with this? Should a define xvacuum=-dx so that the field starts at a constant value at the left boundary and also in in the middle?This ramping generates an additional component of the electric field causing unphysical results.

[mccoys]

If you want a constant field even past the boundary, you should set xvacuum to a large negative value, and let the plateau be much longer than the box. In the middle of the box, there is nothing you can do. It is impossible to have shorter changes then the resolution.

[Tissot11]

Ok. I could get away from the ramp on the left side by defining a python function profile. Is there any way to smoothen the fluctuation electric field that might arise due to this ramp in magnetic field in the middle of the simulation box?

[mccoys]

There are some filters available, but beware it might create unrealistic fields.

[Tissot11]

Ok. I'll check and compare the results. Friedman filter is only available in 2D? I follow the strategy of Plotnikov et al. (2017).

[Tissot11]

I have checked with filters and it made things even worse with unphysical results. I have another question related to the particle injector. I initialized the injector based on the discussions in #271, and I attach the namelist. I define plasma species (electrons and ions) in the simulation box and I'm also injecting these electrons and ions from the left boundary for the continuous flow. The number density is fixed at 1.0, but I see it is not 1.0 (see the figures attached) at the left boundary and shows large fluctuations. Why is this the case? I have played around with different strategies and used trapezoidal fuction with xvacuum to be large negative number. But with the particle injector block on, it always shows unusual behaviour (Ex is also not zero). Am I not understanding how the particle injector block is meant to be used. Please let me know.

I'm running the version of Smilei, I fetched from GitHub since the official tarball releases (4.4 and 4.5) I couldn't compile.

namelist.txt

__

[mccoys]

I am not very familiar with the injector, but I recall that particles are taken from a velocity distribution that may contain negative velocities. Those are discarded. This means there is a fluctuation due to the randomness of this sampling.

[Tissot11]

Thanks for the reply. Your explanation makes sense. The puzzling aspect was that I am initialising a supersonic flow, so i was expecting a lower number of particles with negative velocity. Anyway, is there any way to smooth the initial conditions in Smilei? I know about the current and field filtering but they work for the whole simulation duration. I could be interested in smoothening the initial conditions for the simulation. Do you think if this can be accomplished by using the restart feature of Smilei?

[Tissot11]

To suppress these fluctuations, I tried to do time averaging by choosing time_average= 10 and saved the density in Field and the weight in ParticleBinning diagnostics. The results from two are different and I attach it here. How one can understand this?

[mccoys]

I am not sure how you did your diagnostics here, and I am not sure what is the difference you want to highlight. In general, Field diagnostics obtain the density of particles using a projector like the one that is used to project the currents on the grid (usual PIC stuff). The particleBinning are very different, because they only compute a sum of particle weights in some user-defined binning.

In the future, please use a different issue for each question you have. In this issue you have asked many unrelated things and it is difficult to follow.

[Tissot11]

Sorry about these too many things. Some of these things are related and I'm not sure if one is not causing problems with others. I wanted to check is the profiles functions used to initialise plasma, fields and also use Particle Injector block. To summarize, I see few points that I don't understand.

1. I used your trick with trapezoidal functions that you had explained before for the external field and plasmas. It works well on the left boundary (for external fields but not for the density) but doesn't work on the right boundary for both fields and densities. I write below the lines I used and the figures to explain my point.

ExternalField(
field = "Bx",
profile = constant(0.0358, xvacuum=-2dx)
)
ExternalField(
field = "Bz",
profile = constant(0.140, xvacuum=-2dx)
)
For the ExternalField, I had also tried the trapezoidal function trick which I use below for the densities.

Number density for electrons and ions each

number_density = trapezoidal(3.911, xvacuum=-2dx,xplateau=40.0Lsim,xslope1=0.0,xslope2=0.0),

I also tried

number_density = trapezoidal(3.911, xvacuum=-2dx,xplateau=40.0Lsim,xslope1=0.0,xslope2=100.0),

Please see the figures, below. In my case Lsim=102.4 and dx=0.02. On the right hand boundary, I see fields and density going to zero and I can not control it. For the density I used the Probe diagnostic as it's much faster than Fields. I had used Lsim for defining the corner for the Probe diagnostic to calculate the densities in the figure.

2. For density I thought it was happening due to the fluctuations cause by the Particle Injector block, since even at the left boundary I can't have the correct initialised plasma density (see the figure)

So I tried saving density and fields with Fields and Probe diagnostics and also used ParticleBinning for the density. I had used time averaging to avoid the fluctuation caused by the Particle Injector block. This is what I wrote in the last post. The ParticleBinning diagnostic gives correct value of particle density (1 in my case) but the Field diagnostics does not give correct value of density (~0.10) if I use time_average= 10 for both diagnostics for whole simulation box.

To summarise, I want to know the reasons for the fields going to zero at time t=0 on the right boundary despite profile functions saying otherwise. And problem with the densities on both left and right hand boundaries. I also want to know if I save em fields using Fields and Probe diagnostics then I should get the same results or not.

[Tissot11]

Also does a constant external field is constant in time too?

[mccoys]

• Concerning profiles, have you tried a pure python function, like lambda x, y: 0.140 ?

• I am not sure why the field diag gives 0.1. I will investigate.

• An external field is only applied at the beginning of the simulation. It does not have a time dependence.

[Tissot11]

1. I tried defining some basic Python functions based on conditional expressions. Again on the right side I had the same issue. I can try with lambda function tomorrow. I have also tried the polygonal function that you describe in the documentation.

2. Thanks for looking into it.

3. What I see it that the external applied field shows strange behaviour at the left and right boundaries, which is unphysical. After running simulation for a while, on the left boundary it dies while it grows at the right boundary. In the remaining simulation box, it seems to show the physically meaningful behaviour. It should remain constant like at t=0 at left and right boundaries.

[mccoys]

1. The issue with profiles is actually not related to profiles ! It is due to the non-periodic box, and to the projection scheme. Let me explain. The projection of particle's quantities on the grid assumes some "spread" of each particle (the shape function), related to the order of interpolation (this is classic PIC stuff). When the plasma is uniform, the spreading of each particle will be exactly compensated by the spread of all surrounding particles, leading to uniform grid quantities. But when the plasma is not uniform (say, for instance, a jump in density) there will be a smoothing due to this spread. Now, when the box is not periodic, Smilei must assume there are no particles on the other side of the boundary, meaning it is a jump in density. The interpolation scheme will see it as a ramp in density. Hopefully, this density is only used in diagnostics, it has no effect on the physics. Now, I am not sure how this will appear for the currents as they do have an impact on the physics. I have seen that you have drifting particles in your box. @xxirii @MickaelGrechX do you know how to handle this boundary issue ?

2. Concerning the time_average issue, I just found out that this feature was not done for densities ! Sorry for the inconvenience. I will fix it soon.

3. Concerning external fields I cannot reproduce your problem. Could you provide an example file ?

[xxirii]

Dear @Tissot11, I have run your namelist to test the injection. I don't have the problem at the boundary. Is it something that happens at a given iteration ?

For instance, I have plotted the density from the particle binning diagnostics and the charge from the field diagnostics at timestep 1600:

Note that I agree with @mccoys explanation concerning fluctuations. Increasing the number of particle per cell at injection may improve the situation.

[Tissot11]

1. I understand that drifting plasmas can be tricky. If I remember correctly that I also tried last week periodic boundary conditions but some issues were still there. Anyway for particle density, the drift can be an issue as you explain but not for the external field applied in the simulation. They should not decay at the right boundary.

2. Thanks for the clarification.

3. Last night I ran a simplified set-up of the simulation I was running before. I attach the namelist here and some snapshots. To summarise I see issue with external field (Bz) at the right boundary at t=0 and later on both boundaries. It decays at the left and grows closer to the right boundary. I use Probe and ParticleBinning to sum up the densities of ions and electrons, I see difference in magnitude for in both diagnostics. But strangely this time at the left boundary and the reason of the difference is the ion density I see a big unphysical spike in ion density after t=1000.

First the namelist

namelist.txt

Bz profiles

Densities. First Probe and then ParticleBinning

[Tissot11]

Dear @Tissot11, I have run your namelist to test the injection. I don't have the problem at the boundary. Is it something that happens at a given iteration ?

Probably I had attached a wrong namelist before highlight the problem with time_average=10 with both ParticleBinning and Fields diagnostics. But now this issue has been settled by @mccoys this morning. Please try my namelist that I have just sent and look at the results.

[xxirii]

In fact I have used the namelist you sent 12 days ago. I am only looking at the injection issue. I can see on my computer that the jump can be seen only in the rho diags as you mentioned not in the particle binning. It disappear around iteration 1000 but it may due to the particles coming back from the reflective wall. I think the cause is clearly the interpolation with the boundary behaviour. As soon as the physics is far from it I think it's ok.

[Tissot11]

Dear @xxirii, I have ran at least 30-40 simulations in different setups employing the particle injector block in last two weeks. Sorry that I mixed up namelist before. I've not yet been able to analyze the simulations with a reflecting wall on the right. But please have a look at the namelist that I have just attached 10 minutes before and also have a look at my results.

[mccoys]

Your issue is probably related to the particle boundary conditions, not external fields. The external fields are constant at the beginning so there is no problem. You see a drop on the left side only after a while. The drop at the right side is only caused by the fact that the point x=Lsim is considered outside of the box, so no field is evaluated there (it is always set to 0).

[Tissot11]

Sorry, I'm not able to fully follow you. At t=0 there is a problem on right boundary for the external field even though I have defined xplateau=40*Lsim following your earlier suggestion. After t>500 one starts seeing decline of Bz at the left while it grows at the right boundary, which can't be explained physically. I had also tried using corners=Lsim-dx for the Probe. Though it did help at t=0 but it manifested in a large spike in magnetic field at the right boundary and zero plasma density at later times. After a while density spikes up at the left boundary which is clearly unphysical. I'm using remove boundary conditions for the particle and silver_muller for EM fields.

I'll try to run a simulation without plasmas and Particle Injector block and apply an external field to see what happens.

[mccoys]

• at t=0 : I told you there is no problem at the right boundary. Your last probe point is outside of the box, so it does not record any field.

• at t > 500: You are seeing physical effects. I think that your B-field is transported by the plasma, but this is not caused by the code. You cannot initialize a B-field on particles that have not entered the box yet ! I am sure that this effect is removed if your plasma is not drifting.

[Tissot11]

Dear @Tissot11, I have run your namelist to test the injection. I don't have the problem at the boundary. Is it something that happens at a given iteration ?

For instance, I have plotted the density from the particle binning diagnostics and the charge from the field diagnostics at timestep 1600:

Note that I agree with @mccoys explanation concerning fluctuations. Increasing the number of particle per cell at injection may improve the situation.

Now I looked it these results. Since I haven't put any external fields in this simulation, I only saw that at the left and right boundaries initialization problems. I attach my figures. One can see that at both boundaries the density is around half of what it should be.

[mccoys]

One can see that at both boundaries the density is around half of what it should be.

I answered above this question. It is due to the projection algorithm that assumes there is no particle after the boundary. You have half the value because it makes some sort of smoothing between the plasma and the vacuum outside of the box.

[Tissot11]

• at t=0 : I told you there is no problem at the right boundary. Your last probe point is outside of the box, so it does not record any field.

I have done few more simulations where I didn't use Probe but only Fields to save density and fields and also choosing corners=Lsim-dx. I'll send you figures later after I sort these results out.

How about the spike at the left boundary for ions after t=500? Is this not unphysical?

• at t > 500: You are seeing physical effects. I think that your B-field is transported by the plasma, but this is not caused by the code. You cannot initialize a B-field on particles that have not entered the box yet ! I am sure that this effect is removed if your plasma is not drifting.

Physically, it seems a bit improbable. Because the field is supposed to be initialised at the left boundary and it is constant throughout the simulation at least at the left boundary. In some cases, I can clearly see that it the same effect that you have particle injections. The value initialzed at the left boundary is a bit lower or even zero and somehow because of the drifting plasma, it starts moving to the right. This causes this field to decay. There is no transportation of magnetic field happening that fast on this time scale. I attach some figures to show what I meant.

In this particular simulation, when the Bz becomes stronger at the right boundary, it stops the particle injection on the right boundary after 2e5 iterations. So this is contrary to what you're saying. I attach the namelist and these results on the electron plasma density.

namelist_9.txt

My question would be that if all these projection algorithms can also assume that there is a plasma in ghost cells behind the left boundary. The number of these ghost cells can be increased to get rid of these fluctuations at the boundaries to initialise both particles and external fields with desired values? In other words what one should do to correctly initialise a constant external field and inject a constant flux of particles with a given density at the boundaries during the whole simulation time?

[Tissot11]

Yeah, now it works! Is there some reason behind this particular 2*dx value? I had tried first very large value both for the density and external fields. I noticed some issues with the particle injector, so I choose just 2*dx. Do you have any other suggestions that I should take into account about regarding whatever I have asked till now?

[mccoys]

The default interpolator width (number of ghost cells) is 2.

Concerning the suggestions you ask, I am not sure. Are you talking about the injector setup ? I am really not used to the injection of particles, so I think some others should answer when they return from vacation.

[Tissot11]

Ok. Yeah, I have still few questions about the particle injector but I'll wait until someone replies. Thanks a lot for the help! I did too many runs recently and I have to now understand these better.

[Tissot11]

Sorry again coming back to this external field question. For contact magnetic field it does work. But I wanted to create a jump in magnetic field in the middle of the simulation box. I used both trapezoidal and my own python functions to define this magnetic field Bz. I see that due to this jump there is a numerical instability (both in the middle and at the left boundary) that actually affects the magnetic field applied. Since I don't apply any jump in Bx, it can not show any effect there even for _2*dx_initilization in vacuum. I had hoped this numerical instability arising due to jump i leaves the simulation box after a while since I use silver-muller boundary conditions. This is all in vacuum. Do you have any suggestions on this?

Vacuum-Trial.txt

[mccoys]

A jump of Bz in the x direction makes a curl of the magnetic field. It should create Ey field over time, and then a traveling wave. Do you expect the fields to remain constant ?

[Tissot11]

Yes, your explanation is correct. I expect the Bz field to be constant at and closer to the boundaries. At the location of jump in the simulation box and elsewhere, it will change. I can create a profile for jump in the middle to mitigate this Ey but at the left boundary, I don't know what I should do. I was expecting this Ey to propagate away to the left due to silver-muller boundary condition. It may take a while though.

[mccoys]

I am not sure what is the problem at the left. Your picture shows it constant.

[Tissot11]

The magnitude gets lower. Until x<=40 It should have been 0.035. Due to this instability Ey, Bz becomes much lower even at the left boundary in vacuum. This stair-like structure in the simulation box suggest the formation of two waves at the discontinuity at x=40. Both of them moves in opposite direction. If you run my namelist then you will see that either the wave going to the left due to the discontinuity at _x=40_doesn't respect causality. Or something else happens at the left boundary too and the magnitude of the magnetic field gets lower at the boundary.

[Tissot11]

I guess the solution lies in defining a python function and not using built-in spatial other than constant for the external field? Can time or space resolutions help in suppressing this instability?

[mccoys]

• You forgot to put xvacuum = -3*dx

• This is not an instability. It is correct physics. You simply have an EM wave traveling at the speed of light.

[Tissot11]

I thought it was only needed for the constant function. So this 3*dx initialisation is needed if I use any built-in SMILEI functions for Bz field? For Bx it doesn't matter and even a 2*dx initialisation works because Bx doesn't have any jump. Am I correct in interpreting this?

Would a 3*dx initialisation also be required if I define a Python function like below?

def BzProfile(x):
return 0.035+0.50.105(1+numpy.tanh((x-40.0)))

In the simulations I have done with plasma inlcuded in the simulation, I find the behaviour of the em pulse due to discontinuity in middle according to the expectation. However, in those simulations I see a strange behaviour at the left boundary which is not connected with this discontinuity in the middle. I was using the built-in trapezoidal function to define Bz. So if I try -3*dx initilzation then hopefully I might not see any problem at the left boundary? Also do I need to do the same -3*dx initialisation for defining the plasma density using the trapezoidal function especially for drifting plasma?

[mccoys]

• 3*dx is only needed in functions which implement this xvacuum variable. It means the profile is 0 before that limit.
• For Bx, it works because it does not create a curl of B field.
• 3*dx is not needed in a python function. It does not make sense anyways because there is no such xvacuum variable.
• I tried your namelist, and the problem was caused by Bz=0 at the boundary. You should either use a python profile that does not have 0 at the boundary, or use xvacuum of at least 3*dx

[Tissot11]

Ok. Thanks for tips! Did you try Vacuum-Trial.txt or namelist_9.txt? In Vacuum_trail.txt I had used trapezoidal function without -3*dx. In namelist_9.txt I used polygonal function but fixing the x0 at -2*dx.

[mccoys]

I tried Vacuum_trial

[Tissot11]

Ok. I had done another simulation where I had set xvacuum=-100. This caused also unphysical results. So I should only set xvacuum=-3*dx and not too large value. Sometimes, I also see issue with the left boundary even if I defined a Python function like

_
def BzProfile(x):
if (x <= 40*l0):
return 0.035
else:
return 0.140
_

However, it is inconsistent and seems to depend also on other parameters of the simulation box. Regarding the corners of Probe diagnostic, I had intuitively thought that Lsim should be the corner. But this isn't case as you mentioned. Should I alway set 0.99*Lsim for it to work? Apparently when I increased the resolution of the simulation, choosing corners=Lsim or 0.99*Lsim generate error saying the Probes has no points in the simulation box. Does the numbers of points have anything to do with it and should be smaller than number of points chosen for the simulation grid.

[mccoys]

• setting xvacuum to a large negative value, or using a python function works. If you had a problem, it was due to something else.

• The corners of a Probe can be set to whatever you like. If some points are outside the box, they will be set to 0. If no points are in the box, an error is raised. I don't know why you had an error, but you need to carefully choose your corners and the number of points so that they are in the box. You can have more probe points than grid points if you like, but this is not very helpful usually.

[Tissot11]

1. Yeah the simulation with xvacuum=-100 was the one I spoke few days ago and it could have something to do with the ParticleInjector block. So I’ll wait until colleagues who wrote the ParticleInjector block come back from holidays.

2. I had defined Lsim-dx as the corners and chose Nx=Lsim/dx=102400 to be points of the Probe. Then it showed error about no points being in the domain. Changing the corners to 09*Lsim didn't help either. I then reduced the points to Nx=1000 and the simulation is running fine.

[mccoys]

If you have an issue with probes, please open a new issue with more details on your setup.

[Tissot11]

I have opened the issue and attached the namelist and stdout file.

[mccoys]

Could you recall what was the question with the injector? I am a bit lost

[Tissot11]

So I had two issues.

1. I see that particle injectors seems to inject less or few particles after a certain time on the right boundary but on the left boundary there is a sudden spike. Please have a look at these pictures

2. I suspect that the fluctations in the injected particle density are sometimes causing an issue on the applied magnetic Bz at the left boundary. Like here
   

This was the namelist

namelist_9.txt

I'm running the same simulation again with higher resolution to see if it had any impact on the results. Thanks for looking in this issue!

I forgot to stress that even with 3*dx initilization I see Bz at the left boundary dipping with ParticleInjector block on. But 3*dx initilization somewhat improves the situation a bit.

[mccoys]

• Be careful with profiles at boundaries, as explained before.

• You have positive velocity for species injected on the right. Is this really correct ?

• You have fields and several species. How can you expect them to go smoothly all the way ??? I don't have time to analyse the physics of your system, but you seem to have simply physical effect you need to understand.

[Tissot11]

Yes, I have positive drift for the plasma and I agree that it's a complicated system but I am trying to reproduce results that have been published before using other codes. This is the reason that I suggested that there can be some problems. I have done some test with stationary and drifting plasmas (without particle injections) and external fields and things look reasonable.

[xxirii]

Hi @Tissot11,

I have tested the last version of your namelist on my small computer until time t = 300. There is indeed something strange. I can observe as well a spike growing at the right boundary. Particle injection may not be adapted in your case using an external field near the boundaries. I am not sure injecting particles as we do is correct in this case. Behind the boundaries, the magnetic field do not exist and the plasma is thermalized. In principle, the physics should happen far from the boundaries.

Do you have examples in the literature of people doing the same thing with an injector ?

[Tissot11]

Thanks for the quick answer! I'm also a bit puzzled and I would like to understand the source of this problem. I attach some papers here. My setup is similar but now I'm trying to reproduce these exact results. Please have a look.

Full particle simulation of a perpendicular collisionless shock- A shock-rest-frame model-EarthPlanetsSpace_58_e41.pdf

Dynamics-and-microinstabilities-at-perpendicular-collisionless-shock-A-comparison-of-largescale-twodimensional-full-particle-simulations-with-different-ion-to-electron-mass-ratio2014Physics-of-Plasmas.pdf

[xxirii]

I see, I have to look at the injection methods maybe similar to what was suggested in issue #271. This is a long work and I cannot focus on this for the moment. I will try at least to get more information about the methods when I can.

[Tissot11]

I had tried your workaround that you suggested in #271. I understand that you may have other things on your desk. It would be nice to have an understanding of the issue. Why there is an accumulation of particles on the right and sometime on the left boundaries, and why it seems to affect the external field on the left boundary.

[mccoys]

I am closing this issue as #293 is more specific on the remaining problems.