Field Dependent particle activation #86

Open keithbennett opened 3 years ago

In GitLab by Ernesto Barraza on 05 Sep 2020, 01:56 (GitLab issue GL#2244)

Hi EPOCH developers,

I had a question. I am trying to have a cluster of ions or electrons basically frozen until a strong enough field (in this case a laser) activates them (frees them).

It seems to do this, I have to go in and add things to the code. My question is, to accomplish this is it better to do this within the particle pusher maybe make a subroutine? Or would it be better to make a whole new physics_packages for this? Which would be the quickest dirtiest way to get started and what would be the best permanent solution?

Thanks, Ernesto

In GitLab by Thomas Goffrey (@TomGoffrey) on 05 Sep 2020, 11:53

Are the frozen particles the only particles in the simulation?

In GitLab by Ernesto Barraza on 05 Sep 2020, 15:06

@tom yes. However, I also want to control the energy and momentum distribution as they are activated.

In GitLab by Thomas Goffrey (@TomGoffrey) on 10 Sep 2020, 16:49

If you are triggering the activation by a laser of known temporal profile, is there a reason why particle_tstart can't be used to control the desired effect?

In GitLab by Ernesto Barraza on 10 Sep 2020, 18:08

Interesting. So if I use that I'll just have to time it for when a certain magnitude hits these electron bunches. I could also use the dist_fn and dist_fn_p{x,y,z}_range to control the initial velocity/momentum distribution, correct?

I'll try this out. But I want to have a more physical model closer to electron field emission from a nanotip which requires E-field dependence on emission current density. I feel the simplest way to do this is either have a function like the particle_tstart but rather it dependent on the E-field such as particle_Estart or to have something like the particle injector function but with E-field dependence and not stuck to the boundaries.

In GitLab by Thomas Goffrey (@TomGoffrey) on 10 Sep 2020, 18:35

Yes, essentially the particles won't be updated at all until t > particle_tstart. I realise that's not exactly what you want, but it's worth trying at least.

To truly get what you want I think you'd probably want to modify the particle push so that it tests the magnitude of the electric field at each particle position. So something like:

diff --git a/epoch2d/src/deck/deck_control_block.F90 b/epoch2d/src/deck/deck_control_block.F90
index c1292606d..7b4c120e6 100644
--- a/epoch2d/src/deck/deck_control_block.F90
+++ b/epoch2d/src/deck/deck_control_block.F90
@@ -337,6 +337,9 @@ CONTAINS
     ELSE IF (str_cmp(element, 'particle_tstart')) THEN
       particle_push_start_time = as_real_print(value, element, errcode)

+    ELSE IF (str_cmp(element, 'efrozen')) THEN
+      ecrit = as_real_print(value, element, errcode)
+
     ELSE IF (str_cmp(element, 'use_migration') &
         .OR. str_cmp(element, 'migrate_particles')) THEN
       use_particle_migration = as_logical_print(value, element, errcode)
diff --git a/epoch2d/src/particles.F90 b/epoch2d/src/particles.F90
index de257e55d..7a901478e 100644
--- a/epoch2d/src/particles.F90
+++ b/epoch2d/src/particles.F90
@@ -106,7 +106,7 @@ CONTAINS
     REAL(num) :: fcx, fcy, fcz, fjx, fjy, fjz
     REAL(num) :: root, dtfac, gamma_rel, part_u2, third, igamma
     REAL(num) :: delta_x, delta_y, part_vz
-    REAL(num) :: hy_iy, xfac1, yfac1, yfac2
+    REAL(num) :: hy_iy, xfac1, yfac1, yfac2, emag
     INTEGER :: ispecies, ix, iy, dcellx, dcelly, cx, cy
     INTEGER(i8) :: ipart
 #ifdef WORK_DONE_INTEGRATED
@@ -324,6 +324,9 @@ CONTAINS
 #include "triangle/b_part.inc"
 #endif

+       emag = SQRT(ex_part**2 + ey_part**2 + ez_part**2)
+       IF (emag < efrozen) CYCLE
+
         ! update particle momenta using weighted fields
         uxm = part_ux + cmratio * ex_part
         uym = part_uy + cmratio * ey_part
diff --git a/epoch2d/src/shared_data.F90 b/epoch2d/src/shared_data.F90
index 5fd07c79b..6d79c2650 100644
--- a/epoch2d/src/shared_data.F90
+++ b/epoch2d/src/shared_data.F90
@@ -85,6 +85,9 @@ MODULE shared_data
   ! block of the deck using 'particle_tstart'.
   REAL(num) :: particle_push_start_time = 0.0_num

+  ! Don't move particles below this electric field
+  REAL(num) :: efrozen = 0.0_num
+
   ! Object representing a particle
   ! If you add or remove from this section then you *must* update the
   ! particle pack and unpack routines

I think the above is pretty self explanatory, but in case it's not what I've done is modify the particle push to calculate the magnitude of the electric field at the particle position, and if it is less than efrozen to skip the update of that particle. I've also added efrozen to the control block parser.

Strictly speaking this uses the position the particle would be at if it moved half a timestep. If this isn't accurate enough you'd have to calculate the electric field at an earlier point. I did it this way in an attempt to keep the modifications minimal, whilst guessing this was probably accurate enough.

The only problem with this approach is that if a particle starts to move but then moves into a region where emag < efrozen it will stop moving. The only way I can think of to circumvent this is to add a per-particle flag that tracks if it has started moving, and if the particle is moving to ignore the efrozen test. To do this you'd have to modify (amongst other things) the particle packing and unpacking routines. This would certainly be a more complex modification.

I should stress, I haven't tested any of what I've written above, but I hope this is helpful.

In GitLab by Ernesto Barraza on 10 Sep 2020, 18:42

Awesome. This helps a lot. I'm sure i'll be back with more questions once I get started but this is a good starting point for me.

The other thought I had was modifying the particle injector package to have field dependence and not being bound to the boundaries only. I haven't looked too much into the injector algorithm but I'm guessing there is a reason why it was made to be set only at the boundaries rather than anywhere in space? Would you mind elaborating on this?

In GitLab by Thomas Goffrey (@TomGoffrey) on 10 Sep 2020, 19:01

That's not really something I can comment on with any authority, but I suspect the answer is simply that injection through a boundary was the requested/required use case. From a quick read through the code it certainly looks like the injectors are implemented as flux through a surface, rather than a more general source term as I think you are envisioning. It's possible one of the other dev team who were involved in the implementation of the injectors may want to comment further.

I should add that due to time constraints we aren't able to offer support on modified versions of EPOCH. It's hard enough supporting a code base we know! Best we can do is to offer advice as above. That being said, I hope the advice offered is useful.

In GitLab by Ernesto Barraza on 10 Sep 2020, 19:21

Thomas, thank you! I understand you are all busy and I appreciate your help and responses!

I'm hoping to modify the code slowly and with easier models and build up from there. So before trying to code up a full field emission model, trying to just get frozen bunch of electrons to become active with a field magnitude. Hopefully I'll get a better idea about the inner workings of EPOCH this way. We'll see how it goes and I appreciate your advice.

In GitLab by Christopher Brady (@csbrady-warwick) on 11 Sep 2020, 12:30

Hi, I'm one of the original developers of the injectors so I'll try and expand on @tom 's points. The injectors do work by calculating the flux through a boundary so they aren't really suitable for creating particles in the simulation domain. There are two reasons for this

1) We want to be able to control the statistics of the particles (e.g. particles per cell) and in order to do that we have to know about the flow of particles out of the injector which is much, much easier to do with the flux through a line rather than a general source region.

