d3denergetic/fidasim

A Neutral Beam and Fast-ion Diagnostic Modeling Suite

fidasim's Issues

Refactor rotations

FIDASIM and prefida use the wrong convention for rotations with angles defined to the negative x axis. There is also some device specific routines in fidasim.f90 that relate to rotations. I suggest the following.

1. Right handed coordinate system for beam
2. +x direction should be defined as into the plasma
3. Positive angle == counter-clockwise rotation
4. When going from beam coor. to mach. coor. first rotate about z axis and then about y axis.
5. Replace current rotation matrices with one rotation matrix R = Ry*Rz

Also tilted beam is broken

Calculation switches

@lstagner @bgrierson21 From messing around, I think I found out that if I only want to calculate the thermal D-alpha spectrum due to dcx and halo (to maximise the speed of the calculation I am interested in) I need to set the following switches (if calc_bes isn't selected then the code gets angry):

calc_npa= 0
calc_brems= 0
calc_bes= 1
calc_fida= 0
calc_birth= 0
calc_fida_wght= 0
calc_npa_wght= 0

Is the correct? Is there a dedicated switch for the thermal DCX and halo D-alpha emission or is that done by default if one of hte above switches (i.e calc_bes) is selected?

Should I continue to post questions on github instead of emailing them?

Thanks

NetCDF cannot handle large file sizes

When trying to calculate a the weight functions on a 100x100 energy pitch grid for 41 fibers NetCDF cannot write out the array and prints out the following error.

 NetCDF: One or more variable sizes violate format constraints

I am fairly sure this could be solved by moving over to HDF5 as mentioned in Issue #15

write_spectra fails if NPA chord appear first

Use LAPACK

There is no good reason not to use LAPACK. Its routines are heavily optimized and we will probably get a nice speed bump if we switch over to it.

Distribution Specification

Add distribution_type field to prefida input file.

Same Grid, Multiple Beams

In the case of having multiple beams sharing the same grid (like in NSTXU), it would make sense to allow the specification of multiple beams in some way.

The way I can see this working is that prefida will create input files for each specified beam; appending the beam name onto to runid e.g. <RUNID>_<BEAM_NAME>_inputs.cdf if there is more than one beam.

We can also create a post-processing routine that combines FIDASIM runs or just combines the *_neutrals.cdf files to be used in FIDASIM calculations.

To implement this someone would have to

• Change the nbi structure specification in the device routines to include multiple beams and beam names (also update docs to reflect changes)
• Change prepare_beams prefida subroutine to handle multiple beams
• Modify prefida to loop over the beams when writing files
• Optional: Create post-processing routine that combines FIDASIM runs or just combine the *_neutrals.cdf files

Doppler shift and wavelength axis

This is an issue regarding the definition of the wavelength axis and where on the lambda grid photons are accumulated. In the IDL version of FIDAsim the photons were accumulated on wavelength grid "boundaries", not grid "centers". In order to properly interpret the doppler shift from the halo spectrum the doppler shift had to be calculated with the wavelength axis shifted by "delta lambda" divided by two, where "delta lambda" is the wavelength grid spacing. I have run cases with flat ne, Te, Ti, Zeff profiles and zero rotation. Three runs were executed with wavelength nlambda=768 (default), as well as double and half of this value. If the halo spectrum is turned into a doppler velocity using the output "lambda", then the velocity is finite and scales with the nlambda, which is incorrect. If the halo spectrum is turned into a velocity using lambda+dlambda/2 then the velocity comes out to approximately zero for all three cases. Therefore, the proper wavelength axis for plotting and fitting of FIDASIM output is to define lambda = result.spectra.lambda + (lambda[1]-lambda[0])/2.0. Comments and/or confirmation that this is the proper interpretation is welcome.

Object Oriented based modularization

The current method of modularization is to extend the idl path to include a directory and look for routines that start with "device_name". This is not ideal and could break easily. A better approach is to leverage objects. First create a base Device Class with some standard methods and variables. All other devices, e.g. D3D,NSTX will inherent from this base class. When prefida runs It will instantiate the the proper device class. This will make for code that is harder to break.

Resolving versions of FIDASIM

We have two versions of FIDASIM. The version hosted here and Bens version. Both have diverged and we need to reconcile them. Here is a list of Bens changes that may want to incorporate.

1. Misc. Changes
2. Specific routine for passive fida (includes special background neutral file that can be loaded in)
3. CX - Beam thermalization table for halo calculation (faster and more accurate apparently)
4. Instrumental response function in FIDA weight function calc.
5. Interpolation

We will definitely do 0,2, and 4 but I am unsure about the rest

change reading order in fidasim.f90

from
read(66,) inputs%ab
read(66,) inputs%ai
to
read(66,) inputs%ai
read(66,) inputs%ab

Tables in netCDF format

Binary files are a pain to deal with and everything else in FIDASIM uses netCDF

Version 0.4 Planning Issue

This issue if for planning what features we may want in v0.4.0. This version will break backwards compatibility.

New Features

This is a list new things that would directly affect the user.

• Interpolation of profiles and fbm in fidasim.f90: With this we should be able to calculate the bremsstrahlung and npa attenuation without having the beam grid extend past the separatrix along the sightline. It also opens up the door for orbit calculations.
• load_plasma: We lose the ability to handle 3D fields when using interpolation. This switch would allow for users to load up a file that contains the profiles and fields on the beam grid. It would act similar to load_neutrals switch.
• Use stored velocity vectors for NBI, DCX, Halo calculation. Also would allow for multiple beams. (Issue #39, #19 , and discussion)
• Profile, FBM, Equilibrium, and Tables specification: Users should be able to specify the format of these files so they can be transformed into a form FIDASIM can understand. (Issues #31 #30 #32 #42)
• HDF5 file format: HDF5 has better features and language support. If we are going to break backwards compatibility we might as well destroy it. (Issue #15 #44)
• Improved Documentation: This just needs to be done.
• Apertures for NPA and NBI (Issue #50)
• Dump halo generations (Issue #56)

Interface Routine Improvements

These are things that improve how FIDASIM would interact with other software/databases.

• EFIT/GEQDSK Tools: Ben has improved the current routines. Related to Issue #30
• TRANSP Tools: We have some tools to extract all the information needed for FIDASIM to run (beam geometry, profiles, equilibrium, fbm) but they need to be streamlined. Related to Issue #32 and #21
• FBM Tools: We need routines that can read distributions produced by other codes into a format FIDASIM can understand. Related to Issue #32
• Table Tools: We need routines to make tables for other tables from other databases i.e. ADAS, PREACT, ... Related to #42

Maintenance

Bug fixes and things users wont notice.

• Resolve differences in core analysis routines (Issue #41)
• Tables in new format: It is sometimes necessary to inspect the tables. The current .bin forms is a hindrance in that endeavour. (Issue #44)
• Beam grid and rotations: fidasim.f90 uses the wrong coordinate system (-x into plasma, instead of +x ) because of this the angles alpha and beta variables in the prefida input file are not intuitive. Tilted beams are also broken. (Issue #43)
• Right now prefida checks the inputs to make sure they are valid. With the upcoming changes it would be prudent to do the same for the device routines outputs (Issue #36)

Comments?
@bgrierson21 @geigerb

Backup Geometry file

There should be a backup file that contains the geometry of the beams and diagnostics in case device routines that get the most up to date values are unavailable e.g. when running a ASDEX case at D3D.

Need to launch beam neutral at ion source location, not beam aperture location

The current FIDAsim version launches beam neutrals at beam aperture location, not the ion source location. For a perfectly collimated neutral beam, this is not a problem. However, for any beam with divergence, it induces inaccuracy in beam neutral density profile and magnitude. One example is shown in Fig. 1 and Fig. 2. The beam and halo neutral densities spread much wider and the peak magnitudes are much smaller when launching beam neutrals at ion source location. They will not only change the FIDA light intensity and also modify the FIDA spatial profile. For a NSTX case, it shifts the peak of simulated FIDA spatial profile inward by 5cm, which make the simulation results closer to experimental measurements.

Fig.1 Comparison of beam and halo neutral densities along the NB centerline when launching beam neutrals at beam aperture location (solid lines) and ion source location (dashed lines)

Fig.2 Comparison of beam and halo neutral densities perpendicular to the NB centerline when launching beam neutrals at beam aperture location (solid lines) and ion source location (dashed lines)

RFC: Aperture Specification for NPA and NBI

We currently define a NPA diagnostic by 5 parameters ( 2 points in space, radius of detector and aperture, and the separation between the aperture and detector). This is flawed in a few ways

1. Only supports detectors with circular apertures and detectors
2. Disregards orientation of the detector (i.e. is the detector tilted w.r.t. the midplane)
3. Assumes the aperture and detector are parallel

I propose that we specify an aperture/detector by 3 points in space. Like so...

Using these three points we can also define a set of basis vectors for the aperture/detector plane. With the basis vectors it is easy to go between plane coordinates and machine coordinates.

To determine if a particle goes through the aperture and hits the detector we find the point on the plane where the particle crosses the aperture and detector planes. We then use an indicator function to determine if the point lies within the aperture/detector region. This can all be done in a single function. Here is an example

So the full NPA diagnostic specifications would be

• Three points for the aperture (like pictured)
• Three points for the detector
• Aperture shape e.g. "ellipsoidal" or "rectangular"
• Detector shape

Another advantage with this specification is that with the infrastructure we develop to model the NPA aperture/detector could also be used to model NBI apertures.

Cross section tables specification

It would be nice to be able to specify which tables to use in prefida e.g. "tables":"ADAS" or "tables":"PREACT"

I think this would be done when we switch the table format to be netCDF files.

Need to expand the cross section tables to lower energy and lower plasma temperature

The lower energy limit of CX cross section table is currently set as 0.5keV/amu. It is desirable to expand the cross section tables to lower energy and lower plasma temperature. I have two test runs for a NSTX discharge, one with the standard FIDAsim cross section tables and the other one with new cross section tables expanded to much lower relative energy and plasma temperature (see Fig.1). Figure 2 clearly shows halo neutral density is overestimated when using the standard cross section tables. The beam neutral densities from two runs are the same. The reason is simply because E_rel (in the order of plasma temperature, i.e. 600-800eV in this case) in the calculation of halo neutrals is smaller than the lower energy limit of CX cross section table and sig*v term is overestimated.

Fig.1 Comparison of standard FIDAdim cross section tables and new expanded cross section tables.

Fig.2 Comparison of beam and halo neutral density along the neutral beam centerline. The dashed lines are from the run with standard FIDAsim cross section tables, while the sold line are from the run with expanded cross section tables.

Testing Framework

There needs to be an testing framework to verify that everything is working as it should when changes are made to the source code.

Automatically build and deploy documentation

Lets be honest. When I move on the documentation will never be touched again. So lets just have the robots take over this job.

Table Specification

In the future we may want to be able to switch between different tables e.g. something like "table_type":"ADAS"

Make source code citable via DOI

https://guides.github.com/activities/citable-code/

HDF5

Currently the source code reads and writes netcdf3 files. Netcdf3 files are now depreciated with netcdf4 which is based off HDF5. Netcdf4/HDF5 files allow for grouping within data files which will help with organization. I suggest switching to HDF5 since it is more open then netcdf4 files and has better language support.

JSON input file

Currently the input file for prefida is actually a procedure that is called. This is a very hacky way imputing parameters. This also means that the input files cannot be read by other languages. I suggest switching over to a json input format. This would allow for other languages to use the same input files. Currently json is supported in IDL > 8.2 and is widely supported in other open source languages.

D3D tilted beam

The d3d routines need to be able to handle a tilted beam.

Indicate if fidasim was built using uncommited files

This is just a reminder for myself that I need to add this functionality.

[[ -z $(git status -s) ]]

Add JET Routines

It was requested that we add JET routines. This was already partially implemented in the ls/srs_bench branch in commit 8c2cb0e but everything in that commit is hard coded. A person more familiar with JET will have to make them more general.

New geometry setup

@lstagner @bgrierson21 Previously for the 330lt beams on DIII-D, I/Brian had the geometry as :

origin:[-139.49, 233.77, 0.0]
alpha: -1.46293
beta: 0
gamma: 0
nx,ny,nz: 96,48,50
xmin,xmax:35,95
ymin,ymax:-50,50
zmin,zmax:-70,70

For some reason this causes prefida to say: ERROR: Beam centerline does not intersect grid.

This is what you have in run_tests:

   ;; Non-rotated, Non-tilted Grid
   grid01 = {nx:50,ny:60,nz:70,$
             xmin:-50.d0,xmax:50.d0,$
             ymin:-230.d0,ymax:-110.d0,$
             zmin:-70.d0,zmax:70.d0,$
             alpha:0.d0,beta:0.d0,gamma:0.d0,$
             origin:[0.d0,0.d0,0.d0]}

The obvious differences are the origin, and the symmetry in xmin,xmax. Has something else about the geometry changed between versions? Do you have working examples for any of the DIII-D beams so I can compare?

Thanks

Profile Specification

Prefida assumes that profiles are idl save files and are in a certain location. Prefida input file should include

"profile_dir":"/path/to/profile/dir"
"profile_type":"GAPROFILES"

or

"profile_type":"TRANSP"

Automatically install dependencies

Its a pain to deal with the dependencies. We should include the libraries in the source code and have a script to compile them.

Change how sampling the fast-ion distribution is done for massive speedup

The FIDA and NPA simulations are the slowest part of the code. Many of the grid points are not even seen by the detectors. If only the grid points that contribute to the signal are simulated that would lead to a massive speedup. Right now that is not possible since when the fast ion distribution is sampled it moves the particle a larmor radius away. So a particle from a cell that does not contribute could move into one that does.

The solution is to start at contributing grid points and sample the fast ion distribution a larmor radius away essentially reversing the process.

Need more EFIT routines

The test routines need more EFIT files then we currently have

Rename `species_mix` to `current_fractions`

To remove the ambiguity.

Audit of parallelism of fidasim.f90

Parallelism comes with computational overhead. In some cases running things in parallel slows down the code. I suggest an audit of the parallel portions of the code to determine if the parallelism is actually speeding up the code.

Convert Integers to long in prefida

Specifically the number of spectral channels. Since one could use FIDASIM to simulate a camera which has a lot of channels

Add beam energy index to birth profile

This would make it easier to filter out the full, half, and third energy particles

Feature Request: Option to dump DCX (0th generation halo) neutral density

Here is a snippet from an email conversation describing the issue.

Sorry I had thought you wanted the neutrals from the first generation halo (Deyong wanted this when he was testing FIDASIM's Halo model with TRANSP's)

First lets define our terms (D_0 = Beam Neutral, D_th(+) = Thermal Ion, D_f(+) = Fast Ion)

D_0 + D_f(+) -> D_0(+) + D_f # FIDA Neutral (Produces FIDA light)
D_0 + D_th0(+) -> D_0(+) + D_th0 # DCX (0th Generation Halo) Neutral (Produces Halo light)
D_th0 + D_th1(+) -> D_th0(+) + D_th1 # 1st Generation Halo Neutral (Produces Halo light)
D_th1 + D_th2(+) -> D_th1(+) + D_th2 # 2nd Generation Halo Neutral (Produces Halo light)
...and so on

FIDASIM writes writes out the (0th-Nth) neutrals into the halodens variable. The DCX(0th Generation Halo) neutrals are calculated during the dcx subroutine. Subsequent halo neutral generations are calculated in the halo subroutine. If you just want the DCX(0th Generation Halo) neutral density and spectra then you just have to comment out the call to the halo subroutine located here. If you want the 0th-Nth generation neutrals and spectra then you would change the number of iterations (1-N) the halo subroutine does.

I see no reason not to include an option to do this but this brings up a number of questions that need to be answered.

• Should we merge the functionality of the subroutine dcx with subroutine halo. They both essentially do the same thing (with some minor differences that can be handled by a if statement)
• Should we dump the density for every generation or just the 0th? What about just dumping the ith generation?
• Should I start calling the 0th generation the 1st generation?

@bgrierson21 @shaunhaskey @geigerb

Check device routines outputs

We already do this for the prefida inputs and it has made it easier for users to update to a different version. Since the device routine outputs may change in both variable names and format we should check the outputs to make sure they are still valid

misnaming of {x,y,z}_grid in neutrals/grid?

@lstagner @bgrierson21 I was trying to get all of my old FIDASIM plotting routines to work and track down the differences in the intensities when I noticed that in the neutrals.h5 (and dcx) file the values in grid/x_grid appear to be in machine coordinates (so shouldn't it be named u_grid instead?). Same with y_grid, and z_grid. On the other hand the values in the 1D array grid/x are in FIDASIM co-ordinates, so they seem to be correctly named.

git describe --tags --always --dirty always returns dirty when run in makefile

Intensity issue

@lstagner @bgrierson21 I have been trying to chase down the intensity difference between main ion CER simulations using the old and new FIDASIMs. It seems as though the issue might be related to the beam ion density. I have plotted the n=3 neutrals from old FIDASIM (right plot) and new FIDASIM (left plot). The top row shows the neutrals from fdens, hdens and tdens. The bottom row is DCX+halo, DCX, and Halo only. As far as I can tell the pinj and einj were the same for both simulations. We need to scale the fdens by a factor of about 2.2 to get to the old fdens. The same factor appears to be required for DCX+halo (what gives the main ion CER intensity).

I'll keep checking to make sure the inputs are the same, but unless I have made a mistake, it seems that there might be an issue with scaling of the beam deposition...

Replace IDL with python

IDL is closed source and not everyone has a license for it. I suggest moving over to python. Python is open source and is becoming increasingly popular in science. It also has excellent libraries for plotting and numerical routines (matplotlib,numpy,scipy).

Only write NPA files when we have NPA channels

Error Estimation

We need to be able to propagate errors in the plasma profiles to the calculated spectra. The way to do this is to do the simulation with plasma profiles that are sampled according to their errors. You do this many times and this would build up a posterior distribution. It would take to too long to do this with the Monte Carlo routines but it can be done with the weight functions.

This project should be rather straight forward. I believe this would be a good project for an undergraduate to undertake. The following changes will have to be made

1. Input error estimates for plasma profiles into FIDASIM
2. Change colrad routine to take dne,te,ti ... arguments instead of grid cell index
3. User input switches e.g. calc_error, num_samples, error_type...
4. Routines to handle different error types e.g. Gaussian, Poisson, systematic....(Gaussian probably only one needed)
5. Required changes to the FIDA and NPA weight function routines
6. Output. Should the W.F. be outputted for every sample or just the mean and standard deviation? Perhaps have the output style be chosen by the user. Should we output the sampled parameters as well?
7. How will we propagate the errors in the weight functions to the spectra? If the posterior for each E-p grid point is Gaussian then we can just use standard error propagation. If not then we will have to repeat the above procedure for the spectra calculation.

Once this is in place we can then do a detailed sensitivity study

Routine to create input files from TRANSP run

A routine that pulls all the necessary info from a TRANSP run would be helpful for users that do not have access to a projects resources(devices specific routines, data servers)

Chords use too much memory

Cameras use a lot of chords. When you have too many FIDASIM will stop mainly due to this allocation.

allocate(spec_chords%dlength(spec_chords%nchan, &
                                 beam_grid%nx, & 
                                 beam_grid%ny, &
                                 beam_grid%nz) )

Most of the elements in this array are zero so there is a smarter way to do this.

Negative number of halo particles?

@lstagner @bgrierson21 Negative number of halo particles?

Here is some output:

halo:    16:07:31
     # of markers:    452215
     # of markers:    204975
     # of markers:     80807
     # of markers:     16331
     # of markers:    -18107
     # of markers:    -35284
     # of markers:    -39054

What do the negatives mean?

Have FIDASIM fail gracefully when there are no dcx neutrals

Equilibrium specification

Explicitly state equilibrium file location and type in prefida input file.
Right now prefida makes too many assumptions about the type and naming convention of the file.

"equilibrium_file":"/path/to/file"
"equilbrium_type":"GEQDSK"

In the case of getting equilbrium from server

"equilibrium_file":"EFIT01"
"equilibrium_type":"MDSPLUS"

Bi-linear Interpolation of FBM instead of nearest neighbor averaging

Prefida does this when calculating fast ion density on the grid. It would be best to do this in FIDASIM as well.