Astrophysical Plasmas

Astrophysical Plasmas/Plasma physics

The topic I will be researching is astrophysical plasmas and working on the Vlasov-Maxwell equation code to solve some plasma interactions. I will be working with Bhas and Shaina on this project.

Sources: 

• 1. Here is a textbook about astrophysical plasmas in general. It goes in depth about many different topics and characteristics of plasmas, most of which are very specific and complicated. I doubt they all relate to what I will be studying. 

  2. http://www.sp.ph.imperial.ac.uk/~sjs/APmaster.pdf

1. Phase Space is an important concept used when analyzing plasma dynamics. Below are links about phase space.

   1. http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/phase.html 

   2. http://www.math.ucla.edu/~tao/preprints/phase_space.pdf

   3. Here is a paper entitled "Numerical methods for the Vlasov equations" which talks about astrophysical plasmas and Vlasov systems. It seems complicated, but very relevant to the code I will be working on. 

   4. http://www-m16.ma.tum.de/foswiki/pub/M16/Allgemeines/NumMethVlasov/Num-Meth-Vlasov-Notes.pdf

   5. Astrophysical plasmas are affected by many things and create very complicated phenomena. Here is a source about Turbulence and Magnetic Fields in astrophysical plasmas. 

   6. http://www-thphys.physics.ox.ac.uk/people/AlexanderSchekochihin/mhdbook.pdf

   7. As far as I understand right now, the main use of studying plasmas is to further the goals of achieving nuclear fusion. Here is an article on fusion research and its relations to plasmas. 

   8. http://fire.pppl.gov/nf_50th_6_Meade.pdf

3/8/16 Found more interesting sources 

• 1. This link explains plasma frequency pretty well and describes how it relates to some real world applications. Very interesting. 

  2. http://newton.physics.uiowa.edu/~rmerlino/129Fall12/29_129_Plasma_oscillations.pdf

  3. Another Good link for general plasma stuff 

  4. http://www.iucaa.in/~dipankar/ph217/contrib/plasma.pdf

  5. A good source that talks about formation of plasmas, some plasma parameters, and application of plasmas in words instead of equations.

  6. http://www.britannica.com/science/plasma-state-of-matter

4/5/16

• 1. This is the paper on the Weibel Instability, a specific instability we will be looking at

  2. file:///Users/stephanebissonnette/Downloads/CalifanoWeibel.pdf

Progress:

2/16/16

• • We had our first meeting in AV Williams with James

  • Attempted to download and set up programs that will be needed for our research. 

    • However, I personally ran into problems when downloading SourceTree, the program keeps freezing when opened up

  • Jimmy began talking about what we will do 

  • The Vlasov Equation:

    • It is a very complicated equation that is used to solve plasma behavior. Until recently we did not have the computer power to solve this equation. The computer technology has increased enough now to a point where we can solve this equation if we are clever in the way we approach it.

    • It is so difficult to solve because it is written in terms of phase space, which is essentially 6D where at any moment the position x,y,z and velocity vx,vy,vz are set.

    • What we will look at are instabilities, in which the growth rate becomes nonlinear really fast. However, there is an interval of time in which the growth rate grows according to linear theory, which we can calculate using our code and then run into into the nonlinear part of the instability 

  • There are many instabilities that plasma formations undergo, but jimmy has made it a goal for us to look at one specific instability and try to solve it. Our whole semester project would be a success if we were able to solve this one instability, and if we do, we could move on to others because there is a plethora of them to choose from. 

    • The instability we will be looking at is the tearing mode instability, which has two forms. 

      • Resistive - less math intensive, we can solve this by hand and jimmy will teach us how

      • collision less - more math intensive, 

  • We also explored the topic of magnetic reconnection. 

    • In a plasma, all the particles are charged and are moving, therefore there are many magnetic fields that form within the plasma as a result of these moving particles. 

    • These magnetic field lines experience tension and so when they come in contact, they snap and release energy, which is free energy that is directly released into the plasma

• • This process creates an instability because when the field lines snap and come back together, it leaves an area of low pressure in the middle which makes the plasma move there and then get shot our everywhere with more energy when the field lines snap again. This process is self-sustaining so it goes out of control pretty fast. 

2/23/16

