TJLF

Tglf in Julia Learned from Fortran (TJLF)

How to Run

First you must create a InputTJLF struct or Vector of InputTJLF structs. The best way to do this is to call readInput(directory) where directory points at an input.tglf (not .gen) file. Then, you can call either run(InputTJLF) or runtjlf(InputTJLF)/runtjlf(Vector{InputTJLF}).

run() will return you QLweights, eigenvalue, QLfluxout, fluxspectrum

run_tjlf() will return you the summed fluxes aka the array of actual values you care about with like flux matching (see Indices of Arrays for how these arrays are structured)

Internal to these run() calls, there is a beginning, middle, end like I referred to in my presentation: Beginning: <pre>outputHermite = gausshermite(inputTJLF) satParams = getsatparams(inputTJLF) inputTJLF.KYSPECTRUM .= getkyspectrum(inputTJLF, satParams.gradr0)</pre> Middle: <pre>QLweights, eigenvalue = tjlfTM(inputTJLF, satParams, outputHermite)</pre> End: <pre>QLfluxout, fluxspectrum = sumkyspectrum(inputTJLF, satParams, eigenvalue[:,:,1], QL_weights)</pre>

InputTJLF

The InputTJLF struct can be populated with an InputTGLF struct, the translation is done in the tjlf_modules.jl file, but there are parameters missing between the two (see Deleted Parameters). The structure can also be populated with a input.tglf file (not .gen) if you pass the file directory to the readInput() function.

Currently, InputTJLF structs do NOT have any default values (maybe like one or two exceptions) and should throw an error if it is not properly populated, this is to ensure the user is fully aware of the parameters they are using.

New Parameters

InputTJLF mostly follows the input.tglf file variables, but there are a couple new things

FINDWIDTH::Bool - this is in the input.tglf file, but is used different than in the Fortran, I have repurposed this variable to be a flag to tell the code whether or not we are providing an external width spectrum to avoid extra calls to tjlfmax.jl that solves for width (see Eigenvalue Solver)

WIDTHSPECTRUM::Vector - external width spectrum used if FINDWIDTH is False. If FINDWIDTH is True, this Vector will be populated with widths for each ky grid point. The idea is that you can run with FINDWIDTH set as True, it will automatically save the width spectrum, you can then change the gradients of that particular struct, and then change the FIND_WIDTH flag to False. Or, you can easily save this value to other InputTJLF struct, not sure why you would necessarily, but the information is here.

FINDEIGEN::Bool - this flag tells the code which solver to use when solving for the eigenvalues, if unsure, FINDEIGEN as True is more robust (see Eigenvalue Solver). Should ONLY be False if FINDEIGEN is also False because it would be suspicious if you are confident on the eigenvalues enough to use and not confident on the width values (it would also probably break tjlfmax.jl), there should be an assertion to check this.

EIGENSPECTRUM::Vector - external eigen spectrum used on the first pass if FINDEIGEN is False. This value is set to the sigma parameter of the eigs() solver and acts as an initial guess for the iterative solver. Currently, this is a vector that automatically stores the most unstable eigenvalues from the first pass, but it might be smart to save the eigenvalues for each mode and rework how eigs works in LS.jl

SMALL::Float - value used to rewrite the eigenvalue problem as a system of linear equations like the Fortran does (default 1e-13)

Deleted Parameters

I deleted the parameters in the InputTJLF struct that are old and not used. I left comments and deleted them in commit b9a9c99 and commit 18f810f.

Deleted parameters: USETRANSPORTMODEL,GEOMETRYFLAG,BMODELSA,FTMODELSA,VPARSHEARMODEL,WRITEWAVEFUNCTIONFLAG,VTSSHEAR,VNSSHEAR,VEXB,RMINSA,RMAJSA,QSA,SHATSA,ALPHASA,XWELLSA,THETA0SA,NNMAXERROR

Eigenvalue Solver

Currently, Arpack.jl's eigs(), can be used in tjlfLINEARSOLUTION.jl and tjlfeigensolver.jl as a fast iterative solver. In tjlfLS.jl, eigs() is used to solve for the eigenvector instead of rewriting the generalized eigenvalue problem as a system of linear equations like TGLF, which is potentially a little troublesome. It is also used in tjlfeigensolver.jl if you provide a EIGENSPECTRUM and set FIND_EIGEN to False. eigs() does a shift and inverse iteration to solve the generalized eigenvalue problem, specifying sigma gives the shift, by setting which=:LM you tell eigs() to compute eigenvalues around sigma. I have tried which=:LR (LM is largest magnitude and LR is largest real part), but sometimes the solver will find "fake" eigenvalues with very large real and imaginary parts. I have found using :LM and setting sigma as the most unstable mode has worked, but more testing would be good.