2) There's a fundamental physics problem. If you have particles arriving on the edge of the simulation domain then you are already saying that you're not interested in simulating the area outside the domain so it's reasonable to simply specify the properties of the particles to come through the edge arbitrarily. If you want to put particles inside the domain then you have to ask some questions about where those particles come from, how they were created and what effect that would have on fields etc. at the point of their creation. The physics packages for ionisation and non-linear Compton scattering (and I think bremmstrahlung) all create particles within the domain but they also modifiy other things like the current at the location of creation and the properties of the parent particles. There are lots of models that can work like that but they aren't all the same in the details. It would be mechnically quite easy to write routines to create particles at a specified location with specified properties but it likely wouldn't get the right answer for a general problem and would have an unacceptable risk of giving "plausible" but wrong answers.

In conclusion, if you want to write code to create particles within the domain that isn't too hard (and we're happy give you some help and advice, subject to the limits that @tom mentioned above) but there is a real risk that you might get an unphysical answer unless you are very careful about ensuring the other parts of the simulation are changed to reflect the process that has created those particles.

In GitLab by Ernesto Barraza on 11 Sep 2020, 21:18

@chris Thank you! I think I have enough to get started and get my hands dirty. So looks like after I try some stuff with the particle pusher it might be best to go and modify the ionisation physics_package since it's closer to what I'm hoping to accomplish.

In GitLab by Ernesto Barraza on 24 Sep 2020, 15:36

