Proton in a magnetic field, numerical simulation

[Appelros]

TL;DR

The magnetic field does no work but translating its force to velocity does, the solution is not elegant and inexact.

Hi! I am developing a fusion reactor simulator and when a proton spins in a perpendicular magnetic field the Lorentz force qv×B is applied at a right angle, however whenever you add a value to a vector in this way its magnitude increases, which contradicts that the velocity of the proton should remain constant since the magnetic field does no work. To solve this I simply scaled the vector back to its original size after adding the acceleration. All my units are set to 1 so the proton should spin in a circle with radius 1, which it does, but with two problems.

1. There are now two sets of acceleration vectors in my code, those that need scaling and those who don't. This feels very ad hoc.

2. While the proton spins in a circle of radius 1 it does so from 0.995 to -1.005 and 1.9999 to -1e-5. While the latter is a somewhat acceptable error rate the former is not (I'm assuming it should spin from 1 to -1). It continues to spin within these parameters for many rotations and under different calculation methods so it is not an error accumulation from my numerical algorithms.

Is there a more elegant solution that solves these problems?

 

[berkeman]

Welcome to PF.

when a proton spins in a perpendicular magnetic field the Lorentz force qv×B is applied at a right angle, however whenever you add a value to a vector in this way its magnitude increases,

When you add value to a vector in what way? The centripetal force from the Lorentz Force is perpendicular to the velocity. Can you explain in more detail how you are calculating the changes in the velocity vector at each time step in your simulation?

I am developing a fusion reactor simulator

A fusion simulator will not be making calculations on each proton in the simulation. You should be working more along the lines of a plasma physics simulation, no?

https://www.amazon.com/Introduction...olled-Fusion/dp/3319793918/?tag=pfamazon01-20

 

[Appelros]

Welcome to PF.

Thanks!

When you add value to a vector in what way? The centripetal force from the Lorentz Force is perpendicular to the velocity. Can you explain in more detail how you are calculating the changes in the velocity vector at each time step in your simulation?

Say one vector is along x, if we add along y we form a right angle triangle with the hypotenuse as the new vector, and the hypotenuse is always the longest side in right angle triangles.

In my simulation I first get F from the above Lorentz formula, then divide by m, then multiply with the timestep and add it to the velocity of a particle. Finally I scale it to match the magnitude it had before adding to it.

A fusion simulator will not be making calculations on each proton in the simulation. You should be working more along the lines of a plasma physics simulation, no?

https://www.amazon.com/Introduction...olled-Fusion/dp/3319793918/?tag=pfamazon01-20

Well, fusion reactors use very little fuel and there is development on very small fusion reactors, so running simulations on a particle level isn't too far fetched. Also I want to simulate atomic level reactors that might be buildable in the far future. Plus atomic simulations help collaborate results from fluid/gas approximations.

Thanks for the book recommendation! I'll check it out.

 

[Vanadium 50]

What on earth are you trying to do?

Even at the LHC, where bunches weigh picograms they don't simulate every particle.

Further, nobody is seriously talking about pp fusion reactors. It is far too slow to be useful.

 

[Appelros]

What on earth are you trying to do?

Even at the LHC, where bunches weigh picograms they don't simulate every particle.

As I said above I'd like to simulate futuristic reactors that are extremely small, maybe have a hi score of who can make the smallest one. Also if an approximation of one million particles agrees with the particle by particle simulation of a million particles it helps validate the approximation.

Further, nobody is seriously talking about pp fusion reactors. It is far too slow to be useful.

Yes, obviously D-T is the usual reaction. I use protons for calibration.

 

[Vanadium 50]

This seems like you are putting the cart before the horse. Getting a puff of gas that is only a million particles is not simple. At LHC-quality vacuum, it's a volume of about a quarter of a mm on a side.

At homemade-quality vacuum, this would be even smaller.

And you should absolutely not be messing with tritium.

 

[Appelros]

This seems like you are putting the cart before the horse. Getting a puff of gas that is only a million particles is not simple. At LHC-quality vacuum, it's a volume of about a quarter of a mm on a side.

Loop Quantum Gravity has even harder experimental setups and that is still an interesting field of study.

And you should absolutely not be messing with tritium.

Not even in a computer simulation?

 

[Vanadium 50]

Loop Quantum Gravity has even harder experimental setups and that is still an interesting field of study.

There's a non sequitur.

If you are trying to simulate something that can't be built, well, it's your time to waste. But there is a reason that people who need to simulate real devices don't do it this way.

 

[Appelros]

There's a non sequitur.

No.

If you are trying to simulate something that can't be built, well, it's your time to waste. But there is a reason that people who need to simulate real devices don't do it this way.

If you don't see the value in being able to simulate individual particles for fusion research then I don't know what to tell you.

If you have any input on the problems I listed I'd be glad the hear it.

 

[hutchphd]

TL;DR Summary: The magnetic field does no work but translating its force to velocity does, the solution is not elegant and inexact.

Is there a more elegant solution that solves these problems?

There are far too many to detail here. The motion of a proton in a uniform field is exactly soluble.....is it your intent to put further complications into your model that will demand o simulation??
Your most important error in what you have detailed is how you obtain the new position. You need to use the average velocity during the timestep which will be much more accurate. Please work it out for a single step and show (us) why it is a vast improvement. You are correct that rescaling v each timestep is a very unsatisfactory solution.

 

[Mordred]

If you ever done large particle number N-body simulations you would quickly recognize how quickly a PC can get bogged down in calculations.
If you want your simulation to have any useful processing time using averages and currents would be far quicker.

 

[Appelros]

There are far too many to detail here. The motion of a proton in a uniform field is exactly soluble.....is it your intent to put further complications into your model that will demand o simulation??

Yes.

Your most important error in what you have detailed is how you obtain the new position. You need to use the average velocity during the timestep which will be much more accurate. Please work it out for a single step and show (us) why it is a vast improvement. You are correct that rescaling v each timestep is a very unsatisfactory solution.

Ok, if I understand you correctly this is what you mean:
```
v=[1,0,0]
d=[0,0,0.1]
vavg=(v+(v+d))/2 -> [1,0,0.05]
```
Sure it's an improvement but the magnitude still increases. Scaling coincides exactly with the Boris algorithm and I don't think taking an accuracy loss to be more "elegant" is a good idea.

 

[Appelros]

If you ever done large particle number N-body simulations you would quickly recognize how quickly a PC can get bogged down in calculations.
If you want your simulation to have any useful processing time using averages and currents would be far quicker.

Every plasma course I have taken has covered mechanics of individual particles so it seems like a good place to start. Training a neural network on the physics is also something I see that can speed things up a lot.

 

[Mordred]

Those courses should have also included the related current equations using individual particles when your first learning is useful as it's often easier to understand than a field treatment. It quickly becomes far more practical to apply weighted averages in large multiparticle systems.

 

[Vanadium 50]

If you are simulating one individual particle, there is an analytic solution.
If you are simulating a million, a) it is not enough to capture the bulk physics, and b) it is too many for a PC. There are half a trillion particle-particle pairs that need to be calculated. @Modred is absolutely right here.

