## 3D CuPy-based electric field modeling

- Kuba Anglin
**Posts:**81**Joined:**Fri Dec 19, 2014 4:11 pm**Real name:**Kuba Anglin**Location:**Cambridge, MA

### 3D CuPy-based electric field modeling

Hello everyone,

I thought I would share a project I've been working on. My goal is to help students intuitively understand and predict plasma behavior in fusors.

Here is a short write-up with some relevant graphical outputs. The code is provided at the bottom of the page along with a more comprehensive summary.

The physical accuracy of my program is definitely lacking, as I have not added anything to account for changes in ion densities or other plasma dynamics. However, I feel this approximation good enough for the main purpose of providing an understandable e-field visualization with contour mapping.

Kuba

I thought I would share a project I've been working on. My goal is to help students intuitively understand and predict plasma behavior in fusors.

Here is a short write-up with some relevant graphical outputs. The code is provided at the bottom of the page along with a more comprehensive summary.

The physical accuracy of my program is definitely lacking, as I have not added anything to account for changes in ion densities or other plasma dynamics. However, I feel this approximation good enough for the main purpose of providing an understandable e-field visualization with contour mapping.

Kuba

- Rich Gorski
**Posts:**247**Joined:**Mon Aug 01, 2022 4:34 pm**Real name:**Rich Gorski**Location:**Illinois

### Re: 3D CuPy-based electric field modeling

Nice 3D FD simulation of electric potential and fields around cathode rings. It appears you have the cathode rings at some potential and there is another potential at the four corners of the simulation volume. These four edge potentials (I assume are ground or zero potential) are creating a field around the rings that have a quadrupole component. Did I interpret that right? In a fusor especially those with a spherical chamber the ground or zero potential will be much more continuous (inner surface of a sphere) so the quadrupole field is likely not present. Since you already have the electric field mapped out I assume your next improvement will be the addition of a charged particle trajectory routine using maybe a Runge-Kutta technique (essentially solving Qe = Ma).

Good luck.

Rich G.

Good luck.

Rich G.

- Kuba Anglin
**Posts:**81**Joined:**Fri Dec 19, 2014 4:11 pm**Real name:**Kuba Anglin**Location:**Cambridge, MA

### Re: 3D CuPy-based electric field modeling

That's a great catch that I somehow overlooked haha.

The simulation is set up where there is a spherical anode set at 0 potential. However, I have the actual computation radius set smaller than the anode radius, so the 400x400x400 computation grid never actually contains any of these 0 potential constraints. I did most of my testing with the computation radius set to the anode radius, and only later zoomed in (completely forgetting that this removes the anode from the computation).

The issue is that I want a high resolution view of the cathode, but increasing the grid resolution to account for the larger computation area would lead to obscene computation times even on the best GPU nodes available to me at the PSFC. I've tried to make a dynamic grid resolution function that increases closer to cathode structures, but I failed every time running into so many bugs.

In any case, I should retry these computations with the anode included to see if this is the only thing causing this issue. With a cubic grid size of 400, this will take around 12 hours on my local machine. There are likely a few optimizations I could make to my code to reduce computational cost, but this is definitely a bit outside my area of expertise.

Yes, I do intend to do some particle trajectory simulations. I've previously worked on some matlab simluations which approximated the cathode as a spherical shell. These actually seem to suggest the ideal 1:5 cathode to anode ratio as ideal for fusion. It will certainly be much more accurate using the e-field data generated from a more robust model.

I'll share some media from that project:

model of 2" diameter cathode with 10" diameter anode

matlab simulation results

resulting particle energy distribution

10,000 particle simulation 3D graph

zoomed in

Results from 1" radius cathode (highest efficiency)

========================================

Simulation Summary

========================================

----------------------------------------

Fusion Events:

----------------------------------------

Mean: 4757.58

Standard Deviation: 14.01

Minimum: 4739

Maximum: 4783

Fusion Efficiency: 95.15%

Percent Error: 0.29%

----------------------------------------

Cathode Collisions:

----------------------------------------

Mean: 157.75

Standard Deviation: 19.74

Minimum: 122

Maximum: 188

Percentage of Particles Destroyed by Cathode Collisions: 1.58%

----------------------------------------

Cathode Transparency: 97.05%

----------------------------------------

Simulation Parameters:

----------------------------------------

Simulations: 12

Initial Particle Count: 10000

Simulation Time: 600 ns

Number of Steps per Simulation: 6000

Cathode Radius: 0.0254 m