• • Today we met with Jimmy again to begin learning a few concepts about plasma physics that are going to be useful in our search to solve the instabilities we will be looking at. 

  • The first thing we discussed is what a plasma actually was and a few assumptions we make about this state of matter to study it. 

    • A put together definition of a plasma that we came up with was this: a statistical collection of mobile charged particles

    • A plasma is created when a gas is heated to very high temperatures and the electrons leave their atoms because they are so energized 

    • In a plasma, we assume that a particle's kinetic energy is much greater than its potential energy in relation to any other single particle particle 

      • KE >> U 

    • Because these particles are moving and charged, they are affect by the electromagnetic force

    • Something useful to know is that in plasma physics, temperature is often given in units of energy 

  • Debye Length

  • particles have E fields associated with them, and these E fields have an impact on the surrounding particles 

  • scale that electric field is screened out by plasma

  • based on long derivations involving Gauss’ law, divergence theorem, charge density, etc.

• • 
    

  • Conclusion: In plasma, electric potential and electric field fall off faster than in a vacuum

    • average distance between 2 particles in a plasma with density n

    • from this you can derive the average potential

    • derive average kinetic energy

    • KE>>PE

    • rearrange and substitute in debye length to give plasma parameter>>1

    • number of particles in a debye cube collective effects overpower individual particle interactions

3/1/16

• • Did not meet with Jimmy this week because he had something come up. 

  • Instead we were given a problem to try to solve and present to the class

• • I attempted to solve Part A 

  • To solve for the tCPU, which is the time to numerically simulate the exact evolution of every particles in a single Debye cube (so only one cube not all debye cubes), I reasoned that the formula was going to be (time to simulate two real world particles' time evolution)*(plasma parameter) = tCPU

    • This because we are trying to find how much time it takes to simulate every particle's time evolution and then multiply that by the number of total particles (that's the plasma parameter)

    • To find the time to simulate the time evolution of a single pair of particles, I had to look at the problem more specifically

      • the formula for that was (number of time steps for smallest scale)*(number of floating point operations per time step)*(real world time cost per floating point operation)

    • At first, I did not know what a floating point operation was, but after some clarification, I understood that it was any operation that is performed in formulas calculated by the computer 

    • Jimmy let us know that the equations required to solve the time evolution of any particle were

• • • So I just calculated all the floating point operation in all 3 equations and then I got tCPU = 270*delta*lambda 

3/8/16

• • We went over the problem that was given to us last week. 

  • The time expense does not stem from the floating point of operations when simulating the time evolution of particles, it comes from the number of particles in a debye cube. That number is too large, which is why it would take so long to simulate the exact time evolution of every particle. 

  • If we look at Vlasov equation, we look at the DISTRIBUTION of particles instead of every single particle individually. It is still very complicated to compute, but it is much easier than to find the time evolution of 1E11 particles or so. 

  • We will discuss 3 more plasma parameters 

    • plasma frequency 

    • gyrofrequency (cyclotron frequency)

    • gyroradius (larmor radius)

  • Plasma Frequency

• • • The plasma frequency for electrons is much greater than for the ions because the mass of electrons is much smaller than the other ions. You can think of it as electrons being easier to move around then most massive ions

    • This is the most fundamental time scale of plasmas 

    • It has a wide range of 

• • Gyrofrequency 

• • • Here the velocity is experiencing harmonic motion 

3/22/16

• • Links to look at: JE-14,15,26,28,29

  • Pick one of those links for presentation. 

  • Research the different numerical methods 

    • finite difference 

      • divide up x's into different cells and use the slope between two x's (delta x) to approximate df/dx. Can use taylor series to show that there is an error but has cells get smaller and smaller, error gets smaller. 

    • finite volume 

    • finite element 

  • Jimmy's team utilizes the discontinuous Galerkin finite element method.

3/23/16

• • Looked at JE15

  • Researched Landau Damping 

    • Process in which the oscillations of the number density of electrons exponentially decay 

    • creates a region stability in the parameter space

    • Can be Linear or nonlinear

    • Physical Interpretation 

• • I drew this to describe how the weibel instability basically works. When there is a beam of electrons (Ib in the picture) shot into a plasma, this beam of electrons creates an opposite beam of electrons from the plasma electrons (Ip in the picture) because of the plasmas quasineutrality. In the case where there are no disturbances, this would not do much because the currents would cancel on since they face each other. However, if a disturbance such as a small change in the magnetic field in the z direction (out of the page in the picture) occurs, this would cause a small displacement between the two currents. Now looking at it, the two currents can be thought of as wires that repel and because they repel, they both create magnetic fields in the z direction in the space between them, which then reinforces the original small change in the magnetic field in the z direction. Thus, the magnetic field keeps growing. 