Pooh-pooh others' experience if you want, but there is a reason that people who do this use voxels and/or small numbers of individual particles.

 

[Mordred]

As you haven't described your fusion reactor simulation you may find it also practical to develop your coordinate system with the beam itself.
This is done typically in cyclotrons the beam becomes the reference frame.

For cyclotrons
Frenet-Serret Frame/coordinates is typically employed.

Scatterings however is another detail that will also quickly bog down. Particularly when you try using Breit Wigner. There are methods to simplify scatterings. Don't know how far your taking the simulation

 

[Appelros]

Those courses should have also included the related current equations using individual particles when your first learning is useful as it's often easier to understand than a field treatment. It quickly becomes far more practical to apply weighted averages in large multiparticle systems.

Yes, approximative models will take up a large part of my project. The particle resolution sims are I guess motivated by my mathematical side seeking an "analytical" solution.

 

[Appelros]

If you are simulating one individual particle, there is an analytic solution.
If you are simulating a million, a) it is not enough to capture the bulk physics, and b) it is too many for a PC. There are half a trillion particle-particle pairs that need to be calculated. @Modred is absolutely right here.

You only have to pair particles that are close to each other.

Pooh-pooh others' experience if you want, but there is a reason that people who do this use voxels and/or small numbers of individual particles.