Anode Radius: 0.1270 m

Cathode to Anode Radius Ratio: 0.2000

Cathode Wire Radius: 0.00025 m

Cathode Voltage: -40000.00 V

Peak Deuteron Energy from Polynomial Fit: 1.38 MeV

========================================

Kuba

The simulation is set up where there is a spherical anode set at 0 potential. However, I have the actual computation radius set smaller than the anode radius, so the 400x400x400 computation grid never actually contains any of these 0 potential constraints. I did most of my testing with the computation radius set to the anode radius, and only later zoomed in (completely forgetting that this removes the anode from the computation).

The issue is that I want a high resolution view of the cathode, but increasing the grid resolution to account for the larger computation area would lead to obscene computation times even on the best GPU nodes available to me at the PSFC. I've tried to make a dynamic grid resolution function that increases closer to cathode structures, but I failed every time running into so many bugs.

In any case, I should retry these computations with the anode included to see if this is the only thing causing this issue. With a cubic grid size of 400, this will take around 12 hours on my local machine. There are likely a few optimizations I could make to my code to reduce computational cost, but this is definitely a bit outside my area of expertise.

Yes, I do intend to do some particle trajectory simulations. I've previously worked on some matlab simluations which approximated the cathode as a spherical shell. These actually seem to suggest the ideal 1:5 cathode to anode ratio as ideal for fusion. It will certainly be much more accurate using the e-field data generated from a more robust model.

I'll share some media from that project:

model of 2" diameter cathode with 10" diameter anode

matlab simulation results

resulting particle energy distribution

10,000 particle simulation 3D graph

zoomed in

Results from 1" radius cathode (highest efficiency)

========================================

Simulation Summary

========================================

----------------------------------------

Fusion Events:

----------------------------------------

Mean: 4757.58

Standard Deviation: 14.01

Minimum: 4739

Maximum: 4783

Fusion Efficiency: 95.15%

Percent Error: 0.29%

----------------------------------------

Cathode Collisions:

----------------------------------------

Mean: 157.75

Standard Deviation: 19.74

Minimum: 122

Maximum: 188

Percentage of Particles Destroyed by Cathode Collisions: 1.58%

----------------------------------------

Cathode Transparency: 97.05%

----------------------------------------

Simulation Parameters:

----------------------------------------

Simulations: 12

Initial Particle Count: 10000

Simulation Time: 600 ns

Number of Steps per Simulation: 6000

Cathode Radius: 0.0254 m

Anode Radius: 0.1270 m

Cathode to Anode Radius Ratio: 0.2000

Cathode Wire Radius: 0.00025 m

Cathode Voltage: -40000.00 V

Peak Deuteron Energy from Polynomial Fit: 1.38 MeV

========================================

Kuba

- Liam David
**Posts:**566**Joined:**Sat Jan 25, 2014 5:30 pm**Real name:**Liam David**Location:**PPPL

### Re: 3D CuPy-based electric field modeling

I'd advise against using Runga-Kutta (RK) integrators for these kinds of simulations, especially ones where particles are long-lived, since generic RK is not symplectic (nor time-reversible) and thus does not conserve energy (e.g. https://www.esaim-m2an.org/articles/m2a ... n0876.html). I would suggest the Boris particle-pushing algorithm or in the absence of magnetic fields, the leapfrog algorithm. These are also much faster then RK, requiring fewer field interpolation steps.

Questions and requested clarifications about the particle simulation:

How do you compute the fusion rate? Beam-beam? Beam-gas? I advise strongly against portraying glow-discharge fusors in the typical "converging ion beam fusion" picture.

How is "Fusion Efficiency" defined? Some 4000 fusion events for just 10k particles is rather high.

Are particles initialized uniformly throughout the volume?

Why 600ns duration?

The aliasing of the thin grid wires by the 400^3 grid is probably significant: how sure are you about the cathode collision rate because of this? Or does the particle sim, which I recognize approximates the potential as exactly 1/r, calculate collisions differently?

Regarding the simulation speed: 6000 steps for 10k particles taking 12 hours seems excessive, especially if there's no grid. GPU or CPU? Vectorized?

I suggest looking at:

viewtopic.php?t=14157

viewtopic.php?t=14642

viewtopic.php?p=96483#p96483

Questions and requested clarifications about the particle simulation:

Typo? Also your energy scale bar on the deuteron distribution function goes to 160 keV? In a purely electrostatic simulation without particle collisions and a -40 kV cathode, that will be the maximum possible energy. The general distribution shape is roughly consistent with collisionless simulations I've done.Cathode Voltage: -40000.00 V