gkeyllall

cd builds/gkeyll/par

make

4/8/16

Dowloaded files sent by Jimmy via email to be able to run simulations. Had to copy the files from my computer to stampede cluster with these commands.

#login to stampede

ssh -Y tg833856@stampede.tacc.xsede.org 

# moved to scratch directory

cd $SCRATCH 

# find the path

pwd 

# copy and past the path and then exit 

exit 

#sftp into stampede

sftp tg833856@stampede.tacc.xsede.org 

# cd to scratch directory 

cd $SCRATCH 

# upload the two files and exit 

mput weibel.lua

mput gkeyll_test_trunk.pbs

 exit 

4/12/16

• .pbs extension is the type of file you use to send jobs to the cluster.

• line that says -n 256 may change

  • always greater than 16 and less than 256    

  • -mail-user line is the email the cluster will contact when the job starts and ends

  • will change the code after the -i 

  • suggestion: 

  • mkdir weibel 

    • • every time you change just make a new subfolder and run the simulation from there. 

Input file: 

• • Poly order (1-4)

• vary k0 (.01-10)

• Resolution, time-stepping will be changed 

  • Nx (64-512)

  • NVX (8-64)

  • NVY (8-64)

  • NVZ (8-64) change all velocities spaces simultaneously

  • phaseDecomp 

    • 4 numbers (all should multiply to number of processors you are requesting)

  • C++ is 0 indexed

After you run  a simulation: 

•    log files 

•  h5 -> hdf5

• have access to "my work directory"/gkeyllall/gkeyll/scripts/gkeplot.py

• -y "weibel"_emEnergy--write-history

• gives you a text file 

• plot text file, time vs log of energy 

  • find the slope 

  • plot then gamma vs some varied parameter 

TO QUEUE SIMULATION: sbatch gkeyll_test_trunk.pbs 

4/13/16

Followed the same steps as above to move more files onto stampede. I just had to change the file names after the mput command. 

mput weibel_1X3V.lua 

mput weibel_1X2V.lua 

I changed the gkeyll_test_trunk.pbs file to be able to submit the jobs that I will do in the future. Here is what I changed in the pbs file: 

• change the first line to #SBATCH -A TG-PHY160007

• change the mail user to my email so I get the emails when the simulation starts and ends running 

• change the ibrun line to: 

  • ibrun -n 256 -o 0 /work/03539/tg828380/gkeyll_executable/gkeyll  -i /scratch/04086/tg833856/weibel.lua

  • The last part after the -i is the path to where the simulation is run. It needs the .lua extension is essential.  The code takes a lua file as instructions, so the lua file must exist and that path must point to the lua file you want to run

  • If you do not put the .lua extension, then the .log and .h5 files will not appear.

The right command to get the graph: 

/work/04086/tg833856/gkeyllall/gkeyll/scripts/gkeplot.py -y weibel_mEnergy --write-history 

Our first graph: 

4/19/16

Process to create subdirectories in which the files of each simulation will be stored (the example is with changing polyorder):

• 1. create subdirectory for overall parameter: mkdir testpolyorder

  2. cd into subdirectory: cd polyorder

  3. create another subdirectory for each simulation run (here is just an example of one those subdirectories): mkdir polyorder1 

  4. You need to make sure the .lua file is in every sub-subdirectory in which you send the simulation results: cp /scratch/04086/tg833856/weibel_1X2V.lua . (make sure you are in the sub-subdirectory, the . is for when you are already in the directory you want the file to be copied into)

  5. Rename the .lua file to whatever you are changing in the parameter: mv weibel_1X2V.lua weibel_2Vpolyorder1.lua

  6. When you are in the sub-subdirectory, pwd to find the entire path to be put in the .pbs file: pwd 

  7. go back to the top directory and change the line of code after -i to the path that you had done in the previous step and add the .lua extension

4/20/16

Jimmy fixed the weibel_1X2V.lua file. I had to delete it from my $SCRATCH directory and re-upload it with the procedure above in the logbook. 

•     PS: change the phasedecomp parameters to {16,4,4}

