FACIT -- compare with NEO

[cjperks7]

Currently working on getting a workflow with both FACIT (within AURORA) and NEO (within MIT-IM) to calculate impurity D, V profiles for C-Mod, SPARC shots. Currently isn't working out so here I'm implementing some suggestions from @DanielFajardoJ on how to fix it. This could totally end up being user error on my part, but in case of otherwise I'll make a branch to track changes in FACIT.

Daniel suggested additions to start:

• (tried second) There is a numerical factor in the definition of the collisionality. As you know different people use different factors in the Coulomb collision times. What is in the python FACIT routine comes from the NRL plasma formulary, but the NEO/GACODE has a different factor. Effectively this means multiplying Tauii inside FACIT by *= 1/(3*np.sqrt(2*np.pi)/4).
• (tried first) The qmag that should be used is not directly the safety factor (I call it qprof in the following) but a generalized qmag = qprof * r_min * B0 / d_torfluxov2pi_d_rmin, or equivalently qmag = r_min * B0 / d_polfluxov2pi_d_rmin. This is effectively the same as the use of 'B_unit' instead of the axis B_t in GACODE terms.
• (tried third) The LFS to FSA transformation of the convection is missing a term. One should compute the factor f_conv_fsa = np.gradient(e0imp, amin*rho)/e0imp and then add += f_conv_fsa*self.Dz_PS, etc, to all convection components.
• This doesn't affect your NEO comparisons, but in Aurora, after one calculates the D and V coefficients from FACIT, there is a geometrical transformation missing. The coefficients that go to run_aurora should be multiplied by D *= drvol_drmin**2 and V *= drvol_drmin

[adrianaghiozzi]

@cjperks7 Came upon your issue as I also am trying to do some comparisons between results from NEO and FACIT. I’ve rewritten the standalone FACIT in Julia and see this when benchmarking for a DIII-D shot with carbon and calcium as impurities. Green is the result from NEO, purple is directly calling facit.py, and red is the Julia version I’ve written (interpolated onto a coarser grid):

I tried a few of the quick fixes you suggest above but they seemed to actually worsen the agreement. Have you had any success since the time you opened this issue with getting better agreement between FACIT and NEO and if so could you share what you found?

Also pinging @DanielFajardoJ - echoing the sentiments of @cjperks7 this could be user error on my part but curious to see if you all have updated suggestions! Thanks in advance!

[odstrcilt]

All these changes were already implemented, but they were somehow accidentally overwritten.
It should be here:
#122

[cjperks7]

@adrianaghiozzi So sorry for not keeping the thread up-to-date! I had some more emails with @DanielFajardoJ about this and we settled on at least a temporary method. Indeed, like you said, when I tried to implement some of the changes to facit.py as @odstrcilt illustrates in the changes under #122 I found the changes caused FACIT and NEO to diverge more. So I'm not sure if there's a conclusion to if those changes are necessary. I haven't convinced myself at least if someone else has.

My question for you is if you're running your NEO/FACIT simulations with a pure deuterium plasma or if you have significant impurity dilution? After further emails with Daniel, we found that the largest reason I was getting a disagreement is that when I set up my simulations I make a guess for the impurity mix to account for the measured Zeff. NEO can handle multiple species, but FACIT has so far been coded for a single main ion species and a trace impurity. So what we had settled on as a temporary method is for me in my FACIT simulations to lump together all of the significant species (i.e. D, H, T, B, F, Ar, etc.) into a lumped species of charge Zeff and mass Aeff (i.e. for a pure 50/50 DT plasma Aeff=2.5, Zeff=1, Zeff*n_lumped=n_e). That seemed to work out the best (see plot below for W50+ diffusion/convection coefficients in a diluted SPARC DT plasma, not implementing changes under #122 ) so is kinda what I personally settled on. Anecdotally, I'll say I've found this level of agreement across a wide variety of C-Mod shots I've recently ran (~15 iirc)

Though again I think there needs to be more discussion on improving the comparisons between NEO and FACIT, especially for diluted plasmas. Thoughts @odstrcilt @DanielFajardoJ ?

[odstrcilt]

When Daniel tried it, it was working very well, here is an example:

There is one more correction which I have not included in my commit

qmag = r_min * B0 / d_polfluxov2pi_d_rmin.

it helps a lot at the edge

The red dashed line uses actual q, the blue dashed line uses r_min * B0 / d_polfluxov2pi_d_rmin, and diamonds are NEO.
V/D is unaffected, but D is much closer

This helps a lot! In this case were not included any additional impurities.

I have also tried to use Zeff instead of main ion Z, Meff = sum (Z_i * M_i * n_i) / (Z_i * n_i) instead of mass and ion density = sum ( n_i) instead of bulk ion density. It has helped a bit, but the results was still too far from NEO.

Here it is when I have calculate q_mag from Bunit definition:

[adrianaghiozzi]

Thanks a lot for the detailed response @cjperks, and no worries at all! To give you a little more explanation on how I’ve been handling the multiple impurity species: in the case I’m showing here, based on an ODS of DIII-D shot #190549, there are two impurity species, carbon and calcium. To get the fluxes shown in the plots there, I run FACIT with different combinations of the species - first, deuterium as the main with carbon as the impurity, and that gives the result in the second plot. Then I run once with deuterium as the main and calcium as the impurity, and then another time with carbon as the main and calcium as the impurity, and I sum the fluxes to give the result you see in the first plot.

In the method that you’re describing where you lump everything together (including I guess the background D?), what do you set as the main ion and the impurity when you set up the FACIT run? I suppose the lump of D, H, T, B, F, Ar etc. is the main and then W is the impurity? If I understand correctly, that would map to lumping D+C as the main and leaving Ca as the impurity in my example? Is there any way to coax the particle fluxes for carbon out of FACIT if you do it that way?

[adrianaghiozzi]

Thanks @odstrcilt as well for the response! Are these plots you show coming from running NEO and Aurora/FACIT in OMFIT? And if so could you point me to a shot number/where you’re getting the profiles from so that I could try to reproduce these results with my Julia implementation? Here’s what I see from just changing the qmag, not including any of the other suggested changes you mentioned - it does give a slight improvement near the edge wrt the NEO result, but it diverges everywhere else (purple = FACIT python with the standard definition of q, red = FACIT julia, now with the new definition of qmag, green = NEO). Note that these are the particle fluxes, Flux_z from FACIT. Perhaps I should look more closely at the individual v and D before comparing the fluxes?

[odstrcilt]

@adrianaghiozzi I have pushed changes here: https://github.com/gafusion/OMFIT-source/pull/7375

NEO module in OMFIT is quite easy to use. You can use a standard interface to create input.gacode file from TRANSP. I think that have used TRANSP from run 175860 Z01 at 2s. And you need to use the Aurora from here #122

[odstrcilt]

@adrianaghiozzi have you made any progress in benchmarking FACIT?