search

SmileiPIC / Smilei

Particle-in-cell code for plasma simulation
https://smileipic.github.io/Smilei
344 stars 121 forks source link

Particle Injector and the boundary conditions #293

Open Tissot11 opened 4 years ago

Tissot11 commented 4 years ago

Description

I attach here some results on using the ParticleInjector block in 1D simulations. In these simulations, I define a drifting plasma and ParticleInjector systematically at both boundaries. I always choose remove boundary conditions for particles. The plasma is drifting in +x direction. The only peculiarity is that for xmax boundary, instead of injecting plasma in -x direction, I'm injecting it +x direction which maybe desired for maintaining the momentum flux if one injects these particles in the simulation box from the left side. I discussed some of these issues in #283 but I thought to present simplified results here and I hope that it can help in diagnosing the problem.

Case 1: Initilize plasma and injection at the left boundary.

PI_3.py.txt Density-3

In this picture and below, I plot the plasma density using the Probe diagnostic. One sees that boundary condition remove is working fine on the right but the injection is located around x=0. One would have expected presence of some plasma in the upper part of figure at later times?

Case 2:

Injection at the left boundary and intilize plasma closer to the right boundary.

PI_4.py.txt Density-4

In this case there should not have been any reflection of particles from the right boundary due to remove boundary conditions. Also the spike at x=0 seen in Case 1 disappears now.

Case 3:

Injection from the right boundary in +x direction and plasma is in whole simulation box

PI_2.py.txt Density-2

In this the boundary condition remove doesn't seem to work on either sides. It seems as if it is overridden by the particle injection which is also assuming injection in -x direction from the right boundary though it is defined in +x direction. I had turned off external magnetic fields in this case compared to other cases, but it didn't have much influence on the issue present here.

Case 4. Injection from thr right boundary in +x direction and plasma is closer to the right boundary. PI_5.py.txt

This is similar to Case 2 above but the ParticleInjector is at xmax. Strangely, I can not run this simulation and always see segmentation fault suggesting problems in VectorPatch.

stderr.8870077.txt

It would be great if you can have a look at these results. 

If available, copy-paste faulty code, warnings, errors, etc.

Parameters

xxirii commented 4 years ago

Dear @Tissot11,

Thank you for summarizing your issues.

So far, I only looked at the first case. You have done a mistake in the namelist. The injection is set to xmax instead of xmin. This explains why you do not inject from the left boundary.

xxirii commented 4 years ago

For the case 2, I observe as well the reflection although I cannot explain why. However, it is normal that you do not have the spike at the left and that you do not have particle injection. It is because your injector is related to the number_density of the corresponding Species. The parameter number_density is defined to start at x = 50 meaning that there is not particle behind the left boundary. To inject particle, please put number_density = 1 for instance in the ParticleInjector.

Tissot11 commented 4 years ago

Dear @Tissot11,

Thank you for summarizing your issues.

So far, I only looked at the first case. You have done a mistake in the namelist. The injection is set to xmax instead of xmin. This explains why you do not inject from the left boundary.

Thanks for the reply. But in Case 1, (PI_3.py.txt), I do inject plasma from xmin. Is there a mix-up in the Namelists?

Tissot11 commented 4 years ago

For the case 2, I observe as well the reflection although I cannot explain why. However, it is normal that you do not have the spike at the left and that you do not have particle injection. It is because your injector is related to the number_density of the corresponding Species. The parameter number_density is defined to start at x = 50 meaning that there is not particle behind the left boundary. To inject particle, please put number_density = 1 for instance in the ParticleInjector.

I followed as described in the documentation. Since the _numberdensity parameter is inherited from the species, I didn't define it again. In some cases I had define it also and I had not succeeded. Besides, I want to inject the same density as the plasma species. So can this _numberdensity parameter for ParticleInjector be equal to the density of the plasma species instead one 1?

xxirii commented 4 years ago

Dear @Tissot11, Thank you for summarizing your issues. So far, I only looked at the first case. You have done a mistake in the namelist. The injection is set to xmax instead of xmin. This explains why you do not inject from the left boundary.

Thanks for the reply. But in Case 1, (PI_3.py.txt), I do inject plasma from xmin. Is there a mix-up in the Namelists?

Sorry my mistake. I found why you have no injection.

You should inject the particles at the same positions:

ParticleInjector(
    name      = "Inj_eon",
    species   = "eon",
    box_side  = "xmin",
)
ParticleInjector(
    name      = "Inj_ion",
    position_initialization = 'Inj_eon',
    species   = "ion",
    box_side  = "xmin",
)