Table of the three different combinations of the solvers used in TJLF

+--------------+-----------------------------------------------+
|              |                      FirstPass()              |
+--------------+---------+-----------------+-------------------+
|              |         | ggev!()         | eigs()            |
|              +---------+-----------------+-------------------+
|              | ggev!() |      1.141s     | 522.710 ms        |
|              | eigs()  |     robust,     | robust?           |
|              |         | most "correct", | thread scales??   |
| SecondPass() |         |  thread scales? | best for 1 thread |
| (eigenvalue, +---------+-----------------+-------------------+
| eigenvector) | ggev!() |      1.091s     |                   |
|              | gsev!() |       TGLF      |                   |
|              |         |     robust,     |                   |
|              |         |  thread scales! |                   |
+--------------+---------+-----------------+-------------------+

To run top left, set the InputTJLF SMALL parameter = 0.0, set FINDEIGEN = False To run mid left, use the default InputTJLF SMALL parameter = 1e-13 To run top right, set the InputTJLF SMALL parameter = 0.0, set FINDEIGEN = True

Arpack.jl

NOTE, If you are getting: <pre> Error: XYAUPD_Exception: Maximum number of iterations taken. All possible eigenvalues of OP has been found. │ IPARAM(5) returns the number of wanted converged Ritz values. │ info = 1</pre>

Make sure you are using Arpack v0.5.3 and NOT v0.5.4, the current version does not work. You might have to restart your Julia enviroment, activate TJLF, and build.

Also, Arpack.jl's eigs() is NOT thread safe. I have locks in the code to keep things safe. In the future, GenericArpack.jl should provide a pure Julia version of the Arpack algorithm that is thread safe, but it is still under development and seems to be a ways off.

Indices of Arrays

There are some 3D and 5D arrays where the indices are not obvious. They are specified in the function comments where they appear, but I will repeat them here: <pre>QLweights::Array{5} - [field, species, mode, ky, type] type: (particle, energy, torodial stress, parallel stress, exchange)</pre> <pre>fluxspectrum::Array{5} - [field, species, mode, ky, type] type: (particle, energy, torodial stress, parallel stress, exchange)</pre> <pre>eigenvalue::Array{3} - [mode, ky, type] type: (gamma, frequency)</pre> <pre>QLfluxout::Array{3} - [field, species, type] type: (particle, energy, torodial stress, parallel stress, exchange)</pre> <pre>"Output of run_tjlf()"::Array{2} - [species, type] type: (particle, energy, torodial stress, parallel stress, exchange)</pre> <pre>species: 1 = electron, 2+ = ions</pre> <pre>modes: 1 = most unstable</pre> The order of the indices try to take advantage of Julia's column major memory usage

Multithreading

If you are multithreading, make sure BLAS.setnumthreads(1) is set at some point.

Bash Folder

Wrote a Bash script for testing how number of threads affects speed. Currently it is testing on the seven radial points Tim gave me found in outputs/TIM_test/. You call the executable followed by command line arguements. If you give one number, it gives you the @btime with that number of threads. If you give two numbers, it gives you the @btime looping through those two numbers. If you give three or more numbers, it gives you the @btime for each number in the arguments. I currently don't have any implementation to test the threading speed on other test cases, but you can pretty easily copy my Julia script to meet whatever you need.

Other Notes from DSUN

For SAT0, TM() does not return the value of QLweights (aka "QLfluxspectrumout" in TGLF) and instead returns the value of what TGLF refers to as fluxspectrumout. This is because the saturation rule calculation is done in a different place (LS.jl instead of multiscalespectrum.jl), and otherwise, I would have to pass this extra Array through multiple functions just for SAT0 which is tedious and ugly. Plus SAT0 isn't really used apparently. So, if you are trying to compare the QLweights for a SAT0 run, make sure you keep this in mind.

The name of the SaturationParameters struct is probably not a great struct name since it doesn't really relate to the saturation rules at all and instead are more related to the Miller Geometry, not sure what a better name would be since OutputGeometry is already used.

Currently, InputTJLF's IFLUX parameter tells the code whether it is going to calculate eigenvectors. This is bad and I should probably change it into a local variable instead of a struct parameter. Ran into an issue multithreading with onePass() since it changes IFLUX in the middle of its function. The code is currently fixed with an extra "find_eigenvector" variable, but this should just be the standard IMO. I will try and implement this before I leave and delete the IFLUX parameter, but if this is still in the README I have not done it yet.

Online documentation

For more details, see the online documentation.