TeamAtomECS / AtomECS

Cold atom simulation code
GNU General Public License v3.0
46 stars 12 forks source link

s-wave collisions #14

Open ElliotB256 opened 3 years ago

ElliotB256 commented 3 years ago

To simulate evaporation in optical/magnetic traps, AtomECS requires a way to calculate the effects of short-ranged (eg, s-wave) collisions. Given that the collisions are short-ranged, all methods used should a spatial binning method for speed. Beyond that, there are subtle variations in the implementation details, and we could support a number of approaches:

  1. Monte-Carlo method: determine a number of collisions to occur within a box, eg based on density and mean velocity. Perform collisions between velocities of random atomic pairs.

    • How to do mixtures with this, eg fermions+bosons? Sr87+88?
    • Will be faster, because it's one loop per box with comparatively few collisions. Scales as O(N), where N is average atom number per box.
  2. 'Pair-wise' Monte-Carlo method; Iterate over each pair within a box. Calculate the chance for a collision to occur, then draw a random number to determine if it did.

    • Easier to do mixtures, because you can specify collision lengths between pairs_
    • Will be slower, and scales O(N^2) where N is average atom number per box.

All models should also have the ability to specify some kind of 'macro atom', so that the behaviour of a large number of atoms can be represented by a numerically tractable ensemble. There's a few ways we could do this:

  1. Have a simulation resource which specifies a global scaling parameter.
  2. Specify per-atom MacroAtom components.

I like the global scaling parameter because it is straight forward. Also, I can see how we might be able to simulate multiple orders of magnitude of atom loss by rescaling the parameter as atoms are lost. For instance, this could allow one to start with 1e4 atoms representing 1e10 atoms and end with 1e4 atoms representing 1e4 atoms.

ElliotB256 commented 3 years ago

Worth noting that collision length isn't really a per-atom quantity - it's defined between a pair of atoms - so we can't really have a component for it.

How to store collision lengths for calculations? I would suggest we have a table which holds the collision lengths for all species, and index this table using per-atom collisions::Species { id: byte }. The table should form part of the Resource option which is used to enable/disable s-wave collisions.

ElliotB256 commented 3 years ago

@d-garrick Note that a more general implementation can't just swap velocities, it needs to swap momenta in the COM frame.

Also, it can't just swap momenta in the COM frame either - otherwise that would imply all x-components of velocities are decoupled from all the y-components, which is unphysical. Look up references and prev work to see what model is best to use.

MauriceZeuner commented 3 years ago

Hi @d-garrick ,

I'm Maurice and I'm currently doing the dipole-force branch. Since we are ultimately interested in simulating evaporation (for example in a dipole trap) - I wanted to ask if you are interested in a brief meeting about the collisions. I have implemented a really lame version of them in python beforehand (just to get a feeling for it) and I'd be interested to learn how you implemented them and how to use them. Of course, I can also give you the outline of how the dipole force will work - although it is very similar to the cooling light and forces functionality already in the master branch.

Anyway, I think it would be great to meet so we are aware of what everyone is doing, what do you think?

Greetings, Maurice

d-garrick commented 3 years ago

Hi Maurice,

We can have a meeting about this - the collisions code is still unfinished and requires more testing, but I'm happy to have a chat about it with you. I can probably chat on teams tomorrow (Thursday) or Friday. When are you free?

David