I am not sure about the reason but your ions seem to be quasi-immobile due to the low momentum. Therefore, it is possible that without charge neutrality, they hold back the electrons.

One more thing. I can understand that you though you were injecting at the same location because you have associated injectors to the species but no. Since a species can have several injectors, you have to specify again in the injector to initialize at the same location. I should add this in the doc.

xxirii commented 4 years ago

For the case 2, I observe as well the reflection although I cannot explain why. However, it is normal that you do not have the spike at the left and that you do not have particle injection. It is because your injector is related to the number_density of the corresponding Species. The parameter number_density is defined to start at x = 50 meaning that there is not particle behind the left boundary. To inject particle, please put number_density = 1 for instance in the ParticleInjector.

I followed as described in the documentation. Since the _numberdensity parameter is inherited from the species, I didn't define it again. In some cases I had define it also and I had not succeeded. Besides, I want to inject the same density as the plasma species. So can this _numberdensity parameter for ParticleInjector be equal to the density of the plasma species instead one 1?

If your plasma is not defined until the border of the simulation box, you can still create a parameter called maximum_density that you use to define the internal and injection profile. If you specify a profile that do not cover the boundary, it will not work.

Tissot11 commented 4 years ago

Dear @Tissot11, Thank you for summarizing your issues. So far, I only looked at the first case. You have done a mistake in the namelist. The injection is set to xmax instead of xmin. This explains why you do not inject from the left boundary.

Thanks for the reply. But in Case 1, (PI_3.py.txt), I do inject plasma from xmin. Is there a mix-up in the Namelists?

Sorry my mistake. I found why you have no injection.

You should inject the particles at the same positions:

ParticleInjector(
    name      = "Inj_eon",
    species   = "eon",
    box_side  = "xmin",
)
ParticleInjector(
    name      = "Inj_ion",
    position_initialization = 'Inj_eon',
    species   = "ion",
    box_side  = "xmin",
)

I am not sure about the reason but your ions seem to be quasi-immobile due to the low momentum. Therefore, it is possible that without charge neutrality, they hold back the electrons.

Actually ions should have higher momenta since they are heavier and have the same velocity as electrons in my case.

One more thing. I can understand that you though you were injecting at the same location because you have associated injectors to the species but no. Since a species can have several injectors, you have to specify again in the injector to initialize at the same location. I should add this in the doc.

Thanks for the clarification. Indeed I had thought associating with the species should have assigned the same location as the species itself. Though, I had done some simulations to explicitly assign the initialisation as mentioned in #271. So does the Case 2 works with your workaround? I'll check it tonight.

Tissot11 commented 4 years ago

For the case 2, I observe as well the reflection although I cannot explain why. However, it is normal that you do not have the spike at the left and that you do not have particle injection. It is because your injector is related to the number_density of the corresponding Species. The parameter number_density is defined to start at x = 50 meaning that there is not particle behind the left boundary. To inject particle, please put number_density = 1 for instance in the ParticleInjector.

I followed as described in the documentation. Since the _numberdensity parameter is inherited from the species, I didn't define it again. In some cases I had define it also and I had not succeeded. Besides, I want to inject the same density as the plasma species. So can this _numberdensity parameter for ParticleInjector be equal to the density of the plasma species instead one 1?

If your plasma is not defined until the border of the simulation box, you can still create a parameter called maximum_density that you use to define the internal and injection profile. If you specify a profile that do not cover the boundary, it will not work.

Normally I assign the plasma and the injector to the same boundaries as I mentioned in #283. I wanted to understand the behaviour of ParticleInjector and the boundary conditions because of the problems seen in #283 that's why I created these four cases. In Case 4, I initialise the plasma closer to the right boundary but the injector should inject the plasma in +x direction.

This parameter _maximumdensity isn't mentioned in the documentation. How does it work and where should I define it, in the Species or in the ParticleInjector block? Does it equal to the density of the plasma and density of the injected particle?

xxirii commented 4 years ago

I am still surprised by the behavior of the case 2. I don't understand this reflection and moreover it seems independent from the injection because when I switch it off the reflection still occurs.

Concerning the maximum_density, it's not a Smilei parameter. I was thinking about a Python parameter:

maximum_density = 1.

Species(
...
number_density = trapezoidal(maximum_density , xvacuum=50.0, xplateau=10.0),
...
)

ParticleInjector(
...
number_density = maximum_density ,
...
)