Hi all. I finally got a chance to implement @tom 's code modification with my own code modification in the 1D case. To get around the atoms being stuck after they've been "activated" I used both the critical field activation magnitude and the immobile flag. To do this, I moved the immobile line at the beginning of particle pusher loop (particle.f90) to right after the emag to ecrit comparison with the statement:

#include "triangle/b_part.inc"
#endif

+        emag = SQRT(ey_part**2)
+        IF(emag < ecrit) THEN
+          species_list(ispecies)%immobile = .FALSE.
+        END IF
+        IF (species_list(ispecies)%immobile) CYCLE
        ! update particle momenta using weighted fields
        uxm = part_ux + cmratio * ex_part
        uym = part_uy + cmratio * ey_part
        uzm = part_uz + cmratio * ez_part

I tested it out with all the basic workshop examples and it seems to be working fine. I then tested to see if there was activation of these electrons and indeed the electrons were mobile after seeing an ecrit.

My question would be: is it ok I moved the IF (species_list(ispecies)%immobile) CYCLE line from line 158 at the beginning of the particle pusher loop to the lines after 302? I didn't see anything that could be affected by this. Although, it might be even more helpful if I add another conditional statement to see if ecrit was even initialized. If it wasn't, then use the first original immobile line to cycle. If ecrit was initialized then skip over the first original immobile line to cycle.

Let me know what you all think.

In GitLab by Ernesto Barraza on 25 Sep 2020, 00:38

@tom @chris So I've been trying to see how the physics plays out with this activation. It seems that doing this highlights another problem that I'm not sure what to make of it.

I have run simulation with the original EPOCH with just a localized density of electrons in a simulation cell. Physically, I would expect the electrons begin to spread out to fill the simulation domain and finally escape. Pictures below. However, the electrons continue to stay localized. I'm probably misunderstanding how the fields are generated in a non-uniform plasma.

begin:control
  nx = 200

  # Size of domain
  x_min = -4 * micron
  x_max = -2.0*x_min

  # Final time of simulation
  t_end = 80 * femto

  #stdout_frequency = 10
end:control

begin:boundaries

  bc_x_min = open
  bc_x_max = open
end:boundaries

#begin:laser
#  boundary = x_min
#  intensity_w_cm2 = 1.4e18 * (1/1)^2 * (1)^2
#  lambda = 1 * micron
#  phase = pi / 2
#  t_profile = gauss(time, 2*micron/c, 1*micron/c)
#+gauss(time, 6*micron/c, 1*micron/c)
#  t_end = 3 * micron / c
#end:laser

begin:species

  name = electron
  charge = -1
  mass = 1.0
  temp = 273
  #number_density = 1.e25
  number_density = if( ((x lt 0.0*micron) and (x gt -1.0*micron)), 1e25,0 )
  npart = 5000 * nx

end:species

begin:output
  #name = normal  
  dt_snapshot = 1 * micron / c

  #Properties at particle positions
  particles = always
  px = always

  # Properties on grid
  grid = always
  ey = always
  number_density = always+species

end:output

In GitLab by Ernesto Barraza on 26 Sep 2020, 18:45

Hi all,

So I realized I am misunderstanding the linked list and how species works. It seems that changing the immobile flag was maybe not the right solution.

So I ended up modifying the ionisation code for tunneling in ionise.F90. Specifically, I changed the new and current particle updates of the ionisation process to give me the momentum only in the y-direction and to keep the initial species uncharged. This seems to be working as shown in the pictures below. However, when I add more than one ionisation energy the initial species seems to become ionized (charged) anyways. I'm wondering how that's possible after commenting out the charge subtraction. Additionally, I tried adding a new ionisation feature, basically just copying the tunneling subroutine and adding in my stuff while making sure I had the flags set in the control deck but it did not want to compile. Am I missing something on what to do?

Cheers!

#ifndef PER_SPECIES_WEIGHT
              new%weight = current%weight
#endif
              new%part_pos = current%part_pos
              ! Electron is released without acceleration so simply use momentum
              ! conservation to split the particle
              new%part_p = (/ 0.0_num, ey_part*0.5_num*dt, 0.0_num /) 
        !current%part_p &
              !    * released_mass_fraction(current_state)
              current%part_p = current%part_p !- new%part_p
#ifdef PER_PARTICLE_CHARGE_MASS
              new%charge = species_list( &
                  species_list(current_state)%release_species)%charge
              new%mass = species_list( &
                  species_list(current_state)%release_species)%mass
              current%charge = current%charge !- new%charge
              current%mass = current%mass !- new%mass
#endif

002

004

006

008

In GitLab by Ernesto Barraza on 01 Oct 2020, 17:32

mentioned in issue GL#2257

Warwick-Plasma / epoch

Field Dependent particle activation #86