Python Simulation

This area is for ADVANCED theory and discussions only. If you just joined fusor.net, chances are this is NOT the area for you to be posting.
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Python Simulation

Post by Nicolas Krause »

I've decided to break out the work I'm doing into its own thread. As of right now I've completed porting my initial code over to python, which means that I have a 2d representation of the electric potential in a fusor. My next step is to get a simple particle in cell simulation going with a small number of particles.
2dVoltagepotential.png
Dan Knapp
Posts: 402
Joined: Wed Aug 06, 2008 8:34 am
Real name: Dan Knapp

Re: Python Simulation

Post by Dan Knapp »

You should be aware that a 2D representation of spherical geometry is OK for plotting potential field lines, but 2D will not give you a valid PIC simulation of a spherical system. You have to use full 3D to get valid PIC simulation of spherical geometry. A 2D simulation is only valid when any slice through the system gives the same 2D geometry, e.g. an infinitely long cylindrical geometry.
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

I'm aware it won't be an accurate representation of a spherical system, but it seemed a good way to get my feet wet seeing as this is the first time I've written a numerical simulation. In particular I was inspired by the following work done by a physicist by the name of Matt Lilley, who used a 2d simulation and some nonlinear dynamics to investigate how star mode might form. I'm hoping to first replicate his results, and then see if I can expand upon them. Either by extending the simulation to 3 dimensions or by pursuing some other avenue that appears along the way.
Dan Knapp
Posts: 402
Joined: Wed Aug 06, 2008 8:34 am
Real name: Dan Knapp

Re: Python Simulation

Post by Dan Knapp »

Sounds like a good plan. I just wanted to be sure you were aware of the limitations of 2D. Good luck in your endeavors!
Harald_Consul
Posts: 94
Joined: Sat Oct 20, 2018 6:01 am
Real name: Harald Consul

Re: Python Simulation

Post by Harald_Consul »

If you are not planning a scientific publication of your "electric potential in a Farnsworth fusor simulation", it would be interesting to see the code.

I am not too experienced in Python, but I do a lot of R programming (including simulation tasks, also). Usually a numerical code is easy to understand.
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

Hi Harald,

Definitely not planning to publish anything, but happy to share the code I have so far, it's not very complicated at all. How would you like me to send the file?
Harald_Consul
Posts: 94
Joined: Sat Oct 20, 2018 6:01 am
Real name: Harald Consul

Re: Python Simulation

Post by Harald_Consul »

Just plug it into the thread using the "</>" Button or the BB-code tags

Code: Select all

[code]
PRG listing
[/code]
, so that anyone interested can view it.
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

Here you go! The code uses numpy and matplotlib, two freely available python packages. A matrix called fusor is defined and a finite difference approximation is performed until the matrix converges on a solution. The program is based on some matlab code I found on a professor's website, he has written a nice summary of the general method in the following document.

Code: Select all

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

#fusor is a 2d grid of ~10000 points
Nx=101
Ny=101
fusor = np.zeros((Ny,Nx))

#setup fusor initial voltage
Vsource= -45000

#set the resolution in x and y coordinates
res = 0.0015
hx= res
hy= res
#precalculate kx and ky
kx= (hx**2)/(2*(hx**2+hy**2))
ky= (hy**2)/(2*(hx**2+hy**2))
#initial conditions for our while loop
tol = 0.001
dsum=1
iter=0

while dsum > tol:
	#sum our squared fusor grid
	square1=np.square(fusor)
	sum1=square1.sum()
	print(sum1)

	#iterate through rows
	for ny in range(2,100):

		#iterate through columns
		for nx in range(2,100):
			#source voltages are replaced every loop
			fusor[50,40]=Vsource
			fusor[50,60]=Vsource
			fusor[40,47]=Vsource
			fusor[60,47]=Vsource
			fusor[40,55]=Vsource
			fusor[60,55]=Vsource
			#perform a finite difference approximation
			fusor[ny,nx]=ky*(fusor[ny,nx+1]+fusor[ny,nx-1])+kx*(fusor[ny+1,nx]+fusor[ny-1,nx])
		
			

	#sum our fusor grid again
	square2=np.square(fusor)
	sum2=square2.sum()
	print(sum2)

	#check the absolute value difference between our before and after grid
	dsum = sum2-sum1
	dsum = abs(dsum)
	print(dsum)

	#count iterations
	iter = iter+1


#output fusor grid using a plot
print(iter)
plt.imshow(fusor)
plt.colorbar()
plt.show()