So that you keep the same density everywhere in a simple way.

I did not have a look to 3 and 4 for the moment.

Tissot11 commented 4 years ago

I am still surprised by the behavior of the case 2. I don't understand this reflection and moreover it seems independent from the injection because when I switch it off the reflection still occurs.

Concerning the maximum_density, it's not a Smilei parameter. I was thinking about a Python parameter:

maximum_density = 1.

Species(
...
number_density = trapezoidal(maximum_density , xvacuum=50.0, xplateau=10.0),
...
)

ParticleInjector(
...
number_density = maximum_density ,
...
)

So that you keep the same density everywhere in a simple way.

I did not have a look to 3 and 4 for the moment.

Yes, this reflection is hard to understand. I tried removing magnetic field components and this reflection still persists. Are there some issues with the right boundary? I see some issues while using a moving window and I'll post it in another thread.

mccoys commented 4 years ago

Here is the explanation that I found for this reflection at the right boundary.

As there is some randomness in positions and momenta, electrons and ions will slowly separate, thus creating random electric + magnetic field noise. When the plasma reaches the boundaries, it may happen that an electron is removed, but the ion is still in the box (because the two particles are separated). This creates a local space-charge field that is later compensated by the ion leaving the box. But in the meantime, you have this short-lived noise. When particles come from far from boundary they get more separated over time, thus this noise becomes more and more important. As a consequence, you reach a point where you can add the noise from several particles and it is enough to reflect one other particle. This is why you see particles reflected after a while.

Several ways you can stop this effect:

I know these solutions are not very convenient, but I think the problem is very general in PIC codes. Removing particles is a very tricky business and should be done very carefully. In any case, I believe this is not a problem in Smilei, but a general difficulty in all PIC codes.

phi-kafka commented 4 years ago

Hi mccoys, I have exactly this explanation, just wanted to do some tests. There is another way to avoid the problem, which I was thinking about: When electrons are deleted, they may be stored at the wall. Then their potential will not dissapear, I believe it can work. There is no such option in SMILEI right now, but it worth to try not delete or just stop particles, but stop them and then freeze. What do you think?

mccoys commented 4 years ago

@phi-kafka there is the "stop" boundary condition

MickaelGrech commented 4 years ago

@mccoys I'm not sure the "stop" bc implemented in Smilei will do the work (I explain below). @phi-kafka in theory you are right, we could have a macro-electron stuck there up to the moment an ion gets to it. Unfortunately, I do not see how to do this in Smilei where we use a charge-conserving current deposition (Esirkepov). Hence, a particle that is stopped (not moving) will not deposit current. As a result, I think it just disappears and I'm not quite sure of the difference in that case between stop & remove.

Tissot11 commented 4 years ago

Here is the explanation that I found for this reflection at the right boundary.

As there is some randomness in positions and momenta, electrons and ions will slowly separate, thus creating random electric + magnetic field noise. When the plasma reaches the boundaries, it may happen that an electron is removed, but the ion is still in the box (because the two particles are separated). This creates a local space-charge field that is later compensated by the ion leaving the box. But in the meantime, you have this short-lived noise. When particles come from far from boundary they get more separated over time, thus this noise becomes more and more important. As a consequence, you reach a point where you can add the noise from several particles and it is enough to reflect one other particle. This is why you see particles reflected after a while.

Several ways you can stop this effect:

  • Do not compute fields (see time_fields_frozen). But this is not a PIC code anymore !
  • Remove any randomness to ensure electrons and ions follow each other perfectly : reduce temperature, use "regular spacing".
  • Stay far from boundaries.
  • Have many more particles so that the noise is reduced.

I know these solutions are not very convenient, but I think the problem is very general in PIC codes. Removing particles is a very tricky business and should be done very carefully. In any case, I believe this is not a problem in Smilei, but a general difficulty in all PIC codes.

Thanks for your valuable insights! I'll try your suggestions. I still wanted to clarify few things.

  1. Does usage of regular spacing is as good as random, in a sense that one do not expect any other issues to arise while choosing regular spacing for plasma particles? I have only used random initialisations so-far.

  2. Staying away from the boundary means choosing a larger simulation box can help? Or in a larger box the accumulation of noise will become a problem too?

I'll take more particles for simulating this kind of problems. Thankfully Smilei is super-fast so more particles in a simulation are not an issue on our cluster.

Tissot11 commented 4 years ago