Peak Deuteron Energy from Polynomial Fit: 1.38 MeV

How do you compute the fusion rate? Beam-beam? Beam-gas? I advise strongly against portraying glow-discharge fusors in the typical "converging ion beam fusion" picture.

How is "Fusion Efficiency" defined? Some 4000 fusion events for just 10k particles is rather high.

Are particles initialized uniformly throughout the volume?

Why 600ns duration?

The aliasing of the thin grid wires by the 400^3 grid is probably significant: how sure are you about the cathode collision rate because of this? Or does the particle sim, which I recognize approximates the potential as exactly 1/r, calculate collisions differently?

Regarding the simulation speed: 6000 steps for 10k particles taking 12 hours seems excessive, especially if there's no grid. GPU or CPU? Vectorized?

I suggest looking at:

viewtopic.php?t=14157

viewtopic.php?t=14642

viewtopic.php?p=96483#p96483

- Kuba Anglin
**Posts:**81**Joined:**Fri Dec 19, 2014 4:11 pm**Real name:**Kuba Anglin**Location:**Cambridge, MA

### Re: 3D CuPy-based electric field modeling

Thanks for taking a look at this Liam, I've been checking out your work on computational modeling over the past few months and find it incredibly impressive and thorough. I've only recently started learning to make computational models myself, so I really appreciate the feedback.

Also apologies for the confusion, the matlab particle sim is a completely different program than the python electric field model. The matlab program does not use a grid, e-field computations are absolute without interpolation. The python program uses a set resolution where the cathode is mapped onto the grid points.

The matlab particle simulation was my first project in computational physics, and it definitely needs a lot of work. I intended to make a monte-carlo sim to determine the best cathode anode ratio. For this project, I used parallel CPU processing in my 12 available threads. To give a brief summary of its operation:

-The cathode and anode are approximated by spherical shells

-particles are initially spawned at random locations between the cathode and anode with random initial velocities in the order of 1e5 m/s

-each loop, e-fields are calculated for each particle's position, force is applied, and the particle velocity is updated.

-fusion events occur when two particles are both close enough to each other and have a high enough relative velocity. These thresholds were adjusted until generated error data was meaningful.

The results and values are somewhat meaningless, as the initial conditions and computation variables are pretty arbitrary and not at all physically accurate. This was more of a personal exercise than anything else.

Anyway, the project I am more concerned with is my current python electric field model for mapping contours of various cathode configurations in my first post. I'll check out the Boris particle-pushing algorithm for making a new PIC sim using the exported e-field data from my program. I've added and removed more complex computations such as using the Poisson equation instead of the Laplace equation, determining ideal grid resolution through the Debye length, modeling sheath formation around the cathode rings, but all these basically made my program unusable due to computational cost. This is why I've settled with just using the Laplace equation and solving with the FDM. Not perfect, but I feel it's a good enough approximation for generating cool-looking visuals.

I'm currently designing a research-grade fusor for experimental education. I plan on using an 18" spherical vacuum vessel with options for installing a 10" anode 2" cathode system or a larger cathode without a dedicated anode depending on stability requirements. Various plasma diagnostics will be used such as Langmuir probes and spectroscopic analysis. Later additions might include electron filaments or an ion injector. I'm designing a control console to safely house HV electronics, record all measured data, and regulate pressure and power. The field modeling tool will hopefully give students a more intuitive understanding of how fusors work.

Kuba

Also apologies for the confusion, the matlab particle sim is a completely different program than the python electric field model. The matlab program does not use a grid, e-field computations are absolute without interpolation. The python program uses a set resolution where the cathode is mapped onto the grid points.

The matlab particle simulation was my first project in computational physics, and it definitely needs a lot of work. I intended to make a monte-carlo sim to determine the best cathode anode ratio. For this project, I used parallel CPU processing in my 12 available threads. To give a brief summary of its operation:

-The cathode and anode are approximated by spherical shells

-particles are initially spawned at random locations between the cathode and anode with random initial velocities in the order of 1e5 m/s

-each loop, e-fields are calculated for each particle's position, force is applied, and the particle velocity is updated.

-fusion events occur when two particles are both close enough to each other and have a high enough relative velocity. These thresholds were adjusted until generated error data was meaningful.

The results and values are somewhat meaningless, as the initial conditions and computation variables are pretty arbitrary and not at all physically accurate. This was more of a personal exercise than anything else.