It seems you are projecting, I'm not doing that, you are the one discrediting my decades of mathematical knowledge and coding experience. I'm sure you are very knowledgeable in many fields but please try to be more polite.

 

[Appelros]

As you haven't described your fusion reactor simulation you may find it also practical to develop your coordinate system with the beam itself.
This is done typically in cyclotrons the beam becomes the reference frame.

For cyclotrons
Frenet-Serret Frame/coordinates is typically employed.

Scatterings however is another detail that will also quickly bog down. Particularly when you try using Breit Wigner. There are methods to simplify scatterings.

Funny you should mention cyclotrons, I just today finished a simulation of one.

That is an interesting coordinate system approach, I'll keep that in mind in my further studies.

Don't know how far your taking the simulation

During COVID there was an app that gamified folding molecules to help vaccine research, I'm thinking of something similar.

 

[Mordred]

One other piece of advise though you may already considered it. Use pointers instead of variables you save on clock cycles and hence processing time.

 

[Appelros]

One other piece of advise though you may already considered it. Use pointers instead of variables you save on clock cycles and hence processing time.

Well I use Julia and it doesn't really expose pointers.

 

[Mordred]

Fair enough I would still look into any method to reduce clock cycles. For example in binary operations a bit shift left is divide by 2 bit shift right multiply by two. Far faster than implying the ALU.

 

[hutchphd]

v=[1,0,0]
d=[0,0,0.1]
vavg=(v+(v+d))/2 -> [1,0,0.05]
```
Sure it's an improvement but the magnitude still increases. Scaling coincides exactly with the Boris algorithm and I don't think taking an accuracy loss to be more "elegant" is a good idea.

Are you calculating in 1D? I haven't a clue what this shows. Sorry

 

[Appelros]

Fair enough I would still look into any method to reduce clock cycles. For example in binary operations a bit shift left is divide by 2 bit shift right multiply by two. Far faster than implying the ALU.

Yea optimization is certainly a major consideration, but first I have to get things working. Premature optimization is a waste, for example I recently scrapped almost all of one of my core files.

 

[Appelros]

Are you calculating in 1D? I haven't a clue what this shows. Sorry

No this is 3D vectors as they are written in Julia, python and matlab have similar notation. v is a vector with x value 1 and y/z 0. Why don't you write your solution? Does it have an error less than dt/2?

 

[hutchphd]

This is an exactly solvable problem found in many textbooks. It is called cyclotron motion. I need not recapitulate it here. Please look it up

 

[Appelros]

This is an exactly solvable problem found in many textbooks. It is called cyclotron motion. I need not recapitulate it here. Please look it up

The thread title clearly states this is about numerical simulations. Euler, midpoint and RK, with scaling, match Boris algorithm at an error of dt/2. You said there was a better way without scaling and hinted at averages, then ignored my reply interpreting your solution and showing it was worse, and now you're condescendingly dismissing the entire topic by handwavingly saying complex computer simulations should be done by analytically solving the smallest parts

 

[hutchphd]

You said you wanted accuracy. Toward that end simulations should be done analytically if they can be.
Your microscopic simulation is not time reversal invariant. That should worry you. You introduce that asymmetry artificially by your second difference technique. (Hand-waving finished).
If you wish to communiucate math, please do it in LateX and not Julia.
Remember who is the supplicant here and do not characterise my attempts to help you. Clearly they were not helpful. They have ended. (*******)

 

[Vanadium 50]

my decades of mathematical knowledge

Then why are you shocked that $\sqrt{x^2 + \Delta y^2} > x$? You are approximating dy by Δy in a region where the approximation does not have the property you want and are shocked to find the outcome doesn't have the property you want either? Of course it doesn't!

As has been said multiple times, if you can use an analytic solution, do so.

 

[Vanadium 50]

Funny you should mention cyclotrons, I just today finished a simulation of one.

Really? And you verified that it gives the correct answer? I suspect that if you ran it assuming no driving voltage you'd still see the particle gain energy, for the same reason you see this system gain energy. If you use the same code, you'll get the same result.