D3DEnergetic / FIDASIM

A Neutral Beam and Fast-ion Diagnostic Modeling Suite
http://d3denergetic.github.io/FIDASIM/
Other
27 stars 18 forks source link

ASCOT5 interface #231

Closed mabaquero closed 10 months ago

mabaquero commented 3 years ago

Dear experts, I am looking for a way to correctly format ASCOT5-generated distributions to use as inputs for FIDASIM-2.0.0.

I have up until now used a Python script included in the ASCOT5 release called a5fidasim.py that builds the required dictionaries from the a5 HDF5 output file to feed into fidasim.preprocessing.prefida(). The dictionary corresponding to the fast-ion distribution, dict[], is built upon choosing dist["type"] = 1 and formatting the structure a5.active.dist5d.get_E_xi_dist() (read from the a5 file).

I can run fidasim.preprocessing.prefida() just fine, and I am able to successfully run FIDASIM using the input generated by it. The problem is that, upon comparing synthetic spectra generated for the same shot/time using TRANSP/NUBEAM+FIDASIM, I see a difference by a factor of ~5e5 (the a5 data being smaller).

I am almost sure that the problem comes from a bad scaling that I am doing when manipulating the a5 distribution to format it for prefida.
Do you by any chance have a tested Python script that creates the FIDASIM dist dictionary from a5.active.dist5d.get_E_xi_dist() that you can share with me?

Thanks, Marcelo

lstagner commented 3 years ago

@kgage You are nominally our ASCOT guy any ideas?

lstagner commented 3 years ago

@mabaquero a quick check you can do is to integrate over the input fast-ion distribution F(E,pitch,r,z) to find the total number of fast ions e.g.

image for axisymmetric distributions. This should match what ASCOT gives.

Also I suspect that the factor of 0.5 is due using the wrong pitch units see the following from the comments of read_nubeam.pro which converts a TRANSP/NUBEAM distribution to be compatible with FIDASIM

;;------------Convert d_omega --> pitch
    ;; Fast-ion distribution is given as a function of cm^3, energy
    ;; and d_omega/4Pi. omega is the solid angle in 3D velocity space. In
    ;; order to transform this to a function depending on pitch instead
    ;; of d_omega/4PI, one has to multiply by 0.5!
mabaquero commented 3 years ago

Thanks a lot @lstagner . I will continue debugging my script then.

I will keep the issue open in case anyone else wants to contribute ;)

kgage commented 3 years ago

@mabaquero I'm assuming you're using the version of a5py.postprocessing.a5fidasim.write() in ASCOT's master branch.

I don't see any obvious issues with converting meters and eV into FIDASIM's preferred units, but I'm a bit confused by your volume element Lines 38-43:

    # convert distribution to physical density
    vol = np.pi * (dist_E_xi["r_edges"][1:]**2 - dist_E_xi["r_edges"][:-1]**2) \
        * (dist_E_xi["z_edges"][1] - dist_E_xi["z_edges"][0]) \
        / (dist_E_xi["r_edges"][1] - dist_E_xi["r_edges"][0]) \
        / (dist_E_xi["phi_edges"][1] - dist_E_xi["phi_edges"][0]) \
        / (dist_E_xi["z_edges"][1] - dist_E_xi["z_edges"][0])

It appears you use CodeCogsEqn (1) Isn't a cylindrical volume element CodeCogsEqn (2)

mabaquero commented 3 years ago

Thanks for pointing this out @kgage

My script is indeed based on a5py.postprocessing.a5fidasim.write() in ASCOT's master branch - the part that you mention is identical in my script and ASCOT's master branch.

I will try correcting it and seeing whether it does the trick.

Thanks for your help!

mabaquero commented 2 years ago

Hi @lstagner and @kgage, I finally had some time to work on this issue. I have implemented a transformation of the ASCOT distribution in such a way as to obtain a function of E, p, r and z that fulfills the normalization suggested by Luke above (see entry on Feb. 16). I am testing it by creating the FIDASIM input files using my script, and then using "plot_inputs" in lib/scripts.

I wanted to ask you a couple of questions: 1) Is "plot_inputs" the right script to use for this purpose? 2) Is there another simple script that allows me to access the distribution as handled during the execution of FIDASIM? For example, I am wondering if there is anything that works with types like FastIonDistribution as defined in https://d3denergetic.github.io/FIDASIM/type/fastiondistribution.html . This would allow me to immediately obtain, for example, the total number of fast ions "seen" in FIDASIM to directly compare it with what I have from ASCOT.

The reason why I am asking is that the plots I get using plot_inputs look slightly different from what I get when I integrate the distribution function in my script. I suspect that "plot_inputs" does not include the extra factor of 2pi*r when integrating over r (is this correct?), but I am wondering if I am forgetting something important.

Thanks a lot, Marcelo

mabaquero commented 2 years ago

@lstagner , @kgage Today I noticed that, when running FIDASIM, the info that is printed on the screen during execution may be telling me what I want to know. Indeed, under "Fast-ion distribution settings", there is a field called Ntotal.

Is this the total number of fast ions computed in FIDASIM? In other words, is it the result of the integration of the fast ion distribution function within FIDASIM?

Thanks, Marcelo

lstagner commented 2 years ago

Yes, Ntotal is the total number of fast ions. It is calculated by integrating over the fast-ion distribution.

mabaquero commented 10 months ago

Dear experts,

A few weeks ago I attended the EPS conference on plasma physics in Bordeaux, and had the chance to talk with Bill (Heidbrink) about FIDA in general, and this issue in particular.

I told him that I modified the Python script that performs the preFIDA (preparation of files to run FIDASIM) using ASCOT5 .h5 result files as input, to solve the issue described above.

I have tested the code at SPC (on TCV shots) with good results - for now. Nonetheless, it has not been tested by other groups. Also lacking is a benchmark against TRANSP/NUBEAM + its preFIDA.

Would you be interested in getting the code for further testing?

lstagner commented 10 months ago

Sorry for the late response I just got back from vacation. New code would be appreciated. If you want you can either email me the code or open a pull request. I can also make you a contributer on Github so you can push directly to the repo though if you do that make sure to do it on your own branch.

mabaquero commented 10 months ago

Thanks @lstagner ! I will email you the code.

lstagner commented 10 months ago

A quick question about ascot. Your write function takes in a a5fn variable how is this variable created. Does Ascot produce output files that are read into this variable. Ideally I want to turn this function into a script where I just point it at the ASCOT output files.

mabaquero commented 10 months ago

Hi @lstagner a5fn is just the filename (a string of characters) of the .h5 file created for the ASCOT5 run (and where the ASCOT5 results are stored).

This means that the script gives you upfront the functionality you seek ;)

BTW I usually use this script within a Python script that implements a workflow to run ASCOT5 + FIDASIM 2.0.0 on a specified TCV shot/time. If it is useful for you, I can send it and/or we can do a short Zoom meeting to go over it.

lstagner commented 10 months ago

I haven't run ascot5 myself yet. Could you send me an example output file so I can make sure everything works.

mabaquero commented 10 months ago

Sure. I sent a link to your email so that you can download it (the file is pretty big).

One comment: My script uses some tools in the Python module of ASCOT5, a5py. This includes the routines that open the .h5 and read the data. Therefore, it is necessary to have a5py installed beforehand - I hope it is OK for you.