@mccoys I'm not sure the "stop" bc implemented in Smilei will do the work (I explain below). @phi-kafka in theory you are right, we could have a macro-electron stuck there up to the moment an ion gets to it. Unfortunately, I do not see how to do this in Smilei where we use a charge-conserving current deposition (Esirkepov). Hence, a particle that is stopped (not moving) will not deposit current. As a result, I think it just disappears and I'm not quite sure of the difference in that case between stop & remove.

Just out of curiosity, can applying field filtering around the boundaries help in this situation? Or this is not feasible?

mccoys commented 4 years ago
Tissot11 commented 4 years ago

Thanks! Do the other issues noted in Cases 3 and 4 can also be taken care of if one tries your suggestions?

mccoys commented 4 years ago

Case 3 seems to be just the same without injection, so same as case 2 in my opinion.

Case 4 clearly shows a bug in xmax injection. @xxirii did you notice this bug before.

xxirii commented 4 years ago

I did not study case 4 for the moment but I did not see it before

phi-kafka commented 4 years ago

Hi @mccoys @MickaelGrech , I tried "stop" bc" for idential case, and it shows that what I proposed may work:

When electrons are deleted, they may be stored at the wall. Then their potential will not dissapear, I believe it can work. There is no such option in SMILEI right now, but it worth to try not delete or just stop particles, but stop them and then freeze.

Please, see the comparison: image

  1. For "stop" there is no reflection, but particles accumulate near the boundary.
  2. It is something artificial coming from the left side, because of interferention of injection and "stop" bc
  3. The particle are not frozen at the boundary, when they stop there. That is why they come back when their density increase, due to fluctuations and numerical heating there in high density region.

So, if there would be not just "stop" but "stop"+freeze (switch of pusher for these particles), it would be perfect in this case.

This solution would perfectly work in 1D. In 2D and 3D this would not be perfect, (but still much better). For instance, for non-neutral beam near the boundary it would create artificially stronger fields in the box, than it would be when particles escape the box and go to infinity.

For that 2D and 3D situation it can be more elaborated, For this case it can be (as I understand how it may work, @MickaelGrech ), that particles are not deleted at the boundary, but just frozen with their velocities at the boundary cross. Till they pass, say, several Debye length - this distance may be a parameter in the namelist, their fields do to charge and current densities are adding to the box fields. It is fast because just need to have density and current, then it is analytically integrated.

I understand that it is a serious work, so why don't just add the bc "stop + frozen"?

mccoys commented 4 years ago

As Mickaël said, frozen particles do not participate in any equations in the code. They would have no effect at all on the fields. The Poisson equation is only solved at the beginning of the simulation but never after. This means that charge has no effects. Only current has. Frozen particles have no current.

phi-kafka commented 4 years ago

Yes, indeed, thank you @mccoys, I of course missed this step in my proposition. Yes, it is not enough just to freeze particles, after a "prescribed field" created by these particles should be applied to the box according to the Coulomb law. It is a summation over the surface cells of an analytical expression, seems to be numerically fast. But it is more work with the code, almost the same development work as in the second proposed approach "For that 2D and 3D situation it can be more elaborated...". Just want to mention, that in some setups, it is really important. My example: I look on the effect of plasma surrounding a wire with magnetic field. I have initially dense plasma, but then it goes out from the box, and only ~10^-3 is magnetized. These ions which are not passing the boundary create much more effect on plasma than those who stay near the wire. Increasing the box size is also not easy in 3D.

mccoys commented 4 years ago

Boundary conditions are very often a huge difficulty in all sorts of simulations. What you are asking is basically to create a field that is equal to that of particles exiting the box. If you just make the field of particles stopped at the border, this is NOT equal to the field of escaping particles. Instead, particles that would travel far away would contribute less and less. As a consequence, what you are asking is to do another simulation outside of the box, which is not reasonable.

Generally, the solution would be to make the box much larger and stop the simulation before the perturbation has come back in the center. Otherwise, you have to invent a new boundary condition that would absorb better particles and their fields.

Tissot11 commented 4 years ago

I ran some initial simulations choosing higher number of particles and regular initialization does seem to improve things a bit. If I understood correctly, then these fluctuations might not be an issue if we inject and choose to simulate electron-positron plasmas? Also could choosing smaller grid size and time step also help in easing this sort of fluctuation?

phi-kafka commented 4 years ago