#electric potential has been calclated, so now we need the E-Field, so we can find the force on our charged particles
#have to integrate that somehow, how we gonna do that in python?
#place particles
#move particles
ian_krase
Posts: 636
Joined: Mon Nov 28, 2016 2:48 am
Real name: Ian Krase

Re: Python Simulation

Post by ian_krase »

For those of you who are new to Python:

If you want to use Python on Windows, by far the best way is to install the newest version of Miniconda and then open an Anaconda Command Prompt and type the following:

Code: Select all

conda install ipython numpy scipy matplotlib 
and then follow the prompts.

This works much more cleanly than most other methods to install the very complicated Numpy backend.
Harald_Consul
Posts: 94
Joined: Sat Oct 20, 2018 6:01 am
Real name: Harald Consul

Re: Python Simulation

Post by Harald_Consul »

Another beautiful possibility for numerical programming is R, which can be downloaded at https://cran.r-project.org/ .

The advantages of R are
  • Data formats are all the same through out all R
  • Missings are handled all the same through out all R
  • R installes out of the box an all operating systems
  • R code runs platform independent on all OS
Wyatt Sheffield
Posts: 1
Joined: Thu Aug 29, 2019 11:27 pm
Real name: Wyatt Sheffield
Contact:

Re: Python Simulation

Post by Wyatt Sheffield »

This simulation, while no doubt quite accurate takes a loooong time to compute. If we are wanting to simulate individual particles, I think it might be more effective to calculate the field and forces at the position of each particle only. Of course, if there are n charged particles in the simulation it will run in O(n^2) time - that is, if it takes 1 second to calculate the forces and fields on 1000 particles, it would take 4 seconds for 2000, 9 seconds for 3000, etc. However, given the length of time it took for the given python script to halt I still think this may be a viable alternative.
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

Hi Wyatt,

The current simulation is incomplete, the plan is to calculate a very precise voltage field using a finite difference approximation (which takes a long time) and then use that voltage field as a reference for the rest of the simulation. That way you can do one large chunk of computation and not have to repeat the computation over and over, which I was hoping would help when scaling to large numbers of particles. I've done a bit more work on this over the summer, but while the code compiles I have to tweak some parameters to get it running properly. If you want to speed up the simulation I'd suggest the following.
  • Dial down the sensitivity
  • Use the existing code in a Jupyter notebook:
    This way a single computation could be completed to produce an accurate field, and then you could compute particle paths and collisions on their own, without having to redo the field approximation each and every time.
DR_DANIEL_GOODMAN
Posts: 7
Joined: Wed Dec 27, 2017 11:43 pm
Real name: Daniel Goodman

Re: Python Simulation

Post by DR_DANIEL_GOODMAN »

To all on the simulation thread here:

Purely out of curiosity, I did a search for fuso simulation software, and got this:

https://www.bing.com/search?pc=SUWI&for ... n+software

That then had two apparently relevant SourceForge links:

https://sourceforge.net/projects/fusortools/

And:

https://sourceforge.net/p/fusortools/me ... f/default/

I'll try to see if I can find more; also, anyone here increasing this ever done such simulation in MATLAB at all, I'd just wondered, I'd done a good deal of numerical methods before my disability at various points, and had given thought to trying to do such simukations, I had to do a good deal of such software searching, both freeware/open-source, as well as commercial, at various points, I'll try to send more if I can find it at all, if germane, seemingly, hope the above was of some interest, at least, if nothing ekse, of course.

Many thanks,

Dr. Daniel Goodman.

P.S. I'm also interested in any hobbyist group efforts for chairs on Long Island, as I have an interest in possibly integrating fusors into volunteer STEM tutoring I'm seeking to be allowed to do, if at all possibke, as well, if anyone might know if such groups and/or projects near my wife and myself, many thanks once again, I'd be interested in thoughts about the sites I'd sent in above, and the search also, whenever convenient, no rush, of course, please pardon any typos, I'm on an Android tablet at the moment with not the greatest keyboard and spell checker, I'm afraid.
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

