Major to-do list

[mattsignorelli]

• Add generic element: every element should store basically the same ODE problem but with different params depending on the multipoles
• Decide internal bunch structures: AoS vs SoA. The EnsembleProblem interface gives a very very easy way to do both CPU and GPU parallelization, and appears AoS. Alternatively we could define our ODE on a matrix and attempt to parallelize over that in an SoA layout (discussed as option 1 here) but this may be trickier and possibly slower. EnsembleProblems appears to be recommended for our use case (many initial conditions/parameters but same ODE. I have added a simple example including a new structure for AoS here using the new stucture defined here, however someone should understand this better. We may need different layouts for CPU vs GPU.
• Handling caches/memory pre-allocation: Entering and exiting solve without preallocation can add up over time. Someone should reach out to the SciML developers (Chris?) and ask for advice on how to best do this with our special use case where we have many discontinuous Hamiltonians one after the other. Here is a discourse discussion where someone does this in a loop. Or would it be crazy to have one giant ODE so we stay in solve? These details must be answered and solved carefully or else tracking will be very slow. Another related issue is using the alias_u0 flag in solve - I could not get this to work and it might be a bug or I'm doing something wrong...
• Defining elements: We should check how long the JIT time is for different functions entered into DynamicalODEProblem. If it is costly, then I suggest for simple elements we have a single set of dynamical equations basically includes the solenoid, quad, drift, sextupole, bends fields all together where the params part includes their corresponding field strengths. This will mean only a single JIT for most elements. Higher order multipoles might then need their own functions
• Numerical stability of the dynamics equations: This is an important topic, one issue I solved with the drift is just normalizing to $\beta=0$ and doing some algebra to remove the singularity and also reduce unnecessary calculations. This ensures you never have something close to $0/0$ for very low $\beta$ particles e.g. near a cathode
• Spin tracking: follow the Kepler example
• Array polymorphism: the type definitions should be compatible with as many array types as possible, e.g. StaticArrays and CUArrays. See NonlinearNormalForm.jl or talk to me for inspiration/help