Dylan Langone

Spring 2018: Equilibria and Stability in Magnetic Confinement Fusion

1/25/2018

• I met with Professor Dorland last month before break and he recommended two theses to read as background:

• Up-Down Asymmetric Tokamaks (Justin Ball): https://arxiv.org/pdf/1611.06903.pdf; I'll start with this one and move on to the other in due course.

• The Zero-Turbulence Manifold in Fusion Plasmas (Edmund Highcock): https://arxiv.org/pdf/1207.4419.pdf

• Quick Summary: Turbulence in plasmas can be exploited when magnetic flux surfaces are not mirror symmetric about the midplane of the torus to cause the plasma to spin and become more stable.  More specifically, up-down asymmetry creates turbulent transport from one surface to the other.

1/31/2018

• I set up another meeting with Professor Dorland to specify the goals of the research project

• I have read Ball's thesis until the presentation of Grad-Shafranov.  I'll try to tackle this soon.

2/2/2018

    From Ball's thesis, the main problem in sustaining a fusion reaction is the turbulence that results from the high temperature of the plasma.  Turbulence leads to particles colliding with the walls of the reactor or otherwise reducing the temperature gradient needed to sustain a reaction.  As a result, more energy is put into creating a reaction than energy generated from the reaction itself.  Without at least breaking even, fusion reactors are not viable sources of energy.  The efficiency record is held by the Joint European Torus (JET) experiment, with a Q of 0.67.  Ball's solution to the turbulence problem involves designing a reactor such that the torus is asymmetric about the midplane.  This configuration enables what is known as "intrinsic" rotation, reducing the energy needed to keep the plasma hot enough for fusion.  

2/3/2018

    I met with Professor Dorland today.  We started by discussing the mathematics of vector spaces, function spaces, and the Fourier series for discrete points.  We then discussed more about the project we are going to undertake.  Using a code, we are going to model a plasma for a given reactor geometry and then take infinitesimal perturbations of parts of the plasma to see how they behave at some point in the future.  Given the laws of magnetohydrodynamics (MHD), the displacement of the plasma from its initial condition can get progressively larger in time, in which case the fusion reaction would be more difficult to sustain.  On the other hand, the plasma displacement can dissipate with time like a damped oscillation.  I will study the mechanics the Fourier Series so as to ground my intuition.  

2/6/2018

    I watched some lecture material on the Fourier Series to better ground my intuition.  This technique for decomposing a period function is useful in solving partial differential equations, such as the ones found in MHD that govern the behavior of a plasma system.  Here are my notes from the previous two days.  I am working on formatting mathematics with Google Sites. 

2/9/2018

    I met with Professor Dorland for two hours today.  We discussed complications associated with Fourier analysis.  In particular, he addressed the Gibbs phenomenon, which is when partial sums of the Fourier Series overshoot a piecewise continuous periodic function.  We also talked about aliasing, which is when a higher harmonic is observed even though it is a linear combination of the independent harmonics.  Then we moved on to the details of the project.  I need to use a VPN provided by UMD's OIT to gain remote access to the computer that runs the plasma code.  I also have a paper to read by Edmund Highcock on zero-turbulence manifolds in toroidal plasmas (download link below).  

2/11/2018

    I needed to download and install UMD's VPN in order to establish a secure connection to the computer that runs the MHD code.  I successfully installed the VPN software and I am currently connect to it.  I'll turn it off for now and get started on Highcock's paper.

Main Idea: reducing the magnetic field pitch over inverse aspect ratio control parameter increases the temperature gradient that can be reached before turbulence kicks in.

2/14/2018

    I met with Professor Dorland today for an hour.  We discussed Highcock's paper and the methods he used to calculate the zero-turbulence manifold.  He used several million core hours just to use the MHD code by varying three parameters.  There are about 20 control parameters and performing a search with only two or three of them is computationally very expensive.  This is why we need to be judicious in choosing which control parameters to vary and by how much.  Then we started on some of the theory behind MHD, starting with relativistic electrodynamics.  It turns out that it is impossible to tell whether or not a particle feels a force due to either the electric or magnetic fields alone.  The two must give the same result, implying that they are indistinguishable from each other.  The same applies to gravity, since the math is exactly the same.  There is an analog for the magnetic field for gravity.  He gave me a reading from Griffiths and some problems with electric and magnetic fields.  I'll post key the results later.

2/17/2017

    Today I met with Professor Dorland for 3 hours.  He gave me access to the Kendall computer that performs the MHD calculations.  Then he showed me how to perform a run and analyze the input and output files.  I don't remember everything but he will send me instructions on how to do the analysis.  I'll post some results from runs when I have a better idea of what I am doing in a few days.  

  , in addition to the helical path there is a drift velocity   .   In this case, the drift velocity is in the   direction.  