I figured it was time for an update on this little project. Over the summer I monkeyed around a bit more with the code, it correctly simulates a 2d electrical field from a grid, and I can place a particle in the field and watch it move (if you'd like the completed code I can post it here). However, calculating the field takes a decent amount of time, and one of the core reasons for making this simulation is to be able to test different grid shapes. Already I can get a calculation that takes a few days to complete if I put in a decent level of resolution, if I extend the simulation to 3d I'm quite concerned it'll take too long to be of real use as a tool for testing. However, if I were to build a small computational cluster I'd be able to achieve a significant speedup. So that's what I plan on doing! To start I'm going to port my initial code over to a new language called Julia. It's used extensively in scientific computation and has good support for parallelization. Since I have zero clue at the moment how to program for parallelization I've constructed a small low cost cluster from 4 Raspberry Pi's and a used network switch. I still have to configure the devices and get them setup. I may start a new thread once I actually get around to that!
IMG_1395[1].JPG
User avatar
Nathan Marshall
Posts: 53
Joined: Wed May 08, 2019 8:13 pm
Real name: Nathan Marshall

Re: Python Simulation

Post by Nathan Marshall »

I wrote a numerical solver in Python for the 3D Laplace equation with the relaxation method. Here I set up a two ring grid with a ring in the X-plane and a ring in the Z-plane. The two images are slices of the potential on the X and Y planes so we can see one complete ring in one image and the other just has four points from the two rings. I tried to visualize the E-field by taking the gradient and plotting vector fields/streamplots but it was a real mess... Visualization for the 3D case is much more difficult. I did 10,000 iterations to get these images. I used the numba just-in-time compiler which tremendously speeds up the triple nested for loop needed for 3D. The simulation took about a minute on my laptop.
x_plane.png
y_plane.png
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

Very cool Nathan, as an aside my dad did his undergrad at Kansas State! Yes with 3d I think visualisation is quite difficult, I'm not 100% sure on the tack I'll take yet. From my research the easiest way to do it seems to be to take "slices" of the field at different heights, so you calculate the field say a cm away from the wire, then 2cm, then 3cm etc... and you can view them discretely. Later on I might move to a more complicated method of visualisation once I understand what I'm doing.
User avatar
Nathan Marshall
Posts: 53
Joined: Wed May 08, 2019 8:13 pm
Real name: Nathan Marshall

Re: Python Simulation

Post by Nathan Marshall »

It has been a while since the last post on this thread. Tonight I had some inspiration to continue playing around with fusor simulations and I think I'm getting some interesting results now. I wrote a program that I have jokingly named "SIMION Lite". It takes the 3D potentials I calculated before and converts them to an electric field by taking the gradient. Then I create a deuteron with some initial position and zero velocity. By determining which grid box of my electric field its current position falls within, I can find the force on the deuteron and consequently its acceleration. I then increment its velocity and position with this information. This continues until the particle runs into a boundary like the grid or I reach a predetermined max number of iterations. Here is a plot of 20 different ion tracks using this program:
Two_Ring_Grid_Ion_Paths.png
It seems that the cross of ion jets observed in two ring fusors is clearly showing!

Nicolas, that is really neat that your dad did an undergrad at K-State too! It is a small world.
User avatar
Richard Hull
Moderator
Posts: 14991
Joined: Fri Jun 15, 2001 9:44 am
Real name: Richard Hull

Re: Python Simulation

Post by Richard Hull »

Nice bit of work there. One assumes an infinite mean free path within the field in an ideal vacuum with zero collisions transferring energy to other particles and thereby, having our deuteron becoming a fast neutral. What you show is the absolute ideal, which we in our untrained wisdom back in 1997 believed of the Farnsworth device.... Virtual infinitely looping deuterons.
Now, after many hard won lessons, wakeup calls and being slapped around in the real world...........

Richard Hull
Progress may have been a good thing once, but it just went on too long. - Yogi Berra
Fusion is the energy of the future....and it always will be
The more complex the idea put forward by the poor amateur, the more likely it will never see embodiment
User avatar
Nathan Marshall
Posts: 53
Joined: Wed May 08, 2019 8:13 pm
Real name: Nathan Marshall

Re: Python Simulation

Post by Nathan Marshall »

Thanks Richard! And yes, this is of course an extremely idealized model of a fusor that ignores many factors known to be present in real fusors. Reality is always a rude awakening for the daydreaming theorist.
User avatar
Richard Hull
Moderator
Posts: 14991
Joined: Fri Jun 15, 2001 9:44 am
Real name: Richard Hull

Re: Python Simulation

Post by Richard Hull »

I remain very impressed with the work of others to present to us mathematical simulations in graphical form. I long held, in my mind, the above simplistic example of ion (deuteron) paths from my first learning about the fusor until I actually did fusion in the late 90's. We have long discussed, in many postings, the multitude of flies in the ointment within the fusor environment in which fusion is miraculously taking place in spite of the flies and the many adverse goings-on to doing fusion in the purest sense of the physics involved.

Richard Hull
Progress may have been a good thing once, but it just went on too long. - Yogi Berra
Fusion is the energy of the future....and it always will be
The more complex the idea put forward by the poor amateur, the more likely it will never see embodiment
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

Very nice simulation Nathan. This poster is about the only academic reference I've been able to find that analyzes the possible orbits in a Fusor. It has a very nice perspective from non-linear dynamics about why things like star mode exist in the Fusor.
User avatar
Liam David
Posts: 518
Joined: Sat Jan 25, 2014 5:30 pm
Real name: Liam David
Location: PPPL

Re: Python Simulation

Post by Liam David »

Nicholas, Nathan, your efforts in have inspired me to dig out some simulations that I shelved a year or two ago. With a few more semesters of physics and math under my belt, as well as a functioning fusor (unlike last time), I've decided to take another stab at ion optics and space charge effects.

I've opted to go for a 2D particle-in-cell method as opposed to a direct sum or tree algorithm for a few reasons. It avoids divergences in the forces between particles, allows me to simulate millions of particles at O(N) efficiency, and makes it easier to determine bulk plasma properties. I started on a 3D version but decided to downscale by one dimension due to the matrix sizes needed to solve for the potential. A 3D model discretized into 100^3 cells needs a 1e6 x 1e6 matrix, which even as a sparse matrix takes way more RAM than I have. Sizes up to 70^3 are possible with 32GB, but that's pretty coarse. In 2D, I can easily run a 300^2 simulation with only a few GB.

The potential solver inverts the discrete Laplacian matrix to solve grad^2 (phi) = rho/epsilon0 in the form of a normal matrix equation Ax=b. I take into account boundary conditions, i.e. the grid and chamber, by setting the appropriate cells in b=rho to the boundary values and the corresponding rows of A=grad^2 to the identity. That way, when I compute x = (A^-1)b, I get back the boundary conditions plus the vacuum potential and potential from the particles.

The macroparticles (each represents ~3000 real particles) are binned using a square kernel the size of a grid cell, meaning each can contribute to the properties of up to four cells. This density matrix, with boundary values replaced, is the b in the above equation. In addition to the density, I find the mean velocity of particles in each cell, and subsequently the estimated fusion rate by evaluating the cross-section and f = n1*n2*<sigma*v>*volume. The fusion rate is still very much up-in-the-air, the main uncertainty being whether I should scale my 2D densities to 3D by just adding a 1/length, or by raising it to the 3/2. I could argue either way so perhaps someone has some insights. Particles that hit the walls or cathode are removed.

I found this database of cross-sections for all kinds of particle interactions: https://www-amdis.iaea.org/ALADDIN/collision.html. I downloaded data for three reactions, assuming that values for deuterium are roughly the same as those for hydrogen:
Single charge-exchange: D,D+ --> D+,D
Single ionization: D+,D --> D+,D+,e (discard electron, for now)
Double ionization: D,D --> D+,D+,e
Each time step, the simulation calculates the distance the particles traveled, computes the interaction probability from the cross-sections and density, and adds both ions and fast neutrals when necessary. The particles are given velocities based on the standard elastic collision equations, plus some randomness in what you could call the impact parameter. Spoiler alert, at 10mtorr the mean free path is much smaller than the device size, meaning lots of interactions occur. The vast majority of ions disappear after traveling just 10cm, and the peak fast neutral densities after a couple hundred nanoseconds are comparable to the peak ion densities.

The timestep is 0.1ns such that ions travel across few cells even at maximum speed (~2e6m/s at 50kV). Each step takes ~1 second to compute, increasing somewhat as particles are added, but even so, simulating 1us takes several hours and equilibrium is not achieved.

Space charge from the ions causes a deviation from the vacuum potential of ~2V at ~1us, so negligible. I'd have to run it for a lot longer for this to build up, I'm sure.

Up until today, my initial condition has been a uniform distribution of ions. Currently I am using two, 1mA ion guns.

While I've been focusing on linear cathodes, and specifically ion gunned designs, I decided to model the standard spherical fusor as well.
star mode.PNG
This simulation used 1e6 particles, a 200^2 grid, and ignored charge exchange, neutrals, and anything other than electrostatic interactions, simply because I had not implemented those yet: https://youtu.be/-QoMSFVQbX0. The maximum central density is on the order of 1e8/cm^2, which matches up pretty well with the thesis "Improved Lifetimes and Synchronization Behavior in Multi-grid Inertial Electrostatic Confinement Fusion Devices", which is too large to attach. The star mode also looks pretty nice.

Some more images, specifically of a dual-gunned linear cathode with charge exchange and particle interactions:

Bulk plasma speed for a linear cathode.
Bulk plasma speed for a linear cathode.

Fusion rate over time, assuming ^(3/2) density scaling.
Fusion rate over time, assuming ^(3/2) density scaling.

Ion energy spectrum. x-axis supposed to be labeled keV.
Ion energy spectrum. x-axis supposed to be labeled keV.

I'll post more as I make progress. Perhaps it's worth collaborating to produce a better model.
User avatar
Nicolas Krause
Posts: 230
Joined: Fri Sep 30, 2016 7:36 pm
Real name: Nicolas Krause
Location: Canada
Contact:

Re: Python Simulation

Post by Nicolas Krause »

Great work Liam! I assume you wrote this in Python? Collaboration on a simulation would be fantastic, I've got a ton of ideas that I simply haven't had the time to explore with school. As you've noted with the scaling problems for 3 dimensions, I'd started moving towards some sort of parallelization. I have 4 Raspberry Pi's and a network switch and was hoping to setup a toy beowulf cluster to learn about these things but I've been working on the physical side of the fusor when I've had spare time recently. The other alternative to speed up the simulation is to purchase some computation time on something like Amazon Web Services, I understand it's not too expensive, but I have no idea how to structure code for that. There seem to be a ton of people offering compute time with the rise of AI though.

I think the other big challenge that I simply haven't tackled because I don't know how, is the electron waves that can occur in plasma. I've read Vlasov's original paper, which has some fearsome mathematics in it, but it's pretty apparent that the plasma is a huge strange soup and it doesn't behave the way my simplified model predicts. I hope to use the simulation to provide some clarity about how the fusor actually operates and from there to guide areas of investigation to increase the rate. I'm currently reading some of Irving Langmuir's original papers on electrons in plasma in the hope that they'll be useful.

I have a couple of Pluto notebooks (like jupyter notebooks but for the Julia language) with some more complicated grid shapes that I put together over christmas, but beyond that I think everything I've done is visible here on fusor.net. I suppose the way to collaborate would be to setup a meeting and sort of sketch out some initial ideas?
User avatar
Richard Hull
Moderator
Posts: 14991
Joined: Fri Jun 15, 2001 9:44 am
Real name: Richard Hull

Re: Python Simulation

Post by Richard Hull »

Nicolas, Welcome to the reality club. Plasmas do not obey any real calculable rules until you create a rigidly controlled environment of great simplicity and even then there are gotchas, unexpected, unexplained that throw all of the best mathematical work into a cocked hat. In the case of the fusor, it is a real mess. Wall loading and unloading and the velocity space fusion really mess the thing up. And yet this is just a spherical ball with a central electrode. Seems simple, but it isn't. Precise calculations fail due to the gotchas. Yes you can figure on a lot of generalizations in such a simple system, but high precision and close predictions via mathematics as to fusion, fail.....Even at the hyper advanced amateur level found here. Still, it an interesting and challenging endeavor. Keep at it.

Richard Hull
Progress may have been a good thing once, but it just went on too long. - Yogi Berra
Fusion is the energy of the future....and it always will be
The more complex the idea put forward by the poor amateur, the more likely it will never see embodiment
User avatar
Liam David
Posts: 518
Joined: Sat Jan 25, 2014 5:30 pm
Real name: Liam David
Location: PPPL

Re: Python Simulation

Post by Liam David »

I've been using Matlab for all my simulations. The vector-matrix operations built into the syntax make it comparatively easy to vectorize the code, and it's the language I'm most familiar with because of other work. I've done some stuff in Python, though for fast simulations it isn't ideal, and I've never used Julia but it seems similar to Matlab. The direction I'm leaning is to write overhead processing, visualization, post processing, etc... in a high level language like Matlab, Julia, or Python, but write the actual integrators/solvers in C++. I had some success running CUDA code on an NVIDIA GPU for another project, which could help speed everything up, but my experience there is limited.

In my current simulations I wanted to add electrons, but the particle-in-cell method runs into a problem. My timestep is already 0.1ns, and since electrons are so light, accurately simulating their trajectories would mean a picosecond-scale timestep. They're also relativistic at 50kV, so that adds a layer of complexity to the integrators. I'm not thrilled about waiting years for a few microseconds of total simulation time, so like you, I think the ideal next step is to look at the Vlasov-Maxwell equations. According to this paper I've just started reading (https://www-m16.ma.tum.de/foswiki/pub/M ... -Notes.pdf), the Vlasov equation is fully relativistic and apparently works for non-thermalized plasmas. This paper also goes into some detail (http://www.phys.nthu.edu.tw/~thschang/notes/PP06.pdf). They're both really heavy on the math, and naturally don't include lots of the processes we're interested in like charge exchange, impact ionization, etc... so there's a lot of work ahead.

I think it would be great to discuss general directions and ideas in a meeting. I've got finals till next week but should be good to go after those.
Post Reply

Return to “Advanced Technical Discussion Area”