Next I proceeded to run the simulations with the 4 different polyorder parameters. The runs were finally successful. Here are the graphs

• • The polyorder2 graph is actually the graph with a polyorder of 1 and vice-versa for the polyorder1 graph. 

  • Command to create the graph and text file: 

  • /work/04086/tg833856/gkeyllall/gkeyll/scripts/gkeplot.py -y weibel_2Vpolyorder3_Bz2 --write-history

Now that I have the four graphs of the weibel instability with varying polyorder, I need to transfer the text file data to matlab and create a logEnergy versus time graph in order to calculate alpha (the growth constant), which is just the slope of such a graph. 

The graph of with a polyorder of 1: 

The slope is: 1.098774427637177e-22

The graph with a polyorder of 2: 

The slope is:  1.318429077110039e-22

The graph with a polyorder of 3: 

The slope is:  1.994486976842794e-22

The graph with a polyorder of 4: 

The slope is: 2.777648161063911e-22

Next I plotted the slope (alpha, the growth rate) vs the polyorder number of the simulation. I calculated the slope with a simple slope formula between two points near time = 40 seconds 

4/25/16

I started plotting the variations of k0 for the simulation of the weibel instability. I varied it to be .01, .05., .1, .5, 1, 5 and 10. After running all the simulations and sending the jobs, I got all the text files I needed to create the growth rate graphs. To grab the txt file from my $SCRATCH directory, i used the mget command when you sftp 

sftp tg833856@stampede.tacc.xsede.org 

#then you cd into whatever directory the text file is in

#then you do mget (textfile.txt)

Here are the varyingk graphs: 

k = .01

The slope is 1.654663488277416e-23

k = .05

the slope is  2.110787822455568e-23

k = .1

the slope is  4.222635426849912e-23

k  = .5

slope : 2.182869749817974e-20

k = 1 

slope :   9.768346611665653e-19

k  = 5

slope:  1.556018363995972e-21

k = 10 

slope: 1.442328962553004e-19

Here is the graph of the growth rate alpha vs the varying k0 parameter 

This is the loglog plot: 

BIG ERROR! I have been calculating the slopes of every graph wrong. So I have to redo all the growth rate vs parameter graphs. 

I have decided to restart all of my simulations because the error that I made made everything very confusing, and I want to ensure that my data is correct. 

For the polyorder .txt files: 

/work/04086/tg833856/gkeyllall/gkeyll/scripts/gkeplot.py -y weibel_changingk.05_Bz2 --write-history

sftp tg833856@stampede.tacc.xsede.org 

sbatch gkeyll_test_trunk.pbs 

Slope for polyorder:

polyorder1: 0.140812015016098

polyorder2: 0.179931764069786

polyorder3: 0.176010182212088

polyorder4:  0.222743154092198

Matlab code to plot alpha vs polyorder:

x = [1 2 3 4]

y = [0.140812015016098 0.1799317640697860 .176010182212088 0.222743154092198 ]

plot (x,y)

title('growth rate vs polyorder number')

xlabel('polyorder number')

ylabel('alpha')

Code for Percent change in alpha due to polyorder: 

a =  (0.1799317640697860/0.140812015016098)*100

b = (.176010182212088/0.1799317640697860)*100 

c = ( 0.222743154092198/.176010182212088)*100 

x = [ 1 2 3 ] 

y = [ a b c ] 

scatter (x,y)

title('Percent change in alpha v polyorder step')

xlabel('polyorder step')

ylabel('Percent change in alpha')

Varying k0: 

Graph for k = .01

k = .1

k = .5

k = 5

Slope for k0 parameter

k = .01 -  0.343173420034740

k = .05 - 0.328458704479029

k = .1 - 0.316460292761737

k = .5 - 0.285949203158372

k = 1 - 0.329862958300456

k = 5 - 0.300815419296373

k = 10 -   0.344290769287805

Slopes for k0 with Linear Regression formula:

k .01 - 0.341874781143459

k .05 - 0.324324608952611

k .1 -  0.309622717090245

k .3 - 0.252886508758282

k .4 - 0.267603048078449

k .5 - 0.284309820915140

k .6   0.301115814678802

k .7 - 0.315856048731715

k 1- 0.329597772086216

k 2 - 0.326203428741058

k 3- 0.342775064027326

k 4 -  0.339591311283553

k 5- 0.300207174680946

k 10 - 0.344228381572527