Anyway, the project I am more concerned with is my current python electric field model for mapping contours of various cathode configurations in my first post. I'll check out the Boris particle-pushing algorithm for making a new PIC sim using the exported e-field data from my program. I've added and removed more complex computations such as using the Poisson equation instead of the Laplace equation, determining ideal grid resolution through the Debye length, modeling sheath formation around the cathode rings, but all these basically made my program unusable due to computational cost. This is why I've settled with just using the Laplace equation and solving with the FDM. Not perfect, but I feel it's a good enough approximation for generating cool-looking visuals.

I'm currently designing a research-grade fusor for experimental education. I plan on using an 18" spherical vacuum vessel with options for installing a 10" anode 2" cathode system or a larger cathode without a dedicated anode depending on stability requirements. Various plasma diagnostics will be used such as Langmuir probes and spectroscopic analysis. Later additions might include electron filaments or an ion injector. I'm designing a control console to safely house HV electronics, record all measured data, and regulate pressure and power. The field modeling tool will hopefully give students a more intuitive understanding of how fusors work.

Kuba

- Rich Gorski
**Posts:**247**Joined:**Mon Aug 01, 2022 4:34 pm**Real name:**Rich Gorski**Location:**Illinois

### Re: 3D CuPy-based electric field modeling

I wouldn’t worry about the energy conservation error mode of 4th order RK methods since the flight path in a fusor is so short ie… a few cm MFP length and the particles short life span before collision at fusor pressures. We are not talking about a harmonic oscillator where preservation of energy is important after many oscillations. 4th order RK methods are known to be highly accurate in beam transport situations where only one transit through the acceleration system is required as in transport of an electron from cathode to the phosphor screen of a old fashioned CRT.

Unless you are ignoring collisions in a fusor and just want to allow a particle to circulate indefinitely I wouldn't be so concerned. If collisions are allowed at normal fusor pressures (~10mTorr) and then stopped in the simulation I would just use which ever method, RK, Boris... is easier to code.

Just my limited opinion.

Rich G.

Unless you are ignoring collisions in a fusor and just want to allow a particle to circulate indefinitely I wouldn't be so concerned. If collisions are allowed at normal fusor pressures (~10mTorr) and then stopped in the simulation I would just use which ever method, RK, Boris... is easier to code.

Just my limited opinion.

Rich G.

- Liam David
**Posts:**566**Joined:**Sat Jan 25, 2014 5:30 pm**Real name:**Liam David**Location:**PPPL

### Re: 3D CuPy-based electric field modeling

Sure, in single-pass systems RK is just fine most of the time, but in e.g. a collisionless recirculating fusor, long-term energy conservation is important to avoid artificial losses to the chamber wall and cathode. In a highly collisional (against neutrals only) simulation energy conservation is not quite as important, but RK is more complicated than Boris, requiring more field evaluations, and the former just isn't used in plasma simulations.

I'm glad you're getting into computational physics and I'm sure the fidelity of your models will improve with time. I know how you feel in trying to add additional physics to a simulation only to run into computational limitations or stability problems: I've encountered the same issues in my codes. It sounds like you have quite the work cut out for you preparing the educational fusor. Good luck!

I'm glad you're getting into computational physics and I'm sure the fidelity of your models will improve with time. I know how you feel in trying to add additional physics to a simulation only to run into computational limitations or stability problems: I've encountered the same issues in my codes. It sounds like you have quite the work cut out for you preparing the educational fusor. Good luck!

- Rich Gorski
**Posts:**247**Joined:**Mon Aug 01, 2022 4:34 pm**Real name:**Rich Gorski**Location:**Illinois

### Re: 3D CuPy-based electric field modeling

Liam,

Agreed. Is that a real thing... recirculating, collisionless fusor? It would have to run at very low pressure unlike most fusors described here. My device is probably one of the lower vacuum devices presented here running with ion sources at 1 micron.

Rich G.

Agreed. Is that a real thing... recirculating, collisionless fusor? It would have to run at very low pressure unlike most fusors described here. My device is probably one of the lower vacuum devices presented here running with ion sources at 1 micron.

Rich G.

- Liam David
**Posts:**566**Joined:**Sat Jan 25, 2014 5:30 pm**Real name:**Liam David**Location:**PPPL

### Re: 3D CuPy-based electric field modeling

It's somewhat a real thing, e.g. the SIGFE device in some operating modes or the earlier ion-gunned Hirsch designs, but there certainly aren't any examples here. I was referring only to Kuba's simulation.