@mccoys, sorry, I don't agree that I ask for another simulation outside. There is no need of grid, no need of self-consistency. Knowing particle velocity you just have analytical expression x=x_0+v_0 t, x_0 and v_0 are coordinates and velocities at the boundary. Yes, you are right that particles will contribute less when they far, but only in 2D and 3D. That is what I told before. In 1D it would work perfect even if they stay at the boundary. In 2D and 3D what is written in https://github.com/SmileiPIC/Smilei/issues/293#issuecomment-702408177 can make it more precise. Actually this is what you say "invent a new boundary condition", though not ideal of course. Maybe not considering it now, but please keep that solution in mind.

mccoys commented 4 years ago

@Tissot11 I haven't tried e-p plasmas. But it might help. I don't know if resolutions can help, but maybe.

@phi-kafka Assuming ballistic positions is still somewhat disturbing in a PIC simulation. But it might do things better. Who knows. For now, we cannot dedicate time to this subject.

Tissot11 commented 4 years ago

@mccoys Ok. I'll try e-p plasmas to see. Any thoughts if the thermalize boundary condition can be useful in this type of situation?

xxirii commented 4 years ago

Do not use regular spacing if you have nonrelativistic plasmas. The explanation has been revently added in the doc.

Tissot11 commented 4 years ago

I have always non-relativistic plasma. Following @mccoys earlier suggestion I tried regular initialization for plasma species to improve the fluctuation problems with Particle Injector. I didn't see the explanation yet.

mccoys commented 4 years ago

The issue was unrelated to the injector. It is only boundary condition stuff. I gave an explanation above.

Tissot11 commented 4 years ago

Ok. I think these issues with Particle Injector are not really sorted out yet. Following your suggestions, I ran some simulations with more 60000 particles per cell and regular position initilization. I attach the Namelist of one of the simulations and the figure on the electron density injected into the simulation box. You can see that despite I initialise Particle Injector to be with density 1.0, it doesn't inject that density into the simulation box on the left hand boundary. This was related to Case 1, I had posted in the begining. Also choosing these insanely many particles didn't improve the situation.

Electron density HN-py.txt

I have also done another simulation with that many particles per cell but without the reflecting wall on the right hand side. I'll always see an issue about the buildup of the plasma density and fields on the right hand boundary. I had chosen Case 3 above to mitigate this charge build up on the right hand side, since I was injecting the plasma at the right boundary in +x direction and hoping that it can take care of some of the charge build up at the right boundary. Case 4, as you also said that it a bug, so I hope it can be fixed soon. I'm afraid that currently Particle Injector in SMILEI has few issues to resolve.

xxirii commented 4 years ago

Hi @Tissot11, I do not recommend the use of the regular init for non-maxwellian distributions. This is explained by hand in this link : https://smileipic.github.io/Smilei/particle_injector.html. With the way we inject particles, regular spacing act as a filter that only let the particles with high energy pass the boundary. The velocity of each particle has to be large enough so that the particles reach and cross the boundary. The consequence is a lower density than expected and a truncated distribution. Please, use rather the random init in your case. For the other issues, we have to work on it indeed.

Tissot11 commented 4 years ago

Thanks for your answer and sending the link! I think your explanation on this page make sense! I'm now running the simulations with random initilization and more number of particles. I wanted to ask you two more questions:

  1. If I just inject plasma using the Particle Injector block and keep the density of the species and particles per cell to be zero in the simulation box but I give finite values in the Particle Injector itself. Then I'm injecting only the plasma from the boundaries and I don't have any plasma already present in the simulation box. In this case also, will I see some fluctuations in the injected density at the left boundary for non-relativistic plasmas?

  2. If I assign a Particle Injector on the right boundary that injects the plasma out of the simulation box in +x direction as in cases 3 and 4, should't it help to alleviate some charge build issues at the right boundary?

Thanks again for looking into these issues. I really like Smilei and if I can use Particle Injector without issues, it will be really helpful for my work!

xxirii commented 4 years ago

Dear @Tissot11

Good questions and I don't really know in fact. This injection is really non intuitive for me.

Tissot11 commented 3 years ago

I was wondering if you have any updates on the Particle Injectors issues? I also found that if I use trapezoidal time function for Particle Injector and if I turn off the injection after a fixed time, then the simulation crashes. Additionally, If I freeze the particle motion in the beginning of the simulation, then it's a also a problem with the Particle Injector block. I don't know if it's a bug or is it by design?

xxirii commented 3 years ago

Dear @Tissot11 ,

Concerning the crash situation can you provide a namelist. I will try to investigate soon.

For the moment no more improvement or correction for the injection.

For the freezing, I think the 2 are incompatible. You can not inject a frozen species. However I have never tried and I should add an error message to prevent this.