------ Original Message ------ From: "Maurice Zeuner" @.> To: "TeamAtomECS/AtomECS" @.> Cc: "d-garrick" @.>; "Mention" @.> Sent: Tuesday, 23 Mar, 2021 At 20:52 Subject: Re: [TeamAtomECS/AtomECS] s-wave collisions (#14)

Hi @d-garrick https://github.com/d-garrick , I'm Maurice and I'm currently doing the dipole-force branch. Since we are ultimately interested in simulating evaporation (for example in a dipole trap) - I wanted to ask if you are interested in a brief meeting about the collisions. I have implemented a really lame version of them in python beforehand (just to get a feeling for it) and I'd be interested to learn how you implemented them and how to use them. Of course, I can also give you the outline of how the dipole force will work - although it is very similar to the cooling light and forces functionality already in the master branch. Anyway, I think it would be great to meet so we are aware of what everyone is doing, what do you think? Greetings, Maurice — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/TeamAtomECS/AtomECS/issues/14#issuecomment-805246420 , or unsubscribe https://github.com/notifications/unsubscribe-auth/AST4AXYLLOK6NJQGS55WAVLTFD5QHANCNFSM4XY6IOCA .

MauriceZeuner commented 3 years ago

Hi David,

great, thanks! I am very flexible tomorrow - Friday before noon as well. Would you like to suggest a time that works best for you?

Best wishes, Maurice

d-garrick commented 3 years ago

Shall we say 11am tomorrow on Teams?

David

------ Original Message ------ From: "Maurice Zeuner" @.> To: "TeamAtomECS/AtomECS" @.> Cc: "d-garrick" @.>; "Mention" @.> Sent: Wednesday, 24 Mar, 2021 At 19:30 Subject: Re: [TeamAtomECS/AtomECS] s-wave collisions (#14)

Hi David, great, thanks! I am very flexible tomorrow - Friday before noon as well. Would you like to suggest a time that works best for you? Best wishes, Maurice — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/TeamAtomECS/AtomECS/issues/14#issuecomment-806096898 , or unsubscribe https://github.com/notifications/unsubscribe-auth/AST4AX5F7KXNLWTHEE7L6QLTFI4TTANCNFSM4XY6IOCA .

MauriceZeuner commented 3 years ago

Sounds good - I think we don't have a common teams channel yet. You should be able to find me by my name "Maurice Zeuner" or ----.

See you, Maurice

ElliotB256 commented 3 years ago

image

This graph shows agreement of the MC approach with theoretical rates as a function of collisions/timestep/box and atoms/box. Note that even 30% difference is expected due to inaccuracies in the numerical calculation of the 'expected' rate (I'm using only approximate factors from kinetic theory, there's a few factors missing).

ElliotB256 commented 3 years ago

(The above simulation is for a box trap simulation domain with 10 boxes along each axis for 1000 boxes total.)

ElliotB256 commented 3 years ago

Slightly more correct version of the same: image

ElliotB256 commented 3 years ago

@MauriceZeuner

MauriceZeuner commented 3 years ago

Is there a minimum example that shows how to use the collisions? I'd be very interested to throw some particles into a dipole trap and evaporate.

ElliotB256 commented 3 years ago

@DavidGarrick is still working on it

d-garrick commented 3 years ago

I can add a simple example to the collisions branch tomorrow. It produces accurate collision rates for thermal clouds but haven't tested evaporation yet.

ElliotB256 commented 3 years ago

@DavidGarrick I know we spoke a few times today, could you post a brief update here? :)

d-garrick commented 3 years ago

I've been testing collision rates and evaporative cooling and had some discrepancies but have found the bug - the current version overstimates the collision rate by a factor of the macroparticle value. After fixing this, the s-wave collisions (and magnetic trapping) seem to work well and agree qualitatively with experimental data: image This is a phase space density curve of evaporative cooling in the Time-Orbiting Potential trap, which I've also implemented on the magnetic_force branch. This experimental data is quite old and not perfectly suited to comparison to simulation, so I'm hoping to take some fresh data explicitly for the purposes of comparing to simulation, but the agreement is still reasonably good.

The model for collisions we use assumes that we have identical particles, that the velocity distribution in the collision cell is at least nearly thermal, and that we can model a cell as a uniform density. For most scenarios and a suitably small cell width these last two are probably good assumptions. We may want to modify the existing code to handle mixtures of species.

I also noticed some other interesting bugs. One is that the collisions must be carried out after the velocity verlet integrator has finished, or else the integrator will no longer conserve energy (in quite a major way).

d-garrick commented 3 years ago

It's also worth mentioning that at the moment at least, the user has to spend some time figuring out good parameters to use to ensure accuracy of the simulation. For example, the timestep must be significantly less than the average (or possibly the smallest, in the case of inhomogeneous density) collision time, at least 10 times less. I have also read in the literature that the timestep needs to be much smaller than the mean transit time; that is, the average time taken for a particle to cross a collision box. However, I have not personally noticed any unphysical effects due to this, and I am actually not sure why it would matter on a theoretical basis. The collision box width needs to be large enough to have at least >10 but preferably >20 particles, but too large and the simulation will be slower and correlations between position and velocity will be increasingly coarse-grained. Angular momentum is not conserved in DSMC collision algorithms, so any simulation in which this is of interest is likely to be inaccurate. Small particle number per box is of course fast, but leads to large statistical fluctuations. A time-averaged (or spatially averaged, over many identical collision boxes) collision rate will be accurate, but the fluctuations need to be small enough to average to close to the correct value faster than any relevant dynamics in the trap.

ElliotB256 commented 3 years ago

the timestep must be significantly less than the average ... collision time, at least 10 times less.

Can you produce a 2D plot like the one above to show when it breaks down?

Graph looks great!

ElliotB256 commented 3 years ago

@DavidGarrick discussed in slack the advantages of making a PR for these early collision features, with more configurations exposed/optimisations in a second PR.

ElliotB256 commented 3 years ago

@DavidGarrick ready to pull?

MauriceZeuner commented 3 years ago

Hi David,

I think you mentioned some unit-tests still failing - just to check:

collision_rate() --> fails test_collisions() --> deaktivated/fails if ran

Also, there is a panick that I don't understand where it actually occurs:

<unnamed>' panicked at 'okTried to fetch resource of type 
`CollisionsTracker`[^1] from the `World`, but the resource does not exist.

Is that the currently "correct" behaviour or did I mess something up while merging? Thanks :)

d-garrick commented 3 years ago

The unit tests are very unfinished so I can't comment much on which particular ones are failing, they won't be done until I try and merge collisions (which I'll try and work on soon in between lab work). Let's email about any errors or questions and if we find an actual bug we can mention it here.

ElliotB256 commented 3 years ago

I saw @d-garrick working on this again today

ElliotB256 commented 3 years ago

Just out of interest (and when the experiment is operational again) - did you ever take a more detailed set of data for the evap ramp and try fitting to see what the most realistic cross section for Rb is? eg, ~100 different values for atom number, find the cross section which gives best agreement with experimental data. Then compare that to the expected cross section from theory?

d-garrick commented 3 years ago

I haven't taken that data - might take it once the experiment is working again. Could try, rather than evap, damping rate of a breathing oscillation?

ElliotB256 commented 2 years ago

Any chance to revisit that data now the experiment is operational?

ElliotB256 commented 2 years ago

Tagging @BrianBostwick (hope it lets me)