2/21/2017

    Since I have been given access to the Kendall computer, I have been brushing up on my Linux skills.  I successfully copied a simple C++ program from my account on the CMS cluster to my account on the Kendall computer by using the secure copy (scp) command.  Here is a screenshot:

I can now also copy files from my own computer (a windows machine) using the pscp client that comes with my installation of PuTTY.  It works just like scp in Linux. 

2/26/2017

    I have been working on the Kendall machine.  The key to the analysis is varying certain parameters in the input file that tells the MHD what to simulate.  Various properties of the geometry of the tokamak, plasma, and magnetic field can be varied to search for equilibria and determine stability.  Here is a screenshot of the beginning of the input file.

The parameters are given a comprehensive treatment in the documentation for the gs2 gyrokinetics code online at http://gyrokinetics.sourceforge.net/wiki/index.php/GS2:_A_Guide_for_Beginners

I'll give a brief description of some of the relevant parameters.  

-shat: a sophisticated measurement of the twisting of the magnetic field inside the tokamak, or the magnetic shear.  It is also related to the current density gradient in the direction increasing rho

-rhoc: ratio of r/a

-qinp: ratio of toroidal to poloidal cycle numbers.

The input parameters are described in detail at this link: http://gyrokinetics.sourceforge.net/wiki/index.php/GS2_Input_Parameters

I'll include a list of references at the bottom of the logbook.

3/7/2018

Now I am understanding how the input and output parameters are handled by the code.  Here is a screenshot of my runs directory:

The run I am focused on right now is the one entitled dl01.  The dl01.in file specifies the input parameters.  This is where I will be changing relevant quantities to test the conditions derived by Pierre, one of Professor Dorland's collaborators, for stability.  

There are two output files that I am working with right now.  The first is dl01.fields, pictured below:

Each column is a certain physical quantity.  The leftmost column is the z coordinate of the tokamak (the axis of symmetry).  The column to the left is the wave number of the Fourier mode in the y direction.  The column third from the left is the wave number of the Fourier mode in the x direction, all of which are zero as specified in the input file.  The output cycles through each wave number combination for each cycle through the z coordinate.  The fourth column is the real part of the eigenfunction, the fifth is the imaginary part of the eigenfunction, columns 6 and 7 and are the real and imaginary parts of the direction of the field in the z coordinate, columns 8 and 9 are the real and imaginary parts of the amplitude of the field in the z direction.  Column 10 is usually a repeat of column 1, although this can be changed.

The second relevant output file is dl01.out, a sample of which is given below:

This file gives details about the run at the top of the file, listing the day the run was started and the total time it took to run.  The details of the Fourier modes and eigenfunctions are listed in each column corresponding to a specified time in the leftmost column.  The wave numbers are to the left, the omav column lists the real and imaginary parts of omega in the complex exponential in the formula: 

Here is a graph of the real part of the eigenfunctions vs the z coordinate taken from dl01.fields:

Here is a plot of the imaginary part of the eigenfunction vs z: 

3/11/2018

    I am going to do a run with a different nperiod value in the input file dl02.in with the nohup command.  I ran the following command:

This command runs the gs2 code on the remote computer even after I close the terminal with the nohup command.  -np 4 is the command for using four processors. dl02 specifies the input file used by the gs2 code.  Right now the code is running at 1/6 of real time.  That is 1/6 of a minute is simulated for every minute the codes runs.  

3/12/2018

    The dl02 run is done!  Here is the real part of the eigenfunctions

Here is the imaginary part of the eigenfunction:

I am not sure why the peak values are so different between the dl01 and dl02 plots.  All I changed was nperiod from 5 to 10, which is just the number of 2*pi segments of equilibrium magnetic field to simulate (according to the the input documentation).  I am confused about the physical significance of this, so I will ask Professor Dorland about it tomorrow.

I also did a dl03 run, this time keeping nperiod at 5 and changing naky from .01 to .1 with step size of .01.  I am going to keep track of the runs and their input parameters in a separate spreadsheet to keep myself organized.

Here is Re{Phi} for dl03:

Here is Im{Phi} for dl03:

In the dl03 run, the amplitudes are much smaller than in dl01, a factor 10^5 difference in fact.  I'm guessing this is due to the fact that I used different Fourier modes.  

3/13/2018

    To get a better sense of what is going on, I'll make some plots of the data in dl01/02/03.eigenfunc.   

Re{Eigenfunc}_dl01:

This looks fine, as the eigenfunction is normalized in some basis.

Im{Eigenfunc}_dl01:

If this function is odd I think it's also fine.

