HERMITE-LAGUERRE BASIS FOR MODELLING TURBULENCE
May 18 - I will be working for Dorland over the summer and we will be continuing this project or something similar. Essentially, if you can run the MATLAB code posted below and have access to Dorland's WIP paper (ask him about his work!), you will be able to recreate everything that I have done. Simply by looking at and understanding my code, you will be able to understand all of the major points that my project focused on.
May 8 - I have been in correspondence with Dorland and it seems that my poster is all good. We will print soon and be ready for the poster session on the 11th. Here it is:
May 3 - I met with Dorland and we talked about what we wanted to put on the poster. Additionally we talked about working over the summer. I am going to work for Dorland and we will figure out what we want to work on in the near future. For this project, we essentially confirmed computationally that Dorland's method will work in the larger code, even if it has been proven mathematically before. Now that I have this MATLAB script (listed below) that can calculate Hermite and Laguerre polynomials, projections, and expansions, as well has probability functions and spectral amplitudes, we can quickly check many aspects of Dorland's code in the future. This ability to check and make sure everything is working properly without having to run a massive and complex calculation is invaluable and will help Dorland to make a publishable set of codes quickly.
Apr. 30 - We more or less wrapped up this section of the code checking and are moving on to other sections now. Soon we will need to finish the paper and get to work on a poster as well.
Apr. 27 - Dorland sent me an email and we figured out what I needed to compute. Essentially, Dorland takes this arbitrary function L(x), and defines it to be the roots of the x-th Hermite polynomial, plugged into all of the lower order Hermite polynomials, and then all of those values plugged into some arbitrary Hermite projection function. This is a completely arbitrary way of checking our codes against each other, as these values should be the same no matter what. As my code is definitely right (Dorland and I have looked over it and confirmed that it produces the right things), we will now check that his larger code generates the same values.
L(3) for phi sub 2 = [ -0.242, 0, 0.242, 0.108, 0.154, -0.154]
L(3) for phi sup 2 = [ -1, 0, 1, 2, 3^(1/2), -3^(1/2)]
These are the values that I calculated and it seems that Dorland's code agrees.
Apr. 21 - We talked more about the issues with his code. Now with my code we are able to check different sections of his code and make sure that different sections of it have the right formulas for the Laguerre and Hermite expansion and projection functions. Furthermore, we will soon check a specific set of values that should be the same regardless of how we actually define the functions, as it will depend on the system and not our different definitions of coordinate systems and function properties. This set of values I will talk about in the next update after we have computed them.
Apr. 16 - I finally met with Dorland again after he got back from his conference at MIT. We discussed a lot about my code and what we are going to do next. First off, we developed some plots for the Hermite functions, which was great to visualize and make sure that what I'm doing is correct. Here's a plot of the first six functions;
This works as expected. I plotted these by using MATLAB and generating the first six polynomials. Then, I used the subs() function to find the y value for each polynomial over the desired range and plotted the xs versus the ys. We made some other plots but they were all just variations of this one. After that we looked as his massive turbulence code and he showed me specifically what he wants me to do.
There are two main problems with his code right now. One is that he is not sure that the Laguerre three term recursion formula is right in his code. That can be quickly fixed by looking at my code for that in Matlab. Secondly, there is an issue in the Hermite functions where there is this factor of sqrt(something factorial) in the denominator that gets absolutely massive as the function gets larger. As a result, the factors of G in the input get tiny and tiny, which can have resolution errors. One way to fix that is to put the factor with the factorial in it into the numerator of G by looking at the function and figuring out the relationship. This is one of the problems that I will need to fix by next week. The other problems are as listed in this email I sent the other day;
1. Plot Laguerre Polynomials (0-5);
2. Be able to switch between physicists and probabilist's Hermite functions
3. Three term Laguerre and Probabilist's Hermites functions/polys
4. G' to G and back (meaning Having the L! in the coefficients or in the function itself)
Hopefully I can accomplish all of these tasks by next week. 1-3 are very simple, but 4 may take some more time as I am unsure what the relationship is right now. The project is going very well.
Apr. 9 - Dorland emailed me so that I could add a few more things to my code. I added a three term recursion generator from the equations listed on the wikipedia article for Hermite functions which was very simple to code. Additionally, I tried to get rid of the huge fractions but they won't go away. I think it has to do with how Matlab handles pi, but for the time being I will sit on what I have until I can meet with Dorland again and we can discuss where to go next.
Apr. 4- The code works! I developed Matlab code to generate Laguerre and Hermite functions using symbolic variables and then I was able to transform to g and back again, getting the same coefficients. I will post the code below. It worked perfectly, and even works for complex values of G and g, which is encouraging. Essentially, we are working through these problems so that I can actually see the pitfalls and process of what Dorland is doing in a familiar environment where I can make modifications easily. Next up is to learn how to take my analytic method for solving Dorland's equations to a discrete one. Essentially, I want to be able to take specific values of g and work to get values of G by using the aforementioned Gauss quadrature method that we will tackle soon.
hermiteScript.m
syms vp uB
ell = 2; % these can be changed to whatever value of Hermite/Laguerre
m = 1; % polynomials that you want to calculate
G = zeros(ell,m);%10*rand(ell,m);
G(2,1) = 1; %this is a test case of the values of the amplitudes, G
hermiteSet = sym(zeros(1,ell));
laguerreSet = sym(zeros(1,m));
hermitefSet = sym(zeros(1,ell));
laguerrefSet = sym(zeros(1,m));
hermitepSet = sym(zeros(1,ell));
laguerrepSet = sym(zeros(1,m));
% hermitePoly = (-1)^ell * exp(vp^2 /2) * diff(exp(-vp^2 /2),ell);
% laguerrePoly = (exp(uB)/factorial(m))*diff(uB^m * exp(-uB),m);
% hermiteSubFunction = exp(-vp^2 /2) * hermitePoly / sqrt(2*pi*factorial(ell));
% laguerreSubFunction = (-1)^m * exp(-uB) * laguerrePoly;
% hermiteSupFunction = hermitePoly/sqrt(factorial(ell));
% laguerreSupFunction = (-1)^m * laguerrePoly;
for i = 0:(ell-1)
hermiteSet(i+1) = (-1)^i * exp(vp^2 /2) * diff(exp(-vp^2 /2),i);
hermitefSet(i+1) = exp(-vp^2 /2) * hermiteSet(i+1) / sqrt(2*pi*factorial(i));
hermitepSet(i+1) = hermiteSet(i+1) / sqrt(factorial(i));
end
for j = 0:(m-1)
laguerreSet(j+1) = (exp(uB)/factorial(j))*diff(uB^j * exp(-uB),j);
laguerrefSet(j+1) = (-1)^j * exp(-uB) * laguerreSet(j+1);
laguerrepSet(j+1) = (-1)^j * laguerreSet(j+1);
end
%these create the Hermite and Laguerre polynomials, functions, and
%probability functions.
gsum = sym(zeros(ell,m));
for h = 1:ell
for k = 1:m
gsum(h,k) = hermitefSet(h)*laguerrefSet(k)*G(h,k);
end
end
%this is what we calculate the expressions for g
%be able to plug in vp and uB and make plots
%large numbers
%complex numbers for G
%plots for both real and complex
%three term recursion for probabilists hermite functions
Glm = zeros(ell,m);
for i = 1:ell
for j = 1:m
Glm(i,j) = int(int(hermitepSet(i)*laguerrepSet(j)*gsum(i,j),uB,0,inf),vp,-inf,inf);
end
end
%this is how we get from the values of g to the values of G. These values
%of G should be the same as the values of G at the top.
______________________________________________________________________________________________________________________________________________________________________
Apr. 1- I met with Dorland at his house and we dove into his new paper and talked a lot about the different equations in there. Hermite functions are simply Hermite polynomials that are weighted depending on some factors. In the Wiki page that I have previously linked you can see the relationship. Dorland uses the probabilists' Hermite functions as his basis functions and has developed a series of transformations to go from G's to g in his work. The G's are the different amplitudes of the different basis sets and the function g is a function of v parallel and uB, which are the two main variables in his work that he needs to calculate. By switching bases, Dorland can drastically reduce the computational time for his code, which is huge because the turbulence is the most complex part of a tokamak, and once Dorland finishes this work the global community will be able to simulate entire tokamaks on a massive scale economically. His code is very flexible, which allows it to be used in a variety of different setups. For example, it can be used at low resolution to quickly approximate turbulence in a full-tokamak simulation in conjunction with other simulations where you don't really need high precision. Furthermore, it can also be used to simulate tokamaks at extremely high precision comparable to state of the art codes, but be run in a FAR quicker amount of time. Dorland gave me a project to try and perform these basis changes in Matlab by writing my own code, which I will do over the next few days. Essentially for this, I will need to generate Hermite functions, Laguerre functions (which are analogous to Hermite functions just slightly different. The difference is not important for the log as they both accomplish similar things. Dorland uses Laguerre functions for his variable uB and Hermite functions for his variable v parallel), and G's, and then transform to a simulated function g, and then back again.
Mar. 17- I met with Dorland individually and we talked a lot more about my specific project. The most important part of my project is the Fourier analysis of the turbulence model. Essentially, by splitting apart the turbulence equations into their separate amplitudes, we can look at which specific frequencies most affect the model equations. This is a way of testing the validity of the set of bases that are chosen for that problem. If the model can be described by a very small number of frequencies, then the basis is a good one as it will drastically reduce the computational time. More specifically, we will be looking at the log of the the amplitudes squared, or log(A^2), so that the values are always positive and we can see if there is an exponential decay in the importance of the frequencies. If there is a negative linear description in this plot, then the basis is good and the equations can be described well by the basis. If there is pretty much any other form of these amplitudes, then the basis is probably not very good.
Mar. 16- We met with Dorland last Friday and for the most part we didn't talk about anything new. We mostly just talked about all the different ways of modeling data that we have talked about (polynomial, Haar, Fourier, hard data) and talked about the advantages and disadvantages of each. These should have direct application to our projects.
Polynomials are great for the evaluation of certain types of data sets and can be approximated very quickly. However, they aren't exact and often can't be integrated/differentiated easily. Hard data is another way of looking at things and is entirely local. This means that using data points allows you to very easily evaluate the value of a function at a single point, but it can't really tell you much about the data on a larger scale. For example, you cant integrate a set of data points. You would need to model it in some other way in order to talk about predictions or other analysis. Fourier is the opposite, as they provide great ways to look at the function on a global scale (by looking at specific frequencies) but are unwieldy on a local level. To be exact, one would have to know every function in the Fourier series in order to evaluate the function at a single point, which is incredibly hard to work with. More importantly, Fourier can always be integrated and differentiated, which makes it great for predictions and estimations. Haar is interesting as its main purpose is for condensing data. The main idea behind Haar is that if both you and the person you are sending data to know the process that you're going to perform on the data, then it can be shown that the data can be manipulated such that parts of the data are extraneous and can be removed. For example if you had a straight line you could give the value of two points and tell the other person that the line passes through those two points rather than giving them the values at say 8 points. This compresses the data which can allow for much faster data transmission and data processing.
We should get into individual meetings soon.
Mar. 4- Once again we met with Dorland. This time we worked on a specific polynomial problem using numerical methods. We were given a specific differential equation using a function that was proportional to itself and were asked to find a function that satisfies the equation. We worked on the board and used some Matlab. The basic process was to create a guess function from the information given at the start, which included boundary conditions and the specification of evaluation points, which are points along the function that you would define for your guess function to be equal to the real function at those points. We created a function that we called the residual which was defined to be the difference between our guess function g and the real function f. Essentially this residual function R was defined to be zero at certain points, which allowed us to make a system of linear equations. Essentially, this method uses an N degree guess polynomial evaluated so that the residual is equal to zero at N points. In doing so, you have N unknowns and N equations, so you can fairly easily solve a system of linear equations.
So we solved the system of equations for a0, a1, and a2, because in our case we had a second degree polynomial where the residual was set equal to zero at three evaluation points. This gave us our finalized guess function. The main idea is that the more points that you set the residual equal to zero, the more precise you can make your guess function, as the number of evaluation points determines the number of terms in your guess polynomial. What we talked about after this was how to make this process better, faster, and more efficient. One way that we talked about was to choose your polynomial carefully. By understanding your boundary conditions you can make smart guesses about your polynomial. For example, in our problem we had identical Y values on either boundaries, which meant that our function had to be an even function over the boundary. In doing so, we could have just made our guess function only made up of the even terms (i.e. constant, x^2, x^4, etc.) which would have given us an even more precise guess. Another way to make the guess more precise is to choose your evaluation points more carefully. The way of doing this is called Gaussian Quadrature. I'll have a link to the Wikipedia article down below. This method is great for choosing points and allows you to go through and choose really good points that will make your function approach the real function far far quicker than using random points. By approach, I mean that it will closely model the real function with fewer polynomial terms.
These pictures describe what we did for the most part. We looked at a function and solved it using numerical methods.
Here is the code that we wrote to solve the system of linear equations and here is the plot of our polynomial function versus the actual, exponential function that we didnt know actually made the model. There are many ways of going about this so this just served as a lesson really.
A short explanation of a Gauss Quadrature
Feb. 25- We met again with Dorland. The main thing that we went over is the idea of basis functions. A big part of what Dorland does is transferring data or computations from one basis to another, which makes it easier to compute, and then switching back again. These basis functions can also be described as coordinate transformations. Dorland takes complicated data or computations for them, finds a way to convert them to a simpler coordinate system, then executes some computation, then switches to the next basis. It turns out that doing this is actually way faster than just solving each issue in one single coordinate system. The way Dorland does this is through matrix multiplication. He uses complex linear algebra to go through and transfer from basis to basis until all his computing is done then transfers to the basis that he wants his data to be in.
What he had us do is take a system of numbers from sin(pi*x/4) from integers [0-7] and try to change the basis so that it was far easier to model than a sine wave.
However, we are still having trouble with it and I will have to update later with what the purpose of all of this is. We are not completely sure what we are doing and we also have not started on our projects yet.
Feb. 21 - We finally met with Dorland and we essentially got introduced to our projects. Dorland models something to do with astrophysical plasma which he models specifically so that he can look at plasma within a fusion reactor, the main design of which is a Tokamak. So Dorland has two main areas of interest - fusion and space plasma. What he specifically investigates for both of these environments is the turbulence associated with the moving plasma. The idea is that turbulence is turbulence everywhere, and that by studying the turbulence of plasma in outer space (which is very easy relatively to collect data for and to analyze) you can make inferences and gain insight into modeling turbulence within a fusion reactor.
Turbulence within a fusion reactor is a very debilitating thing. When the plasma touches the wall of the Tokamak it cools down tremendously. This cooling in turn drastically reduces the efficiency of the reactor as a whole. For nuclear fusion reactors today, the main issue is with efficiency. There are operational reactors today but they aren't efficient enough to use them to power anything substantial. Fusion offers some absolutely incredible benefits for humanity, so if Dorland can reduce the turbulence he can do a lot of good.
Dorland's library of code encompasses four main groups; old fusion, new fusion, old space, and new space. He has two sets of code (in FORTRAN) that can simulate these environments and give a lot of data. However, one problem with these codes is that they are incredibly slow, and require large computing grids and multiple days to actually run. These codes are called GS2 (for fusion reactors) and AstroGK (or Astrogeek, which is for space plasma). He also has some new code written in CUDA + C known as Gandalf which simulates astrophysical plasma with way less processing power (like on a desktop in an hour power) than the old code. He is working on bringing code over from FORTRAN to CUDA + C that will simulate fusion reactors. These new codes take advantage of GPUs in gaming computers to simulate things as they can work incredibly well.
Specifically the project I will be working on will be running diagnostics on some calculations made using AstroGK. He has calculations, he has work done, but he wants to run the code a few more times using different parameter values to see what works and check that he is right. I will be slightly modifying his code to perform these diagnostics and will be using a little bit of Python to go ahead and do some analysis on the output of his code.
Feb. 9 - We contacted Bill Dorland and hopefully he will get back to us soon. In the meantime, I am trying to better understand Dorland's research. I will be going over this website http://accelconf.web.cern.ch/Accelconf/p05/PAPERS/TPAT041.PDF to try and have some idea of what Vlasov-Maxwell equations are, which seem to be really important to Dorland's work.
REFERENCES:
http://www.umiacs.umd.edu/~ramani/cmsc460/Lecture16_integration.pdf - Gauss quadrature
http://mathworld.wolfram.com/HermitePolynomial.html - Hermite polynomials
https://en.wikipedia.org/wiki/Fourier_transform - Fourier Transforms
https://en.wikipedia.org/wiki/Haar_wavelet - Haar Wavelets
https://en.wikipedia.org/wiki/Sturm%E2%80%93Liouville_theory
https://en.wikipedia.org/wiki/Bessel_function
I have kept a personal logbook throughout the semester and I have been periodically updating this one with my information.