Re{Eigenfunc}_dl02:

Also looks fine.

Im{Eigenfunc}_dl02:

Same story as last time

Re{Eigenfunc}_dl03:

This looks scary because the end behavior is not very nice.  This might just be because the functions are odd and they were forced to be one at z = 0 because of the normalization, but I am not sure so I will consult with Professor Dorland about this.  

Im{Eigenfunc}_dl03:

This one looks weirder because the amplitudes are not dying away with distance like the real part.  I am not sure what is happening at all here.  

For reference, here is what I did in gnuplot:

3/14/2018

    I met with Professor Dorland yesterday for an hour and we discussed the runs I did and what was actually going on.  I'll get into that in just a bit.  First, I did a derivation of the amplitude (Phi) of the eigenfunctions.  

It turns out that those weird eigenfunctions that didn't decay away with distance have a very low growth rate, as seen in the plot below of gamma versus ky.

For values of ky between .01 and .1, the growth rate is small, even negative at the very end.  The fastest growing modes occur at ky = .4 or so.  They go negative around ky = .7.  This is perhaps the most important aspect of a stability calculation.

3/20/2018

    I am working on plotting the real part of Phi in python on my machine so I can better automate the whole process.  My goal is to write a script on Dorland's computer that will take the output files, send them to my computer, then run my python code to plot all the relevant quantities.  Here is a graph of RePhi from the dl01 run:

3/26/2018

    I have been working on automating the file transfer and plotting process for the last few days.  I can plot a new file with python from the command prompt of my windows machine in only five steps.  

The first command loads my macros into memory.  Here is what the myMacros.txt file containing all my macros looks like.  

I then run the copyK macro to copy a file from Dorland's computer to my machine.  Then I adjust the plotting parameters in the PlotFields.py program I wrote to plot the specified quantities with matplotlib (a python library for object oriented plotting).  Finally, I run the PlotFields.py program to produce the plot I want.

3/27/2018

    I finished the a program to parse the .out files so I can plot the growth rates.  This one is called PlotOut.py, and both are provided at the bottom of the page.  

4/6/2018

    Dorland and I met yesterday to discuss my plotting code and talk about where this project is headed.  He suggested that I format the code so that I don't have to go into the code each time I want to plot from a new file or plot a different variable from a file.  I fixed that and allowed to user to specify the parameters to be plotted from the command line directly.  Here is an example:

I can now even plot the growthrates from one file on the same graph.  The only issue is that I have to know which files contain the smallest wave numbers, because the files have to be listed in order from smallest to largest wave numbers.  Here is the improved plot with titles, no zero compression, and a grid:

4/11/2018

    During out meeting today Dorland showed me how to calculate what is known as "balloon" instability.  The code for evaluating stability uses an MHD rather than a gyrokinetics approach.  With the MHD approach, only perturbations of wavelengths much shorter than the radius of the tokamak device can be accurately evaluated for stability.  The code takes input parameters specified in an input file.  The input parameters can be changed such that one can evaluate the stability for a given radius (a flux surface really), with specified values of beta' and s_hat, the pressure gradient and current density gradient respectively.  Setting the parameter bishop = 4 allows for s_hat to be changed at will and beta' can also be changed to a starting and ending value with a specified step size.  The input file is shown below:

s_hat is specified by changing the value assigned to the variable s_hat_input.  The radius is changed by changing the value assigned assigned to the variable rhoc.  Here is an example of the output of ballooning code, run with the simple command ./ball :

4/25/2018

It turns out that zooming on a region of the s_hat -10 to 10 plot such that -beta' runs from 0 to 2 looks like this:

The region to the left of the boundary is unstable and the region to the right is stable.  This matches how Dorland predicted the boundary should look like.

5/4/2018

    I want to do a point analysis to test the microstability of points that are MHD stable and MHD unstable.  I cannot seem to change s_hat and beta_prime for the gyrokinetic input file with ./eiktest.

5/5/2018

    I have found that the ./eiktest command runs a program that only tests a file with the name "eik.in".  I have found another file that specifies which values of s_hat and beta' are used for any input file.  I can now vary s_hat and beta' at will to do such a point analysis.  

5/6/2018

    I performed the gyrokinetic analysis for microstability on the following two points on the stability plot:

Notice how in the MHD regime the growth rates are all negative, so the amplitudes are decaying in time.  Notice also how the point is MHD stable but experiences gyrokinetic instabilities because the growth rates in the gyrokinetic regime are positive, indicating growing amplitudes at shorter wavelengths.  The conclusion we can draw is that even MHD stable points can by gyrokinetically unstable, so all instabilities have to be checked.

5/7/2018

    I finished my poster and posted it at the bottom of this page.